# A Trifecta of Modelling Tools: Navigating COMPAS, Shark and Bayesian
Inference for Binary Black Hole Model Selection

Liana Rauf1,2, Cullan Howlett1,2, Simon Stevenson2,3, Jeff Riley2,4
1School of Mathematics and Physics, The University of Queensland, Brisbane,
QLD 4072, Australia
2OzGrav: The ARC Centre of Excellence for Gravitational Wave Discovery,
Australia
3Centre for Astrophysics and Supercomputing, Swinburne University of
Technology, Mail H39, PO Box 218, VIC 3122, Australia
4School of Physics and Astronomy, Monash University, Clayton VIC 3800,
Australia E-mail<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

Gravitational waves (GWs) have revealed surprising properties of binary black
hole (BBH) populations, but there is still mystery surrounding how these
compact objects evolve. We apply Bayesian inference and an efficient method to
calculate the BBH merger rates in the Shark host galaxies, to determine the
combination of COMPAS parameters that outputs a population most like the GW
sources from the LVK transient catalogue. For our COMPAS models, we calculate
the likelihood with and without the dependence on the predicted number of BBH
merger events. We find strong correlations between hyper-parameters governing
the specific angular momentum (AM) of mass lost during mass transfer, the
mass-loss rates of Wolf–Rayet stars via winds and the chemically homogeneous
evolution (CHE) formation channel. For the AM loss, we also find opposing
trends in the likelihood with and without the dependence on the number of GW
events. We conclude that analysing the marginalised and unmarginalised
likelihood is a good indicator of whether the population parameters
distribution and number of observed events reflect the LVK data. In doing so,
we see that the majority of the models preferred in terms of the population-
level parameters of the BBHs greatly overpredict the number of events we
should have observed to date. Looking at the smaller number of models which
perform well with both likelihoods, we find that those with no CHE, AM loss
occurring closer to the donor during the first mass-transfer event, and/or
higher rates of mass-loss from Wolf–Rayet winds are generally preferred by
current data. We find these conclusions to be robust to our choice of
selection criteria. However, we do identify several notable events that all of
our models struggle to reproduce, and that, even when incorporating selection
effects into our models to match our selection of LVK data, COMPAS tends to
underpredict the proportion of events with chirp mass $>50M_{\odot}$ and
$z>0.5$, leaving avenues for future investigations.

###### keywords:

gravitational waves – black hole physics – stars: evolution – galaxies: star
formation – methods: analytical – methods: numerical – methods: statistical

††pagerange: A Trifecta of Modelling Tools: Navigating COMPAS, Shark and
Bayesian Inference for Binary Black Hole Model Selection–C

##  1 Introduction

The detection of gravitational waves (GWs) have prompted the astrophysical
community to revise the models for binary black hole (BBH) formation and
evolution. There are currently two main formation channels proposed:

  1. 1.

Isolated: The formation of merging BBHs through the evolution of isolated
binaries. These stars are gravitationally bound since birth and are not
perturbed by other stars or compact objects.

  2. 2.

Dynamical: BBH mergers occurring in dense environments via three-body
encounters followed by energy exchange and hardening of the binary. Some
examples of dense environments are young star, nuclear or globular clusters,
or in accretion disks of the central supermassive BHs in galaxies.

GWs can be used to gain insight into the poorly constrained aspects of massive
stellar binary evolution and help in distinguishing between binary formation
channels. The detailed origins and formation channels for binary systems are
still unknown, particularly what fraction of these systems come from isolated
and dynamical channels.

There are also aspects of stellar evolution that are poorly understood, such
as the common envelope phase, mass transfer efficiency and angular momentum
loss. This information is critical for understanding the population of BBHs
observed by LIGO, Virgo and KAGRA (LVK), including their mass and redshift
distributions, and to inform us on what is required for future detectors
(Mapelli, 2019, 2020; Mandel & Broekgaarden, 2022; Arca Sedda et al., 2023).
Combining the data from GWTC-1 (Abbott et al., 2019a), GWTC-2 (Abbott et al.,
2021a) and GWTC-3 (Abbott et al., 2023b), there is a total of 90 confident
($p_{\rm astro}>0.5$) GW event candidates observed to date, including 85 BBHs,
2 binary neutron stars (BNSs), and 3 black-hole-neutron-stars (BHNSs). This
growing population of GWs, including those to be released from the ongoing
fourth LVK observing run (O4), provides key information on the source masses,
spins and local merger rates (Abbott et al., 2023a) of binary compact objects,
which can be used to answer some of these questions.

However, there are many formation pathways that could lead to the same
individual source properties. Simulations of ensembles of events, produced by
varying the physical assumptions, can provide better insight into how these
binary systems evolve and merge over time. When combined with the full set of
data observed to date, they can potentially break the degeneracy between
uncertainties in stellar evolution parameters and BBH populations. They also
enable us to produce merger rates that can be compared to observations. Mandel
& Broekgaarden (2021) provide a comprehensive review of the local merger rates
for the isolated and dynamical formation channels. The studies in the review
combine population synthesis codes with different assumptions for the
metallicity dependent star formation rate (SFR) density, which plays a crucial
role in the evolution of massive binaries (Neijssel et al., 2019; Chruślińska,
2024). Broekgaarden et al. (2022a) utilises COMPAS (Stevenson et al., 2017;
Vigna-Gómez et al., 2018; Team COMPAS: Riley, J. et al., 2022) to explore the
impact of stellar evolution parameters and metallicity specific star
formation, $\mathcal{S}(z,Z)$, on the merger rate, chirp mass and coalescence
time distributions. They found that uncertainties in the binary evolution
models only led to variations in the BBH merger rate of up to a factor 2, with
the BBH merger rates most sensitive to $\mathcal{S}(z,Z)$. Santoliquido et al.
(2021) and Neijssel et al. (2019) also come to similar conclusions. Many of
these studies choose to explore the variations in binary evolution models over
a grid of metallicities, and generate populations for each grid point. Some
studies have shown that population III stars can also contribute to the total
BBH merger rate density (Santoliquido et al., 2024; Tanikawa, 2024). These
studies emphasise the need to further study the assumptions in binary
evolution models to understand the impact on the merger rate, specifically the
combined impact of model variations (Stevenson & Clarke, 2022).

There have also been various studies proposing methods for constraining
different binary evolution models. Wong et al. (2023) applies “backward"
modelling of each GW event to the progenitors to avoid fixed astrophysical
assumptions, by allowing variations in the hyper-parameters. However, they
note that their methods become less efficient as the dimensionality of the
parameter space increases. Barrett et al. (2018) uses Fisher information
matrices to quantitatively analyse how the observed BBH population will
constrain the COMPAS models, given the chirp mass distribution and accounting
for the metallicity specific cosmic SFR, selection effects and measurement
uncertainty. Riley & Mandel (2023) uses artificial neural networks to generate
an interpolant to calculate the merger detection rates for double compact
objects for varying prescriptions of the metallicity specific SFR. This is
used to create GW populations for each model, where the chirp mass and
redshift distribution are compared using Bayesian analysis to the LVK data.
Delfavero et al. (2023) calculates a joint likelihood for binary evolution
models. They fit normal distribution fits to the likelihood grids and
implement Gaussian process regression to interpolate between models, the
detected number of events and the likelihood. This is also used to quantify
the systematic uncertainty in the parameters. However, they emphasise that the
inference on a limited subspace of models could be biased. These studies
highlight the high computational costs for multi-dimensional parameter
studies. However, as more observations are made, the uncertain astrophysics
will produce correlated impacts in the binary population, so there is a
crucial need for more parameters covering a large range in model inference.
Realistically, all formation channels (within the isolated and dynamical
channels) must be considered to obtain a parameter space representative of the
true BBH population (Zevin et al., 2021; Cheng et al., 2023a)

To build on this body of work, we develop here a new method of combining a
more accurate, non-parametric star-formation and metallicity history of the
Universe with a wider range of different massive binary evolution models. We
then use these models alongside data from GWTC-1 (Abbott et al., 2019a),
GWTC-2 (Abbott et al., 2021a) and GWTC-3 (Abbott et al., 2023b) to identify
the most likely avenues for generating binary compact objects, and to rule out
those that are unable to produce the distributions of objects seen to date. We
do this by improving on Rauf et al. (2023), where we combined runs from the
COMPAS rapid population synthesis code, for formation and evolution of BBH
formed through isolated binary evolution, with a star formation rate and
metallicity history for the observed host galaxies predicted by the semi-
analytic galaxy evolution code Shark (Lagos et al., 2018; Lagos et al., 2019).
In our original work, we predicted (among other things) the BBH merger rate as
a function of redshift for four COMPAS runs, where we varied the remnant mass
prescription. We found the merger rate ranged over several magnitudes and we
concluded that the Mandel & Müller (2020) remnant mass prescription was the
best fit out of the four models to the local merger rate currently predicted
using the BBH mergers from GWTC-3. However, we also found that cross-matching
between the Shark galaxies and the COMPAS population was prohibitively time-
consuming for exploring larger portions of our parameter space. Since we
explored only rates, we were also not able to confidently rule out any models
that do not align with the GW population-level parameters, and did not test
the impact of varying sub-parameters of these models such as those covering
other aspects of binary evolution rather than just the mass of the remnant
compact objects.

So, in this paper, we substantially extend on our work in Rauf et al. (2023)
and explore an alternative and efficient method of predicting the BBH merger
rate and population level parameters such as the binary masses and redshift
distributions. We do so by rescaling a high-fidelity single run of our
original simulation, varying some of the most relevant COMPAS input
parameters. This allows us to test in this work $>300$ different models for
population synthesis, varying ingredients/model assumptions such as the common
envelope efficiency, mass transfer accretion efficiency, mass-loss rates of
Wolf–Rayet stars, chemically homogeneous evolution model and specific angular
momentum loss. For each of these we evaluate the ability of the model to
reproduce the observed distribution of black hole masses and redshifts, using
a Bayesian framework to account for selection effects. We also explore the
impact on our findings with and without including the merger rate as a fitting
parameter (i.e., including a comparison of the model-predicted and observed
merger rates in evaluating our model posteriors). Compared to previous studies
this work also has the improvement of incorporating a basic estimate of
sampling error in the simulated model merger rates and posteriors. This allows
us to set a threshold for what we consider to be tangible differences in model
performance, and also demonstrates more generally the importance of further
investigation into the uncertainties in the typical ‘simulation-based
inference’ used in GW population studies.

This paper is structured as follows. In Section 2, we summarise recent
progress in our understanding of various massive binary evolution models and
parameters, motivating the range of models we explore in this work. In Section
3, we show that starting with a fiducial model, where the merger rate in a set
of simulated galaxies has been evaluated by carefully cross-matching between
COMPAS and Shark, we can rescale the simulation to produce the time and
metallicity dependent GW merger rate for another model where the COMPAS
parameters have been altered. In Section 4, we implement Bayesian inference to
determine which of these COMPAS runs output a BBH population that best matches
the current observed population. We do so by comparing the chirp mass, mass
ratio and redshift distributions from each model with the LVK posterior sample
distributions. We then validate the total number of theoretical and detected
GW events expected for the models and the models posteriors by comparing to
the values obtained using the four ‘brute-force’ evaluations from Rauf et al.
(2023). We repeat the method for a wide range of new COMPAS runs in Section 5,
identifying the best and worst fit models to the data and studying the impact
of different selection effects and analysis choices on the conclusions we can
draw about stellar evolution. We conclude in Section 6. In this work, we use
the Shark cosmology $H_{0},\Omega_{m}=\\{67.8,0.308\\}$ where necessary.

##  2 COMPAS: Overview of Models

In the following section, we detail the physics and the implementation in
COMPAS of the common envelope, mass transfer, mass loss via winds, angular
momentum loss and formation channels such as the isolated and chemically
homogeneous evolution. We review the literature on the assumptions and
parameters associated with these processes that can impact the evolution of
the binary stars and, hence, the number and source properties of GWs for each
COMPAS model. We conclude with a summary of the big questions surrounding
these prescriptions and their limitations in the context of BBH population
modelling.

###  2.1 Isolated formation channels

In this paper, we focus on the isolated formation channel. This is the most
commonly studied channel. It involves massive stars formed on the zero-age-
main-sequence (ZAMS). We describe, one of the formation channels, the
classical isolated binary evolution via a common envelope phase as follows:

  1. 1.

The primary star in the binary reaches the end of main sequence (MS)
evolution. Hydrogen burning in the core causes the outer stellar envelope to
expand past the Roche lobe and transfer mass to the secondary star.

  2. 2.

Mass transfer from the primary stops and leaves a Wolf-Rayet star. This then
reaches the end of its life, undergoes a supernova and forms a compact object.
The system is now an X-ray binary.

  3. 3.

The secondary evolves and transfers mass back to the primary (compact) object.
This process is dynamically unstable and forms a common envelope, which
encompasses the compact object.

  4. 4.

Orbital energy released due to the compact object rapidly in-spiralling can
remove the common envelope.

  5. 5.

The secondary star evolves and eventually also forms a compact object.

  6. 6.

The double compact object (DCO) system slowly in-spirals due to orbital energy
loss via GWs.

However, there are many stages within this formation channel where the binary
may not survive, or the events can subtly differ from the above prescription
leading to dramatically different outcomes. The system may not survive the
first supernova explosion or the common envelope may be not be ejected
resulting in the compact object merging with the secondary core before a DCO
is formed. Alternatively, the second mass transfer could be stable and avoid
the common envelope (e.g., van den Heuvel et al., 2017; Neijssel et al., 2019;
van Son et al., 2022a; Gallegos-Garcia et al., 2021; Marchant et al., 2021).
The wind mass-loss rates and remnant prescription can also impact the orbital
separation, as well as the final masses. There are various circumstances and
factors that can affect the separation of DCO systems, which is strongly
correlated with the delay time and the merger rate. On one hand, this makes GW
sources unique probes of cosmic history (Neijssel et al., 2019; Mapelli, 2021;
Mandel & Farmer, 2022; Costa et al., 2023), but on the other, leads to
difficulties in interpreting or identifying the most important/critical points
in the evolution of the system.

As an example, Stevenson et al. (2017) used COMPAS to show that the first
three GW events detected by LIGO: GW150914, GW151226 and GW151012, can be
explained by the classical isolated channel with mass transfer and common
envelope phase, and characterised their progenitor metallicity and masses.
However, Dorozsmai & Toonen (2022) found that the stable channel, consisting
of two stable mass transfer phases, dominates the GW sources — formation of GW
sources via a radiative common envelope channel is not possible for low mass
transfer efficiency ($f_{\rm MT}\leq 0.3$) due to the mass ratio being below
the critical mass ratio ($q_{\rm crit}=M_{d}/M_{a}$, where $M_{a}$ is the
acretor mass and $M_{d}$ is the donor mass). Picco et al. (2023) similarly
found that stable mass transfer, where mass and angular momentum is lost
through the L2 point, is a possible formation channel for BBHs. Comparing the
two channels, van Son et al. (2022a) found that the common envelope channel
tends to dominate BBH mergers with primary masses $<30M_{\odot}$ binaries with
short delay times ($<1$Gyr). They also found that the common envelope phase is
rare for massive stars due to heavy mass loss preventing them from filling
their Roche lobe. As such, the stable mass transfer channel, which widens the
binary orbit as mass is lost from the system during secondary mass transfer,
forms primary BH masses $>30M_{\odot}$ and systems with long delay times
($>1$Gyr).

####  2.1.1 Common envelope evolution ($\alpha_{\rm CE}$) and Mass transfer
efficiency ($f_{\rm MT}$)

Given the above, in the formation of isolated binary stars, uncertainties
remain in the common envelope (CE) phase and mass transfer process, which play
a key role in the final component and orbital properties.

The mass transfer efficiency is the fraction of mass lost by the donor star
and accreted by the companion star. Team COMPAS: Riley, J. et al. (2022) uses
a simplified approach for the mass transfer. They parameterise the fraction of
mass lost by the donor to the accretor with $f_{\rm MT}$,

$\displaystyle\dot{M}_{a}=-f_{\rm MT}\dot{M}_{d},$ (1)

