[]supp1a []supp1b

# Exploring compact binary merger host galaxies and environments with zELDA

S. Mandhai1,2, G. P. Lamb1, N. R. Tanvir1, J. Bray3,4, C. J. Nixon1, R. A. J.
Eyles-Ferris1, A. J. Levan5,6, B. P. Gompertz7.
1 School of Physics and Astronomy, University of Leicester, University Road,
Leicester, LE1 7RH, United Kingdom
2 Jodrell Bank Centre for Astrophysics, Department of Physics and Astronomy,
The University of Manchester, Manchester, M13 9PY, United Kingdom
3 Department of Physical Sciences, Robert Hooke Building, Kents Hill, Milton
Keynes MK7 6AA, United Kingdom
4 Department of Physics, The University of Auckland, Private Bag 92019,
Auckland, New Zealand
5 Department of Astrophysics/IMAPP, Radboud University, Nijmegen, Netherlands
6 Department of Physics, University of Warwick, Coventry, CV4 7AL, UK
7 School of Physics and Astronomy, University of Birmingham, Birmingham, B15
2TT, UK Contact e-mail<EMAIL_ADDRESS><EMAIL_ADDRESS>e-mail<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

Compact binaries such as double neutron stars or a neutron star paired with a
black-hole, are strong sources of gravitational waves during coalescence and
also the likely progenitors of various electromagnetic phenomena, notably
short-duration gamma-ray bursts (SGRBs), and kilonovae. In this work, we
generate populations of synthetic binaries and place them in galaxies from the
large-scale hydrodynamical galaxy evolution simulation EAGLE. With our zELDA
code, binaries are seeded in proportion to star formation rate, and we follow
their evolution to merger using both the BPASS and COSMIC binary stellar
evolution codes. We track their dynamical evolution within their host galaxy
potential, to estimate the galactocentric distance at the time of the merger.
Finally, we apply observational selection criteria to allow comparison of this
model population with the legacy sample of SGRBs. We find a reasonable
agreement with the redshift distribution (peaking at $0.5<z<1$), host
morphologies and projected galactocentric offsets (modal impact parameter
$\lesssim 10$ kpc). Depending on the binary simulation used, we predict $\sim
16-35\%$ of SGRB events would appear “host-less", i.e. sources that merge with
high impact parameters or have hosts fainter than the detection limit
($H>26$).

###### keywords:

neutron star mergers, black hole - neutron star mergers, binaries: close,
galaxies: general, gamma-ray burst: general

††pubyear: 2021††pagerange: Exploring compact binary merger host galaxies and
environments with zELDA–LABEL:lastpage

##  1 Introduction

The orbits of compact binaries (containing neutron stars or stellar mass black
holes) decay via gravitational wave emission, with inspiral times ranging from
a few Myr to many Hubble times. The final mergers have become of increasing
interest following the first gravitational wave (GW) detections (Abbott et
al., 2016). Where one or both of the components are a neutron star, such a
coalescence can also produce detectable electromagnetic counterparts, opening
the route to multi-messenger astrophysics (Branchesi, 2016; Abbott et al.,
2017b; Murase & Bartos, 2019; Huerta et al., 2019; Mészáros et al., 2019).
This was exemplified by the historic discovery of GW170817, and the
accompanying emission from both a kilonova and relativistic jetted material
(Abbott et al., 2017c; Tanaka et al., 2017; Levan et al., 2017; Im et al.,
2017; Tanvir et al., 2017; Evans et al., 2017; Pian et al., 2017; Alexander et
al., 2017; Arcavi et al., 2017; Chornock et al., 2017; Coulter et al., 2017;
Cowperthwaite et al., 2017; Haggard et al., 2017; Hallinan et al., 2017;
Kasliwal et al., 2017; McCully et al., 2017; Margutti et al., 2017; Nicholl et
al., 2017; Shappee et al., 2017; Smartt et al., 2017; Soares-Santos et al.,
2017; Troja et al., 2017; Villar et al., 2017; Corsi et al., 2018; Lyman et
al., 2018; Resmi et al., 2018; Mooley et al., 2018; Nynka et al., 2018; Hajela
et al., 2019; Piro et al., 2019; Troja et al., 2019a; Lamb et al., 2019a;
Troja et al., 2020).

The kilonova is powered by the radioactive decay of newly synthesised heavy
elements from ejected neutron star material, and may be the dominant site of
r-process nucleosynthesis in the universe. The kilonova accompanying GW170817
(AT2017gfo) was extensively studied from ultraviolet to infrared wavelengths,
thanks to its relative proximity at $\approx 40$ Mpc (Abbott et al., 2017c).
However, more typically, the current generation of GW detectors 111Third
Observing Run (O3) 222https://gracedb.ligo.org/superevents/public/O3/ finds
neutron-star containing mergers at $\sim 200$ Mpc (Abbott et al., 2018), which
combined with the large (10s or 100s of sq-deg) error regions, presents a
significant challenge to counterpart searches (Mandhai et al., 2018; Dichiara
et al., 2020).

The production of relativistic jets during mergers is less well understood,
but is believed to be responsible for the short-duration gamma-ray burst
(SGRB) emission; these events have typical durations, $T_{90}<2\,\rm{s}$
(Mazets & Golenetskii, 1981; Kouveliotou et al., 1993). This emission is much
brighter and can be seen at cosmological distances, indeed, apart from GRB
170817A accompanying GW170817, the nearest SGRB with a secure redshift is GRB
080905A at $z=0.1218$ (Rowlinson et al., 2010) whilst the most distant SGRBs
have been found at $z\gtrsim 2$ (e.g. Levan et al., 2006; Fong et al., 2017;
Selsing et al., 2018). The association of SGRBs with compact binary mergers
confirmed both by analysis of the non-thermal emission from AT2017gfo (Perego
et al., 2017; Wang et al., 2017), and the evidence of kilonova emission
following other SGRBs (e.g. Tanvir et al., 2013; Jin et al., 2016; Gompertz et
al., 2018; Lamb et al., 2019b; Troja et al., 2019b; Jin et al., 2020).

SGRBs were first localised to arcsec precision in the Neil Gehrels Swift
Observatory (hereafter Swift) era, which resulted in the identifications of
host galaxies and hence redshifts (e.g. Gehrels et al., 2005; Bloom et al.,
2006; Fox et al., 2005; Hjorth et al., 2005). The typical redshifts for SGRBs
at the the depth reached by Swift/BAT is $z<1$ (Lien et al., 2016). These
events originate from a wide range of host galaxies and environments,
including some with little or no ongoing star formation (e.g. Fong et al.,
2013; Fong & Berger, 2013; Berger, 2014). The Swift rate of SGRB discovery is
$\sim 10$ yr-1, in contrast, detectors with larger fields-of-view such as the
Fermi/GBM (Gamma-ray Burst Monitor) have a rate of $\sim 40$ yr-1, although
the Fermi localisations are much worse (typically a $\sim\rm{few}\times
10\,\rm{deg^{2}}$) , so redshifts are rarely obtained. Since the afterglows of
SGRBs are also typically fainter than those from long-GRBs
($T_{90}>2\,\rm{s}$) , the net result is that the total sample of SGRBs with
better constrained ($\sim$arcsec–arcmin) sky-positions is $\sim 140$, and only
25–30% of these have confident redshifts (e.g. Fong et al., 2015a; Lien et
al., 2016). It is worth cautioning that the assignation of a given GRB to the
short-hard class is also uncertain in some cases (Bromberg et al., 2013), and
also that a small proportion of the candidate host galaxies may be chance
alignments rather than the real host (Levan et al., 2007). The existing sample
should be understood to be incomplete and inhomogeneous, but still sufficient
to show general trends.

Based on this observed population of SGRBs, a non-negligible fraction ($\sim
15\%$; Fong et al., 2013) appear to have no obvious host galaxy. These events
are considered “host-less" and are typically seen at distances $30-75\,$kpc in
projection from the nearest host galaxy candidate (Berger, 2010; Tunnicliffe
et al., 2014; Fong & Berger, 2013). Studying the orbital evolution of the
progenitor binary pairs relative to their host galaxy, before they merge, may
allow us to investigate the frequency of progenitors migrating to distances
such that they are classed host-less. Establishing constraints on the expected
projected offsets of SGRBs, and the fraction of host-less events, may also aid
in the efforts to identify the original host of the progenitor binary where it
is uncertain, and by extension an improved understanding of the distribution
of these compact binary mergers over the cosmological history. This
combination of properties supports the compact binary merger hypothesis, in
which the binary systems can have long lifetimes before merger, and can
receive large kicks during the supernovae in which the compact remnants were
formed (Hobbs et al., 2005; Nakar, 2007).

In order to understand the relationship between the offsets of detected SGRBs
and their host galaxies, we created a population of synthetic binaries and
evolved them within hydrodynamically simulated galaxies. The compact binaries
of interest discussed in this paper are primarily NSNS and BHNS as these are
the most likely progenitors of SGRBs (Narayan et al., 1992, 2001; Katz &
Canel, 1996; Lee & Ramirez-Ruiz, 2007). Previous studies have also explored
the relationship between galaxies and the offsets of GRBs within static
potentials (e.g. Bloom et al., 1999; Bulik et al., 1999; Bloom et al., 2002;
Belczynski et al., 2006; Church et al., 2011). In this study, we use updated
stellar evolution prescriptions and evolving galaxy potentials to contrast
with and predict observational SGRB offsets. The search for host galaxies and
EM counterparts to GW events can be challenging as seen with the follow up of
GW190814 (Ackley et al., 2020; Morgan et al., 2020; Andreoni et al., 2020).
Studies such as the aforementioned and the one detailed in this paper are
beneficial in improving the efficiency of follow-up campaigns by enabling
prediction of the most likely host populations and offsets. We also comment on
the merger rates predicted by our analysis, scaled by LIGO/VIRGO gravitational
wave constraints for NSNS systems.

In this work, Section 2 introduces our code, zELDA, that is used for the
majority of the data analysis; in Section 3, we outline the results; in
Section 4, we compare our results with observations of SGRBs; in Section 5, we
discuss the relationship of the host galaxies and SGRB localisations; and
finally, in Section 6, we summarise our key findings.

##  2 zELDA: Redshift Electromagnetic Localisation & Deduction Algorithm

Figure 1: A schematic breakdown of the Python based zELDA coding suite used
within this study. Package dependencies include: NumPy (Oliphant, 2006),
Matplotlib (Hunter, 2007), SciPy (Virtanen et al., 2020), CosmoloPy, Pandas
(McKinney, 2010), AstroPy (Astropy Collaboration et al., 2013; Price-Whelan et
al., 2018), and GalPy (Bovy, 2015).

The Redshift Electromagnetic Localisation & Deduction Algorithm
(zELDA)333Repository: https://github.com/hirizon/zELDA-COMBIN is a collection
of scripts designed to process and evolve a population of compact binaries
within synthetic galaxies. An overview of the zELDA code schematic is shown in
Figure 1 444CosmoloPy: http://roban.github.io/CosmoloPy/. The objective of
zELDA is to determine the orbital and stellar evolution of binaries relative
to their host galaxy environment.

Throughout, we have adopted the same flat ($k=0$), $\Lambda$ Friedmann
cosmology as used for the EAGLE simulation. Specifically the density
parameters for matter and dark energy are $\Omega_{M0}=0.307$ and
$\Omega_{\Lambda 0}=0.693$ respectively, and the value for the Hubble constant
used is $H_{0}\,=67.77\,\text{km\,s${}^{-1}\,$Mpc${}^{-1}$}$ (Planck
Collaboration et al., 2014; Furlong et al., 2015).

###  2.1 Compact Binaries

Figure 2: The normalised (by population) merger time distributions (i.e.
$N_{\rm bin}/N_{\rm tot}$) for the seeded NSNS and eBHNS type binaries
produced using different binary population synthesis configurations/packages.
The top, middle, and bottom panels correspond to binaries with metallicities
of Z $=0.004$, $0.008$, and $0.02$, respectively. The numbers in brackets
correspond to the total number of EM-bright unique models that are used for
the binary seeding. We note an increase in systematic noise resulting from the
relative rarity of eBHNS systems within the binary simulations used.

Figure 3: The final systematic velocity distributions corresponding to the
seeded NSNS [upper] and eBHNS [lower] binaries for each simulation used.

####  2.1.1 Formation channels

The evolution of binary stars with masses $>8\ \textrm{M}_{\odot}$ may result
in the formation of compact binary systems that consist of a neutron star or a
black hole paired with a neutron star, i.e. neutron star-neutron star (NSNS),
or black hole-neutron star (BHNS).

For this study, we have simulated a population of synthetic compact binaries
that are formed from the end-point evolution of regular, isolated, stellar
binary systems. Compact binaries may also form through dynamical channels
(e.g. Rodriguez et al., 2015; Rodriguez et al., 2016; Rodriguez & Loeb, 2018;
Fragione & Kocsis, 2018; Samsing & D’Orazio, 2018; Choksi et al., 2019;
Fernández & Kobayashi, 2019; Andrews & Mandel, 2019), for example in globular
or nuclear clusters, or during galactic mergers (Palmese et al., 2017),
although the merger rate of such binaries has been estimated to be small when
compared to isolated binaries (Belczynski et al., 2018; Ye et al., 2019).

Recently Santoliquido et al. (2020) have argued that compact binaries formed
in young stellar clusters will have a more significant contribution to global
merger rates, increasing the rate of NSNS mergers by $\sim 50$% compared to
isolated binaries alone. However, this conclusion also depends on their
assumption of typically low natal-kicks for neutron stars, which itself seems
inconsistent with the frequent large impact parameters seen for SGRBs.
Therefore, we do not explicitly include dynamically formed systems within our
calculations.

####  2.1.2 Simulated binary samples

Separate simulations were performed using two independent stellar evolution
and population synthesis codes. For the first case, the binaries were produced
by the Binary Population and Spectral Synthesis (Eldridge et al., 2017, BPASS
v2.1) code555Whilst BPASS v2.2 is available (Stanway & Eldridge, 2018), we use
an older version for consistency with the REAPER code. and processed through
the Remnant Ejecta and Progenitor Explosion Relationship (REAPER) population
synthesis code (Bray & Eldridge, 2016). We use two configurations of compact
remnant binaries from BPASS which correspond to the prescriptions highlighted
in Bray & Eldridge (2018). Henceforth, these are referred to as
BPASS/Bray666The demonstrative binary sample used for the figures and
numerical values in the main text. and BPASS/Hobbs, respectively. In short,
both configurations use the detailed BPASS stellar evolution models, however,
the net systematic kick velocities for each compact object in BPASS/Bray are
calculated from a kick prescription based on the mass ratio of the supernova
ejecta mass and the final mass of the compact object remnant while for the
BPASS/Hobbs dataset, the kick was randomly selected from the Maxwell-Boltzmann
distribution outlined in Hobbs et al. (2005). In each case, the kicks are
applied in random directions relative to the orbital plane of the binaries.
The relative number of compact binaries that are non-disrupted is dependent on
the kick prescription used.

