# Forward Modelling the O3(a+b) GW transient mass distributions with bpass by
varying compact remnant mass and SNe kick prescriptions

Sohan Ghodla1, Wouter G. J. van Zeist1, J. J. Eldridge1, Héloïse F. Stevance1,
Elizabeth R. Stanway2
1Department of Physics, University of Auckland, Private Bag 92019, Auckland,
New Zealand
2Department of Physics, University of Warwick, Gibbet Hill Road, Coventry, CV4
7AL, UK<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

We present forward modelling from the bpass code suite of the population of
observed gravitational wave (GW) transients reported by the LIGO/VIRGO
consortium (LVC) during their third observing run, O3(a+b). Specifically, we
predict the expected chirp mass and mass ratio distributions for GW
transients, taking account of detector sensitivity to determine how many
events should have been detected by the current detector network in O3(a+b).
We investigate how these predictions change by alternating between four
different remnant mass estimation schemes and two supernovae (SNe) kick
prescriptions. We find that none of the model populations resulting from these
variations accurately match the whole O3(a+b) GW transient catalog. However,
agreement from some models to part of the catalog suggests ways to achieve a
more complete fit. These include reducing the number of low mass black holes
(BHs) close to the mass gap, while also increasing the number of higher mass
BHs below the pair-instability SN limit. Finally, we find that the interaction
between the value of the remnant mass from a stellar model and the choice of
SN kick is complex and different kick prescriptions may be required depending
on whether a neutron star or BH is formed.

###### keywords:

stars: evolution – binaries: general – supernovae: general – methods:
numerical

††pubyear: 2021††pagerange: Forward Modelling the O3(a+b) GW transient mass
distributions with bpass by varying compact remnant mass and SNe kick
prescriptions–A

##  1 Introduction

In the span of only almost half a decade, the detection of GW transients has
become almost routine. After LVC O3(a+b), there are now 90111Overall, 11
events from O1+O2, 44 events from O3a and 35 events from O3b known GW
transients events with a probability of being a real astrophysical signal
($p_{\rm astro})>50\%$ (Abbott et al., 2020a; The LIGO Scientific
Collaboration et al., 2021b). Each event in itself reveals much about the
evolution of binary stars (e.g. for individual studies see Abbott et al.
(2016, 2017, 2020b) and for evolutionary implications for binaries, see
Eldridge & Stanway (2016); Belczynski et al. (2016); Stevenson et al. (2017);
Mapelli et al. (2017); Vigna-Gómez et al. (2018); Eldridge et al. (2019);
Lipunov et al. (2021); Broekgaarden et al. (2021) and references therein).
However with a growing understanding of the distribution of masses, we are now
beginning to constrain the population of GW transients both through parametric
models (e.g. Abbott et al. 2019, 2020a; The LIGO Scientific Collaboration et
al. 2021b) and using the technique of population synthesis to understand the
physical process of stellar evolution that gives rise to the mass distribution
(e.g. Belczynski et al. (2008); Stevenson et al. (2017); Eldridge et al.
(2017); Giacobbo et al. (2018) and references therein).

While binary population synthesis is a powerful tool, there are many
uncertainties in the physics incorporated in these models (e.g. De Marco &
Izzard 2017) as well as in the prescriptions used to model the star formation
history and chemical evolution of the Universe (Belczynski et al. 2016;
Mapelli et al. 2017; Chruslinska et al. 2019; Tang et al. 2020; Boco et al.
2021; Broekgaarden et al. 2021). In particular, apart from the uncertainties
in evolution of massive single stars (eg. mass-loss rate, rotation,
overshooting, etc.) the understanding of the mass-transfer and common envelope
phase in binaries is also incomplete (e.g., Ivanova et al. 2013 and references
therein). In addition, though the delaytime distribution is believed to mostly
lie within 0-4.5 Gyrs (e.g. Fishbach & Kalogera 2021), a subset of GW mergers
has also been found to come from very long time delays (e.g., Eldridge &
Stanway 2016; Mapelli et al. 2017). As such, the wide distribution of possible
delaytime intervals between the formation of the progenitor stellar binary of
a GW transient and the inspiral event introduces sensitivity to past star
formation conditions. Another required assumption is the remnant mass
prescription used to determine the compact remnant mass once the stellar model
reaches its final state, which often does not proceed all the way to core-
collapse for numerical reasons. To this end, several remnant mass
prescriptions with different features also exist in the literature, e.g. see,
Patton & Sukhbold (2020).

In this letter, we utilize bpass (Binary Population and Spectral Synthesis)
code results to forward model the population of GW transients resulting from
their underlying stellar population up till their chirp event, specifically
predicting the distribution of expected chirp masses and mass ratios for them
during LVC O3(a+b). Section 2 outlines the numerical model, highlighting the
differences between this and earlier bpass work. Specifically, here we vary
two key uncertainties in our analysis, (i) the nature of the kick that
remnants receive while undergoing a SN explosion and (ii) the prescription for
calculating the resulting remnant masses. Section 3 presents our results on
the number of expected detections and the predicted GW transient object mass
distribution. We then compare these to the LVC O3(a+b) catalog (The LIGO
Scientific Collaboration et al., 2021b) to see which of the prescriptions
produce a mass distribution that closely matches their results, presenting
discussion and concluding remarks in section 4.