with $0\leq f_{\rm MT}\leq 1$. This fraction depends on the properties of the
binary prior to mass transfer and varies over time. We use a fixed fraction
prescription in COMPAS, so the fraction of mass accreted is assumed to be
constant over the mass transfer phase. $f_{\rm MT}=1$ corresponds to fully
conservative mass transfer, so all the mass lost from the donor is accreted.
$f_{\rm MT}<1$ is non-conservative and $f_{\rm MT}=0$ is fully non-
conservative. In the default setting, the mass transfer efficiency is
determined by the ratio of the Kelvin–Helmholtz thermal timescale of the two
stars,

$\displaystyle\tau_{\rm KH}=\frac{E_{\rm int}}{L}.$ (2)

This is the time taken for the star’s internal energy, $E_{\rm int}$, to be
radiated at it’s current luminosity, $L$. The maximum accretion rate onto a
stellar accretor is,

$\displaystyle\dot{M}_{a}=\frac{CM_{a}}{\tau_{\rm KH,a}},$ (3)

where the default factor $C=10$ is assumed to account for the expansion of the
star due to mass loss (Paczyński & Sienkiewicz, 1972; Neo et al., 1977; Hurley
et al., 2002; Schneider et al., 2015). This is calculated for every mass
transfer event and is non-constant. Accretion onto the compact objects, such
as black holes and neutron stars, is assumed to be Eddington limited such
that,

$\displaystyle\dot{M}_{a}=\dot{M}_{\rm edd}=1.5\times
10^{-8}\frac{R_{*}}{10\rm km}\frac{M_{\odot}}{\rm yr},$ (4)

where the black hole radius, $R_{*}$, is the Schwarzschild radius. In the
unstable mass transfer scenario, for simplicity, the CE phase is parameterised
by COMPAS in the context of energy balance with the $\alpha_{\rm CE}$
formalism (Xu & Li, 2010). The CE efficiency, $\alpha_{\rm CE}$, describes the
efficiency of converting the orbital energy into kinetic energy to unbind the
common envelope. This unbinding energy is defined as,

$\displaystyle E_{\rm bind}=\alpha_{\rm
CE}\left(-\frac{GM_{d}M_{2}}{2a_{i}}+\frac{GM_{\rm
core}M_{2}}{2a_{f}}\right),$ (5)

where G is the gravitational constant, $M_{d}$ is the mass of the pre-CE
donor, $M_{2}$ is the mass of the secondary, $M_{\rm core}$ is the mass of the
donor’s core, $a_{i}$ and $a_{f}$ are the pre- and post-CE orbital
separations. Increasing $\alpha_{\rm CE}$ increases the unbinding energy, so
the CE can be ejected more efficiently at wider binary separations. Too small
values of $\alpha_{\rm CE}$ would require the binary to shrink in order to
obtain enough unbinding energy (Tutukov & Yungelson, 1979; Webbink, 1984;
Ivanova et al., 2013; Hirai & Mandel, 2022; Stevenson & Clarke, 2022). de Kool
(1990) proposed an alternative version of the binding energy, introducing a
new parameter, $\lambda$, to characterise the central concentration of the
donor’s envelope,

$\displaystyle E_{\rm bind}=-\frac{GM_{1}M_{\rm env}}{\lambda a_{i}r_{L}},$
(6)

where $M_{\rm env}=M_{1}-M_{\rm core}$ is the mass of the envelope and
$r_{L}=R_{L}/a_{i}$ is the ratio of the Roche lobe radius and pre-CE orbital
separation. The post-CE orbital separation is determined by substituting Eq 6
into Eq 5,

$\displaystyle\frac{a_{f}}{a_{i}}=\frac{M_{\rm
core}M_{2}}{M_{1}}\frac{1}{M_{2}+2M_{\rm env}/\alpha_{\rm CE}\lambda r_{L}}.$
(7)

Due to lack of understanding and convenience in population synthesis
calculations, $\alpha_{\rm CE}$ and $\lambda$ are set to constants, but should
be varied depending the stellar and binary parameters (Xu & Li, 2010).

Progenitors of GW sources provide an alternative path for comparison between
models and observations. Zuo & Li (2014) uses the binary population synthesis
code by Hurley et al. (2000) and Hurley et al. (2002) to study high-mass X-ray
binaries (HMXB) and the correlation between their apparent luminosity and
displacement from star clusters. A high value of $\alpha_{\rm CE}=0.8-1.0$ is
preferred when comparing to observations by Kaaret et al. (2004). However, the
ambiguity in the definition of the donor core mass boundary introduces a large
uncertainty (almost two orders of magnitude), which translates directly to the
expected value of $\alpha_{\rm CE}$. Romero-Shaw et al. (2023) uses COMPAS to
study black hole high-mass X-ray binaries (BH-HMXB). They found for the case
where $f_{\rm MT}\leq 0.1$ (and $f_{\rm WR}=0.1$, explained in Section 2.1.2),
that $8\%$ of BH-HMXBs end up as BBHs or BHNSs. de Mink et al. (2007) attempts
to fit the mass transfer efficiency fraction for massive close binaries. They
find a wide spread in $f_{\rm MT}$ and conclude that a single value is not
valid. Therefore, more electromagnetic observations are required to constrain
modelling degeneracies of massive binaries.

There have been many studies on the mass transfer and common envelope phase.
Many of these have explored how varying these parameters affect the BBH merger
rate, the BH mass distribution range and which binary formation channels
dominate these BBH mergers (Santoliquido et al., 2021; Grichener, 2023; van
Son et al., 2020, 2022a, 2022b; Marchant et al., 2021). To infer the formation
channels of BBH mergers, many studies further explore the connection between
the models and data (Cheng et al., 2023b). Bouffanais et al. (2021) find the
preferred parameter range for the common envelope and mass transfer efficiency
is $[4,7]$ and $[0.7,0.8]$ respectively. They also showed that increasing the
common envelope efficiency increases the uncertainty on the local BBH merger
rate density. However studies such as Bavera et al. (2021),Klencki et al.
(2021) and Cui & Li (2023) find there is a preference for lower $\alpha_{\rm
CE}$ ($0.2<\alpha_{\rm CE}<1$, $\alpha_{\rm CE}\leq 0.7$ and $\alpha_{\rm
CE}=0.5,1$ respectively). Klencki et al. (2021) further questions the validity
of the CE channel for producing massive BHs. They find that successful
ejection of the CE is possible if the donor is a massive convective-envelope
giant.

We briefly mention the various studies that explore the uncertainties in
parameters. But another uncertainty to note is in how these parameters are
implemented in population synthesis codes. Gallegos-Garcia et al. (2021)
highlights the uncertainty in stellar evolution tracks by comparing the mass
transfer and common envelope treatments in the population synthesis codes MESA
and COSMIC. The COSMIC models primarily form BBH mergers through the common
envelope phase, while the stable mass transfer formation channel is preferred
with MESA. Increasing $\alpha_{\rm CE}$ also increases the variation in the
BBH merger rates between the models, due to variations in merger times.
Agrawal et al. (2022a) and Romagnolo et al. (2023) have compared population
synthesis codes to explore the uncertainty in radial expansion of massive
stars and the variation in the remnant mass.

####  2.1.2 Wolf–Rayet Winds ($f_{\rm WR}$)

Stellar winds are outflows of gas from the atmosphere of a star. Wolf–Rayet
stars are evolved stars with high mass loss rates, where the hydrogen envelope
of the stars are removed. To account for uncertainties in the mass-loss rates
of Wolf–Rayet stars, in COMPAS, the Wolf–Rayet multiplier, $f_{\rm WR}$, is
used to parameterise the mass loss rate (Belczynski et al., 2010),

$\displaystyle\frac{dM}{dt}=f_{\rm WR}\times
10^{-13}L^{1.5}\left(\frac{Z}{Z_{\odot}}\right)^{m}\rm M_{\odot}yr^{-1},$ (8)

where $M$ is the mass of the helium star, $L$ is the luminosity, $Z$ is the
metallicity, $Z_{\odot}=0.014$ (Asplund et al., 2009) and $m=0.86$ (Vink & de
Koter, 2005). Wolf–Rayet stars have been studied to understand their explosion
mechanisms and determine if they are progenitors of gamma ray bursts (Detmers
et al., 2008). More recently, Wolf–Rayet stars have been used to study GW
sources (Vink & Harries, 2017; Bogomazov et al., 2018). Decreasing $f_{\rm
WR}$ leads to an increase in the BH masses (Dorozsmai & Toonen, 2022; Romero-
Shaw et al., 2023).

HMXBs are usually ideal candidates to study and constrain the effects of
stellar mass loss in interacting binaries, as radiation driven winds are
expected for giant stars. Schrøder et al. (2021) finds that the orbital
evolution of HMXBs depends on the wind velocity and mass ratio. Slow winds and
high mass objects results in the binary shrinking, due to drag forces. In the
case where there is no drag, the binary separation widens. GWs can provide an
indirect approach to study the separation of binaries, as their progenitors
include HMXBs. However, Liotine et al. (2023) found that no currently observed
HMXB is predicted to form a merging BBH with high probability. They use COSMIC
(Breivik et al., 2020) to conclude that the BHs in the detectable HMXBs are
not likely to have masses $>35M_{\odot}$, the detectable HMXBs will not merge
as BBHs within the Hubble time and the distribution of redshifts, masses and
metallicities of detectable HMXBs and BBHs are different. Fishbach & Kalogera
(2022) argue that discrepancies in the BH masses in BBHs and HMXBs are due to
GW observational selection effects.

Neijssel et al. (2021) and Romero-Shaw et al. (2023) find a preference for low
$f_{\rm WR}$ ($<0.2$) for HMXB systems such as Cygnus X-1 (Miller-Jones et
al., 2021). Stevenson & Clarke (2022) explores the correlated impact of model
parameters on the merger rate. They identified a correlation between the Wolf-
Rayet factor and metallicity, as increasing both parameters increases mass
loss via stellar winds. They found most models that match the current
observations have $f_{\rm WR}>1$. This contradicts the previously mentioned
studies due to COMPAS over-predicting the rates and masses of BBHs formed
through chemically homogeneous evolution (see Section 2.2).

Riley et al. (2021) studies the CHE binaries as a formation channel for
various $f_{\rm WR}$ with COMPAS. They found a correlation between the delay
time and $f_{\rm WR}$, as increasing mass loss leads to wider binaries. This
results in the merger rate decreasing with increasing $f_{\rm WR}$. For
$f_{\rm WR}=1$, the CHE channel contributes $20$ $\rm Gpc^{-3}yr^{-1}$ to the
total rate of $50$ $\rm Gpc^{-3}yr^{-1}$, which is an overestimation compared
to GWTC-3. Due to the short delay times, the BBH merger rate peaks at higher
redshifts ($z\approx 3-4$) compared to the peak SFR.

####  2.1.3 Specific angular momentum loss ($f_{\gamma}$)

Mass transfer from a massive donor to a lower-mass accretor adds specific
angular momentum to the transferred mass, so conservation of angular momentum
requires that binary’s orbit shrinks. If mass transfer is non-conservative,
angular momentum, $J_{\rm orb}$, is carried away by mass lost from the system.
The change in angular momentum is parameterised by $\gamma$,

$\displaystyle\frac{\dot{J}_{\rm orb}}{\dot{M}_{d}}=\gamma\frac{J_{\rm
orb}}{M_{d}}.$ (9)

The orbital evolution of the binary is given by,

$\displaystyle\frac{\dot{a}}{a}$
$\displaystyle=-2\frac{\dot{M}_{d}}{M_{d}}\left[1-f_{\rm
MT}\frac{M_{d}}{M_{a}}-(1-f_{\rm
MT})\left(\gamma+\frac{1}{2}\right)\frac{M_{d}}{M_{d}+M_{a}}\right]$
$\displaystyle=-2\frac{\dot{M}_{d}}{M_{d}}\left[1-\frac{f_{\rm
MT}}{q}-(1-f_{\rm MT})\left(\gamma+\frac{1}{2}\right)\frac{1}{1+q}\right].$
(10)

Willcox et al. (2023) introduces a new parametrisation, $a_{\gamma}$,

$\displaystyle a_{\gamma}=a_{\rm acc}+f_{\gamma}(a_{\rm L2}-a_{\rm acc}),$
(11)

which is the effective decoupling radius between the accretor, $a_{\rm acc}$,
and the L2 point, $a_{\rm L2}$. COMPAS default model assumes isotropic re-
emission of matter during non-conservative stable mass transfer
($f_{\gamma}=0$). In this model, mass lost from the donor is transported near
the accretor and ejected via fast isotropic winds. Hydrodynamical simulations
indicate that this mass is lost between the accretor and L2 point (MacLeod et
al., 2018a, b). $f_{\gamma}$ determines the point where this angular momentum
is lost. If the donor star is growing rapidly with respect to its Roche Lobe,
the system will lose mass through the outer Lagrangian point (L2), which
causes the orbit to shrink. This corresponds to $f_{\gamma}=1$. Angular
momentum loss directly from the accretor corresponds to $f_{\gamma}=0$ (Costa
et al., 2023; Romero-Shaw et al., 2023).

Lu et al. (2023) and Marchant et al. (2021) found that their models predicted
high mass transfer rates, which lead to mass loss through the L2 point.
MacLeod et al. (2018a); MacLeod & Loeb (2020) argue that the specific angular
momentum is ejected between the accretor and the L2 point. Willcox et al.
(2023) uses COMPAS to study the stable mass transfer channel. They found that
increasing $f_{\gamma}$ leads to a reduction in the number of systems that
undergo stable mass transfer. Increasing $f_{\rm MT}$ while decreasing
$f_{\gamma}$ increases the likelihood of stable mass transfer. For
$f_{\gamma}\gtrsim 0.7$, nearly all systems experience common envelope
evolution.

van Son et al. (2022b) studies the impact of $f_{\gamma}$ on the local BBH
merger rate and the primary BH mass distribution. for $f_{\gamma}=0.75$ and
$f_{\gamma}=1$, there are no BBHs formed through the stable mass transfer
channel. The merger rate drops from $141\rm Gpc^{-3}yr^{-1}$ with
$f_{\gamma}=0.25$ to $16\rm Gpc^{-3}yr^{-1}$ with $f_{\gamma}=0.5$, which
implies the stable mass transfer contributes strongly to the merger rate for
$f_{\gamma}<0.5$. As the low mass end of the primary BH mass distribution is
dominated by the stable mass transfer channel, as $f_{\gamma}$ increases, the
minimum primary BH mass increases.

###  2.2 Chemically homogeneous evolution

In the chemically homogeneous evolution (CHE) channel, the stars do not
expand. They are metal poor and rapidly rotating in tight binaries. The
rotation efficiently transports hydrogen to the core and helium to the
envelope. This evolves and forms a Wolf–Rayet star, which contracts rather
than expands. Wind-driven mass loss does not widen the binary significantly
for low metallicity systems and there is no mass transfer, with the exception
of massive over-contact binaries (Marchant et al., 2016). Highly spinning DCOs
would produce GWs in the age of the Universe (Nathaniel, 2017; Mandel &
Farmer, 2022; Costa et al., 2023). de Mink et al. (2009, 2010) introduces the
concept of CHE to justify short period massive binaries such as the HMXB M33
X-7 (Orosz et al., 2007).

Mandel & de Mink (2016) found that the CHE channel significantly contributes
to the BBH merger rate ($10.5\pm 0.5$ $\rm Gpc^{-3}yr^{-1}$ at redshift 0).
These systems typically have a total mass ranging between $50-110$
$M_{\odot}$. Marchant et al. (2016) explores the parameter space for forming
CHE binaries. They found that at low metallicity, BBHs will merge within the
age of the Universe due to weaker winds ($t_{\rm delay,min}=0.4$ Gyr for
$Z_{\odot}/50$). There is strong evidence for increasing chirp masses with
decreasing metallicity, with masses above the pair instability SN gap. If
these mergers were detected, they could probe the evolution of massive stars
in the early Universe. In addition to calculating the merger rate and mass
distribution for the CHE channel, du Buisson et al. (2020) finds that the
cosmological merger rate ($5.8$ $\rm Gpc^{-3}yr^{-1}$ at redshift 0) can vary
up to 40% for extreme deviations in the metallicity specific SFR.