The third set of simulated binaries were produced using the COSMIC v3.2.0/BSE
code (Breivik et al., 2020; Hurley et al., 2002). Binaries produced with BPASS
are formed using detailed stellar evolution models and prescriptions, whilst
COSMIC uses analytical scalings and is best suited for rapid-population
synthesis of binaries. Including binaries from COSMIC in our analysis allows
us to more densely populate the parameter space, and reduce simulation noise
by increasing the number of models sampled over, in comparison to the BPASS
binaries used. Use of several codes to describe the stellar and binary
evolution also allows us to compare and contrast results stemming from varied
approaches. The net kick velocities for the COSMIC compact objects are
obtained from a Maxwellian distribution with a velocity dispersion of
$265\,\rm{kms^{-1}}$ which is reduced by the factor of $1-f_{\rm fb}$, where
$f_{\rm fb}$ is the fraction of supernova ejecta falling back onto the remnant
object (Fryer et al., 2012).

The stellar evolution codes provide outputs for a range of discrete
metallicities. We use binaries corresponding to a sub-solar, metal poor
($Z=0.004$) stellar population, a moderately metal abundant ($Z=0.008$)
population, and a roughly solar, metal rich ($Z=0.02$) population. The role of
metallicity in the formation of the binaries is discussed in Section 2.3.1.
Noticeably, the COSMIC binaries exhibit lower merger times than the BPASS
systems likely owing to the smaller separations between the compact objects
within the sample. This can be cumulatively attributed to the differences in
the adopted stellar evolutionary prescriptions between both simulations.
COSMIC binaries initially have elliptical orbits and incorporate fallback on
to the compact objects from ejected supernova material. Conversely for BPASS,
the orbits of the binaries are circular and the supernova ejecta fallback is
not considered. Other differences in simulation input fractions are also
present and are detailed in Stanway & Eldridge (BPASS 2018) and Breivik et al.
(COSMIC 2020). Within the work presented in this paper, BPASS/Bray is used as
a primary fiducial model.

Whilst we expect the majority of NSNS binaries to produce an EM counterpart
such as an SGRB, only BHNS systems with a comparatively small mass ratio are
expected to. The formation of an EM counterpart depends on the ejecta mass
released from the neutron star when merging with its companion black-hole
(Tanaka et al., 2014; Kawaguchi et al., 2016; Hinderer et al., 2019; Foucart
et al., 2019). For more massive black holes the NS can pass the event horizon
before it is tidally disrupted, therefore, stifling the formation of an EM
counterpart. The mass ratio of the black-hole to neutron star,
$Q=M_{BH}/M_{NS}$ is, therefore, crucial in governing the mass remaining
following the coalescence of the binary (Shibata & Taniguchi, 2011; Ruiz et
al., 2020). However, a higher black-hole spin increases the value of $Q$ for
which significant NS ejecta is expected (Barbieri et al., 2019). For this
study, we use a $Q<3$ threshold for the mass ratio (Hayashi et al., 2020). The
notation of “eBHNS" within this work refers to the potentially EM bright sub-
population of systems that satisfy our condition. This applies to $\sim
30-37\%$ and $\lesssim\,1\%$ of the BPASS and COSMIC BHNS binaries,
respectively. Notably, COSMIC imposes a mass-gap between $3-5\,\rm{M_{\odot}}$
which results in a smaller fraction of BHNS binaries within the $Q$ threshold.
The formation of these systems relative to the total BHNS population is rare.
The number of unique eBHNS models produced by BPASS and COSMIC are relatively
smaller than the simulated NSNS binaries. Regardless, we have included these
binaries to account for their contribution towards the formation of SGRBs,
although it does lead to a more sparse sampling of the merger time and kick
velocity distributions for eBHNS systems, as is evident in Figures 2 and 3.

###  2.2 Cosmological simulations

We place the simulated binaries into environments representative of the real
universe. To achieve this, we allocate these systems to host galaxies produced
by a cosmological simulation. We use the publicly available output from the
Evolution and Assembly of GaLaxies and their Environments (EAGLE )
project777http://icc.dur.ac.uk/Eagle/database.php (Schaye et al., 2015; Crain
et al., 2015; McAlpine et al., 2016), specifically the RefL0100N1504 sample;
the box size of $(100\,$Mpc$)^{3}$ gives us a large, cosmologically
representative volume to work with. Galaxies with stellar masses in the range,
$\,\textrm{M}_{*}/\textrm{M}_{\odot}>10^{8}$ , are considered well resolved
and successfully reproduce the observed stellar galaxy mass function to
$\lesssim 0.2\,$dex within this mass range (Furlong et al., 2015). We also
include poorly resolved galaxies with masses close to the baryonic particle
mass used in the simulation, $\sim 10^{6}\,\textrm{M}_{\odot}$ (Crain et al.,
2015), to enable us to gauge the fraction of binaries that are seeded (see
Section 2.3.1) into very small host galaxies, while recognising that
individual galaxy properties are not well determined in these cases.

The simulation is split into 28 non-uniformly spaced epochs spanning from
$z\sim 15$ to $z\sim 0$. The use of EAGLE cosmological simulations allows
binaries to be seeded in galaxies over cosmic time. These galaxies will
undergo mergers and produce typically larger galaxies of varied morphologies
(Qu et al., 2017) before the binaries coalesce. Thus, by tracing the evolution
of the binary systems within the evolving potential of the host, we can
extract the galaxy properties when the compact binary coalesces.

###  2.3 Algorithm functionality

####  2.3.1 Binary seeding in EAGLE galaxies

For each cosmological epoch available in our selected EAGLE sample, we seed
the simulated binaries in galaxies in proportion to their star formation rates
and the averaged global metallicity. As the binaries form in discrete
snapshots, a birth time offset is applied to create a continuous distribution
of binary formation redshifts. Previous studies have used similar techniques
to predict the morphology of binary merger host galaxies (e.g. Mapelli et al.,
2018; Artale et al., 2019; Jiang et al., 2020).

Within zELDA, the binaries are seeded within their host galaxies at a distance
that is weighted by the mass profile of the stellar disc. The evolution of the
binary is traced from the birth of the progenitor massive stars to the point
at which the final compact binary coalesces. The stellar evolutionary lifetime
of each star is typically short ($\sim 10$ Myr), and usually ends with a core-
collapse supernova, leaving behind a compact remnant. Each supernova provides
the system with a natal-kick velocity (Belczyński & Bulik, 1999; Fragos et
al., 2009; Giacobbo & Mapelli, 2018). In some cases, the imparted impulse may
disrupt the system entirely, and therefore, a compact binary is never created.
For surviving compact binaries, the direction of the kick is applied to the
entire system randomly with respect to the galactic plane. For sufficiently
large kicks and long lifetimes, systems may migrate to great distances from
their hosts before merging, consistent with the prevalence of SGRBs that are
either apparently host-less or at large offsets on the sky from their likely
hosts (e.g. Fong et al., 2013; Tunnicliffe et al., 2014; Berger, 2014; Zevin
et al., 2019; Gompertz et al., 2020).

The overall star formation rate density within the EAGLE simulation is noted
to be $0.2-0.5\,$dex lower than the observed rate density (Furlong et al.,
2015). For the purpose of this study, we assume the star formation in
different galaxy types reflects the real universe (Madau & Dickinson, 2014),
despite the offset in the overall normalisation. We stack the ordered EAGLE
galactic star formation rates within a snapshot to obtain the cumulative star
formation rate (CSFR). For each binary, a random value between zero and the
maximum CSFR is then selected, and the corresponding galaxy is taken as the
host. Galaxies with higher star formation rates, are therefore, more likely to
host binaries. We use three sub-samples for the synthesised binaries. These
consist of a low-metallicity ($Z=0.004$), an intermediate-metallicity
($Z=0.008$), and a high-metallicity ($Z=0.02$) population (see Section 2.1.2).
The binaries are seeded into galaxies based on the averaged star forming gas
metallicity, $Z_{\rm gal}$. This quantity is allocated to bins corresponding
to the three binary metallicity populations used. An associated binary from
these populations is then selected and seeded into a binned host galaxy, i.e.,
low-met binaries are seeded in galaxies with $Z_{\rm gal}<0.006$, intermediate
binaries in galaxies with $0.006\leq Z_{\rm gal}\leq 0.014$, and high
metallicity binaries in $Z_{\rm gal}>0.014$. The binaries are drawn according
to a Kroupa IMF (Kroupa, 2001) profile for the primary stellar component.

Figure 2 shows the population density of the seeded binaries as a function of
$\tau_{\rm merge}$, the total system lifetime (i.e. stellar evolution plus in-
spiral lifetimes), for all metallicities. The relative fraction of binaries
merging with a given lifespan varies between the type (i.e. NSNS or eBHNS) and
the metallicity of the binary. The total natal-kick velocity distribution of
the seeded NSNS and eBHNS binaries is shown in Figure 3.

For the different binary evolution codes, the rate of merging compact binaries
as a function of redshift is shown in Figure 4. These rates have been
normalised using the LIGO/Virgo O3 estimate of
$320^{+490}_{-240}\,\rm{Gpc^{-3}\,yr^{-1}}$ for the merger rate density of
NSNS systems (Abbott et al., 2021a). The NSNS merger rate density always
exceeds that of EM-bright eBHNS, but the ratio varies significantly between
our three simulations, as demonstrated by the cyan lines in the figure.

Figure 4: Comoving volumetric rate of NSNS (blue) and eBHNS (orange) compact
binary mergers normalised by the local LIGO/Virgo gravitational wave rate,
$320^{+490}_{-240}\,\rm{Gpc^{-3}\,yr^{-1}}$ for NSNS systems (Abbott et al.,
2021a). The shaded regions correspond to the rates contained within the upper
and lower bound of the LIGO/Virgo estimate used. The cyan lines in the bottom
panels show the relative ratio of NSNS against eBHNS systems. The plotted
points correspond to the local volumetric rate estimates for NSNS, and BHNS
($45^{+75}_{-33}\,\rm{Gpc^{-3}\,yr^{-1}}$) (Abbott et al., 2021b) mergers.

####  2.3.2 Orbital evolution of kicked compact binaries within their host
potentials

For binaries that remain intact, we trace their orbital evolution within their
host galaxies. The morphology and potential of the galaxy evolves according to
its descendant in the following EAGLE snapshots encountered during the
lifetime of the orbit. The gravitational potential is modelled as a static
well until the next snapshot, where it is updated if the host evolves or
merges.

In our study, the global potential is composed of three contributors, a
Navarro-Frenk-White (NFW; Navarro et al., 1996) profile for the dark matter
halo, and two Miyamoto-Nagai (MN) (Miyamoto & Nagai, 1975) potentials for the
stellar and gas bodies. The prescriptions for the galactic potential are
scaled accordingly by the component masses for the EAGLE galaxies. The
potential for the NFW profile is described by Equation 1. The EAGLE simulation
provides particle information corresponding to the galaxies which can be used
to construct the gravitational potentials. However, this would be
significantly more computationally expensive with little practical gain. As
such, we do not expect a significant difference between the two methods when
averaged over the total population of galaxies.

$\Phi(r)=-\frac{4\pi G\rho_{0}R_{s}^{3}}{r}{\rm
ln}\Bigg{(}1+\frac{r}{R_{s}}\Bigg{)},$ (1)

where $G$ is the gravitational constant, $\rho_{0}$ is the central density,
$R_{s}$ the scale radius, and $r$ is the radius from the central axis of the
galaxy. $R_{s}$ is fixed by requiring that the half mass radius for the halo
matches that reported in EAGLE. The MN potential on the other hand, is
constructed with the following general prescription,

$\Phi(r,z)=\frac{-GM}{[r^{2}+(a+(z^{2}+b^{2})^{1/2})^{2}]^{1/2}},$ (2)

where $z$ is the height relative to the galactic plane, $a$ is the scale
length, $b$ is the scale height. Increasing the ratio of $b/a$ dictates
whether the geometry of the galaxy is more disc-like ($b<a$) or more
spheroidal ($b\geq a$).

The MN potential (see Equation 2) for the stellar distributions are created
assuming a non-flat potential. We can approximate the scale length, $a$, by
applying a similar minimisation for the half mass of the component as seen
with the estimates of $R_{s}$. The scale height, $b$, is approximated using
the relation, $b=a[(1-\epsilon)/\epsilon]^{2}$, where $\epsilon$ is the
stellar ellipticity for the EAGLE galaxy. A Milky Way like galaxy has an
approximate scale length ratio of $b/a\sim 0.1$.

The morphological structure can be quantified using the
$\alpha_{m}\,=\,(\epsilon^{2}+1-T)/2$ diagnostic described in Thob et al.
(2019), where $T$ is the triaxiality of the galaxy. We broadly distinguish
early-/late-type morphologies based on a fiducial value of $\alpha$. Galaxies
with $\alpha_{m}\leq 0.4$, are assumed to be early-type galaxies, with late-
types assumed to have $\alpha_{m}>0.4$. We note that this is an essentially
structural definition, and at higher redshifts other typical characteristics
of early-type galaxies, such as low specific star formation rate, will not
apply. For galaxies where $\alpha_{m}$ can not be determined due to the
absence of the required shape parameters, we use the red/blue sequence
criterion discussed in Correa et al. (2017). This affects $\lesssim 5\%$ of
the early-type galaxy classifications made.

The effects of galaxy mergers on the binary are approximated within zELDA. If
the host merges with a larger galactic dark matter halo, we apply a pseudo-
kick equivalent to the maximum circular velocity of the major halo in a random
direction to approximately account for the velocity offset entailed in
switching to the new origin. Tidal interactions with the central super-massive
black-holes of galaxies can also result in the disruption of binaries (Amaro-
Seoane et al., 2012; Mainetti et al., 2016). Similar effects may also arise
from interactions with nearby gravitating bodies such as satellite galaxies.
However, in the analysis underlined in this paper, we do not consider the
effects of the central black-hole or nearby neighbouring galaxies. As such,
the resulting gravitational intricacies such as these are beyond the scope of
this paper.

The orbital and galaxy dynamic prescriptions that are used in zELDA
incorporate functionality from the GalPy 888http://github.com/jobovy/galpy
Python package, (Bovy, 2015) using the relevant mass, radius, and shape
properties from the EAGLE simulation.

##  3 Host Galaxies of Compact Binary Mergers

Table 1 shows the host galaxy demographics broken down for each binary
simulation. Figures 5, 6, 8–14 and Tables 2, 3 correspond to BPASS/Bray
binaries. The equivalent results for BPASS/Hobbs and COSMIC are found in
Appendices999The appendices are available online on the journal website. A and
B, respectively.

###  3.1 Galaxy host characteristics of coalescing binaries

Figure 5: [BPASS/Bray - NSNS] The star formation rate ($\phi$) against stellar
mass ($\rm{M_{*}}$) evolution of NSNS merger hosts. In each panel, the point
colours represent the binary merger time. The vertical dashed line indicates
the resolved threshold for EAGLE galaxies (see Section 2.2). The hosts of
ejected binaries are represented in grey. Galaxies that are considered
observationally faint (apparent magnitude $H>26$) are depicted as purple
triangles, the indicated hatched region within $1\sigma$ deviation from the
mean $\phi$ and $\rm{M_{*}}$ of this population. The purple crosses indicate
the host properties of localised SGRBs and their $\phi$ and $\rm{M_{*}}$
uncertainties, (Berger, 2014; Leibler & Berger, 2010). The cyan point
corresponds to the galaxy (NGC 4993) that is the known of host of SGRB 170817A
(Levan et al., 2017; Im et al., 2017) \- the EM counterpart of the binary
neutron star GW detection, GW170817. In each panel, the percentage of
predicted late/early merger hosts are indicated. The estimated total
population of host-less binaries (see Section 3.2.2) and faint hosts are also
stated.

