# Beyond Gaps and Bumps: Spectral Siren Cosmology with Non-Parametric
Population Models

Ignacio Magaña Hernandez McWilliams Center for Cosmology, Department of
Physics, Carnegie Mellon University, Pittsburgh, PA 15213, USA
<EMAIL_ADDRESS>Anarya Ray University of Wisconsin-Milwaukee,
Milwaukee, WI 53201, USA<EMAIL_ADDRESS>
###### Abstract

Gravitational wave standard sirens typically require electromagnetic (EM) data
to obtain redshift information to constrain cosmology. Difficult to find EM
counterparts for bright sirens and galaxy survey systematics for dark sirens
make cosmological constraints with spectral sirens, a gravitational wave data-
only approach, extremely appealing. In this work, we use the GWTC-3 BBH
detections as spectral sirens to constrain the BBH population and the
underlying cosmological expansion with a flexible model for the black hole
mass spectrum. We use a binned Gaussian process to model the BBH mass
distribution in the source frame without any astrophysical assumptions on the
shape and or inclusion (or lack of) features that drive the cosmological
constraints as the redshifted detector frame masses become consistent with the
underlying astrophysical mass distribution features. For GWTC-3 we find a
measurement on the Hubble constant of $H_{0}=73.0^{+13.3}_{-7.7}$ $\mbox{km
s}^{-1}\,\mbox{Mpc}^{-1}$ at $68\%$ C.L. when combined with that obtained from
the bright standard siren analysis with GW170817 and its associated host
galaxy NGC 4993. We find an improved estimate for the Hubble constant of
around a factor of 1.4 times better than the GW170817 measurement alone. We
validate our nonparametric spectral siren approach with simulations and
benchmark its scalability and constraining performance when compared with
parametric methods.

gravitational-waves

##  1 Introduction

With the recent release of the third Gravitational-wave Transient Catalog
(Abbott et al., 2021b, GWTC-3) from the LIGO Scientific, Virgo, and KAGRA
Collaborations (LVK, Aasi et al., 2015; Acernese et al., 2015; Akutsu et al.,
2021; Abbott et al., 2018), the number of confident gravitational wave (GW)
detections from compact binary mergers is around 70, most being from binary
black hole (BBH) mergers. The increasing size of GW catalogs has enabled the
study of the BBH population (Abbott et al., 2021c, d), the cosmic expansion
history (Abbott et al., 2021a, e) as well as allowing for extensive tests of
general relativity (GR) within the strong field regime (Abbott et al., 2023a,
2021f, b). With half of the LVK’s fourth observing run (O4) having concluded,
the number of additional BBH merger detections has already doubled, and by its
conclusion, we expect to have $\mathcal{O}(300)$ additional BBH mergers and
$\mathcal{O}(10)$ mergers that contain at least one neutron star with
$\mathcal{O}(1)$ being a multimessenger event.

Gravitational waves have become a promising avenue to study the cosmic
expansion history of the universe. Measurements of the Hubble constant
($H_{0}$) with local probes such as type 1A supernovae standard candles (Riess
et al., 2022) are in tension with early universe constraints from the cosmic
microwave background (CMB) (Planck Collaboration, 2018). These two state-of-
the-art measurements provide an independent way of determining the Hubble
constant at the percent level, however, there is currently a larger than
5-sigma tension between these constraints, implying either new physics or
unaccounted sources of systematics (Riess et al., 2022; Di Valentino et al.,
2021).

To potentially resolve this tension, one ideally needs a third independent
cosmological probe. Gravitational wave sources are so-called standard sirens
(Schutz, 1986; Holz & Hughes, 2005), that is, they provide an absolute
measurement for the luminosity distance to the source without the need for the
cosmic distance ladder as a calibrating step. To do cosmology with GW sources,
we require an independent measurement of the source redshift.

For _bright sirens_ , such as binary neutron star mergers we can attempt to
detect the associated electromagnetic counterpart (EM) e.g., a kilonovae, and
measure its redshift directly (Abbott et al., 2017a, b). For _dark sirens_ ,
such as binary black hole mergers or bright sirens without a detectable (or
missed) EM counterpart, we can determine the redshift statistically using
galaxy surveys as prior information on the potential host galaxies for these
sources (Schutz, 1986; Del Pozzo, 2012; Nair et al., 2018; Chen et al., 2018;
Fishbach et al., 2019; Gray et al., 2019; Soares-Santos et al., 2019; Abbott
et al., 2021a; Palmese et al., 2020; Mukherjee et al., 2021a; Diaz &
Mukherjee, 2022; Ghosh et al., 2023; Mukherjee et al., 2021b; Gray et al.,
2023, 2022). The growing catalog of BBH mergers is thus critical for
cosmological studies since multimessenger events have proven difficult to
find.

However, even without accessible electromagnetic information, either as EM
counterparts for direct $H_{0}$ measurements or complete enough galaxy
surveys, one can still make a statistical measurement of redshift using the
features of the population distribution of compact binary mergers. Since we
measure redshifted detector-frame masses, $m^{\rm det}=m(1+z)$, one can model
the expected source frame mass distribution to estimate the redshift $z$. So
called _spectral sirens_ (Ezquiaga, 2021; Ezquiaga & Holz, 2022) therefore
allow for a direct measurement of the cosmic expansion history with
gravitational wave data alone.

The mass spectrum of LIGO–Virgo–KAGRA events introduces at least five
independent mass features (Abbott et al., 2023a): the upper and lower edges of
the pair-instability supernova (PISN) gap (Woosley et al., 2002; Heger &
Woosley, 2002; Heger et al., 2003; Woosley, 2017, 2019), the upper and lower
edges of the neutron star–black hole gap (Ye & Fishbach, 2022), and the
minimum neutron star mass (Suwa et al., 2018). However, the location of these
features is still uncertain, and some BBH formation channels might form events
with masses in the gaps, e.g., hierarchical mergers (Gerosa & Fishbach, 2021)
leading to a more complex mass spectrum. By using the full mass distribution,
degeneracies between mass evolution and cosmological evolution can be broken.
This self-calibrating _spectral siren_ method has the potential to provide
precision constraints of both cosmology and the potential evolution of the
mass distribution. For recent works on cosmological inference that relies on
the astrophysical mass distribution of CBCs see Mastrogiovanni et al. (2021,
2023, 2024); Gray et al. (2023)