In COMPAS, there are two CHE modes, ‘optimistic’ and ‘pessimistic’, which
refer to how a CH star evolves on the MS (Riley et al., 2021). In the
‘optimistic’ mode, if the star begins the main sequence with fast enough
rotation to evolve as a CHE star, it is assumed the star continues rotating
fast enough to remain CH throughout the MS and the spin is not checked. For
the CHE mode ‘pessimistic’ the spin is checked at every timestep, and if it
has slowed below the CH metallicity-dependent rotational frequency threshold,
it is immediately converted to a regular MS star and evolve to the end of the
MS (Team COMPAS: Riley, J. et al., 2022).

###  2.3 Remnant mass prescription

GWs can be used to infer the underlying mass distributions of BBHs. This
remnant mass distributions depend on the metallicity and the type of explosion
mechanism. Their formation can be described in three stages; stellar collapse
and bounce, convective engine, and post-explosion fallback. Massive stars
undergo a series of successful burning stages, building an iron core until
electron degeneracy pressure can no longer support this. Electron capture and
dissociation of the core elements into alpha particles accelerates the core
collapse. This collapse is halted when the core reaches nuclear densities and
nuclear forces cause the core to bounce and emit a shock outwards. This
results in a mixing of the hot and cold layers in the star, known as
convection. The convective engine provides the means to convert potential
energy released in the collapse into explosion energy in the form of a
supernovae (Fryer et al., 2012; Fryer et al., 2021). The commonly accepted
supernova explosion mechanism is neutrino driven, where some of the neutrinos
emitted during the core collapse are reabsorbed, powering an explosion (Janka
et al., 2007; Marek & Janka, 2009; Janka, 2012, 2017).

Fryer et al. (2012) provide analytic prescriptions for estimating the remnant
mass in population synthesis codes, based on at the time recent understanding
of supernova and gamma-ray burst explosions. In their model/fit, the remnant
mass depends on the core, fallback mass and explosion mechanism. They explore
two neutrino-driven explosion mechanisms– the rapid convection, which leads to
explosion in the 250 ms after bounce, and the delayed convection, which can
occur over a longer timescale. Fryer et al. (2022) altered the remnant mass
prescription from Fryer et al. (2012) to study the growth rate of convection.
They introduce a new parameter in their models denoted $f_{\rm mix}$, the
mixing growth time for convective instabilities. A fast growth time, $f_{\rm
mix}=4.0$, and a slow growth time, $f_{\rm mix}=0.5$, correspond to the rapid
and delayed models respectively from Fryer et al. (2012). Whether a star
undergoes a successful or failed supernovae provides clues on the final
remnant object. However, there are uncertainties surrounding these mechanisms
(Mapelli, 2021). GW detection from a core collapse SN can be distinguish
between explosion mechanisms (Powell et al., 2016, 2023)

Mandel & Müller (2020) develop probabilistic prescriptions for the remnant
mass based on results from 3D supernovae simulations in Müller et al. (2016).
Instead of a single value prescription, the BH mass branches over a range of
core masses, depending on the fallback strength. If the BH is not formed by
complete fallback, the mass follows as normal distribution. These models are
more realistic as they capture the stochasticity of stellar evolution.

Schneider et al. (2021) also applies the Mandel & Müller (2020) semi-analytic
SN prescription, but the resulting remnant mass prescription differs. The
model focuses more on the removal of the hydrogen-rich common envelope and the
evolution steps before core collapse. They also explore the amount of
fallback. Their models overestimate the BH masses by 10% if the BHs are not
formed by direct collapse. They do not account for any convective boundary
mixing beyond the core helium burning.

Schneider et al. (2023) uses the models in Schneider et al. (2021) and the
stellar evolution code MESA (Paxton et al., 2011) to explore the black hole
mass distribution and chirp mass distribution of BBH mergers. The ability of
stars to explode in the delayed neutrino-driven SN mechanism is characterised
by the compactness of the progenitor cores (O’Connor & Ott, 2011; Ertl et al.,
2016). This depends on the mass and radius, and a high compactness indicates
unsuccessful explosions and collapse into a BH. Thus, the chirp mass
distribution is bimodal and the model predicts a lack of chirp masses in the
range $9-13M_{\odot}$. Mass accretion can influence the final outcome of core
collapse and final remnant mass. Schneider et al. (2024) studies models with
MESA that allow for multiple mergers. They find BH masses in the pair
instability SN gap, which can explain events such as GW190521. They note that
for a more complete picture of stellar mergers, further detailed multi-
dimensional, hydrodynamic simulations and a systematic exploration of the
relevant parameter space is needed.

###  2.4 Summary

We outline the parameters and model assumptions varied in COMPAS in Table 1.
We note that the initial conditions and distributions for the binary systems
are set to default (except for the metallicity where we use Shark), as shown
in Table 2. Future studies may choose to vary these parameters.

Parameters | Value | Default
---|---|---
Remnant mass prescription | [MULLERMANDEL, SCHNEIDER2020, FRYER2012, FRYER2022] | MULLERMANDEL
Chemical homogeneous evolution | [NONE, PESSIMISTIC,OPTIMISTIC] | PESSIMISTIC
Mass transfer fraction ($f_{\rm MT}$) | [0.1,0.2,0.5,1] | 0.5
Common envelope efficiency ($\alpha_{\rm CE}$) | [0.1,1,10,100] | 1.0
L2 point angular momentum loss linear fraction ($f_{\gamma}$) | [0,0.25,0.5,0.75,1] | 0.5
Wolf-Rayet winds multiplier ($f_{\rm WR}$) | [0.1,0.5,1,5,10] | 1.0
Table 1: COMPAS hyper-parameters explored in this paper. Note that some
prescriptions are only used in conjunction with others. The mass transfer
fraction is only applied if the mass transfer accretion efficiency
prescription is set to FIXED_FRACTION. The default is THERMAL. Similarly, the
angular momentum loss linear fraction is only applied when the mass transfer
angular momentum loss prescription is MACLEOD_LINEAR. The default is ISOTROPIC
($f_{\gamma}=0$). Refer to Tables 4,5 and 6 in Appendix C for all relevant
prescriptions used for each model.

As discussed, we are aware of the limitations and approximations made in
COMPAS and other population synthesis codes, that may cause discrepancies in
the outputs. We summarize these limitations and the key questions below:

  1. 1.

The binding energy formalism utilises a fixed $\alpha_{\rm CE}$, but is this a
valid approach? Hirai & Mandel (2022) proposes a two-step CE process, which
depends on the evolutionary stage of the donor and the companion’s mass. They
find that the predicted post-CE separation is wider compared to the
$\alpha_{\rm CE}$ formalism. This may pose a problem for BBHs merging within
Hubble time. Picker et al. (2024) has recently implemented fits for the
convective envelope mass and associated values of $\lambda$ to be used in
conjunction with the two-step CE formalism. They find the effective combined
value of $\alpha_{\rm CE}\lambda$ can range from $\sim 10^{-3}$ to $10^{3}$
when converted to the $\alpha$-formalism.

  2. 2.

Are the masses in COMPAS overestimated and how does this impact the CHE
binaries? CHE stars evolve from the hydrogen-ZAMS directly to the helium-ZAMS.
When a CHE star is close to the end of the MS, it is already rich in helium,
and may start to resemble a helium/ Wolf–Rayet star. However, COMPAS
implements the same wind loss prescription throughout the MS, instead of a
specific Wolf–Rayet wind prescription with higher mass loss rates. This change
is applied to the MESA CHE modelling (Marchant et al., 2016; du Buisson et
al., 2020; Marchant et al., 2023).

  3. 3.

Mass loss via radiation-driven winds is still an uncertainty in massive star
progenitor models and accurate predictions are necessary for predicting their
evolutionary paths. Björklund et al. (2023) derive radiation-driven mass-loss
from hot, massive stars depending on their fundamental stellar parameters and
implement this in MESA. They find that lower mass-loss rates potentially allow
for the creation of high-mass black holes even at higher metallicities
compared to previous models.

  4. 4.

At masses $>50M_{\odot}$, Hurley et al. (2000) smoothly extrapolate the Pols
et al. (1998) model grid up to $150M_{\odot}$, which is not well tested
outside of the prescribed mass range. Significant improvement on this
extrapolation method is far from trivial, as the evolution of high mass
binaries is still very uncertain. This is reflected by the discrepancies
between population synthesis codes in their stellar tracks (Agrawal et al.,
2020, 2022a, 2022b).

  5. 5.

There are many uncertainties in the mass transfer efficiency and angular
momentum loss, which influence the Roche lobe response to mass loss.
Uncertainties in how a donor star responds to mass loss, combined with the
Roche lobe response uncertainties contributes to the lack of understanding
around the mass transfer stability boundary as a function of donor mass and
radius (Ge et al., 2010; Ge et al., 2015, 2020; Temmink et al., 2023).

  6. 6.

Is the angular momentum loss occurring closer to the accretor or the L2 point?
And should the implementation of $f_{\gamma}$ vary for degenerate and non-
degenerate accretors? When performing the analysis in this paper, split
functionality for AM loss from degenerate and non-degenerate accretors was yet
to be implemented in COMPAS. Physically, it is more accurate for $f_{\gamma}$
to be implemented differently for the primary and secondary star, particular
when the accretor is a compact object. It is not clear how this would impact
our current results, where we have only considered the angular momentum loss
near the non-degenerate object. Gallegos-Garcia et al. (2023) derive an
analytic model for mass loss from the binary that is informed by observations
of X-ray binary accretion disks and hydrodynamical simulations, which
considers angular momentum loss from a disk wind around a compact object
accretor. They find their model predicts more angular momentum loss compared
to the standard case of isotropic re-emission and the orbital separation of a
binary can be considerably affected by mass loss via a disk wind.

  7. 7.

COMPAS has two modes; single stellar evolution (SSE) and binary stellar
evolution (BSE). The Hurley et al. (2000, 2002) SSE algorithms are implemented
in the BSE mode at all evolutionary phases and provides the stellar attributes
for both components of the binary. However, the fitting formulae are not well
behaved outside of the metallicity range $Z\in[10^{-4},0.03]$.

  8. 8.

The SSE tracks in COMPAS are based on non-rotating stars and do not account
for effects on stellar evolution due to mild rotation. Thus, the modelling of
rotational mixing, core overshooting (convection beyond the Schwarzschild
criterion) and donor stripping is not accurate. Many studies had investigated
the impact of these parameters on stellar structures using MESA and PARSEC
(Claret & Torres, 2018; Costa et al., 2019; Laplace et al., 2020; Temaj et
al., 2024).

##  3 Predicting the GW merger rate for different COMPAS models

In this section we demonstrate our method for rescaling the merger rate from a
fiducial, carefully cross-matched, simulation of galaxy formation (via Shark)
and binary population synthesis (from COMPAS). This allows us to produce
predictions for the rate and properties of an ensemble of GW events for any
model without the large computational cost of redoing the cross-matching from
scratch. From Rauf et al. (2023), the number of GW events per Shark galaxy
that occur at time $t_{j}$ from stars born at time $t_{i}$ using COMPAS model
$\alpha$ was written as

$\displaystyle\frac{dN_{\rm GW}^{\alpha}}{dN_{\rm gal}}(Z,t_{i},t_{j})$
$\displaystyle=\frac{\psi(Z,t_{i})}{M_{\rm tot}^{\alpha}}\Delta
t_{i}\sum_{n}^{N_{\mathrm{binary}}}\vartheta(t_{i}+\tau_{n}^{\alpha}(Z),t_{j},\Delta
t_{j}),$ (12)

where

$\displaystyle\Delta t_{j}$ $\displaystyle=t_{j}-t_{j-1},$ (13)
$\displaystyle\vartheta(t,t_{j},\Delta t_{j})$
$\displaystyle=\left\\{\begin{array}[]{cc}1&\hskip 4.2679pt\text{if
$t_{j}-\frac{\Delta t_{j}}{2}<t<t_{j}+\frac{\Delta t_{j}}{2}$,}\\\ 0&\hskip
4.2679pt\text{otherwise.}\\\ \end{array}\right\\}$ (16)

and $\psi(Z,t_{i})$ is the metallicity ($Z$) dependent star formation rate in
the galaxy when the progenitor stars were born. $M_{\rm tot}^{\alpha}$ is a
normalisation constant for a given model $\alpha$ (see Eq 10 in Rauf et al.
2023), which normalises the star formation rate based on the chosen stellar
initial mass function, in such a way that it cancels out the dependence of the
number of GW events on the number of binaries that we simulate
$N_{\mathrm{binary}}$, and any lower or upper mass cuts we apply when running
COMPAS. $\vartheta(t,t_{j},\Delta t_{j})$ is a counter for the number of BBH
mergers in Shark snapshot, $j$. $\tau_{n}^{\alpha}(Z)$ is the coalescence time
of each binary outputted from the model $\alpha$. If the merger time,
$t=t_{i}+\tau_{n}^{\alpha}(Z)$, falls in the $j^{\mathrm{th}}$ bin, a count is
added to that bin, otherwise it does not contribute. Overall, this equation
encapsulates that the rate consists of two contributions, one from the Shark
host galaxy and the other from population synthesis (in our case from COMPAS).

###  3.1 The effective factor for different models

Cross-matching the star formation history of every galaxy in a simulation with
a COMPAS run is computationally expensive. To simplify this, we now
demonstrate how we can rescale the rate. If instead of considering GW events
on a per galaxy basis we consider them in terms of number per metallicity and
time $t_{i}$ and $t_{j}$ bins, we can write

$\displaystyle\frac{d^{3}N_{\rm GW}^{\alpha}}{dZdt_{i}dt_{j}}$
$\displaystyle=\frac{d^{3}N_{\rm gal}}{dZdt_{i}dt_{j}}\frac{dN_{\rm
GW}^{\alpha}}{dN_{\rm gal}}$ $\displaystyle=\frac{d^{3}N_{\rm
gal}}{dZdt_{i}dt_{j}}\frac{\overline{\psi}(Z,t_{i})}{M_{\rm
tot}^{\alpha}}\Delta
t_{i}\overline{\sum_{n}\theta(t_{i}+\tau_{n}^{\alpha}(Z,t_{i}),t_{j},\Delta
t_{j})}.$ (17)

In the second line we have substituted in Eq. 12 and use $\overline{\psi}$,
$\overline{\theta}$ to acknowledge that this is the average of the star-
formation history and Heaviside function over all galaxies that fall into the
the bin $dZdt_{i}dt_{j}$ (i.e., we assume that all galaxies in this bin
produce the same number of gravitational waves, which is strictly true only
for infinitesimally small bins). We now can introduce a new population
synthesis model, $\beta$, where the merger rate is unknown. If we have the
merger rate from model $\alpha$, we can take a ratio of the two models and
rearrange,

$\displaystyle\frac{d^{3}N_{\rm GW,k}^{\beta}}{dZdt_{i}dt_{j}}$
$\displaystyle=\frac{M_{\rm tot}^{\alpha}}{M_{\rm
tot}^{\beta}}\frac{d^{3}N_{\rm
GW}^{\alpha}}{dZdt_{i}dt_{j}}\frac{\overline{\sum_{n}\theta(t_{i}+\tau_{n}^{\beta}(Z,t_{i}),t_{j},\Delta
t_{j})}}{\overline{\sum_{n}\theta(t_{i}+\tau_{n}^{\alpha}(Z,t_{i}),t_{j},\Delta
t_{j})}}$ $\displaystyle=\frac{M_{\rm tot}^{\alpha}}{M_{\rm
tot}^{\beta}}\frac{d^{3}N_{\rm
GW}^{\alpha}}{dZdt_{i}dt_{j}}f^{\alpha\beta}_{\rm eff}(Z,t_{i},t_{j},\Delta
t_{j}).$ (18)