##  2 Methods and Prescriptions, Simulations and Observations

Our aim here is to forward model the population of GW transients that may have
been detected in O3(a+b) (henceforth O3), making predictions independent of
the O3 observations, and explore some of the stellar evolutionary
uncertainties that might impact these results. We use four different schemes
to estimate the remnant masses resulting from stellar evolution and two
different kicks to investigate how they collectively affect the populations of
the predicted GW mergers.

Figure 1: The detection distance versus chirp mass as a function of detection
probability and mass ratio. The solid, dotted and dashed lines represent the
mass ratio values $q=1$, $q=0.5$ and $q=0.1$. The black lines are for a
detection probability of 0.9, red 0.5 and blue 0.1.

###  2.1 Population synthesis

We use bpass v2.2.1 stellar models (Eldridge et al., 2017; Stanway & Eldridge,
2018) and build upon our earlier work on predicting the GW transient rates. We
refer the reader to those papers for full details of our model and method
(Eldridge & Stanway, 2016; Eldridge et al., 2019; Tang et al., 2020). We use
the same underlying model of the cosmological star formation history and
metallicity evolution as in Eldridge et al. (2019). However, two aspects of
our previous population synthesis have been improved.

Firstly, we have updated the treatment of rejuvenation during mass transfer
onto the secondary stars. bpass has two principal sets of models: primary
models where the more massive star’s full structure is computed in detail and
the secondary’s evolution is approximated using the equations of Hurley et al.
(2000); and secondary models where the secondary is evolved in detail as
either a single star or a star in orbit around a compact object – depending on
the outcome of the first SN. Many different primaries map onto the same
secondary model. As a result, in earlier bpass calculations we used an average
rejuvenation age (i.e. time step at which mass transfer leads to rotational
mixing, causing the accretor to become chemically homogeneous eg., Maeder
1987) for each secondary model. In this work, we use a distribution of
rejuvenation ages for each secondary model. This has the effect of smoothing
out our delaytime distribution for transient events driven by GW emission.

Secondly, we now compute the GW merger time by calculating the orbital
evolution in full as outlined by Peters (1964) rather than by simple
interpolation between the two analytic forms for circular and highly eccentric
orbits. This is motivated by the fact that many of the orbits lie in the mid-
range between circular and eccentric after the second supernova, and the
updated calculations offer a more accurate prediction of the merger times.

###  2.2 Detection horizons

Detection of predicted mergers will depend on instrumental and alignment
effects. Suites of codes to evaluate this do exist, such as LALSuite (LIGO
Scientific Collaboration, 2018), but we implement our own Python package to do
this because of our specific focus on population synthesis. The package
Riroriro222A riroriro is a bird native to Aotearoa New Zealand that is
frequently heard but seldom seen. has been created with the aim of determining
the fraction of predicted GW transients that can be detected (van Zeist et
al., 2021). Currently, it only deals with GW transients caused by the merger
of NSNS, NS-BH, and BHBH. In the future, we plan to extend its scope to all
areas of the GW population and the resulting GW spectral synthesis.

For each chirp mass and mass ratio pair Riroriro creates a synthetic merger
gravitational waveform (including inspiral and ringdown phases). The effects
of source orientation and distance are then taken into account to calculate an
effective detection volume and the evolution of detection fraction of these
events with luminosity distance. Our calculations span from 1 Mpc to 10 Gpc in
steps of 0.1 dex. The package is inspired from Buskirk & Babiuc-Hamilton
(2019) who created a Mathematica routine to calculate the expected GW waveform
for a GW transient event. They used the results of Huerta et al. (2017) to
model the inspiral and ringdown emission. We use the same methods they outline
as well as using information from the Findchirp algorithm (Allen et al., 2012)
to improve our modeling of the inspiral waveform and calculation of the
signal-to-noise ratio (SNR). It is often difficult to match the inspiral and
ringdown waveforms and we switch between them when the frequency and time
derivative are closest between the two models.