However, these studies (Farr, 2019; Ezquiaga & Holz, 2022; Chen et al., 2024)
have assumed _a priori_ some knowledge about the underlying astrophysical
processes to forward model the CBC population with simple parametric
population models. For example, modeling expected features in the mass
spectrum which are then used to measure the cosmological expansion. However,
Pierra et al. (2023) has shown that incorrect assumptions regarding the shape
of the mass-spectrum and its redshift evolution can lead to significant
($3\sigma$) biases in the inferred cosmological parameters. For more details
regarding the redshift evolution of the mass spectrum and its impact on the
associated cosmological inference, see Mukherjee (2022); Karathanasis et al.
(2023).

Therefore, given the significant uncertainty regarding known BBH formation
models and the new features found in the BBH mass spectrum (Abbott et al.,
2023a), there is a need for a flexible and non-parametric approach to spectral
siren cosmology. Examples of non-parametric population models include:
autoregressive processes (Callister & Farr, 2023a), splines (Edelman et al.,
2022a, b), Gaussian mixture models (Tiwari, 2021; Tiwari & Fairhurst, 2021;
Tiwari, 2022), adaptive kernel density estimation Sadiq et al. (2022), maximum
population likelihood (Payne & Thrane, 2023), Dirichlet processes (Rinaldi &
Del Pozzo, 2021) and binned Gaussian processes (GPs, Ray et al. (2023); Mohite
(2022); Abbott et al. (2023a)).

In this work, we employ, for the first time, a flexible population model on
the BBH population distribution and perform hierarchical Bayesian inference on
the GWTC-3 BBH detections to measure the cosmic expansion history using the
spectral siren methodology. We model the BBH population using the binned
gaussian process (BGP) approach (Ray et al., 2023, 2024; Mohite, 2022; Abbott
et al., 2023a) to place constraints on the Hubble constant without assuming a
particular shape on the mass distribution of BBH mergers. For a discussion
regarding how our method compares with a parallel investigation on non-
parametric spectral sirens (Farah et al., 2024), see Sec. 5.

This paper is organized as follows. In Section 2, we summarize the
hierarchical Bayesian framework describing the BGP spectral siren model used
in our analysis. In Section 3, we present the main results of this paper by
analyzing the GWTC-3 BBH observations with our framework. In Section 4, we
test our methodology with a simulated population that mimicks the GWTC-3
population and current GW catalog sizes. Finally, in Section 5 we provide a
summary of our work and discuss future directions.

##  2 Methods

In this section, we describe the hierarchical inference framework used to
simultaneously infer the cosmological parameters with the shape of the CBC
mass spectrum.

###  2.1 Flexible population model

Following Ray et al. (2023); Mohite (2022), we construct our flexible
population model in the source frame as a piece-wise binned function in masses
and redshift, such that,

$\displaystyle\frac{dN}{dm_{1}dm_{2}dz}(m_{1},m_{2},z|n^{\gamma},n^{\alpha})=\frac{n^{\gamma}_{m}n^{\alpha}_{z}}{m_{1}m_{2}}\frac{dV}{dz}T_{\rm{obs}}(1+z)^{-1},$
(1)

where $n^{\gamma}_{m}n^{\alpha}_{z}$ is the merger rate density of CBCs per
log component mass per co-moving volume per unit time in the $\gamma$th mass
bin and the $\alpha$th redshift bin. While the model in Equation 1 is more
flexible than the one used in Mohite (2022) which assumes no evolution of the
merger rate, it is less so than the correlated model of Ray et al. (2023)
where the joint distribution of masses and redshifts were modelled using a
single piece-wise binned function.

Keeping in mind that the primary objective of measuring the Hubble parameter
independent of the uncertainties regarding CBC formation channels, and the
large measurement uncertainties in the mass-redshift correlations reported by
Ray et al. (2023), we expect the un-correlated mass-redshift model given by
Equation 1 to cause no significant biases in our results given the size of
current datasets. However, as catalogs continue to grow, we hope to relax this
restriction on the existence of mass-redshift correlations in future studies
for which we will use the generalized correlated model of Ray et al. (2023).

To fit a population model such as Equation 1 to GW observations of CBC
systems, it is necessary to convert detector frame observables such as
redshifted masses and luminosity distances into their corresponding source-
frame counterparts using a fixed cosmological model. In the next subsection,
following the spectral siren method of Farr (2019); Ezquiaga & Holz (2022), we
describe how to vary and infer the cosmological parameters simultaneously with
the population distribution in the context of our flexible model.

###  2.2 Simultaneous cosmological inference

Given a cosmological model characterized by parameters $\Omega$, it is
possible to express the population model described in the previous section as
a function of detector frame CBC observables instead,

$\displaystyle\frac{dN}{dm_{1}dm_{2}dz}(m^{d}_{1},m^{d}_{2},D_{L}|n^{\gamma}_{m},n^{\alpha}_{z},\Omega)$
$\displaystyle=\frac{n^{\gamma}_{m}n^{\alpha}_{z}}{m^{d}_{1}m^{d}_{2}}\frac{dV_{c}}{dz}\left(D_{L}|\Omega\right)T_{\rm{obs}}\left(1+z(D_{L}|\Omega)\right),$
(2)

where $m_{1}^{d}$, $m_{2}^{d}$ and $D_{L}$ are the observed detector frame
masses and luminosity distances of CBCs and $z(D_{L}|\Omega)$ is the redshift
of the source as a function of luminosity distance given a set of cosmological
parameters. As usual, $dV_{c}/dz$ is the differential uniform-in-comoving
volume element, $T_{\mathrm{obs}}$ is the total observation time, and the
extra factor of $1/(1+z)$ converts source-frame time to detector-frame time.

By modeling the occurrence of CBCs as realizations of an inhomogeneous Poisson
process, it is possible to infer hyper-parameters characterizing the spectral
siren model of Equation 2.2 from the posterior samples of masses and
luminosity distances of each observed CBC using Bayesian hierarchical
inference (Mandel et al., 2019; Loredo, 2004; Vitale et al., 2020). Within
this framework, the likelihood of population and cosmological hyper-parameters
given GW data ($d$) from all events in an observed catalog can be constructed
as,