We hence find that retaining the information from Shark is not required as the
star formation history, $\overline{\psi}(Z,t_{i})$, and the galaxy
distribution are independent of the population synthesis model. In practice we
also expect that for most use cases $\frac{M_{\rm tot}^{\alpha}}{M_{\rm
tot}^{\beta}}\approx 1$ since changes to the initial mass function between
different literature values primarily affect stars that are too small to
become binary black holes or neutron stars, although we retain it here for
completeness.111Note this factor also needs to be included if one were to
rescale between models where different lower or upper mass limits have been
used to draw stars from the Initial Mass Function (IMF), where a different
mass ratio distribution for the initial binaries is assumed, or the overall
binary fraction is changed The only caveat then is on the free choice of
binning used to evaluate what we call the ‘effective factor’
$f^{\alpha\beta}_{\mathrm{eff}}$ between the two models — although the
dependence on galaxy simulation has been absorbed into our fiducial
$\frac{d^{3}N_{\rm GW}^{\alpha}}{dZdt_{i}dt_{j}}$, there is still the implicit
assumption that the rate in each galaxy in the bin is the same. We find that
this approximation remains valid over a range of metallicity binning choices,
even for the quite wide bins we adopt as our default (100 un-evenly spaced
bins for $\log Z\in[-7.0,-0.3]$). We show an example distribution of $f_{\rm
eff}$ in these metallicity and merger time bins in Fig. 1 for a single choice
of binary formation redshift of 2.5 (corresponding to $t_{i}\approx 2.62\rm
Gyr$).

Figure 1: An example of the effective factor between two COMPAS models with
different remnant mass prescriptions $\alpha=$ Fryer et al. (2012) and
$\beta=$ Mandel & Müller (2020) in bins of merger time and metallicity for
binaries forming at redshift 2.5 (corresponding to $t_{i}\approx 2.62\rm
Gyr$). In this plot and the rest of our work, an un-equal metallicity binning
is chosen such that the number of Shark galaxies in each metallicity bin in
the same. Unshaded bins indicate regions of the parameter space where there
are no gravitational waves forming in either one or both of the models.

###  3.2 Generating the merger rate

Here we give a qualitative overview of our algorithm for predicting the GW
merger rate for our different COMPAS models.

  1. 1.

The starting point of our algorithm is a high-fidelity set of GW events
simulated from a combination of Shark and COMPAS model $\alpha$. The detailed
steps used for generating these are covered in Rauf et al. (2023) and so will
not be repeated here. Briefly, this is done by cross-matching the galaxy
metallicities to the binary population metallicities outputted from COMPAS. We
inject this population into the Shark galaxies and track their evolution until
they merge using Eq 16.

  2. 2.

From this, $\frac{d^{3}N_{\rm GW}^{\alpha}}{dZdt_{i}dt_{j}}$ is calculated by
taking the fiducial merger rate outputs and binning in terms of metallicity
and birth and merger time. For the metallicity we use quantiles to determine
the bin edges such that there is an even number of galaxies in each bin. For
the birth and merger time, we simply use the 170 snapshots from Shark to
define the 170 bins.

  3. 3.

Having established our metallicity, birth time and merger time bins, we take
our full list COMPAS binaries, read in the metallicity and evolution time,
which is added to the birth time ($t_{\rm birth}+t_{\rm BBH}+t_{\rm
coalescence}$), and add these to the corresponding bins. We apply this to both
the fiducial and new model, take the ratio and multiply by the merger rate in
(ii) to estimate the new merger rate, shown in Eq 18.

Figure 2: Volumetric BBH merger rate with remnant mass prescription from
Mandel & Müller (2020) as a function of redshift. The blue line is the true
merger rate from running COMPAS and Shark. The red dashed line is the merger
rate rescaled with the effective factor calculated using the full COMPAS run
containing $170\times 10^{6}$ binaries. The grey lines are the rescaled merger
rate where the effective factor has been calculated using a unique sub-sample
of $10^{6}$ binaries from the full COMPAS run. The green line is the average
over 170 of these subsamples.

Once we have our estimate for the number of gravitational waves in our new
model $\beta$ we evaluate the volumetric merger rate at time $t_{j}$ by simply
summing over the metallicity bins and birth time bins,

$\displaystyle R_{\rm GW}^{\beta}(t_{\rm
j})=\frac{1}{V_{\mathrm{sim}}}\sum_{i=0}^{170}\sum_{k=0}^{100}\frac{d^{3}N_{\rm
GW}^{\beta}}{dZ_{k}dt_{i}dt_{j}},$ (19)

where $V_{\mathrm{sim}}$ is the volume of our original Shark simulation.