For each resultant waveform, we use the method of Barrett et al. (2018) to
calculate the SNR of the merger model at a given redshift using the
corresponding luminosity distance and redshifting the waveform. We
subsequently calculate a total SNR assuming a triple-detection by combining
the SNR calculated for the LIGO Hanford, LIGO Livingston, and VIRGO detectors
during O3 in quadrature. We also calculate the average expected SNR of a GW
observation with a single, non-specific detector operating. These SNR values
assumed optimal alignment of detector and source. To take into account the
effect of arbitrary orientation we multiplied these with a projection function
(Finn, 1996; Belczynski et al., 2013, 2014), giving a probability distribution
of the expected SNR and the probability that this would exceed the commonly-
used detectability threshold value of 8, i.e. the detection probability. The
variation in the detection probability of a GW transient depending on the
chirp mass, mass ratio, and distance of the binary is illustrated in Fig. 1.
For each chirp mass, there is a maximum distance beyond which no mergers can
be detected but there is a smoothly decreasing detection probability up to
this horizon. Higher chirp masses can be detected to larger distances.
However, at the highest chirp masses a hard upper limit is introduced due to
the decreasing maximum frequency of these mergers. This upper limit is most
sensitive to the mass ratio, as shown in Fig. 1, and the lowest mass ratio has
a significantly lower maximum chirp mass for the peak comoving distance to the
source. We combine the detection probabilities as a function of chirp mass,
mass ratio and distance with the population statistics from bpass, which
dictate the frequency with which systems of such parameters form, to create a
forward model for the LVC detection rates. We integrate the rates in events
per year per unit volume over the observable Universe to obtain rates in
events per year. We use the method of Hogg (1999) to convert the luminosity
distances used in the previous steps to comoving distances which are employed
in the volume calculations333We assume a standard $\Lambda$CDM cosmology with
$H_{0}=100$h kms-1 Mpc-3 with h = 0.696, $\Omega_{\mathrm{M}}$ = 0.286,
$\Omega_{\Lambda}$ = 0.714. We reduce our rates by a factor of (177.3 + 142.0)
days / 365.25 days) to account for the elapsed duration during O3a + O3b where
at least one detector was active and taking data.

###  2.3 Remnant mass estimates and supernova kick velocities

Within bpass the standard method of calculating remnant masses at the end of
stellar evolution for massive stars is to compare the binding energy of the
stellar envelope to the typical SN explosion energy of $10^{51}\,{\rm ergs}$.
The ejecta mass is considered to be the amount of mass with this much binding
energy, removed from the star, and the remainder goes into the remnant mass
(Eldridge & Tout, 2004). Here we refer to this method as the Standard remnant
mass scheme. However, to explore the effect of varying the remnant mass
prescription on our predictions, we look into three other schemes.

(i) In AlwaysNS scheme, motivated by the Chandrasekhar mass limit (and to save
computation time), we force all the models to form a neutron star (NS) of 1.4
M⊙ at core-collapse. Forcing all the remnants to form a NS at birth (albeit of
the same mass) allows us to observe the influence this could have on the NS
mergers rates. (ii) For most stars the final mass of the carbon-oxygen (CO)
core would be a fair estimate of the maximum possible remnant mass, hence in
MCO,final remnant mass scheme, the final CO core mass of the progenitor model
is set equal to the remnant mass (i.e. only the hydrogen/helium envelope will
be ejected, independent of SN energy). (iii) For comparative purpose we employ
a commonly-used rapid synthesis scheme from Fryer et al. (2012), here referred
as FryerRapid. This scheme gives a remnant mass that depends on MCO,final (and
partly on Mfinal) of the stellar model. One interesting feature is that it
predicts a mass-gap in the BH mass distribution with no BHs less massive than
5 M⊙. Consequently, using the FryerRapid scheme allows us to determine the
impact of allowing for a mass-gap. Finally, in addition to using the above
remnant mass prescriptions, we also use two different SN kick schemes namely
the Hobbs et al. (2005) and Bray & Eldridge (2018) kick velocity relations.
The former is described by a Maxwell-Boltzmann distribution with 1-dimensional
root mean square speed of 265 kms-1 while the latter depends strongly by
construction on the masses of the ejecta and the remnant and thus we expect
the kick and remnant mass distributions to interact in a complex and non-
linear way. Mathematically, the Bray kick relation is as follows

$v_{\mathrm{kick}}=100\left(\frac{M_{\mathrm{ejecta}}}{M_{\mathrm{remnant}}}\right)-170\left(\frac{1.4}{M_{\mathrm{remnant}}}\right)$
(1)

where $v_{\rm kick}$ is the velocity, in km s-1, imparted to the compact
remnant as a result of the SN, $M_{\mathrm{ejecta}}$ is the SN ejecta mass and
$M_{\mathrm{remnant}}$ is the remnant mass (all masses in [M⊙] units).

##  3 Results

Below we show some ways in which one could constrain or verify the form of the
SNe kick and the remnant mass relation using the GW transient data from O3.

Figure 2: The Chirp mass distributions of cumulative events observed by LVC
O3(a+b) (black line) and those predicted by various schemes under the bpass
model populations for the same duration of observation. The green ticks show
the observed values of Mchirp and the solid grey area represents the
uncertainty between the upper and lower Mchirp bounds for these observed
events (Abbott et al. 2020a; The LIGO Scientific Collaboration et al. 2021b).
The black curve is the cumulative observed number of events over the observed
values of Mchirp. Solid lines represent Bray kick while dashed lines Hobbs
kick.

###  3.1 Some implications based on merger rates