$\displaystyle p(\vec{d}|\vec{n}_{m},\vec{n}_{z},\Omega)$
$\displaystyle=e^{-N_{\rm{det}}(\vec{n}_{m},\vec{n}_{z},\Omega)}\prod_{i}^{N_{\rm{obs}}}\left\langle\frac{\frac{dN}{dm_{1}dm_{2}dz}(\vec{n}_{m},\vec{n_{z}},\Omega)}{p_{\text{PE}}(m_{1},m_{2},z)}\right\rangle_{\text{samples},i},$
(3)

where $\left\langle\cdot\right\rangle_{{\rm{samples}},i}$ represents a Monte
Carlo sum over posterior samples from the ith event reweighted by the prior
used in its parameter estimation and $N_{\rm{det}}$ is the expected number of
detections as a function of hyper-parameters.

Since the criteria for detection takes the form of a threshold imposed on some
statistic such as false alarm rate, the estimation of $N_{\rm{det}}$ has to
account for Malmquist biases resulting from the imposition of the detection
criteria (Mandel et al., 2019; Loredo, 2004). This is implemented by
simulating a fiducial population of CBCs and injecting them into detector
noise realizations Essick & Farr (2022); Farr (2019). The set of simulated
sources that are detected above thresholds can be re-weighted to the
population model of Equation 2.2 to estimate
$N_{\rm{det}}(\vec{n}_{m},\vec{n}_{z},\Omega)$ in the following way:

$\displaystyle N_{\rm{det}}(\vec{n}_{m},\vec{n}_{z},\Omega)$
$\displaystyle=\frac{K_{\rm{det}}}{K_{\rm{draw}}}\left\langle\frac{\frac{dN}{dm_{1}dm_{2}dz}(\vec{n}_{m},\vec{n}_{z},\Omega)}{p_{\rm{draw}}(m_{1},m_{2},z)}\right\rangle_{\text{samples,det}}$
(4)

where similarly, $\left\langle\cdot\right\rangle_{\text{samples,det}}$
represents a Monte Carlo sum over the parameters of detectable events re-
weighted by the population used to generate the mock simulations. For a
discussion regarding the convergence of Monte Carlo sums used to evaluate
various terms in the likelihood see Appendix A

The likelihood defined in Equation 3 can then be used to inform the shape of
the CBC mass spectrum, the redshift evolution of the merger rate and the
cosmological parameters self-consistently given a catalog of GW events.

###  2.3 Gaussian process prior and HMC sampling

To infer the shape of the mass distribution using the likelihood of Eq. (3),
we draw the rate densities in each mass bin from a Gaussian process prior (Ray
et al., 2023; Mohite, 2022) so that the posterior on the rate densities is
given by,

$\displaystyle
p(\vec{n}_{m},\vec{n}_{z},\vec{\mu},\vec{\sigma},\vec{l}|\vec{d},\Omega)\propto
p(\vec{d}|\vec{n}_{m},\vec{n}_{z},\Omega)$ (5) $\displaystyle\times
p(\vec{\mu},\vec{\sigma},\vec{l})p(\vec{n}_{m}|\vec{\mu}_{m},\sigma_{m},l_{m})$
$\displaystyle\times p(\vec{n}_{z}|\vec{\mu}_{z},\sigma_{z},l_{z})$

where $p(\vec{n}|\vec{\mu},\sigma,\lambda)$ is the GP prior and $\vec{\mu}$ is
the mean function of the GP while $l,\sigma$ are parameters that control the
correlation length and amplitudes of the GP’s covariance matrix.

The prior on the GP hyperparameters is defined by $p(\vec{\mu},\sigma,l)$.
Following previous implementations of the binned Gaussian process population
analysis, we chose the priors on the GP hyper-parameters
($\vec{\mu},l,\sigma$) to be Normal, Log Normal, and Half Normal respectively
while modeling the covariance matrix as an exponential quadratic function.

Along with the rate-densities, we simultaneously draw the cosmological
parameters uniform priors and use Equation 3 to sample the posterior
distribution,

$\displaystyle p(\vec{n}_{m},\vec{n}_{z},\Omega|\vec{d})\propto
p(\vec{n}_{m},\vec{n}_{z},\vec{\mu},\vec{\sigma},\vec{l}|\vec{d},\Omega)p(\vec{\Omega})$
(6)

The posterior samples of $\vec{\Omega}$ represent measurements of the
cosmological expansion that have been marginalized over the uncertainties
regarding the shape of the CBC mass spectrum without making any astrophysical
assumptions about any features present in the mass spectrum. Simultaneously,
the samples of $\vec{n}$ can then be used to reconstruct the shape of the CBC
mass distribution in a data-driven way, while being free of any biases that
might result from the uncertainties in the measurements of cosmological
parameters.

##  3 Results for GWTC-3

In this section, we present our BGP spectral siren constraints using GWTC-3
BBH mergers. Following the LVK GWTC-3 population analysis choices Abbott et
al. (2021d), we make use only of the BBH events that pass a 1 per year IFAR
threshold. That is, we analyze a total of 69 confidently detected BBH mergers
and exclude GW190814 as a population outlier. For the events incorporated in
our analysis, the parameter estimation samples of detector frame masses and
luminosity distances are used to compute the posterior weights of Equation 3.
We also use the publicly available LVK’s GWTC-3 sensitivity estimate injection
campaigns to compute Equation 4 (The LIGO Scientific Collaboration et al.,
2021, 2023).

We infer the mass rate densities as well as the redshift evolution of the BBH
merger rate jointly with the cosmological parameters $H_{0},\Omega_{m}$. In
Figure 1, we show the inferred marginal posterior distributions on both the
primary mass $m_{1}$ and secondary mass $m_{2}$ along with the redshift
evolution of the combined merger rate. We also show the corresponding LVK
marginal posterior distributions on the masses for the Powerlaw+Peak model for
reference Talbot & Thrane (2018); Abbott et al. (2021d). Our results show
broad consistency with the LVK Powerlaw+Peak posterior. Our results are also
consistent with the BGP results presented in Ray et al. (2023) which assume a
fixed cosmological model.

