11institutetext:  1 INFN, Sezione di Roma, I-00185 Roma, Italy
2 Université Lyon, Université Claude Bernard Lyon 1, CNRS, IP2I Lyon/IN2P3,UMR
5822, F-69622 Villeurbanne, France
3 Laboratoire des 2 Infinis - Toulouse (L2IT-IN2P3), Université de Toulouse,
CNRS, UPS, F-31062 Toulouse Cedex 9, France
4 Institut de Física d’Altes Energies (IFAE), Barcelona Institute of Science
and Technology, Barcelona, Spain
5 Department of Physics and Astronomy, Ghent University, Proeftuinstraat 86,
9000 Ghent, Belgium
6 SUPA, University of Glasgow, Glasgow, G12 8QQ, United Kingdom
7 Department of Physics & Astronomy, Queen Mary University of London, Mile End
Road, London, E1 4NS, United Kingdom
8 Université Paris Cité, CNRS, Astroparticule et Cosmologie, F-75013 Paris,
France

#  icarogw: A python package for inference of astrophysical population
properties of noisy, heterogeneous and incomplete observations

Simone Mastrogiovanni1<EMAIL_ADDRESS>Grégoire Pierra2 Stéphane
Perriès2 Danny Laghi3 Giada Caneva Santoro4 Archisman Ghosh5 Rachel Gray6,7
Christos Karathanasis4 Konstantin Leyde8

We present icarogw 2.0 , a pure CPU/GPU python code developed to infer
astrophysical and cosmological population properties of noisy, heterogeneous,
and incomplete observations. icarogw 2.0 is mainly developed for compact
binary coalescence (CBC) population inference with gravitational wave (GW)
observations. The code contains several models for masses, spins, and redshift
of CBC distributions, and is able to infer population distributions as well as
the cosmological parameters and possible general relativity deviations at
cosmological scales. We present the theoretical and computational foundations
of icarogw, and we describe how the code can be employed for population and
cosmological inference using (i) only GWs, (ii) GWs and galaxy surveys and
(iii) GWs with electromagnetic counterparts. Although icarogw 2.0 has been
developed for GW science, we also describe how the code can be used for any
physical and astrophysical problem involving observations from noisy data in
presence of selection biases. With this paper, we also release tutorials on
Zenodo (Mastrogiovanni et al., 2023).

###### Key Words.:

Methods: data analysis, Methods: statistical, Cosmology: cosmological
parameters,Cosmology: observations, Gravitational waves

###### Contents

  1. 1 Introduction
  2. 2 Inhomogeneous Poisson process and Bayesian inference
  3. 3 Structure of the code
    1. 3.1 A custom example: Estimation of the mean of a gaussian process in presence of selection biases
  4. 4 Application to a compact binary coalescence case
    1. 4.1 Spectral sirens merger rates
    2. 4.2 Galaxy catalog merger rates
    3. 4.3 Multi-messenger sources merger rates
  5. 5 Sanity checks for CBC populations
  6. 6 Conclusions and future development
  7. A Cosmological and GR deviation models
    1. A.1 Cosmological background models
    2. A.2 Beyond-GR models
      1. A.2.1 The $\Xi_{0}$ model
      2. A.2.2 The phenomenological $\log$ parametrization
      3. A.2.3 Extra-dimensions
      4. A.2.4 The $c_{M}$ parametrization
  8. B CBC Population models
    1. B.1 CBC redshift rate evolution models
      1. B.1.1 Power Law
      2. B.1.2 Madau-Dickinson
    2. B.2 Mass models
      1. B.2.1 Truncated Power-Law
      2. B.2.2 Power-Law + Peak
      3. B.2.3 Broken Power Law
      4. B.2.4 Multi-Peak
    3. B.3 Spin models
      1. B.3.1 Default spin model
      2. B.3.2 Gaussian spin model
  9. C GPU implementation

##  1 Introduction

Inferring Cosmology and AstRophysics with Observations of Gravitational Waves
(icarogw) 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 long-standing
problem shared among several research topics. With almost 100 GW observations
from the last runs of the LIGO, Virgo, and KAGRA (LVK) collaboration (Abbott
et al., 2021), 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 firstly used for the LVK
analysis of Abbott et al. (2021a) and released on the LVK official GitLab page
. The first distribution of icarogw was also used in independent studies for
population code validations (Karathanasis et al., 2022; Turski et al., 2023),
beyond General Relativity (GR) (Leyde et al., 2022), astrophysical processes
(Karathanasis et al., 2022), and primordial black holes models (Liu et al.,
2023; Zheng et al., 2023).

In general, in astrophysics population inference involves the correction of a
_selection bias_ , or Malmquist bias (Malmquist, 1922), that prevents the
observation of a particular class of astrophysical processes, thus biasing the
population analysis. Selection biases are very common in many astrophysical
observations involving neutrino physics (Loredo & Lamb, 2002), exoplanets
(Foreman-Mackey et al., 2014; Winn & Fabrycky, 2015), galaxy and clusters
observations (Gonzalez & Faber, 1997), and $\gamma$-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
binary parameters such as masses and redshift.

Besides the correction of selection biases, a population inference for GW
signals should also account for the fact that the source parameters are
generated according to non-trivial distribution (data is _heterogeneous_) and
they are not perfectly measured (presence of noise in the data). Therefore,
when reconstructing the population properties of GW signals, we need to deal
with noisy, heterogeneous, and incomplete observations that require specific
statistical frameworks. Bovy et al. (2011) also refers to this type of
analysis as _extreme deconvolution_. One might be tempted to think that, even
in absence of selection biases, by “stacking” measures with their errors of a
given quantity, _e.g._ the CBC chirp mass $\mathcal{M}_{c}$, it is possible to
reconstruct the generative chirp mass distribution. However, this is not true
and it will lead the analysis to reconstruct a biased generative model for the
chirp mass. Let us provide a concrete example: if a random gaussian process
generates samples of $x$ which are then measured with gaussian uncertainties
$\bar{x}_{i}\pm\sigma$, then summing gaussians centered in $\bar{x}_{i}$ with
standard deviation (s.t.d.) $\sigma$ will not reconstruct the original
distribution. We, therefore, need a proper statistical framework to
reconstruct population properties. Current techniques for population inference
include the use of _hierarchical Bayesian inference_ (Mandel et al., 2019;
Vitale et al., 2022). icarogw 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 observations, (ii)
a set of simulations, or injections, to calculate the explorable volume in
parameter space and (iii) a rate, or parameter population distributions
problem.

In this paper, we present the version 2.0 of icarogw. The new icarogw updates
include:

  * •

A more user-friendly environment that is able to accommodate a quick
implementation of new population models without the necessity of entering most
of the code details.

  * •

An interface to calculate various numerical stability estimators for Bayesian
hierarchical inference.

  * •

Two models for the cosmological background expansion, four models for GR
modifications, 4 models for binary black hole (BBH) source mass distributions,
2 models for BBH spins distributions, and 2 models for BBH merger rate
evolution.

  * •

Three different methodologies based on electromagnetic (EM) counterparts, GWs
alone, or galaxy catalogs that can combine all the population models available
in icarogw.

This paper is organized as follows. In Sec. 2 we summarise the theoretical and
code implementation basics of hierarchical Bayesian inference. In Sec. 3 we
provide a general overview of icarogw’s structure and also discuss how it can
be used with a non-GW-related example of population inference. In Sec. 4 we
enter into the details on how icarogw can be used for population inference
using GWs alone (spectral sirens), galaxy catalogs, and electromagnetic
counterparts. In Sec. 6 we draw our conclusions and discuss prospects for
future development. Finally, in App. A.1 and App. A.2 we describe implemented
cosmological background and beyond-GR models, in App. B mass, redshift, and
spin models for CBCs, and in App. C the GPU (Graphical Processing Unit)
implementation of the code.

##  2 Inhomogeneous Poisson process and Bayesian inference

The main application of icarogw is to infer the population parameters
$\Lambda$ that describe the production rate of events in terms of parameters
$\theta$, namely $\frac{{\rm d}N}{{\rm d}t{\rm d}\theta}(\Lambda)$. For
instance, in Sec. 4 we will consider the rate of CBC in terms of source frame
masses. The hierarchical likelihood of obtaining $N_{\rm obs}$ observations,
each described by some parameters $\theta$, in a data collection $\\{x\\}$ for
a given observing time $T_{\rm obs}$ 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):

$\begin{split}\mathcal{L}(\\{x\\}|\Lambda)&\propto e^{-N_{\rm
exp}(\Lambda)}\prod_{i=1}^{N_{\rm
obs}}\int\mathcal{L}(x_{i}|\theta,\Lambda)\frac{{\rm d}N}{{\rm d}t{\rm
d}\theta}(\Lambda){\rm d}t{\rm d}\theta\\\ &\propto e^{-N_{\rm
exp}(\Lambda)}\prod_{i=1}^{N_{\rm obs}}T_{\rm
obs}\int\mathcal{L}(x_{i}|\theta,\Lambda)\frac{{\rm d}N}{{\rm d}t{\rm
d}\theta}(\Lambda){\rm d}\theta.\end{split}$ (1)

Eq. (1) is often referred to as hierarchical likelihood. By assuming a ”scale-
free” (i.e., neglecting information about the rate) prior $\pi(N_{\rm
exp})\propto 1/N_{\rm exp}$ on the expected number of detections, an
equivalent form of Eq. (1) can be derived:

$\mathcal{L}(x|\Lambda)\propto\prod_{i=1}^{N_{\rm
obs}}\frac{\int\mathcal{L}(x_{i}|\theta,\Lambda)\frac{{\rm d}N}{{\rm d}t{\rm
d}\theta}{\rm d}\theta}{\int p_{\rm det}(\theta,\Lambda)\frac{{\rm d}N}{{\rm
d}t{\rm d}\theta}{\rm d}\theta}.$ (2)

Both the hierarchical likelihoods in Eq. (1) and Eq. (2) contain several
crucial quantities to the inference problem. In what follows, we will explain
how icarogw numerically computes these quantities and uses them to compute the
full hierarchical likelihood defined in Eq. (1).

The first central quantity is the single-event likelihood
$\mathcal{L}(x_{i}|\theta,\Lambda)$: this term tells us how well we are able
to measure the parameters $\theta$. Typically, we are provided with $N_{\rm
s,i}$ PE posterior samples drawn from
$p(\theta|x_{i},\Lambda)\propto\mathcal{L}(x_{i}|\theta,\Lambda)\pi_{\rm
PE}(\theta|\Lambda)$, where $\pi_{\rm PE}(\theta|\Lambda)$ is the prior used
to generate the samples. It has to be noted that $\theta$ represents the
single-event parameters that we believe can be true. icarogw evaluates
numerically the likelihood integral in Eq. (1) and in the numerator of Eq. (2)
via Monte Carlo integration by summing over PE samples:

$\int\mathcal{L}(x_{i}|\theta,\Lambda)\frac{{\rm d}N}{{\rm d}t{\rm
d}\theta}(\Lambda){\rm d}\theta\approx\frac{1}{N_{{\rm
s},i}}\sum_{j=1}^{N_{{\rm s},i}}\frac{1}{\pi_{\rm
PE}(\theta_{i,j}|\Lambda)}\frac{dN}{dtd\theta}(\Lambda)\bigg{|}_{i,j}\equiv\frac{1}{N_{{\rm
s},i}}\sum_{j=1}^{N_{{\rm s},i}}w_{i,j},$ (3)

where the index $i$ refers to the event and the index $j$ to the posterior
samples of the events. We have also defined a weight $w_{i,j}$ of dimension
equal to the number of events generated per unit of time. As Eq. (3) is
evaluated with a finite sum over posterior samples, we introduce a numerical
stability estimator that is the effective number of posterior samples per
event $i$ as Farr (2019):