Table 1:  Predicted number of observable detections during O3(a+b) for the three expected mergers types, employing various methods for calculating SN kicks and remnant masses. Rates in parenthesis are local intrinsic merger rates (in units of events Gpc-3 yr-1) calculated before the consideration of detection probabilities. Last row states the LVC O3(a+b) observed and intrinsic rates. (Abbott et al., 2020a; The LIGO Scientific Collaboration et al., 2021b). The detectable rates account for the duration of 177.3 + 142.0 days of LVC O3 during which at least one detector was active. Kick | Remnant | NSNS | BHNS | BHBH
---|---|---|---|---
Hobbs | Standard | 4.9 (208) | 12.6 (131) | 45 (31)
Hobbs | M${}_{\mathrm{CO,\;final}}$ | 1.6 (43) | 46.8 (417) | 420 (873)
Hobbs | FryerRapid | 3.8 (176) | 27 (208) | 171 (169)
Hobbs | AlwaysNS | 4.9 (223) | 0.4 (8.7) | 0.0
Bray | Standard | 17.1 (745) | 13.3 (180) | 27.0 (13)
Bray | M${}_{\mathrm{CO,\,final}}$ | 5.6 (179) | 43.2 (498) | 212.5 (569)
Bray | FryerRapid | 13.7 (677) | 16.0 (157) | 64.8 (70)
Bray | AlwaysNS | 14.4 (708) | 2 (48) | 0.0
LVC O3(a+b) | 1$\pm$1 | 5$\pm$2.2 | 69$\pm$8.3
Intrinsic [Gpc-3 yr-1] | (13–1900) | (7.4–320) | (17.3–45)

Our results shown in Table 1 (also see Fig. 4, 5 in the Appendix), presents
the expected number of events that should have been detected in O3 along with
the predicted local intrinsic rates under each set of prescriptions. On
considering these under the light of O3 merger rate data, we observe patterns
with interesting implications:

  1. 1.

The Bray kick consistently gives a low BHBH intrinsic and detectable merger
rate, as compared to the Hobbs kick. The strength of these differences varies
between the remnant mass schemes but a general trend holds. This is acceptable
as in the Bray kick, vkick is inversely proportional to the Mremnant. Thus
more massive remnants receive weaker kicks and so have less eccentric orbits,
longer merger times and a lower BHBH merger rate. In particular, we find that
for very massive BHs the eccentricity of the orbit is independent of the kick
but for intermediate mass BHs, Hobbs kick produces more eccentric orbits and
therefore more mergers. This behavior can be seen in Fig. 2 where the
cumulative observable count for both kicks closely follow each other up till
$\sim$ 1.35 ($\sim$ 22 M⊙) from the left, after which Hobbs kick starts to
produce more observable detection for the BBH systems.

  2. 2.

In MCO,final, the BHBH and BHNS intrinsic and detected rates are significantly
higher than those seen in O3 implying that under our current settings, stellar
remnants must be consistently less massive than the MCO,final of their
progenitors. Fig. 9 (discussed later) also demonstrates this bias of MCO,final
toward higher mass remnants.

  3. 3.

The Standard scheme over-predicts the number of detectable NSNS and BHNS
mergers while under-predicting the BHBH mergers (fig 4), more so using the
Bray kick. On the other hand, a better agreement is seen once their intrinsic
merger rates are compared with the expected local intrinsic rate.
Additionally, the O3 BHBH merger count is reasonably reproduced by the Bray
kick \+ FryerRapid remnant scheme. However, this produces too many mergers
involving binary NSs. Though none of the schemes fully satisfy either the
detectability or intrinsic criteria set by LVC O3, the Standard \+ Hobbs and
Fryer Rapid + Bray offer the closest fit to the data in Table 1. This can also
be seen in figure 6 which provides the local BHBH merger rate density
distribution as a function of the primary BH’s mass. Both the Standard \+
Hobbs and Fryer Rapid + Bray produce mergers owing only to the isolated binary
evolution scenario and their rates mostly fall below the 90% credibility
interval. Other relations also show similar trends but the above two schemes
along with Fryer Rapid + Hobbs come the closest to the PowerLaw+Peak
distribution. We note that though Fig. 6 helps us to perform a check on the
validity of the chosen schemes it is not yet possible to reach a definitive
conclusion.

  4. 4.

The AlwaysNS scheme does have some BHNS mergers despite no BHs being formed in
core-collapse. These BHs are the result of the subsequent accretion-induced
collapse (due to accretion from the companion) of the already formed NS. This
gives an indication of the relative importance of such a channel compared to
normal BHNS mergers where one of the BHs is formed directly at core-collapse.

  5. 5.

On comparing Fig. 4 and 5 we find that most NSNS and BHNS rates fail to fall
within the allowed detectability range. On the other hand, a better agreement
is seen once their intrinsic merger rates are compared with the expected local
intrinsic rates (Fig. 5). Also, we find that there are some variations in the
relative magnitudes of the intrinsic rates in comparison to the detectable
rates for the same underlying remnant mass and SNe kick scheme. At present, we
think that more detections (especially of NSNS and BHNS mergers) need to be
performed before any substantive conclusion can be drawn from these trends.

###  3.2 Implications based on the chirp mass distribution