Figure 6: [BPASS/Bray - eBHNS] The star formation rate ($\phi$) against
stellar mass ($\rm{M_{*}}$) evolution of eBHNS merger hosts. In each panel,
the point colours represent the binary merger time. The vertical dashed line
indicates the resolved threshold for EAGLE galaxies (see Section 2.2). The
hosts of ejected binaries are represented in grey. Galaxies that are
considered observationally faint (apparent magnitude $H>26$) are depicted as
purple triangles, the indicated hatched region within $1\sigma$ deviation from
the mean $\phi$ and $\rm{M_{*}}$ of this population. The purple crosses
indicate the host properties of localised SGRBs and their $\phi$ and
$\rm{M_{*}}$ uncertainties, (Berger, 2014; Leibler & Berger, 2010). In each
panel, the percentage of predicted late/early merger hosts are indicated. The
estimated total population of host-less binaries (see Section 3.2.2) and faint
hosts are also stated.

The majority of binaries seeded within the EAGLE galaxies are born within
$0.5\lesssim z\lesssim 3$ and merge within $0\lesssim z\lesssim 1$. The
distribution of the merger time varies between different prescriptions for
binary formation. This enables us to gauge the typical age and relative binary
population as a function of redshift (Chruslinska et al., 2019). The
understanding of binaries merging within their host galaxy is also pivotal in
providing insight into the chemical enrichment of the host galaxy (Shen et
al., 2015; van de Voort et al., 2020). Similarly, ejected systems have
relevance for our understanding of the enrichment of the intergalactic medium.

The successful association of a compact binary and its host galaxy can further
yield constraints on the system’s merger timescale, if the star formation
history of the host is known (McCarthy et al., 2020). In Figures 5 (NSNS), and
6 (eBHNS), we show the relationship between the host galaxy mass and SFR at
various redshift ranges for merging binaries. Colour coding shows the binary
lifetimes. which range from $10$s to $1000$s of Myr. At lower redshift
($z\,<\,1$), a large fraction of BPASS/Bray and BPASS/Hobbs merging binaries
have lifespans of $\tau_{\rm merge}>1000\,$Myr. These systems formed at
earlier epochs ($z>1$) within actively star forming galaxies, but ultimately
merge when the host is in a more quiescent phase.

Coalescing binaries within $z<0.4$ have host galaxies with a median stellar
mass between $M_{*}\sim 1.5-3\times 10^{10}\,\textrm{M}_{\odot}$, and a star
formation rate within the range $\phi\sim
1.3-2.1\,\textrm{M}_{\odot}\,\textrm{yr}^{-1}$. These galaxies tend to be
massive, star-forming, and disk dominated (e.g. Calvi et al., 2018). For
comparison, the Milky Way has a stellar mass, $M_{*}\sim 6\times
10^{10}\,\textrm{M}_{\odot}$, and a star formation rate of, $\phi\sim
1.6\,\textrm{M}_{\odot}\,\textrm{yr}^{-1}$ (Licquia & Newman, 2015). The
purple crosses marked in Figures 5 and 6 indicate the properties of galaxies
associated with localised SGRB detections, for the subset of cases for which
these host parameters have been estimated. Generally, the majority of these
galaxies are consistent with the expected hosts of merging compact objects
predicted by zELDA. For the lowest redshift slice, the over-plotted properties
of the host galaxies reside below the peak of the host distributions. We note
that intrinsically very faint hosts may be missed, or unidentified, in follow-
up observations, so their absence in the observed sample may not be
surprising. The properties of the NSNS merger host associated with GW170817,
NGC 4993, falls within a region where very few galaxies are expected to host
coalescing NSNS systems. An event such as this is rare based on the zELDA
predictions.

We distinguish between early and late-type host galaxies at the point of
coalescence based on the morphological criteria described in Section 2.3.2.
Approximately $50-70\%$ of mergers for BPASS compact binaries appear to
originate from late-type galaxies at all redshifts in each simulation. For the
COSMIC sample, the fraction of binary mergers originating from late-type hosts
can range up to $\sim 80\%$. We note that Artale et al. (2020a) performed a
similar simulation based on EAGLE, and concluded that at low redshifts
($z<0.1$) early-type galaxies dominate the host population. However, the
discrepancy between our studies is due to a different approach in
distinguishing galaxy morphologies. Specifically, they define early-type
galaxies to have specific star formation rates, ${\rm
SFR}/M_{*}<10^{-10}\,{\rm yr}^{-1}$. This means that many galaxies that we
classify as late-type disks at low redshifts are considered early-types within
their work based on their star formation rate criterion.

A non-negligible fraction of these binaries are shown to migrate to distances
where they are considered to be host-less (elaborated on in Section 3.2.2);
the host galaxies for these systems are depicted as grey squares in Figures 5
and 6. The figures also indicate galaxies that could be missed with an
observational follow-up, as they are faint (H$>$26; purple triangles). Small
spiral/dwarf galaxies, and more distant intrinsically brighter galaxies, are
likely to fall into this category due to their low apparent luminosity. The
overall galaxy demographics including the contribution of host-less systems is
shown in Table 1.

###  3.2 Projected galactocentric distance at merger

####  3.2.1 The impact parameter

By following the evolution of binaries from their birth to their coalescence,
we are able to observe the overall migration distance during their lifetimes.
In our simulations, these are turned into projected on-sky distances by
choosing random orientations for each host galaxy and calculating the
corresponding impact parameter.

Figure 7:  The population normalised ($N/N_{\rm bin}$) migrated distance
distributions over all redshifts, for the BPASS/Bray, BPASS/Hobbs, and COSMIC
simulated NSNS (blue) and eBHNS (orange) compact binaries. The fainter red
lines correspond to the total population of seeded binaries, inclusive of non-
bright EM systems (refer to Section 2.1.2). Relative to the COSMIC binary
models sampled over, fewer unique BPASS models are used within our simulation.
The resolution of the binned population is therefore limited*.
*Secondary and tertiary peaks are artefacts of noise. The green dashed line corresponds to a distance of $10\,$kpc.

Figure 7 shows the global distribution of the migration distances that the
NSNS and eBHNS mergers exhibit for the various population synthesis codes
used. We find $\sim 60-70\%$ of NSNS and eBHNS mergers occur within $\lesssim
10\,\textrm{kpc}$ (indicated by the vertical dashed line) from the centre of
their host galaxies, for BPASS/Bray binaries. For BPASS/Hobbs binaries, $\sim
50\%$ of binaries are retained within this distance. For COSMIC binaries,
$\sim 70-80\%$ of systems merge within this distance. In extreme cases, we
find $\lesssim 1\%$ of these binaries merging at galactocentric distances
$>1\,\rm{Mpc}$ for BPASS/Bray and COSMIC binaries, and $\lesssim 5\%$ for the
BPASS/Hobbs population.

Figures 8 and 9 show the impact parameter distributions, where the impact
parameter corresponds to the projected on-sky separations, split by the host
galaxy characteristics (stellar mass and star formation rate) and redshift.
Generally, binaries with high velocities have a higher likelihood of
overcoming the gravitational potential of their host galaxy and thus
travelling further before coalescing. For low mass galaxies such as small disk
or dwarf galaxies, the likelihood increases as the gravitational potential is
weaker and therefore, easier to escape. As such, these galaxies are expected
to retain very few binaries. We see this for both NSNS and eBHNS systems in
the top panels of both figures.

For binaries merging in larger galaxies such as large spirals (middle and
bottom panels) or ellipticals (bottom panel), the coalescence is likely to
occur within $\sim 10\,\rm{kpc}$ and therefore, close to or within the host
galaxy. Tables 2 and 3, show the fractions of binaries merging in each panel
per redshift slice, corresponding to the aforementioned figures.

####  3.2.2 Chance alignment of binaries and their host galaxies

A binary is considered host-less if the offset from its host is likely to be
observationally registered as a chance alignment, with a probability, $P_{\rm
chance}>0.05$ (e.g. Fong et al., 2013) and/or if the host is fainter than
$H>26$. The $P_{\rm chance}$ quantity is calculated using the prescription for
chance alignments from Bloom et al. (2002) and defined as

$P_{\rm chance}=1-e^{-\eta},$ (3)

where

$\eta=\pi r_{i}^{2}\sigma(\leq m).$ (4)

Here, $r_{i}$ is an angular distance that is determined based on the offset of
a binary merger relative to the centre of its host, and the size of the
galaxy. Specifically, if a binary merger occurs at a projected radius less
than twice the projected half mass stellar radius ($R_{*}$) of its host
galaxy, it is assumed to be “within" the host and $r_{i}$ is set to
$r_{i}=2R_{*}$. For mergers occurring further out from their host galaxy
centre $r_{i}=\sqrt{R_{0}^{2}+4R_{*}^{2}}$, where $R_{0}$ is the angular
offset between the centre of the host and the site of the merger. $\sigma(<m)$
is the mean surface density of galaxies brighter than the $H$-band magnitude
($m$) of the host galaxy.

Our procedure differs slightly from the original Bloom et al. (2002)
prescription in that we use the projected half-mass stellar radius in place of
the half-light radius; and we use $H$-band rather than $r$-band magnitudes as
early-type galaxies in particular tend to become rapidly fainter in the
$r$-band at $z\gtrsim 1$ due to the redshifting of the 4000 Å break through
the filter passband. $H$-band number counts for field galaxies are taken from
Metcalfe et al. (2006); Frith et al. (2006), and we use a polynomial
interpolation between the values to extract the mean surface density for any
$H$-band magnitude. In the context of our simulations, we use the given
apparent magnitude of the host galaxies reported by EAGLE.

The $H>26$ criterion accounts for hosts that appear faint (particularly those
at higher redshift) and therefore are likely to be missed when following up
SGRBs, resulting in them being classified as host-less by observers.
Recognising the existence of these cases is essential when accounting for
potential observational biases that can arise with SGRB samples.

Table 1: A population breakdown of EM bright compact binary mergers according
to the probable observed host classification, for the binary simulations used.
The bottom segment of the table compares the total simulated “BAT-detectable"
SGRB progenitors against the demographics breakdown of observed SGRBs from
Fong et al. (2013). The square brackets indicate the demographics of systems
expected to exhibit bright afterglows. A further breakdown of the host-less
population is highlighted beneath the divisions for the “NSNS", “eBHNS", and
“BAT-detectable SGRB Progenitors" segments. Faint galaxies are defined as
sources with a magnitude, $H>26$, while binaries are considered remote if
$P_{\rm chance}>0.05$.
* The demographics fractions corresponding to a sub-population of SGRB progenitors merging in afterglow-bright environments (see Section 4.2).

Binary Source | BPASS/Bray ($\%$) | BPASS/Hobbs ($\%$) | COSMIC ($\%$) |  Fong et al. (2013) ($\%$)
---|---|---|---|---
NSNS
Early Type Galaxies | 26.5 | 20.8 | 26.5 | -
Late Type Galaxies | 50.7 | 40.5 | 51.5 | -
Total Obs. Host-less | 22.9 | 38.7 | 22.1 | -
Faint Galaxies | 0.9 | 1.4 | 5.8 | -
Remote Binaries | 14.1 | 25.7 | 9.1 | -
Faint Host + Remote | 7.8 | 11.7 | 7.2 | -
eBHNS
Early Type Galaxies | 29.6 | 25.0 | 19.3 | -
Late Type Galaxies | 38.0 | 28.8 | 63.9 | -
Total Obs. Host-less | 32.4 | 46.2 | 16.8 | -
Faint Galaxies | 1.5 | 2.9 | 1.8 | -
Remote Binaries | 24.5 | 31.7 | 13.1 | -
Faint Host + Remote | 6.4 | 11.6 | 1.8 | -
BAT-detectable SGRB Progenitors [*]
Early Type Galaxies | 30.2 [31.4] | 23.8 [30.0] | 23.6 [25.1] |  $\gtrapprox$17
Late Type Galaxies | 50.6 [63.1] | 41.5 [59.7] | 61.0 [67.6] | $\approx 47$
Total Obs. Host-less | 19.2 [5.6] | 34.8 [10.3] | 15.5 [7.3] |  $\gtrapprox$17
Faint Galaxies | 0.3 [0.3] | 0.6 [0.6] | 2.7 [1.6] | -
Remote Binaries | 15.1 [4.5] | 28.8 [8.6] | 10.1 [4.3] | -
Faint Host + Remote | 3.7 [0.8] | 5.4 [1.2] | 2.7 [1.5] | -

The total fraction of NSNS and eBHNS that will appear to be host-less is shown
in Table 1 for each binary simulation used. A breakdown of the fraction of
binaries that fulfil our host-less criteria is shown for each panel of Figures
8 and 9 within the second sections of Tables 2 and 3.

Figure 8: [BPASS/Bray - NSNS] For varied redshift slices, each panel
corresponds to a breakdown of the impact parameter distributions for different
galaxy masses (row), and star formation rates (column). Each histogram is
composed of the number of binaries with a given impact parameter normalised
over the total summation of merging systems within the respective panel and
redshift. The dashed vertical line corresponds to the galaxy averaged impact
parameter corresponding to the $P_{\rm chance}$ limit used to identify systems
that may be potentially classed as host-less. A numerical breakdown of each
panel is shown in Table 2.  Table 2: [BPASS/Bray - NSNS] A breakdown of the
relative binary mergers occurring within each panel of Figure 8. The top
division refer to the relative fractions of binaries coalescing for a given
redshift slice. The last column indicates the population of the binaries
contained within each slice as a percentage of the cumulative population. The
second section of the table shows the estimated host-less population per panel
[overall fraction corresponding to the redshift slice] based on the criteria
defined in Section 3.2.

z | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | Total Population (%)
---|---|---|---|---|---|---|---|---|---|---
Binary Fraction (%)
0.0<z<0.4 | 4.6 | 0.1 | 0.0 | 9.9 | 33.9 | 0.0 | 5.5 | 38.0 | 8.0 | 10.6
0.4<z<1.0 | 5.4 | 0.5 | 0.0 | 5.5 | 39.8 | 0.9 | 2.2 | 28.8 | 17.0 | 42.4
1.0<z<2.0 | 7.8 | 2.2 | 0.0 | 2.1 | 43.8 | 7.6 | 0.2 | 11.8 | 24.5 | 36.5
2.0<z<3.6 | 11.4 | 6.0 | 0.0 | 0.9 | 41.4 | 22.4 | 0.0 | 2.4 | 15.6 | 10.5
| Host-less Fraction/Panel [Relative to z-slice] (%) | Total Host-less Fraction (%)
0.0<z<0.4 | 86.8 [4.0] | 67.9 [0.1] | 0.0 [0.0] | 34.2 [3.4] | 14.9 [5.0] | 13.3 [0.0] | 14.8 [0.8] | 5.3 [2.0] | 2.8 [0.2] | 15.5
0.4<z<1.0 | 90.5 [4.9] | 65.9 [0.3] | 0.0 [0.0] | 40.8 [2.2] | 18.9 [7.5] | 9.4 [0.1] | 12.1 [0.3] | 6.3 [1.8] | 5.2 [0.9] | 18.0
1.0<z<2.0 | 93.1 [7.2] | 76.8 [1.7] | 0.0 [0.0] | 51.3 [1.1] | 25.9 [11.4] | 7.8 [0.6] | 6.0 [0.0] | 7.3 [0.9] | 6.5 [1.6] | 24.4
2.0<z<3.6 | 96.5 [11.0] | 86.2 [5.1] | 25.0 [0.0] | 67.9 [0.6] | 37.3 [15.4] | 8.5 [1.9] | 0.0 [0.0] | 6.8 [0.2] | 8.7 [1.4] | 35.6