$N_{{\rm eff},i}=\frac{(\sum_{j}^{N_{{\rm
s},i}}w_{i,j})^{2}}{\sum_{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 Abbott
et al. (2021a), it is required to have at least an effective number of
posterior samples equal to 20 for each event and population model supported by
the analysis (although this is a flag that can be set in the code). 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 $N_{\rm
exp}(\Lambda)$, which is related to the selection bias and can be evaluated
as:

$N_{\rm exp}(\Lambda)=T_{\rm obs}\int p_{\rm
det}(\theta,\Lambda)\frac{dN}{{\rm d}t{\rm d}\theta}{\rm d}\theta,$ (5)

where $p_{\rm det}(\theta,\Lambda)$ is a detection probability that can be
calculated as:

$p_{\rm
det}(\theta,\Lambda)=\int_{x\in\rm{detectable}}\mathcal{L}(x_{i}|\theta,\Lambda){\rm
d}x.$ (6)

Typically, we do not have access to an analytical form of the detection
probability (see Gair et al. (2022) for an introductory example in the context
of GW cosmology with galaxy catalogs). 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 volume that can be explored in the
parameter space and correct for selection biases. Therefore, their occurrence
is proportional to $p_{\rm det}(\theta,\Lambda)$ and the population model used
to generate them. icarogw takes in input a set of $N_{\rm det}$ detected
injections out of $N_{\rm gen}$ total injections generated from a prior
$\pi_{\rm inj}(\theta)$ to calculate the integral in Eq. (5) using Monte Carlo
integration:

$N_{\rm exp}\approx\frac{T_{\rm obs}}{N_{\rm gen}}\sum_{j=1}^{N_{\rm
det}}\frac{1}{\pi_{\rm
inj}(\theta_{j})}\frac{dN}{dtd\theta}\bigg{|}_{j}\equiv\frac{T_{\rm
obs}}{N_{\rm gen}}\sum_{j=1}^{N_{\rm det}}s_{j}.$ (7)

Here we have again defined a weight $s_{j}$ with the dimension of a rate of
events. Note that there is one fundamental difference with Eq. (3). The
injection prior $\pi_{\rm inj}(\theta)$ must be properly normalized to obtain
a reasonable value of $N_{\rm exp}$, while a wrong normalization of $\pi_{\rm
PE}(\theta)$ (which is used in Eq. (3)) will only result in an overall
normalization factor to the overall hierarchical likelihood. Following Farr
(2019), also for Eq. (7) we can define a numerical stability estimator, the
effective number of injections:

$N_{\rm eff,inj}=\frac{\left[\sum_{j}^{N_{\rm
det}}s_{j}\right]^{2}}{\left[\sum_{j}^{N_{\rm det}}s_{j}^{2}-N_{\rm
gen}^{-1}(\sum_{j}^{N_{\rm det}}s_{j})^{2}\right]}.$ (8)

A typical value for numerical stability is $N_{\rm eff,inj}>4N_{\rm obs}$.

To summarise:

  * •

icarogw computes numerically the hierarchical likelihood in Eq. (1) as:

$\ln[\mathcal{L}(\\{x\\}|\Lambda)]\approx-\frac{T_{\rm obs}}{N_{\rm
gen}}\sum_{j=1}^{N_{\rm det}}s_{j}+\sum_{i}^{N_{\rm obs}}\ln\left[\frac{T_{\rm
obs}}{N_{{\rm s},i}}\sum_{j=1}^{N_{{\rm s},i}}w_{i,j}\right]\,.$ (9)

  * •

In the case that a scale-free likelihood is used, the hierarchical likelihood
in Eq. (2) is:

$\ln[\mathcal{L}(\\{x\\}|\Lambda)]\approx-N_{\rm obs}\ln\left[\frac{1}{N_{\rm
gen}}\sum_{j=1}^{N_{\rm det}}s_{j}\right]+\sum_{i}^{N_{\rm
obs}}\ln\left[\frac{1}{N_{{\rm s},i}}\sum_{j=1}^{N_{{\rm
s},i}}w_{i,j}\right]\,.$ (10)

  * •

For each population model, icarogw calculates two numerical stability
estimators, the effective number of posterior samples for each event in Eq.
(4) and the effective number of injections in Eq. (8). If at least one of
these numerical estimators is below the threshold set by the user, icarogw
returns a $\ln[\mathcal{L}(\\{x\\}|\Lambda)]=-\infty$. This prevents the
population model to be chosen.

##  3 Structure of the code

icarogw contains several python modules used for population inference. icarogw
modules are divided in core modules and support modules. The core modules are
required to calculate the hierarchical likelihood and are necessary even when
using icarogw for a non-GW application (see Sec. 3.1 for an example). The
support modules provide utility functions that could be required by the core
modules. Most of the support modules contain functions used for GW
applications. In Fig. 1 we display a schematic view of icarogw core modules in
red boxes and the support modules in the yellow box.

Figure 1: icarogw modules structure. Core modules used for the general
hierarchical Bayesian inference are colored in red, while support modules
mostly used for GW applications are colored in yellow. The green boxes
highlight the computation of the rate weights which are needed for the
calculation of the hierarchical likelihood given in Eq. (9).

Here we provide a general description of icarogw’s modules, while we refer the
reader to Sec. 3.1 and Sec. 4 for more practical examples. The wrappers.py
module contains python classes for the events production rate $\frac{{\rm
d}N}{{\rm d}t{\rm d}\theta}$. Each rate class should specify what are the
event-level parameters $\theta$ used to calculate $\frac{{\rm d}N}{{\rm
d}t{\rm d}\theta}$ and what are the population-level parameters $\Lambda$ used
to calculate the rate model. Each rate class should also contain instructions
on how to update the rate model from the population parameters $\Lambda$. The
posterior_samples.py module contains classes that are used to allocate PE
samples that represent the possible true values of the event parameters. The
injections.py module contains a class used to allocate the injections that are
used to evaluate the selection bias. Both the single-event PE samples and
injections should be provided with the priors samples used to generate them.
The likelihood.py module contains the hierarchical likelihood class used by
the python package bilby (Ashton et al., 2019b, a; Ashton & Talbot, 2021). The
other modules are used by the core modules as utilities. The icarogw’s general
logic works as follows:

  * 1./

Injections to evaluate selection biases and PE samples from the observed
events are allocated in the classes present in injections.py and
posterior_samples.py.

  * 2./

Injections and PE samples are passed to the population models from the
wrappers.py module to evaluate the $w_{i,j}$ and $s_{j}$ rate coefficients
defined in Eq. (3) and Eq. (7). Support modules can be used if necessary,
_e.g._ we might want to use cosmological models or a particular set of
probabilistic distributions. The most invoked support routines come from the
priors.py module (which contains probabilistic models) and the cosmology.py
module which contains routines for cosmological calculations.

  * 3./

The rate coefficients $w_{i,j}$ and $s_{j}$ are passed in likelihood.py to
calculate the hierarchical likelihood in Eq. (9) or (10).

###  3.1 A custom example: Estimation of the mean of a gaussian process in
presence of selection biases

In this section, we show how to work with icarogw using a single custom
example. We will discuss how to generate mock data consistently with the
statistical model and how to use icarogw to infer population properties of
observations subject to a selection bias. Below we will describe step-by-step
how to set up the problem, generate the mock data, and perform the
hierarchical Bayesian inference with icarogw.

##### The problem

A gaussian process with mean $\mu_{*}=0$ and variance $\sigma^{2}_{*}=1$
generates events $x$ with a rate $R_{x}$ of 100 samples per second. The sample
rate per unit time can be written as:

$\frac{{\rm d}N}{{\rm d}x{\rm
d}t}=R_{x}\mathcal{N}(x|\mu_{*}=0,\sigma^{2}_{*}=1),$ (11)

where $\mathcal{N}$ indicates the normal distribution. Samples of $x$ are
recorded by an experimental apparatus that is contaminated by a noise process.
Due to the presence of noise, the detector does not record directly $x$ values
but records $y$ values. The likelihood of obtaining $y$ given $x$ is given by:

$\mathcal{L}_{\rm noise}(y|x)=\mathcal{N}(y|\mu=x,\sigma^{2}_{n}=1).$ (12)

Moreover, the detector is able to only record values of $y\in[-1,1]$, which is
an equivalent of the selection bias.

The user’s task is to estimate the mean $\mu_{*}$ of the original gaussian
process, assuming $\sigma_{n}$ is known, from the observed values of
$\\{y\\}=\\{y_{1},y_{2},\ldots,y_{N_{\rm obs}}\\}$.

##### Generation of the observations:

The first step is to generate a set of observations, or events $\\{y\\}$ to
which we want to apply the hierarchical inference. For each detected $y_{i}$,
we will also need to generate a set of PE samples $x_{i,j}$ that we believe
could be the true value of $x_{i}$ from which $y_{i}$ was generated. This
simulation step is required to be statistically consistent to not obtain
biased results.

⬇

1import numpy as np

2import bilby

3T_obs = 1 # Observation time in seconds

4Rx = 100 # Rate of samples per second

5x_true = np.random.randn(int(Rx*T_obs)) # Generates values of x according to
the rate

6

7y_measured = x_true + np.random.randn(len(x_true)) # Generate a measure of y

8idx = (y_measured>=-1) & (y_measured<=1) # Find detected samples

9y_measured = y_measured[idx] # Apply selection bias

10x_PE = []

11for i in range(len(y_measured)):

12 gsampler = bilby.core.prior.Normal(mu=y_measured[i], sigma=1.)

13 x_PE.append(gsampler.sample(5000)) # Generate gaussian distributed
posterior samples

Listing 1: This python code shows how to generate mock observations $\\{y\\}$
with inferred PE samples on the true values $x_{i}$ consistently with the
statistical model and presence of selection biases.

The Lst. 1 shows how to generate a set of observed events $\\{y\\}$ with
associated posterior samples of the true values $\\{x\\}$. Given an observed
value $y_{i}$, and a uniform prior $\pi_{\rm PE}(x)\propto{\rm constant}$, we
can generate posterior samples on $x$. According to the noise likelihood model
in Eq. (12) the PE samples on $x$ are distributed according to:

$p(x|y_{i})\propto\mathcal{L}_{\rm noise}(y_{i}|x)\pi_{\rm PE}(x),$ (13)

that is a gaussian distribution centered in $y_{i}$ with $\sigma^{2}_{n}=1$.
In Fig. 2, left panel, we display the histogram of the events’ true $x$ values
generated and detected in 1 second of observation. In one second of
observation, we generate 100 events and we are able to detect only 56 events
due to the presence of the selection bias. It is important to notice that we
are actually able to detect events $x$ outside the detection area of
$y\in[-1,1]$. This is because our selection threshold is on $y$ and not on $x$
itself. In other words, the noise process in Eq. (12) can generate a sample of
$y$ in the detection area even if $x\notin[-1,1]$.

Figure 2: Left panel: Histogram of the 100 generated (blue line) and 56
detected (orange line) samples in one second of data taking. Right panel:
Histogram of the $10^{5}$ generated (blue line) and 16763 detected (orange
line) injections used to evaluate the selection bias. The gray areas in the
plots display the detection region identified in the $y$ space.

##### Generation of injections to evaluate the selection bias:

We now need to generate a set of injections to evaluate the selection bias,
namely Eq. (7). The procedure to generate a set of detected injections is
logically equivalent to the procedure of generating a set of observed events.
However, there are some different choices we need to do.

⬇

1Ngen = 100000 # Ngen for the injections

2x_injections = np.random.uniform(-6, 6, Ngen) # Generate the x samples from a
distribution which we think can cover all the detectable cases

3y_measured = x_injections + np.random.randn(len(x_injections)) # Generate a
measured value of y

4idx = (y_measured>-1) & (y_measured<1) # Check what we can detect

5y_measured = y_measured[idx]

Listing 2: python code showing how to generate a set of injections to evaluate
selection biases for the gaussian custom problem.

First, the injection set should be generated from a prior $\pi_{\rm inj}(x)$
wide enough so that we are sure that our injection set captures _all_ the
values of $x$ that we might be able to detect even with noise fluctuations.
The code written in Lst. 2 shows how to generate the injections set. The logic
is the following: a set of $N_{\rm gen}$ injections (events) is generated from
a flat prior $\pi_{\rm inj}(x)$. For each injection, a value of $y_{i}$ is
generated according to the noise likelihood model, then the $N_{\rm det}$
injections with $y_{i}\in[-1,1]$ are saved. The right panel in Fig. 2 depicts
this first requirement. The original simulated injections are distributed in
$x$ according to a uniform distribution between $[-6,6]$, namely $\pi_{\rm
inj}(x)=1/12$. We can see that the detected injections are distributed with a
different profile than the original one (this is due to noise fluctuations in
Eq. (12)), and most importantly they smoothly transition to 0 (not detected)
around $x=\pm 4$. The smooth transition to 0 is indicative of the fact that
our prior range $\pi_{\rm inj}(x)$ is wide enough to capture all the
observable $x$. The list of detected $\\{x\\}$ injections are saved together
with the prior $\pi_{\rm inj}(x)$ used for their generation i.e. the values of
the prior evaluated at the detected $\\{x\\}$.

##### Set up and run of icarogw analysis:

We now have all the material, namely a set of observations $\\{y\\}$ with
posterior samples $x_{i,j}$ and a set of injections for selection bias, needed
to run icarogw. The only remaining item is the rate model. icarogw does not
contain a rate model in wrappers.py for this custom problem, so we will need
to create a new one. Rate models in icarogw are standardized, thus we will use
this example to see how you can build a custom rate model in icarogw.

We need to create a python class for the rate model that corresponds to the
following rate:

$\frac{{\rm d}N}{{\rm d}x{\rm
d}t}=R_{x}\mathcal{N}(x|\mu_{*},\sigma^{2}_{*}=1).$ (14)

The code snippet in Lst. 3 shows how to create an icarogw rate for this
problem.

⬇

1class my_gaussian_rate(object):

2 ’’’

3 A class for the rate of the gaussian example

4 ’’’

5 def __init__(self, scale_free=False):

6 self.scale_free = scale_free # The class needs to have a scale_free flag

7 if scale_free:

8 self.population_parameters = [’mu_star’] # Population parameters if we want
the scale-free version

9 else:

10 self.population_parameters = [’mu_star’, ’Rx’] # Population parameters

11 self.event_parameters = [’x’] # Events parameters, in this case we just
have x

12

13 def update(self, **kwargs): # We need a method that updates the population
model

14 self.gmod = bilby.core.prior.Normal(mu=kwargs[’mu_star’], sigma=1.) # Save
a gaussian

15 if not self.scale_free:

16 self.Rx = kwargs[’Rx’] # Save Rx

17

18 def log_rate_PE(self, prior, **kwargs): # Tell us how to calculate the log
of the rate for the PE samples

19 # prior is the prior that you applied to generate the PE samples, written
in the variable x

20 log_weights = self.gmod.ln_prob(kwargs[’x’]) - np.log(prior)

21 if not self.scale_free:

22 log_out = log_weights + np.log(self.Rx)

23 else:

24 log_out = log_weights

25 return log_out

26

27 def log_rate_injections(self,prior, **kwargs): # Tells us how to calculate
the log of the rate for the injections

28 # prior is the prior that you applied to generate injections, written in
the variable x

29 return self.log_rate_PE(prior, **kwargs)

30

31myrate = my_gaussian_rate() # Initialize the rate model

Listing 3: python example on how to code in icarogw a rate model for a custom
problem.

The rate requires a  init() method, where the user must specify the population
parameters (the rate $R_{x}$ and mean $\mu_{*}$ in this example) and the
single-level event variables ($x$) for which we calculate the rate. The user
should also specify if the rate is scale-free or not: in this case, and in
general, $R_{x}$ is not a variable that we consider. The rate should also have
an  update() method, that specifies how to update the rate model with a new
set of population parameters $R_{x},\mu_{*}$. In addition, the  log_rate_PE()
and  log_rate_injections() methods tell respectively how to calculate the
$\log w_{i,j}$ for the posterior samples and the $\log s_{j}$ for the
injections. Note that both the  log_rate_PE() and  log_rate_injections()
methods calculate the rate coefficients using the single event parameters (in
our example, $x$), and they both require the priors used to generate the PE
and injections $\pi_{\rm PE/inj}(x)$ written in terms of $x$. Also note that
for this custom problem,  log_rate_PE() and  log_rate_injections() perfectly
coincide as it should be according to the theoretical basis in Sec. 2.
However, for some particular applications (see Sec. 4) we might want to treat
injections slightly differently from PE samples for numerical reasons.

At this point, we can start preparing the inputs for the icarogw analysis. The
first step is to allocate the observed events and their PE samples on $x$ in
the  posterior_samples_catalog() class from the posterior_samples.py module.
This can be done as in the code snippet below (Lst. 4).

⬇

1 from icarogw.posterior_samples import posterior_samples,
posterior_samples_catalog

2# Icarogw wants you to store all the PE samples in a dict of
posterior_samples classes

3

4posterior_dict = {}

5for i in range(len(x_PE)): # Loop over all the events that you detect

6 # You need to provide the PE samples in a dict, IMPORTANT field of the dict,
it is equal to the variable needed to evaluate the rate, i.e. x.

7 posterior_dict[str(i)] = posterior_samples({’x’: x_PE[i]},

8 prior=np.ones(5000)) # PE samples have been generated with uniform prior in
x, so here we pass a constant. NOT IMPORTANT if not normalized

9posterior_dict = posterior_samples_catalog(posterior_dict)

Listing 4: python example on how to allocate the observed events and their PE
samples in the  posterior_samples_catalog() class of icarogw.

It is important to notice that we also need to pass the prior $\pi_{\rm
PE}(x)$ used to generate the PE samples. In this case, since we have used a
constant prior to generate the PE samples, we are going to pass constant
values (Note: It is not important that the PE prior is correctly normalized in
$x$). We can now allocate the injections in the  injections() class from the
injections.py module. This can be done with the code snippet in Lst. 5.

⬇

1# Icarogw wants the injections to evaluate the selection bias to be stored

2injections = icarogw.injections.injections({’x’: x_injections}, # Same logic
as the posterior samples class

3 prior = np.ones_like(x_injections)/12, # Prior used for the injections,
IMPORTANT need to be correctly normalized

4 ntotal = Ngen, # Number of generated injections

5 Tobs = 1.) # Observation time in units compatible with Rx)

6injections.update_cut(idx) # We apply the cut on the detected injections
stored above

7myrate.update(Rx=100., mu_star=0.) # We update the rate model, note that this
is the true model

8injections.update_weights(myrate) # We update the weights of the injections

9Nexp = injections.expected_number_detections() # We calculate the expected
number of detections

10

11print(’Example, for the true model you would expect {:f} detections in 1
second’.format(Nexp))

12print(’The effective number of injections used to calculate this number is
{:f}’.format(injections.effective_injections_number()))

Listing 5: python example on how to allocate the injections for evaluation of
selection biases in the  injections class.()

The definition of the injection class is very similar to the definition of the
PE samples class. There are, however, a few extra inputs and steps that need
to be given and done. The injection class requires the prior $\pi_{\rm
inj}(x)$ to be correctly normalized. Moreover, it requires the observation
time (same units of the rate class, in this case, seconds) and the total
number of injections generated (even the non-detected ones). Moreover, we need
to tell the injection class which injections have been detected. This can be
done with the  update_cut() method. The code in Lst. 5 also shows how the
injections can be used with the rate model to calculate the expected number of
detections $N_{\rm exp}$.

We are now ready to start the icarogw analysis. With the rate model handled by
the class  my_gaussian_rate(), the posterior samples of the observed events,
and the injection set for selection biases, we can initialize the hierarchical
likelihood class from the likelihood.py module. The Lst. 6 shows how to define
the hierarchical likelihood in icarogw.

⬇

1from icarogw.likelihood import hierarchical_likelihood

2# We define the hierarchical likelihood. This is a classical Bilby likelihood

3likelihood = hierarchical_likelihood(posterior_dict, #Posterior samples dict

4 injections, # Injections class

5 myrate, # Rate model

6 nparallel=100, # Number of PE to use per event

7 neffPE=10, # Number of effective PE for numerical stability

8 neffINJ=None) # Number of effective injections

9# Priors

10prior = {’mu_star’: bilby.core.prior.Uniform(-2, 2),

11 ’Rx’: bilby.core.prior.LogUniform(1e-1, 1e3)}

12

13# Run bilby

14myres = bilby.run_sampler(likelihood, priors=prior, save=False)

Listing 6: python example showing how to initialize in icarogw an hierarchical
likelihood class. The bottom line of the code snippet shows how to run the
nested sampling algorithm dynesty from bilby.

The hierarchical likelihood takes as input the PE samples and injections
classes together with the rate class. Moreover, the hierarchical likelihood
requires some technical flags:  nparallel is the number of posterior samples
that will be used per event to compute $w_{i,j}$,  neffPE is the effective
number of PE samples per event in Eq. (4) and  neffINJ the effective number of
injections required by the hierarchical likelihood. If  neffINJ=None, a
default threshold of  neffINJ=4*Nobs will be used. The last line in Lst. 6
runs the dynesty nested sampling algorithm from the bilby package (Ashton et
al., 2019b; Ashton & Talbot, 2021) to generate posterior samples on the
population parameters. The results are going to be stored in the myres object.

Figure 3:  Joint and marginal posterior distributions on the population
parameters $R_{x}$ and $\mu_{*}$ for the custom example. The shaded areas in
the 2-D panel correspond to the $68.3\%$ and $95\%$ credible intervals, while
the marginalised 1-D panels show the the $68.3\%$ credible intervals. The
black dashed lines mark the true values used for the simulations.

Fig. 3 shows the reconstructed posterior distribution on the population
parameters. As we can see from the 2-D plot, the population parameters are
included in the 2$\sigma$ credible intervals.

##  4 Application to a compact binary coalescence case

Despite that icarogw could 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. As such, icarogw contains already
coded CBC merger rate models in the module wrappers.py. In this section, we
show several examples of how to use icarogw for CBC and cosmological
population inference.

From GW detections we estimate the luminosity distance of the source $d_{L}$,
the detector masses $m_{1,d},m_{2,d}$, the sky position $\Omega$, and also
parameters related to the spins of the compact objects (we indicate them
generally with $\chi$ although two parameterizations are available, see App.
B.3). All the rate models implemented in icarogw work with
$d_{L},m_{1,d},m_{2,d}$, and $\chi$, namely with a rate:

$\frac{{\rm d}N}{{\rm d}d_{L}{\rm d}\Omega{\rm d}m_{1,d}{\rm d}m_{2,d}{\rm
d}\chi{\rm d}t_{d}},$ (15)

where $t_{d}$ is the detector unit time. Eq. (15) is the CBC detector rate
model. However, in GW astronomy we are mostly interested in inferring
properties for source masses $m_{1/2,s}$ and redshift $z$, namely we want to
infer a CBC source rate model. As a reminder, the source masses are related to
the detector masses from the relation $m_{1/2,d}=m_{1/2,s}(1+z)$, and the
luminosity distance is related to the redshift by the choice of a cosmological
model and possible GR deviations at cosmological scales. So, without loss of
generality, icarogw works with a detector rate written in terms of “detector
frame” variables $\theta_{\rm D}$ and aims at inferring properties in “source
frame” variables $\theta_{\rm S}$. In icarogw, the detection rate is written
in terms of source rate as:

$\frac{{\rm d}N}{{\rm d}\theta_{\rm D}{\rm d}t_{d}}=\frac{{\rm d}N}{{\rm
d}\theta_{\rm S}{\rm d}t_{s}}\frac{{\rm d}t_{s}}{{\rm
d}t_{d}}\frac{1}{\det{J_{\rm D\rightarrow S}}}=\frac{{\rm d}N}{{\rm
d}\theta_{\rm S}{\rm d}t_{s}}\frac{1}{1+z}\frac{1}{\det{J_{\rm D\rightarrow
S}}}.$ (16)

In the equation above, the factor $1/1+z$ comes from the difference between
source-frame and detector-frame time, $\det{J_{\rm D\rightarrow S}}$ is the
determinant of the Jacobian from the change of variables $\theta_{\rm
D}\rightarrow\theta_{\rm S}$. In the next sections we will see how the
detector rate model is parameterized for the different types of analyses
available in icarogw.

Let us clarify an important aspect. The right-hand side of Eq. (16) is
evaluated with the following procedure.

  * •

Given a certain cosmological model and/or GR deviation model, we convert
detector-frame quantities to source-frame quantities, _i.e._ $\theta_{\rm
D}\rightarrow\theta_{\rm S}$.

  * •

Using the source-frame quantities, we calculate the source rate model and the
Jacobian. Both the Jacobian and source rate model, therefore, depends
typically on the chosen cosmological model which is not fixed but varies as we
explore the parameter space.

  * •

We then calculate the detector rate model using Eq. (16) which is used for the
hierarchical inference.

  * •

All the previous actions are done in the rate classes available in the
wrappers.py module, more specifically in the methods  log_rate_PE() and
log_rate_injections().

###  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 to infer 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
$\theta_{\rm D}$ are $(d_{L},m_{1,d},m_{2,d},\chi)$, _i.e._ luminosity
distance, detector masses, and spin parameters. The source event parameters
are $\theta_{\rm S}=(z,m_{1,s},m_{2,s},\chi)$, _i.e._ the redshift, two source
masses, and the spin parameters. The injections and posterior samples classes
_must_ be passed with values (and associated priors) written in terms of
$d_{L},m_{1,d},m_{2,d}$, and $\chi$. The Jacobian between detector and source
event parameters is:

$\det{J_{\rm D\rightarrow S}}=\frac{\partial d_{L}}{\partial z}(1+z)^{2}.$
(17)

In Eq. (17), the factor $(1+z)^{2}$ comes from the transformation from source-
frame to detector-frame of both $m_{1}$ and $m_{2}$ (Note : There is no term
associated to the spin part in the Jacobian, since we assume that they are
independent of the redshift, hence they are the same in the source and
detector frame). The expression of the differential of the luminosity distance
can be found in App. A.1 for standard cosmological models and in App. A.2 for
modified gravity models. In the end, the detector rate for the spectral sirens
analysis is parameterized as:

$\frac{{\rm d}N}{{\rm d}d_{L}{\rm d}m_{1,d}{\rm d}m_{2,d}{\rm d}\chi{\rm
d}t_{d}}=R_{0}\Psi(z;\Lambda)p_{\rm pop}(m_{1,s},m_{2,s}|\Lambda)p_{\rm
pop}(\chi|\Lambda)\frac{{\rm d}V_{c}}{{\rm
d}z}\frac{1}{1+z}\frac{1}{\det{J_{\rm D\rightarrow S}}},$ (18)

where $R_{0}$ is the CBC merger rate per comoving volume per year (in $\rm
Gpc^{-3}yr^{-1}$), $\Psi(z;\Lambda)$ is a function parametrizing the rate
evolution in redshift such that $\Psi(z=0;\Lambda)=1$, $p_{\rm
pop}(m_{1,s},m_{2,s}|\Lambda)$ is a prior distribution describing the
production of source masses, and $p_{\rm pop}(\chi|\Lambda)$ is a prior
distribution for the production of spin parameters. Finally, $V_{c}$ is the
comoving volume. Note that all these quantities depend on a given set of
population parameters $\Lambda$.

The rate in Eq. (18) is handled by the  CBC_vanilla_rate() class in the
wrappers.py module. An example of how to initialize this rate model is
provided in Lst. 7.

⬇

1import icarogw

2cw = icarogw.wrappers.FlatLambdaCDM_wrap(zmax=10.) #Initialize a wrapper for
the cosmological model

3mw = icarogw.wrappers.massprior_PowerLawPeak() #Initialize a Power Law plus
peak model

4rw = icarogw.wrappers.rateevolution_Madau() # Initialize a Madau model for
the rate

5rate_model = icarogw.wrappers.CBC_vanilla_rate(cw, mw, rw, scale_free=True) #
Collate everything in one single rate, note the flag scale_free that tells we
don’t want to use R0

6print(’The population parameters are ’, rate_model.population_parameters)

7print(’The event parameters you need to pass are ’,
rate_model.event_parameters)

8

9’’’

10The population parameters are [’H0’, ’Om0’, ’alpha’, ’beta’, ’mmin’, ’mmax’,
’delta_m’, ’mu_g’, ’sigma_g’, ’lambda_peak’, ’gamma’, ’kappa’, ’zp’]

11The event parameters you need to pass are [’mass_1’, ’mass_2’,
’luminosity_distance’]

12’’’

Listing 7: python code example of how to initialize a CBC rate model for
spectral sirens.

To be initialized, the spectral siren rate model requires: (i) a cosmological
model with optional GR deviations, (ii) a model for the source mass
distribution $p_{\rm pop}(m_{1,s},m_{2,s}|\Lambda)$, and (iii) a model for the
rate evolution function $\Psi(z;\Lambda)$. Optionally we can pass a model for
the spin parameters $p_{\rm pop}(\chi|\Lambda)$. The list of population,
cosmology, and GR deviation models available in icarogw is provided in App. B.
For the example in Lst. 7 we initialized a source rate model composed by (i) a
standard flat $\Lambda$CDM cosmology, (ii) source masses distributed according
to a power law + peak model, and (iii) a Madau-Dickinson like rate evolution
model. No spin model was used for this example.

In analogy to the custom example presented in Sec. 3.1, now we need to load
the PE samples and injections that are required for the analysis. Regarding
the injections needed to evaluate the selection biases, we will use a set of
injections publicly released by the LVK collaboration (Aasi et al., 2015;
Acernese et al., 2015; Akutsu et al., 2021) which were used in (Abbott et al.,
2021a, a, 2023b) and are also available in Abbott et al. (2021c) and the data
distribution related to that paper. As PE samples, we will be using the 42 BBH
events used in (Abbott et al., 2021a) which have been collected in a pickle
file for the data release of that paper.

⬇

1import pickle

2data=pickle.load(open(’../injections/GWTC3_cosmo_paper_injections.p’, ’rb’))
# Injections set used for the O3 cosmology paper

3time_O3 = data[’Tobs’] # Extract observation time in yr

4prior = data[’prior’] # Extract the prior; the prior is already written in
detector frame masses and luminosity distance [Mpc]

5injections_dict = {’mass_1’: data[’mass_1’],

6 ’mass_2’: data[’mass_2’],

7 ’luminosity_distance’:data[’distance’]} # Initialize the dictionary

8

9# Here we initialize the injection set

10inj = icarogw.injections.injections(injections_dict,

11 prior=prior,

12 ntotal=data[’ntotal’], #Total number of injections generated, even the
missed one

13 Tobs=time_O3)

14inj.update_cut(data[’snr’]>=11.) # Update the injection set

15

16posterior_dict = {}

17poss = pickle.load(open(’../GW_events/GWTC3_BBH_SNR_11_IFAR_4.p’, ’rb’)) #
Save the posterior samples of GW events from the O3 cosmology paper together

18for ev in poss.keys():

19 pos_dict = {’mass_1’: poss[ev][’mass_1_det’],

20 ’mass_2’: poss[ev][’mass_2_det’],

21 ’luminosity_distance’: poss[ev][’distance’]} # Extract detector frame
masses and luminosity distance

22 posterior_dict[ev] = icarogw.posterior_samples.posterior_samples(pos_dict,

23 prior=np.power(poss[ev][’distance’], 2.)) # Only a dl2 prior was applied on
PE

24posterior_dict =
icarogw.posterior_samples.posterior_samples_catalog(posterior_dict)

Listing 8: python example showing how to load the PE samples and injections
used for the spectral siren analysis.

The Lst. 8 loads the set of PE samples and injections and allocates them in
the  posterior_samples_catalog() and the  injections() classes. In the code
example, the PE samples and injections are both passed alongside their
generation priors written in terms of luminosity distance and masses.

⬇

1from icarogw.likelihood import hierarchical_likelihood

2H0array = np.linspace(20, 120, 200)

3

4# We initialize the likelihood passing everything that is needed

5likelihood = hierarchical_likelihood(posterior_dict, inj, rate_model,
nparallel=2048, # PE to use for each event

6 neffINJ=None, # Effective number of injections will be set to 4 N_obs

7 neffPE=20) # Number of PE to use

8

9posterior = np.zeros_like(H0array)

10# On grid inference

11for i, H0 in enumerate(H0array):

12 likelihood.parameters = {’Om0’: 0.308, ’alpha’: 3.78, ’beta’: 0.81, ’mmin’:
4.98, ’mmax’: 112.5, ’delta_m’: 4.8, ’mu_g’: 32.27, ’sigma_g’: 3.88,
’lambda_peak’: 0.03, ’gamma’: 4.59, ’kappa’: 2.86, ’zp’: 2.47, ’H0’: H0}

13 posterior[i] = likelihood.log_likelihood()

14

15posterior -= posterior.max()

16posterior = np.exp(posterior)

17posterior = posterior/np.trapz(posterior, H0array)

18

19# MCMC inference

20priors = {’Om0’: 0.308, ’alpha’: 3.78, ’beta’: 0.81, ’mmin’: 4.98, ’mmax’:
112.5, ’delta_m’: 4.8, ’mu_g’: 32.27, ’sigma_g’: 3.88, ’lambda_peak’: 0.03,
’gamma’: 4.59, ’kappa’: 2.86, ’zp’: 2.47, ’H0’: bilby.prior.Uniform(20, 140)}

21my_res = bilby.run_sampler(likelihood, priors=priors, save=False,
sampler=’emcee’, nsteps=600, nwalkers=4)

Listing 9: python example showing how to calculate the hierarchical Hubble
constant posterior on a grid and how to calculate it with emcee sampling.

Now we can initialize the hierarchical likelihood using the PE samples,
injections, and rate models. The Lst. 9 shows how to initialize and use the
hierarchical likelihood for this example. In the code above, we calculate the
posterior on $H_{0}$ using two methods. In the first one, we fix the
population parameters and calculate the posterior on a $H_{0}$ grid, while in
the second one we run a Markov Chain Monte Carlo algorithm to sample $H_{0}$.
The first method is not computationally doable when we want to estimate other
population parameters at the same time. In this example, we expect the two
methods to provide the same result. In Fig 4 (left plot), we show the $H_{0}$
posterior computed with these two methods.

Figure 4: Left plot: Posterior probability density distribution on the Hubble
constant $H_{0}$ using 42 spectral sirens from GWTC-3 and icarogw. Right plot:
Posterior on the modified gravity parameter $\Xi_{0}$ using the same set of
events and injections.

icarogw can also consider modified gravity models (see App. A.2 for more
details). Typically, in these modified gravity models the GW luminosity
distance is different from the standard luminosity distance. For instance, the
$\Xi_{0}$ model (Belgacem et al., 2019) prescripts the following
parametrisation:

$d_{L}^{\mathrm{GW}}=d_{L}^{\mathrm{EM}}\left(\Xi_{0}+\frac{1-\Xi_{0}}{(1+z)^{n}}\right).$
(19)

The logic to infer modification of gravity with icarogw is always the same, we
need to define a detector rate model. This can be done with just the third
line of Lst. 10. The rest of the example code calculates the posterior on
$\Xi_{0}$ using the grid and MCMC methods. The PE samples and injections
inputs required for icarogw are the same as in Lst. 8.

⬇

1import icarogw

2# This is the only difference, the cosmology wrapper needs to be modified by
a modified gravity wrapper.

3cw =
icarogw.wrappers.Xi0_mod_wrap(icarogw.wrappers.FlatLambdaCDM_wrap(zmax=20.))

4

5mw = icarogw.wrappers.massprior_PowerLawPeak()

6rw = icarogw.wrappers.rateevolution_Madau()

7rate_model = icarogw.wrappers.CBC_vanilla_rate(cw, mw, rw, scale_free=True)

8

9Xi0array = np.linspace(0.5, 5, 200)

10likelihood = icarogw.likelihood.hierarchical_likelihood(posterior_dict, inj,
rate_model, nparallel=2048, neffINJ=None, neffPE=20)

11

12posterior = np.zeros_like(Xi0array)

13effn = np.zeros_like(Xi0array)

14effPE = np.zeros_like(Xi0array)

15

16for i, Xi0 in enumerate(Xi0array):

17 likelihood.parameters={’Om0’: 0.308, ’alpha’: 3.78, ’beta’: 0.81, ’mmin’:
4.98, ’mmax’: 112.5, ’delta_m’: 4.8, ’mu_g’: 32.27, ’sigma_g’: 3.88,
’lambda_peak’: 0.03, ’gamma’: 4.59, ’kappa’: 2.86, ’zp’: 2.47, ’H0’: 67.7,
’Xi0’: Xi0, ’n’: 2.}

18 posterior[i] = likelihood.log_likelihood()

19 effn[i] = likelihood.injections.effective_injections_number() # Injections
are updated in the likelihood class

20 effPE[i] =
likelihood.posterior_samples_dict.get_effective_number_of_PE().min()

21

22posterior -= posterior.max()

23posterior = np.exp(posterior)

24posterior = posterior/np.trapz(posterior, Xi0array)

25

26priors = {’Om0’: 0.308, ’alpha’: 3.78, ’beta’: 0.81, ’mmin’: 4.98, ’mmax’:
112.5, ’delta_m’: 4.8, ’mu_g’: 32.27, ’sigma_g’: 3.88, ’lambda_peak’: 0.03,
’gamma’: 4.59, ’kappa’: 2.86, ’zp’: 2.47, ’H0’: 67.7, ’Xi0’:
bilby.prior.Uniform(0.5, 5), ’n’: 2.}

27my_res = bilby.run_sampler(likelihood, priors=priors, save=False,
sampler=’emcee’, nsteps=600, nwalkers=4, outdir=’out2’)

Listing 10: python example showing how to initialize and run an analysis for
the modified gravity model $\Xi_{0}$. The modified gravity model is loaded on
top of the standard cosmological background model.

The only difference with the standard case is that we need to change the
cosmology background model to a modified gravity model. Fig. 4 (right plot)
shows a comparison of the $\Xi_{0}$ posterior obtained from this example.

Finally, as mentioned in Sec. 1, icarogw also offers some functionalities that
return numerical stability estimators such as the effective number of PE
samples per event and the effective number of injections for a given rate
model, that are not specific to modified gravity models. The methods
get_effective_number_of_PE() from the  posterior_samples_catalog() class and
effective_injections_number() from the  injections() class can return these
two quantities after that the classes are updated with the rate model. It is
important that the rate model is updated with its  update() method each time a
new value of $\Xi_{0}$ for this specific case is used. In Fig. 5 we display
the effective number of PE samples and the effective number of injections as a
function of $\Xi_{0}$. As we can see from the plots, the numerical stability
estimators are always above the threshold for numerical stability.

Figure 5: Left plot: Minimum over events of the effective number of PE samples
(each event has an effective number of PE samples) computed with Eq. (4) as a
function of $\Xi_{0}$. Right plot: Effective number of injections computed
with Eq. (8) as a function of $\Xi_{0}$. In both plots the black dashed lines
indicate the thresholds required by icarogw for numerical stability.

###  4.2 Galaxy catalog merger rates

The “galaxy catalog” analysis adds information on the GW event redshift from
galaxy surveys (Schutz, 1986; Del Pozzo, 2012; Gray et al., 2020, 2022). A
detailed description of the method employed in icarogw is given in
Mastrogiovanni (2023). Also in 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, this time with extra information
coming from the galaxy catalog.

In this method, the detector event parameters $\theta_{\rm
D}=(d_{L},m_{1,d},m_{2,d},\Omega,\chi)$ 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 $\theta_{\rm
S}=(z,m_{1,s},m_{2,s},\Omega,\chi)$, that is, the redshift and the two source
masses, sky direction, and spins. The Jacobian between detector and source
event parameters is still given by Eq. (17).

The detector rate for the galaxy catalog analysis is parameterized as:

$\displaystyle\frac{{\rm d}N}{{\rm d}d_{L}{\rm d}m_{1,d}{\rm d}m_{2,d}{\rm
d}\Omega{\rm d}\chi{\rm d}t_{d}}$ $\displaystyle=R^{*}_{\rm
gal,0}\Psi(z;\Lambda)p_{\rm pop}(m_{1,s},m_{2,s}|\Lambda)p_{\rm
pop}(\chi|\Lambda)\frac{1}{1+z}\frac{1}{\det{J_{\rm D\rightarrow S}}}\times$
$\displaystyle\times\left[\frac{dV_{c}}{dzd\Omega}\phi_{*}(H_{0})\Gamma_{\rm
inc}(\alpha+\epsilon+1,x_{\rm max}(M_{\rm thr}),x_{\rm
min})+\sum_{j=1}^{N_{\rm
gal}(\Omega)}\frac{f_{L}(M(m_{j},z);\epsilon)p(z|z^{j}_{\rm
obs},\sigma^{j}_{\rm z,obs})}{\Delta\Omega}\right],$

where $R^{*}_{\rm gal,0}$ is the local CBC merger rate per galaxy per year (in
$\rm 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 (2023) for more details). The first term is the completeness
correction, _i.e._ it accounts for missing galaxies. It depends on the
absolute magnitude threshold of galaxy detection $M_{\rm thr}$, on how likely
more luminous galaxies can emit GW events (through the $\epsilon$ parameter),
and on the Schecter luminosity function and its parameters $\phi_{*}$ and
$\alpha$, $x_{\rm min/max}$ are defined in (Mastrogiovanni, 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 $f_{L}(M(m_{j},z);\epsilon)$ quantifies how likely
luminous galaxies emit GW events, while $p(z|z^{j}_{\rm obs},\sigma^{j}_{\rm
z,obs})$ is the probability of having a certain value of $z$ given observed
values of galaxy redshift inside the catalog.

The rate in Eq. (18) is handled by the  CBC_catalog_vanilla_rate() class in
the wrappers.py module. An example of how to initialize this rate model is
provided in Lst. 11. The initialization of the merger rate for the catalog
analysis is similar to the spectral siren case, but it requires few additional
inputs.

⬇

1from icarogw.wrappers import CBC_catalog_vanilla_rate

2

3# Wrappers definition

4cosmo_wrap =
icarogw.wrappers.Xi0_mod_wrap(icarogw.wrappers.FlatLambdaCDM_wrap(zmax=20.))

5mass_wrap = icarogw.wrappers.massprior_PowerLawPeak()

6rate_wrap = icarogw.wrappers.rateevolution_Madau()

7

8# Rate definition

9rate_model = CBC_catalog_vanilla_rate(cat, cosmo_wrap,

10 mass_wrap, rate_wrap,

11 average=True, # This flag tells you that you want to use a sky-averaged
detection probability to evaluate selection biases

12 scale_free=True)

Listing 11: python example showing how to initialize a CBC merger model that
accounts for a galaxy catalog.

The  cat() class is the precompiled galaxy catalog (see later), while the
average flag can be set to True if the user wishes to calculate $N_{\rm exp}$
using a sky-averaged111The sky averaged detection probability in icarogw is
computed by defining a source rate averaged over the sky:  $\frac{{\rm
d}N_{\rm CBC}}{{\rm d}z{\rm d}\theta_{\rm S}{\rm d}t_{s}}=N_{\rm
pix}\Delta\Omega\left\langle\frac{{\rm d}N_{\rm CBC}}{{\rm d}z{\rm
d}\Omega{\rm d}\theta_{\rm S}{\rm
d}t_{s}}(\Lambda)\right\rangle_{\rm\Omega}=4\pi\left\langle\frac{{\rm d}N_{\rm
CBC}}{{\rm d}z{\rm d}\Omega{\rm d}\theta_{\rm S}{\rm
d}t_{s}}(\Lambda)\right\rangle_{\rm\Omega}.$ (21) galaxy number density rather
than sky-depend galaxy number density. Typically, we set  average to True
since there are only a few injections that fall in a given sky pixel. The code
snippet in Lst. 12 shows how to construct a compiled galaxy catalog for
icarogw.

⬇

1nside = 64 # Same nside as the posterior

2

3# We need to build the galaxy catalog with a reference cosmology

4# The only cosmological parameter that impacts the catalog construction is
Om0, no problem for H0

5cosmo_ref = icarogw.cosmology.astropycosmology(zmax=10.)

6cosmo_ref.build_cosmology(FlatLambdaCDM(H0=67.7, Om0=0.308))

7

8

9data=h5py.File(’../icarogwCAT/glade+.hdf5’) # Load the glade+ galaxy catalog

10

11cat = icarogw.catalog.galaxy_catalog() # Initialize class

12

13# Below we extract information from glade and save everything in a
dictionary to pass to icarogw

14cat_data={}

15for key in [’ra’, ’dec’, ’z’, ’sigmaz’, ’m_K’]:

16 if key == ’m_K’:

17 cat_data[’m’] = data[key][:]

18 else:

19 cat_data[key] = data[key][:]

20

21# Create the HDF5 file

22cat.create_hdf5(filename=’glade+_kband_BBH.hdf5’,

23 cat_data=cat_data,

24 band=’K’, # Select band

25 nside=nside) # Nside

26

27cat.calculate_mthr(50) # Calculate the apparent magnitude threshold using
the 50% percentile (median) of galaxies in each bin

28

29# Calculate the interpolant, it will take a long time

30cat.calc_dN_by_dzdOmega_interpolant(cosmo_ref, # Reference cosmology

31 1., # Luminosity weight

32 Nintegration=20, # Bin resolution for each gaussian redshift localization

33 zcut=0.5, # Where to cut the catalog, completeness will be 0 after that
value

34 type=’gaussian’, # Type of redshift localization, options available:
‘uniform‘, ‘gaussian‘

35 Numsigma=3.) # Number of sigmas for the redshift uncertainties.

Listing 12: python example showing how to create a galaxy catalog class
starting from a list of galaxy right ascensions, declinations, redshifts,
apparent magnitudes, and redshift uncertainties.

The galaxy catalog is created with the following steps:

  * 1./

We allocate the  galaxy_catalog() class from the catalog.py module.

  * 2./

The galaxy catalog requires inputting a list of galaxies with right ascension
and declination (in radians), redshift, redshift uncertainties, and apparent
magnitudes. The  create_hdf5() method creates an hdf5 file where it will
allocate all the data required by icarogw.

  * 3./

The  calculate_mthr() method calculates the apparent magnitude threshold in
each sky pixel. The apparent magnitude threshold is defined as the apparent
magnitude percentile (defined by the user) of all the galaxies reported in
each pixel.

  * 4./

The method  calc_dN_by_dzdOmega_interpolant() creates an interpolant for the
catalog-based galaxy number density in Eq. (LABEL:eq:ratecat) (second term in
the square brackets). The interpolant is used to quickly evaluate the galaxy
number density as a function of sky position and redshift. The creation of the
interpolant requires several technical flags which are described in the code
snippet.

  * 5./

The pre-compiled galaxy catalog is stored in the created hdf5 and can be
loaded for later use by the method  load_hdf5() without the need to regenerate
it.

The hierarchical analysis for the galaxy catalog can be run in the exact same
way of the spectral siren case (see Lst. 9). The only difference is that when
loading the PE samples and the injection set (see Lst 8), the user must recall
to give to the PE samples and injections classes the right ascension (flag
right_ascension) and declination (flag  declination) in radians. The user must
also correctly normalize the prior to take into account the event sky area.
Typically an isotropic prior is used to create PE and injections; the user
must remember to add a prior term $\pi_{\rm inj,PE}(\Omega)=1/4\pi$. Finally,
the user should remember to divide in sky pixels the PE samples and
injections, using the same resolution used for the galaxy catalog. This can be
done by invoking the  pixelize() method from the PE samples and injection
classes.

Besides the usual hierarchical Bayesian inference, when working with galaxy
catalogs the user might want to check the galaxy density profile in the sky
localization area of the GW event. This quantity plays a fundamental role in
Eq. (LABEL:eq:ratecat) when estimating possible rates for the GW event. The
galaxy_catalog() class method  effective_galaxy_number_interpolant() returns
the galaxy density profile (and its completeness correction) for a given sky
position and redshift. This function can be used to generate figures such as
Fig. 6, where it is possible to visualize the galaxy number density profile in
the area of a given GW event.

Figure 6: Top row: Effective number of CBC emitters per redshift and steradian
in-catalog (solid lines) and out-catalog (dashed line). The first column is
generated using a CBC emission rate independent on galaxy luminosity
($\epsilon=0$) while the second is generated using a rate linearly
proportional to the galaxy luminosity ($\epsilon=1$). Second row: CBC
completeness calculated. The shaded areas indicate the contours identified by
the 90% credible sky area of GW190814 while the lines correspond to the median
values. Reproduction from Mastrogiovanni (2023).

###  4.3 Multi-messenger sources merger rates

In icarogw, there is also the implementation of a third methodology based on
GW events with an associated EM counterpart. The new rate model now takes into
account additional constraints on sky position and redshift from the EM
counterpart. In this case, the hierarchical likelihood needs to be modified
(let us use Eq. (2)) to include the information on the EM counterpart. By
assuming that the GW measure is independent of the EM measure, the overall
likelihood term is now $\mathcal{L}_{\rm
EM+GW}(x_{i}|z,\Omega,m_{1,s},m_{2,s})$, which describes the measure of
$z,\Omega,m_{1,s},m_{2,s}$ from EM and GW data. Here we assume that the EM
data measures $\Omega$ and $z$, while the GW data can measure
$\Omega,z,m_{1,s},m_{2,s}$ independently, so that:

$\mathcal{L}_{\rm
EM+GW}(x_{i}|z,\Omega,m_{1,s},m_{2,s},\chi)\propto\mathcal{L}_{\rm
EM}(x_{i}|z,\Omega)\mathcal{L}_{\rm GW}(x_{i}|z,\Omega,m_{1,s},m_{2,s},\chi).$
(22)

The integral of the numerator in Eq. (2) now becomes:

$I=\int\mathcal{L}_{\rm EM}(x_{i}|z,\Omega)\mathcal{L}_{\rm
GW}(x_{i}|z,\Omega,m_{1,s},m_{2,s},\chi)\frac{{\rm d}N}{{\rm d}z{\rm
d}\Omega{\rm d}m_{1,s}{\rm d}m_{2,s}{\rm d}\chi{\rm d}t_{s}}\frac{1}{1+z}{\rm
d}m_{1,s},m_{2,s}{\rm d}\chi{\rm d}z{\rm d}\Omega.$ (23)

To perform this integral, icarogw defines a function of redshift $F(z)$
starting from a collection of GW PE samples. The function is:

$F(z)=\left[\frac{1}{N_{s}^{\rm EM}}\sum_{i}^{N_{s}^{\rm EM}}w_{i}\right]{\rm
KDE}[z_{i},{\rm weights}=w_{i}],$ (24)

where KDE is a kernel density estimate performed on the GW redshift samples
$z_{i}$ with weights $w_{i}$. The weights are given by:

$w_{i}=\frac{1}{\pi_{\rm EM}(z^{i})\pi_{\rm
PE}(z^{i},\@vec{m}^{i},\@vec{\chi}^{i})}\frac{{\rm d}N}{{\rm d}z{\rm
d}m_{1,s}{\rm d}m_{2,s}{\rm d}\chi{\rm
d}t_{s}}\biggr{|}_{i}\frac{1}{1+z^{i}},$ (25)

where $\pi_{\rm EM}(z)$ is the prior used by the EM experiment to provide a
redshift measure. The integral $I$ is evaluated by summing over posterior
samples of $z$ given by the EM counterpart (that lives in a very narrow $z$
region is compared to where $F(z)$ is defined). Numerically:

$I(z)\approx\frac{1}{N_{\rm s,EM}}\sum_{i}^{N_{\rm sEM}}F(z_{i}).$ (26)

icarogw assumes that the detection probability is dominated $p_{\rm
det}^{EM+GW}(\cdot)=p_{\rm det}^{GW}(\cdot)$, _i.e._ if we are able to detect
the GW, then we are able to detect for sure the EM counterpart. We made this
assumption as EM selection biases could be very model dependent.

⬇

1 from astropy import constants

2

3# Load GW170817 data as usual, we need information on the sky location of
course

4

5GW170817 = h5py.File(’../GW_events/GW170817.hdf5’)

6ppd = {’mass_1’:
GW170817[’IMRPhenomPv2NRT_lowSpin_posterior’][’m1_detector_frame_Msun’],

7
’mass_2’:GW170817[’IMRPhenomPv2NRT_lowSpin_posterior’][’m2_detector_frame_Msun’],

8
’luminosity_distance’:GW170817[’IMRPhenomPv2NRT_lowSpin_posterior’][’luminosity_distance_Mpc’],

9
’right_ascension’:GW170817[’IMRPhenomPv2NRT_lowSpin_posterior’][’right_ascension’],

10 ’declination’:GW170817[’IMRPhenomPv2NRT_lowSpin_posterior’][’declination’]}

11

12# Initialize the posterior dict

13posterior_dict = {’GW170817’:
icarogw.posterior_samples.posterior_samples(ppd,

14
prior=np.power(GW170817[’IMRPhenomPv2NRT_lowSpin_posterior’][’luminosity_distance_Mpc’],
2.))}

15

16# Pixelize GW170817

17posterior_dict[’GW170817’].pixelize(64)

18

19# Now, we add an EM counterpart to GW170817. To add an EM counterpart you
need to know perfectly RA and DEC in radians

20# And you can pass posterior samples on the redshift inferred from EM

21

22# Values taken from Nature volume 551 (2017). These need to be recessional
velocities due to the Hubble flow.

23zcenter = 3017/constants.c.to(’km/s’).value

24sigma = 166/constants.c.to(’km/s’).value

25

26# Add the EM counterpart. The code will also tell you how many PE you have
in the pixel of the EM counterpart

27posterior_dict[’GW170817’].add_counterpart(np.random.randn(10000)*sigma +
zcenter, ppd[’right_ascension’][0], ppd[’declination’][0])

28

29posterior_dict =
icarogw.posterior_samples.posterior_samples_catalog(posterior_dict)

Listing 13: python example showing how to add an EM counterpart to the
posterior events class.

Lst. 13 shows an example of how to prepare the data inputs for the EM
counterpart method. In order to use the electromagnetic counterpart method, we
need to load and pixelize a set of PE samples. Then we can invoke the
add_counterpart() method to provide a set of possible cosmological redshifts
that were inferred from the EM counterpart signal. The  add_counterpart()
method will select the PE samples falling in the sky direction of the EM
counterpart and will print a message on how many PE samples fall in that area.
The definition of the detector rate, the loading of injection set, and the
running of the analysis is equivalent with respect to what was discussed in
Lst. 8 and Lst. 9.

##  5 Sanity checks for CBC populations

We have performed several tests reproducing population results which are
consistent with previous analyses.

##### Spectral siren analyses

: We used the same 42 BBHs and injection set used in Abbott et al. (2021a) to
reproduce several populations, cosmological, and beyond-GR analyses generated
with icarogw 1.0. We reproduced constraints on the mass and redshift
distribution, together with the cosmological parameters $H_{0},\Omega_{m}$,
and $w_{0}$ obtained in Abbott et al. (2021a). Fig. 7 shows the results of one
of our tests for some of the population parameters. The posterior probability
density distribution was obtained by icarogw 2.0 in comparison to the ones
obtained by icarogw 1.0. The two posterior distributions are in perfect
agreement with each other. We also reproduced the constraints on modified
gravity models and the population of BBHs generated with icarogw 1.0 in Leyde
et al. (2022). Fig. 8 the results of the test for the running Planck mass-
modified gravity model (see Sec. A.2 for more details) for a subset of the
population parameters. Also in this case the posteriors generated with icarogw
1.0 and icarogw 2.0 are in a perfect match between each other. Finally, we
replicated the constraints on the default and gaussian spin models obtained by
Abbott et al. (2023a) to validate the new spin implementation in icarogw.
Using 60 GW events from the third observing run, we jointly inferred the spin
parameters with other population parameter. For this analysis, the
cosmological parameters were fixed to the Astropy Planck15 values, in order to
replicate the same setup as in Abbott et al. (2023a). We found that both the
Default spin model and the gaussian spin model results are in very good
agreement, as shown in Fig. 7 and Fig. 8.

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. (2021) (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).

Figure 7: Comparison of the posteriors on the mass, redshift, and cosmological
background population parameters obtained by icarogw 1.0 (blue line) and
icarogw 2.0 (green line) using 42 BBHs. We used a powerlaw+peak model for the
source masses, a Madau-Dickinson-like model for the merger rate, and a flat
$w$CDM cosmology. The plot labels are set as the code flag names present in
icarogw 2.0, see App. A.1-B.1-B.2, for more details on population models and
the code flags. Figure 8: Comparison of the posteriors on the mass, redshift,
cosmological background, and beyond-GR population parameters obtained by
icarogw 1.0 (blue line) and icarogw 2.0 (green line) using 42 BBHs. We used a
powerlaw+peak model for the source masses, a Madau-Dickinson-like model for
the merger rate, and a flat $w$CDM cosmology. The plot labels are set as the
code flag names present in icarogw 2.0, see App. A.1-A.2-B.1-B.2, for more
details on population models and the code flags.

##### Galaxy catalog analysis:

We tested icarogw against the results generated by gwcosmo (Gray et al., 2020,
2022) in Abbott et al. (2021a) using the glade+ (Dálya et al., 2022) galaxy
catalog with the infrared K-band. For this particular type of test, we have
fixed the source mass and redshift population model to the same used by
gwcosmo in Abbott et al. (2021a). We run only the inference for the Hubble
constant $H_{0}$. Results are shown in Fig. 9. We find 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., 2021a). However, 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
obtain posteriors that are in agreement. We have also tested the catalog
analysis on GW190814, one of the best localized and close-by dark sirens. For
this event, we find results in perfect agreement in the case that the Galaxy
catalog is not used. However, we find a slightly different $H_{0}$ posterior
(more “peaked” in the $H_{0}$ tension region), with respect to the one
obtained by gwcosmo in Abbott et al. (2021a). In Mastrogiovanni (2023) we
describe in more detail this comparison.

Figure 9: Posterior probability density distributions obtained by icarogw 2.0
(blue line) from 42 BBHs used in Abbott et al. (2021a) in comparison with
gwcosmo (orange line).

We note that the posteriors we obtain for the catalog analysis are also
consistent with the ones generated in Finke et al. (2021) from the code
DarkSirensStat. Note, however, that these results are generated with a
different choice of BBHs and galaxy catalog descriptions with respect to the
analyses performed in Abbott et al. (2021a).

##### Electromagnetic counterpart analysis:

To test the electromagnetic counterpart method, we infer the Hubble constant
using the BNS merger GW170817 and its EM counterpart. We use the low spin
IMRPhenom PE samples from Abbott et al. (2021b) and we assume a Hubble flow
recessional velocity of $v_{H}=3017$ km/s with uncertainty $\sigma_{v}=166$
km/s. As an injection set, we used the BNS injection set released for O3
sensitivity in (Abbott et al., 2023a).

The analysis done by the code gwcosmo in (Abbott et al., 2017) used the
following assumptions to describe the CBC merger rate. First, the luminosity
distance was approximated using linear cosmology. This has a set of
consequences, namely:

$d_{L}=d_{c}=\frac{cz}{H_{0}},$ (27)

thus implying the following relations for the comoving volume:

$\displaystyle\frac{\partial d_{L}}{\partial z}$ $\displaystyle=$
$\displaystyle\frac{c}{H_{0}},$ (28) $\displaystyle\frac{\partial z}{\partial
d_{L}}$ $\displaystyle=$ $\displaystyle\frac{H_{0}}{c},$ (29)
$\displaystyle\frac{\partial V_{c}}{\partial z}$ $\displaystyle=$
$\displaystyle 4\pi\frac{c^{3}z^{2}}{H_{0}^{3}}.$ (30)

The CBC merger rate model for GW170817 used in the analysis was:

$p(m_{1,d},m_{2,d})=\frac{\Theta(m_{2,d}<m_{1,d})}{2(m_{\rm d,max}-m_{\rm
d,min})^{2}},$ (31)

where the $\Theta$ function ensures that the detector secondary mass is
lighter than the primary one. Moreover, the analysis neglected the $1/1+z$
factor coming from the difference between source and detector frames. The
overall merger rate was:

$\frac{{\rm d}N_{\rm CBC}}{{\rm d}t_{d}{\rm d}d_{L}{\rm d}m_{1}{\rm
d}m_{2}}=\frac{{\rm d}N_{\rm CBC}}{{\rm d}t_{d}{\rm
d}d_{L}}p(m_{1,d},m_{2,d})=R_{0}\frac{\partial V_{c}}{\partial
z}\frac{\partial z}{\partial
d_{L}}p(m_{1,d},m_{2,d})=R_{0}4\pi\frac{c^{2}z^{2}}{H_{0}^{2}}\frac{\Theta(m_{2,d}<m_{1,d})}{2(m_{\rm
d,max}-m_{\rm d,min})^{2}}.$ (32)

We remark that the aforementioned assumptions on cosmology and rate model are
not expected to provide a noticeable difference when calculating the weights
$w_{i,j}$ for GW170817 PE samples. This is because GW170817 is a very close-by
GW event, and even for extreme values of $H_{0}$, 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 in the calculation of the selection bias. With O3 sensitivities, BNSs
are detected up to a luminosity distance of $\sim 300-400$ Mpc, where the
linear cosmology approximation can fail (especially if a high value of $H_{0}$
is chosen). So, even when reproducing GW170817 it is important to consider the
rate model assumed in Abbott et al. (2021d). Lst. 14 shows how to define the
aforementioned rate model in icarogw consistently with the EM method presented
in Sec. 4.

⬇

1SPEED_L = constants.c.to(’km/s’).value

2class my_170817_rate(object):

3 ’’’

4 Note that in the rate below we do not put any mass prior term because it is
constant as we said.

5 ’’’

6

7 def __init__(self):

8 self.scale_free = True

9 self.population_parameters = [’H0’] # The only pop parameter is H0

10 self.event_parameters = [’luminosity_distance’, ’z_EM’]

11 def update(self, **kwargs):

12 self.H0 = kwargs[’H0’] # Update H0

13 def log_rate_PE(self, prior, **kwargs):

14 if len(kwargs[’luminosity_distance’].shape) != 2:

15 raise ValueError(’The EM counterpart rate wants N_ev x N_samples arrays’)

16

17 z = kwargs[’luminosity_distance’]*self.H0/SPEED_L # Calculation of redshift

18

19 # Differential of the comoving volume

20 log_dVc_dz = xp.log(4*xp.pi) + 3*xp.log(SPEED_L) + 2*xp.log(z) -
3*xp.log(self.H0)

21

22 # Compute the weights. The first two terms are the rate term, the second
removes the prior used to generate PE as usual

23 log_weights = log_dVc_dz - xp.log(SPEED_L/self.H0) - xp.log(prior)

24

25 n_ev = kwargs[’luminosity_distance’].shape[0]

26 lwtot = xp.empty(kwargs[’z_EM’].shape)

27 for i in range(n_ev):

28 ww = xp.exp(log_weights[i,:])

29 kde_fit = gaussian_kde(z[i,:], weights=ww/ww.sum())

30 lwtot[i,:] = logsumexp(log_weights[i,:]) -
xp.log(kwargs[’luminosity_distance’].shape[1]) +
kde_fit.logpdf(kwargs[’z_EM’][i,:])

31 log_out = lwtot

32 return log_out

33

34 def log_rate_injections(self, prior, **kwargs):

35

36 # Same rate model for the injections. Just remember that here we neglect
the EM seleciton bias.

37 z = kwargs[’luminosity_distance’]*self.H0/SPEED_L

38 log_dVc_dz = xp.log(4*xp.pi) + 3*xp.log(SPEED_L) + 2*xp.log(z) -
3*xp.log(self.H0)

39 # Sum over posterior samples

40 log_weights = log_dVc_dz - xp.log(prior) - xp.log(SPEED_L/self.H0)

41 log_out = log_weights

42 return log_out

Listing 14: python snippet showing how to define a new rate model for
GW170817.

Fig. 10 shows the posterior that we obtain with icarogw 2.0 and the method
highlighted in Sec. 4 for GW170817 in comparison with gwcosmo. The posteriors
are in good agreement with each other.

Figure 10: Posterior distributions for $H_{0}$ obtained from GW170817 with
icarogw 2.0 (blue line) in comparison with gwcosmo (orange dotted line). The
posterior of gwcosmo is taken from Abbott et al. (2021b).

##  6 Conclusions and future development

In this paper we have presented icarogw 2.0, a python software for population
properties inference in presence of selection biases. We have provided several
working examples that are available in a GitHub repository . We described 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). Moreover, the checks on
the spin parameters estimation also coincide with previous studies with
studies present in the literature. At last, the galaxy catalog and the EM
counterpart methods agree with previous studies.

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 correlation 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
repository.

## Software packages

icarogw uses the public python packages astropy (Astropy Collaboration et al.,
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).

###### 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. RG
was supported by ERC starting grant SHADE 949572 and STFC grant ST/V005634/1.

##  Appendix A Cosmological and GR deviation models

Cosmological models and GR deviation models are handled by the classes and
functions in cosmology.py and are organized in high-level wrappers for quick
use in the wrappers.py module. They can be passed to
CBC_vanilla_EM_counterpart(),  CBC_vanilla_rate(), and
CBC_catalog_vanilla_rate() to construct an overall merger rate.

In Tab. LABEL:tab:cosmomod we provide an overview of all the cosmological and
GR models available. 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 next
sections) while leaving the other cosmological quantities unchanged.