For efficiency, we choose a sample $\tilde{N}_{\beta}^{\rm GW}$ of $10^{6}$
binaries from our $\beta$ COMPAS run to estimate the effective factor, while
our fiducial model $\alpha$ has $170\times 10^{6}$.222Note that we draw only
$\sim 10^{6}$ unique masses for this process from the IMF, but each one is
randomly assigned a different metallicity from Shark, resulting in $\sim
170\times 10^{6}$ unique combinations of initial masses and metallicities, and
so $\sim 170\times 10^{6}$ unique binaries. The merger rate is reduced by a
factor of 170 but this is compensated for by the normalisation $\frac{M_{\rm
tot}^{\alpha}}{M_{\rm tot}^{\beta}}\sim 170$. Thus, this formalism removes the
dependence on the sample size. We test the validity of our rescaling procedure
and using only $10^{6}$ binaries for model $\beta$ by splitting one of our
fiducial COMPAS runs from Rauf et al. (2023) (using the remnant mass model of
Mandel & Müller 2020) into 170 by $10^{6}$ sub-samples. We then use these to
estimate the rescaled merger rate for this model given $\alpha=\,$Fryer et al.
(2012), and look at the spread in the merger rate. We show in Fig. 2 that
overall the rescaling method using equal numbers of binaries for models
$\alpha$ and $\beta$ recovers the true merger rate from our high-fidelity
cross-matching, but also that using only $10^{6}$ binaries for model $\beta$
returns, on average, an unbiased estimate of the merger rate. The use of a
smaller number of binaries in the rescaling procedure does introduce
additional Poisson scatter in the merger rate, which we propagate through into
our model comparison/evaluation (see Section 4.

In Fig. 3, we show the volumetric BBH merger rates in the Shark simulations
using the effective factor approximation for four different models. In all
cases we find that we are able to recover rates from Rauf et al. (2023) very
well. We note that the plot is noisier than in Rauf et al. (2023) due to the
smaller sample of BBHs taken from each COMPAS run, but again this noise arises
only from Poisson sampling error.

Figure 3: Volumetric rate plot as a function of redshift from Rauf et al.
(2023) recreated with the effective factor scaling method. The error on the
plots are the range of merger rates from sampling from each COMPAS run 170
times.

##  4 Bayesian inference and model comparison

In this section we outline our procedure for comparing different COMPAS models
to data from GWTC-1 (Abbott et al., 2019a), GWTC-2 (Abbott et al., 2021a) and
GWTC-3 (Abbott et al., 2023b) including our rescaling methodology. We validate
this by method by comparing the likelihoods given the data for the four models
from Rauf et al. (2023) for which we have both ‘brute-force’ measurements of
the GW rate and our rescaled calculations.

###  4.1 Selection criteria

The Bayesian analysis will depend strongly on $N_{\rm obs}$, which is the
number of confident GW detections. How confident we are in the GW signals and
the validity of the sources is determined by the following ranking statistics:

  1. 1.

Signal-to-noise ratio (SNR): This is proportional to the signal amplitude
divided by the noise amplitude. The higher the SNR, the stronger the signal
compared to the noise, and more likely that the signal will be detected. This
is the main ranking statistic when the signal shape is known from theory, such
as the signal for a BBH merger. We use the network SNR, which is the square
root of the sum of the squares of SNRs from the triggers which form the GW
event (only from the detectors that produce triggers contributing to the
event; Cabourn Davies & Harry 2022). The SNR, $\rho$, in a single detector is
defined as,

$\displaystyle\rho^{2}=\langle h|h\rangle=4\mathcal{R}\int_{f_{\rm
min}}^{f_{\rm max}}\frac{h^{*}(f)h(f)}{S(f)}df,$ (20)

where $h(f)$ is the waveform measured by the detector, $S(f)$ is the power
spectrum density and $f_{\rm min}$ and $f_{\rm max}$ are the minimum and
maximum limits for the considered frequency range (Farr et al., 2015; Barrett
et al., 2018).

  2. 2.

False alarm rate (FAR): Once the threshold for the SNR is set, the FAR can be
determined, providing an alternative way to estimate the significance of the
event (Zheng et al., 2021). The FAR is the number of noise events with a
signal equal to or greater than the GW event, divided by the total duration of
the background data. The false alarm rate of single-detector candidates given
$\rho$ is,

$\displaystyle FAR(\rho)=\int_{\rho}^{\rho_{\rm
max}}\Lambda_{n}p_{n}(\rho^{\prime})d\rho^{\prime},$ (21)

where $\Lambda_{n}$ is the mean Poisson rate of noise triggers and
$p_{n}(\rho^{\prime})$ is the probability densities describing the
distribution of $\rho$ (Callister et al., 2017).

  3. 3.

$p_{\rm astro}$: The probability that a particular event is of astrophysical
origin, as opposed to noise fluctuations or terrestrial, is defined by $p_{\rm
astro}$. It is measured by comparing the GW event rate to the background event
rate for a fixed SNR and FAR. The probability of a trigger, $x$, being
astrophysical is,

$\displaystyle p_{\rm
astro}(x)=\frac{p(x|\mathcal{S})\pi(\mathcal{S})}{p(x|\mathcal{S})\pi(\mathcal{S})+p(x|\mathcal{N})\pi(\mathcal{N})},$
(22)

where $\mathcal{S}$ and $\mathcal{N}$ describe the signal and noise model, as
well as the distribution of signal and noise sources respectively (Callister
et al., 2017; Broekgaarden et al., 2023).

For our most rigorous selection, we combine $\rm FAR\leq 0.25\,yr^{-1}$ and
$p_{\rm astro}\geq 0.9$, as the FAR alone will not provide information on the
sources being astrophysical, but rather the behaviour of the signal
background. The $p_{\rm astro}$ reveals more events in the population rich
areas of the parameter space, and is a more realistic metric for setting
threshold for electromagnetic follow-up (Lynch et al., 2018; Andres et al.,
2022). The caveat of $p_{\rm astro}$ is that that independent pipelines will
obtain different $p_{\rm astro}$ due to a difference in data reduction
techniques. Banagiri et al. (2023) proposes an unified $p_{\rm astro}$ by
combining triggers from all pipelines, which incorporates extra information
about a signal. Galaudage et al. (2020) also emphasises the need for a
pipeline independent $p_{\rm astro}$, as they show that the posterior
distribution of the population hyper-parameters can change with an
astrophysically motivated prior. Given these caveats, we also investigate a
simpler, broader selection using only an $SNR\geq 8$ cut, which matches that
used in evaluating our detection probability function $p_{\rm det}(\theta)$
for our models and provides us with a larger samples of events. However, the
downside of this is that it potentially allows events into our sample that may
not be astrophysical. In these cases, we allow our Bayesian framework to
decide if these events should be ruled out.

Our supersample of GW events for querying the PDF for each COMPAS model
consists of confirmed BBHs and BH-mass-gap events such as GW190521 (Abbott et
al., 2020a) and GW190814 (Abbott et al., 2020c). This gives a total of 93
catalog GW detections. Of these, the following GW events are not included in
our analysis for the following reasons:

  * •

GW170817: Confirmed BNS (Abbott et al., 2017)

  * •

GW190425: Confirmed BNS (Abbott et al., 2020b)

  * •

GW190426_152155: Not astrophysical (Abbott et al., 2021a)

  * •

GW190909_114149: Not astrophysical (Abbott et al., 2021a)

  * •

GW191219_163120: Confirmed NSBH (Abbott et al., 2023b)

  * •

GW200115_042309: Confirmed NSBH (Abbott et al., 2021b)

We then apply our two sets of selection criteria. We have 79 GW events with
$\rm SNR\geq 8$333We use the network SNR provided in
https://gwosc.org/eventapi/html/GWTC/?page=1&pagesize=100. and 60 GW events
with $p_{\rm astro}\geq 0.9$ and $\rm FAR\leq 0.25yr^{-1}$. Tab 3 lists the
events used with each selection criteria. Out of the 87 BBH GW events we could
have included (93 less the 6 listed above) the following GW events do not fit
either criteria:

  * •

GW190403_051519

  * •

GW190514_065416

  * •

GW191113_071753

  * •

GW200208_222617

  * •

GW200220_061928

  * •

GW200306_093714

  * •

GW200308_173609

  * •

GW200322_091133

###  4.2 Fitting LVK posterior for various COMPAS models

Our basic fitting methodology is based on the widely-used Bayesian framework
(e.g., Stevenson et al., 2015; Bouffanais et al., 2019; Mandel et al., 2019;
Vitale et al., 2022; Périgois et al., 2023). Given the BBH GW data from some
observing runs, $\mathcal{H}$, the posterior for a model specified via hyper-
parameters $\lambda$ (for instance a certain mass transfer efficiency $f_{\rm
MT}$ or specific angular momentum loss $f_{\gamma}$), takes the form of an
inhomogeneous Poisson distribution,

$\mathcal{L}(\lambda,N_{\lambda}|\mathcal{H})=\pi(\lambda)e^{-\mu_{\lambda}}N_{\lambda}^{N_{\rm
obs}}\prod_{k}^{N_{\rm obs}}\mathcal{I}^{k}.$ (23)

$N_{\lambda}$ is the number of GW events predicted by the model in a given
time-frame, while $\mu_{\lambda}$ is the number the model suggests we would
actually detect. For convenience we can use these to define some detection
efficiency $\beta(\lambda)=\frac{\mu_{\lambda}}{N_{\lambda}}$. $N_{\rm obs}$
is the actual number of detected events in our dataset. We assume that a
priori all hyper-parameter combinations are equally likely, so
$\pi(\lambda)=1$. Finally, $\mathcal{I}^{k}$ denotes the likelihood of
observing GW event $h^{k}$ given the hyper-parameters $\lambda$, written as a
marginalisation over the binary properties predicted by that model,

$\displaystyle\mathcal{I}^{k}=\int\mathcal{L}(h^{k}|\theta)p_{\rm
sim}(\theta|\lambda)d\theta=\frac{1}{N_{s}^{k}}\sum_{i=1}^{N_{s}^{k}}\frac{p_{\rm
sim}(\theta_{i}^{k}|\lambda)}{\pi^{k}(\theta_{i}^{k})},$ (24)

where $\theta=\\{\mathcal{M}_{c},q,z\\}$ are the GW event parameters chirp
mass, mass ratio and redshift respectively and $p_{\rm sim}(\theta|\lambda)$
is the probability distribution function for the event-level parameters
generated by COMPAS using the Shark metallicity as input. The latter sum
arises from the use of the $N_{s}^{k}$ publicly available posterior samples to
represent the likelihood of event $k$, $\mathcal{L}(h^{k}|\theta)$. Formally
$p_{\rm sim}(\theta|\lambda)$ is obtained by marginalising over all the
metallicities and birth times used that can give rise to an event with
parameters $\theta$:

$p_{\rm sim}(\theta|\lambda)=\int_{0}^{\infty}\int_{0}^{t_{j}}p_{\rm
sim}(\mathcal{M}_{c},q,z(t_{j})|Z,t_{i},\lambda)\,dt_{i}dZ.$ (25)

$N_{\lambda}$ is the number of GW events in the Universe within an observing
time that is expected from each COMPAS model. We can evaluate this given our
rescaled merger rate density (expressed in individual metallicity and birth
time bins, and in units of $\rm Mpc^{-3}yr^{-1}$) and $p_{\rm sim}$,

$\displaystyle N_{\lambda}$ $\displaystyle=4\pi T_{\rm obs}\int_{\theta\in
V_{\theta}}\int_{0}^{\infty}\int_{0}^{t_{j}(z)}p_{\rm
sim}(\theta|Z,t_{i},\lambda)$ $\displaystyle\qquad\qquad R_{\rm
GW}^{\lambda}(Z,t_{i},t_{j}(z))\,\chi^{2}(z)\frac{c}{H(z)}\frac{1}{1+z}dt_{i}dZd\theta,$
(26)

where $\chi(z)$ is the comoving distance to the redshift $z$, $c$ is the speed
of light, $H(z)$ is the Hubble parameter, and $T_{\rm obs}$ is the total LVK
observing time (duty cycle). This is 118.0, 138.8 and 142.0 days for O1, O2
and O3 respectively, which we converted to years.

We can extend this to the predicted detected number of GW events using a
detection probability for the event-level parameters $p_{\rm
det}(\mathcal{M}_{c},q,z)$ given a particular detector configuration and
signal-to-noise threshold.

$\displaystyle\mu_{\lambda}$ $\displaystyle=4\pi T_{\rm obs}\int_{\theta\in
V_{\theta}}\int_{0}^{\infty}\int_{0}^{t_{j}(z)}p_{\rm
sim}(\theta|Z,t_{i},\lambda)p_{\rm det}(\theta)$ $\displaystyle\qquad\qquad
R_{\rm
GW}^{\lambda}(Z,t_{i},t_{j}(z))\,\chi^{2}(z)\frac{c}{H(z)}\frac{1}{1+z}dt_{i}dZd\theta.$
(27)

For Eq. 23, we are assuming that $N_{\lambda}$ events are all formed via the
isolated formation channel. There is some evidence for some events in GWTC-3
(Abbott et al., 2023b) forming via the dynamical channel (Arca Sedda, 2021;
Romero-Shaw et al., 2022). To avoid any assumptions on the formation channel
and discount information on the predicted vs. observed GW rate, we can
marginalise over $N_{\lambda}$ to remove its dependence by multiplying with a
prior $\pi(N_{\lambda})=\frac{1}{N_{\lambda}}$ and integrating (Fishbach et
al., 2018).

$\displaystyle\mathcal{L}(\lambda|\mathcal{H})$
$\displaystyle=\int\frac{1}{N_{\lambda}}e^{-N_{\lambda}\beta(\lambda)}N_{\lambda}^{N_{\rm
obs}}\prod_{k}^{N_{\rm obs}}\mathcal{I}^{k}dN_{\lambda}$
$\displaystyle=(N_{\rm obs}-1)!\prod_{k}^{N_{\rm
obs}}\frac{\mathcal{I}^{k}}{\beta(\lambda)},$ (28)

We make model inferences using both of these posteriors in our work and
compare how our conclusions change when $N_{\lambda}$ is included as
information.

In order to evaluate our posteriors, we first calculate the model volumetric
rate using a variant of Eq. 19 but kept in separate metallicity and birth time
bins. Given $T_{\rm obs}$, the merger rate and the appropriate volumetric
factors, we then find the number of GW events that would fall in each
metallicity and time bin, $N_{\lambda,i,j}$, and draw that many events
(characterised by a chirp mass and mass ratio) from our full list of $10^{6}$
COMPAS binaries.444We find $N_{\rm COMPAS,i,j}$ in each metallicity and
snapshot bins from the COMPAS run. We calculate a weight defined as
$N_{\lambda,i,j}/N_{\rm COMPAS,i,j}$ to avoid repeatedly drawing from the
COMPAS run. We apply this weight to the PDF. When calculating the detected
PDF, we multiply this weight by $p_{\rm det}$. We store the chirp mass, mass
ratio and merger redshift of the events across all metallicity and birth time
bins before then building the PDF $p_{\rm\textsc{sim}}(\theta|\lambda)$ using
Kernel Density Estimation, in particular using the fast KDEpy.FFTKDE package.
The bandwidth for our KDE is calculated using Silverman’s rule of thumb
(Silverman, 1986) for each hyper-parameter distribution and averaged to obtain
the single band-width for the full 3D distribution. As a by-product of method
we can also estimate $N_{\lambda}$ and $\mu_{\lambda}$ at the same time by
simply keeping track of how many COMPAS binaries we are drawing.

For $p(\theta_{i}^{k}|\lambda)$, we query the KDE-built PDF using the chirp
mass, mass ratio and redshift from the public LVK posterior samples. We
construct $\pi^{k}(\theta_{i}^{k})$ by building another KDE with the public
LVK prior samples and querying this again with the posterior samples at each
point $\theta_{i}^{k}$. In the case where $p_{\rm sim}(\theta|\lambda)=0$ for
all posterior samples, this implies that the COMPAS model cannot reproduce the
BBHs in the posterior samples, as they fall outside of the KDE. We flag these
events and ultimately remove them from the analysis for that specific model
comparison. For each of the $N_{\lambda,i,j}$, we calculate the detection
probability, $p_{\rm det}$, based on the masses, redshift and luminosity
distance. The detection probability comes from a selection effect function in
COMPAS, which takes the BH masses, redshift and luminosity distance as inputs
(Barrett et al., 2018; Team COMPAS: Riley, J. et al., 2022). The SNR is
calculated by computing the source waveform using the waveform model from the
LAL software suite (LIGO Scientific Collaboration et al., 2018).

We apply the same procedure also to estimate marginalised, observed 1D
distributions $p_{\rm\textsc{sim\\_det}}(\mathcal{M}_{c}|\lambda)$,
$p_{\rm\textsc{sim\\_det}}(q|\lambda)$ and
$p_{\rm\textsc{sim\\_det}}(z|\lambda).$555$p_{\rm\textsc{sim\\_det}}$ refers
to the detected PDF, which is marginalised over the other two fitting
parameters after the detection probability is applied. We show these 1D
distributions for our four fiducial models in Fig. 4 alongside the GWTC data,
with our main observations being:

  * •

The models follow the chirp mass distribution reasonably well except for the
chirp mass range 100–130 M⊙. These models predict high mass BHs that either
may have not been detected due to the rarity of finding these types of events
with current detectors, as they may be associated with noise in the low
frequency band, or as a result of inaccuracies in the models resulting in an
incorrect representation of the observed population.

  * •

The mass ratio distribution is quite noisy but overall relatively flat,
similar to the observed distribution. However for Mandel & Müller (2020) and
Schneider et al. (2021) there is a peak in the distribution after $q\approx
0.8$ indicating a preference for equal mass binaries. Mandel & Müller (2020)
with $f_{\rm WR}=0.2$ has a preference for $q\sim 0.6$ and Fryer et al. (2012)
for $q\sim 0.7-0.8$.

  * •

The shape of the redshift distributions for the models fit the observed
distribution quite well. Overall, the Fryer et al. (2012) model is the best
fit, while the other models have a wider distribution. This indicates BBHs
from these models are more likely to be detected at higher redshifts.

Figure 4: Distribution of $\mathcal{M}_{c}$, $q$ and $z$ for different models.
Difference colours refer to the different remnant mass prescriptions
investigated in Rauf et al. (2023). The grey shaded regions are the inferred
distribution from the LVK collaborations, calculated using the posterior
samples for the $N_{\rm obs}=79$ GW events based on the SNR $\geq$ 8 cut.

###  4.3 Quantifying the uncertainty in $N_{\lambda}$,
$p(\lambda,N_{\lambda}|\mathcal{H})$ and $p(\lambda|\mathcal{H})$

For our four models in Rauf et al. (2023) we were able to split COMPAS outputs
into 170 sub-samples to compute 170 effective factors and produce 170
theoretical, $N_{\lambda}$, and detected, $\mu_{\lambda}$, numbers of GW
events. We show the distribution of the detected number of events in Fig. 5.
For the error on the total number of events, we can use the standard deviation
across our 170 subsamples.

Figure 5: Normalised histograms of the detected number of BBH merger events
calculated using Eq. 27 over the LVK observing time for the last three
observing runs. The coloured distributions are obtained from the predictions
from the model for 170 sub-samples. The grey band represents the uncertainty
in the observed number of BBH mergers depending on the selection criteria
adopted (see Section 4.1). The dotted lines are the number of events predicted
using the true merger rate from Rauf et al. (2023).

However, if we have only one sample of $N_{\rm sample}$ binaries for a COMPAS,
we require an alternative approach to assign an error for these
models.666Given the number of models we test in this work, the alternatives of
simply using more samples or bootstrapping using the $N_{\rm sample}=10^{6}$
samples would be computationally challenging. Furthermore, subsampling or
bootstrapping to generate an error would fail to capture the impact of
progenitors of rare events that are not represented in the $N_{\rm
sample}=10^{6}$, but which are present in the distributions obtained from our
four models with $N_{\rm sample}=170\times 10^{6}$, and from which our errors
are extrapolated. We find that the ratio of the standard deviation to the mean
number of events as a function of $\log N_{\lambda}$ for the four models has a
linear trend such that the relative error increases as the total number of GW
events generated by the model decreases, in a roughly Poissonion fashion. In
Fig. 6 we show the linear fits for the standard deviation in $N_{\lambda}$ and
$\mu_{\lambda}$, which we show in the following equations:

$\displaystyle\frac{\Delta N_{\lambda}}{N_{\lambda}}$
$\displaystyle\approx-0.04\log N_{\lambda}+0.55$ (29)
$\displaystyle\frac{\Delta\mu_{\lambda}}{\mu_{\lambda}}$
$\displaystyle\approx-0.03\log N_{\lambda}+0.43$ (30)

Figure 6: Relative error on the numbers of generated and detected events
versus the average $\log N_{\lambda}$ for the four Rauf et al. (2023) models
from 170 sub-samples. The crosses correspond to the relative error from the
models and the dashed lines are linear fits to the points. The black points
and fit are for $N_{\lambda}$ and the red points and fits are for
$\mu_{\lambda}$.

We can now apply this to single samples, assuming similar Gaussian
distributions for other COMPAS models. To justify this, we find that the mean
$\mu_{\lambda}$ over the 170 sub-samples in Fig. 5 is equivalent to the
$\mu_{\lambda}$ we obtain using the true merger rate from the full COMPAS
output (one sample) for all four models that we test.

Now to address the uncertainty in $\mathcal{I}^{k}$; We find another linear
trend for the standard deviation in $\sum_{k}^{N_{\rm
obs}}\log\mathcal{I}^{k}$ over the 170 sub-samples as a function of the mean
$\log N_{\lambda}$ for both our selection criteria, as shown in Fig. 7.

Figure 7: The standard deviation in the sum of the log-likelihood per event,
$\log\mathcal{I}/N_{\rm obs}$, versus the mean $\log N_{\lambda}$ over the 170
sub-samples for the four Rauf et al. (2023) models. We show the linear fit
(green dashed line) that applies to both sets of selection criteria; $\rm
FAR\leq 0.25yr^{-1}$ and $p_{\rm astro}\geq 0.9$ (black) and $\rm SNR\geq 8$
(red).

The linear fits for the standard deviation in the likelihood sums are given
by:

$\displaystyle\frac{1}{N_{\rm obs}}\Delta\sum_{k}^{N_{\rm
obs}}\log\mathcal{I}^{k}\approx-0.05\log N_{\lambda}+0.70$ (31)

For new COMPAS models shown in Section 5, we are able to use Eq. 31 and to
calculate the error on $\sum_{k}^{N_{\rm obs}}\log\mathcal{I}^{k}$ given
$N_{\lambda}$ without needing to resample. We can now propagate the error in
$N_{\lambda}$, $\mu_{\lambda}$ and $\sum_{k}^{N_{\rm obs}}\log\mathcal{I}^{k}$
into the error on the log-likelihood.

$\displaystyle\Delta\log\mathcal{L}(\lambda|\mathcal{H})$
$\displaystyle=\sqrt{\left(N_{\rm
obs}\frac{\Delta\mu_{\lambda}}{\mu_{\lambda}}\right)^{2}+\left(N_{\rm
obs}\frac{\Delta
N_{\lambda}}{N_{\lambda}}\right)^{2}+\left(\Delta\sum_{k}^{N_{\rm
obs}}\log\mathcal{I}^{k}\right)^{2}}$
$\displaystyle\Delta\log\mathcal{L}(\lambda,N_{\lambda}|\mathcal{H})$
$\displaystyle=\sqrt{(\Delta\mu_{\lambda})^{2}+\left(N_{\rm obs}\frac{\Delta
N_{\lambda}}{N_{\lambda}}\right)^{2}+\left(\Delta\sum_{k}^{N_{\rm
obs}}\log\mathcal{I}^{k}\right)^{2}}$ (32)

Although only rudimentary given the small number of models available to us to
generate fitting formulae, even adopting this simple uncertainty on the
likelihood allows us to consider the statistical variations in model
evaluation using population synthesis. As such it enables us to draw more
confident conclusions on which hyper-parameters cause major changes in the
likelihood, and whether these changes are physically significant.

###  4.4 Comparing posteriors with and without rescaling merger rate

Figure 8: $\mathcal{I}^{k}$ term for each GW event with $\rm SNR\geq 8$. Each
colour represents a different remnant mass prescription from COMPAS and the
error bar represents the sampling error over 170 samples when building the
simulated PDF of $\mathcal{M}_{c}$, $q$ and $z$ given the model.

Fig. 8 shows the $I^{k}$ term for individual GW events. This is a practical
approach to determine the probability of individual GW events occurring in the
populations generated by the COMPAS models. The Schneider et al. (2021) model
generally returns high probabilities when querying the PDF, while the Fryer et
al. (2012) model is disfavoured. We find some events are outliers and are not
well reproduced by some or all of the models relative to the other events. As
an example, Mandel & Müller (2020) with $f_{\rm WR}=0.2$ and Schneider et al.
(2021) are less likely to output the progenitors of GW190814, GW190917_114630
and GW200210_092254 compared to Fryer et al. (2012) and Mandel & Müller
(2020). For GW190426_190642, GW190602_175927 and GW191109_010717 Schneider et
al. (2021) is distinctly preferred over the other models, although in these
cases, and especially that of GW190521, the models are generally less able to
produce the required event-level properties. This is expected for GW190521 as
both masses, $\sim 85M_{\odot}$ and $\sim 66M_{\odot}$, most likely lie in the
pair-instability supernovae mass gap, making it atypical of binaries formed
via the isolated formation channel (Abbott et al., 2020a).

Figure 9: Normalised histograms of the log likelihoods obtained using 170 sub-
samples for the models explored in Rauf et al. (2023) with and without the
merger rate dependence. The top two panels are the likelihoods comparing the
models to the GW events that fit the $\rm FAR\leq 0.25yr^{-1}$ and $p_{\rm
astro}\geq 0.9$ criteria, while the bottom two panels compares to the $\rm
SNR\geq 8$ criteria. The dotted lines are the likelihoods using the true
merger rates calculated from Rauf et al. (2023). In all cases we see that the
true likelihood is well represented within our distribution of 170 subsamples,
indicating our use of a relatively smaller number of $10^{6}$ binaries is not
expected to bias our results.

For a more comprehensive assessment of the models, we can analyse the
likelihoods. In Fig. 9 we now incorporate the detection probability and show
the likelihoods with and without $N_{\lambda}$ dependence. The Mandel & Müller
(2020) appears to the best fit to the observed population. This can be
justified with Fig. 4. The Mandel & Müller (2020) model has two peaks in the
chirp mass distribution and a relatively flat mass ratio distribution, similar
to the LVK population distribution. If we assume that the models output
accurate populations, then incorporating the number of GW events makes the
likelihood worse – Fryer et al. (2012) greatly overestimates the number of
detected events, while Schneider et al. (2021) underestimates it. In general
the latter is preferred as it leaves ‘room’ for some of the LVK observed
events to arise from the dynamical formation channel. The exception is the
Mandel & Müller (2020) model, which performs well since the model predicts the
same number of GW events within the observed range. The true number of events
predicted is $\sim 60$, comparable to $N_{\rm obs}$ for the $p_{\rm astro}\geq
0.9$ and $\rm FAR\leq 0.25yr^{-1}$ selection criteria. For $\rm SNR\geq 8$, we
see an increase in the likelihood for all models, but the distributions remain
the same. The biggest difference in the unmarginalised likelihood for our two
selections occurs for the Fryer et al. (2012) model, which produces the
largest number of detected events out of the four models. This is due to
$N_{\rm obs}$ increasing to 79, such that the difference between the predicted
$\mu_{\lambda}$ and the number of observed LVK events is slightly reduced.

##  5 Results

Having validated our methodology using the four models analysed in Rauf et al.
(2023), this section discusses the outputs from our new set of COMPAS models
where we vary the range of hyper-parameters described previously in Section 2
(see also Table 1. We start with describing the trends in the number of
detected events when varying the COMPAS hyper-parameters. We calculate the
Bayes factor to identify groups of models that are more preferred. We then
show our likelihood analysis, where we discuss the best and worst models when
compared to the observed population parameters. We follow up this analysis and
look at specific events and their $\mathcal{I}^{k}$ value to determine which
models are more and less likely to generate the progenitors of these events.
See Tables 4, 5 and 6 in Appendix C for a detailed description of the hyper-
parameters for all 304 COMPAS models in this analysis.

###  5.1 $\mu_{\rm GW}$ versus COMPAS parameters

Figure 10: Detected number of BBH merger events versus the mass transfer
fraction for no CHE at fixed $f_{\gamma}$. Left: We fix the common envelope
evolution parameter $\alpha_{\rm CE}=0.1$ and show results for different
specific angular momentum loss parameters. Right: We show the models for any
$f_{\gamma}$ and $\alpha_{\rm CE}$ which output $\mu_{\lambda}$ close to the
number of observed LVK events given by the grey band.

Since we explore the number of GW events as a parameter in our posterior
analysis, in this section we comment briefly on the trends between the
detected number of GW events, $\mu_{\lambda}$, and the COMPAS hyper-
parameters.

In the left plot of Fig. 10, we see $\mu_{\rm GW}$ increase with mass transfer
fraction, $f_{\rm MT}$, when the corresponding angular momentum is lost
directly from the accretor ($f_{\gamma}=0$). Increasing the mass transfer can
lead to more equal mass objects, which we find most of our COMPAS models
producing. This trend disappears for higher $f_{\gamma}$. We credit this to
non-trivial evolutionary effects such as the stability prescriptions, that
determine any future interactions within the binary. This preference for
intermediate values of $f_{\gamma}$ agrees with the findings by Santoliquido
et al. (2021). Irrespective of $f_{\rm MT}$, $\mu_{\rm GW}$ peaks for moderate
values of $f_{\gamma}=0.5$ but drops again for high $f_{\gamma}$ when angular
momentum loss occurs closer to the L2 point. At this point, mass is lost at
higher velocities so more angular momentum is lost from the system, which in
turn can induce CE evolution and merger during the mass transfer stage,
stopping the formation of a BBH. These results are in agreement with Lu et al.
(2023); van Son et al. (2022b) that the number of stable mass channel objects
increases for decreasing $f_{\gamma}$, as opposed to the objects formed by CE
evolution. Marchant et al. (2021) also found that the stable channel dominates
the BBH merger rate. We note that when $f_{\rm MT}=1$, $f_{\gamma}$ should
have no impact. Any variations in the number of events at $f_{\rm MT}=1$ is
hence due to Eddington-limited accretion occurring after the first mass-
transfer event(s), once the accretor has become a BH. 777As a result of this
finding, dual functionality in the AM loss was later implemented in COMPAS, as
discussed in Section 2.4. If specific AM loss is set at the compact object as
well as prior to the formation of the first compact object we expect the
number of events will no longer vary for $f_{\rm MT}=1$, although at time of
writing this functionality is yet to be comprehensively tested.

We show the models that produce the closest $\mu_{\lambda}$ to the observed
LVK events in the right plot of Fig. 10. We find that for $\mu_{\rm GW}$ to be
below 100 events (closer to the number of detected events from all observing
runs), a CE efficiency $\alpha_{\rm CE}=10,100$ and $f_{\gamma}=0.75,1.0$ with
no CHE is required. This implies that the CE is ejected earlier in the
evolution and the orbit tightens at a wider separation, which means longer
merger times and less events in the same observed time. This result validates
the studies by Santoliquido et al. (2021), Bouffanais et al. (2021) and
Grichener (2023). With the exception of these models, the CE efficiency has no
significant impact on the populations, as the range for the number of events
does not vary significantly with $\alpha_{\rm CE}$.

Figure 11: Detected number of BBH merger events versus the Wolf–Rayet
multiplier for fixed CHE modes and the Schneider et al. (2021) remnant mass
prescription. The grey band is the observed number of LVK events.

In Figure 11, $\mu_{\rm GW}$ decreases with increasing $f_{\rm WR}$. We find
this for any CHE mode and remnant mass prescription, but we show this for
Schneider et al. (2021). This is due to increasing the mass loss rates leading
to smaller BHs and wider separations, so less mergers occur in the same
observing time. These predictions agree with the work by Riley et al. (2021)
and Stevenson & Clarke (2022) using COMPAS. We find that for no CHE, the rate
of change in $\mu_{\lambda}$ across $f_{\rm WR}$ is smaller compared to the
pessimistic and optimistic cases.

The pessimistic and optimistic CHE modes produced $>1000$ events due to their
tight orbits. This suggests the CHE channel is not ideal for generating the
current LVK detected population. We verify this with the likelihoods in
Section 5.2.

###  5.2 Bayes Factor and Likelihood Analysis

To identify groups of models that are more preferred than others, we use the
Bayes factor, $\mathcal{B}$. Since the priors for all models are the same, the
Bayes factor is the ratio of the likelihoods between two models, here denoted
$\alpha$ and $\kappa$:

$\mathcal{B(\lambda|\mathcal{H})}=\frac{\mathcal{L}(\lambda_{\alpha}|\mathcal{H})}{\mathcal{L}(\lambda_{\kappa}|\mathcal{H})};\qquad\mathcal{B(\lambda,N_{\lambda}|\mathcal{H})}=\frac{\mathcal{L}(\lambda_{\alpha},N_{\lambda_{\alpha}}|\mathcal{H})}{\mathcal{L}(\lambda_{\kappa},N_{\lambda_{\kappa}}|\mathcal{H})}.$
(33)

In Fig. 12, we show the Bayes factor with the likelihoods using the GW events
with the $\rm SNR\geq 8$ cut.

Figure 12: Comparison of models with $\rm SNR\geq 8$ using the Bayes factor.
Each grid is coloured by $\log_{10}(\mathcal{B})$. For
$\log_{10}(\mathcal{B})>0$, the model in the row (y-axis) is favoured over
than the model in the column (x-axis), while for $\log_{10}(\mathcal{B})<0$
the model in the row is disfavoured over the model in the column. Left plot is
the Bayes factor for the marginalised likelihood and right plot is the Bayes
factor for the unmarginalised likelihood.

We note that the Bayes factor indicates there is strong evidence that some
groups of models are preferred more than others, which we will further discuss
in the likelihood analysis. We emphasise these with highlighted boxes in Fig.
12 labeled A-G. We can understand the range of Bayes factor by considering the
different components that enter Eq 33. The difference in $\log\beta(\lambda)$
between different models ranges from $\sim 0.05-3.44$. However, we find that
the $I^{k}$ terms varies over 50 orders of magnitude (see bottom plot in Fig.
18) and therefore dominates the Bayes factor. This confirms that some of the
models are distinguishable and strongly preferred (more likely to produce
progenitors of particular GW events) over others purely based on the masses
and redshift distributions they predict. When including $N_{\lambda}$ as a
parameter in the likelihood, the $-\mu_{\lambda}+N_{\rm
obs}\log(\mu_{\lambda})$ terms this introduces into the log-likelihood has
significant weight, and comes to dominate the Bayes factor, meaning that a
comparison of two models using this likelihood strongly favours those that
produce a number of events close to the observed number of events in our data
sample, even if the event-level distributions are less well matched. Hence, we
carefully consider the results using both of our likelihoods.

Having established groups of models that are strongly preferred over other
models, we now discuss the likelihoods for individual cases and the trends
within these groups. The following results discussed use the $\rm SNR\geq 8$
selection cut. We see an increase in the likelihood across all models, but the
trends don’t change substantially with the selection criteria and we find the
differences are made more distinct when using the $\rm SNR\geq 8$ sample. We
show this in Tables 4, 5 and 6 with Bayes factors between our various models.

We see an increase in the likelihood from $f_{\gamma}=0$ to $f_{\gamma}=0.5$
and we find that the likelihoods are similar with no CHE and pessimistic and
optimistic CHE. In Fig. 13, we show the likelihoods for fixed $f_{\gamma}$. It
is clear that $f_{\gamma}=0.5$ is preferred over $f_{\gamma}=1.0$. Less BBH
mergers are produced with $f_{\gamma}=1.0$. Models with $f_{\gamma}=0.5$
produce more equal mass binaries, which is more consistent with the LVK
population (Abbott et al., 2019b). However, $f_{\gamma}=0.75$ is unique in
that for low mass transfer $f_{\rm MT}<0.5$ the population is disfavoured
similar to $f_{\gamma}=1.0$, but for high mass transfer $f_{\rm MT}>0.5$ the
model is more favourable. In all cases, we see little impact from the choice
of common envelope efficiency. This highlights the degeneracy between stellar
modelling assumptions and populations. Ultimately, we find that in order to
reproduce the observed events from LVK, we need either momentum loss to occur
close to the accretor, or, if occurring further away, for the mass transfer
fraction to be high. These trends are also highlighted by boxes A and B in
Fig. 12, where moderate values of $f_{\gamma}$ are strongly favoured when
comparing population parameters.

Fig. 14 then explores different combinations of the remnant mass prescription,
CHE mode and Wolf–Rayet multiplier. Overall, we see weaker trends here than
with the mass transfer fraction. However, the remnant mass prescription can be
seen to impact the likelihood, with the Fryer et al. (2012) or Mandel & Müller
(2020) prescriptions preferred. It is worth reiterating that, of these two,
the Fryer et al. (2012) was seen to greatly over-predict the number of
observed events to date. We find a weak trend of increasing likelihood in the
absence of CHE, most evident for the Schneider et al. (2021) model, but
overall less significant an impact than was seen for changes in the mass
transfer and AM loss parameters. We also see a weak preference for larger
$f_{\rm WR}$ depending on the remnant mass model.

Figure 13: Likelihoods for models with no CHE. Each band corresponds to a
fixed $f_{\gamma}$. Within each band, excluding the Rauf et al. (2023) models
in blue, the mass transfer fraction and common envelope efficiency are
increasing with increasing model number, within the ranges in Table 1, with
$f_{\rm MT}$ varying more rapidly. Figure 14: Marginalised likelihood for
varying CHE and remnant mass prescription. In each band, the Wolf–Rayet
multiplier is increased with increasing model number, across the values given
in Table 1.

####  5.2.1 Likelihood including $N_{\lambda}$

When including the dependence on the number of GW events, we are able to
identify more distinct groups of trends, particularly with the CHE modes as
shown in Fig. 15. This is in contrast to the likelihoods without the number of
events, indicating that the CHE mode more strongly affects just the number of
BBH mergers, rather than the underlying distributions of mass or redshift. The
likelihood generally reduces when switching from no CHE to Pessimistic or
Optimistic, which is not as notable for the likelihood without $N_{\lambda}$
dependence. We also see this impact in Fig. 12, where box C and D correspond
to no CHE being strongly preferred over Pessimistic and Optimistic CHE and box
E implies Pessimistic is strongly preferred over Optimistic CHE. Allowing for
CHE would mean tighter orbits and more mergers within Hubble time, which would
over-predict the number of observed events when compared to $N_{\rm obs}$. In
addition, the likelihood now decreases from $f_{\gamma}=0.25$ to
$f_{\gamma}=0.5$. This is also due to the increased number of GW events,
shifting them further from the observed number of LVK events. It then
increases again with further increasing $f_{\gamma}=0.75,1.0$ due to the
reduced numbers of events. Within each choice of CHE, we see a preference for
low mass transfer efficiency for $f_{\gamma}=0$. This is required to keep the
predicted number of events low or comparable to the actual number of LVK
events, as these produce unequal mass wide binaries with long merger times.

Although no CHE is preferred, this can depend on $f_{\gamma}$. We show a
subset of no CHE with varying $f_{\gamma}$ in Fig. 15. Again we note the
degeneracy in the hyper-parameters and populations. Low values of $f_{\rm MT}$
and $f_{\gamma}=0.25$ produce similar populations to $f_{\gamma}=0$, while
models with higher values of $f_{\rm MT}$ at $f_{\gamma}=0.25$ appear similar
to Pessimistic CHE with $f_{\gamma}=0$ populations. Conversely, further
increases to $f_{\gamma}=0.5$ and high $f_{\rm MT}$ produce high-likelihood
populations similar to $f_{\gamma}=0$, with now a decreasing likelihood as
$f_{\rm MT}$ decreases. The likelihood further increases with
$f_{\gamma}=0.75$ particularly again for high $f_{\rm MT}$. There is hence a
flip in the behaviour of increasing mass transfer fraction depending on
whether the momentum loss is occurring closer to the accretor or L2 point.
These changes in likelihood closely reflect the trends seen in the number of
predicted events seen in Fig. 10, but further understanding of this complex
behaviour warrants further study.

Figure 15: Model versus unmarginalised likelihood, where each coloured band
corresponds to the a fixed CHE mode and $f_{\gamma}$. Within each band,
excluding the Rauf et al. (2023) models, the mass transfer fraction and common
envelope efficiency are increasing with increasing model number, within the
ranges in Table 1, with $f_{\rm MT}$ varying more rapidly.

Finally, in Fig. 16, we find a strong case for likelihood increasing with
$f_{\rm WR}$, with the exception of Schneider et al. (2021) with no CHE. This
is also shown in Fig. 12 by boxes F and G, where $f_{\gamma}=0$ and thermal
timescale mass transfer rates are strongly preferred over moderate values of
$f_{\gamma}$ and fixed $f_{\rm MT}$. For fixed CHE and remnant mass
prescription, there is a preference for models with $f_{\rm WR}=10$. This is
because mass loss via strong winds can widen binaries. This would increase the
merger time and decrease the number of GW events closer to the LVK detected
number of events. Overall, again no CHE is preferred over pessimistic and
optimistic CHE — in the presence of CHE, one needs a rather finely-tuned set
of: high Wolf–Rayet factor, high $f_{\gamma}$ and low mass transfer efficiency
to recover a reasonable number of mergers.

Figure 16: Unmarginalised likelihood for No CHE (Left) and pessimistic CHE
(Right) for fixed remnant mass prescriptions. Within each band we vary the
Wolf–Rayet multiplier with values given in Table 1, increasing the value with
increasing model number.

Given the main assumption in this paper being all GW events originate from the
isolated binary formation channel, the observed number of events in the LVK
data places an upper bound on $\mu_{\lambda}$. Any models that output
$\lesssim N_{\rm obs}$ cannot be ruled out. However, even in the case where
the number of detected events is close to $N_{\rm obs}$, this is not a good
indicator of if the model population will match the observed population. We
propose checking the marginalised and unmarginalised likelihood simultaneously
to find if there is overlap in models. This will determine the models with
comparable populations and detected number of events to the LVK data. We check
the 30 largest likelihoods out of the 300 new models for both these cases
separately, and identify the common models and their hyper-parameters. We
again find a clear preference for no CHE for $f_{\rm WR}\geq 1$ and
$f_{\gamma}=0$, in that these models often return representative numbers of
detected events, and populations of masses and redshifts, but the choice of
remnant mass prescription varies between Mandel & Müller (2020), Fryer et al.
(2012) and Fryer et al. (2022).

####  5.2.2 Comparisons of the best individual models to data

We now compare the visual population parameter distributions to the GW events
with $\rm SNR\geq 8$ posteriors for the 300 new models. Model 44 has the best
likelihood ($\mathcal{L}(\lambda|\mathcal{H})$) which implies that the model
population best matches the observed LVK population. As shown in Fig. 17, we
find that the chirp mass and mass ratio distribution falls within the
2$\sigma$ contour of the GW events. The LVK chirp mass distribution has three
distinct peaks. The best model does not replicate the first peak but it does
for the second and third and generally matches the overall chirp mass
distribution. That said, the mass ratio distribution is uniform for the LVK
data but the model does prefer high mass ratios ($q>0.6$). Model 79 outputs
the worst likelihood. The chirp mass and mass ratio distributions have very
narrow peaks corresponding to the second peak in the LVK mass distribution
around $\mathcal{M}_{c}=20M_{\odot}$ and at $q=0.5$. The redshift distribution
range is also smaller compared to Model 44 and the data. However, the number
of detected events, $\mu_{\lambda}$ is closer to $N_{\rm obs}=79$. This
demonstrates the intrinsic difficulty in population synthesis of returning the
correct event-level properties and overall GW rate.

The best likelihood with a dependence on the number of events
($\mathcal{L}(\lambda,N_{\lambda}|\mathcal{H})$) is Model 274. This model has
a peak in the chirp mass distribution corresponding to the first peak in the
LVK mass distribution and recovers the higher mass peaks reasonably well too,
albeit with absence in the very high mass tail. The mass ratio distribution is
not uniform, there is a preference for mass ratios $>0.5$, but $<0.5$ are also
possible such that the model and data distributions agree reasonably well.
This model outputs $\mu_{\lambda}\approx 83$ events. The worst likelihood is
model 199. This model has has a rather high peak chirp mass around
$40M_{\odot}$ and a strong preference for high mass ratios ($q>0.8$). It is
likely that as a result of producing higher mass BHs with smaller separations,
this models outputs the largest detected number of events, $\sim 3805$, which
negatively impacts the likelihood.

Figure 17: $\mathcal{M}_{c}$, $q$ and $z$ $3\sigma$ distributions for the
models with the best (green contour) and worst (red contour) likelihood. On
the left we show the likelihood without the $N_{\lambda}$ dependence and the
right plot shows the likelihood with the $N_{\lambda}$ dependence. The blue
contour combines the posteriors for the $N_{\rm obs}=79$ GW events based on
the $\rm SNR\geq 8$ cut.

In all cases, our models seem to underpredict the high chirp mass and high
mass ratio events with $z\geq 0.5$ when compared to the LVK distribution. It
is possible that these events are not produced frequently by COMPAS. This may
be a limitation in the SSE implementation, as discussed in Section 2.4. It is
likely that these systems are not produced efficiently through the isolated
channel, which could indicate that some of our LVK events progenitors form
dynamically. And finally, these events may just be rare or ‘exceptional’ in
any formation channel/ environment. In this case, approaches more like
Passenger et al. (2024) — which uses the “normalised evidence” within a
statistical framework to determine if hierarchical mergers in active galactic
nuclei or globular clusters can justify exceptionally large events such as
GW190521 — may be more appropriate.

###  5.3 Which data are the most informative?

Given our range of models and the form with which the likelihood is evaluated
by summing over each GW event, we can also ask the question of which events
allow for the greatest distinction between our various COMPAS models. For each
GW event, we evaluate the distribution of $\log(\mathcal{I}^{k})$ terms for
all 300 models. We use the range and mean of the distributions to determine if
the progenitors of the GW events are distinguishable between models. We also
look for skewed distributions, which indicate that some models are more or
less likely to produce progenitors of particular GW events. We show an example
of 20 of these — the top with the least variance between models, and the ten
with the most — in Fig. 18.

Figure 18: Distribution of $\log(\mathcal{I}^{k})$ terms for the 300 new
COMPAS models. The top plot shows the GW events where the distribution are
constrained to high $\log\mathcal{I}^{k}$ values and the bottom plot shows GW
events with broad distributions. Out of the distinguishable events, the pink
distributions correspond to events with one component in the lower mass gap
and the green distributions correspond to events with at least one component
in the upper mass gap.

The top plot of Fig. 18 shows the GW events where their $\mathcal{I}^{k}$
terms are tightly constrained to high probabilities. This implies that all the
models are very likely to produce the progenitors of these events. The
distributions appear Gaussian. However, upon focusing on the
$\log(\mathcal{I}^{k})<-5$ region they appear to have a right skewed
distribution, implying there a few models in this tail that are less preferred
to form the progenitors of these events.

The bottom plot shows the GW events where the progenitors populations are more
clearly distinguishable between models. In particular, GW190814 and GW190521
have wide bimodal distributions, apparently separating our COMPAS models into
those that can, and those that cannot, reproduce the characteristics of these
events. In Figure 19 we show these events with low probability of their
progenitors forming in some models. They appear to be close to, or in, the
pair-instability mass gap $\sim 50-130M_{\odot}$ (Woosley & Heger, 2021).
GW190426_190642 has a total mass that supersedes the total mass of GW190521,
and it could be argued that they are from the same formation channel. However,
GW190426_190642 is less certain to be a true event, having $p_{\rm
astro}=0.75$ and $\rm FAR=4.1yr^{-1}$ (Banerjee, 2022). GW190602_175927 and
GW191109_010717 have at least one component in the mass gap. The spin-orbit
misalignment of GW191109_010717 provides further evidence it has a dynamical
origin (Zhang et al., 2023). However, GW191109_010717 has been flagged as
problematic due to data quality issues, which may lead to unreliable inference
of the source parameters such as the precession (Tong et al., 2022).
GW190917_114630 and GW200210_092254 have a secondary component in the lower
mass gap between the most massive neutron star and smallest BH, like GW190814.
The low probability of generating these events by some models indicate that
the isolated formation channel is not ideal for forming unequal mass binaries.
However, the ability of some models to form the progenitors of these events
suggests that this mass gap might be narrower or non-existent. More
observations of similar events is paramount in confirming this (Zevin et al.,
2020).

Overall, looking at the models that perform well for these events, they tend
to produce a dearth of higher mass BHs, with chirp mass distributions that
peak around $40M_{\odot}$ and mass ratio distributions that peak around
$q\approx 0.9$. However, this often comes at the expense of producing enough
lower mass or more extreme mass ratio binaries, to explain the mass and
redshift distribution of the other GW events.

This is particularly emphasised for GW190521 where we found that the posterior
samples often did not fall at all within the chirp mass or mass ratio range
output from COMPAS. In these cases we assign a probability of zero to these
events such that they make no difference to the sum in Eq 23 and 28, but are
also hence not represented in Fig. 18. Model choices that fail to produce any
progenitor close to GW190521 include:

  * •

No CHE (i.e., models 1-80 and 241-250), implying that formation channels with
mass transfer and CE evolution produce smaller BH masses.

  * •

The remnant mass prescriptions Mandel & Müller (2020) and Schneider et al.
(2021), which tend to produce smaller BH masses. The exception being model
271-280, which use Fryer et al. (2012) and Fryer et al. (2021).

  * •

Models 221 and 270, where CHE is set to Optimistic, but $f_{\gamma}=0.75$
(Model 221) and/or $f_{\rm WR}=10$ (Model 270), which may be contributing to
extreme mass loss from the progenitors again leading to a population of
particularly low-mass BHs.

Overall, the difficulty in simultaneously producing both these larger events
and the bulk of LVK data even across the large range of models we explore in
this work lends credence to the idea of multiple formation pathways for these
events.

Figure 19: Posterior distributions for ‘outlier’ events, (identified by the
green distributions in Fig. 18) that are not well reproduced by the majority
of our COMPAS models. The blue contour is the full posterior distribution for
LVK data with SNR $\geq$ 8, which clearly indicates the majority of the data
are at both lower mass and lower redshift. We show GW190426_190642 as an
unfilled contour, as it does not fit the $p_{\rm astro}\geq 0.9$ and $\rm
FAR\leq 0.25yr^{-1}$ criteria.

##  6 Conclusion

In this work we present a thorough analysis of the ability of over 300 binary
population synthesis models to reproduce the GW events observed so far within
Gravitational Wave Transients Catalogues 1-3. To enable this, we developed an
efficient method for predicting the BBH merger rate in the isolated channel
for various COMPAS models (Section 3), validated based on previous work from
Rauf et al. (2023). We then used this within a Bayesian inference method to
calculate the likelihood of these models with and without the dependence on
the number of GW events predicted by each model, $N_{\lambda}$ (Section 4).

The COMPAS population synthesis code contains a number of free parameters — we
introduce some of the most salient of these in Section 2. Our results fitting
these models to data in Section 5 demonstrate there is a preference for
subsets of models dictated by the type of isolated formation channel and
specific angular momentum loss. However, we frequently find opposing trends in
the likelihoods with and without $N_{\lambda}$ dependence — models that
produce realistic event-level masses and redshifts frequently over-predict the
number of events we should have observed to date. Looking at models which
perform well with both likelihoods, we find that models with no Chemically
Homogeneous Evolution, low to moderate values of $f_{\gamma}$ (such that
angular momentum loss occurs closer to the donor) and higher rates of mass-
loss from Wolf–Rayet winds (larger $f_{\rm WR}$) are generally preferred by
current data.

We find that the selection effects impact our likelihoods, as the $\rm SNR\geq
8$ cut on the GW events increases the likelihood for all models. This may vary
with future GW detections, particularly if the BH masses are in the mass gaps.
We also use the individual event likelihoods $I^{k}$ across all our models
term to distinguish individual events that are easy or more difficult to
reproduce by our models. We identify the usual suspects GW190521 and GW190814
as typically having low probability of existing in our new models, as well as
some outlier events in the the lower and upper mass gap. This motivates future
studies to consider joint population inference and model selection studies
across multiple formation channels (Gompertz et al., 2022; Gayathri et al.,
2023).

Our method for efficiently calculating the merger rate can also be applied to
other population synthesis codes. However, we note that it is only applicable
in the case where the metallicity dependent SFR is the same for the fiducial
and new model. We argue that using a different SFR may not impact the PDF
significantly. van Son et al. (2023) finds that variation in the SFR do not
impact the features in the primary BH mass distribution and have minimal
impact on the low mass end of the distribution.

We acknowledge this work relies on various assumptions and simplifications.
For future work, we consider including the effective spin as a population
parameter in the analysis. The covariance between the mass ratio and effective
spin can arise due to common envelope physics, mass transfer efficiency and
Eddington limited accretion (Mandel & Fragos, 2020; Bavera et al., 2021; Zevin
& Bavera, 2022; Broekgaarden et al., 2022b; Adamcewicz et al., 2023; Olejak et
al., 2024). This could be used to further distinguish between COMPAS models.
Alternatively, we can build the PDF only with $\mathcal{M}_{c}$, mass ratio or
redshift to test the validity of the model outputs. We could also consider
including natal kick velocity as one of the hyper-parameters to vary in COMPAS
in the future, as this is not well constrained for BHs (Callister et al.,
2021; Cui & Li, 2023; Vigna-Gómez et al., 2024), or more fundamental
parameters such as the binary fraction. Varying the selection function cuts
applied to current data to restrict the analysis to only the strongest
detections may further inform us on trends but given we found consistent
results with our two different selection functions this may not necessarily
provide improved constraints on the hyper-parameters.

Future GW detectors will observe GW events at a wider range of frequencies and
redshifts, and will provide data with increasing precision (Broekgaarden et
al., 2023). This will answer a variety of astrophysical questions and provide
a clearer picture of the BBH population distribution. In the absence of
electromagnetic counterparts, BBH merger populations with the correct
modelling are crucial for an unbiased inference of cosmological parameters
(Mastrogiovanni et al., 2021; Gray et al., 2023). With deeper catalogues and
higher SNR detections, our selection criteria will be more refined, which will
hopefully allow us to further constrain the uncertainties in massive binary
evolution models. However, as always, with more GW events comes greater change
of observing rare exceptional events to challenge these models. Until then, in
the case of limited number of GW events we must rely on observations of
massive binary stars to constrain the different assumptions for massive binary
evolution. Observations of X-ray binaries, red supergiants, luminous blue
variables and Wolf–Rayet stars can help constrain binaries at different
evolutionary stages (Eldridge & Stanway, 2022; Marchant & Bodensteiner, 2023;
Li et al., 2024).

## Acknowledgements

We thank Carole Perigois and Michella Mapelli for useful discussion on their
paper that motivated this work. We thank Ilya Mandel, Eric Thrane and Paul
Lasky for useful discussions on the methodology and results. We also thank
Reinhold Wilcox and Floor Broekgaarden for useful discussions on the COMPAS
hyper-parameters and relevant literature.

The authors acknowledge support from the Australian Research Council (ARC)
Centre of Excellence for Gravitational Wave Discovery (OzGrav), through
project numbers CE170100004 and CE230100016. LR and CH acknowledge support
through ARC Discovery Project DP20220101395. SS is a recipient of an ARC
Discovery Early Career Research Award (DE220100241).

## Data Availability

Prior and posterior samples of GW events were provided by Carole Perigois.
Simulations in this paper made use of the COMPAS rapid binary population
synthesis code (version v02.38.08), which is freely available at this
repository.

## References

  * Abbott et al. (2017) Abbott B. P., et al., 2017, Phys. Rev. Lett., 119, 161101
  * Abbott et al. (2019a) Abbott B. P., et al., 2019a, Physical Review X, 9, 031040
  * Abbott et al. (2019b) Abbott B. P., et al., 2019b, ApJ, 882, L24
  * Abbott et al. (2020a) Abbott R., et al., 2020a, Phys. Rev. Lett., 125, 101102
  * Abbott et al. (2020b) Abbott B. P., et al., 2020b, ApJ, 892, L3
  * Abbott et al. (2020c) Abbott R., et al., 2020c, ApJ, 896, L44
  * Abbott et al. (2021a) Abbott R., et al., 2021a, Physical Review X, 11, 021053
  * Abbott et al. (2021b) Abbott R., et al., 2021b, ApJ, 915, L5
  * Abbott et al. (2023a) Abbott R., et al., 2023a, Physical Review X, 13, 011048
  * Abbott et al. (2023b) Abbott R., et al., 2023b, Physical Review X, 13, 041039
  * Adamcewicz et al. (2023) Adamcewicz C., Lasky P. D., Thrane E., 2023, ApJ, 958, 13
  * Agrawal et al. (2020) Agrawal P., Hurley J., Stevenson S., Szécsi D., Flynn C., 2020, MNRAS, 497, 4549
  * Agrawal et al. (2022a) Agrawal P., Szécsi D., Stevenson S., Eldridge J. J., Hurley J., 2022a, MNRAS, 512, 5717
  * Agrawal et al. (2022b) Agrawal P., Stevenson S., Szécsi D., Hurley J., 2022b, A&A, 668, A90
  * Andres et al. (2022) Andres N., et al., 2022, Classical and Quantum Gravity, 39, 055002
  * Arca Sedda (2021) Arca Sedda M., 2021, ApJ, 908, L38
  * Arca Sedda et al. (2023) Arca Sedda M., Mapelli M., Benacquista M., Spera M., 2023, MNRAS, 520, 5259
  * Asplund et al. (2009) Asplund M., Grevesse N., Sauval A. J., Scott P., 2009, ARA&A, 47, 481
  * Banagiri et al. (2023) Banagiri S., Berry C. P. L., Davies G. S. C., Tsukada L., Doctor Z., 2023, Phys. Rev. D, 108, 083043
  * Banerjee (2022) Banerjee S., 2022, A&A, 665, A20
  * Barrett et al. (2018) Barrett J. W., Gaebel S. M., Neijssel C. J., Vigna-Gómez A., Stevenson S., Berry C. P. L., Farr W. M., Mandel I., 2018, MNRAS, 477, 4685
  * Bavera et al. (2021) Bavera S. S., et al., 2021, A&A, 647, A153
  * Belczynski et al. (2010) Belczynski K., Bulik T., Fryer C. L., Ruiter A., Valsecchi F., Vink J. S., Hurley J. R., 2010, ApJ, 714, 1217
  * Björklund et al. (2023) Björklund R., Sundqvist J. O., Singh S. M., Puls J., Najarro F., 2023, A&A, 676, A109
  * Bogomazov et al. (2018) Bogomazov A. I., Cherepashchuk A. M., Lipunov V. M., Tutukov A. V., 2018, New Astron., 58, 33
  * Bouffanais et al. (2019) Bouffanais Y., Mapelli M., Gerosa D., Di Carlo U. N., Giacobbo N., Berti E., Baibhav V., 2019, ApJ, 886, 25
  * Bouffanais et al. (2021) Bouffanais Y., Mapelli M., Santoliquido F., Giacobbo N., Iorio G., Costa G., 2021, MNRAS, 505, 3873
  * Breivik et al. (2020) Breivik K., et al., 2020, ApJ, 898, 71
  * Broekgaarden et al. (2022a) Broekgaarden F. S., et al., 2022a, MNRAS, 516, 5737
  * Broekgaarden et al. (2022b) Broekgaarden F. S., Stevenson S., Thrane E., 2022b, ApJ, 938, 45
  * Broekgaarden et al. (2023) Broekgaarden F. S., Banagiri S., Payne E., 2023, arXiv e-prints, p. arXiv:2303.17628
  * Cabourn Davies & Harry (2022) Cabourn Davies G. S., Harry I. W., 2022, Classical and Quantum Gravity, 39, 215012
  * Callister et al. (2017) Callister T. A., Kanner J. B., Massinger T. J., Dhurandhar S., Weinstein A. J., 2017, Classical and Quantum Gravity, 34, 155007
  * Callister et al. (2021) Callister T. A., Farr W. M., Renzo M., 2021, ApJ, 920, 157
  * Cheng et al. (2023a) Cheng A. Q., Zevin M., Vitale S., 2023a, ApJ, 955, 127
  * Cheng et al. (2023b) Cheng A. Q., Zevin M., Vitale S., 2023b, Astrophys. J., 955, 127
  * Chruślińska (2024) Chruślińska M., 2024, Annalen der Physik, 536, 2200170
  * Claret & Torres (2018) Claret A., Torres G., 2018, ApJ, 859, 100
  * Costa et al. (2019) Costa G., Girardi L., Bressan A., Marigo P., Rodrigues T. S., Chen Y., Lanza A., Goudfrooij P., 2019, MNRAS, 485, 4641
  * Costa et al. (2023) Costa G., Chruślińska M., Klencki J., Broekgaarden F. S., Rodriguez C. L., Joseph T. D., Saracino S., 2023, arXiv e-prints, p. arXiv:2311.15778
  * Cui & Li (2023) Cui Z., Li X.-D., 2023, MNRAS, 523, 5565
  * Delfavero et al. (2023) Delfavero V., O’Shaughnessy R., Belczynski K., Drozda P., Wysocki D., 2023, Phys. Rev. D, 108, 043023
  * Detmers et al. (2008) Detmers R. G., Langer N., Podsiadlowski P., Izzard R. G., 2008, A&A, 484, 831
  * Dorozsmai & Toonen (2022) Dorozsmai A., Toonen S., 2022, arXiv e-prints, p. arXiv:2207.08837
  * Eldridge & Stanway (2022) Eldridge J. J., Stanway E. R., 2022, ARA&A, 60, 455
  * Ertl et al. (2016) Ertl T., Janka H. T., Woosley S. E., Sukhbold T., Ugliano M., 2016, ApJ, 818, 124
  * Farr et al. (2015) Farr W. M., Gair J. R., Mandel I., Cutler C., 2015, Phys. Rev. D, 91, 023005
  * Fishbach & Kalogera (2022) Fishbach M., Kalogera V., 2022, ApJ, 929, L26
  * Fishbach et al. (2018) Fishbach M., Holz D. E., Farr W. M., 2018, ApJ, 863, L41
  * Fryer et al. (2012) Fryer C. L., Belczynski K., Wiktorowicz G., Dominik M., Kalogera V., Holz D. E., 2012, The Astrophysical Journal, 749, 91
  * Fryer et al. (2021) Fryer C. L., Karpov P., Livescu D., 2021, Astronomy Reports, 65, 937
  * Fryer et al. (2022) Fryer C. L., Olejak A., Belczynski K., 2022, ApJ, 931, 94
  * Galaudage et al. (2020) Galaudage S., Talbot C., Thrane E., 2020, Phys. Rev. D, 102, 083026
  * Gallegos-Garcia et al. (2021) Gallegos-Garcia M., Berry C. P. L., Marchant P., Kalogera V., 2021, ApJ, 922, 110
  * Gallegos-Garcia et al. (2023) Gallegos-Garcia M., Jacquemin-Ide J., Kalogera V., 2023, arXiv e-prints, p. arXiv:2308.13146
  * Gayathri et al. (2023) Gayathri V., Bartos I., Rosswog S., Miller M. C., Veske D., Lu W., Marka S., 2023, arXiv e-prints, p. arXiv:2307.09097
  * Ge et al. (2010) Ge H., Hjellming M. S., Webbink R. F., Chen X., Han Z., 2010, ApJ, 717, 724
  * Ge et al. (2015) Ge H., Webbink R. F., Chen X., Han Z., 2015, ApJ, 812, 40
  * Ge et al. (2020) Ge H., Webbink R. F., Chen X., Han Z., 2020, ApJ, 899, 132
  * Gompertz et al. (2022) Gompertz B. P., Nicholl M., Schmidt P., Pratten G., Vecchio A., 2022, MNRAS, 511, 1454
  * Gray et al. (2023) Gray R., et al., 2023, J. Cosmology Astropart. Phys., 2023, 023
  * Grichener (2023) Grichener A., 2023, MNRAS, 523, 221
  * Hirai & Mandel (2022) Hirai R., Mandel I., 2022, ApJ, 937, L42
  * Hurley et al. (2000) Hurley J. R., Pols O. R., Tout C. A., 2000, MNRAS, 315, 543
  * Hurley et al. (2002) Hurley J. R., Tout C. A., Pols O. R., 2002, MNRAS, 329, 897
  * Ivanova et al. (2013) Ivanova N., et al., 2013, A&ARv, 21, 59
  * Janka (2012) Janka H.-T., 2012, Annual Review of Nuclear and Particle Science, 62, 407
  * Janka (2017) Janka H.-T., 2017, in Alsabti A. W., Murdin P., eds, , Handbook of Supernovae. p. 1095, doi:10.1007/978-3-319-21846-5_109
  * Janka et al. (2007) Janka H. T., Langanke K., Marek A., Martínez-Pinedo G., Müller B., 2007, Phys. Rep., 442, 38
  * Kaaret et al. (2004) Kaaret P., Alonso-Herrero A., Gallagher J. S., Fabbiano G., Zezas A., Rieke M. J., 2004, MNRAS, 348, L28
  * Klencki et al. (2021) Klencki J., Nelemans G., Istrate A. G., Chruslinska M., 2021, A&A, 645, A54
  * Kroupa (2002) Kroupa P., 2002, Science, 295, 82
  * LIGO Scientific Collaboration et al. (2018) LIGO Scientific Collaboration Virgo Collaboration KAGRA Collaboration 2018, LVK Algorithm Library - LALSuite, Free software (GPL), doi:10.7935/GT1W-FZ16
  * Lagos et al. (2018) Lagos C. d. P., Tobar R. J., Robotham A. S. G., Obreschkow D., Mitchell P. D., Power C., Elahi P. J., 2018, MNRAS, 481, 3573
  * Lagos et al. (2019) Lagos C. d. P., et al., 2019, MNRAS, 489, 4196
  * Laplace et al. (2020) Laplace E., Götberg Y., de Mink S. E., Justham S., Farmer R., 2020, A&A, 637, A6
  * Li et al. (2024) Li Z., Zhu C., Lü G., Li L., Liu H., Guo S., Yu J., Lu X., 2024, arXiv e-prints, p. arXiv:2405.11571
  * Liotine et al. (2023) Liotine C., Zevin M., Berry C. P. L., Doctor Z., Kalogera V., 2023, ApJ, 946, 4
  * Lu et al. (2023) Lu W., Fuller J., Quataert E., Bonnerot C., 2023, MNRAS, 519, 1409
  * Lynch et al. (2018) Lynch R., Coughlin M., Vitale S., Stubbs C. W., Katsavounidis E., 2018, ApJ, 861, L24
  * MacLeod & Loeb (2020) MacLeod M., Loeb A., 2020, ApJ, 895, 29
  * MacLeod et al. (2018a) MacLeod M., Ostriker E. C., Stone J. M., 2018a, ApJ, 863, 5
  * MacLeod et al. (2018b) MacLeod M., Ostriker E. C., Stone J. M., 2018b, ApJ, 868, 136
  * Mandel & Broekgaarden (2021) Mandel I., Broekgaarden F. S., 2021, arXiv preprint arXiv:2107.14239
  * Mandel & Broekgaarden (2022) Mandel I., Broekgaarden F. S., 2022, Living Reviews in Relativity, 25, 1
  * Mandel & Farmer (2022) Mandel I., Farmer A., 2022, Phys. Rep., 955, 1
  * Mandel & Fragos (2020) Mandel I., Fragos T., 2020, ApJ, 895, L28
  * Mandel & Müller (2020) Mandel I., Müller B., 2020, Monthly Notices of the Royal Astronomical Society, 499, 3214
  * Mandel & de Mink (2016) Mandel I., de Mink S. E., 2016, MNRAS, 458, 2634
  * Mandel et al. (2019) Mandel I., Farr W. M., Gair J. R., 2019, MNRAS, 486, 1086
  * Mapelli (2019) Mapelli M., 2019, in Oskinova L. M., Bozzo E., Bulik T., Gies D. R., eds, , Vol. 346, High-mass X-ray Binaries: Illuminating the Passage from Massive Binaries to Merging Compact Objects. pp 397–416, doi:10.1017/S1743921319001868
  * Mapelli (2020) Mapelli M., 2020, Frontiers in Astronomy and Space Sciences, 7, 38
  * Mapelli (2021) Mapelli M., 2021, in , Handbook of Gravitational Wave Astronomy. p. 16, doi:10.1007/978-981-15-4702-7_16-1
  * Marchant & Bodensteiner (2023) Marchant P., Bodensteiner J., 2023, arXiv e-prints, p. arXiv:2311.01865
  * Marchant et al. (2016) Marchant P., Langer N., Podsiadlowski P., Tauris T. M., Moriya T. J., 2016, A&A, 588, A50
  * Marchant et al. (2021) Marchant P., Pappas K. M. W., Gallegos-Garcia M., Berry C. P. L., Taam R. E., Kalogera V., Podsiadlowski P., 2021, A&A, 650, A107
  * Marchant et al. (2023) Marchant P., Podsiadlowski P., Mandel I., 2023, arXiv e-prints, p. arXiv:2311.14041
  * Marek & Janka (2009) Marek A., Janka H. T., 2009, ApJ, 694, 664
  * Mastrogiovanni et al. (2021) Mastrogiovanni S., et al., 2021, Phys. Rev. D, 104, 062009
  * Miller-Jones et al. (2021) Miller-Jones J. C. A., et al., 2021, Science, 371, 1046
  * Müller et al. (2016) Müller B., Heger A., Liptai D., Cameron J. B., 2016, MNRAS, 460, 742
  * Nathaniel (2017) Nathaniel K., 2017
  * Neijssel et al. (2019) Neijssel C. J., et al., 2019, MNRAS, 490, 3740
  * Neijssel et al. (2021) Neijssel C. J., Vinciguerra S., Vigna-Gómez A., Hirai R., Miller-Jones J. C. A., Bahramian A., Maccarone T. J., Mandel I., 2021, ApJ, 908, 118
  * Neo et al. (1977) Neo S., Miyaji S., Nomoto K., Sugimoto D., 1977, PASJ, 29, 249
  * O’Connor & Ott (2011) O’Connor E., Ott C. D., 2011, ApJ, 730, 70
  * Olejak et al. (2024) Olejak A., Klencki J., Xu X.-T., Wang C., Belczynski K., Lasota J.-P., 2024, arXiv e-prints, p. arXiv:2404.12426
  * Orosz et al. (2007) Orosz J. A., et al., 2007, Nature, 449, 872
  * Paczyński & Sienkiewicz (1972) Paczyński B., Sienkiewicz R., 1972, Acta Astron., 22, 73
  * Passenger et al. (2024) Passenger L., Thrane E., Lasky P. D., Payne E., Stevenson S., Farr B., 2024, arXiv e-prints, p. arXiv:2405.09739
  * Paxton et al. (2011) Paxton B., Bildsten L., Dotter A., Herwig F., Lesaffre P., Timmes F., 2011, ApJS, 192, 3
  * Périgois et al. (2023) Périgois C., Mapelli M., Santoliquido F., Bouffanais Y., Rufolo R., 2023, arXiv e-prints, p. arXiv:2301.01312
  * Picco et al. (2023) Picco A., Marchant P., Sana H., Nelemans G., 2023, arXiv e-prints, p. arXiv:2309.05736
  * Picker et al. (2024) Picker L., Hirai R., Mandel I., 2024, arXiv e-prints, p. arXiv:2402.13180
  * Pols et al. (1998) Pols O. R., Schröder K.-P., Hurley J. R., Tout C. A., Eggleton P. P., 1998, MNRAS, 298, 525
  * Powell et al. (2016) Powell J., Gossan S. E., Logue J., Heng I. S., 2016, Phys. Rev. D, 94, 123012
  * Powell et al. (2023) Powell J., Iess A., Llorens-Monteagudo M., Obergaulinger M., Müller B., Torres-Forné A., Cuoco E., Font J. A., 2023, arXiv e-prints, p. arXiv:2311.18221
  * Rauf et al. (2023) Rauf L., Howlett C., Davis T. M., Lagos C. D. P., 2023, Monthly Notices of the Royal Astronomical Society, 523, 5719
  * Riley & Mandel (2023) Riley J., Mandel I., 2023, ApJ, 950, 80
  * Riley et al. (2021) Riley J., Mandel I., Marchant P., Butler E., Nathaniel K., Neijssel C., Shortt S., Vigna-Gómez A., 2021, MNRAS, 505, 663
  * Romagnolo et al. (2023) Romagnolo A., Belczynski K., Klencki J., Agrawal P., Shenar T., Szécsi D., 2023, MNRAS, 525, 706
  * Romero-Shaw et al. (2022) Romero-Shaw I., Lasky P. D., Thrane E., 2022, ApJ, 940, 171
  * Romero-Shaw et al. (2023) Romero-Shaw I., Hirai R., Bahramian A., Willcox R., Mandel I., 2023, MNRAS, 524, 245
  * Sana et al. (2012) Sana H., et al., 2012, Science, 337, 444
  * Santoliquido et al. (2021) Santoliquido F., Mapelli M., Giacobbo N., Bouffanais Y., Artale M. C., 2021, MNRAS, 502, 4877
  * Santoliquido et al. (2024) Santoliquido F., Mapelli M., Iorio G., Costa G., Glover S. C. O., Hartwig T., Klessen R. S., Merli L., 2024, MNRAS, 528, 954
  * Schneider et al. (2015) Schneider F. R. N., Izzard R. G., Langer N., de Mink S. E., 2015, ApJ, 805, 20
  * Schneider et al. (2021) Schneider F., Podsiadlowski P., Müller B., 2021, Astronomy & Astrophysics, 645, A5
  * Schneider et al. (2023) Schneider F. R. N., Podsiadlowski P., Laplace E., 2023, ApJ, 950, L9
  * Schneider et al. (2024) Schneider F. R. N., Podsiadlowski P., Laplace E., 2024, arXiv e-prints, p. arXiv:2403.03984
  * Schrøder et al. (2021) Schrøder S. L., MacLeod M., Ramirez-Ruiz E., Mandel I., Fragos T., Loeb A., Everson R. W., 2021, arXiv e-prints, p. arXiv:2107.09675
  * Silverman (1986) Silverman B. W., 1986, Density estimation for statistics and data analysis
  * Stevenson & Clarke (2022) Stevenson S., Clarke T. A., 2022, MNRAS, 517, 4034
  * Stevenson et al. (2015) Stevenson S., Ohme F., Fairhurst S., 2015, ApJ, 810, 58
  * Stevenson et al. (2017) Stevenson S., Vigna-Gómez A., Mandel I., Barrett J. W., Neijssel C. J., Perkins D., de Mink S. E., 2017, Nature Communications, 8, 14906
  * Tanikawa (2024) Tanikawa A., 2024, arXiv e-prints, p. arXiv:2403.04389
  * Team COMPAS: Riley, J. et al. (2022) Team COMPAS: Riley, J. et al., 2022, ApJS, 258, 34
  * Temaj et al. (2024) Temaj D., Schneider F. R. N., Laplace E., Wei D., Podsiadlowski P., 2024, A&A, 682, A123
  * Temmink et al. (2023) Temmink K. D., Pols O. R., Justham S., Istrate A. G., Toonen S., 2023, A&A, 669, A45
  * Tong et al. (2022) Tong H., Galaudage S., Thrane E., 2022, Phys. Rev. D, 106, 103019