Figure 9: [BPASS/Bray - eBHNS] For varied redshift slices, each panel
corresponds to a breakdown of the impact parameter distributions for different
galaxy masses (row), and star formation rates (column). Each histogram is
composed of the number of binaries with a given impact parameter normalised
over the total summation of merging systems within the respective panel and
redshift. The dashed vertical line corresponds to the galaxy averaged impact
parameter corresponding to the $P_{\rm chance}$ limit used to identify systems
that may be potentially classed as host-less. A numerical breakdown of each
panel is shown in Table 3.  Table 3: [BPASS/Bray - eBHNS] A breakdown of the
relative binary mergers occurring within each panel of Figure 9. The top
division refer to the relative fractions of binaries coalescing for a given
redshift slice. The last column indicates the population of the binaries
contained within each slice as a percentage of the cumulative population. The
second section of the table shows the estimated host-less population per panel
[overall fraction corresponding to the redshift slice] based on the criteria
defined in Section 3.2.

z | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | Total Population (%)
---|---|---|---|---|---|---|---|---|---|---
Binary Fraction (%)
0.0<z<0.4 | 5.5 | 0.2 | 0.0 | 12.3 | 31.4 | 0.0 | 5.3 | 35.2 | 10.1 | 8.5
0.4<z<1.0 | 7.3 | 1.2 | 0.0 | 5.8 | 42.4 | 0.8 | 1.9 | 24.0 | 16.6 | 31.5
1.0<z<2.0 | 5.4 | 2.9 | 0.0 | 1.5 | 54.2 | 8.9 | 0.1 | 6.9 | 20.1 | 44.1
2.0<z<3.6 | 2.5 | 2.7 | 0.0 | 0.4 | 43.1 | 32.6 | 0.0 | 1.7 | 17.0 | 15.9
| Host-less Fraction/Panel [Relative to z-slice] (%) | Total Host-less Fraction (%)
0.0<z<0.4 | 93.9 [5.2] | 56.2 [0.1] | 0.0 [0.0] | 55.9 [6.9] | 31.8 [10.0] | 0.0 [0.0] | 10.8 [0.6] | 8.9 [3.1] | 3.6 [0.4] | 26.2
0.4<z<1.0 | 90.0 [6.6] | 58.9 [0.7] | 0.0 [0.0] | 59.5 [3.5] | 41.4 [17.5] | 18.6 [0.1] | 15.7 [0.3] | 12.1 [2.9] | 8.7 [1.4] | 33.1
1.0<z<2.0 | 87.7 [4.8] | 69.2 [2.0] | 0.0 [0.0] | 69.6 [1.0] | 38.7 [21.0] | 14.5 [1.3] | 17.9 [0.0] | 13.7 [0.9] | 11.7 [2.4] | 33.4
2.0<z<3.6 | 91.5 [2.3] | 68.2 [1.9] | 0.0 [0.0] | 64.5 [0.3] | 29.6 [12.8] | 8.3 [2.7] | 0.0 [0.0] | 8.1 [0.1] | 13.5 [2.3] | 22.3

##  4 Simulating Observational Selection Effects

###  4.1 Luminosity function and selection

We produce a simulated population of observable SGRBs assuming the broken
power-law intrinsic luminosity function for SGRB peak luminosities from Paul
(2018),

$\Phi(L)\propto\begin{cases}\left(\frac{L}{10^{52.18}{\rm
erg/s}}\right)^{-0.5}&L\leq 10^{52.18}{\rm erg/s},\\\
\left(\frac{L}{10^{52.18}{\rm erg/s}}\right)^{-1.86}&L>10^{52.18}{\rm
erg/s}.\end{cases}$ (5)

A broken power-law is consistent with the findings in Wanderman & Piran
(2015); Ghirlanda et al. (2016); alternatively, see Yonetoku et al. (2014);
Sun et al. (2015) for examples of single power-law luminosity functions. As
the distribution of intrinsically low-luminosity SGRBs is poorly constrained,
we fix the cut-off at $10^{49}\,\rm{erg\ s^{-1}}$ (Mogushi et al., 2019). As
such, the luminosities, $L$, are drawn from a distribution within the range
$10^{49}\leq L\leq 10^{55}\,\rm{erg\ s^{-1}}$. The ascribed SGRB luminosity
has no dependence on the binary characteristics, since we lack any theoretical
or observational mapping that would motivate this. Whilst some SGRBs, such as
GRB 170817A, exhibit luminosities below our minimum, these events are expected
to account for a small fraction of the cumulative observed population (Tan &
Yu, 2020). We use a sensitivity of $0.2$ ph s-1 cm-2 to estimate the fraction
of GRBs that would be detectable by Swift/BAT. We apply a Band function (Band
et al., 1993) with an index $\alpha=0.5$ and $\beta=2.25$ for the distribution
below and above $E_{P}$ respectively, where $E_{p}$ is the spectral peak
energy given by the relation for SGRBs found by (Tsutsui et al., 2013).

###  4.2 Afterglow dependence on environment

GRB afterglows are thought to arise when the relativistic jet impacts and
shocks the external ambient medium, producing long-lived emission from radio
to X-ray (Zhang et al., 2006, 2007). Identification of candidate host
galaxies, and hence redshifts and impact parameters, relies on the
$\sim$arcsec localisations that come from X-ray and optical/nIR detections,
rather than the $\sim$arcmin positions that are produced by the Swift/BAT. For
some SGRBs, particularly those with “extended soft emission” (Mészáros & Rees,
1997; Sari et al., 1998), the X-ray localisation may be found independently of
afterglow detection, however, more commonly observations of the afterglow is
necessary. Thus there is the possibility that SGRBs in low density
environments may be selectively lost from the observational sample we consider
here. In this section we outline an approach to introduce this into the
simulations, which may indicate a plausible maximum size of this selection
bias.

The average gas density, $n$, surrounding a galaxy will decrease with
increasing distance from the galaxy core, therefore, the afterglows from
bursts at large separations from their host will typically have fainter
afterglow emission (Perna & Belczynski, 2002). The maximum luminosity
distance, $D_{\rm max}$, for the peak afterglow emissions to be above a given
threshold depends on the ambient density as $D_{\rm max}\propto n^{(p+1)/8}$
assuming slow-cooling and an observed frequency above the characteristic
synchrotron frequency at the peak, where $p$ is the relativistic electron
distribution index in the shock – typically $p\sim 2.4$ (Fong et al., 2015a).
The afterglows of SGRBs can be used to put constraints on the ambient density,
where the earliest afterglow detection can be used to constrain the minimum
density. Using this minimum ambient density and an assumed galaxy gas density
profile, O’Connor et al. (2020) infer a maximum impact parameter distribution
for the population of observed SGRBs.

For binaries capable of producing an SGRB e.g. NSNS, and eBHNS systems, (see
Section 2.1.2), we estimate the fraction that merge in a very low-density
environment, defined by Fong et al. (2015a) as $n<10^{-4}\,\rm{cm^{-3}}$.
Based on the location of the merger with respect to its host galaxy, we
estimate the ambient density by considering the gas mass contained within an
EAGLE aperture. These apertures are composed of spheres with sizes of $1$,
$3$, $5$, $10$, $20$, $30$, $40$, $50$, $70$, and $100$ kpc. We split these
spheres into shells based on the radii of neighbouring apertures and calculate
the spherically averaged gas density. We added an additional hot gas
component, similar to the approach described in O’Connor et al. (2020). This
secondary component is described by the Maller & Bullock gas model (Maller &
Bullock, 2004; Fang et al., 2013),

$\rho^{MB}_{g}(x)=\rho_{v}\bigg{[}1+\frac{3.7}{x}\rm{ln}(1+x)-\frac{3.7}{c_{v}}\rm{ln}(1+c_{v})\bigg{]}^{3/2},$
(6)

where $\rho_{v}$ is the density at the virial radius, $r_{v}$ (Carlberg et
al., 1997), corresponding to an overdensity parameter of $\Delta_{v}=200$; $x$
is the ratio of $r/R_{s}$ ; and $c_{v}$ is the concentration parameter, (e.g.
Bullock et al., 2001; Comerford & Natarajan, 2007). Using this prescription,
we identify the fraction of the SGRB progenitor mergers that occur in
environments with $n<10^{-4}$ cm-3 as afterglow-faint. Mergers occurring in
denser environments are considered to be afterglow-bright.

The total population of the simulated SGRB progenitors merging in environments
with low ambient densities is shown in the first half of Table 4. This is
further split according to whether the merger fulfils our criteria for being
host-less or not. The proportion of mergers occurring in low density media is
therefore somewhat increased over the whole population. Notably, the
proportion of host-less BAT-detectable SGRBs formed by COSMIC and BPASS/Bray
binaries merging within afterglow-faint environments are broadly consistent,
within errors, with the estimated fraction of SGRBs occurring in a
$n<10^{-4}\,{\rm cm}^{-3}$ ambient density medium by O’Connor et al. (2020),
whereas the $\approx 35$% fraction of SGRBs in low density environments
predicted by the BPASS/Hobbs simulation is rather in excess of the
observations.

In Figure 10 we compare the redshift distribution of observed SGRBs against
our predictions for the BAT-detectable mergers with host associations. The
agreement is generally reasonable, although the finding of several SGRBs at
$z\gtrsim 1.7$ is harder to reconcile with the small numbers predicted at
those redshifts. This may suggest that the bright-end slope of our adopted
luminosity function (Eqn. 5) is rather too steep.

Figure 10: [BPASS/Bray] The distribution of the NSNS and eBHNS binary mergers
from our study that are likely to produce SGRBs detectable by Swift/BAT,
excluding host-less systems, refer to Section 3.2.2. The magenta histogram
outlines the total population of BAT-detectable SGRBs with identifiable hosts.
The lime histograms corresponds to the subset of events merging in afterglow-
bright environments, ($n>10^{-4}\,\rm{cm^{-3}}$). The orange points correspond
to the population of Swift observed SGRBs within redshift bin sizes of
$z=0.2$, with a Poisson uncertainty.  Table 4: Summary of the fractions of
SGRB progenitors merging in environments with ambient densities
$<10^{-4}\,\rm{cm^{-3}}$, for each binary simulation used. Systems are
considered host-less if the host galaxy has a magnitude, $H>26$, and/or if the
offset of the binary is such that $P_{\rm chance}>0.05$.

Binary Simulation | All SGRB Progenitors ($\%$) | BAT-detectable ($\%$)
---|---|---
With Host | Host-less | Total | With Host | Host-less | Total
BPASS/Bray | 4.9 | 14.9 | 19.9 | 9.3 | 14.8 | 24.2
BPASS/Hobbs | 2.9 | 26.1 | 29.0 | 7.0 | 28.1 | 35.0
COSMIC | 1.7 | 8.9 | 10.6 | 3.6 | 9.1 | 12.7

##  5 Identifying the host galaxies of SGRBs

###  5.1 Comparing study predictions against observations

Binaries capable of producing SGRBs within the sensitivity threshold of
Swift/BAT are shown in Figure 11. We have overlaid offset and redshift
information for SGRBs noted in literature. There is a relatively high
concentration of zELDA processed binaries at redshift $\sim 0.5$ with
relatively low impact parameters, $\sim 10\,$kpc. The observed localised SGRBs
tend to cluster at low redshift ($z\lesssim 1$) with spatial offsets
$<100\,$kpc from their likely hosts. This is found to be consistent with the
peak concentration predicted by the binary simulations used.

The identification of the parent galaxy associated with an SGRB is often non-
trivial, especially for events with no bright galaxies in or near their
localisation uncertainty region, or in cases where no arcsec level
localisation was obtained in the first place. Campaigns aiming to successfully
isolate and identify the hosts of these events have been ongoing throughout
the Swift era (e.g. Levan et al., 2007; Berger, 2009; Fong et al., 2013;
Berger, 2014). A host that has a low-apparent luminosity may easily be missed
when searching within or around the localisation uncertainty of a given SGRB
detection. Apart from the rare cases where the SGRB is located on or close to
an obvious bright (typically low redshift) galaxy, the first requirement is to
conduct a deep-imaging of the field of sky surrounding the region of the
burst. This is often provided by follow-up imaging of the afterglow following
the SGRB. Ultimately, associations made between the burst and the host galaxy
is probabilistic given that the galaxy could be within the localised field by
chance, see Section 3.2.2.

Fong et al. (2013) analysed the host demographics for 36 well associated SGRBs
(i.e. events with convincing host galaxies backed by afterglow observations),
finding that $\sim 47$% have late-type hosts, $\sim 17$% early-type, $\sim
17$% host-less, and $\sim 19$% are “inconclusive" owing to their uncertain
association with a nearby galaxy. Comparing our predictions (refer to the
final segment of Table 1) with these values we find reasonably consistent
fractions for SGRBs arising from late-type galaxies (e.g. BPASS/Bray: $\sim
51\%$). For the early-type hosts we predict a higher fraction (e.g.
BPASS/Bray: $\sim 30\%$), although given the limited statistics, and the
possibility that some of the “inconclusive" cases could in fact be early-type
hosts, we do not think this is a serious discrepancy. Similarly, the host-less
population fraction of BAT-detectable SGRBs only slightly varies from the Fong
et al. (2013) numbers for the BPASS/Bray and COSMIC simulations, which could
easily be accounted for if some “inconclusive" cases are ultimately considered
host-less, combined with the observational selection bias against localising
afterglow-faint SGRBs since they are also more likely to be host-less (cf.
Table 4). For BPASS/Hobbs, the $\sim 35$% host-less fraction is only
consistent with the observed population if the majority of the “inconclusive”
population from Fong et al. (2013) is considered host-less.

Figure 11: [BPASS/Bray] A density map of the impact parameter and redshift for
binaries that are expected to produce a SGRB within Swift/BAT’s detection
sensitivity range and are not considered host-less (see Section 3.2.2). In
both cases, offsets and redshifts corresponding to SGRB-galaxy from literature
are plotted with their relevant uncertainties, refer to Table 5. The contour
levels shown, proceeding outwards, indicate the $\sim\,25,\,50,\,75,\,90$ and
$98\%$ enclosed synthetic merging binary population. The total binary
population is comprised of both NSNS and eBHNS binaries.

###  5.2 Implications for host searches for GW detected events

Figure 12: The ratio of NSNS (solid) and eBHNS (dashed) mergers ($N$) per unit
luminosity ($L$) in bins of absolute magnitude, for all galaxies in the
simulation volume (restricted to redshift $z<0.05$, and the BPASS/Bray model).
Shown are results for the $g$/$z$/$H$-band: of these, the $g$-band is
generally flatter, indicating that merger rate is, on average, more nearly
directly proportional to luminosity in this band. Note that the luminosity
scale here is arbitrary, so there is no particular significance to the
relative values of $N/L$.