Table 1: Summary table for all the background cosmology and models available in wrappers.py of icarogw. More details on the models can be found in Sec. A.1. Model name | Wrapper name | Population parameters
---|---|---
Symbol | Code flag | Description
Flat $\Lambda$CDM | FlatLambdaCDM_wrap() | $H_{0}$ |  H_0 | Hubble constant in [km/s/Mpc]
$\Omega_{m}$ |  Om0 | Matter energy density
Flat $w$CDM | FlatwCDM_wrap() | $H_{0}$ |  H_0 | Hubble constant in [km/s/Mpc]
$\Omega_{m}$ |  Om0 | Matter energy density
$w_{0}$ |  w0 | Dark Energy equation of state parameter
Table 2: Summary table for all the background cosmology and beyond-GR models in wrappers.py of icarogw. More details on the models can be found in Sec. A.2. Model name | Wrapper name | Population parameters
---|---|---
Symbol | Code flag | Description
$\Xi_{0}$ model | Xi0_mod_wrap() | $\Xi_{0}$ | Xi_0 |  See Eq. (39)
$n$ | n |  See Eq. (39)
Running planck mass | cM_mod_wrap() | $c_{M}$ | c_M |  See Eq. (48)
Extra dimensions | extraD_mod_wrap() | $D$ | D |  # spacetime dimensions, See Eq. (44)
$n$ | n |  Scaling parameter, See Eq. (44)
$R_{c}$ | Rc | GR deviation scale in Mpc
$\alpha$-log | alphalog_mod_wrap() | $\alpha_{1}$ | alphalog_1 |  See Eq. (41)
$\alpha_{2}$ | alphalog_2 |  See Eq. (41)
$\alpha_{3}$ | alphalog_3 |  See Eq. (41)

