Open Access
Issue
A&A
Volume 682, February 2024
Article Number A167
Number of page(s) 16
Section Numerical methods and codes
DOI https://doi.org/10.1051/0004-6361/202347007
Published online 23 February 2024

© The Authors 2024

Licence Creative CommonsOpen Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

This article is published in open access under the Subscribe to Open model. Subscribe to A&A to support open access publication.

1 Introduction

Inferring Cosmology and AstRophysics with Observations of Gravitational Waves (ICAROGW1) is a pure python code developed to infer the population properties of compact binary coalescences (CBCs) observed with gravitational waves (GWs). The problem of inferring the population properties from a sample of observations of astrophysical sources is a very common and longstanding problem across a range of research topics. With almost 100 GW observations from the last runs of the LIGO, Virgo, and KAGRA (LVK) collaboration (Abbott et al. 2023a), GW sources are moving rapidly to the “population-inference” domain. The first distribution of ICAROGW was presented in the context of GW cosmology in Mastrogiovanni et al. (2021) and first used in the LVK analysis of Abbott et al. (2023b) and released on the LVK official GitLab page 2. The first version of ICAROGW was also used in independent studies for population code validations (Karathanasis et al. 2023b; Turski et al. 2023) beyond general relativity (GR; Leyde et al. 2022), astrophysical processes (Karathanasis et al. 2023a), and primordial black hole models (Liu et al. 2023; Zheng et al. 2023).

For astrophysics, population inference involves the correction of a selection bias, namely, the Malmquist bias (Malmquist 1922) that prevents the observation of a particular class of astrophysical processes, thus biasing the population analysis when not properly taken into account. Selection biases appear in many astrophysical observations involving neutrino physics (Loredo & Lamb 2002), exoplanets (Foreman-Mackey et al. 2014; Winn & Fabrycky 2015), galaxies and galaxy clusters (Gonzalez & Faber 1997), and γ-rays (Loredo & Wasserman 1998). For GW observations, selection biases are introduced by the sensitivity of the detector as a function of the GW frequency, which is also related to the binary parameters such as masses and redshift.

Besides the correction for the selection bias, GW population inference also has to account for possible uncertainties in the measurement of the source parameters, since the GW signal is always observed in addition to noise in the data (heterogeneous data). Therefore, when reconstructing the population properties of GW signals, we face noisy, heterogeneous, and incomplete observations that require a specific statistical framework. Bovy et al. (2011) also referred to this type of analysis as extreme deconvolution. Current techniques for population inference include the use of hierarchical Bayesian inference (Mandel et al. 2019; Vitale et al. 2022). The code provides a user-friendly python environment to work with hierarchical Bayesian inference. As we will discuss later, only a few “ingredients” are required for population inference: (i) a set of parameter estimation (PE) samples from the finite number of GW observations, (ii) a set of simulations (or injections) to calculate the detectable volume in parameter space and (iii) a rate, or parameter population distribution.

With respect to its first release, ICAROGW 2.0 includes several useful improvements for all types of population studies. The code incorporates a user-friendly environment that allows the implementation of new population models without entering into the code details and contains utilities to study the numerical stability of the Bayesian hierarchical inference. For GW population studies, ICAROGW 2.0 contains several cosmological and population models that are commonly used in previous literature. In Appendices A and B, we describe the CBC population models for spins, masses, and redshift as well as cosmological expansion and modified gravity models. It also contains three methodologies that can be used to perform GW cosmological studies, namely: (i) the spectral siren (ii) galaxy catalogue, and (iii) the electromagnetic (EM) counterpart method.

This paper is organized as follows. In Sect. 2 we summarise the theoretical and code implementation basics of hierarchical Bayesian inference. In Sect. 3, we briefly discuss the infrastructure of the code and its performance on GPUs. In Sect. 4, we describe our process of validating the code by reproducing the population and cosmological results published in literature from the third Gravitational Wave Transient Catalogue (GWTC-3; Abbott et al. 2023a). In Sect. 6, we draw our conclusions and discuss prospects for future development.

2 Inhomogeneous Poisson process and Bayesian inference

The main application of ICAROGW is to infer the population parameters, A, that describe the production rate of events in terms of the GW parameters, θ, namely, dN dt dθ(Λ)${{{\rm{d}}N} \over {{\rm{d}}t{\rm{d}}\theta }}(\Lambda )$. For instance, in Sect. 4 we will consider the rate of CBC in terms of the source frame masses. The hierarchical likelihood of obtaining Nobs observations, each described by some parameters, θ, in a data collection {x} for a given observing time, Tobs, from a population of events, with a constant rate and in presence of selection biases is given by (see Mandel et al. 2019; Vitale et al. 2022 for a detailed derivation): ({x}Λ)eNexp(Λ)i=1Nobs (xiθ,Λ)dN dtdθ(Λ)dtdθ        eNexp(Λ)i=1Nobs Tobs (xiθ,Λ)dNdtdθ(Λ)dθ.$\eqalign{ & {\cal L}(\{ x\} \mid \Lambda ) \propto {e^{ - {N_{\exp }}(\Lambda )}}\prod\limits_{i = 1}^{{N_{{\rm{obs }}}}} {\int {\cal L} } \left( {{x_i}\mid \theta ,\Lambda } \right){{{\rm{d}}N} \over {{\rm{d}}t{\rm{d}}\theta }}(\Lambda ){\rm{d}}t{\rm{d}}\theta \cr & & \,\,\,\,\,\,\, \propto {e^{ - {N_{\exp }}(\Lambda )}}\prod\limits_{i = 1}^{{N_{{\rm{obs }}}}} {{T_{{\rm{obs }}}}} \int {\cal L} \left( {{x_i}\mid \theta ,\Lambda } \right){{{\rm{d}}N} \over {{\rm{d}}t{\rm{d}}\theta }}(\Lambda ){\rm{d}}\theta . \cr} $(1)

Equation (1) is often referred to as a “hierarchical likelihood”. By assuming a “scale-free” (i.e., neglecting information about the rate) prior π(Nexp) ∝ 1/Nexp on the expected number of detections, an equivalent form of Eq. (1) can be derived: (xΛ)i=1Nobs (xiθ,Λ)dN dtdθdθpdet(θ,Λ)dN dtdθdθ.${\cal L}(x\mid \Lambda ) \propto \prod\limits_{i = 1}^{{N_{{\rm{obs }}}}} {{{\int {\cal L} \left( {{x_i}\mid \theta ,\Lambda } \right){{{\rm{d}}N} \over {{\rm{d}}t{\rm{d}}\theta }}{\rm{d}}\theta } \over {\int {{p_{\det }}} (\theta ,\Lambda ){{{\rm{d}}N} \over {{\rm{d}}t{\rm{d}}\theta }}{\rm{d}}\theta }}} $(2)

Both the hierarchical likelihoods in Eqs. (1) and (2) contain several crucial quantities to the inference problem. In what follows, we will explain how ICAROGW numerically computes these quantities and obtains a numerical evaluation of the full hierarchical likelihood defined in Eq. (1).

The first central quantity is the single-event likelihood ℒ(xiθ, Λ), it quantifies the measurement uncertainties of the GW parameters θ. The single-event likelihood is usually not directly accessible, instead we are provided with Ns,i (PE) posterior samples drawn from p(θ|xi, Λ) ∝ ℒ(xi|θ, Λ)πPE(θ|Λ), where πPE(θ|Λ) is the prior used to generate the samples. The PE samples represent the source parameters that we believe given the observation of data x. ICAROGW numerically evaluates the likelihood integral in Eq. (1) and in the numerator of Eq. (2) via Monte Carlo integration by summing over PE samples: (xiθ,Λ)dN dtdθ(Λ)dθ 1Ns,ij=1Ns,i1πPE(θi,jΛ)dNdtdθ(Λ)|i,j1Ns,ij=1Ns,iwi,j,$\eqalign{ & \int {\cal L} \left( {{x_i}\mid \theta ,\Lambda } \right){{{\rm{d}}N} \over {{\rm{d}}t{\rm{d}}\theta }}(\Lambda ){\rm{d}}\theta \approx \cr & & & {\left. {{1 \over {{N_{{\rm{s}},i}}}}\sum\limits_{j = 1}^{{N_{{\rm{s}},i}}} {{1 \over {{\pi _{{\rm{PE}}}}\left( {{\theta _{i,j}}\mid \Lambda } \right)}}} {{dN} \over {dtd\theta }}(\Lambda )} \right|_{i,j}} \equiv {1 \over {{N_{{\rm{s}},i}}}}\sum\limits_{j = 1}^{{N_{{\rm{s}},i}}} {{w_{i,j}}} , \cr} $(3)

where the index, i, refers to the event and the index, j, to the posterior samples of the events. The weight wi,j is the number of CBC mergers happening per unit of time. As Eq. (3) is evaluated with a finite sum over posterior samples, we introduce the numerical stability estimator called “effective number of posterior samples” per event i, as introduced in Talbot & Golomb (2023): Neff,i=(jNs,iwi,j)2jNs,iwi,j2.${N_{{\rm{eff}},i}} = {{{{\left( {\sum\limits_j^{{N_{{\rm{s}},i}}} {{w_{i,j}}} } \right)}^2}} \over {\sum\limits_j^{{N_{{\rm{s}},i}}} {w_{i,j}^2} }}.$(4)

This estimator quantifies how many samples per event are contributing to the evaluation of the integral. Typically, in population analyses such as the one presented in Abbott et al. (2023b), it is necessary to have at least an effective number of posterior samples equal to 20 for each event and population model supported by the analysis. In case this requirement is not satisfied, ICAROGW will artificially associate a null likelihood to the population model, as the model cannot be trusted.

The second central quantity is the “expected number of events”, Nexp(Λ), which is related to the selection bias and can be evaluated as: Nexp (Λ)=Tobs pdet (θ,Λ)dN dtdθdθ,${N_{{\rm{exp }}}}(\Lambda ) = {T_{{\rm{obs }}}}\int {{p_{{\rm{det }}}}} (\theta ,\Lambda ){{{\rm{d}}N} \over {{\rm{d}}t{\rm{d}}\theta }}{\rm{d}}\theta ,$(5)

where pdet(θ, Λ) is a detection probability that can be calculated as: pdet(θ,Λ)=xdetectable(xiθ,Λ)dx${p_{{\rm{det }}}}(\theta ,\Lambda ) = \int_{x \in {\rm{ detectable }}} {\cal L} \left( {{x_i}\mid \theta ,\Lambda } \right){\rm{d}}x$(6)

However, we do not have access to an analytical form of the detection probability, unless some simplifying assumptions are made; we refer to Gair et al. (2023) for an introductory example in the context of GW cosmology with galaxy catalogues. The current approach to evaluate selection biases is to use Monte Carlo simulations of injected and detected events (Abbott et al. 2021b), often shortly referred to as “injections”. The injections are used to evaluate the signal detectable volume that can be explored in the parameter space and correct for selection biases. The occurrence of detected injections, proportional to pdet(θ, Λ), and the population model used to generate them can be used to evaluate the selection bias. ICAROGW takes as it input a set of Ndet detected signals out of Ngen total injections generated from

a prior πinj (θ) to calculate the integral in Eq. (5), using Monte Carlo integration as follows: NexpTobsNgenj=1Ndet1πinj(θj)dNdtdθ|jTobsNgenj=1Ndetsj.${\left. {{N_{{\rm{exp }}}} \approx {{{T_{{\rm{obs }}}}} \over {{N_{{\rm{gen }}}}}}\sum\limits_{j = 1}^{{N_{{\rm{det }}}}} {{1 \over {{\pi _{{\rm{inj }}}}\left( {{\theta _j}} \right)}}} {{{\rm{d}}N} \over {{\rm{d}}t{\rm{d}}\theta }}} \right|_j} \equiv {{{T_{{\rm{obs }}}}} \over {{N_{{\rm{gen }}}}}}\sum\limits_{j = 1}^{{N_{{\rm{det }}}}} {{s_j}} $(7)

Here, we have again defined a weight, Sj, with the dimension of a CBC merger rate per detector time. The injection prior πinj(θ) must be properly normalized to obtain a reasonable value of Nexp, while a wrong normalization of πPE(θ) (used in Eq. (3)) will only result in an overall normalization factor to the overall hierarchical likelihood. Following Farr (2019), for Eq. (7), we can also define a numerical stability estimator, namely, the “effective number of injections”, as follows: Neff, inj =[ jNdet sj ]2[ jNdet sj2Ngen 1(jNdet sj)2 ],${N_{{\rm{eff}},{\rm{ inj }}}} = {{{{\left[ {\sum\limits_j^{{N_{{\rm{det }}}}} {{s_j}} } \right]}^2}} \over {\left[ {\sum\limits_j^{{N_{{\rm{det }}}}} {s_j^2} - N_{{\rm{gen }}}^{ - 1}{{\left( {\sum\limits_j^{{N_{{\rm{det }}}}} {{s_j}} } \right)}^2}} \right]}},$(8)