The rich science returns from observations of neutron-star compact binary
mergers in both gravitational waves and electromagnetic radiation has been
exemplified by the case of GW170817 (e.g. Abbott et al., 2017a, b). However,
even with the current three advanced GW detectors, positional error regions
are typically 10s or 100s of sq-deg in size, while the distance range for such
events is $\sim$ few $\times 10^{2}\,$Mpc (Abbott et al., 2018). This is
challenging both in terms of sky area and depth required to locate EM
counterparts, apart from the rare cases of on-axis SGRBs. To aid in mapping
these error regions, particularly with narrow-field optical, near-IR or X-ray
instruments, various strategies have been suggested to prioritise observations
based on the known positions and properties of potential host galaxies (e.g.
Gehrels et al., 2016; Evans et al., 2016; Artale et al., 2020b).

A question arises as to how to weight the galaxies in such schemes, most
simply by making use of their absolute luminosities in some pass-band (since
that information is more widely available than, for example, more detailed
characterisation of their star formation histories). With this in mind, in
Figure 12 we plot the ratio of the number of mergers (restricted to $z<0.05$)
occurring per unit luminosity (averaged over all galaxies in the simulation
box), as a function of host absolute luminosity, in various pass-bands. A
situation in which binary merger rate was exactly proportional to luminosity
would result in completely flat curves in this diagram. As can be seen, this
is most nearly the case for the $g$-band (both NSNS and eBHNS), which may
reflect the balance between an enhanced rate of mergers from relatively young
(and therefore bluer) stellar populations and the overall dominance of
intermediate and older (redder) populations. In contrast, the (near-IR)
$H$-band magnitude does not provide such a good tracer of merger likelihood,
presumably due to large early-type (elliptical/lenticular) galaxies being
bright in that band, but producing comparatively few mergers.

###  5.3 Implications for SGRB host identification

Figure 13: [BPASS/Bray] The distribution of host-less (see Section 3.2.2)
binary merger offsets against redshift, $z$, corresponding to different ranges
of photon flux, $\Phi$. The colour of the points indicate the brightness of
the host galaxy within the $H$-band, with redder points indicating brighter
host galaxies than the bluer points.

Figure 14: [BPASS/Bray] The impact parameter distribution for the population
of SGRB compatible NSNS and eBHNS binaries expected to be observable by
Swift/BAT for i) all systems (lilac dashed), ii) systems with hosts only
(orange dotted). The blue line corresponds to a sample of SGRBs with known
redshifts and offsets (blue) inclusive of non-Swift detected events. Refer to
Table 5 for sources used.  Table 5: Redshift and impact parameter values for
the SGRBs used within this study. The Eiso is calculated from Swift/BAT’s
$15-350\,\rm{keV}$ energy band.
*Calculated using Fermi/GBM’s $10-10000\,\rm{keV}$ energy band.

GRB | z | IP [kpc] | Eiso [$\times 10^{48}$erg] | Refs.
---|---|---|---|---
050509B | 0.225 | 63.40 | 7.42 | [1–5]
050709 | 0.161 | 3.64 | 29.6* | [1–4, 6]
050724 | 0.258 | 2.57${}_{-0.8}^{+0.8}$ | 83.9 | [1–4, 7]
050813 | 0.719 | 35.50 | 219 | [1, 3, 8]
051210 |  $>$1.400 |  $>$24.90 | 2480 | [1, 3]
051221A | 0.550 | 0.76${}_{-0.25}^{+0.25}$ | 2680 | [1–4, 9–10]
060502B | 0.287 | 73.30 | 49.2 | [1–4, 11]
060801 | 1.130 | 17.60 | 2290 | [1–4, 12–13]
061006 | 0.436 | 3.50 | 1070 | [1–3, 14]
061210 | 0.410 | 15.60 | 552 | [1–4, 15–17]
070429B | 0.902 | 4.70 | 460 | [1–4, 15, 18–20]
070707 |  $<$3.6 | - | 178000* | [1, 18, 21]
070714B | 0.923 | 3.10 | 3340 | [1–4, 15, 18, 22]
070724A | 0.457 | 5.46 ${}_{-0.14}^{+0.14}$ | 39.7 | [1–4, 15, 18, 23–24]
070729 | 0.800 ${}_{-0.1}^{+0.1}$ | 3.21 ${}_{-18.80}^{+18.80}$ | 897 | [3]
070809 | 0.473 | 35.40 | 140 | [1–2, 15, 18]
071227 | 0.381 | 15.50${}_{-0.24}^{+0.24}$ | 204 | [1–3, 14–15, 18]
080503 |  $<$4.000 | - | 16200 | [1, 18]
080905A | 0.122 | 18.50 | 16.1 | [1–4, 18, 25]
090305A |  $<$4.100 | - | 36000 | [18]
090510 | 0.903 | 9.40 | 4980 | [1–3, 15, 18, 26]
100117A | 0.915 | 1.32${}_{-0.33}^{+0.33}$ | 1100 | [1–3, 15, 18]
100206A | 0.407 | 21.90${}_{-18.10}^{+18.10}$ | 263 | [1–2, 27–28]
100625A | 0.452 |  $<$19.80 | 525 | [1–2, 27]
100816A | 0.804 | 10.07 | 7730 | [29–30]
101219A | 0.718 | 3.50 | 3500 | [1–2, 27, 31–32]
111117A | 2.211 | 10.6${}_{-1.70}^{+1.70}$ | 19800 | [2, 33–34]
120804A | 1.300 | 2.20${}_{-1.20}^{+1.20}$ | 18200 | [2, 27, 35]
130603B | 0.356 | 5.21${}_{-0.17}^{+0.17}$ | 803 | [2, 18, 27, 36–37]
140903A | 0.351 | 0.5${}_{-0.2}^{+0.2}$ | 90.7 | [2, 38–40]
150101B | 0.134 | 7.35${}_{-0.07}^{+0.07}$ | 5.84 | [2, 41–43]
150424A | 0.300 | 22.50 | 1080 | [2, 44]
160410A | 1.717 | $\sim 0$ | 1590 | [45–46]
160624A | 0.483 | <11.91 | 192 | [2, 47–48]
160821B | 0.162 | 16.40 | 13.9 | [2, 49–51]
170428A | 0.454 | 6.90${}_{-1.70}^{+1.70}$ | 663 | [2, 52]
170817A | 0.001 | 2.13 | 0.031* | [2, 53]
181123B | 1.754 | 5.08${}_{-1.38}^{+1.38}$ | 10400 | [54]
200522A | 0.554 | 0.93${}_{-0.19}^{+0.19}$ | 287 | [48, 55]
200826A | 0.748 | - | 13000* | [56]
201221D | 1.046 | - | 6660* | [57]

References: [1] Tunnicliffe et al. (2014) [2] Gompertz et al. (2020) [3] Kann
et al. (2011) [4] Church et al. (2012) [5] Bloom et al. (2006) [6] Fox et al.
(2005) [7] Prochaska et al. (2005) [8] Foley et al. (2005) [9] Berger &
Soderberg (2005) [10] Soderberg et al. (2006) [11] Bloom et al. (2007) [12]
Cucchiara et al. (2006) [13] Racusin et al. (2006) [14] D’Avanzo et al. (2009)
[15] Leibler & Berger (2010) [16] Ziaeepour et al. (2006) [17] Berger (2006)
[18] Fong & Berger (2013) [19] Perley et al. (2007) [20] Holland et al. (2007)
[21] Piranomonte et al. (2008) [22] Graham et al. (2009) [23] Berger et al.
(2009) [24] Kocevski et al. (2010) [25] Rowlinson et al. (2010) [26] Rau et
al. (2009) [27] Berger (2014) [28] Cenko et al. (2010) [29] Tanvir et al.
(2010) [30] Im et al. (2010) [31] Kuin & Gelbord (2011) [32] Chornock & Berger
(2011) [33] Selsing et al. (2018) [34] Margutti et al. (2012) [35] Berger et
al. (2013) [36] Thone et al. (2013) [37] Levan et al. (2013) [38] Levan et al.
(2014) [39] Cucchiara et al. (2014) [40] Troja et al. (2016) [41] Levan et al.
(2015) [42] Fong et al. (2015b) [43] Fong et al. (2016) [44] Castro-Tirado et
al. (2015) [45] Selsing et al. (2016) [46] Selsing et al. (2019) [47]
Cucchiara & Levan (2016) [48] O’Connor et al. (2021) [49] Lamb et al. (2019b)
[50] Troja et al. (2019b) [51] Acciari et al. (2021) [52] Izzo et al. (2017)
[53] Abbott et al. (2017c) [54] Paterson et al. (2020) [55] Fong et al. (2020)
[56] Rothberg et al. (2020) [57] de Ugarte Postigo et al. (2020).

Table 6: Two sample Kolmogorov-Smirnov test rejection levels for the simulated BAT-detectable SGRB progenitors impact parameters as compared against observed SGRBs. The second half of data includes the KS-test rejection levels for the subset of systems that merge in afterglow-bright environments, refer to Section 4.2.  BAT-Detectable SGRBs | Kolmogorov-Smirnov Test Rejection ($\%)$
---|---
BPASS/Bray | BPASS/Hobbs | COSMIC
All | 18 | 95 | 81
With Host only | 92 | 64 | 98
| Afterglow Bright Events
All | 90 | 42 | 98
With Host only | 95 | 79 | 99

Successfully ejected binaries with long lifespans are likely to merge further
away from their hosts due to prolonged migration. The greater their
galactocentric radius, the lower the likelihood of a successful association
with their parent galaxy, which, as we have seen in previous sections, will
bias SGRB redshift samples, for example. It is therefore interesting to ask,
if a SGRB detected by Swift/BAT is found to be host-less, what can we say
about its likely redshift and offset from its parent galaxy, based solely on
the peak gamma-ray luminosity?

Figure 13 shows the distribution of the resulting offsets with redshift for
host-less merging systems. Each panel corresponds to a range of peak photon
fluxes, $\Phi$, that would be detected by BAT. The colour of each point
corresponds to the $H$-band magnitude of the host galaxy. The majority of
events with high photon fluxes $\Phi\gtrsim 100\,\rm{cm^{-2}s^{-1}}$ occur
between $0.1\lesssim z\lesssim 0.5$. At all fluxes there are some mergers that
are low redshift and have large offsets ($>100^{\prime\prime}$) on the sky
from their hosts, making association very challenging. For SGRBs with low
photon fluxes ($<1$ ph s-1 cm-2; refer to the final panel of Figure 13) close
to the BAT detection threshold, the fraction of these events with faint host
galaxies ($H>26$) is found to be $\sim 23\%$, $\sim 19\%$, and $\sim 38\%$ for
BPASS/Bray, BPASS/Hobbs, and COSMIC simulations respectively.

###  5.4 Impact Parameters of SGRBs

In Figure 14 we compare the impact parameters for an observed sample of $32$
SGRBs with well-constrained offsets and redshifts (see Table 5) against our
simulated BPASS/Bray SGRB progenitor systems that are detectable by Swift/BAT
(for the same comparison with BPASS/Hobbs, and COSMIC binaries see Figures 22,
30). The median of the observed SGRB impact parameters is $\sim 5.5\,$kpc. In
comparison, for all simulated BAT-detectable SGRBs, we find median impact
parameters of $6.1$, $11.6$, and $4.2\,$kpc for BPASS/Bray, BPASS/Hobbs, and
COSMIC respectively. When the host-less population of progenitors are
excluded, these values are revised to $4.7$, $6.2$, and $3.8\,$kpc.

In Table 6, we present the rejection levels for the two-sample
Kolmogorov–Smirnov (KS) tests comparing impact parameters for the observed
SGRB sample and each binary simulation. We conduct this test for the whole
BAT-detectable SGRB samples as well as the sub-samples of afterglow-bright
events (see Section 4.2). The most relevant comparison is with the “with host"
samples: COSMIC simulations perform poorly in these tests, due to the
typically short life-times and consequently low impact parameters that it
produces. Otherwise BPASS/Hobbs matches well, and BPASS/Bray is marginally
consistent with the observed sample. The results do not vary significantly if
restricted to “afterglow-bright" systems. Finally, we note that all
simulations are reasonably consistent with the $<10\,$kpc median SGRB
galactocentric offset estimated by O’Connor et al. (2020).

##  6 Summary

Using our new zELDA code suite along with binary population synthesis codes,
specifically BPASS and COSMIC, in addition to the cosmological simulation,
EAGLE, we have traced the orbital evolution and host properties of a
population of simulated isolated compact binaries from birth to coalescence.
We use three binary simulated populations with varied prescriptions which are
referred to as, BPASS/Bray, BPASS/Hobbs, and COSMIC. The key difference
between the samples is the employed natal-kick prescriptions and merger time
distributions. Most seeded BPASS/Bray binaries merge on relatively long
timescales of 1–$10\,\rm{Gyr}$ with the majority of natal-kick velocities
ranging between $\sim 100-300\,\rm{km\,s^{-1}}$; BPASS/Hobbs systems exhibit a
broader range of lifespans with higher kick velocities, with the modal
velocity at $v\sim 350\,\rm{km\,s^{-1}}$; COSMIC systems merge with
comparatively short lifespans, $<100\,\rm{Myr}$ and natal-kick velocities that
are typically $<100\,\rm{km\,s^{-1}}$. Despite these differences, we find some
conclusions are common to all the models, specifically:

  * •

For NSNS and eBHNS SGRB progenitors at merger, we find $\sim 50-70\%$ of the
hosts to be late-type galaxies, at all redshifts. For COSMIC binaries, the
upper limit on this demographics estimate can extend to $\sim 80\%$.

  * •

A consequence of this is that, for low redshift events $z<0.05$ (i.e. within
the current range of GW detection), we find the majority of mergers to arise
from relatively blue, star forming hosts. Thus, we find the $g$ band to be
best suited for predicting compact binary merger rates in low redshift
galaxies, and hence useful for weighting/prioritising galaxies for
gravitational wave electromagnetic follow-up observations.

  * •

The modal average of the radial offset from the centre of the host galaxy at
the time of merger is $\lesssim 10\,\rm{kpc}$. For each binary population
used, these values range between $\sim 5-10\,\rm{kpc}$ for BPASS/Hobbs systems
and eBHNS mergers from the BPASS/Bray sample. For BPASS/Bray NSNS and COSMIC
systems, this range varies between $\sim 2.5-5\,\rm{kpc}$. Long-lived systems
with $\tau_{\rm merge}\gtrapprox 1\,\rm{Gyr}$ can potentially travel well
beyond $100\,\rm{kpc}$ from their host galaxy provided that they have natal-
kick velocities that are sufficient for them to escape the gravitational
potential of their host galaxy.

To gauge the consistency of the zELDA processed binaries and their hosts we
compare against the population of well-localised short-duration gamma-ray
bursts (SGRBs) (i.e. events with convincing host galaxies and offsets). We
achieve this by simulating Swift/BAT observations of NSNS and eBHNS mergers
from our study, using an assumed luminosity function consistent with
literature estimates.