###  A.1 Cosmological background models

In principle icarogw is able 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:

$d_{L}=\frac{c(1+z)}{H_{0}}\int_{0}^{z}\frac{{\rm
d}z^{\prime}}{E(z^{\prime})},$ (33)

where $H(z)=H_{0}E(z)$, which is the same as the EM luminosity distance
assuming GR, as in this section. The differential of the luminosity distance
is:

$\frac{\partial d_{L}}{\partial
z}=\frac{d_{L}(z)}{1+z}+\frac{c(1+z)}{H_{0}}\frac{1}{E(z)}.$ (34)

The comoving volume is:

$V_{c}=\int_{0}^{z}{\rm d}\Omega{\rm d}z^{\prime}\frac{{\rm d}V_{c}}{{\rm
d}z^{\prime}{\rm d}\Omega},$ (35)

and the differential of the comoving volume is:

$\frac{{\rm d}V_{c}}{{\rm d}z{\rm
d}\Omega}=\frac{c^{3}}{H^{3}_{0}}\frac{1}{E(z)}\left[\int_{0}^{z}\frac{{\rm
d}z^{\prime}}{E(z^{\prime})}\right]^{2}.$ (36)

The function $E(z)$ depends on the cosmological model assumed. For the Flat
$\Lambda$CDM model:

$E^{2}(z)=\Omega_{m}(1+z)^{3}+(1-\Omega_{m}),$ (37)

while for the Flat $w_{0}$CDM model:

$E^{2}(z)=\Omega_{m}(1+z)^{3}+(1-\Omega_{m})(1+z)^{3(1+w_{0})}.$ (38)

###  A.2 Beyond-GR models

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

Figure 11: Left panel: Luminosity distance as a function of redshift for the
modified gravity models. Right panel: Differential of the luminosity distance
as a function of redshift for the modified gravity models.

####  A.2.1 The $\Xi_{0}$ model

The luminosity distance is given by (see Eq. 2.31 of Belgacem et al. (2019)):

$d_{L}^{\mathrm{GW}}=d_{L}^{\mathrm{EM}}\left(\Xi_{0}+\frac{1-\Xi_{0}}{(1+z)^{n}}\right)\,.$
(39)

The Jacobian is given by:

$\frac{\mathrm{d}d_{L}^{\mathrm{GW}}}{\mathrm{d}z}=\frac{\mathrm{d}d_{L}^{\mathrm{EM}}}{\mathrm{d}z}\left(\Xi_{0}+\frac{1-\Xi_{0}}{(1+z)^{n}}\right)-d_{L}^{\mathrm{EM}}\frac{n(1-\Xi_{0})}{(1+z)^{n+1}}\,.$
(40)