where typical value for numerical stability is Neff,inj > 4Nobs.

In summary, the hierarchical likelihoods defined in Eqs. (1) and (2) are calculated using the approximations of the Monte Carlo integrals in Eqs. (3) and (7). The Monte Carlo integrated version of Eq. (1) is then: ln[({x}Λ)]TobsNgenj=1Noctsj+iNobsln[ TobsNs,ij=1Ns,iwi,j ],$\ln [{\cal L}(\{ x\} \mid \Lambda )] \approx - {{{T_{{\rm{obs}}}}} \over {{N_{{\rm{gen}}}}}}\sum\limits_{j = 1}^{{N_{{\rm{oct}}}}} {{s_j}} + \sum\limits_i^{{N_{{\rm{obs}}}}} {\ln } \left[ {{{{T_{{\rm{obs}}}}} \over {{N_{{\rm{s}},i}}}}\sum\limits_{j = 1}^{{N_{{\rm{s}},i}}} {{w_{i,j}}} } \right]$(9)

while for the scale-free version, the Monte Carlo version for Eq. (2) is: ln[({x}Λ)]Nobsln[ 1Ngenj=1Ndet sj ]+iNobsln[ 1Ns,ij=1Ns,iwi,j].$\ln [{\cal L}(\{ x\} \mid \Lambda )] \approx - {N_{{\rm{obs}}}}\ln \left[ {{1 \over {{N_{{\rm{gen}}}}}}\sum\limits_{j = 1}^{{N_{{\rm{det }}}}} {{s_j}} } \right] + \sum\limits_i^{{N_{{\rm{obs}}}}} {\ln } \left[ {{1 \over {{N_{{\rm{s}},i}}}}\sum\limits_{j = 1}^{{N_{{\rm{s}},i}}} {{w_{i,j}}} } \right]$(10)

For each population model, ICAROGW calculates one of the two likelihoods and the numerical stability estimators. If at least one of these numerical estimators is below the threshold set by the user, ICAROGW returns a result of ln[ℒ({x}|Λ)] = −∞. This keeps the population model from being chosen.

thumbnail Fig. 1

ICAROGW logical structure. The orange boxes highlight the computation of the rate weights which are needed for the calculation of the hierarchical likelihood given in Eq. (9).

3 Code structure and GPU performance

ICAROGW contains several python modules used for population inference. in Fig. 1, we display a schematic view of ICAROGW. The PE samples and injections are stored in two different objects, although they interact with common functions to calculate the CBC merger rate. The CBC rate models are the central objects of the code, they are python classes that contain the production rate dNdtdθ${{{\rm{d}}N} \over {{\rm{d}}t{\rm{d}}\theta }}$. Each one specifies what are the event-level parameters θ used to calculate dNdtdθ${{{\rm{d}}N} \over {{\rm{d}}t{\rm{d}}\theta }}$ and what are the population-level parameters, Λ, used to calculate the rate model. The rate class also contains instructions on how to update the rate model from the population parameters Λ. As we describe more in detail in Sect. 4, for GW cosmology, ICAROGW distinguishes between detector frame frame parameters, θD, and source frame parameters, θS. Detector frame parameters are the CBC physical parameters observed on Earth, while the source frame ones are the actual parameters of the source. The Gw parameters θD that are necessary to calculate the CBC merger rate are passed to a CBC rate class that performs the conversion to source frame quantities, θS, and calculates the weights, wi,j and Sj, that are to be used for the calculation of the hierarchical likelihood and numerical stability estimators. Finally, the hierarchical likelihood computation is handled by a function embedded in the python package BILBY (Ashton et al. 2019b,a; Ashton & Talbot 2021) for Bayesian inference.

Almost all the functions of ICAROGW can also be executed on Graphical Processing Units (GPUs). A GPU is a specialized processor devoted to light and parallel computations. ICAROGW uses cupy (Okuta et al. 2017) to parallelize on GPU the computation of the weights, wi,j, sk, and all the quantities that are computed from them.

The use of the GPU greatly enhances the computation time needed for the calculation of the hierarchical likelihood with respect to the time taken by Central Processing Units (CPUs). in Fig. 2, we show a comparison between the timing of the hierarchical likelihood in Eq. (1) computation with CPU and GPU using double float precision. The figure reports two figures of merit, the total time for the likelihood computation and the time taken per sample, in terms of the number of samples (PE+injections) used in the code. The figures of merit have been generated on a computing machine employed to obtain the results of Mastrogiovanni et al. (2023).

The total time of the likelihood is a proxy for the time required to obtain a population fit as a function of the data size.

From Fig. 2, we can observe three regimes. For a number of samples (<104), the CPU’s and GPU’s timing coincide, meaning that the computing time is not dominated by the sample size but by the internal python routines. For a number of samples between 104 and 106, the CPU’s time starts to be dominated by the sample size. In this middle regime, the GPU timing can be further improved by acting on routines that do not use data samples. The last regime is achieved for a number of samples higher than 106, where both the CPU and GPU timings are dominated by the data samples. In agreement with Talbot et al. (2019), we find that the GPU enhance the computation time of a factor ~100 for population inference.

In Mastrogiovanni et al. (2023), about 105 samples were used on a CPU to obtain a population fit in 1 week of computation time. If 106 samples on CPU are used (i.e., ten times more Gw events are considered), then the total time required would be equal to as many as ten weeks of computation. Indeed, ten weeks of computation would not be strictly computationally prohibitive, but this choice could certainly impact the number of models and exploratory studies that could be performed.

The efficiency per sample is another way of quantifying when the CPU’s and GPU’s timings will start scaling with the number of samples. The efficiency per sample is defined as the data sample size divided by the total computing time. Figure 2 shows three regimes analogous to the total likelihood time. For a number of samples (<104), the efficiency improves for both the CPU and GPU and the computing time is not dominated by the sample size. For a number of samples higher than 104, the CPU has already saturated its efficiency to a value of 106 samples per second and the CPU’s total time will start scaling linearly with the number of samples. Finally, for a total number of samples higher than 106, the GPU’s efficiency has saturated to a value of 108 samples per second per sample and the total likelihood time will start scaling linearly as the samples pass.

Indeed, we expect that different hardware configurations will display distinct values for the total likelihood time, efficiency, and for the transition between the three regimes. However, we expect the considerations discussed above to apply in general to all hardware facilities including GPU and CPU. Therefore, GPU computing will play a crucial role for future populations of GW events.

thumbnail Fig. 2

ICAROGW computational benchmarking on CPU and GPU. Top: average elapsed time taken by ICAROGW to evaluate the hierarchical likelihood in Eq. (1). Bottom: efficiency, defined as number of samples processed per second to evaluate the hierarchical likelihood in Eq. (1). For both panels, the CPU is reported with blue pentagons and GPU with orange squares. The two panels are reported as a function of a number of samples, namely, the sum of the number of PE samples and injections used for the computation. The CPU markers are generated up to 106 samples since at 105 samples the efficiency is already saturated. Figures generated with double float precision on CPU Intel Core i9-11950H (8 cores HT, 2.6–5.0 GHz Turbo) and GPU NVIDIA GeForce RTX3080 (16Gb GDDR6, 6144 cores CUDA).

4 Application to a compact binary coalescence case

Although ICAROGW can be adapted to any problem involving an inhomogeneous Poisson process in the presence of selection bias, the code is mostly developed for GW science with CBCs. From the Gw detections, we estimated the luminosity distance of the source, dL, the detector masses, m1,d, m2,d, and the sky position, Ω, as well as a number of parameters related to the spins of the compact objects (we indicate them generally with χ, although two parameterizations are available; see Appendix B.3). Therefore, the rate of CBC mergers as seen in detectors' data must be modelled as a function of dL, m1,d, m2,d, and χ, namely with a rate: dN ddLdΩdm1,d dm2,d dχdtd,${{{\rm{d}}N} \over {{\rm{d}}{d_{\rm{L}}}{\rm{d}}\Omega {\rm{d}}{m_{1,d}}{\rm{d}}{m_{2,d}}{\rm{d}}\chi {\rm{d}}{t_{\rm{d}}}}},$(11)

where td is the detector unit time. Equation (11) is the CBC detector rate model. The rate in Eq. (11) does not strictly include any information about the Universe expansion model and its parameters. The information about cosmology can be included by rewriting the CBC merger rate in terms of source-frame parameters such as the source masses and redshift. The source masses are related to the detector masses through the relation: m1/2,d=m1/2,s(1+z),${m_{1/2,d}} = {m_{1/2,s}}(1 + z),$(12)

and the luminosity distance is related to the redshift by the choice of a cosmological model and possible GR deviations at cosmo-logical scales. ICARoGw contains several parameterizations for the luminosity distance that are reported explicitly with their “population parameters” in Appendix A.1 for GR models and in Appendix A.2 for modified gravity models.

The CBC merger rate at the detector can be written as a function of source-frame parameters by performing a change of variables that involves the choice of the cosmological parameters, namely: dN dθDdtd=dN dθSdtsdtsdtd1detJDS=dN dθSdts11+z1detJDS.${{{\rm{d}}N} \over {{\rm{d}}{\theta _{\rm{D}}}{\rm{d}}{t_{\rm{d}}}}} = {{{\rm{d}}N} \over {{\rm{d}}{\theta _{\rm{S}}}{\rm{d}}{t_{\rm{s}}}}}{{{\rm{d}}{t_{\rm{s}}}} \over {{\rm{d}}{t_{\rm{d}}}}}{1 \over {\det {J_{{\rm{D}} \to {\rm{S}}}}}} = {{{\rm{d}}N} \over {{\rm{d}}{\theta _{\rm{S}}}{\rm{d}}{t_{\rm{s}}}}}{1 \over {1 + z}}{1 \over {\det {J_{{\rm{D}} \to {\rm{S}}}}}}.$(13)

In the equation above, the factor of 1/1 + z comes from the difference between source-frame and detector-frame time, and detJDS=dLz(1+z)2$\det {J_{{\rm{D}} \to {\rm{S}}}} = {{\partial {d_{\rm{L}}}} \over {\partial z}}{(1 + z)^2}$(14)

is the Jacobian for the change of variables between the detector and the source frames. The expression of the differential of the luminosity distance can be found in Appendix A.1 for standard cosmological models and in Appendix A.2 for modified gravity models. The term dN dθSdts${{{\rm{d}}N} \over {{\rm{d}}{\theta _{\rm{S}}}{\rm{d}}{t_{\rm{s}}}}}$

is the CBC merger rate in the source frame and it is described by the distributions of CBCs in terms of redshift, spins and source frame masses. Several parameterizations of the CBC merger rate in the source frame are possible depending on the observational data available. We cover all these cases in the subsections below.

4.1 Spectral sirens merger rates

The first case that we discuss is the “spectral siren” analysis (Ezquiaga & Holz 2022). For this case, we are interested in inferring population properties of the source rate model, as well as cosmology and GR deviations from a population of GW events alone. For this model, the detector event parameters θD are (dL, m1,d, m2,d,χ), namely, luminosity distance, detector masses, and spin parameters. The source event parameters are θs = (z, m1,s, m2,s,χ), namely, the redshift, two source masses, and the spin parameters. The detector rate for the spectral sirens analysis is parameterized as: dNddLdm1,ddm2,ddχdtd=R0Ψ(z;Λ)×ppop(m1,s,m2,sΛ)ppop(χΛ)dVcdz11+z1detJDS,$\eqalign{ & {{{\rm{d}}N} \over {{\rm{d}}{d_{\rm{L}}}{\rm{d}}{m_{1,d}}{\rm{d}}{m_{2,d}}{\rm{d}}\chi {\rm{d}}{t_{\rm{d}}}}} = {R_0}\Psi (z;\Lambda ) \cr & \times {p_{{\rm{pop }}}}\left( {{m_{1,s}},{m_{2,s}}\mid \Lambda } \right){p_{{\rm{pop }}}}(\chi \mid \Lambda ){{{\rm{d}}{V_{\rm{c}}}} \over {{\rm{d}}z}}{1 \over {1 + z}}{1 \over {\det {J_{{\rm{D}} \to {\rm{S}}}}}}, \cr} $(15)