Observationally, some SGRBs have been found to be “host-less", in the sense
that there is no clear association with a galaxy, presumably due to either
large offsets at merger and/or because the host is very faint. GRBs without a
confident host association will typically have no redshift and impact
parameter estimates, biases of these cases must be considered in our
comparison. To mimic this, we define simulated binaries to be "host-less" if
either the parent galaxy has a greater than 5% probability of being simply a
chance alignment based on the impact parameter (i.e. projected radial offset
from galaxy centre) and apparent magnitude of the galaxy, and/or has a host
galaxy with a $H$-band magnitude of $H>26$ (i.e. requiring deep follow-up
observations, e.g. with HST to detect). Our findings are:

  * •

For both the BPASS/Bray and COSMIC simulations, we find the hosts of the
simulated SGRBs to be $\sim$50–60% late-type, and $\sim$15–20% host-less, in
reasonable agreement with the demographic distribution of observed hosts (Fong
et al., 2013). In contrast, the BPASS/Hobbs sample predicts a much higher
proportion of host-less events ($\sim 35$%), due largely to the typically
higher kick velocities, which is likely in excess of the host-less fraction
observed.

  * •

The impact parameter distribution of the simulated BAT-detectable SGRBs for
BPASS/Bray and BPASS/Hobbs are reasonably consistent with the observed SGRB
offsets. In this comparison, however, the COSMIC simulation under-predicts the
offsets that are observed, which is a consequence of the lower kick velocities
and shorter lifetimes.

It is interesting to ask the question, if a SGRB is found to be host-less, can
our simulations give a prior expectation for its most likely redshift, host
and projected offset?

  * •

For the brightest host-less simulated SGRBs (those with peak fluxes $>10^{2}$
photons cm-2 s-1), the parent galaxies are found to be always at redshifts
$z<1$ and offsets are typically tens of arcsec, but can range up to in excess
of $\sim 10^{3}$ arcsec (as expected, the BPASS/Hobbs simulation typically
results in rather larger offsets).

  * •

At the other extreme, for bursts close to the Swift/BAT detection threshold,
we find that the majority of events (for all simulations) to be $z<1$, but
with a significant tail to higher redshifts. Notably, for these faint bursts,
a sizable fraction of events ($\sim 23\%$, $\sim 19\%$, and $\sim 38\%$ for
BPASS/Bray, BPASS/Hobbs, and COSMIC populations, respectively) are found to be
host-less by virtue of their hosts being faint ($H>26$) rather than having
large impact parameters.

As a final step, we also considered an additional selection effect, namely
that bursts arising in a low ambient density medium (defined here as
$n<10^{-4}\,\rm{cm^{-3}}$) are expected to have on average fainter afterglows,
and so would be selected against in observed samples that rely on detection of
X-ray or optical counterparts for precise localisation.

  * •

We estimate the fraction of BAT-detectable SGRBs arising from such low density
environments. For the BPASS/Bray, BPASS/Hobbs, and COSMIC samples these
fractions are found to be $\sim 24\%$, $\sim 35\%$, and $\sim 13\%$,
respectively.

Overall, we conclude that of the three population synthesis models considered,
BPASS/Bray gives results most closely aligning with the available
observational constraints i.e. host demographics, and impact parameter
distributions. This may be expected as it includes the most sophisticated
binary evolution modelling, although this conclusion reflects most strongly
the differences in kick velocity than merger time distributions.

Our study does not include binaries that are formed through dynamical
processes, for example in globular or open clusters. In comparison to isolated
systems, the merger rate of dynamically formed binaries has generally been
estimated to be significantly lower. As such, isolated EM-bright compact
binaries are expected to be the dominant channel for producing the observed
population of SGRBs. If these systems are found to be more common, for
instance due to their formation in young stellar clusters, then it will be
important to account for their contribution. With further host associations of
SGRBs and/or with gravitational waves produced by NSNS and BHNS compact binary
mergers, we will be able to further refine and improve constraints on the
demographic and offset estimates derived and discussed in this work.
Similarly, with a greater sample of observed SGRBs, predictions for the
physical input parameters such as the total net-kick velocity and merger time
distributions can be tested from the comparison with associated impact
parameters.

## Acknowledgements

S.F.M. is supported by a PhD studentship funded by the College of Science and
Engineering at the University of Leicester; G.P.L. is supported by STFC grant
ST/S000453/1; N.R.T. acknowledge support from ERC Grant 725246 TEDE; C.J.N is
supported by the Science and Technology Facilities Council grant number
ST/M005917/1 and has received funding from the European Union’s Horizon 2020
research and innovation program under the Marie Skłodowska-Curie grant
agreement No 823823 (Dustbusters RISE project); R.A.J.E-F. is supported by an
STFC PhD studentship. Furthermore, this research used the ALICE High
Performance Computing Facility at the University of Leicester. We acknowledge
the Virgo Consortium for making their simulation data available. The EAGLE
simulations were performed using the DiRAC-2 facility at Durham, managed by
the ICC, and the PRACE facility Curie based in France at TGCC, CEA, Bruyéres-
le-Châtel.

We would like to extend our gratitude to the numerous colleagues who have
provided valuable comments and meaningful discussion over the course of the
work encapsulated within this paper. We would like to extend our gratitude to
Rob Crain and Adrian Thob for their useful discussions regarding EAGLE. We
would also like to thank Kim Page for discussions regarding Swift/XRT
observations.

## Data Availability

The cosmological simulation data underlying this article are available at
EAGLE, for the galaxies used. The BPASS/Bray and BPASS/Hobbs datasets are
derived using the BPASS models from BPASS 2.1. The COSMIC binaries used are
produced using the COSMIC python package available at COSMIC Homepage.