Figure 3: Contour plots of the expected distributions of detected events in O3
in Mchirp and mass ratio space. The upper panel shows the Hobbs kick and the
Standard remnant mass scheme, the lower panel shows the Bray kick and the
FryerRapid remnant mass scheme. Darker contours represent a higher
probability. The black dotted line distinguishes the regions where mergers
occur with component masses both either below 3M⊙, one above and one below or
both above 3M⊙ (i.e, the expected regions of NSNS, BHNS, and BHBH mergers).
The dark blue symbols are the observed events from O1, the blue symbols are
from O2, the cyan symbols from O3a and the yellow triangles from O3b.

Using rates alone is a limited approach to constrain the models. We must
consider the number of events expected in chirp mass and mass ratio space as
well. To this end, we first consider the predicted chirp mass distribution of
our model GW transient populations. Fig. 2 shows the cumulative number of GW
transients versus chirp mass from the high-mass end of the chirp mass
distribution. We see that most of the model predictions reasonably match the
observations at the high-mass end down to log 1.4 (i.e. $\sim$ 25 M⊙),
although the Standard scheme (both kicks) tends to predict a higher count of
massive mergers than observed. This can be seen more clearly in Fig. 7 which
provides the differential numbers for the same. Below $\sim$ 25 M⊙, the
predicted distributions diverge significantly, with most schemes predicting
more GW transients of lower chirp mass than observed. In the mass range of log
1.4 - log 0.5 ($\sim 25-5M_{\odot}$) Hobbs kick produces more observable
mergers but at the lowest chirp masses (spanning BHNS/NSNS mergers),
populations arising from Bray kick have an excess of systems. Cumulatively,
the Hobbs kick tends to produce more detectable systems. However, among the
schemes presented here, Fig. 2 and 7 makes it possible to identify those that
offer the closest match for the observed cumulative and differential events
count over the Mchirp range.

###  3.3 Implications based on chirp mass - mass ratio distribution

We now further the analysis by comparing the GW transient distributions in
chirp mass and mass ratio ($M_{\rm chirp}$ vs $q$) space. Fig. 3, shows the
schemes that offer the closest fit to the total event rates in Table 1 (i.e.,
Standard \+ Hobbs kick and FryerRapid \+ Bray kick). While the Standard scheme
spans nearly all combinations of remnant mass, however, in the FryerRapid
scheme the mass gap (Bailyn et al., 1998; Özel et al., 2010) becomes apparent
due to a clear lack of BH systems near the BHNS/BHBH dividing lines.

For log(Mchirp/M⊙) in the neighborhood of log 0.7 (i.e. $\sim$ 5 M⊙) - close
to the mass-gap end, - the Standard scheme overpredicts the number of expected
mergers over most parts of the $q$ interval. The FryerRapid scheme in Fig. 3
does achieve a reasonably accurate peak near $M_{\rm chirp}\sim 1.5$ and
$q\sim 0.7$. However, it predicts a second, much larger peak for $M_{\rm
chirp}$ and $q\sim 0.8$ and a stronger BHNS peak. In Fig. 3 the peaks in the
predicted and observed number of events for larger $M_{\rm chirp}$ lie at
similar but not identical $M_{\rm chirp}$ vs $q$ ranges. This suggests a need
for further investigation at the upper end of the BH mass distribution, close
to the mass range for pair-instability SNe (e.g. Woosley 2017; Stevenson et
al. 2019). Additionally, in Fig. 8, the Standard \+ Bray kick scheme
underpredicts the BHBH and over-predicts the NSNS events, consistent with the
findings in fig 2, 4 and 7. The FryerRapid \+ Hobbs kick scheme in Fig. 8 also
fails to provide a reliable match for the total events over the $M_{\rm
chirp}$ vs $q$ range. The $M_{\rm CO,\;final}$ scheme in Fig. 9 produces too
many mergers nearly over the entire $q$ range at the lower and intermediate
Mchirp values. In our simulation, these mergers could be avoided if the
remnant mass is consistently less than $M_{\rm CO,\;final}$ in this $M_{\rm
chirp}$ range which leads to a stronger kick and wider orbits. Overall, on
comparing the distribution of the observed GW transients to our forward-model
predictions, we find that the small number statistics still limit a meaningful
comparison.

##  4 Discussion & Conclusions