where R0 is the CBC merger rate per comoving volume per year (in Gpc−3 yr−1), Ψ(z; Λ) is a function parametrizing the rate evolution in redshift such that Ψ(z=0; Λ) = 1, ppop(m1,s, m2,s∣Λ) is a prior distribution describing the production of source masses, and ppop(χ∣Λ) is a prior distribution for the production of spin parameters. ICAROGW implements the prior and rate distributions that are typically used for CBC population studies and motivated by population syntheses studies of CBCs (Abbott et al. 2023b). Appendix B.2 lists and describes the models implemented for the masses, Appendix B.3 the models for spins and Appendix B.1 the models for the merger rates as a function of redshift. Finally, Vc is the comoving volume, that also depends on the cosmological parameters.

4.2 Galaxy catalog merger rates

The “galaxy catalog” parametrization adds information on the GW event redshift from galaxy surveys (Schutz 1986; Del Pozzo 2012; Gray et al. 2020, 2022). Also in this case, we are interested in inferring population properties of the source rate model, as well as the cosmology and GR deviations from a population of GW events, this time with extra information coming from the galaxy catalog.

In this method, the detector parameters θD = (dL, m1,d, m2,d, Ω,χ) are the luminosity distance and detector masses, sky direction pixel, and spins. The sky direction pixel area is measured in squared radians. The source event parameters are θS = (z, m1,s, m2,s, Ω χ), that is, the redshift and the two source masses, sky direction, and spins. Here, we directly provide the parametrization of the CBC merger rate for this method, for a more detailed derivation, see Mastrogiovanni et al. (2023). The detector rate for the galaxy catalog analysis is parameterized as: dN ddLdm1,d dm2,d dΩdχdtd=Rgal,0*Ψ(z;Λ)1+zppop(m1,s,m2,sΛ)detJDSppop(Λ) [ dVcdzdΩϕ*(H0)Γinc(α+ϵ+1,xmax(Mthr),xmin) +j=1Ngal(Ω)fL(M(mj,z);ϵ)p(zzobsj,σz,obsj)ΔΩ ]$\eqalign{ & {{{\rm{d}}N} \over {{\rm{d}}{d_{\rm{L}}}{\rm{d}}{m_{1,d}}{\rm{d}}{m_{2,d}}{\rm{d}}\Omega {\rm{d}}\chi {\rm{d}}{t_{\rm{d}}}}} = R_{{\rm{gal}},0}^*{{\Psi (z;\Lambda )} \over {1 + z}}{{{p_{{\rm{pop}}}}\left( {{m_{1,s}},{m_{2,s}}\mid \Lambda } \right)} \over {\det {J_{{\rm{D}} \to {\rm{S}}}}}} \cr & {p_{{\rm{pop}}}}(\mid \Lambda )\left[ {{{d{V_{\rm{c}}}} \over {dzd\Omega }}{\phi _*}\left( {{H_0}} \right){\Gamma _{{\rm{inc}}}}\left( {\alpha + + 1,{x_{\max }}\left( {{M_{{\rm{thr}}}}} \right),{x_{\min }}} \right)} \right. \cr & \left. { + \sum\limits_{j = 1}^{{N_{{\rm{gal}}}}(\Omega )} {{{{f_L}\left( {M\left( {{m_j},z} \right);} \right)p\left( {z\mid z_{{\rm{obs}}}^j,\sigma _{{\rm{z}},{\rm{obs}}}^j} \right)} \over {\Delta \Omega }}} } \right] \cr} $(16)

where Rgal,0*$R_{{\rm{gal}},0}^*$ is the local CBC merger rate per galaxy per year (in yr−1). The sum of the two terms in the square brackets represents the galaxy number density in redshift and sky area that could host GW sources (see Mastrogiovanni et al. 2023 for more details). The first term is the completeness correction, that is, it accounts for missing galaxies. It depends on the absolute magnitude threshold of galaxy detection, Mthr, on how likely more luminous galaxies can emit GW events (through the e parameter), and on the Schecter luminosity function and its parameters ϕ*, and α, xmin/max are defined in Mastrogiovanni et al. (2023) and are related to the minimum and maximum of the Schecter function. The second term in the square brackets is given by the galaxy distribution reported in the catalog. The function fL(M(mj,z); ) quantifies how likely luminous galaxies emit GW events, while p(zzobsJ,σz,obsJ)$p\left( {z\mid z_{{\rm{obs}}}^J,\sigma _{{\rm{z}},{\rm{obs}}}^J} \right)$ is the probability of having a certain value of z given observed values of galaxy redshift inside the catalog.

4.3 Sources with observed EM counterparts

A third methodology to infer population and cosmological properties for GW events with an associated EM counterpart is also implemented. The model now takes into account additional constraints on sky position and redshift from a putative EM counterpart. In this case, the CBC merger rate can still be parametrized as done in the previous sections, but the overall hierarchical likelihood should be modified to include the additional EM information. We might be in two types of observational cases that we describe in the following.

4.3.1 High-latency analyses

In the case that we perform a “high latency” population analysis for the GW and EM observation, we expect GW studies to have produced PE samples on all the GW parameters such as sky position, redshift, and masses. In addition, the EM counterpart will provide us with an accurate sky localization of the source and a redshift localization with small uncertainties due to the peculiar velocities of the host galaxy.

By assuming that the GW measure is independent of the EM measure, the overall likelihood term is now ℒEM+GW(xi|z, Ω, m1,s, m2,s), which describes the measure of z, Ω, m1,s, m2,s from EM and GW data. Here we assume that the EM data measures Ω and z, while the GW data can measure Ω, z, m1,s, m2,s independently, so that: EM+GW(xiz,Ω,m1,s,m2,s,χ)EM(xiz,Ω)GW(xiz,Ω,m1,s,m2,s,χ).$\eqalign{ & {{\cal L}_{{\rm{EM}} + {\rm{GW}}}}\left( {{x_i}\mid z,\Omega ,{m_{1,s}},{m_{2,s}},\chi } \right) \propto \cr & \quad {{\cal L}_{{\rm{EM}}}}\left( {{x_i}\mid z,\Omega } \right){{\cal L}_{{\rm{GW}}}}\left( {{x_i}\mid z,\Omega ,{m_{1,s}},{m_{2,s}},\chi } \right). \cr} $(17)

The integral of the numerator in Eq. (2) becomes I=EM(xiz,Ω)GW(xiz,Ω,m1,s,m2,s,χ)      dN dz dΩdm1,s dm2,s dχdts11+z dm1,s,m2,s dχdz .$\eqalign{ & I = \int {{{\cal L}_{{\rm{EM}}}}} \left( {{x_i}\mid z,\Omega } \right){{\cal L}_{{\rm{GW}}}}\left( {{x_i}\mid z,\Omega ,{m_{1,s}},{m_{2,s}},\chi } \right) \cr & \,\,\,\,\,\,{{{\rm{d}}N} \over {{\rm{d}}z{\rm{d}}\Omega {\rm{d}}{m_{1,s}}{\rm{d}}{m_{2,s}}{\rm{d}}\chi {\rm{d}}{t_{\rm{s}}}}}{1 \over {1 + z}}{\rm{d}}{m_{1,s}},{m_{2,s}}{\rm{d}}\chi {\rm{d}}z{\rm{d}}\Omega . \cr} $(18)

Regarding the denominator of Eq. (2), namely the selection bias, following the implementation of the same methodology in Gray et al. (2020), we assume that the probability of detecting an EM counterpart (if present) given the GW detection is always 1. This assumption allows us to calculate the selection bias in Eq. (2) only using GW detectable signals, without the need to fold in a detection probability model for the EM counterpart. This assumption is physically valid only when the detection range of GW observatories is significantly lower than the detection range of EM observatories. As GW detectors improve their sensitivity this assumption should be revised as it could introduce a systematic bias on the estimation of H0 (Chen 2020).

The integral in Eq. (18) can not be directly evaluated by calculating CBC rate weights wi,j and slj as described in Sect. 2. Some additional steps need to be taken to account for the information from the EM counterpart. To perform the integral in Eq. (18), ICAROGW firstly identifies all the GW PE samples falling in the sky pixel in which we believe the EM counterpart is located. The PE selection procedure practically corresponds to performing the integral on dΩ in Eq. (18). Next, a function of redshift F(z) is defined from the selected GW PE samples. The function is F(z)=[ 1NsEMiNsEMwi ]KDE[ zi, weights =wi ],$F(z) = \left[ {{1 \over {N_s^{{\rm{EM}}}}}\sum\limits_i^{N_s^{{\rm{EM}}}} {{w_i}} } \right]{\rm{KDE}}\left[ {{z_i},{\rm{ weights }} = {w_i}} \right],$(19)

where KDE is a kernel density estimate performed on the GW redshift samples zi with weights wi. The weights are given by: wi=1πEM(zi)πPE(zi,mi,χi)dNdzdm1,sdm2,sdχdts|i11+zi,${w_i} = {\left. {{1 \over {{\pi _{{\rm{EM}}}}\left( {{z^i}} \right){\pi _{{\rm{PE}}}}\left( {{z^i},{{\bf{m}}^i},{{\bf{\chi }}^i}} \right)}}{{{\rm{d}}N} \over {{\rm{d}}z{\rm{d}}{m_{1,s}}{\rm{d}}{m_{2,s}}{\rm{d}}\chi {\rm{d}}{t_{\rm{s}}}}}} \right|_i}{1 \over {1 + {z^i}}},$(20)

where πEM(z) is the prior used by the EM experiment to provide a redshift measure. The definition of F(z), and the fact that we selected GW PE in the EM counterpart localization area, allows us to rewrite Eq. (18) as IpEM(zxi,Ω)F(z)dz,$I \propto \int {{p_{{\rm{EM}}}}} \left( {z\mid {x_i},\Omega } \right)F(z){\rm{d}}z$(21)

where pEM(z|xi, Ω) is EM posterior that localizes the EM counterpart in redshift given a certain sky area. The above integral can be evaluated with the usual Monte Carlo sum, by assuming a set of PE samples on z drawn from the EM posterior as: I1Ns,EMiNsEMF(zi).$I \approx {1 \over {{N_{{\rm{s}},{\rm{EM}}}}}}\sum\limits_i^{{N_{{\rm{sEM}}}}} F \left( {{z_i}} \right).$(22)

4.3.2 Low-latency analyses

In some cases, we might want to produce a “low latency” inference of population parameters (e.g. Hubble constant) given a GW and EM observation. in this case, we are typically provided with low-latency EM and GW low-latency skymap. The low-latency GW skymap can be thought of as a “posterior” PGW(dL, Ω|xi) of the luminosity distance and sky position given the observed GW event. The low-latency EM skymap can be thought of as a posterior pEM(z, Ω|xi) for redshift and sky position. Also for this case, we make the assumption that the selection bias is dominated by the GW detection process and therefore no modelling for the EM selection bias is needed.

The GW-EM likelihood is a function only of sky position and redshift and can be written as EM+GW(xiz,Ω)EM(xiz,Ω)GW(xiz,Ω).${{\cal L}_{{\rm{EM}} + {\rm{GW}}}}\left( {{x_i}\mid z,\Omega } \right) \propto {{\cal L}_{{\rm{EM}}}}\left( {{x_i}\mid z,\Omega } \right){{\cal L}_{{\rm{GW}}}}\left( {{x_i}\mid z,\Omega } \right).$(23)

The integral in the numerator of the hierarchical Likelihood in Eq. (2) becomes I=EM(xiz,Ω)GW(xiz,Ω)dNdzdΩdts11+zdz,$I = \int {{{\cal L}_{{\rm{EM}}}}} \left( {{x_i}\mid z,\Omega } \right){{\cal L}_{{\rm{GW}}}}\left( {{x_i}\mid z,\Omega } \right){{{\rm{d}}N} \over {{\rm{d}}z{\rm{d}}\Omega {\rm{d}}{t_{\rm{s}}}}}{1 \over {1 + z}}{\rm{d}}z{\rm{d}}\Omega ,$(24)