In Figure 2, we show the corner plot for the joint posterior distribution on
the Hubble constant $H_{0}$, the matter density parameter $\Omega_{m}$, as
well as the GP kernel length scales $l_{m},l_{z}$. We find a BGP spectral
siren BBH measurement on the Hubble constant of $H_{0}=82.2^{+51.7}_{-33.7}$
$\mbox{km s}^{-1}\,\mbox{Mpc}^{-1}$ while for the matter density parameter, we
find $\Omega_{m}=0.3^{+0.1}_{-0.1}$ both at $68\%$ C.L. We measure length
scales $l_{m}=2.1^{+1.0}_{-1.1}$, $l_{z}=2.1^{+1.0}_{-1.1}$ at $68\%$ C.L.

In Figure 3, we show joint constraints for our BGP spectral siren measurement
on the Hubble constant with the 69 BBHs from GWTC-3 only when combined with
the bright siren $H_{0}$ measurement from GW170817 and its EM counterpart NGC
4993 (Abbott et al., 2017a, b). We find a joint $H_{0}$ measurement of
$H_{0}=73.0^{+13.3}_{-7.7}$ $\mbox{km s}^{-1}\,\mbox{Mpc}^{-1}$ at $68\%$ C.L.
which is an improvement with a factor of 1.5 times better relative to the
GW180817 with NCG 4993 measurement alone ($H_{0}=71^{+23}_{-8}$ $\mbox{km
s}^{-1}\,\mbox{Mpc}^{-1}$ at $68\%$ C.L.). In all of our reported $H_{0}$
measurements, we used a uniform prior in the range $[20,140]\ $ $\mbox{km
s}^{-1}\,\mbox{Mpc}^{-1}$.

Figure 1:  Inferred constraints for both the primary and secondary mass
distributions using the 69 BBH mergers from GWTC-3 and our BGP spectral siren
model. We show our results as the binned piecewise function plotted in purple
and for comparison, we also show the corresponding results obtained by the LVK
using the Powerlaw + Peak model in Abbott et al. (2021d) as the gray band
depicting the 95% posterior credible region. Figure 2:  Corner plot showing
the joint marginalized posterior distribution on the Hubble constant $H_{0}$,
the matter density parameter $\Omega_{m}$, and the inferred length scales
$l_{m},l_{z}$ from our BGP spectral siren model and the 69 GWTC-3 BBH events.
Figure 3:  Marginalized posterior distributions on the Hubble constant
$H_{0}$. In green, we show the posterior from the 69 GWTC-3 BBH events we have
considered in the spectral siren analysis using our BGP model. In gray, we
show the GW170817 $H_{0}$ posterior obtained with its host galaxy NGC4993. In
blue, we show the resulting combined GW170817 and GWTC-3 BBH spectral siren
posterior to obtain a joint measurement on $H_{0}$. For reference, in the
orange band we show measurements on $H_{0}$ from the CMB (Aghanim et al.,
2020) and in the pink shaded band from standard candle type 1A SN measurements
(Riess et al., 2019).

##  4 Validation of Results

We validate our inference on simulated GW catalogs made from events drawn from
a known fiducial population and cosmological model. We draw simulated BBH
observations from the following models: for the primary mass, we draw from the
Powerlaw + Peak model (Talbot & Thrane, 2018), mass ratios are drawn from a
powerlaw in $q$ with slope $\beta$ (Fishbach & Holz, 2020) and we allow for
the merger rate to evolve with redshift (Fishbach et al., 2018a). The fiducial
hyperparameter values for the simulated population correspond to
$m_{\rm{min}}=5M_{\odot}$, $m_{\rm{max}}=65M_{\odot}$, $\alpha=3.14$,
$\beta=1.4$, $m_{{\rm peak}}=35M_{\odot}$, $\sigma_{{\rm peak}}=5M_{\odot}$,
$f_{{\rm peak}}=0.01$, and $\kappa=3$. For explicit expressions, we refer the
reader to Appendix 1 of X and to Figure 4 for the marginal posterior
distributions over $m_{1}$ and $m_{2}$ that we use in our simulations. We use
a spatially flat $\Lambda$CDM cosmological model with assumed Planck 2015
cosmological parameters, that is, with $H_{0}=67.8\ \rm{km\ s^{-1}\
{Mpc}^{-1}}$ and $\Omega_{m}=0.308$ in our simulations.

The drawn masses and redshifts from the fiducial population are then assigned
corresponding observed values calibrated to the expected measurement
uncertainties for BBH mergers at advanced LIGO design sensitivity (Chen et
al., 2018; Abbott et al., 2018; Fishbach et al., 2018b) and selected so that
only a subset is observable, i.e., simulating the selection effects seen in GW
observations. The details of this procedure are described in Farah et al.
(2023); Fishbach et al. (2020, 2018b). A corresponding set of simulated
detectable injections is generated from a broad distribution to take into
account the estimation of selection effects, a requirement for estimating
$N_{\rm{det}}$ through Equation 4.

Our catalog is comprised of 144 observed BBH events drawn from the fiducial
population model and underlying cosmological model. As with our analysis of
the GWTC-3 BBH events, we employ the BGP spectral siren model to fit our
simulated population with the model defined in Equation 2.2. In Figure 4 we
show the inferred marginal posterior distributions on both the primary mass
$m_{1}$ and secondary mass $m_{2}$ using our simulated data. We can see that
the BGP model is sufficiently accurate at recovering the shape of the fiducial
mass population with the 121 events that we consider. Although there is some
discrepancy in the inferred shape of the primary mass population at high
masses, we are still broadly consistent with the simulated population. We note
that this discrepancy is due to a small number of events out of the total in
our catalog being at the higher end of the mass spectrum as well as due to the
increased uncertainty in their corresponding single-event parameter estimation
posterior distributions.

We infer the BBH mass distribution as well as the redshift evolution of the
BBH merger rate jointly with the Hubble constant $H_{0}$. In Figure 4, we show
the inferred marginal posterior distributions on both the primary mass
$m_{1}$, secondary mass $m_{2}$, and redshift. We also show the fiducial
populations for comparison. In figure 5, we show the posterior distributions
of the cosmological parameters. The BGP spectral siren method can accurately
reconstruct the shape of the injected population as well as the fiducial
values of the cosmological parameters.