Here, we have presented theoretical predictions of the expected number of
events (and their distribution in chirp mass and mass ratio space) detected
within bpass under the LVC O3 sensitivity criteria. We find that none of the
remnant mass schemes and SN kick combinations presented fully describe the
observed population of GW transients in the O3 catalogs. Although most
population models do reasonably well in reproducing the chirp mass
distribution at masses above $\sim$ 25 M⊙ but diverge considerably at lower
masses. Our results suggest that the Standard \+ Hobbs kick scheme performs
the best overall but requires modification at the upper and lower ends of the
chirp mass distribution. The values of the most massive chirp masses predicted
by this scheme need to be reduced by 0.1-0.2 dex, while the mass-gap needs to
be included in some form at the lower end of the BH mass distribution.
However, since the event GW190814 (Abbott et al., 2020b) is not predicted in
the FryerRapid scheme, this suggests that a prescription making low mass BHs
significantly rarer, rather than preventing them from forming entirely within
the mass-gap, may be required (see also, Zevin et al. 2020). The Standard \+
Hobbs kick somewhat underpredicts the observed BHBH events in Table 1, (though
it satisfies the intrinsic merger rate constraint) but these are the outcome
of isolated binary evolution only (e.g. see, Zevin et al. 2021). We note that
here we only vary the SNe remnant mass and SNe kick schemes while keeping the
other uncertain parameters fixed (see, Eldridge et al. 2017). Hence the
predicted values have an implicit uncertainty built into them. Moreover, our
models do not yet allow for tidally-induced chemical homogeneous evolution
which (see de Mink et al. 2009; de Mink & Mandel 2016; Marchant et al. 2016)
could result in a higher BHBH merger rate. Additionally, BHBH mergers could
also come from dynamical channels in young stellar or globular clusters (e.g.,
Kulkarni et al. 1993; Portegies Zwart & McMillan 2000) or at the disks of
active galactic nuclei (e.g., McKernan et al. 2012) and some BHBH systems
could be of primordial origin, arising in the early Universe (Carr & Hawking,
1974). Tang et al. (2020) found that the GW transient rates involving BHs are
highly dependent on the rate of chemical enrichment of the Universe and the
number of BHNS and BHBH mergers can be reduced by altering the cosmic
metallicity evolution. However, rapid enrichment would be required for our
highest predicted rates to match those observed in O3. Lastly, all the
predictions made here, assume fixed values for the rest of the uncertain
parameters that go into the binary stellar evolution (see Eldridge et al. 2017
for more information). Changing these parameters could affect the results of
this paper and hence any reliable prediction can only be made when more
constraints (for the other unknowns) are available from other independent
sources.

To conclude, we find some of our model populations produce a reasonable
agreement between the predictions and observed distributions. Nonetheless,
discrepancies remain, highlighting remaining uncertainties in our
understanding of high mass compact binary formation. Our results indicate that
no one model completely matches the observed distribution, although the
Standard remnant mass scheme + Hobbs kick and FryerRapid scheme + Bray kick
offer the closest match. The key features which are difficult to match for
these model sets are the smaller number of low chirp mass transients and a
significant peak around log 1.48 ($\sim$ 30 M⊙). In the future, we might
consider using a prescription where both the kick and remnant mass are
stochastically sampled, similar to Mandel & Müller (2020). Additionally, at
the high mass end the physics of pair-instability and pulsational pair-
instability systems also need to be taken into account (as in Stevenson et al.
2019). In Fig. 4 and 5, we find variations in the relative magnitudes of the
intrinsic rates in comparison to the detectable rates for the same underlying
remnant mass and SNe kick scheme. Given that supernovae remnant masses and
kicks are interrelated, higher number statistics are needed to fully explore
their effects. This should be possible as GW transient catalogs grow in size
and completeness.

## Acknowledgements

We thank the anonymous referee for their helpful feedback and suggestions. SG
and WGJvZ acknowledge support from the University of Auckland. JJE and HFS
acknowledge the support of the Marsden Fund Council managed through the Royal
Society of New Zealand Te Apārangi. ERS received support from United Kingdom
Science and Technology Facilities Council (STFC) grant ST/T000406/1.

## Data Availability

The data underlying this article will be shared on reasonable request to the
corresponding author.