and the CBC merger rate is parameterized only in terms of red-shift and sky position. We can now use the Bayes theorem to rewrite the above equation as: IpEM(z,Ωxi)πEM(z,Ω)pGW(dL(z,Λ),Ωxi)πGW(dL(z,Λ),Ω)dN dz dΩdts11+z dz ,$I \propto \int {{{{p_{{\rm{EM}}}}\left( {z,\Omega \mid {x_i}} \right)} \over {{\pi _{{\rm{EM}}}}(z,\Omega )}}} {{{p_{{\rm{GW}}}}\left( {{d_{\rm{L}}}(z,\Lambda ),\Omega \mid {x_i}} \right)} \over {{\pi _{{\rm{GW}}}}\left( {{d_{\rm{L}}}(z,\Lambda ),\Omega } \right)}}{{{\rm{d}}N} \over {{\rm{d}}z{\rm{d}}\Omega {\rm{d}}{t_{\rm{s}}}}}{1 \over {1 + z}}{\rm{d}}z{\rm{d}}\Omega ,$(25)

where with Λ we indicate any population parameter required to calculate the luminosity distance from redshift, while with πPE and πGW the priors used to generate the low-latency EM and GW skymaps. The integral now explicitly depends on the EM and GW low-latency skymaps. Eq. (25) can be evaluated analytically by drawing NEM samples on z and Ω the EM skymap, namely, I1NEMj=0NEM1πEM(zj,Ωj)pGW(dL(zj,Λ),Ωjxi)πGW(dL(zj,Λ),Ωj)dN dz dΩdts|j11+zj.${\left. {I \approx {1 \over {{N_{{\rm{EM}}}}}}\sum\limits_{j = 0}^{{N_{{\rm{EM}}}}} {{1 \over {{\pi _{{\rm{EM}}}}\left( {{z_j},{\Omega _j}} \right)}}} {{{p_{{\rm{GW}}}}\left( {{d_{\rm{L}}}\left( {{z_j},\Lambda } \right),{\Omega _j}\mid {x_i}} \right)} \over {{\pi _{{\rm{GW}}}}\left( {{d_{\rm{L}}}\left( {{z_j},\Lambda } \right),{\Omega _j}} \right)}}{{{\rm{d}}N} \over {{\rm{d}}z{\rm{d}}\Omega {\rm{d}}{t_{\rm{s}}}}}} \right|_j}{1 \over {1 + {z_j}}}.$(26)

The motivation for which the integral is evaluated by drawing samples from the EM skymap is that typically the EM source is better localized than the GW source. By drawing samples from the EM skymaps, we have a higher probability for the weights in Eq. (26) to be different from zero. Consequently, the numerical evaluation of I is more stable.

5 Code validation

We performed several tests reproducing population results on CBC that are consistent with previous analyses. in all the cases, we found that ICAROGW 2.0 is able to reproduce the results of previous works in the literature.

5.1 Spectral siren analyses

We used the same 42 BBHs used in Abbott et al. (2023b) to reproduce populations, cosmological, and beyond-GR analyses generated with ICAROGW 1.0. For the benchmarking, we employed a POWERLAW + PEAK model (see Appendix B.2) for the source mass distribution of BBHs and a MADAU-DICKISON merger rate model (see Appendix B.1). To evaluate the selection biases, we used the procedure described in Sect. 2, based on a set of detectable GW signals released with Abbott et al. (2023b).

The first analysis that we reproduced is the one presented in Abbott et al. (2023b), which conjointly estimates the population parameters for the BBH distribution together with the cosmo-logical parameters H0, Ωm, and w0 (see Appendix A.1 for more details). The first column of Fig. 3 shows the quantile-quantile plots (QQ-plots) for the marginal posteriors on the population parameters obtained with ICAROGW 1.0 and ICAROGW 2.0. When the two posterior are in agreement, the QQ-plot displays a bisector. For all the population parameters, we obtained an excellent agreement between the two posteriors, as shown in the first column of Fig. 3.

We also reproduced the constraints on modified gravity models and the population of BBHs generated with ICAROGW 1.0 in Leyde et al. (2022). The second, third, and fourth columns of Fig. 3 show the QQ-plots for three modified gravity models (see Sect. A.2 for more details). Also, in this case, the posteriors generated with ICAROGW 1.0 and ICAROGW 2.0 are in a good match between each other. The tests with spectral sirens are also in good agreement with results obtained from the independent code MGcosmoPop (Mancarella & Genoud-Prachex 2022) in Mancarella et al. (2022) (for the same set of BBH events). They are also in agreement with the results of Ezquiaga (2021) using events from O3a modified gravity propagation and with the population-only results generated with the code gwpopulation (Talbot et al. 2019).

thumbnail Fig. 3

Quantile-quantile plots of the marginal posterior distributions of the CBC merger rate population parameters inferred from the 42 BBHs in GWTC-3. The plots columns correspond to the model used to calculate the GW luminosity distance as a function of redshift, see Appendices A.1A.2. The legend and different lines report the population parameters. All the benchmark runs use the POWER LAW + PEAK model for the source mass (cf. Appendix B.2) and the MADAU-DICKINSON merger rate model (cf. Appendix B.1).

5.2 Galaxy catalog analysis

We tested ICAROGW against the results generated by gwcosmo (Gray et al. 2020, 2022) in Abbott et al. (2023b) using the GLADE+ (Dalya et al. 2022) galaxy catalog with the infrared K-band.

The procedure to construct a redshift distribution of possible host galaxies from the catalog is described in detail in Mastrogiovanni et al. (2023). In summary, all the galaxies reported in GLADE+ are divided into equal-sized sky pixels large 53 deg2. For each pixel, an apparent magnitude threshold is defined as the median of the apparent magnitude reported for all the galaxies. Galaxies fainter than the apparent magnitude threshold are removed from the catalog. Then, galaxies are subdivided into smaller pixels of 13 deg2 that are the ones used to run the analysis. The completeness correction is applied following the procedure described in Mastrogiovanni et al. (2023), which consists of assuming a Schecter function model for the galaxy distribution and counting all the galaxies missing due to the apparent magnitude threshold cut. The Schechter function model assumed for the K-band galaxies is specified by the parameters Mmin = −27.85, Mmax = −19.84, α = −1.09, ϕ* = 0.03 Mpc−3 (for H0 = 67.7 hu). For this benchmark analysis, following Abbott et al. (2023b), we assume that for each galaxy reported in the catalog, the probability of hosting a GW event is proportional to the galaxy’s intrinsic luminosity. The source mass and redshift population model is fixed to the same used by gwcosmo in Abbott et al. (2023b).

We only considered the inference for the Hubble constant, H0, to make a direct comparison with Abbott et al. (2023b). The results that we obtained are shown in Fig. 4. We found a good agreement for almost all the 42 BBHs present in the dataset. This is not surprising since most of the results are dominated by the assumption of the source mass distribution (Abbott et al. 2023b) and, as discussed in the previous section, ICARoGw 2.0 agrees well with previous literature when using the spectral siren method. However, we note that the most close-by events such as GW150914, GW170814, and GW190824_021846 are only partially dominated by population assumption on masses and even for these events, we obtained posteriors that are in agreement.

We note that the posteriors we obtained for the catalog analysis are also consistent with the ones generated in Finke et al. (2021) from the DarkSirensStat code3. However, the results in Finke et al. (2021) are generated with a different choice of BBHs and galaxy catalogue descriptions compared to the analyses performed in Abbott et al. (2023b).

We also tested the H0-catalog analysis on GW190814, one of the best localized and close-by dark sirens. For this event, we compare with the new version of GWCOSMO (Gray et al. 2023).

We performed two analyses, the first one using only the GW event and the spectral siren method, the second using galaxies reported in the K-band from GLADE+. In contrast to the previous analysis using the galaxy catalog, for GW190814, we used a pixel size (after computing the apparent magnitude cut) of 3 deg2 because GW190814 localization area is a few tens of square degrees. Following Abbott et al. (2023b), we assumed that GW190814 is an NSBH and we used a source mass POWER LAW + PEAK source mass model with the same population parameters used in Abbott et al. (2023b). In Fig. 5, we show the comparison between GWCOSMO and ICAROGW for the spectral and catalog analyses. The two results are in excellent agreement, thus indicating the robustness of our analysis.

thumbnail Fig. 4

Posterior probability density distributions obtained by ICAROGW 2.0 (blue line) from 42 BBHs used in Abbott et al. (2023b) in comparison with gwcosmo (orange dashed line).

5.3 Sources with observed EM counterparts

To test the electromagnetic counterpart method, we inferred the Hubble constant using the BNS merger GW170817 and its EM counterpart. We performed two types of analyses, one using a “high-latency” approach and the other using a “low-latency” approach. For both cases, the injection set used to evaluate the selection biases is given by the BNS injection set released for 03 sensitivity in Abbott et al. (2023c).

For the high-latency approach, we used IMRPhenom PE samples from Abbott et al. (2021b), consistently with the analysis of GWCOSMO, while for the low-latency approach we use the GW170817 skymap publicly released for GWTC-24. For both analyses, we used a sky position for the EM counterpart at right ascension 197.6 deg and declination of −23.4 deg. Moreover, we assumed a recessional velocity in the Hubble flow of vH = 3017 km s−1 with an uncertainty of σv = 166 km s−1 (Abbott et al. 2021b).

The analysis done with the code gwcosmo in Abbott et al. (2021b) used slightly different population and cosmological models to what it is currently implemented in ICAROGW. Therefore, to make a fair comparison, we implemented in ICAROGW the population and cosmological models used in Abbott et al. (2021b). First, the luminosity distance was approximated using linear cosmology, namely: dL=dc=czH0.${d_{\rm{L}}} = {d_c} = {{cz} \over {{H_0}}}.$(27)

The assumption corresponds to the following relations for the jacobians between source and detector frame and the differential of the comoving volume: dLz=cH0,${{\partial {d_{\rm{L}}}} \over {\partial z}} = {c \over {{H_0}}},$(28) zdL=H0c,${{\partial z} \over {\partial {d_{\rm{L}}}}} = {{{H_0}} \over c},$(29) Vcz=4πc3z2H03.${{\partial {V_{\rm{c}}}} \over {\partial z}} = 4\pi {{{c^3}{z^2}} \over {H_0^3}}.$(30)

Second, The CBC merger rate model for GW170817 used in the analysis was uniform in detector masses, namely p(m1,d,m2,d)=Θ(m2,d<m1,d)2(md,maxmd,min)2,$p\left( {{m_{1,d}},{m_{2,d}}} \right) = {{\Theta \left( {{m_{2,d}} < {m_{1,d}}} \right)} \over {2{{\left( {{m_{{\rm{d}},\max }} - {m_{{\rm{d}},\min }}} \right)}^2}}},$(31)

where the Θ function ensures that the detector’s secondary mass is lighter than the primary one. The analysis also neglected the 1/1 + z factor from the difference between source and detector frames. The overall merger rate was: dNCBCdtdddLdm1 dm2=dNCBCdtdddLp(m1,d,m2,d)       =R0VczzdLp(m1,d,m2,d)=R04πc2z2H02Θ(m2,d<m1,d)2(md,maxmd,min)2.$\matrix{ {{{{\rm{d}}{N_{{\rm{CBC}}}}} \over {{\rm{d}}{t_{\rm{d}}}{\rm{d}}{d_{\rm{L}}}{\rm{d}}{m_1}{\rm{d}}{m_2}}} = {{{\rm{d}}{N_{{\rm{CBC}}}}} \over {{\rm{d}}{t_{\rm{d}}}{\rm{d}}{d_{\rm{L}}}}}p\left( {{m_{1,d}},{m_{2,d}}} \right)} \hfill \cr {\,\,\,\,\,\,\quad = {R_0}{{\partial {V_{\rm{c}}}} \over {\partial z}}{{\partial z} \over {\partial {d_{\rm{L}}}}}p\left( {{m_{1,d}},{m_{2,d}}} \right) = {R_0}4\pi {{{c^2}{z^2}} \over {H_0^2}}{{\Theta \left( {{m_{2,d}} < {m_{1,d}}} \right)} \over {2{{\left( {{m_{{\rm{d}},\max }} - {m_{{\rm{d}},\min }}} \right)}^2}}}.} \hfill \cr } $(32)

We remark that the aforementioned assumptions on cosmology and rate model are not expected to provide a noticeable difference when calculating the weights wi,j for the GW170817 PE samples. This is because GW 170817 is a very close-by GW event, and even for extreme values of H0, it remains at low redshift where the linear cosmology approximation is enough. Moreover, assumptions about the masses for GW170817 are not expected to strongly bias the result in the presence of an EM counterpart, as shown in Mastrogiovanni et al. (2021). However, both masses and cosmological assumptions are expected to have an impact on the calculation of the selection bias. With 03 sensitivities, BNSs are detected up to a luminosity distance of ~300–400 Mpc, where the linear cosmology approximation can fail (especially if a high value of H0 is chosen). So, even when reproducing GW170817, it is important to consider the rate model assumed in Abbott et al. (2023c). Thus, we implemented the rate model in Eq. (32) to study this comparison.