####  A.2.2 The phenomenological $\log$ parametrization

The luminosity distance is given by:

$d_{L}^{\mathrm{GW}}=d_{L}^{\mathrm{EM}}\left[1+\sum_{\nu=1}^{n=3}\alpha_{\nu}\log^{\nu}(1+z)\right]\,.$
(41)

The Jacobian is given by:

$\frac{\mathrm{d}d_{L}^{\mathrm{GW}}}{\mathrm{d}z}=\frac{\mathrm{d}d_{L}^{\mathrm{EM}}}{\mathrm{d}z}\frac{d_{L}^{\mathrm{GW}}}{d_{L}^{\mathrm{EM}}}+d_{L}^{\mathrm{EM}}\left[\sum_{\nu=1}^{n=3}\alpha_{\nu}\nu\frac{\log^{\nu-1}(1+z)}{1+z}\right]\,.$
(42)

####  A.2.3 Extra-dimensions

In the extra-dimensions model, the luminosity distance is given by (see Eq.
2.22 in Corman et al. (2022)):

$d_{L}^{\rm GW}=d_{L}^{\rm EM}\left[1+\left(\frac{d_{L}^{\rm
EM}}{(1+z)R_{c}}\right)^{n}\right]^{\frac{D-4}{2n}}\,.$ (43)

Let’s define the following function:

$\mathcal{A}=\left[1+\left(\frac{d_{L}^{\rm
EM}}{(1+z)R_{c}}\right)^{n}\right]\,.$ (44)

We also define the exponential:

$\mathcal{E}=\frac{D-4}{2n}\,.$ (45)

We can write the previous equation as:

$\frac{\mathrm{d}d_{L}^{\mathrm{GW}}}{\mathrm{d}z}=(\mathcal{A})^{\mathcal{E}}\left[\frac{\mathrm{d}d_{L}^{\mathrm{EM}}}{\mathrm{d}z}+\frac{n\mathcal{E}}{\mathcal{A}}\left(\frac{d_{L}^{\rm
EM}}{R_{c}}\right)^{n}\left(\frac{\mathrm{d}d_{L}^{\rm
EM}}{\mathrm{d}z}\frac{1}{(1+z)^{n}}-\frac{d_{L}^{\rm
EM}}{(1+z)^{1+n}}\right)\right]\,.$ (46)

####  A.2.4 The $c_{M}$ parametrization

Lastly, we consider a model with a running Planck mass (Lagos et al., 2019):

$d_{L}^{\rm{GW}}=d_{L}^{\rm{EM}}{\rm{exp}}\left[\frac{c_{M}}{2}\int_{0}^{z}\frac{1}{(1+z^{\prime})E^{2}(z^{\prime})}dz^{\prime}\right]\equiv
d_{L}^{\rm{EM}}{\rm{exp}}\left[{\frac{c_{M}}{2}I(z)}\right],$ (47)

which defines $I(z)$. In a flat $\Lambda$CDM model, $I(z)$ can be calculated
analytically and the result is (Eq. 19 in Lagos et al. (2019)):

$d_{L}^{\rm
GW}=d_{L}^{\rm{EM}}{\rm{exp}}\left[\frac{c_{M}}{2\Omega_{\Lambda,0}}\ln\frac{1+z}{(\Omega_{\rm
m,0}(1+z)^{3}+\Omega_{\Lambda,0})^{1/3}}\right]\,,$ (48)

otherwise it needs to be calculated numerically. In any cosmology, the
Jacobian is given by:

$\frac{\mathrm{d}d_{L}^{\mathrm{GW}}}{\mathrm{d}z}=\frac{\mathrm{d}d_{L}^{\mathrm{EM}}}{\mathrm{d}z}\frac{d_{L}^{\mathrm{GW}}}{d_{L}^{\mathrm{EM}}}+d_{L}^{\mathrm{GW}}\cdot\frac{c_{M}}{2}\cdot\frac{1}{(1+z)E^{2}(z)}.$
(49)

##  Appendix B CBC Population models

Population models for the mass, redshift, and spins for CBC are hosted in
wrappers.py and usually make use of probability density distributions defined
in the priors.py module. All the models currently available in icarogw are not
conditionally dependent from each other, _i.e._ the probability distributions
of redshift, source masses, and spins are independent from each other. The
redshift, mass, and spin models provided in icarogw can be passed to the
CBC_vanilla_EM_counterpart(),  CBC_vanilla_rate(), and
CBC_catalog_vanilla_rate() classes to construct a full CBC rate model for
hierarchical inference.

We provide a list of models for the CBC merger rate in redshift in Sec. B.1,
for source masses in Sec. B.2, and for spins in Sec. B.3.

Table 3: Summary table for all the merger rate models available in wrappers.py of icarogw. More details on the models can be found in Sec. B.1. Model name | Wrapper name | Population parameters
---|---|---
Symbol | Code flag | Description
Power Law rate | rateevolution_PowerLaw() | $\gamma$ | gamma | See Eq. (50)
Madau Dickinson | rateevolution_Madau() | $\gamma$ | gamma | See Eq. (51)
$k$ | kappa | See Eq. (51)
$z_{p}$ | zp | Redshift peak of merger rate

(a) (a) Power Law rate model with $\gamma=2$.

(b) (b) Madau-Dickinson rate model with $\gamma=2.7,k=6,z_{p}=0.5$.

Figure 12: Sample of rate models implemented in icarogw. Table 4: Summary table for all the background cosmology and models available in wrappers.py of icarogw. The same mass models are also available for NSBH binaries. For NSBH, the models have an additional string at the end of the name specifying  _NSBH(). NSBH models also include three extra parameters:  mmin_NS,  mmax_NS, and  delta_m_NS 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 Sec. B.2. Model name | Wrapper name | Population parameters
---|---|---
Symbol | Code flag | Description
Power Law | massprior_PowerLaw() | $-\alpha$ | alpha | Power Law index primary mass
$\beta$ | beta | Power Law index secondary mass
$m_{\rm min}$ | mmin |  Minimum source mass [$M_{\odot}$]
$m_{\rm max}$ | max |  Maximum source mass [$M_{\odot}$]
Power Law + Peak | massprior_PowerLawPeak() | $-\alpha$ | alpha | Power Law index primary mass
$\beta$ | beta | Power Law index secondary mass
$m_{\rm min}$ | mmin |  Minimum source mass [$M_{\odot}$]
$m_{\rm max}$ | max |  Maximum source mass [$M_{\odot}$]
$\delta_{\rm m}$ | delta_m |  Smoothing parameter [$M_{\odot}$] Eq. (70)
$\mu_{g}$ | mu_g |  Peak of the gaussian [$M_{\odot}$]
$\sigma_{g}$ | sigma_g |  s.t.d of the gaussian [$M_{\odot}$]
$\lambda_{\rm peak}$ | lambda_peak |  Fraction of events in gaussian $\in[0,1]$
Broken Power Law | massprior_BrokenPowerLaw() | $-\alpha_{1}$ | alpha_1 | First Power Law index primary mass
$-\alpha_{2}$ | alpha_2 | Second Power Law index primary mass
$\beta$ | beta | Power Law index secondary mass
$m_{\rm min}$ | mmin |  Minimum source mass [$M_{\odot}$]
$m_{\rm max}$ | max |  Maximum source mass [$M_{\odot}$]
$\delta_{\rm m}$ | delta_m |  Smoothing parameter [$M_{\odot}$] Eq. (70)
$b$ | b |  Defines $m_{\rm break}=b(m_{\rm max}-m_{\rm min})$
Multi Peak | massprior_MultiPeak() | $-\alpha$ | alpha | Power Law index primary mass
$\beta$ | beta | Power Law index secondary mass
$m_{\rm min}$ | mmin |  Minimum source mass [$M_{\odot}$]
$m_{\rm max}$ | max |  Maximum source mass [$M_{\odot}$]
$\delta_{\rm m}$ | delta_m |  Smoothing parameter [$M_{\odot}$] Eq. (70)
$\mu_{g,\rm low}$ | mu_g_low |  Peak lower mass gaussian [$M_{\odot}$]
$\sigma_{g,\rm low}$ | sigma_g_low |  s.t.d. of the lower mass gaussian in [$M_{\odot}$]
$\mu_{g,\rm high}$ | mu_g_high |  Peak higher mass gaussian [$M_{\odot}$]
$\sigma_{g,\rm high}$ | sigma_g_high |  s.t.d higher mass gaussian in [$M_{\odot}$]
$\lambda_{\rm g}$ | lambda_g |  Events in gaussians $\in[0,1]$
$\lambda_{\rm g,low}$ | lambda_g_low |  Events in lower gaussian $\in[0,1]$
Table 5: Summary table for all the merger spin models available in wrappers.py of icarogw. More details on the models can be found in Sec. B.3. Model name | Wrapper name | Population parameters
---|---|---
Symbol | Code flag | Description
Default | spinprior_default() | $\alpha_{\chi}$ |  alpha_chi |  $\beta$-distribution parameter, see Eq. (89)
$\beta_{\chi}$ |  beta_chi |  $\beta$-distribution parameter, see Eq. (90)
$\sigma_{t}$ |  sigma_t |  Aligned spins dispersion, see Eq. (91)
$\xi$ |  csi_spin | Fraction of events with semi-aligned spins
gaussian | spinprior_gaussian() | $\mu_{\xi,{\rm eff}}$ |  mu_chi_eff |  Mean of $\chi_{\rm eff}$
$\mu_{\xi,{\rm p}}$ |  mu_chi_p |  Mean of $\chi_{\rm p}$
$\sigma_{\xi,{\rm eff}}$ |  sigma_chi_eff |  s.t.d of $\chi_{\rm eff}$
$\sigma_{\xi,{\rm p}}$ |  sigma_chi_p |  s.t.d of $\chi_{\rm p}$
$\rho$ |  rho |  $\chi_{\rm eff}-\chi_{\rm p}$ correlation term