Figure 4:  Inferred constraints for both the primary and secondary mass
distributions from the simulated population described in Section 4 with the
BGP spectral siren model. We show our results as the binned piecewise function
plotted in purple. The simulated population is shown as the solid green lines.
Figure 5:  Corner plot showing the joint posterior distributions on the Hubble
constant $H_{0}$, the matter density parameter $\Omega_{m}$, and the inferred
GP length scales from our BGP Spectral Siren model on the simulated population
described in Section 4. We plot the simulated values for both $H_{0}=67.8\
\rm{km\ s^{-1}\ {Mpc}^{-1}}$, $\Omega_{m}=0.308$ for reference in gray.

##  5 Discussion and Conclusion

In this work, we have explored how the spectral siren framework can be used
jointly with nonparametric approaches in the modeling of the BBH population.
In particular, we used a binned Gaussian process to model the black hole mass
distribution without assuming any particular shape inclusive (or exclusive) of
features expected from the astrophysics of BBH formation channels. Our proof
of principle demonstration made use of the GWTC-3 confident BBH detections to
place the first nonparametric constraints on the Hubble constant. We have
validated our approach with GWTC-3 like populations and have made use of
current catalog sizes to demonstrate the constraining power of our
nonparametric method relative to traditional parametric approaches to spectral
siren cosmology.

In a parallel investigation, Farah et al. (2024) have also attempted to use
Gaussian-process-based mass models on a simulated dataset for a proof of
concept demonstration of astrophysics agnostic spectral siren cosmology. Even
though their non-parametric population model is not restricted by bin
resolutions, the mock data analysis presented in Farah et al. (2024) models
_only_ the primary mass distribution of CBCs with gaussian processes. They
further assume that all BBHs are equal mass systems and use a parametric model
for the redshift evolution of the merger rate. By simulating a large
population of equal mass BBHs with a known functional form of the redshift
evolution of the merger rate, they demonstrate that their method can implement
accurate spectral siren cosmology in a simplistic Universe while being free of
the necessity to bin up the space of BBH parameters. We note that while an
unbinned population model in both component masses and redshifts could in
principle be more flexible than ours, Farah et al. (2024) have not
demonstrated whether their method can be scalably generalized and applied to
more realistic datasets wherein BBH masses are expected to follow more
complicated pairing functions (Abbott et al., 2023a) than a delta function in
mass ratio centered around one. Furthermore, given the strong correlations
between the measurement uncertainties of the Hubble parameter and the redshift
evolution parameters (Abbott et al., 2021e), previously unmodelled features in
the redshift evolution function can serve to equally bias the cosmological
parameters. On the other hand, our model can reconstruct features in the joint
distribution of component masses and the redshift evolution function
completely driven by the data, up to the resolution limit imposed by our
choice of binning, variations of which are found to yield consistent results.

As the size of GW catalogs grow, we expect to expand our nonparametric data-
driven approach to take advantage of potential correlations between mass and
redshift. In this work, we have assumed an uncorrelated mass-redshift
distribution for simplicity. This choice can be relaxed as was shown in Ray et
al. (2023) where a single three dimensional BGP infers the shape of the joint
mass-redshift distribution. However, as it was shown in Ray et al. (2023),
most of the structure is washed out due to the increased dimensionality of the
problem and thus would require larger GW catalogs to measure such
correlations. This is important, if say for examples, the BBH mass
distribution evolves at roughly the same rate as the cosmic expansion history
as argued by Ezquiaga & Holz (2022). We leave these simulated and systematic
studies as future work.

Another caveat of our analysis is the fact that we only use BBH mergers to
infer the cosmic expansion rate. Our methodology is CBC agnostic, that is, we
fit for the merger rate density in mass bins. Including the whole population
of CBC mergers (such as the inclusion of BNS and NSBH mergers) should make use
of the full mass distribution of CBCs, a more powerful and event category
agnostic spectral siren probe as explored by Ezquiaga & Holz (2022).

Our method can also be expanded to be model agnostic concerning the underlying
functional of the background cosmological model. In our work, we have assumed
a spatially flat $\Lambda$CDM cosmological model which is sufficient given the
number of observations that we have analyzed and their associated single event
parameter uncertainties. Fitting the luminosity distance and redshift
relation, $D_{L}(z)$, with a GP may provide a way to test the underlying
cosmological model with future networks of GW detectors.

Finally, nonparametric methods are data-driven by nature and therefore require
large amounts of observations for meaningful inference. As such, the
scalability of such methods is critical. However, for the next generation of
GW detectors (XG/3G), the development of scalable and potentially approximate
inference will likely be necessary to make use of the
$\mathcal{O}(10^{4}-10^{6})$ expected GW detections catalog sizes (Chen et
al., 2024).

## Acknowledgements

The authors would like to thank Jolien Creighton for his useful comments and
feedback throughout this work. IMH is supported by a McWilliams postdoctoral
fellowship at Carnegie Mellon University. AR is supported by the National
Science Foundation award PHY-2207728. The authors are grateful for
computational resources provided by the LIGO Laboratory and supported by
National Science Foundation Grants PHY-0757058 and PHY-0823459, PHY-0823459,
PHY-1626190, and PHY-1700765. This material is based upon work supported by
NSF’s LIGO Laboratory which is a major facility fully funded by the National
Science Foundation.

##  Appendix A Convergence of Monte-Carlo integrals

As mentioned in section 2.2, the Monte-Carlo sums used to compute various
terms in Eq. (3) are susceptible to uncertainties which might raise concerns
regarding their convergence. Farr (2019) shows that marginalizing over Monte
Carlo uncertainties renders sampling of the marginalized posterior unable to
find the maxima in the space of hyper-parameters unless sparsely populated
regions of that space are penalized by means of a likelihood. Specifically,
the effective number of single-event posterior samples and detectable
injections are required to be high enough compared to the total number of
observations for unbiased inference. Essick & Farr (2022); Farr (2019) outline
how to impose this requirement in the context of an hierarchical framework
such as the one employed in this study.

Following the implementations of the methods of Essick & Farr (2022); Farr
(2019) in the context of the binned model, as elaborated in Ray et al. (2023),
we compute the means and standard deviations of the Monte Carlo sums used
evaluate the likelihood of Eq. (3). In particular, we compute the mean and
deviations of the event specific posterior-weights and the detectable time-
volumes in each bin as a function of the cosmological parameters, which take
the following forms:

$\displaystyle\sum_{\gamma,\alpha}n^{\gamma}_{m}n^{\alpha}_{z}\mu_{w,i}^{\gamma\alpha}(\Omega)$
$\displaystyle=$
$\displaystyle\left\langle\frac{\frac{dN}{dm_{1}dm_{2}dz}(\vec{n}_{m},\vec{n}_{z})}{p_{\text{PE}}(m_{1},m_{2},z)}\right\rangle_{\text{samples},i}$
(A1)
$\displaystyle\sum_{\gamma,\alpha}\left(n^{\gamma}_{m}n^{\alpha}_{z}\sigma_{w,i}^{\gamma\alpha}(\Omega)\right)^{2}$
$\displaystyle=$
$\displaystyle\left\langle\left(\frac{\frac{dN}{dm_{1}dm_{2}dz}(\vec{n}_{m},\vec{n_{z}})}{p_{\text{PE}}(m_{1},m_{2},z)}\right)^{2}\right\rangle_{\text{samples},i}-\frac{1}{N_{\text{samples},i}}\left(\sum_{\gamma,\alpha}n^{\gamma}_{m}n^{\alpha}_{z}\mu_{w,i}^{\gamma\alpha}(\Omega)\right)^{2}$
(A2)
$\displaystyle\sum_{\gamma,\alpha}n^{\gamma}_{m}n^{\alpha}_{z}\mu_{\rm{VT}}^{\gamma\alpha}(\Omega)$
$\displaystyle=$
$\displaystyle\frac{K_{\text{det}}}{K_{\text{draw}}}\left\langle\frac{\frac{dN}{dm_{1}dm_{2}dz}(\vec{n}_{m},\vec{n_{z}})}{p_{\text{draw}}(m_{1},m_{2},z)}\right\rangle_{\text{samples,det}}$
(A3)
$\displaystyle\sum_{\gamma,\alpha}\left(n^{\gamma}_{m}n^{\alpha}_{z}\sigma_{\rm{VT},i}^{\gamma}(\Omega)\right)^{2}$
$\displaystyle=$
$\displaystyle\left(\frac{K_{\text{det}}}{K_{\text{draw}}}\right)^{2}\left\langle\left(\frac{\frac{dN}{dm_{1}dm_{2}dz}(\vec{n}_{m},\vec{n_{z}})}{p_{\text{draw}}(m_{1},m_{2},z)}\right)^{2}\right\rangle_{\text{samples,
det}}-\left(\frac{K_{\text{det}}}{K_{\text{draw}}}\right)\left(\sum_{\gamma,\alpha}n^{\gamma}_{m}n^{\alpha}_{z}\mu_{\rm{VT}}^{\gamma\alpha}(\Omega)\right)^{2}$
(A4)

. The effective number of samples for the detectable injections and event-
specific posteriors can be computed from these means and, deviations to be
$N_{\text{eff}}=\left(\frac{\mu}{\sigma}\right)^{2}$ and are required to be
large enough to avoid biases in the inferred hyper-parameters arising from
Monte Carlo uncertainties. In particular, following Ray et al. (2023);
Callister & Farr (2023b), we demand the following conditions:

$\displaystyle\frac{N_{\rm{eff}}^{\rm{VT}}(\Omega)}{N_{\rm{det}}(\vec{n}_{m},\vec{n_{z}},\Omega)}=\frac{\sum_{\gamma,\alpha}\left(\frac{\mu_{\rm{VT}}^{\gamma\alpha}}{\sigma_{\rm{VT}}^{\gamma\alpha}}\right)^{2}}{\sum_{\gamma,\alpha}n^{\gamma}_{m}n^{\alpha}_{z}\mu_{\rm{VT}}^{\gamma\alpha}}$
$\displaystyle\geq$ $\displaystyle 2$ (A5)
$\displaystyle\min_{i}N_{\rm{eff}}^{w,i}(\vec{n}_{m},\vec{n}_{z},\Omega)=\min_{i}\left\\{\frac{\left(\sum_{\gamma,\alpha}n^{\gamma}_{m}n^{\alpha}_{z}\mu_{w,i}^{\gamma\alpha}\right)^{2}}{\sum_{\gamma,\alpha}(n^{\gamma}_{m}n^{\alpha}_{z}\sigma_{w,i}^{\gamma\alpha})^{2}}\right\\}$
$\displaystyle\geq$ $\displaystyle 10^{0.6}$ (A6)

be used on all of the hyperparameter samples that are accepted at each HMC
step. Given the bin resolution and sample-sizes used in our studies we expect
these conditions to be satisfied automatically. Hence, instead of imposing the
above conditions in the form of a likelihood penalization, we verify post
sampling that all of the hyper-parameter samples satisfy them.