## References

  * Abbott et al. (2016) Abbott B. P., et al., 2016, Phys. Rev. Lett., 116, 241102
  * Abbott et al. (2017) Abbott B. P., et al., 2017, Phys. Rev. Lett., 119, 161101
  * Abbott et al. (2019) Abbott B., et al., 2019, Physical Review X, 9, 031040
  * Abbott et al. (2020a) Abbott R., et al., 2020a, arXiv e-prints, p. arXiv:2010.14527
  * Abbott et al. (2020b) Abbott R., et al., 2020b, ApJ, 896, L44
  * Allen et al. (2012) Allen B., Anderson W. G., Brady P. R., Brown D. A., Creighton J. D. E., 2012, Phys. Rev. D, 85, 122006
  * Bailyn et al. (1998) Bailyn C. D., Jain R. K., Coppi P., Orosz J. A., 1998, ApJ, 499, 367
  * 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
  * Belczynski et al. (2008) Belczynski K., Kalogera V., Rasio F. A., Taam R. E., Zezas A., Bulik T., Maccarone T. J., Ivanova N., 2008, ApJS, 174, 223
  * Belczynski et al. (2013) Belczynski K., Bulik T., Mandel I., Sathyaprakash B. S., Zdziarski A. A., Mikołajewska J., 2013, ApJ, 764, 96
  * Belczynski et al. (2014) Belczynski K., Buonanno A., Cantiello M., Fryer C. L., Holz D. E., Mandel I., Miller M. C., Walczak M., 2014, ApJ, 789, 120
  * Belczynski et al. (2016) Belczynski K., Holz D. E., Bulik T., O’Shaughnessy R., 2016, Nature, 534, 512
  * Boco et al. (2021) Boco L., Lapi A., Chruslinska M., Donevski D., Sicilia A., Danese L., 2021, ApJ, 907, 110
  * Bray & Eldridge (2018) Bray J. C., Eldridge J. J., 2018, MNRAS, 480, 5657
  * Broekgaarden et al. (2021) Broekgaarden F. S., et al., 2021, arXiv e-prints, p. arXiv:2103.02608
  * Buskirk & Babiuc-Hamilton (2019) Buskirk D., Babiuc-Hamilton M. C., 2019, EJP, 40, 025603
  * Carr & Hawking (1974) Carr B. J., Hawking S. W., 1974, MNRAS, 168, 399
  * Chruslinska et al. (2019) Chruslinska M., Nelemans G., Belczynski K., 2019, MNRAS, 482, 5012
  * De Marco & Izzard (2017) De Marco O., Izzard R. G., 2017, Publ. Astron. Soc. Australia, 34, e001
  * Eldridge & Stanway (2016) Eldridge J. J., Stanway E. R., 2016, MNRAS, 462, 3302
  * Eldridge & Tout (2004) Eldridge J. J., Tout C. A., 2004, MNRAS, 353, 87
  * Eldridge et al. (2017) Eldridge J. J., Stanway E. R., Xiao L., McClelland L. A. S., Taylor G., Ng M., Greis S. M. L., Bray J. C., 2017, PASA, 34, e058
  * Eldridge et al. (2019) Eldridge J. J., Stanway E. R., Tang P. N., 2019, MNRAS, 482, 870
  * Finn (1996) Finn L. S., 1996, Phys. Rev. D, 53, 2878
  * Fishbach & Kalogera (2021) Fishbach M., Kalogera V., 2021, arXiv e-prints, p. arXiv:2105.06491
  * Fryer et al. (2012) Fryer C. L., Belczynski K., Wiktorowicz G., Dominik M., Kalogera V., Holz D. E., 2012, ApJ, 749, 91
  * Giacobbo et al. (2018) Giacobbo N., Mapelli M., Spera M., 2018, MNRAS, 474, 2959
  * Hobbs et al. (2005) Hobbs G., Lorimer D. R., Lyne A. G., Kramer M., 2005, MNRAS, 360, 974
  * Hogg (1999) Hogg D. W., 1999, arXiv e-prints, pp astro–ph/9905116
  * Huerta et al. (2017) Huerta E. A., et al., 2017, Phys. Rev. D, 95, 024038
  * Hurley et al. (2000) Hurley J. R., Pols O. R., Tout C. A., 2000, MNRAS, 315, 543
  * Ivanova et al. (2013) Ivanova N., et al., 2013, A&ARv, 21, 59
  * Kulkarni et al. (1993) Kulkarni S. R., Hut P., McMillan S., 1993, Nature, 364, 421
  * LIGO Scientific Collaboration (2018) LIGO Scientific Collaboration 2018, LIGO Algorithm Library - LALSuite, Free software (GPL), doi:10.7935/GT1W-FZ16
  * Lipunov et al. (2021) Lipunov V., Gorbovskoy E., Grinshpun V., Vlasenko D., 2021, MNRAS, 502, 1925
  * Maeder (1987) Maeder A., 1987, A&A, 178, 159
  * Mandel & Müller (2020) Mandel I., Müller B., 2020, MNRAS, 499, 3214
  * Mapelli et al. (2017) Mapelli M., Giacobbo N., Ripamonti E., Spera M., 2017, MNRAS, 472, 2422
  * Marchant et al. (2016) Marchant P., Langer N., Podsiadlowski P., Tauris T. M., Moriya T. J., 2016, A&A, 588, A50
  * McKernan et al. (2012) McKernan B., Ford K. E. S., Lyra W., Perets H. B., 2012, MNRAS, 425, 460
  * Özel et al. (2010) Özel F., Psaltis D., Narayan R., McClintock J. E., 2010, ApJ, 725, 1918
  * Patton & Sukhbold (2020) Patton R. A., Sukhbold T., 2020, MNRAS, 499, 2803
  * Peters (1964) Peters P. C., 1964, Physical Review, 136, 1224
  * Portegies Zwart & McMillan (2000) Portegies Zwart S. F., McMillan S. L. W., 2000, ApJ, 528, L17
  * Stanway & Eldridge (2018) Stanway E. R., Eldridge J. J., 2018, MNRAS, 479, 75
  * 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
  * Stevenson et al. (2019) Stevenson S., Sampson M., Powell J., Vigna-Gómez A., Neijssel C. J., Szécsi D., Mandel I., 2019, ApJ, 882, 121
  * Tang et al. (2020) Tang P. N., Eldridge J. J., Stanway E. R., Bray J. C., 2020, MNRAS, 493, L6
  * The LIGO Scientific Collaboration et al. (2021a) The LIGO Scientific Collaboration the Virgo Collaboration the KAGRA Collaborationn 2021a, doi:10.5281/zenodo.5655785, https://doi.org/10.5281/zenodo.5655785
  * The LIGO Scientific Collaboration et al. (2021b) The LIGO Scientific Collaboration et al., 2021b, arXiv e-prints, p. arXiv:2111.03606
  * Vigna-Gómez et al. (2018) Vigna-Gómez A., et al., 2018, MNRAS, 481, 4009
  * Woosley (2017) Woosley S. E., 2017, ApJ, 836, 244
  * Woosley et al. (2007) Woosley S. E., Blinnikov S., Heger A., 2007, Nature, 450, 390
  * Zevin et al. (2020) Zevin M., Spera M., Berry C. P. L., Kalogera V., 2020, ApJ, 899, L1
  * Zevin et al. (2021) Zevin M., et al., 2021, ApJ, 910, 152
  * de Mink & Mandel (2016) de Mink S. E., Mandel I., 2016, MNRAS, 460, 3545
  * de Mink et al. (2009) de Mink S. E., Cantiello M., Langer N., Pols O. R., Brott I., Yoon S. C., 2009, A&A, 497, 243
  * van Zeist et al. (2021) van Zeist W. G. J., Stevance H. F., Eldridge J. J., 2021, The Journal of Open Source Software, 6, 2968