Figure 6 shows the posteriors that we obtain with ICAROGW 2.0 and the method highlighted in Sect. 4 for GW170817 in comparison with gwcosmo. The posteriors for the high-latency analyses are in good agreement with each other. The main difference that we have found is that the low-latency posterior based on the GW skymap supports higher values of H0. The result is motivated as follows. By fixing the sky localization of GW170817 EM counterpart, the GW PE samples and skymap have different support for the luminosity distance location. In Fig. 7, we show that the skymap has a posterior distribution approximated as a Gaussian and it peaks at lower distance values if compared to the PE samples. At a fixed value of redshift for the EM counterpart, a lower luminosity distance can be fit using a higher value of H0. This is the motivation for which the low-latency analysis supports slightly higher values of H0. Concerning the high-latency analyses, the ICAROGW analysis excellently agrees with the GWCOSMO analysis of Abbott et al. (2021b).

thumbnail Fig. 5

Hubble constant posteriors generated for the well-localized event GW190814. The plot reports results generated with ICAROGW and the GWCOSMO code for two cases: the spectral siren and galaxy catalog cases (see legend).

thumbnail Fig. 6

Posterior distributions for H0 obtained from GW170817 with ICAROGW 2.0 with the high-latency approach (blue line) and low-latency approach (green dotted line) in comparison with gwcosmo using the high-latency approach (orange dash line). The posterior of gwcosmo is taken from Abbott et al. (2021b).

thumbnail Fig. 7

Luminosity distance posterior probability density function conditioned on the EM counterpart line-of-sight for GW 170817 using PE samples (orange histogram) and the low-latency skymap (blue line).

6 Conclusions and future developments

In this paper, we present ICAROGW 2.0, a python software for population properties inference in the presence of selection biases. We describe some of the tests performed to check the validity of the CBC rate models implemented. We show that the results obtained by ICAROGW 2.0 with the spectral sirens method are consistent with the previous version of the code (with and without the use of modified gravity models). We also show that the galaxy catalog analysis matches results generated by the code GWCOSMO. Finally, we show how EM counterparts can be employed with the CBC rate models to infer the vale of the Hubble constant. We discuss two approaches from the EM counterpart case that are able to produce consistent results with GWCOSMO.

ICAROGW 2.0 can be easily adapted to any custom population inference problem involving the presence of noisy measurements and selection biases. Future development plans in GW science for ICAROGW include more realistic models for CBCs that might include correlations among different variables (e.g. mass and redshift), the inclusion of more beyond-GR models, and time-delay models. The latest version of ICAROGW 2.0 is available to use in a public GitHub 5 repository.

Acknowledgements

The authors are grateful for computational resources provided by the LIGO Laboratory and supported by National Science Foundation Grants PHY-0757058 and PHY-0823459. This research has made use of data or software obtained from the Gravitational Wave Open Science Center (gwosc.org), a service of LIGO Laboratory, the LIGO Scientific Collaboration, the Virgo Collaboration, and KAGRA. LIGO Laboratory and Advanced LIGO are funded by the United States National Science Foundation (NSF) as well as the Science and Technology Facilities Council (STFC) of the United Kingdom, the Max-Planck-Society (MPS), and the State of Niedersachsen/Germany for support of the construction of Advanced LIGO and construction and operation of the GEO600 detector. Additional support for Advanced LIGO was provided by the Australian Research Council. Virgo is funded, through the European Gravitational Observatory (EGO), by the French Centre National de Recherche Scientifique (CNRS), the Italian Istituto Nazionale di Fisica Nucleare (INFN) and the Dutch Nikhef, with contributions by institutions from Belgium, Germany, Greece, Hungary, Ireland, Japan, Monaco, Poland, Portugal, Spain. KAGRA is supported by Ministry of Education, Culture, Sports, Science and Technology (MEXT), Japan Society for the Promotion of Science (JSPS) in Japan; National Research Foundation (NRF) and Ministry of Science and ICT (MSIT) in Korea; Academia Sinica (AS) and National Science and Technology Council (NSTC) in Taiwan. R.G. was supported by ERC starting grant SHADE 949572 and STFC grant ST/V005634/1. Software packages: ICAROGW uses the public python packages astropy (Astropy Collaboration 2022), bilby (Ashton et al. 2019b; Ashton & Talbot 2021), cupy (Okuta et al. 2017), h5py (Collette 2013), healpy (Górski et al. 2005; Zonca et al. 2019), numpy (Harris et al. 2020), pickle (Van Rossum 2020), and scipy (Virtanen et al. 2020) and their dependencies. This paper has used plotting utilities from the python packages chainconsumer (Hinton 2016) and matplotlib (Hunter 2007).

Appendix A Cosmological and GR deviation models

Cosmological models and GR deviation models are handled by several classes and functions and are organized in high-level wrappers for quick use. They can be passed to the various CBC rate models.

In Table A.1 we provide an overview of all the cosmo-logical and GR models available. We note that GR deviation models are extensions of cosmological models, with beyond-GR population parameters on top of the cosmological background population parameters. GR deviation models only override the way in which the GW luminosity distance is computed (see the sections below), while leaving the other cosmological quantities unchanged.

Table A.1

Summary table for all the background cosmology and models available in ICAROGW. More details on the models can be found in Sect. A.1.

Table A.2

Summary table for all the background cosmology and beyond-GR models in ICAROGW. More details on the models can be found in Sect. A.2.

Appendix A.1 Cosmological background models

In principle, ICAROGW is equipped to use all the cosmologies included in ASTROPY. However, for hierarchical inference, we have implemented only the models listed in the next subsections. For all the models, we calculate the GW luminosity distance: dL=c(1+z)H00zdzE(z),${d_{\rm{L}}} = {{c(1 + z)} \over {{H_0}}}\mathop \smallint \limits_0^z {{{\rm{d}}z'} \over {E\left( {z'} \right)}},$(A.1)

where H(z) = H0E(z), which is the same as the EM luminosity distance assuming GR (as described in this section). The differential of the luminosity distance is: dLz=dL(z)1+z+c(1+z)H01E(z)${{\partial {d_{\rm{L}}}} \over {\partial z}} = {{{d_{\rm{L}}}(z)} \over {1 + z}} + {{c(1 + z)} \over {{H_0}}}{1 \over {E(z)}}$(A.2)

The comoving volume is: Vc=0zdΩdzdVcdz,${V_{\rm{c}}} = \mathop \smallint \limits_0^z {\rm{d}}\Omega {\rm{d}}z'{{{\rm{d}}{V_{\rm{c}}}} \over {{\rm{d}}z'{\rm{d}}\Omega }},$(A.3)

and the differential of the comoving volume is: dVcdz=c3H031E(z)[ 0zdzE(z) ]2${{{\rm{d}}{V_{\rm{c}}}} \over {{\rm{d}}z{\rm{d}}\Omega }} = {{{c^3}} \over {H_0^3}}{1 \over {E(z)}}{\left[ {\mathop \smallint \limits_0^z {{{\rm{d}}z'} \over {E\left( {z'} \right)}}} \right]^2}$(A.4)

thumbnail Fig. A.1

Luminosity distance and its differential in redshift for various models available in ICAROGW. Top panel: Luminosity distance as a function of redshift for the modified gravity models. Bottom panel: Differential of the luminosity distance as a function of redshift for the modified gravity models (right). The lines indicate the models used to calculate the GW luminosity distance from the redshift.

The function E(z) depends on the cosmological model assumed.

For the Flat ACDM model: E2(z)=Ωm(1+z)3+(1Ωm),${E^2}(z) = {\Omega _m}{(1 + z)^3} + \left( {1 - {\Omega _m}} \right),$(A.5)

while for the Flat w0CDM model: E2(z)=Ωm(1+z)3+(1Ωm)(1+z)3(1+w0).${E^2}(z) = {\Omega _m}{(1 + z)^3} + \left( {1 - {\Omega _m}} \right){(1 + z)^{3\left( {1 + {w_0}} \right)}}.$(A.6)

Appendix A.2 Beyond-GR models

All the beyond-GR models implemented modify the luminosity distance, which we now refer to as dLGW$d_L^{{\rm{GW}}}$ (and its differential), while leaving the comoving volume untouched. We refer to the standard luminosity distance as dLEM$d_L^{{\rm{EM}}}$. In Fig. A.1, we show how the luminosity distance and its differential with respect to redshift are modified for the models described below.

Appendix A.2.1 The Ξ0 model

The luminosity distance is given by (see Eq. 2.31 of Belgacem et al. 2019): dLGW=dLEM(Ξ0+1Ξ0(1+z)n).$d_L^{{\rm{GW}}} = d_L^{{\rm{EM}}}\left( {{\Xi _0} + {{1 - {\Xi _0}} \over {{{(1 + z)}^n}}}} \right).$(A.7)

The Jacobian is given by: ddLGWdz=ddLEMdz(Ξ0+1Ξ0(1+z)n)dLEMn(1Ξ0)(1+z)n+1${{{\rm{d}}d_L^{{\rm{GW}}}} \over {{\rm{d}}z}} = {{{\rm{d}}d_L^{{\rm{EM}}}} \over {{\rm{d}}z}}\left( {{\Xi _0} + {{1 - {\Xi _0}} \over {{{(1 + z)}^n}}}} \right) - d_L^{{\rm{EM}}}{{n\left( {1 - {\Xi _0}} \right)} \over {{{(1 + z)}^{n + 1}}}}$(A.8)

Appendix A.2.2 The phenomenological log parametrization

The luminosity distance is given by: dLGW=dLEM[ 1+v=1n=3αvlogv(1+z) ].$d_L^{{\rm{GW}}} = d_L^{{\rm{EM}}}\left[ {1 + \mathop \sum \limits_{v = 1}^{n = 3} {\alpha _v}{{\log }^v}(1 + z)} \right]$(A.9)

The Jacobian is given by: ddLGWdz=ddLEMdzdLGWdLEM+dLEM[ ν=1n=3αvvlogν1(1+z)1+z].${{{\rm{d}}d_L^{{\rm{GW}}}} \over {{\rm{d}}z}} = {{{\rm{d}}d_L^{{\rm{EM}}}} \over {{\rm{d}}z}}{{d_L^{{\rm{GW}}}} \over {d_L^{{\rm{EM}}}}} + d_L^{{\rm{EM}}}\left[ {\mathop \sum \limits_{\nu = 1}^{n = 3} {\alpha _v}v{{{{\log }^{\nu - 1}}(1 + z)} \over {1 + z}}} \right]$(A.10)

Appendix A.2.3 Extra-dimensions

In the extra-dimensional model, the luminosity distance is given by (see Eq. 2.22 in Corman et al. 2022): dLGW=dLEM[ 1+(dLEM(1+z)Rc)n ]D42π.$d_{\rm{L}}^{{\rm{GW}}} = d_{\rm{L}}^{{\rm{EM}}}{\left[ {1 + {{\left( {{{d_{\rm{L}}^{{\rm{EM}}}} \over {(1 + z){R_c}}}} \right)}^n}} \right]^{{{D - 4} \over {2\pi }}}}$(A.11)

We define the following function: A=[ 1+(dLEM(1+z)Rc)n ].${\cal A} = \left[ {1 + {{\left( {{{d_{\rm{L}}^{{\rm{EM}}}} \over {(1 + z){R_c}}}} \right)}^n}} \right].$(A.12)

We also define the exponential: =D42n.${\cal E} = {{D - 4} \over {2n}}.$(A.13)

We can write the previous equation as: ddLGWdz=(A)ε[ ddLEMdz+nA(dLEMRc)n(ddLEMdz1(1+z)ndLEM(1+z)1+n) ].${{{\rm{d}}d_L^{{\rm{GW}}}} \over {{\rm{d}}z}} = {({\cal A})^\varepsilon }\left[ {{{{\rm{d}}d_L^{{\rm{EM}}}} \over {{\rm{d}}z}} + {{n{\cal E}} \over {\cal A}}{{\left( {{{d_{\rm{L}}^{{\rm{EM}}}} \over {{R_c}}}} \right)}^n}\left( {{{{\rm{d}}d_L^{{\rm{EM}}}} \over {{\rm{d}}z}}{1 \over {{{(1 + z)}^n}}} - {{d_{\rm{L}}^{{\rm{EM}}}} \over {{{(1 + z)}^{1 + n}}}}} \right)} \right].$(A.14)