## References

  * Aasi et al. (2015) Aasi, J., et al. 2015, Class. Quant. Grav., 32, 074001, doi: 10.1088/0264-9381/32/7/074001
  * Abbott et al. (2017a) Abbott, B. P., et al. 2017a, Phys. Rev. Lett., 119, 161101, doi: 10.1103/PhysRevLett.119.161101
  * Abbott et al. (2017b) —. 2017b, Nature, 551, 85, doi: 10.1038/nature24471
  * Abbott et al. (2018) —. 2018, Living Rev. Rel., 21, 3, doi: 10.1007/s41114-020-00026-9
  * Abbott et al. (2021a) —. 2021a, Astrophys. J., 909, 218, doi: 10.3847/1538-4357/abdcb7
  * Abbott et al. (2021b) Abbott, R., et al. 2021b.  https://arxiv.org/abs/2111.03606
  * Abbott et al. (2021c) —. 2021c, Astrophys. J. Lett., 913, L7, doi: 10.3847/2041-8213/abe949
  * Abbott et al. (2021d) —. 2021d.  https://arxiv.org/abs/2111.03634
  * Abbott et al. (2021e) —. 2021e.  https://arxiv.org/abs/2111.03604
  * Abbott et al. (2021f) —. 2021f.  https://arxiv.org/abs/2112.06861
  * Abbott et al. (2023a) —. 2023a, Phys. Rev. X, 13, 011048, doi: 10.1103/PhysRevX.13.011048
  * Abbott et al. (2023b) —. 2023b.  https://arxiv.org/abs/2304.08393
  * Acernese et al. (2015) Acernese, F., et al. 2015, Class. Quant. Grav., 32, 024001, doi: 10.1088/0264-9381/32/2/024001
  * Aghanim et al. (2020) Aghanim, N., Akrami, Y., Ashdown, M., et al. 2020, Astronomy & Astrophysics, 641, A6
  * Akutsu et al. (2021) Akutsu, T., et al. 2021, PTEP, 2021, 05A102, doi: 10.1093/ptep/ptab018
  * Callister & Farr (2023a) Callister, T. A., & Farr, W. M. 2023a, A Parameter-Free Tour of the Binary Black Hole Population
  * Callister & Farr (2023b) —. 2023b.  https://arxiv.org/abs/2302.07289
  * Chen et al. (2024) Chen, H.-Y., Ezquiaga, J. M., & Gupta, I. 2024.  https://arxiv.org/abs/2402.03120
  * Chen et al. (2018) Chen, H.-Y., Fishbach, M., & Holz, D. E. 2018, Nature, 562, 545, doi: 10.1038/s41586-018-0606-0
  * Del Pozzo (2012) Del Pozzo, W. 2012, Phys. Rev. D, 86, 043011, doi: 10.1103/PhysRevD.86.043011
  * Di Valentino et al. (2021) Di Valentino, E., Mena, O., Pan, S., et al. 2021, Class. Quant. Grav., 38, 153001, doi: 10.1088/1361-6382/ac086d
  * Diaz & Mukherjee (2022) Diaz, C. C., & Mukherjee, S. 2022, Mon. Not. Roy. Astron. Soc., 511, 2782, doi: 10.1093/mnras/stac208
  * Edelman et al. (2022a) Edelman, B., Doctor, Z., Godfrey, J., & Farr, B. 2022a, The Astrophysical Journal, 924, 101, doi: 10.3847/1538-4357/ac3667
  * Edelman et al. (2022b) Edelman, B., Farr, B., & Doctor, Z. 2022b, Cover Your Basis: Comprehensive Data-Driven Characterization of the Binary Black Hole Population
  * Essick & Farr (2022) Essick, R., & Farr, W. 2022.  https://arxiv.org/abs/2204.00461
  * Ezquiaga (2021) Ezquiaga, J. M. 2021, Phys. Lett. B, 822, 136665, doi: 10.1016/j.physletb.2021.136665
  * Ezquiaga & Holz (2022) Ezquiaga, J. M., & Holz, D. E. 2022, Phys. Rev. Lett., 129, 061102, doi: 10.1103/PhysRevLett.129.061102
  * Farah et al. (2024) Farah, A. M., Callister, T. A., Ezquiaga, J. M., Zevin, M., & Holz, D. E. 2024.  https://arxiv.org/abs/2404.02210
  * Farah et al. (2023) Farah, A. M., Edelman, B., Zevin, M., et al. 2023, Things that might go bump in the night: Assessing structure in the binary black hole mass spectrum
  * Farr (2019) Farr, W. M. 2019, Research Notes of the AAS, 3, 66, doi: 10.3847/2515-5172/ab1d5f
  * Farr (2019) Farr, W. M. 2019, Research Notes of the American Astronomical Society, 3, 66, doi: 10.3847/2515-5172/ab1d5f
  * Fishbach et al. (2020) Fishbach, M., Farr, W. M., & Holz, D. E. 2020, The Astrophysical Journal Letters, 891, L31, doi: 10.3847/2041-8213/ab77c9
  * Fishbach & Holz (2020) Fishbach, M., & Holz, D. E. 2020, Astrophys. J. Lett., 891, L27, doi: 10.3847/2041-8213/ab7247
  * Fishbach et al. (2018a) Fishbach, M., Holz, D. E., & Farr, W. M. 2018a, The Astrophysical Journal Letters, 863, L41, doi: 10.3847/2041-8213/aad800
  * Fishbach et al. (2018b) —. 2018b, Astrophys. J., 863, L41, doi: 10.3847/2041-8213/aad800
  * Fishbach et al. (2019) Fishbach, M., et al. 2019, Astrophys. J. Lett., 871, L13, doi: 10.3847/2041-8213/aaf96e
  * Gerosa & Fishbach (2021) Gerosa, D., & Fishbach, M. 2021, Nature Astron., 5, 749, doi: 10.1038/s41550-021-01398-w
  * Ghosh et al. (2023) Ghosh, T., More, S., Bera, S., & Bose, S. 2023.  https://arxiv.org/abs/2312.16305
  * Gray et al. (2019) Gray, R., Magaña Hernandez, I., Qi, H., Sur, A., et al. 2019.  https://arxiv.org/abs/1908.06050
  * Gray et al. (2022) Gray, R., Messenger, C., & Veitch, J. 2022, Mon. Not. Roy. Astron. Soc., 512, 1127, doi: 10.1093/mnras/stac366
  * Gray et al. (2023) Gray, R., et al. 2023, JCAP, 12, 023, doi: 10.1088/1475-7516/2023/12/023
  * Heger et al. (2003) Heger, A., Fryer, C. L., Woosley, S. E., Langer, N., & Hartmann, D. H. 2003, Astrophys. J., 591, 288, doi: 10.1086/375341
  * Heger & Woosley (2002) Heger, A., & Woosley, S. E. 2002, Astrophys. J., 567, 532, doi: 10.1086/338487
  * Holz & Hughes (2005) Holz, D. E., & Hughes, S. A. 2005, Astrophys. J., 629, 15, doi: 10.1086/431341
  * Karathanasis et al. (2023) Karathanasis, C., Mukherjee, S., & Mastrogiovanni, S. 2023, Mon. Not. Roy. Astron. Soc., 523, 4539, doi: 10.1093/mnras/stad1373
  * Loredo (2004) Loredo, T. J. 2004, AIP Conf. Proc., 735, 195, doi: 10.1063/1.1835214
  * Mandel et al. (2019) Mandel, I., Farr, W. M., & Gair, J. R. 2019, Mon. Not. Roy. Astron. Soc., 486, 1086, doi: 10.1093/mnras/stz896
  * Mastrogiovanni et al. (2021) Mastrogiovanni, S., Leyde, K., Karathanasis, C., et al. 2021, Phys. Rev. D, 104, 062009, doi: 10.1103/PhysRevD.104.062009
  * Mastrogiovanni et al. (2023) Mastrogiovanni, S., Laghi, D., Gray, R., et al. 2023, Phys. Rev. D, 108, 042002, doi: 10.1103/PhysRevD.108.042002
  * Mastrogiovanni et al. (2024) Mastrogiovanni, S., Pierra, G., Perriès, S., et al. 2024, Astron. Astrophys., 682, A167, doi: 10.1051/0004-6361/202347007
  * Mohite (2022) Mohite, S. 2022, PhD thesis, Wisconsin U., Milwaukee
  * Mukherjee (2022) Mukherjee, S. 2022, Mon. Not. Roy. Astron. Soc., 515, 5495, doi: 10.1093/mnras/stac2152
  * Mukherjee et al. (2021a) Mukherjee, S., Wandelt, B. D., Nissanke, S. M., & Silvestri, A. 2021a, Phys. Rev. D, 103, 043520, doi: 10.1103/PhysRevD.103.043520
  * Mukherjee et al. (2021b) Mukherjee, S., Wandelt, B. D., & Silk, J. 2021b, Mon. Not. Roy. Astron. Soc., 502, 1136, doi: 10.1093/mnras/stab001
  * Nair et al. (2018) Nair, R., Bose, S., & Saini, T. D. 2018, Phys. Rev., D98, 023502, doi: 10.1103/PhysRevD.98.023502
  * Palmese et al. (2020) Palmese, A., et al. 2020, Astrophys. J. Lett., 900, L33, doi: 10.3847/2041-8213/abaeff
  * Payne & Thrane (2023) Payne, E., & Thrane, E. 2023, Phys. Rev. Res., 5, 023013, doi: 10.1103/PhysRevResearch.5.023013
  * Pierra et al. (2023) Pierra, G., Mastrogiovanni, S., Perriès, S., & Mapelli, M. 2023.  https://arxiv.org/abs/2312.11627
  * Planck Collaboration (2018) Planck Collaboration. 2018.  https://arxiv.org/abs/1807.06209
  * Ray et al. (2024) Ray, A., Magaña Hernandez, I., Breivik, K., & Creighton, J. 2024.  https://arxiv.org/abs/2404.03166
  * Ray et al. (2023) Ray, A., Magaña Hernandez, I., Mohite, S., Creighton, J., & Kapadia, S. 2023, Astrophys. J., 957, 37, doi: 10.3847/1538-4357/acf452
  * Riess et al. (2019) Riess, A. G., Casertano, S., Yuan, W., Macri, L. M., & Scolnic, D. 2019, Astrophys. J., 876, 85, doi: 10.3847/1538-4357/ab1422
  * Riess et al. (2022) Riess, A. G., Wenlong, Y., Macri Lucas, M., et al. 2022, Astrophys. J. Lett, 934, L7
  * Rinaldi & Del Pozzo (2021) Rinaldi, S., & Del Pozzo, W. 2021, Monthly Notices of the Royal Astronomical Society, 509, 5454, doi: 10.1093/mnras/stab3224
  * Sadiq et al. (2022) Sadiq, J., Dent, T., & Wysocki, D. 2022, Phys. Rev. D, 105, 123014, doi: 10.1103/PhysRevD.105.123014
  * Schutz (1986) Schutz, B. F. 1986, Nature, 323, 310, doi: 10.1038/323310a0
  * Soares-Santos et al. (2019) Soares-Santos, M., et al. 2019, Astrophys. J. Lett., 876, L7, doi: 10.3847/2041-8213/ab14f1
  * Suwa et al. (2018) Suwa, Y., Yoshida, T., Shibata, M., Umeda, H., & Takahashi, K. 2018, Mon. Not. Roy. Astron. Soc., 481, 3305, doi: 10.1093/mnras/sty2460
  * Talbot & Thrane (2018) Talbot, C., & Thrane, E. 2018, Astrophys. J., 856, 173, doi: 10.3847/1538-4357/aab34c
  * The LIGO Scientific Collaboration et al. (2021) The LIGO Scientific Collaboration, The Virgo Collaboration, & The KAGRA Collaboration. 2021, GWTC-3: Compact Binary Coalescences Observed by LIGO and Virgo During the Second Part of the Third Observing Run — Parameter estimation data release, Zenodo, doi: 10.5281/zenodo.5546663
  * The LIGO Scientific Collaboration et al. (2023) The LIGO Scientific Collaboration, Virgo Collaboration, & KAGRA Collaboration. 2023, GWTC-3: Compact Binary Coalescences Observed by LIGO and Virgo During the Second Part of the Third Observing Run — O1+O2+O3 Search Sensitivity Estimates, Zenodo, doi: 10.5281/zenodo.7890398
  * Tiwari (2021) Tiwari, V. 2021, Classical and Quantum Gravity, 38, 155007, doi: 10.1088/1361-6382/ac0b54
  * Tiwari (2022) —. 2022, The Astrophysical Journal, 928, 155, doi: 10.3847/1538-4357/ac589a
  * Tiwari & Fairhurst (2021) Tiwari, V., & Fairhurst, S. 2021, The Astrophysical Journal Letters, 913, L19, doi: 10.3847/2041-8213/abfbe7
  * Vitale et al. (2020) Vitale, S., Gerosa, D., Farr, W. M., & Taylor, S. R. 2020, doi: 10.1007/978-981-15-4702-7_45-1
  * Woosley (2017) Woosley, S. E. 2017, Astrophys. J., 836, 244, doi: 10.3847/1538-4357/836/2/244
  * Woosley (2019) Woosley, S. E. 2019, ApJ, 878, 49, doi: 10.3847/1538-4357/ab1b41
  * Woosley et al. (2002) Woosley, S. E., Heger, A., & Weaver, T. A. 2002, Rev. Mod. Phys., 74, 1015, doi: 10.1103/RevModPhys.74.1015
  * Ye & Fishbach (2022) Ye, C., & Fishbach, M. 2022, Astrophys. J., 937, 73, doi: 10.3847/1538-4357/ac7f99