###  B.1 CBC redshift rate evolution models

icarogw contains two models for the redshift evolution of the merger rate, see
Eqs. (18)-(LABEL:eq:ratecat). Table LABEL:tab:ratemod summarises the merger
rate redshift models, while Fig. 12 provides some examples of the models for
specific values of the parameters.

####  B.1.1 Power Law

The rate is parametrized as:

$\psi(z;\gamma)=(1+z)^{\gamma}\,.$ (50)

####  B.1.2 Madau-Dickinson

The rate is parametrized following the star formation rate of Madau &
Dickinson (2014) as:

$\psi(z;\gamma)=[1+(1+z_{p})^{-\gamma-k}]\frac{(1+z)^{\gamma}}{1+\left(\frac{1+z}{1+z_{p}}\right)^{\gamma+k}}\,.$
(51)

###  B.2 Mass models

Most of the mass models are composed of gaussian and power law distributions
which we report in the following. The simple truncated power law distribution
is given by

$\displaystyle\mathcal{P}(x|a,b,\alpha)=\bigg{\\{}\begin{array}[]{lr}\frac{1}{N}x^{\alpha},&\quad\left(a<x<b\right)\\\
0,&\quad\left({\rm{otherwise}}\right)\end{array}$ (54)

and the normalization factor is given by:

$\displaystyle
N=\bigg{\\{}\begin{array}[]{lr}\dfrac{1}{\alpha+1}[b^{\alpha+1}-a^{\alpha+1}],&\quad{\rm{if\,\alpha\neq-1}}\\\
{\ln}{\dfrac{b}{a}},&\quad\rm{if\,\alpha=-1}\end{array}$ (57)

The truncated gaussian distribution is given by:

$\displaystyle\mathcal{G}_{[a,b]}(x|\mu,\sigma)=\bigg{\\{}\begin{array}[]{lr}\dfrac{1}{N}\dfrac{1}{\sigma\sqrt{2\pi}}e^{-\frac{(x-\mu)^{2}}{2\sigma^{2}}},&\quad
a<x<b\\\ 0,&\quad\text{otherwise}\end{array}$ (60)

The normalization factor is expressed through the error function. Then the
normalization factor is:

$N=\int_{a}^{b}\frac{1}{\sigma\sqrt{2\pi}}e^{-\frac{(x-\mu)^{2}}{2\sigma^{2}}}dx=\int_{(a-\mu)/(\sigma\sqrt{2})}^{(b-\mu)/(\sigma\sqrt{2})}\frac{1}{\sqrt{\pi}}e^{-t^{2}}dt.$
(61)

Using the symmetry of the integrand around $x=\mu$ ($t=0$) and the definition
of erf function222 Scipy function.:

$\mathrm{erf}(z)=\frac{2}{\sqrt{\pi}}\int_{0}^{z}e^{-t^{2}}dt,$ (62)

it follows that:

$N=\frac{1}{2}\left({\rm{erf}}\left[\frac{b-\mu}{\sigma\sqrt{2}}\right]-{\rm{erf}}\left[\frac{a-\mu}{\sigma\sqrt{2}}\right]\right).$
(63)

In icarogw, we factorize the prior on mass as:

$\pi(m_{1,s},m_{2,s}|\Lambda)=\pi(m_{1,s}|\Lambda)\pi(m_{2,s}|m_{1,s},\Lambda)\,.$
(64)

When dealing with a NSBH, the neutron star is assigned to $m_{2,s}$ and the
distribution of $m_{2,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=m_{\rm min}$:

$\pi(m_{1,s},m_{2,s}|\Lambda)=[\pi(m_{1,s}|\Lambda)\pi(m_{2,s}|m_{1,s},\Lambda)]S(m_{1}|m_{\rm
min},\delta_{m})S(m_{2}|m_{\rm min},\delta_{m}),$ (65)

where $S$ is a sigmoid-like window function (Eqs. B6-B7 of Abbott et al.
(2020)):

$\displaystyle S(m_{1,s}|m_{\rm
min},\delta_{m})=\bigg{\\{}\begin{array}[]{lr}0,&\left(m<m_{\rm min}\right)\\\
\left[f(m-m_{\rm min},\delta_{m})+1\right]^{-1},&\quad\left(m_{\rm min}\leq
m<m_{\rm min}+\delta_{m}\right)\\\ 1,&\quad\left(m\geq m_{\rm
min}+\delta_{m}\right)\end{array}$ (69)

with

$f(m^{\prime},\delta_{m})=\exp\left(\frac{\delta_{m}}{m^{\prime}}+\frac{\delta_{m}}{m^{\prime}-\delta_{m}}\right).$
(70)

When we apply this window, the priors are numerically renormalized.

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

(a) (a) Truncated Power Law model with parameters $\alpha=0$,$\beta=1$,$m_{\rm
min}=5M_{\odot}$,$m_{\rm max}=100M_{\odot}$.

(b) (b) Power Law + peak model with parameters $\alpha=2$,$\beta=1$,$m_{\rm
min}=5M_{\odot}$,$m_{\rm max}=100M_{\odot}$,
$\mu_{g}=35M_{\odot},\sigma_{g}=5M_{\odot},\lambda_{p}=0.1$

(c) (c) Broken power law model with parameters $\beta=1$,$m_{\rm
min}=5M_{\odot}$,$m_{\rm
max}=100M_{\odot}$,$\alpha_{1}=0$,$\alpha_{2}=1$,$b=0.5$,$\delta_{m}=5M_{\odot}$.

(d) (d) multi peak model with parameters $\alpha=2,\beta=1,m_{\rm
min}=5M_{\odot},m_{\rm max}=100M_{\odot},\delta_{m}=5M_{\odot},\mu_{\rm
g,low}=9M_{\odot},\sigma_{\rm g,low}=1M_{\odot},\lambda_{\rm
g,low}=0.05,\mu_{\rm g,high}=35M_{\odot},\sigma_{\rm
g,high}=5M_{\odot},\lambda_{g}=0.5$

Figure 13: Sample of mass models implemented in icarogw.

####  B.2.1 Truncated Power-Law

The Truncated Power Law model is given by Eq. (64) with:

$\displaystyle\pi(m_{1,s}|m_{\rm min},m_{\rm max},\alpha)$ $\displaystyle=$
$\displaystyle\mathcal{P}(m_{1,s}|m_{\rm min},m_{\rm max},-\alpha)\,,$ (71)
$\displaystyle\pi(m_{2,s}|m_{\rm min},m_{1,s},\beta)$ $\displaystyle=$
$\displaystyle\mathcal{P}(m_{2,s}|m_{\rm min},m_{1,s},\beta)\,,$ (72)

where $\mathcal{P}$ is defined in Eq. (54).

####  B.2.2 Power-Law + Peak

This model was proposed in Talbot et al. (2019) and it is given by Eq. (65)
with:

$\displaystyle\pi(m_{1,s}|m_{\rm min},m_{\rm max},\alpha)$ $\displaystyle=$
$\displaystyle(1-\lambda)\mathcal{P}(m_{1,s}|m_{\rm min},m_{\rm
max},-\alpha)+\lambda\mathcal{G}(m_{1,s}|\mu_{g},\sigma)\,,\quad(0\leq\lambda\leq
1)$ (73) $\displaystyle\pi(m_{2,s}|m_{\rm min},m_{1,s},\beta)$
$\displaystyle=$ $\displaystyle\mathcal{P}(m_{2,s}|m_{\rm
min},m_{1,s},\beta)\,.$ (74)

####  B.2.3 Broken Power Law

This model is based on Eq. (65) and consists basically of two truncated power-
law distributions attached at the point $b$:

$b=m_{\rm min}+(m_{\rm max}-m_{\rm min})f,$ (75)

where $f$ is a scalar in $[0,1]$, so $b=m_{\mathrm{min}}$ for $f=0$. This
model was proposed in Abbott et al. (2020). The priors are the following:

$\displaystyle\pi(m_{1,s}|m_{\rm min},m_{\rm max},\alpha)$ $\displaystyle=$
$\displaystyle\frac{1}{N}[\mathcal{P}(m_{1,s}|m_{\rm
min},b,-\alpha_{1})+\frac{\mathcal{P}(b|m_{\rm
min},b,-\alpha_{1})}{\mathcal{P}(b|b,m_{\rm
max},-\alpha_{2}))}\mathcal{P}(m_{1,s}|b,m_{\rm max},-\alpha_{2})]\,,$ (76)
$\displaystyle\pi(m_{2,s}|m_{\rm min},m_{1,s},\beta)$ $\displaystyle=$
$\displaystyle\mathcal{P}(m_{2,s}|m_{\rm min},m_{1,s},\beta)\,,$ (77)

where the new normalization factor $N$ here is:

$N=1+\frac{\mathcal{P}(b|m_{\rm min},b,-\alpha_{1})}{\mathcal{P}(b|b,m_{\rm
max},-\alpha_{2}))}.$ (78)

####  B.2.4 Multi-Peak

This model is based on Eq. (65) and consists of one power-law + two gaussian
models with:

$\displaystyle\pi(m_{1,s}|m_{\rm min},m_{\rm max},\alpha)$ $\displaystyle=$
$\displaystyle\Bigl{[}(1-\lambda)\mathcal{P}(m_{1,s}|m_{\rm min},m_{\rm
max},-\alpha)+\lambda\lambda_{\rm low}\mathcal{G}(m_{1,s}|\mu_{g,\rm
low},\sigma_{\rm low})+$ $\displaystyle+\lambda(1-\lambda_{\rm
low})\mathcal{G}(m_{1,s}|\mu_{g,\rm high},\sigma_{\rm high})\Bigr{]}\,,$
$\displaystyle\pi(m_{2,s}|m_{\rm min},m_{1,s},\beta)$ $\displaystyle=$
$\displaystyle\mathcal{P}(m_{2,s}|m_{\rm min},m_{1,s},\beta)\,.$ (80)

Note that since $\mathcal{G}$ and $\mathcal{P}$ are already normalized,
$\pi(m_{1,s}|m_{\rm min},m_{\rm max},\alpha)$ is automatically normalised. The
model was used in Abbott et al. (2020).

###  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. 14, we provide a definition for the spin parameters
typically employed in GW studies.

Figure 14: Representation of the spin components for a compact binary
coalescence.

By definition, the $z$ axis of a binary is aligned to the instantaneous
orbital angular momentum ${\@vec{L}}$. The (normalized) spin amplitudes
$\chi_{1,2}$, defined from the Cartesian components of the spin vectors, are:

$\displaystyle\chi_{1}$ $\displaystyle=$
$\displaystyle\sqrt{s_{1,x}^{2}+s_{1,y}^{2}+s_{1,z}^{2}}\,,$ (81)
$\displaystyle\chi_{2}$ $\displaystyle=$
$\displaystyle\sqrt{s_{2,x}^{2}+s_{2,y}^{2}+s_{2,z}^{2}}.$ (82)

The tilt angles $\theta_{1,2}$ are defined as the angle between the BH spins
and the orbital angular momentum, namely:

$\displaystyle\cos\theta_{1}$ $\displaystyle=$
$\displaystyle\frac{s_{1,z}}{\chi_{1}}\,,$ (83) $\displaystyle\cos\theta_{2}$
$\displaystyle=$ $\displaystyle\frac{s_{2,z}}{\chi_{2}}.$ (84)

The effective spin parameter $\chi_{\rm eff}$ and precession spin parameter
$\chi_{\rm p}$ are defined by (Abbott et al. (2021d)):

$\displaystyle\chi_{\rm eff}$ $\displaystyle=$
$\displaystyle\frac{\chi_{1}\cos\theta_{1}+q\chi_{2}\cos\theta_{2}}{1+q}=\frac{s_{1,z}+qs_{2,z}}{1+q}\,,$
(85) $\displaystyle\chi_{\rm p}$ $\displaystyle=$ $\displaystyle{\rm
max}\left[\chi_{1}\sin\theta_{1};\,\left(\frac{4q+3}{3q+4}\right)q\chi_{2}\sin\theta_{2}\right],$
(86)

where the mass ratio $q$ is:

$q=\frac{m_{2}}{m_{1}}\,,\qquad(q\leq 1).$ (87)

(The factors of $q$ appearing in the expression for $\chi_{\rm p}$ come from
the leading order PN equation for $\dot{\@vec{L}}$. Note that the 4 in-plane
components of the perpendicular components of $\@vec{s}$ have been replaced by
one scalar $\chi_{\rm p}$ which is an averaged quantity, see Eq. (3.1) of
Schmidt et al. (2015)).

The parameter $\chi_{\rm eff}$ accounts for the amount of spin aligned with
the orbital angular momentum, as well as the magnitude of the BH spins. Since
$\chi_{\rm eff}$ is bounded between $[-1,1]$, values close to $\chi_{\rm
eff}=1$ correspond to highly spinning BHs with aligned spins, whereas values
close to $\chi_{\rm eff}=-1$ support highly spinning BHs with anti-aligned
spins, and $\chi_{\rm eff}=0$ is consistent with non spinning BHs. The
precession spin parameter $\chi_{\rm 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. 15.

(a) (a) default spins model with
$\alpha_{\chi}=3,\beta_{\chi}=5,\sigma_{t}=0.3,\xi=0.5$

(b) (b) gaussian spins model with $\mu_{\chi_{\rm eff}}=0.06,\sigma_{\chi_{\rm
eff}}=0.12,\mu_{\chi_{\rm p}}=0.5,\sigma_{\chi_{\rm p}}=0.05,\rho=0.5$.

Figure 15: Sample of spin models implemented in icarogw.

####  B.3.1 Default spin model

The model was used in Abbott et al. (2020) and it is proposed after Talbot et
al. (2019); Wysocki et al. (2019). The model works with spins parameterized
using the two spin magnitudes $\chi_{1},\chi_{2}$ and the two cosine of
inclination angles $\cos t_{1},\cos t_{2}$ with respect to the orbital angular
momentum. Note: 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
$\phi_{1}$ and $\phi_{2}$ are not considered here and supposed uniform.

The population distribution is given by:

$\pi(\chi_{1},\chi_{2},\cos\theta_{1},\cos\theta_{2})={\rm
Beta}(\chi_{1}|\alpha,\beta)\pi(\cos\theta_{1}|\xi,\sigma_{t}){\rm
Beta}(\chi_{2}|\alpha,\beta)\pi(\cos\theta_{2}|\xi,\sigma_{t}),$ (88)

namely it is factored into two parts. Above “Beta” is the beta distribution,
calculated with parameters $\alpha$ and $\beta$ defined by:

$\displaystyle\alpha$ $\displaystyle=$
$\displaystyle\left(\frac{1-\mu_{\chi}}{\sigma^{2}_{\chi}}-\frac{1}{\mu_{\chi}}\right)\mu_{\chi}^{2}\geq
1\,,$ (89) $\displaystyle\beta$ $\displaystyle=$
$\displaystyle\alpha\left(\frac{1}{\mu_{\chi}}-1\right)\geq 1.$ (90)

The condition $(\alpha,\beta)\geq 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 (see Eq. (14) in Abbott et al.
(2021d)):

$\pi(\cos\theta_{1,2}|\zeta,\sigma_{t})=\xi\mathcal{G}_{[-1,1]}(\cos\theta_{1,2}|1,\sigma_{t})+\frac{1-\xi}{2},$
(91)

where $\mathcal{G}_{[-1,1]}(\cos\theta_{i}|1,\sigma_{t})$ (see Eq.(60)) is a
truncated gaussian between $-1$ and $1$ on $\cos\theta_{i}$, with mean $1$ and
standard deviation $\sigma_{t}$.
Note: the parameter $\xi$ 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 $\xi$ is the mixing parameter.

####  B.3.2 Gaussian spin model

The Gaussian spin model seeks to measure the joint distribution of $\chi_{\rm
eff}$ and $\chi_{p}$. It was proposed in Miller et al. (2020) and it depends
on 5 parameters that are $\mu_{\chi_{\rm eff}},\sigma_{\chi_{\rm
eff}},\mu_{\chi_{\rm p}},\sigma_{\chi_{\rm p}}$, and $\rho$. The population
probability on $\chi_{\rm eff},\chi_{\rm p}$ is a bivariate gaussian truncated
between [-1,1] for $\chi_{\rm eff}$ and between $[0,1]$ for $\chi_{\rm p}$.
The covariance of the bivariate gaussian is ${\rm cov}_{[\chi_{\rm
eff},\chi_{\rm p}]}=\rho\sigma_{\chi_{\rm eff}}\sigma_{\chi_{\rm p}}$. In
icarogw, this bivariate gaussian is factorized as:

$\pi(\chi_{\rm eff},\chi_{\rm p}|\mu_{\chi_{\rm eff}},\sigma_{\chi_{\rm
eff}},\mu_{\chi_{\rm p}},\sigma_{\chi_{\rm
p}},\rho)=\mathcal{G}_{[-1,1]}(\chi_{\rm eff}|\mu_{\chi_{\rm
eff}},\sigma_{\chi_{\rm eff}})\mathcal{G}_{[0,1]}(\chi_{\rm
p}|\mu_{*},\sigma_{*}),$ (92)

where

$\displaystyle\mu_{*}$ $\displaystyle=$ $\displaystyle\mu_{\chi_{\rm
p}}+\frac{{\rm cov}_{[\chi_{\rm eff},\chi_{\rm p}]}}{\sigma^{2}_{\chi_{\rm
eff}}}(\chi_{\rm eff}-\mu_{\chi_{\rm eff}})\,,$ (93) $\displaystyle\sigma_{*}$
$\displaystyle=$ $\displaystyle\frac{\sigma_{\chi_{\rm p}}{\rm
cov}_{[\chi_{\rm eff},\chi_{\rm p}]}}{\sigma^{2}_{\chi_{\rm eff}}}\,.$ (94)

This factorization is equivalent to a bivariate gaussian distribution.

Note: Typically GW priors for PE samples are not applied in $\chi_{\rm eff}$
and $\chi_{\rm p}$ but in spin magnitudes and inclination angles. This
includes non-trivial priors and jacobians to account for, for further details
see Callister (2021). In the conversions.py, icarogw implements a code from
the GitHub project “effective-spin-priors” T. Callister  to obtain $\chi_{\rm
eff}$ and $\chi_{\rm p}$ from priors isotropic or aligned in spin directions.

##  Appendix C GPU implementation

icarogw is compatible for parallel computation on GPUs by implementing cupy
(Okuta et al., 2017). The GPU/CPU interaction is handled in the cupy_pal.py
module. The code is able to automatically recognize if a cupy compatible GPU
is available for the system. If the GPU is available, icarogw will
automatically import cupy, while if this is not the case, icarogw will import
numpy. In both cases, when importing icarogw a message will be printed to
indicate what package has been loaded. It is also possible to set icarogw to
load numpy by default. To do so, if the user does not wish to use the GPU, add
a config.py file under the working directory with a global variable
CUPY=False.

In Fig. 16 we show a comparison between the timing of the hierarchical
likelihood in Eq. 1 computation with CPU and GPU. The CPU for the test was an
Intel Core i9-11950H (8 cores HT, 2.6-5.0GHz Turbo) and the GPU was an NVIDIA
GeForce RTX3080 (16Gb GDDR6, 6144 cores CUDA). With less than 100 GW events,
the GPU gains a factor of 3-5 in the computation of the hierarchical
likelihood. With almost a thousand events (about 10 million total PE samples
used), the GPU is able to gain more than an order of magnitude in speed.

Figure 16: Average elapsed time taken by icarogw to evaluate the hierarchical
likelihood in Eq. 1 with CPU (blue diamonds) and GPU (orange filled circles).
For each event we use 2048 PE samples to evaluate the hierarchical likelihood.
Each point is calculated as the median elapsed time over 500 computations of
the hierarchical likelihood computation. CPU: Intel Core i9-11950H (8 cores
HT, 2.6-5.0GHz Turbo). GPU: NVIDIA GeForce RTX3080 (16Gb GDDR6, 6144 cores
CUDA).

## References

  * Aasi et al. (2015) Aasi, J. et al. 2015, Classical and Quantum Gravity, 32, 074001
  * Abbott et al. (2021a) Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2021a, Constraints on the cosmic expansion history from GWTC-3
  * Abbott et al. (2017) Abbott, B. P. et al. 2017, Nature, 551, 85
  * Abbott et al. (2021b) Abbott, B. P. et al. 2021b, Astrophys. J., 909, 218
  * Abbott et al. (2021c) Abbott, B. P. et al. 2021c, Data distribution of Constraints on the cosmic expansion history from the GWTC-3, 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.
  * Abbott et al. (2021d) Abbott, R., Abbott, T., Acernese, F., et al. 2021d, arXiv preprint arXiv:2111.03634
  * Abbott et al. (2023a) Abbott, R., Abbott, T., Acernese, F., et al. 2023a, Physical Review X, 13, 011048
  * Abbott et al. (2021) Abbott, R., Abbott, T. D., Acernese, F., et al. 2021, arXiv, arXiv:2111.03606
  * Abbott et al. (2020) Abbott, R. et al. 2020 [arXiv:2010.14533]
  * Abbott et al. (2021a) Abbott, R. et al. 2021a, SoftwareX, 13, 100658
  * Abbott et al. (2021b) Abbott, R. et al. 2021b, Astrophys. J. Lett., 913, L7
  * Abbott et al. (2023b) Abbott, R. et al. 2023b [arXiv:2302.03676]
  * Acernese et al. (2015) Acernese, F. et al. 2015, Class. Quant. Grav., 32, 024001
  * Akutsu et al. (2021) Akutsu, T. et al. 2021, PTEP, 2021, 05A101
  * Ashton et al. (2019a) Ashton, G., Hübner, M., Lasky, P., & Talbot, C. 2019a, Bilby: A User-Friendly Bayesian Inference Library, Zenodo
  * Ashton et al. (2019b) Ashton, G., Hübner, M., Lasky, P. D., et al. 2019b, ApJS, 241, 27
  * Ashton & Talbot (2021) Ashton, G. & Talbot, C. 2021, MNRAS, 507, 2037
  * Astropy Collaboration et al. (2022) Astropy Collaboration, Price-Whelan, A. M., Lim, P. L., et al. 2022, apj, 935, 167
  * Belgacem et al. (2019) Belgacem, E. et al. 2019, JCAP, 07, 024
  * Bovy et al. (2011) Bovy, J., Hogg, D. W., & Roweis, S. T. 2011, AnApS, 5, 1657
  * Callister (2021) Callister, T. 2021 [arXiv:2104.09508]
  * Collette (2013) Collette, A. 2013, Python and HDF5 (O’Reilly)
  * Corman et al. (2022) Corman, M., Ghosh, A., Escamilla-Rivera, C., et al. 2022, Phys. Rev. D, 105, 064061
  * Dálya et al. (2022) Dálya, G. et al. 2022, Mon. Not. Roy. Astron. Soc., 514, 1403
  * Del Pozzo (2012) Del Pozzo, W. 2012, PhRvD, 86, 043011
  * Ezquiaga (2021) Ezquiaga, J. M. 2021, PhLB, 822, 136665
  * Ezquiaga & Holz (2022) Ezquiaga, J. M. & Holz, D. E. 2022, Phys. Rev. Lett., 129, 061102
  * Farr (2019) Farr, W. M. 2019, RNAAS, 3, 66
  * Finke et al. (2021) Finke, A., Foffa, S., Iacovelli, F., Maggiore, M., & Mancarella, M. 2021, JCAP, 2021, 026
  * Foreman-Mackey et al. (2014) Foreman-Mackey, D., Hogg, D. W., & Morton, T. D. 2014, ApJ, 795, 64
  * Gair et al. (2022) Gair, J. R. et al. 2022 [arXiv:2212.08694]
  * Gonzalez & Faber (1997) Gonzalez, A. H. & Faber, S. M. 1997, The Astrophysical Journal, 485, 80
  * Górski et al. (2005) Górski, K. M., Hivon, E., Banday, A. J., et al. 2005, ApJ, 622, 759
  * Gray et al. (2020) Gray, R., Hernandez, I. M., Qi, H., et al. 2020, Phys. Rev. D, 101, 122001
  * Gray et al. (2022) Gray, R., Messenger, C., & Veitch, J. 2022, Monthly Notices of the Royal Astronomical Society, 512, 1127
  * Harris et al. (2020) Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357
  * Hinton (2016) Hinton, S. R. 2016, The Journal of Open Source Software, 1, 00045
  * Hunter (2007) Hunter, J. D. 2007, Computing in Science & Engineering, 9, 90
  * Karathanasis et al. (2022) Karathanasis, C., Mukherjee, S., & Mastrogiovanni, S. 2022, arXiv, arXiv:2204.13495
  * Karathanasis et al. (2022) Karathanasis, C., Revenu, B., Mukherjee, S., & Stachurski, F. 2022 [arXiv:2210.05724]
  * Lagos et al. (2019) Lagos, M., Fishbach, M., Landry, P., & Holz, D. E. 2019, Phys. Rev. D, 99, 083504
  * Leyde et al. (2022) Leyde, K., Mastrogiovanni, S., Steer, D. A., Chassande-Mottin, E., & Karathanasis, C. 2022, arXiv, arXiv:2202.00025
  * Liu et al. (2023) Liu, B., Li, Z., Zhao, S., Zhou, H., & Gao, H. 2023, ApJ, 943, 29
  * Loredo & Lamb (2002) Loredo, T. J. & Lamb, D. Q. 2002, Phys. Rev. D, 65, 063002
  * Loredo & Wasserman (1998) Loredo, T. J. & Wasserman, I. M. 1998, ApJ, 502, 75
  * Madau & Dickinson (2014) Madau, P. & Dickinson, M. 2014, Ann. Rev. Astron. Astrophys., 52, 415
  * Malmquist (1922) Malmquist, K. G. 1922, MeLuF, 100, 1
  * Mancarella & Genoud-Prachex (2022) Mancarella, M. & Genoud-Prachex, E. 2022, CosmoStatGW/MGCosmoPop: v1.0.0
  * Mancarella et al. (2021) Mancarella, M., Genoud-Prachex, E., & Maggiore, M. 2021 [arXiv:2112.05728]
  * Mandel et al. (2019) Mandel, I., Farr, W. M., & Gair, J. R. 2019, MNRAS, 486, 1086
  * Mastrogiovanni et al. (2021) Mastrogiovanni, S., Leyde, K., Karathanasis, C., et al. 2021, PhRvD, 104, 062009
  * Mastrogiovanni et al. (2023) Mastrogiovanni, S. et al. 2023, Tutorials associated with ”Icarogw: A python package for inference of population properties of noisy, heterogeneous and incomplete observations”
  * Mastrogiovanni (2023) Mastrogiovanni, S. e. a. 2023, In preparation
  * Miller et al. (2020) Miller, S., Callister, T. A., & Farr, W. M. 2020, The Astrophysical Journal, 895, 128
  * Okuta et al. (2017) 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)
  * Schmidt et al. (2015) Schmidt, P., Ohme, F., & Hannam, M. 2015, Phys. Rev. D, 91, 024043
  * Schutz (1986) Schutz, B. F. 1986, Nature, 323, 310
  * Talbot et al. (2019) Talbot, C., Smith, R., Thrane, E., & Poole, G. B. 2019, Phys. Rev. D, 100, 043030
  * Talbot et al. (2019) Talbot, C., Smith, R., Thrane, E., & Poole, G. B. 2019, Phys. Rev. D, 100, 043030
  * Turski et al. (2023) Turski, C., Bilicki, M., Dálya, G., Gray, R., & Ghosh, A. 2023, arXiv, arXiv:2302.12037
  * Van Rossum (2020) Van Rossum, G. 2020, The Python Library Reference, release 3.8.2 (Python Software Foundation)
  * Virtanen et al. (2020) Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nature Methods, 17, 261
  * Vitale et al. (2022) Vitale, S., Gerosa, D., Farr, W. M., & Taylor, S. R. 2022, in Handbook of Gravitational Wave Astronomy, 45
  * Winn & Fabrycky (2015) Winn, J. N. & Fabrycky, D. C. 2015, ARA&A, 53, 409
  * Wysocki et al. (2019) Wysocki, D., Lange, J., & O’Shaughnessy, R. 2019, Phys. Rev. D, 100, 043012
  * Zheng et al. (2023) Zheng, L.-M., Li, Z., Chen, Z.-C., Zhou, H., & Zhu, Z.-H. 2023, PhLB, 838, 137720
  * Zonca et al. (2019) Zonca, A., Singer, L., Lenz, D., et al. 2019, Journal of Open Source Software, 4, 1298