Appendix A.2.4 The cM parametrization

Lastly, we consider a model with a running Planck mass (Lagos et al. 2019): dLGW=dLEMexp[ cM20z1(1+z)E2(z)dz ]dLEMexp[ cM2I(z) ],$d_{\rm{L}}^{{\rm{GW}}} = d_{\rm{L}}^{{\rm{EM}}}\exp \left[ {{{{c_M}} \over 2}\mathop \smallint \limits_0^z {1 \over {\left( {1 + z'} \right){E^2}\left( {z'} \right)}}dz'} \right] \equiv d_{\rm{L}}^{{\rm{EM}}}\exp \left[ {{{{c_M}} \over 2}I(z)} \right],$(A.15)

which defines I(z). In a flat ACDM model, I(z) can be calculated analytically and the result is (Eq. 19 in Lagos et al. 2019): dLGW=dLEMexp[ cM2ΩΛ,0ln1+z(Ωm,0(1+z)3+ΩΛ,0)1/3 ]$d_{\rm{L}}^{{\rm{GW}}} = d_{\rm{L}}^{{\rm{EM}}}\exp \left[ {{{{c_M}} \over {2{\Omega _{\Lambda ,0}}}}\ln {{1 + z} \over {{{\left( {{\Omega _{{\rm{m}},0}}{{(1 + z)}^3} + {\Omega _{\Lambda ,0}}} \right)}^{1/3}}}}} \right]{\rm{,}}$(A.16)

otherwise it needs to be calculated numerically. In any cosmology, the Jacobian is given by: ddLGWdz=ddLEMdzdLGWdLEM+dLGWcM21(1+z)E2(z).${{{\rm{d}}d_L^{{\rm{GW}}}} \over {{\rm{d}}z}} = {{{\rm{d}}d_L^{{\rm{EM}}}} \over {{\rm{d}}z}}{{d_L^{{\rm{GW}}}} \over {d_L^{{\rm{EM}}}}} + d_L^{{\rm{GW}}} \cdot {{{c_M}} \over 2} \cdot {1 \over {(1 + z){E^2}(z)}}.$(A.17)

Appendix B CBC Population models

Population models for the mass, redshift, and spins for CBC usually make use of probability density distributions. All the models currently available in ICAROGW are not conditionally dependent on each other, namely, the probability distributions of redshift, source masses, and spins are independent from each other. The redshift, mass, and spin models provided can be used to construct a CBC rate model.

We provide a list of models for the CBC merger rate in red-shift in Sect. B.1, for source masses in Sect. B.2, and those for spins in Sect. B.3.

Table B.1

Summary table for all the merger rate models available in ICAROGW. More details on the models can be found in Sect. B.1.

Table B.2

Summary table for all the background cosmology and models available in ICAROGW. The same mass models are also available for NSBH binaries. NSBH models also include three extra parameters: that are the minimum and maximum mass of the NS and the smoothing window for the lower end of the NS mass spectrum, respectively. More details on the models can be found in Sect. B.2.

Table B.3

Summary table for all the merger spin models available in ICAROGW. More details on the models can be found in Sect. B.3.

thumbnail Fig. B.1

Sample of rate models implemented in ICAROGW as a function of redshift.

Appendix B.1 CBC redshift rate evolution models

ICAROGW contains two models for the redshift evolution of the merger rate, as shown by Eqs. (15)-(16). Table B.1 summarises the merger rate redshift models, while Fig. B.1 provides some examples of the models for specific values of the parameters.

Appendix B.1.1 Power law

The rate is parametrized as: ψ(z;γ)=(1+z)γ.$\psi (z;\gamma ) = {(1 + z)^\gamma }.$(B.1)

Appendix B.1.2 Madau-Dickinson

The rate is parametrized following the star formation rate of Madau & Dickinson (2014) as: ψ(z;γ)=[ 1+(1+zp)γk ](1+z)γ1+(1+z1+zp)γ+k.$\psi (z;\gamma ) = \left[ {1 + {{\left( {1 + {z_p}} \right)}^{ - \gamma - k}}} \right]{{{{(1 + z)}^\gamma }} \over {1 + {{\left( {{{1 + z} \over {1 + {z_p}}}} \right)}^{\gamma + k}}}}.$(B.2)

Appendix B.2 Mass models

Most of the mass models are composed of the Gaussian and power law distributions, which we report in the following. The simple truncated power law distribution is given by 𝒫(xa,b,α)={ 1Nxα,(a<x<b)0, (otherwise)  $(x\mid a,b,\alpha ) = \{ \matrix{ {{1 \over N}{x^\alpha },} \hfill & {(a < x < b)} \hfill \cr {0,} \hfill & {{\rm{(otherwise)}}} \hfill \cr } $(B.3)

and the normalization factor is given by: N={ 1α+1[ bα+1aα+1 ], if α1lnba, if α=1 $N = \{ \matrix{ {{1 \over {\alpha + 1}}\left[ {{b^{\alpha + 1}} - {a^{\alpha + 1}}} \right],} \hfill & {{\rm{if}}\alpha \ne - 1} \hfill \cr {\ln {b \over a},} \hfill & {{\rm{if}}\alpha = - 1} \hfill \cr } $(B.4)

The truncated Gaussian distribution is given by: G[a,b](xμ,σ)={ 1N1σ2πe(xμ)22σ2,a<x<b0, otherwise  ${{\cal G}_{[a,b]}}(x\mid \mu ,\sigma ) = \{ \matrix{ {{1 \over N}{1 \over {\sigma \sqrt {2\pi } }}{e^{ - {{{{(x - \mu )}^2}} \over {2{\sigma ^2}}}}},} \hfill & {a < x < b} \hfill \cr {0,} \hfill & {{\rm{otherwise}}} \hfill \cr } $(B.5)

The normalization factor is expressed through the error function. Then the normalization factor is: N=ab1σ2πe(xμ)22σ2dx=(aμ)/(σ2)(bμ)/(σ2)1πet2dt.$N = \mathop \smallint \limits_a^b {1 \over {\sigma \sqrt {2\pi } }}{e^{ - {{{{(x - \mu )}^2}} \over {2{\sigma ^2}}}}}dx = \mathop \smallint \limits_{(a - \mu )/(\sigma \sqrt 2 )}^{(b - \mu )/(\sigma \sqrt 2 )} {1 \over {\sqrt \pi }}{e^{ - {t^2}}}dt.$(B.6)

Using the symmetry of the integrand around x = μ (t = 0) and the definition of erf function6: erf(z)=2π0zet2dt,${\rm{erf}}(z) = {2 \over {\sqrt \pi }}\mathop \smallint \limits_0^z {e^{ - {t^2}}}dt,$(B.7)

it follows that: N=12(erf[ bμσ2 ]erf[ aμσ2 ]).$N = {1 \over 2}\left( {{\rm{erf}}\left[ {{{b - \mu } \over {\sigma \sqrt 2 }}} \right] - {\rm{erf}}\left[ {{{a - \mu } \over {\sigma \sqrt 2 }}} \right]} \right).$(B.8)

In ICAROGW, we factorize the prior on mass as: π(m1,s,m2,sΛ)=π(m1,sΛ)π(m2,sm1,s,Λ).$\pi \left( {{m_{1,s}},{m_{2,s}}\mid \Lambda } \right) = \pi \left( {{m_{1,s}}\mid \Lambda } \right)\pi \left( {{m_{2,s}}\mid {m_{1,s}},\Lambda } \right).$(B.9)

When dealing with a NSBH, the neutron star is assigned to m2,s and the distribution of m2,s will be a simple power law defined between a minimum and a maximum mass, which are different from the ones assumed for the black hole.

In some of the models, we also apply a smoothing factor to the lower end of the mass distribution at m = mmin: π(m1,s,m2,sΛ)=[ π(m1,sΛ)π(m2,sm1,s,Λ) ]×S(m1mmin,δm)S(m2mmin,δm),$\matrix{ {{\rm{\pi }}\left( {{m_{1,s}},{m_{2,s}}\mid \Lambda } \right) = \left[ {\pi \left( {{m_{1,s}}\mid \Lambda } \right)\pi \left( {{m_{2,s}}\mid {m_{1,s}},\Lambda } \right)} \right] \times } \hfill \cr {\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,S\left( {{m_1}\mid {m_{\min }},{\delta _m}} \right)S\left( {{m_2}\mid {m_{\min }},{\delta _m}} \right),} \hfill \cr } $(B.10)

where S is a sigmoid-like window function (Eqs. B6–B7 of Abbott et al. 2021b): S(m1,smmin,δm)=                ={ 0,(m<mmin)[ f(mmmin,δm)+1 ]1,(mminm<mmin+δm)1,(mmmin+δm) $\matrix{ {S\left( {{m_{1,s}}\mid {m_{\min }},{\delta _m}} \right) = } \hfill \cr {\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\, = \{ \matrix{ {0,} & {\left( {m < {m_{\min }}} \right)} \cr {{{\left[ {f\left( {m - {m_{\min }},{\delta _m}} \right) + 1} \right]}^{ - 1}},} & {\left( {{m_{\min }} \le m < {m_{\min }} + {\delta _m}} \right)} \cr {1,} & {\left( {m \ge {m_{\min }} + {\delta _m}} \right)} \cr } } \hfill \cr } $(B.11)

with f(m,δm)=exp(δmm+δmmδm).$f\left( {m',{\delta _m}} \right) = \exp \left( {{{{\delta _m}} \over {m'}} + {{{\delta _m}} \over {m' - {\delta _m}}}} \right).$(B.12)

When we apply this window, the priors are numerically renor-malized.

All the mass models presented in this section can be visualized in Fig. B.2.

thumbnail Fig. B.2

Sample of mass models implemented in ICAROGW as distributions of primary and secondary source masses. The histograms display samples generated from the mass models and the blue solid lines the analytical marginal distributions on the primary and secondary source masses.

Appendix B.2.1 Truncated power law

The truncated power law model is given by Eq. (B.9) with: π(m1,smmin,mmax,α)=𝒫(m1,smmin,mmax,α),${\rm{\pi }}\left( {{m_{1,s}}\mid {m_{\min }},{m_{\max }},\alpha } \right) = \left( {{m_{1,s}}\mid {m_{\min }},{m_{\max }}, - \alpha } \right),$(B.13) π(m2,smmin,m1,s,β)=𝒫(m2,smmin,m1,s,β),${\rm{\pi }}\left( {{m_{2,s}}\mid {m_{\min }},{m_{1,s}},\beta } \right) = \left( {{m_{2,s}}\mid {m_{\min }},{m_{1,s}},\beta } \right),$(B.14)

where 𝒫 is defined in Eq. (B.3).

Appendix B.2.2 Power law + peak

This model was proposed in Talbot et al. (2019) and it is given by Eq. (B.10) with: π(m1,smmin,mmax,α)=(1λ)𝒫(m1,smmin,mmax,α)+λG(m1,sμg,σ),(0λ1)$\matrix{ {\pi \left( {{m_{1,s}}\mid {m_{\min }},{m_{\max }},\alpha } \right) = (1 - \lambda )\left( {{m_{1,s}}\mid {m_{\min }},{m_{\max }}, - \alpha } \right) + } \hfill \cr {\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\lambda {\cal G}\left( {{m_{1,s}}\mid {\mu _g},\sigma } \right),\quad (0 \le \lambda \le 1)} \hfill \cr } $(B.15) π(m2,smmin,m1,s,β)=𝒫(m2,smmin,m1,s,β).$\pi \left( {{m_{2,s}}\mid {m_{\min }},{m_{1,s}},\beta } \right) = \left( {{m_{2,s}}\mid {m_{\min }},{m_{1,s}},\beta } \right).$(B.16)

thumbnail Fig. B.3

Representation of the spin components for a compact binary coalescence.

Appendix B.2.3 Broken power law

This model is based on Eq. (B.10) and basically consists of two truncated power-law distributions attached at the point b: b=mmin+(mmaxmmin)f,$b = {m_{\min }} + \left( {{m_{\max }} - {m_{\min }}} \right)f,$(B.17)

where f is a scalar in [0, 1], so b = mmin for f = 0. This model was proposed in Abbott et al. (2021b). The priors are the following: π(m1,smmin,mmax,α)=1N [𝒫(m1,smmin,b,α1)+ 𝒫(bmmin,b,α1)𝒫(bb,mmax,α2)𝒫(m1,sb,mmax,α2) ],π(m2,smmin,m1,s,β)=𝒫(m2,smmin,m1,s,β),$\matrix{ {\pi \left( {{m_{1,s}}\mid {m_{\min }},{m_{\max }},\alpha } \right) = {1 \over N}[\left( {{m_{1,s}}\mid {m_{\min }},b, - {\alpha _1}} \right) + } \hfill \cr {{{\left( {b\mid {m_{\min }},b, - {\alpha _1}} \right)} \over {\left( {b\mid b,{m_{\max }}, - {\alpha _2}} \right)}}\left( {{m_{1,s}}\mid b,{m_{\max }}, - {\alpha _2}} \right)],} \hfill \cr {\pi \left( {{m_{2,s}}\mid {m_{\min }},{m_{1,s}},\beta } \right) = \left( {{m_{2,s}}\mid {m_{\min }},{m_{1,s}},\beta } \right),} \hfill \cr } $(B.18)

where the new normalization factor, N, here is: N=1+𝒫(bmmin,b,α1)𝒫(bb,mmax,α2)$N = 1 + {{\left( {b\mid {m_{\min }},b, - {\alpha _1}} \right)} \over {\left( {b\mid b,{m_{\max }}, - {\alpha _2}} \right)}}$(B.19)

Appendix B.2.4 Multi-peak model

This model is based on Eq. (B.10) and consists of one power-law + two Gaussian models with: π(m1,smmin,mmax,α)= [(1λ)𝒫(m1,smmin,mmax,α)+λλlow G(m1,sμg, low ,σlow )++λ(1λlow )G(m1,sμg, high ,σhigh ) ],$\eqalign{ & \pi \left( {{m_{1,s}}\mid {m_{\min }},{m_{\max }},\alpha } \right) = [(1 - \lambda ){\cal P}\left( {{m_{1,s}}\mid {m_{\min }},{m_{\max }}, - \alpha } \right) + \cr & \,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\lambda {\lambda _{{\rm{low }}}}{\cal G}\left( {{m_{1,s}}\mid {\mu _{g,{\rm{ low }}}},{\sigma _{{\rm{low }}}}} \right) + \cr & \,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\left. { + \lambda \left( {1 - {\lambda _{{\rm{low }}}}} \right){\cal G}\left( {{m_{1,s}}\mid {\mu _{g,{\rm{ high }}}},{\sigma _{{\rm{high }}}}} \right)} \right], \cr} $(B.20) π(m2,smmin,m1,s,β)=𝒫(m2,smmin,m1,s,β).$\pi \left( {{m_{2,s}}\mid {m_{\min }},{m_{1,s}},\beta } \right) = \left( {{m_{2,s}}\mid {m_{\min }},{m_{1,s}},\beta } \right).$(B.21)

We note that since 𝒢 and 𝒫 are already normalized, π(m1,s |mmin, mmax, α) is then automatically normalized. This model was used in Abbott et al. (2021b).

Appendix B.3 Spin models

In ICAROGW we implemented two models for the CBC spins. The two models are based on two different parametrizations for the spin parameters of a binary. Referring to Fig. B.3, we provide a definition for the spin parameters typically employed in GW studies.

By definition, the z axis of a binary is aligned to the instantaneous orbital angular momentum L. The (normalized) spin

amplitudes χ1,2, defined from the Cartesian components of the spin vectors, are expressed as: χ1=s1,x2+s1,y2+s1,z2${\chi _1} = \sqrt {s_{1,x}^2 + s_{1,y}^2 + s_{1,z}^2} {\rm{,}}$(B.22) χ2=s2,x2+s2,y2+s2,z2.${\chi _2} = \sqrt {s_{2,x}^2 + s_{2,y}^2 + s_{2,z}^2} .$(B.23)

The tilt angles θ1,2 are defined as the angle between the BH spins and the orbital angular momentum, namely: cosθ1=s1,zχ1$\cos {\theta _1} = {{{s_{1,z}}} \over {{\chi _1}}}$(B.24) cosθ2=s2,zχ2.$\cos {\theta _2} = {{{s_{2,z}}} \over {{\chi _2}}}.$(B.25)

The effective spin parameter χeff and precession spin parameter χp are defined by (Abbott et al. 2023c): χeff=χ1cosθ1+qχ2cosθ21+q=s1,z+qs2,z1+q,${\chi _{{\rm{eff}}}} = {{{\chi _1}\cos {\theta _1} + q{\chi _2}\cos {\theta _2}} \over {1 + q}} = {{{s_{1,z}} + q{s_{2,z}}} \over {1 + q}},$(B.26) χp=max[ χ1sinθ1;(4q+33q+4)qχ2sinθ2 ],${\chi _{\rm{p}}} = \max \left[ {{\chi _1}\sin {\theta _1};\left( {{{4q + 3} \over {3q + 4}}} \right)q{\chi _2}\sin {\theta _2}} \right],$(B.27)

where the mass ratio q is: q=m2m1,(q1).$q = {{{m_2}} \over {{m_1}}},\quad (q \le 1).$(B.28)

(The factors of q appearing in the expression for χp come from the leading order PN equation for L˙$\dot L$. We note that the four in-plane components of the perpendicular components of s have been replaced by one scalar χp, which is an averaged quantity, as shown in Eq. (3.1) of Schmidt et al. (2015).

The parameter χeff accounts for the amount of spin aligned with the orbital angular momentum, as well as the magnitude of the BH spins. Since χeff is bounded between [−1, 1], values close to χeff = 1 correspond to highly spinning BHs with aligned spins, whereas values close to χeff = −1 support highly spinning BHs with anti-aligned spins and χeff = 0 is consistent with non spinning BHs. The precession spin parameter χp, bounded between [0, 1], quantifies the amount of spin perpendicular to the angular momentum. All the spin models described below can be visualized in Fig. B.4.

Appendix B.3.1 DEFAULT spin model

This model was used in Abbott et al. (2021b) and it is proposed after Talbot et al. (2019); Wysocki et al. (2019). The model works with spins parameterized using the two spin magnitudes, χ1, χ2, and the two cosine of inclination angles cos, t1, cos t2, with respect to the orbital angular momentum. The total number of degrees of freedom (d.o.f) of a BBH system in terms of spin is 6. The last two remaining d.o.f are the azimuthal angles ϕ1 and ϕ2 are not considered here and supposed uniform. The population distribution is given by: π(χ1,χ2,cosθ1,cosθ2)=Beta(χ1α,β)π(cosθ1ξ,σt)×Beta(χ2α,β)π(cosθ2ξ,σt),$\matrix{ {\pi \left( {{\chi _1},{\chi _2},\cos {\theta _1},\cos {\theta _2}} \right) = {\rm{Beta}}\left( {{\chi _1}\mid \alpha ,\beta } \right)\pi \left( {\cos {\theta _1}\mid \xi ,{\sigma _t}} \right) \times } \hfill \cr {\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,{\rm{Beta}}\left( {{\chi _2}\mid \alpha ,\beta } \right)\pi \left( {\cos {\theta _2}\mid \xi ,{\sigma _t}} \right),} \hfill \cr } $(B.29)

namely, it is factored into two parts. Above, "beta" refers to the beta distribution, calculated with parameters α and β defined by: α=(1μχσχ21μχ)μχ21,$\alpha = \left( {{{1 - {\mu _\chi }} \over {\sigma _\chi ^2}} - {1 \over {{\mu _\chi }}}} \right)\mu _\chi ^2 \ge 1,$(B.30) β=α(1μχ1)1. $\beta = \alpha \left( {{1 \over {{\mu _\chi }}} - 1} \right) \ge 1.{\rm{}}$(B.31)

thumbnail Fig. B.4

Sample of spin models implemented in ICAROGW as a function of different spin parameters.

The condition (α,β) ≥ 1 is imposed to avoid any non-singular asymptotic behavior of the Beta distribution. The probability density function for the angle distribution is given by (as in Eq. (14) in Abbott et al. 2023c): π(cosθ1,2ζ,σt)=ξG[1,1](cosθ1,21,σt)+1ξ2,$\pi \left( {\cos {\theta _{1,2}}\mid \zeta ,{\sigma _t}} \right) = \xi {{\cal G}_{[ - 1,1]}}\left( {\cos {\theta _{1,2}}\mid 1,{\sigma _t}} \right) + {{1 - \xi } \over 2},$(B.32)

where 𝒢[−1, 1](cos θi| 1, σt) (see Eq.(B.5)) is a truncated Gaussian between −1 and 1 on cos θi, with mean 1 and standard deviation σt. We note that the parameter ξ is bounded between [0, 1]: the form of the angle distribution is a mixed model between a truncated Gaussian and a uniform distribution between −1 and 1, where ξ is the mixing parameter.

Appendix B.3.2 Gaussian spin model

The Gaussian spin model seeks to measure the joint distribution of χeff and χp. It was proposed in Miller et al. (2020) and it depends on five parameters: μχeff,σχeff,μχp,σχp${\mu _{{\chi _{{\rm{eff}}}}}},{\sigma _{{\chi _{{\rm{eff}}}}}},{\mu _{{\chi _{\rm{p}}}}},{\sigma _{{\chi _{\rm{p}}}}}$, and ρ. The population probability on χeff,χp is a bivariate Gaussian truncated between [−1, 1] for χeff and between [0, 1] for χp. The covariance of the bivariate Gaussian is cov[ χeff,χp ]=ρσχeffσχp${{\mathop{\rm cov}} _{\left[ {{\chi _{{\rm{eff}}}},{\chi _{\rm{p}}}} \right]}} = \rho {\sigma _{{\chi _{{\rm{eff}}}}}}{\sigma _{{\chi _{\rm{p}}}}}$. In ICAROGW, this bivariate Gaussian is factorized as: π( χeff,χp μχeff,σχeff,μχp,σχp,ρ )=                       G[1,1](χeffμχeff,σχeff)G[0,1](χpμ*,σ*),$\matrix{ {\pi ({\chi _{{\rm{eff}}}},{\chi _{\rm{p}}}\mid {\mu _{{\chi _{{\rm{eff}}}}}},{\sigma _{{\chi _{{\rm{eff}}}}}},{\mu _{{\chi _{\rm{p}}}}},{\sigma _{{\chi _{\rm{p}}}}},\rho ) = } \hfill \cr {\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,{{\cal G}_{[ - 1,1]}}\left( {{\chi _{{\rm{eff}}}}\mid {\mu _{{\chi _{{\rm{eff}}}}}},{\sigma _{{\chi _{{\rm{eff}}}}}}} \right){{\cal G}_{[0,1]}}\left( {{\chi _{\rm{p}}}\mid {\mu _*},{\sigma _*}} \right),} \hfill \cr } $(B.33)

where μ*=μχp+cov[ χeff,χp ]σχeff2(χeffμχeff),${\mu _*} = {\mu _{{\chi _{\rm{p}}}}} + {{{{{\mathop{\rm cov}} }_{\left[ {{\chi _{{\rm{eff}}}},{\chi _{\rm{p}}}} \right]}}} \over {\sigma _{{\chi _{{\rm{eff}}}}}^2}}\left( {{\chi _{{\rm{eff}}}} - {\mu _{{\chi _{{\rm{eff}}}}}}} \right),$(B.34) σ*=σχPcov[ χeff,χp ]σχeff2${\sigma _*} = {{{\sigma _{{\chi _{\rm{P}}}}}{{{\mathop{\rm cov}} }_{\left[ {{\chi _{{\rm{eff}}}},{\chi _{\rm{p}}}} \right]}}} \over {\sigma _{{\chi _{{\rm{eff}}}}}^2}}$(B.35)

This factorization is equivalent to a bivariate Gaussian distribution.

References

  1. Abbott, R., Abbott, T. D., Abraham, S., et al. 2021b, ApJ, 913, L7 [NASA ADS] [CrossRef] [Google Scholar]
  2. Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2021b, ApJ, 909, 218 [NASA ADS] [CrossRef] [Google Scholar]
  3. Abbott, R., Abbott, T. D., Acernese, F., et al. 2023a, Phys. Rev. X, 13, 041039 [Google Scholar]
  4. Abbott, R., Abe, H., Acernese, F., et al. 2023b, ApJ, 949, 76 [NASA ADS] [CrossRef] [Google Scholar]
  5. Abbott, R., Abbott, T., Acernese, F., et al. 2023c, Phys. Rev. X, 13, 011048 [NASA ADS] [Google Scholar]
  6. Ashton, G., & Talbot, C. 2021, MNRAS, 507, 2037 [NASA ADS] [CrossRef] [Google Scholar]
  7. Ashton, G., Hübner, M., Lasky, P., & Talbot, C. 2019a, https://doi.org/10.5281/zenodo.2602178 [Google Scholar]
  8. Ashton, G., Hübner, M., Lasky, P. D., et al. 2019b, ApJS, 241, 27 [NASA ADS] [CrossRef] [Google Scholar]
  9. Astropy Collaboration (Price-Whelan, A. M., et al.) 2022, ApJ, 935, 167 [NASA ADS] [CrossRef] [Google Scholar]
  10. Belgacem, E., Calgani, G., Crisostomi, M., et al. 2019, J. Cosmol. Astropart. Phys., 07, 024 [CrossRef] [Google Scholar]
  11. Bovy, J., Hogg, D. W., & Roweis, S. T. 2011, Ann. Appl. Stat., 5, 1657 [Google Scholar]
  12. Chen, H.-Y. 2020, Phys. Rev. Lett, 125, 201301 [NASA ADS] [CrossRef] [Google Scholar]
  13. Collette, A. 2013, Python and HDF5 (O’Reilly Media, Inc.) [Google Scholar]
  14. Corman, M., Ghosh, A., Escamilla-Rivera, C, et al. 2022, Phys. Rev. D, 105, 064061 [NASA ADS] [CrossRef] [Google Scholar]
  15. Dalya, G., Díaz, R., Bouchet, F. R., et al. 2022, MNRAS, 514, 1403 [CrossRef] [Google Scholar]
  16. Del Pozzo, W. 2012, Phys. Rev. D, 86, 043011 [NASA ADS] [CrossRef] [Google Scholar]
  17. Ezquiaga, J. M. 2021, Phys. Lett. B, 822, 136665 [NASA ADS] [CrossRef] [Google Scholar]
  18. Ezquiaga, J. M., & Holz, D. E. 2022, Phys. Rev. Lett, 129, 061102 [NASA ADS] [CrossRef] [Google Scholar]
  19. Farr, W. M. 2019, Res. Notes AAS, 3, 66 [NASA ADS] [CrossRef] [Google Scholar]
  20. Finke, A., Foffa, S., Iacovelli, F., Maggiore, M., & Mancarella, M. 2021, J. Cosmol. Astropart. Phys., 2021, 026 [CrossRef] [Google Scholar]
  21. Foreman-Mackey, D., Hogg, D. W., & Morton, T. D. 2014, ApJ, 795, 64 [Google Scholar]
  22. Gair, J. R., Ghosh, A., Gray, R., et al. 2023, AJ, 166, 22 [NASA ADS] [CrossRef] [Google Scholar]
  23. Gonzalez, A. H., & Faber, S. M. 1997, ApJ, 485, 80 [NASA ADS] [CrossRef] [Google Scholar]
  24. Górski, K. M., Hivon, E., Banday, A. J., et al. 2005, ApJ, 622, 759 [Google Scholar]
  25. Gray, R., Hernandez, I. M., Qi, H, et al. 2020, Phys. Rev. D, 101, 122001 [NASA ADS] [CrossRef] [Google Scholar]
  26. Gray, R., Messenger, C, & Veitch, J. 2022, MNRAS, 512, 1127 [NASA ADS] [CrossRef] [Google Scholar]
  27. Gray, R., Beirnaert, F., Karathanasis, C, et al. 2023, J. Cosmol. Astropart. Phys., 2023, 023 [CrossRef] [Google Scholar]
  28. Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357 [NASA ADS] [CrossRef] [Google Scholar]
  29. Hinton, S. R. 2016, J. Open Source Softw., 1, 00045 [NASA ADS] [CrossRef] [Google Scholar]
  30. Hunter, J. D. 2007, Comput. Sci. Eng., 9, 90 [NASA ADS] [CrossRef] [Google Scholar]
  31. Karathanasis, C, Mukherjee, S., & Mastrogiovanni, S. 2023a, MNRAS, 523, 4539 [NASA ADS] [Google Scholar]
  32. Karathanasis, C, Revenu, B., Mukherjee, S., & Stachurski, F. 2023b, A&A, 677, A124 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  33. Lagos, M., Fishbach, M., Landry, P., & Holz, D. E. 2019, Phys. Rev. D, 99, 083504 [NASA ADS] [CrossRef] [Google Scholar]
  34. Leyde, K., Mastrogiovanni, S., Steer, D. A., Chassande-Mottin, E., & Karathanasis, C. 2022, J. Cosmol. Astropart. Phys., 09, 012 [CrossRef] [Google Scholar]
  35. Liu, B., Li, Z., Zhao, S., Zhou, H., & Gao, H. 2023, ApJ, 943, 29 [NASA ADS] [CrossRef] [Google Scholar]
  36. Loredo, T. J., & Lamb, D. Q. 2002, Phys. Rev. D, 65, 063002 [NASA ADS] [CrossRef] [Google Scholar]
  37. Loredo, T. J., & Wasserman, I. M. 1998, ApJ, 502, 75 [NASA ADS] [CrossRef] [Google Scholar]
  38. Madau, P., & Dickinson, M. 2014, ARA&A, 52, 415 [Google Scholar]
  39. Malmquist, K. G. 1922, MeLuF, 100, 1 [NASA ADS] [Google Scholar]
  40. Mancarella, M., & Genoud-Prachex, E. 2022, https://doi.org/10.5281/zenodo.6323173 [Google Scholar]
  41. Mancarella, M., Genoud-Prachex, E., & Maggiore, M. 2022, Phys. Rev. D, 105, 064030 [NASA ADS] [CrossRef] [Google Scholar]
  42. Mandel, I., Farr, W. M., & Gair, J. R. 2019, MNRAS, 486, 1086 [Google Scholar]
  43. Mastrogiovanni, S., Leyde, K., Karathanasis, C, et al. 2021, Phys. Rev. D, 104, 062009 [NASA ADS] [CrossRef] [Google Scholar]
  44. Mastrogiovanni, S., Laghi, D., Gray, R., et al. 2023, Phys. Rev. D, 108, 042002 [NASA ADS] [CrossRef] [Google Scholar]
  45. Miller, S., Callister, T. A., & Farr, W. M. 2020, ApJ, 895, 128 [NASA ADS] [CrossRef] [Google Scholar]
  46. Okuta, R., Unno, Y., Nishino, D., Hido, S., & Loomis, C. 2017, in Proceedings of Workshop on Machine Learning Systems (LearningSys) in The Thirty-first Annual Conference on Neural Information Processing Systems (NIPS) [Google Scholar]
  47. Schmidt, P., Ohme, F., & Hannam, M. 2015, Phys. Rev. D, 91, 024043 [NASA ADS] [CrossRef] [Google Scholar]
  48. Schutz, B. F. 1986, Nature, 323, 310 [Google Scholar]
  49. Talbot, C, & Golomb, J. 2023, MNRAS, 526, 3495 [NASA ADS] [CrossRef] [Google Scholar]
  50. Talbot, C, Smith, R., Thrane, E., & Poole, G. B. 2019, Phys. Rev. D, 100, 043030 [NASA ADS] [CrossRef] [Google Scholar]
  51. Turski, C, Bilicki, M., Dálya, G., Gray, R., & Ghosh, A. 2023, MNRAS, 526, 6224 [NASA ADS] [CrossRef] [Google Scholar]
  52. Van Rossum, G. 2020, The Python Library Reference, release 3.8.2 (Python Software Foundation) [Google Scholar]
  53. Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nat. Methods, 17, 261 [Google Scholar]
  54. Vitale, S., Gerosa, D., Farr, W. M., & Taylor, S. R. 2022, in Handbook of Gravitational Wave Astronomy (Berlin: Springer), 45 [Google Scholar]
  55. Winn, J. N., & Fabrycky, D. C. 2015, ARA&A, 53, 409 [Google Scholar]
  56. Wysocki, D., Lange, J., & O'Shaughnessy, R. 2019, Phys. Rev. D, 100, 043012 [NASA ADS] [CrossRef] [Google Scholar]
  57. Zheng, L.-M., Li, Z., Chen, Z.-C, Zhou, H., & Zhu, Z.-H. 2023, Phys. Lett. B, 838, 137720 [NASA ADS] [CrossRef] [Google Scholar]
  58. Zonca, A., Singer, L., Lenz, D., et al. 2019, J. Open Source Softw., 4, 1298 [Google Scholar]

6

Scipy function.

All Tables

Table A.1

Summary table for all the background cosmology and models available in ICAROGW. More details on the models can be found in Sect. A.1.

Table A.2

Summary table for all the background cosmology and beyond-GR models in ICAROGW. More details on the models can be found in Sect. A.2.

Table B.1

Summary table for all the merger rate models available in ICAROGW. More details on the models can be found in Sect. B.1.

Table B.2

Summary table for all the background cosmology and models available in ICAROGW. The same mass models are also available for NSBH binaries. NSBH models also include three extra parameters: that are the minimum and maximum mass of the NS and the smoothing window for the lower end of the NS mass spectrum, respectively. More details on the models can be found in Sect. B.2.

Table B.3

Summary table for all the merger spin models available in ICAROGW. More details on the models can be found in Sect. B.3.

All Figures

thumbnail Fig. 1

ICAROGW logical structure. The orange boxes highlight the computation of the rate weights which are needed for the calculation of the hierarchical likelihood given in Eq. (9).

In the text
thumbnail Fig. 2

ICAROGW computational benchmarking on CPU and GPU. Top: average elapsed time taken by ICAROGW to evaluate the hierarchical likelihood in Eq. (1). Bottom: efficiency, defined as number of samples processed per second to evaluate the hierarchical likelihood in Eq. (1). For both panels, the CPU is reported with blue pentagons and GPU with orange squares. The two panels are reported as a function of a number of samples, namely, the sum of the number of PE samples and injections used for the computation. The CPU markers are generated up to 106 samples since at 105 samples the efficiency is already saturated. Figures generated with double float precision on CPU Intel Core i9-11950H (8 cores HT, 2.6–5.0 GHz Turbo) and GPU NVIDIA GeForce RTX3080 (16Gb GDDR6, 6144 cores CUDA).

In the text
thumbnail Fig. 3

Quantile-quantile plots of the marginal posterior distributions of the CBC merger rate population parameters inferred from the 42 BBHs in GWTC-3. The plots columns correspond to the model used to calculate the GW luminosity distance as a function of redshift, see Appendices A.1A.2. The legend and different lines report the population parameters. All the benchmark runs use the POWER LAW + PEAK model for the source mass (cf. Appendix B.2) and the MADAU-DICKINSON merger rate model (cf. Appendix B.1).

In the text
thumbnail Fig. 4

Posterior probability density distributions obtained by ICAROGW 2.0 (blue line) from 42 BBHs used in Abbott et al. (2023b) in comparison with gwcosmo (orange dashed line).

In the text
thumbnail Fig. 5

Hubble constant posteriors generated for the well-localized event GW190814. The plot reports results generated with ICAROGW and the GWCOSMO code for two cases: the spectral siren and galaxy catalog cases (see legend).

In the text
thumbnail Fig. 6

Posterior distributions for H0 obtained from GW170817 with ICAROGW 2.0 with the high-latency approach (blue line) and low-latency approach (green dotted line) in comparison with gwcosmo using the high-latency approach (orange dash line). The posterior of gwcosmo is taken from Abbott et al. (2021b).

In the text
thumbnail Fig. 7

Luminosity distance posterior probability density function conditioned on the EM counterpart line-of-sight for GW 170817 using PE samples (orange histogram) and the low-latency skymap (blue line).

In the text
thumbnail Fig. A.1

Luminosity distance and its differential in redshift for various models available in ICAROGW. Top panel: Luminosity distance as a function of redshift for the modified gravity models. Bottom panel: Differential of the luminosity distance as a function of redshift for the modified gravity models (right). The lines indicate the models used to calculate the GW luminosity distance from the redshift.

In the text
thumbnail Fig. B.1

Sample of rate models implemented in ICAROGW as a function of redshift.

In the text
thumbnail Fig. B.2

Sample of mass models implemented in ICAROGW as distributions of primary and secondary source masses. The histograms display samples generated from the mass models and the blue solid lines the analytical marginal distributions on the primary and secondary source masses.

In the text
thumbnail Fig. B.3

Representation of the spin components for a compact binary coalescence.

In the text
thumbnail Fig. B.4

Sample of spin models implemented in ICAROGW as a function of different spin parameters.

In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.