## References

  * Abbott et al. (2016) Abbott B. P., et al., 2016, Phys. Rev. Lett., 116, 061102
  * Abbott et al. (2017a) Abbott B. P., et al., 2017a, Nature, 551, 85
  * Abbott et al. (2017b) Abbott B. P., et al., 2017b, ApJ, 848, L12
  * Abbott et al. (2017c) Abbott B. P., et al., 2017c, ApJ, 848, L13
  * Abbott et al. (2018) Abbott B. P., et al., 2018, Living Reviews in Relativity, 21, 3
  * 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
  * Acciari et al. (2021) Acciari V. A., et al., 2021, ApJ, 908, 90
  * Ackley et al. (2020) Ackley K., et al., 2020, arXiv e-prints, p. arXiv:2002.01950
  * Alexander et al. (2017) Alexander K. D., et al., 2017, ApJ, 848, L21
  * Amaro-Seoane et al. (2012) Amaro-Seoane P., Miller M. C., Kennedy G. F., 2012, MNRAS, 425, 2401
  * Andreoni et al. (2020) Andreoni I., et al., 2020, ApJ, 890, 131
  * Andrews & Mandel (2019) Andrews J. J., Mandel I., 2019, ApJ, 880, L8
  * Arcavi et al. (2017) Arcavi I., et al., 2017, ApJ, 848, L33
  * Artale et al. (2019) Artale M. C., Mapelli M., Giacobbo N., Sabha N. B., Spera M., Santoliquido F., Bressan A., 2019, MNRAS, 487, 1675
  * Artale et al. (2020a) Artale M. C., Mapelli M., Bouffanais Y., Giacobbo N., Pasquato M., Spera M., 2020a, MNRAS, 491, 3419
  * Artale et al. (2020b) Artale M. C., Bouffanais Y., Mapelli M., Giacobbo N., Sabha N. B., Santoliquido F., Pasquato M., Spera M., 2020b, MNRAS, 495, 1841
  * Astropy Collaboration et al. (2013) Astropy Collaboration et al., 2013, A&A, 558, A33
  * Band et al. (1993) Band D., et al., 1993, ApJ, 413, 281
  * Barbieri et al. (2019) Barbieri C., Salafia O. S., Perego A., Colpi M., Ghirlanda G., 2019, A&A, 625, A152
  * Belczyński & Bulik (1999) Belczyński K., Bulik T., 1999, A&A, 346, 91
  * Belczynski et al. (2006) Belczynski K., Perna R., Bulik T., Kalogera V., Ivanova N., Lamb D. Q., 2006, ApJ, 648, 1110
  * Belczynski et al. (2018) Belczynski K., et al., 2018, Astronomy & Astrophysics, 615, A91
  * Berger (2006) Berger E., 2006, GRB Coordinates Network, 5965, 1
  * Berger (2009) Berger E., 2009, ApJ, 690, 231
  * Berger (2010) Berger E., 2010, ApJ, 722, 1946
  * Berger (2014) Berger E., 2014, ARA&A, 52, 43
  * Berger & Soderberg (2005) Berger E., Soderberg A. M., 2005, GRB Coordinates Network, 4384, 1
  * Berger et al. (2009) Berger E., Cenko S. B., Fox D. B., Cucchiara A., 2009, ApJ, 704, 877
  * Berger et al. (2013) Berger E., et al., 2013, ApJ, 765, 121
  * Bloom et al. (1999) Bloom J. S., Sigurdsson S., Pols O. R., 1999, MNRAS, 305, 763
  * Bloom et al. (2002) Bloom J. S., Kulkarni S. R., Djorgovski S. G., 2002, AJ, 123, 1111
  * Bloom et al. (2006) Bloom J. S., et al., 2006, ApJ, 638, 354
  * Bloom et al. (2007) Bloom J. S., et al., 2007, ApJ, 654, 878
  * Bovy (2015) Bovy J., 2015, ApJS, 216, 29
  * Branchesi (2016) Branchesi M., 2016, in Journal of Physics Conference Series. p. 022004, doi:10.1088/1742-6596/718/2/022004
  * Bray & Eldridge (2016) Bray J. C., Eldridge J. J., 2016, MNRAS, 461, 3747
  * Bray & Eldridge (2018) Bray J. C., Eldridge J. J., 2018, Monthly Notices of the Royal Astronomical Society, 480, 5657
  * Breivik et al. (2020) Breivik K., et al., 2020, ApJ, 898, 71
  * Bromberg et al. (2013) Bromberg O., Nakar E., Piran T., Sari R., 2013, ApJ, 764, 179
  * Bulik et al. (1999) Bulik T., Belczyński K., Zbijewski W., 1999, MNRAS, 309, 629
  * Bullock et al. (2001) Bullock J. S., Kolatt T. S., Sigad Y., Somerville R. S., Kravtsov A. V., Klypin A. A., Primack J. R., Dekel A., 2001, MNRAS, 321, 559
  * Calvi et al. (2018) Calvi R., Vulcani B., Poggianti B. M., Moretti A., Fritz J., Fasano G., 2018, MNRAS, 481, 3456
  * Carlberg et al. (1997) Carlberg R. G., et al., 1997, ApJ, 485, L13
  * Castro-Tirado et al. (2015) Castro-Tirado A. J., Sanchez-Ramirez R., Lombardi G., Rivero M. A., 2015, GRB Coordinates Network, 17758, 1
  * Cenko et al. (2010) Cenko S. B., Bloom J. S., Perley D. A., Cobb B. E., Morgan A. N., Miller A. A., Modjaz M., James B., 2010, GRB Coordinates Network, 10389, 1
  * Choksi et al. (2019) Choksi N., Volonteri M., Colpi M., Gnedin O. Y., Li H., 2019, ApJ, 873, 100
  * Chornock & Berger (2011) Chornock R., Berger E., 2011, GRB Coordinates Network, 11518, 1
  * Chornock et al. (2017) Chornock R., et al., 2017, ApJ, 848, L19
  * Chruslinska et al. (2019) Chruslinska M., Nelemans G., Belczynski K., 2019, MNRAS, 482, 5012
  * Church et al. (2011) Church R. P., Levan A. J., Davies M. B., Tanvir N., 2011, MNRAS, 413, 2004
  * Church et al. (2012) Church R. P., Levan A. J., Davies M. B., Tanvir N., 2012, Memorie della Societa Astronomica Italiana Supplementi, 21, 104
  * Comerford & Natarajan (2007) Comerford J. M., Natarajan P., 2007, MNRAS, 379, 190
  * Correa et al. (2017) Correa C. A., Schaye J., Clauwens B., Bower R. G., Crain R. A., Schaller M., Theuns T., Thob A. C. R., 2017, MNRAS, 472, L45
  * Corsi et al. (2018) Corsi A., et al., 2018, ApJ, 861, L10
  * Coulter et al. (2017) Coulter D. A., et al., 2017, Science, 358, 1556
  * Cowperthwaite et al. (2017) Cowperthwaite P. S., et al., 2017, ApJ, 848, L17
  * Crain et al. (2015) Crain R. A., et al., 2015, MNRAS, 450, 1937
  * Cucchiara & Levan (2016) Cucchiara A., Levan A. J., 2016, GRB Coordinates Network, 19565, 1
  * Cucchiara et al. (2006) Cucchiara A., Fox D., Berger E., P.A. P., 2006, GRB Coordinates Network, 5470, 1
  * Cucchiara et al. (2014) Cucchiara A., Cenko S. B., Perley D. A., Capone J., Toy V., 2014, GRB Coordinates Network, 16774, 1
  * D’Avanzo et al. (2009) D’Avanzo P., et al., 2009, A&A, 498, 711
  * Dichiara et al. (2020) Dichiara S., Troja E., O’Connor B., Marshall F. E., Beniamini P., Cannizzo J. K., Lien A. Y., Sakamoto T., 2020, MNRAS, 492, 5011
  * 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, Publ. Astron. Soc. Australia, 34, e058
  * Evans et al. (2016) Evans P. A., et al., 2016, MNRAS, 455, 1522
  * Evans et al. (2017) Evans P. A., et al., 2017, Science, 358, 1565
  * Fang et al. (2013) Fang T., Bullock J., Boylan-Kolchin M., 2013, ApJ, 762, 20
  * Fernández & Kobayashi (2019) Fernández J. J., Kobayashi S., 2019, MNRAS, 487, 1200
  * Foley et al. (2005) Foley R. J., Bloom J. S., Chen H. W., 2005, GRB Coordinates Network, 3808, 1
  * Fong & Berger (2013) Fong W., Berger E., 2013, ApJ, 776, 18
  * Fong et al. (2013) Fong W., et al., 2013, ApJ, 769, 56
  * Fong et al. (2015a) Fong W., Berger E., Margutti R., Zauderer B. A., 2015a, ApJ, 815, 102
  * Fong et al. (2015b) Fong W., Berger E., Fox D., Shappee B. J., 2015b, GRB Coordinates Network, 17333, 1
  * Fong et al. (2016) Fong W., et al., 2016, ApJ, 833, 151
  * Fong et al. (2017) Fong W., et al., 2017, ApJ, 848, L23
  * Fong et al. (2020) Fong W., et al., 2020, arXiv e-prints, p. arXiv:2008.08593
  * Foucart et al. (2019) Foucart F., Duez M. D., Kidder L. E., Nissanke S. M., Pfeiffer H. P., Scheel M. A., 2019, Phys. Rev. D, 99, 103025
  * Fox et al. (2005) Fox D. B., et al., 2005, Nature, 437, 845
  * Fragione & Kocsis (2018) Fragione G., Kocsis B., 2018, Phys. Rev. Lett., 121, 161103
  * Fragos et al. (2009) Fragos T., Willems B., Kalogera V., Ivanova N., Rockefeller G., Fryer C. L., Young P. A., 2009, ApJ, 697, 1057
  * Frith et al. (2006) Frith W. J., Metcalfe N., Shanks T., 2006, MNRAS, 371, 1601
  * Fryer et al. (2012) Fryer C. L., Belczynski K., Wiktorowicz G., Dominik M., Kalogera V., Holz D. E., 2012, ApJ, 749, 91
  * Furlong et al. (2015) Furlong M., et al., 2015, MNRAS, 450, 4486
  * Gehrels et al. (2005) Gehrels N., et al., 2005, Nature, 437, 851
  * Gehrels et al. (2016) Gehrels N., Cannizzo J. K., Kanner J., Kasliwal M. M., Nissanke S., Singer L. P., 2016, ApJ, 820, 136
  * Ghirlanda et al. (2016) Ghirlanda G., et al., 2016, A&A, 594, A84
  * Giacobbo & Mapelli (2018) Giacobbo N., Mapelli M., 2018, MNRAS, 480, 2011
  * Gompertz et al. (2018) Gompertz B. P., et al., 2018, ApJ, 860, 62
  * Gompertz et al. (2020) Gompertz B. P., Levan A. J., Tanvir N. R., 2020, ApJ, 895, 58
  * Graham et al. (2009) Graham J. F., et al., 2009, ApJ, 698, 1620
  * Haggard et al. (2017) Haggard D., Nynka M., Ruan J. J., Kalogera V., Cenko S. B., Evans P., Kennea J. A., 2017, ApJ, 848, L25
  * Hajela et al. (2019) Hajela A., et al., 2019, ApJ, 886, L17
  * Hallinan et al. (2017) Hallinan G., et al., 2017, Science, 358, 1579
  * Hayashi et al. (2020) Hayashi K., Kawaguchi K., Kiuchi K., Kyutoku K., Shibata M., 2020, arXiv e-prints, p. arXiv:2010.02563
  * Hinderer et al. (2019) Hinderer T., et al., 2019, Phys. Rev. D, 100, 063021
  * Hjorth et al. (2005) Hjorth J., et al., 2005, ApJ, 630, L117
  * Hobbs et al. (2005) Hobbs G., Lorimer D. R., Lyne A. G., Kramer M., 2005, MNRAS, 360, 974
  * Holland et al. (2007) Holland S. T., de Pasquale M., Markwardt C. B., 2007, GRB Coordinates Network, 7145, 1
  * Huerta et al. (2019) Huerta E. A., et al., 2019, Nature Reviews Physics, 1, 600
  * Hunter (2007) Hunter J. D., 2007, Computing in Science & Engineering, 9, 90
  * Hurley et al. (2002) Hurley J. R., Tout C. A., Pols O. R., 2002, MNRAS, 329, 897
  * Im et al. (2010) Im M., Park W. K., Pak S., Jeong H., Kim E., Kim J., 2010, GRB Coordinates Network, 11108, 1
  * Im et al. (2017) Im M., et al., 2017, ApJ, 849, L16
  * Izzo et al. (2017) Izzo L., Cano Z., de Ugarte Postig A., Kann D. A., Thoene C., Geier S., 2017, GRB Coordinates Network, 21059, 1
  * Jiang et al. (2020) Jiang Z., et al., 2020, MNRAS, 498, 926
  * Jin et al. (2016) Jin Z.-P., et al., 2016, Nature Communications, 7, 12898
  * Jin et al. (2020) Jin Z.-P., Covino S., Liao N.-H., Li X., D’Avanzo P., Fan Y.-Z., Wei D.-M., 2020, Nature Astronomy, 4, 77
  * Kann et al. (2011) Kann D. A., et al., 2011, ApJ, 734, 96
  * Kasliwal et al. (2017) Kasliwal M. M., et al., 2017, Science, 358, 1559
  * Katz & Canel (1996) Katz J. I., Canel L. M., 1996, ApJ, 471, 915
  * Kawaguchi et al. (2016) Kawaguchi K., Kyutoku K., Shibata M., Tanaka M., 2016, ApJ, 825, 52
  * Kocevski et al. (2010) Kocevski D., et al., 2010, MNRAS, 404, 963
  * Kouveliotou et al. (1993) Kouveliotou C., Meegan C. A., Fishman G. J., Bhat N. P., Briggs M. S., Koshut T. M., Paciesas W. S., Pendleton G. N., 1993, ApJ, 413, L101
  * Kroupa (2001) Kroupa P., 2001, MNRAS, 322, 231
  * Kuin & Gelbord (2011) Kuin N. P. M., Gelbord J. M., 2011, GRB Coordinates Network, 11575, 1
  * Lamb et al. (2019a) Lamb G. P., et al., 2019a, ApJ, 870, L15
  * Lamb et al. (2019b) Lamb G. P., et al., 2019b, ApJ, 883, 48
  * Lee & Ramirez-Ruiz (2007) Lee W. H., Ramirez-Ruiz E., 2007, New Journal of Physics, 9, 17
  * Leibler & Berger (2010) Leibler C. N., Berger E., 2010, ApJ, 725, 1202
  * Levan et al. (2006) Levan A. J., et al., 2006, ApJ, 648, L9
  * Levan et al. (2007) Levan A. J., et al., 2007, MNRAS, 378, 1439
  * Levan et al. (2013) Levan A. J., Tanvir N. R., Wiersema K., Hartoog O., Kolle K., Mendez J., Kupfer T., 2013, GRB Coordinates Network, 14742, 1
  * Levan et al. (2014) Levan A. J., Cenko S. B., Cucchiara A., Perley D. A., 2014, GRB Coordinates Network, 16784, 1
  * Levan et al. (2015) Levan A. J., Hjorth J., Wiersema K., Tanvir N. R., 2015, GRB Coordinates Network, 17281, 1
  * Levan et al. (2017) Levan A. J., et al., 2017, ApJ, 848, L28
  * Licquia & Newman (2015) Licquia T. C., Newman J. A., 2015, ApJ, 806, 96
  * Lien et al. (2016) Lien A., et al., 2016, ApJ, 829, 7
  * Lyman et al. (2018) Lyman J. D., et al., 2018, Nature Astronomy, 2, 751
  * Madau & Dickinson (2014) Madau P., Dickinson M., 2014, Annual Review of Astronomy and Astrophysics, 52, 415–486
  * Mainetti et al. (2016) Mainetti D., Lupi A., Campana S., Colpi M., 2016, MNRAS, 457, 2516
  * Maller & Bullock (2004) Maller A. H., Bullock J. S., 2004, MNRAS, 355, 694
  * Mandhai et al. (2018) Mandhai S., Tanvir N., Lamb G., Levan A., Tsang D., 2018, Galaxies, 6, 130
  * Mapelli et al. (2018) Mapelli M., Giacobbo N., Toffano M., Ripamonti E., Bressan A., Spera M., Branchesi M., 2018, MNRAS, 481, 5324
  * Margutti et al. (2012) Margutti R., et al., 2012, ApJ, 756, 63
  * Margutti et al. (2017) Margutti R., et al., 2017, ApJ, 848, L20
  * Mazets & Golenetskii (1981) Mazets E. P., Golenetskii S. V., 1981, Ap&SS, 75, 47
  * McAlpine et al. (2016) McAlpine S., et al., 2016, Astronomy and Computing, 15, 72
  * McCarthy et al. (2020) McCarthy K. S., Zheng Z., Ramirez-Ruiz E., 2020, arXiv e-prints, p. arXiv:2007.15024
  * McCully et al. (2017) McCully C., et al., 2017, ApJ, 848, L32
  * McKinney (2010) McKinney W., 2010, in van der Walt S., Millman J., eds, Proceedings of the 9th Python in Science Conference. pp 51 – 56
  * Mészáros & Rees (1997) Mészáros P., Rees M. J., 1997, ApJ, 476, 232
  * Mészáros et al. (2019) Mészáros P., Fox D. B., Hanna C., Murase K., 2019, Nature Reviews Physics, 1, 585
  * Metcalfe et al. (2006) Metcalfe N., Shanks T., Weilbacher P. M., McCracken H. J., Fong R., Thompson D., 2006, MNRAS, 370, 1257
  * Miyamoto & Nagai (1975) Miyamoto M., Nagai R., 1975, PASJ, 27, 533
  * Mogushi et al. (2019) Mogushi K., Cavaglià M., Siellez K., 2019, ApJ, 880, 55
  * Mooley et al. (2018) Mooley K. P., et al., 2018, Nature, 561, 355
  * Morgan et al. (2020) Morgan R., et al., 2020, ApJ, 901, 83
  * Murase & Bartos (2019) Murase K., Bartos I., 2019, Annual Review of Nuclear and Particle Science, 69, 477
  * Nakar (2007) Nakar E., 2007, Phys. Rep., 442, 166
  * Narayan et al. (1992) Narayan R., Paczynski B., Piran T., 1992, ApJ, 395, L83
  * Narayan et al. (2001) Narayan R., Piran T., Kumar P., 2001, ApJ, 557, 949
  * Navarro et al. (1996) Navarro J. F., Frenk C. S., White S. D. M., 1996, ApJ, 462, 563
  * Nicholl et al. (2017) Nicholl M., et al., 2017, ApJ, 848, L18
  * Nynka et al. (2018) Nynka M., Ruan J. J., Haggard D., Evans P. A., 2018, ApJ, 862, L19
  * O’Connor et al. (2020) O’Connor B., Beniamini P., Kouveliotou C., 2020, MNRAS, 495, 4782
  * O’Connor et al. (2021) O’Connor B., et al., 2021, MNRAS, 502, 1279
  * Oliphant (2006) Oliphant T. E., 2006, A guide to NumPy. Vol. 1, Trelgol Publishing USA
  * Palmese et al. (2017) Palmese A., et al., 2017, ApJ, 849, L34
  * Paterson et al. (2020) Paterson K., et al., 2020, ApJ, 898, L32
  * Paul (2018) Paul D., 2018, MNRAS, 477, 4275
  * Perego et al. (2017) Perego A., Radice D., Bernuzzi S., 2017, ApJ, 850, L37
  * Perley et al. (2007) Perley D. A., Bloom J. S., Modjaz M., Poznanski D., Thoene C. C., 2007, GRB Coordinates Network, 7140, 1
  * Perna & Belczynski (2002) Perna R., Belczynski K., 2002, ApJ, 570, 252
  * Pian et al. (2017) Pian E., et al., 2017, Nature, 551, 67
  * Piranomonte et al. (2008) Piranomonte S., et al., 2008, A&A, 491, 183
  * Piro et al. (2019) Piro L., et al., 2019, MNRAS, 483, 1912
  * Planck Collaboration et al. (2014) Planck Collaboration et al., 2014, A&A, 571, A16
  * Price-Whelan et al. (2018) Price-Whelan A. M., et al., 2018, AJ, 156, 123
  * Prochaska et al. (2005) Prochaska J. X., Bloom J. S., Chen H. W., Hansen B., Kalirai J., Rich M., Richer H., 2005, GRB Coordinates Network, 3700, 1
  * Qu et al. (2017) Qu Y., et al., 2017, MNRAS, 464, 1659
  * Racusin et al. (2006) Racusin J. L., Burrows D. N., Gehrels N., 2006, GRB Coordinates Network, 5773, 1
  * Rau et al. (2009) Rau A., McBreen S., Kruehler T., 2009, GRB Coordinates Network, 9353, 1
  * Resmi et al. (2018) Resmi L., et al., 2018, ApJ, 867, 57
  * Rodriguez & Loeb (2018) Rodriguez C. L., Loeb A., 2018, ApJ, 866, L5
  * Rodriguez et al. (2015) Rodriguez C. L., Morscher M., Pattabiraman B., Chatterjee S., Haster C.-J., Rasio F. A., 2015, Phys. Rev. Lett., 115, 051101
  * Rodriguez et al. (2016) Rodriguez C. L., Chatterjee S., Rasio F. A., 2016, Phys. Rev. D, 93, 084029
  * Rothberg et al. (2020) Rothberg B., Kuhn O., Veillet C., Allanson S., 2020, GRB Coordinates Network, 28319, 1
  * Rowlinson et al. (2010) Rowlinson A., et al., 2010, MNRAS, 408, 383
  * Ruiz et al. (2020) Ruiz M., Paschalidis V., Tsokaros A., Shapiro S. L., 2020, Phys. Rev. D, 102, 124077
  * Samsing & D’Orazio (2018) Samsing J., D’Orazio D. J., 2018, MNRAS, 481, 5445
  * Santoliquido et al. (2020) Santoliquido F., Mapelli M., Bouffanais Y., Giacobbo N., Di Carlo U. N., Rastello S., Artale M. C., Ballone A., 2020, arXiv e-prints, p. arXiv:2004.09533
  * Sari et al. (1998) Sari R., Piran T., Narayan R., 1998, ApJ, 497, L17
  * Schaye et al. (2015) Schaye J., et al., 2015, MNRAS, 446, 521
  * Selsing et al. (2016) Selsing J., et al., 2016, GRB Coordinates Network, 19274, 1
  * Selsing et al. (2018) Selsing J., et al., 2018, A&A, 616, A48
  * Selsing et al. (2019) Selsing J., et al., 2019, A&A, 623, A92
  * Shappee et al. (2017) Shappee B. J., et al., 2017, Science, 358, 1574
  * Shen et al. (2015) Shen S., Cooke R. J., Ramirez-Ruiz E., Madau P., Mayer L., Guedes J., 2015, ApJ, 807, 115
  * Shibata & Taniguchi (2011) Shibata M., Taniguchi K., 2011, Living Reviews in Relativity, 14, 6
  * Smartt et al. (2017) Smartt S. J., et al., 2017, Nature, 551, 75
  * Soares-Santos et al. (2017) Soares-Santos M., et al., 2017, ApJ, 848, L16
  * Soderberg et al. (2006) Soderberg A. M., et al., 2006, ApJ, 650, 261
  * Stanway & Eldridge (2018) Stanway E. R., Eldridge J. J., 2018, MNRAS, 479, 75
  * Sun et al. (2015) Sun H., Zhang B., Li Z., 2015, ApJ, 812, 33
  * Tan & Yu (2020) Tan W.-W., Yu Y.-W., 2020, arXiv e-prints, p. arXiv:2006.02060
  * Tanaka et al. (2014) Tanaka M., Hotokezaka K., Kyutoku K., Wanajo S., Kiuchi K., Sekiguchi Y., Shibata M., 2014, ApJ, 780, 31
  * Tanaka et al. (2017) Tanaka M., et al., 2017, PASJ, 69, 102
  * Tanvir et al. (2010) Tanvir N. R., et al., 2010, GRB Coordinates Network, 11123, 1
  * Tanvir et al. (2013) Tanvir N. R., Levan A. J., Fruchter A. S., Hjorth J., Hounsell R. A., Wiersema K., Tunnicliffe R. L., 2013, Nature, 500, 547
  * Tanvir et al. (2017) Tanvir N. R., et al., 2017, ApJ, 848, L27
  * Thob et al. (2019) Thob A. C. R., et al., 2019, MNRAS, 485, 972
  * Thone et al. (2013) Thone C. C., de Ugarte Postigo A., Gorosabel J., Tanvir N., Fynbo J. P. U., 2013, GRB Coordinates Network, 14744, 1
  * Troja et al. (2016) Troja E., et al., 2016, ApJ, 827, 102
  * Troja et al. (2017) Troja E., et al., 2017, Nature, 551, 71
  * Troja et al. (2019a) Troja E., et al., 2019a, MNRAS, 489, 1919
  * Troja et al. (2019b) Troja E., et al., 2019b, MNRAS, 489, 2104
  * Troja et al. (2020) Troja E., et al., 2020, MNRAS, 498, 5643
  * Tsutsui et al. (2013) Tsutsui R., Yonetoku D., Nakamura T., Takahashi K., Morihara Y., 2013, MNRAS, 431, 1398
  * Tunnicliffe et al. (2014) Tunnicliffe R. L., et al., 2014, MNRAS, 437, 1495
  * Villar et al. (2017) Villar V. A., et al., 2017, ApJ, 851, L21
  * Virtanen et al. (2020) Virtanen P., et al., 2020, Nature Methods, 17, 261
  * Wanderman & Piran (2015) Wanderman D., Piran T., 2015, MNRAS, 448, 3026
  * Wang et al. (2017) Wang H., et al., 2017, ApJ, 851, L18
  * Ye et al. (2019) Ye C. S., Fong W.-f., Kremer K., Rodriguez C. L., Chatterjee S., Fragione G., Rasio F. A., 2019, arXiv e-prints, p. arXiv:1910.10740
  * Yonetoku et al. (2014) Yonetoku D., Nakamura T., Sawano T., Takahashi K., Toyanago A., 2014, ApJ, 789, 65
  * Zevin et al. (2019) Zevin M., Kelley L. Z., Nugent A., Fong W.-f., Berry C. P. L., Kalogera V., 2019, arXiv e-prints, p. arXiv:1910.03598
  * Zhang et al. (2006) Zhang B., Fan Y. Z., Dyks J., Kobayashi S., Mészáros P., Burrows D. N., Nousek J. A., Gehrels N., 2006, ApJ, 642, 354
  * Zhang et al. (2007) Zhang B., et al., 2007, ApJ, 655, 989
  * Ziaeepour et al. (2006) Ziaeepour H., et al., 2006, GRB Coordinates Network, 5948, 1
  * de Ugarte Postigo et al. (2020) de Ugarte Postigo A., Kann D. A., Izzo L., Thoene C. C., Blazek M., Agui Fernandez J. F., Lombardi G., 2020, GRB Coordinates Network, 29132, 1
  * van de Voort et al. (2020) van de Voort F., Pakmor R., Grand R. J. J., Springel V., Gómez F. A., Marinacci F., 2020, MNRAS, 494, 4867