xxxx

##  Appendix A Supplementary Figures for GW event distribution

This appendix presents a visual plot (Fig. 4, 5) of the data shown in Table 1,
the local merger rate density distribution as a function of the primary BH’s
mass (Fig. 6) for the various schemes, the differential rates (Fig. 7) for
data shown in Fig. 2 and also the chirp mass - mass ratio distributions
associated with the additional supernova kick and remnant mass prescriptions
as discussed in the text (Fig. 8, 9).

Figure 4: To assist in reading, here we plot the predicted number of
detectable mergers in O3 as shown in Table 1 by the various schemes. Pink band
represents the range of allowed values for BHBH, cyan for BHNS and grey for
NSNS merger rates as obtained by LVC O3, (Abbott et al., 2020a; The LIGO
Scientific Collaboration et al., 2021b). The red crosses are for NSNS, blue
for BHNS and green for BHBH merger rates. The BHBH merger rate for MCO,final
(both kicks) and FryerRapid (Hobbs kick) are not shown due to their large
values. Also, the BHBH merger rate for AlwaysNS (both kicks) is zero. Figure
5: Here we plot the intrinsic local merger rates in Gpc-3 yr-1 as shown in
Table 1 by the various schemes. Pink band represents the range of allowed
values for BHBH, cyan for BHNS and grey for NSNS merger rates as obtained by
LVC O3, (Abbott et al., 2020a; The LIGO Scientific Collaboration et al.,
2021b) using their assumptions for the underlying remnant mass distribution
(which differs from those explored here). The overlapping cyan and grey region
represents the overlap in the allowed ranges for the NSNS and BHNS mergers
rates. The crosses have the same meaning. The BHBH merger rate for AlwaysNS
(both kicks) is zero. Note the change of scale on y-axis. Figure 6: The local
merger rate density (per unit solar mass) distribution as a function of the
primary BH’s mass in the BBH mergers for the various schemes. The blue curve
is the PowerLaw+Peak posterior population distribution for the same with its
90% credible interval, represented by the shaded extension (the PowerLaw+Peak
curve has been reproduced from The LIGO Scientific Collaboration et al. 2021a
which is based on the work of The LIGO Scientific Collaboration et al. 2021b).
The slope ($\alpha$) of the powerlaw in log-log space is
$\alpha=3.4_{-0.49}^{+0.58}$ supplemented with a Gaussian peak at
$34_{-3.8}^{+2.3}M_{\odot}$. The unshaded histogram is for Hobbs while the
shaded is for Bray kick respectively. The bins are logarithmically spaced and
the minimum birth mass of BH is set to 2.5M⊙ resulting in systems on the LHS
of 3M⊙ as well. The presence of systems in the pair-instability region is due
to subsequent accretion by the BH from the companion star. Moreover, the data
presented here does not yet account for the pulsational pair-instability SNe
(e.g. Woosley et al. 2007) and therefore does not feature a clear peak around
34 M⊙. Figure 7: The predicted number of observations by various schemes as a
function of chirp mass. The red histogram represents the number of observed
events in O3 in their respective mass bin. The filled histogram represent Bray
kick while the empty ones represent Hobbs kick. Note that the y-scale varies
among the plots. Apart from the Standard + Hobbs and FryerRapid + Bray, others
schemes predict large (or tiny) peaks at various Mchirp values that are highly
inconsistent with the observed data.

Figure 8: Details as in Fig. 3 but here the LHS panel shows the Standard
remnant mass prescription and the Bray kick while the RHS shows the FryerRapid
remnant mass prescription and the Hobbs kick

Figure 9: Details as in Fig. 3 but here the remnant mass is estimated from the
M${}_{\mathrm{CO,\;final}}$ remnant mass prescription with the LHS showing the
Hobbs kick and the RHS the Bray kick. The mergers are biased towards the mid-
range of the Mchirp values, providing a poor match to the observed data.