##  Appendix A BPASS/Hobbs Binaries

Here we show the equivalent figures and tables from the main text, for
BPASS/Hobbs simulated binaries.

Figure 15: [BPASS/Hobbs - NSNS] Similar to Figure 5.

Figure 16: [BPASS/Hobbs - eBHNS] Similar to Figure 6.  Figure 17: [BPASS/Hobbs
- NSNS] For varied redshift slices, each panel corresponds to a breakdown of
the impact parameter distributions for different galaxy masses (row), and star
formation rates (column). Each histogram is composed of the number of binaries
with a given impact parameter normalised over the total summation of merging
systems within the respective panel and redshift. The dashed vertical line
corresponds to the galaxy averaged impact parameter corresponding to the
$P_{\rm chance}$ limit used to identify systems that may be potentially
classed as "host-less". A numerical breakdown of each panel is shown in Table
7.  Table 7: [BPASS/Hobbs - NSNS] A breakdown of the relative binary mergers
occurring within each panel of Figure 17. The top division refer to the
relative fractions of binaries coalescing for a given redshift slice. The last
column indicates the population of the binaries contained within each slice as
a percentage of the cumulative population. The second section of the table
shows the estimated host-less population per panel [overall fraction
corresponding to the redshift slice] based on the criteria defined in Section
3.2.

z | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | Total Population (%)
---|---|---|---|---|---|---|---|---|---|---
Binary Fraction (%)
0.0<z<0.4 | 4.9 | 0.3 | 0.0 | 9.0 | 34.8 | 0.0 | 5.2 | 37.4 | 8.3 | 6.9
0.4<z<1.0 | 6.6 | 1.2 | 0.0 | 4.3 | 45.0 | 1.3 | 1.6 | 23.2 | 16.8 | 31.6
1.0<z<2.0 | 7.9 | 3.5 | 0.0 | 1.5 | 47.0 | 10.4 | 0.1 | 7.4 | 22.2 | 43.4
2.0<z<3.6 | 11.6 | 9.5 | 0.0 | 0.6 | 41.1 | 21.9 | 0.0 | 1.7 | 13.5 | 18.1
| Host-less Fraction/Panel [Relative to z-slice] (%) | Total Host-less Fraction (%)
0.0<z<0.4 | 92.2 [4.5] | 67.3 [0.2] | 0.0 [0.0] | 53.6 [4.8] | 36.1 [12.6] | 35.3 [0.0] | 22.8 [1.2] | 12.6 [4.7] | 6.2 [0.5] | 28.6
0.4<z<1.0 | 92.0 [6.1] | 72.0 [0.8] | 50.0 [0.0] | 60.7 [2.6] | 39.7 [17.9] | 16.8 [0.2] | 19.0 [0.3] | 14.9 [3.4] | 12.0 [2.0] | 33.4
1.0<z<2.0 | 93.2 [7.3] | 79.6 [2.8] | 40.0 [0.0] | 68.3 [1.0] | 42.9 [20.2] | 16.1 [1.7] | 18.7 [0.0] | 16.0 [1.2] | 10.5 [2.3] | 36.5
2.0<z<3.6 | 95.1 [11.0] | 87.5 [8.3] | 33.3 [0.0] | 73.5 [0.5] | 56.6 [23.3] | 20.2 [4.4] | 0.0 [0.0] | 12.5 [0.2] | 12.6 [1.7] | 49.4

Figure 18: [BPASS/Hobbs - eBHNS] For varied redshift slices, each panel
corresponds to a breakdown of the impact parameter distributions for different
galaxy masses (row), and star formation rates (column). Each histogram is
composed of the number of binaries with a given impact parameter normalised
over the total summation of merging systems within the respective panel and
redshift. The dashed vertical line corresponds to the galaxy averaged impact
parameter corresponding to the $P_{\rm chance}$ limit used to identify systems
that may be potentially classed as "host-less". A numerical breakdown of each
panel is shown in Table 8.  Table 8: [BPASS/Hobbs - eBHNS] A breakdown of the
relative binary mergers occurring within each panel of Figure 18. The top
division refer to the relative fractions of binaries coalescing for a given
redshift slice. The last column indicates the population of the binaries
contained within each slice as a percentage of the cumulative population. The
second section of the table shows the estimated host-less population per panel
[overall fraction corresponding to the redshift slice] based on the criteria
defined in Section 3.2.

z | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | Total Population (%)
---|---|---|---|---|---|---|---|---|---|---
Binary Fraction (%)
0.0<z<0.4 | 9.9 | 0.7 | 0.0 | 10.2 | 32.0 | 0.1 | 4.9 | 34.3 | 8.0 | 3.4
0.4<z<1.0 | 11.2 | 2.0 | 0.0 | 4.9 | 48.2 | 1.0 | 0.9 | 17.7 | 14.0 | 24.1
1.0<z<2.0 | 7.9 | 4.2 | 0.0 | 1.2 | 54.0 | 8.4 | 0.1 | 5.7 | 18.5 | 48.0
2.0<z<3.6 | 7.1 | 7.1 | 0.0 | 0.3 | 44.6 | 27.1 | 0.0 | 1.1 | 12.6 | 24.5
| Host-less Fraction/Panel [Relative to z-slice] (%) | Total Host-less Fraction (%)
0.0<z<0.4 | 88.5 [8.7] | 45.9 [0.3] | 0.0 [0.0] | 75.8 [7.7] | 57.8 [18.5] | 25.0 [0.0] | 38.5 [1.9] | 28.6 [9.8] | 13.7 [1.1] | 48.1
0.4<z<1.0 | 84.9 [9.5] | 54.8 [1.1] | 0.0 [0.0] | 83.6 [4.1] | 57.2 [27.6] | 36.9 [0.4] | 40.8 [0.4] | 34.2 [6.0] | 21.6 [3.0] | 52.1
1.0<z<2.0 | 84.3 [6.7] | 64.1 [2.7] | 0.0 [0.0] | 81.8 [1.0] | 43.7 [23.6] | 21.1 [1.8] | 46.2 [0.0] | 28.6 [1.6] | 19.5 [3.6] | 41.0
2.0<z<3.6 | 91.6 [6.5] | 78.2 [5.6] | 25.0 [0.0] | 75.4 [0.2] | 44.7 [19.9] | 15.9 [4.3] | 0.0 [0.0] | 13.5 [0.1] | 15.9 [2.0] | 38.7

Figure 19: [BPASS/Hobbs] Similar to Figure 13.

Figure 20: [BPASS/Hobbs] Similar to Figure 11.

Figure 21:  [BPASS/Hobbs] Similar to Figure 10.

Figure 22: [BPASS/Hobbs] Similar to Figure 14. Refer to Table 5 for sources
used.

##  Appendix B COSMIC Binaries

Here we show the equivalent figures and tables from the main text, for COSMIC
simulated binaries.

Figure 23: [COSMIC - NSNS] Similar to Figure 5.

Figure 24: [COSMIC - eBHNS] Similar to Figure 6.  Figure 25: [COSMIC - NSNS]
For varied redshift slices, each panel corresponds to a breakdown of the
impact parameter distributions for different galaxy masses (row), and star
formation rates (column). Each histogram is composed of the number of binaries
with a given impact parameter normalised over the total summation of merging
systems within the respective panel and redshift. The dashed vertical line
corresponds to the galaxy averaged impact parameter corresponding to the
$P_{\rm chance}$ limit used to identify systems that may be potentially
classed as "host-less". A numerical breakdown of each panel is shown in Table
9.  Table 9: [COSMIC - NSNS] A breakdown of the relative binary mergers
occurring within each panel of Figure 25. The top division refer to the
relative fractions of binaries coalescing for a given redshift slice. The last
column indicates the population of the binaries contained within each slice as
a percentage of the cumulative population. The second section of the table
shows the estimated host-less population per panel [overall fraction
corresponding to the redshift slice] based on the criteria defined in Section
3.2.

z | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | Total Population (%)
---|---|---|---|---|---|---|---|---|---|---
Binary Fraction (%)
0.0<z<0.4 | 4.7 | 1.3 | 0.0 | 3.1 | 53.9 | 0.1 | 0.7 | 26.7 | 9.6 | 3.3
0.4<z<1.0 | 4.6 | 2.4 | 0.0 | 1.2 | 55.9 | 2.5 | 0.1 | 14.4 | 19.0 | 26.5
1.0<z<2.0 | 5.2 | 4.8 | 0.0 | 0.4 | 48.1 | 13.8 | 0.0 | 4.3 | 23.4 | 46.4
2.0<z<3.6 | 7.6 | 11.0 | 0.0 | 0.1 | 39.7 | 26.9 | 0.0 | 0.9 | 13.7 | 23.8
| Host-less Fraction/Panel [Relative to z-slice] (%) | Total Host-less Fraction (%)
0.0<z<0.4 | 48.5 [2.3] | 31.8 [0.4] | 0.0 [0.0] | 31.4 [1.0] | 10.2 [5.5] | 0.0 [0.0] | 63.6 [0.4] | 3.9 [1.0] | 4.3 [0.4] | 11.0
0.4<z<1.0 | 48.2 [2.2] | 29.3 [0.7] | 20.0 [0.0] | 40.4 [0.5] | 11.3 [6.3] | 4.7 [0.1] | 26.8 [0.0] | 6.6 [0.9] | 10.2 [1.9] | 12.7
1.0<z<2.0 | 56.1 [2.9] | 38.9 [1.9] | 14.3 [0.0] | 49.5 [0.2] | 13.4 [6.4] | 4.0 [0.6] | 8.0 [0.0] | 5.3 [0.2] | 5.6 [1.3] | 13.5
2.0<z<3.6 | 70.6 [5.4] | 59.8 [6.6] | 6.5 [0.0] | 57.0 [0.1] | 28.6 [11.4] | 6.1 [1.7] | 0.0 [0.0] | 5.3 [0.0] | 9.0 [1.2] | 26.3

Figure 26: [COSMIC - eBHNS] For varied redshift slices, each panel corresponds
to a breakdown of the impact parameter distributions for different galaxy
masses (row), and star formation rates (column). Each histogram is composed of
the number of binaries with a given impact parameter normalised over the total
summation of merging systems within the respective panel and redshift. The
dashed vertical line corresponds to the galaxy averaged impact parameter
corresponding to the $P_{\rm chance}$ limit used to identify systems that may
be potentially classed as "host-less". A numerical breakdown of each panel is
shown in Table 10.  Table 10: [COSMIC - eBHNS] A breakdown of the relative
binary mergers occurring within each panel of Figure 26. The top division
refer to the relative fractions of binaries coalescing for a given redshift
slice. The last column indicates the population of the binaries contained
within each slice as a percentage of the cumulative population. The second
section of the table shows the estimated host-less population per panel
[overall fraction corresponding to the redshift slice] based on the criteria
defined in Section 3.2.

z | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | Total Population (%)
---|---|---|---|---|---|---|---|---|---|---
Binary Fraction (%)
0.0<z<0.4 | 2.6 | 0.8 | 0.0 | 3.5 | 52.5 | 0.2 | 1.5 | 29.5 | 9.4 | 5.7
0.4<z<1.0 | 2.0 | 1.6 | 0.0 | 1.6 | 55.0 | 2.5 | 0.3 | 16.8 | 20.3 | 39.4
1.0<z<2.0 | 1.3 | 2.5 | 0.0 | 0.5 | 41.5 | 18.6 | 0.1 | 5.5 | 30.1 | 45.8
2.0<z<3.6 | 3.0 | 5.7 | 0.0 | 0.2 | 20.2 | 40.7 | 0.0 | 1.4 | 28.9 | 9.0
| Host-less Fraction/Panel [Relative to z-slice] (%) | Total Host-less Fraction (%)
0.0<z<0.4 | 52.0 [1.3] | 30.4 [0.2] | 0.0 [0.0] | 60.4 [2.1] | 19.8 [10.4] | 0.0 [0.0] | 55.6 [0.9] | 12.0 [3.5] | 8.1 [0.8] | 19.2
0.4<z<1.0 | 49.2 [1.0] | 19.9 [0.3] | 0.0 [0.0] | 48.1 [0.8] | 18.1 [10.0] | 8.3 [0.2] | 28.1 [0.1] | 14.8 [2.5] | 14.9 [3.0] | 17.8
1.0<z<2.0 | 64.9 [0.8] | 25.2 [0.6] | 0.0 [0.0] | 48.6 [0.2] | 15.0 [6.2] | 7.3 [1.4] | 38.5 [0.0] | 15.8 [0.9] | 5.8 [1.7] | 11.9
2.0<z<3.6 | 77.0 [2.3] | 54.6 [3.1] | 0.0 [0.0] | 62.5 [0.1] | 26.9 [5.4] | 6.0 [2.4] | 0.0 [0.0] | 19.4 [0.3] | 4.0 [1.2] | 14.8

Figure 27: [COSMIC] Similar to Figure 13.

Figure 28: [COSMIC] Similar to Figure 11.

Figure 29: [COSMIC] Similar to Figure 10.

Figure 30: [COSMIC] Similar to Figure 14. Refer to Table 5 for sources used.