# Gravitational waves from the remnants of the first stars in nuclear star
clusters

Boyuan Liu 1 and Volker Bromm1
1Department of Astronomy, University of Texas, Austin, TX 78712, USA
E-mail<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

We study Population III (Pop III) binary remnant mergers in nuclear star
clusters (NSCs) with a semi-analytical approach for early structure formation.
Within this framework, we keep track of the dynamics of Pop III binary
(compact object) remnants during cosmic structure formation, and construct the
population of Pop III binary remnants that fall into NSCs by dynamical
friction of field stars. The subsequent evolution within NSCs is then derived
from three-body encounters and gravitational-wave (GW) emission. We find that
7.5% of Pop III binary remnants will fall into the centres ($<3\ \rm pc$) of
galaxies. About $5-50\%$ of these binaries will merge at $z>0$ in NSCs,
including those with very large initial separations (up to 1 pc). The merger
rate density (MRD) peaks at $z\sim 5-7$ with $\sim 0.4-10\ \rm yr^{-1}\ \rm
Gpc^{-3}$, leading to a promising detection rate $\sim 170-2700\ \rm yr^{-1}$
for 3rd-generation GW detectors that can reach $z\sim 10$. Low-mass ($\lesssim
10^{6}\ \rm M_{\odot}$) NSCs formed at high redshifts ($z\gtrsim 4.5$) host
most ($\gtrsim 90$%) of our mergers, which mainly consist of black holes (BHs)
with masses $\sim 40-85\ \rm M_{\odot}$, similar to the most massive BHs found
in LIGO events. Particularly, our model can produce events like GW190521
involving BHs in the standard mass gap for pulsational pair-instability
supernovae with a MRD $\sim 0.01-0.09\ \rm yr^{-1}\ Gpc^{-3}$ at $z\sim 1$,
consistent with that inferred by LIGO.

###### keywords:

early Universe – dark ages, reionization, first stars – gravitational waves –
galaxies: nuclei

††pubyear: 2021††pagerange: Gravitational waves from the remnants of the first
stars in nuclear star clusters–Gravitational waves from the remnants of the
first stars in nuclear star clusters

##  1 Introduction

The detection of gravitational waves (GWs) from mergers of compact objects,
such as black holes (BHs) and neutron stars (NSs), has opened a new
observational window in astrophysics, cosmology and fundamental physics
(reviewed by, e.g. Barack et al. 2019). These ripples of spacetime carry
valuable information for a variety of astrophysical processes, such as
formation and evolution of compact object binaries, BH mass function, cosmic
expansion and structure/star formation (e.g. Dvorkin et al. 2016; Fishbach et
al. 2018; Mapelli et al. 2019; Vitale et al. 2019; Perna et al. 2019;
Safarzadeh & Berger 2019; Farr et al. 2019; Adhikari et al. 2020; Tang et al.
2020; Safarzadeh 2020; Bouffanais et al. 2021; De Luca et al. 2021; Fishbach
et al. 2021; Mastrogiovanni et al. 2021; Wang et al. 2021). The GW window
ideally complements the electromagnetic (EM) window with two advantages that
are particularly important at high redshifts: (i) EM signals decay rapidly
with (luminosity) distance ($\propto d_{L}^{-2}$) such that EM observations at
high-$z$ are significantly biased towards massive/luminous systems, while the
amplitudes of GWs decay slower ($\propto d_{L}^{-1}$), making it easier for GW
observations to reach high-$z$. (ii) Compact object binaries formed in the
early Universe can merge at later times such that imprints of high-$z$
processes can be inferred from GW events detected in the local Universe.

These unique features of GWs make them a promising probe to the first stars,
the so-called Population III (Pop III) with extremely low or zero metallicity
($Z\lesssim 10^{-6}-10^{-4}\ \rm Z_{\odot}$, reviewed by e.g. Bromm et al.
2009; Bromm 2013; Haemmerlé et al. 2020), for which direct observations in the
EM window will be challenging even with the James Webb Space Telescope (JWST)
(Gardner et al., 2006; Schauer et al., 2020). Compared with Population I/II
(Pop I/II) stars formed in metal-enriched gas, Pop III stars produce compact
object remnants much more efficiently because of their lack of strong mass
loss and top-heavy initial mass function (IMF) resulting from the insufficient
cooling of primordial gas. This makes them ideal progenitors of compact object
mergers that can potentially account for a significant fraction of the GW
events detected at $z\lesssim 1$ so far by LIGO (Nitz et al., 2021), and those
to be discovered at higher redshifts. The recent detection of the special
event GW190521 with unusual BH masses $85_{-14}^{+21}\ \rm M_{\odot}$ and
$66_{-18}^{+17}$ (Abbott et al., 2020a, b) further highlights the importance
of Pop III remnants. These BH masses are forbidden for Pop I/II by standard
pulsational pair-instability supernova (PPISN) models (e.g. Heger et al. 2003;
Belczynski et al. 2016; Woosley 2017; Marchant et al. 2019). Pop III stars, on
the other hand, can produce BHs in this mass range by retaining the hydrogen
envelope before collapse and avoiding the PPISN regime due to their
compactness and strongly reduced mass loss (Farrell et al., 2021; Tanikawa et
al., 2021a). Furthermore, the observed merger rate density (MRD) for
GW190521-like events can also be reproduced considering a Pop III origin
(Kinugawa et al., 2021a; Liu & Bromm, 2020c).

The GW signals from Pop III binary remnant mergers have been intensively
studied for the classical binary stellar evolution (BSE) channel (e.g.
Kinugawa et al. 2014, 2015; Hartwig et al. 2016; Belczynski et al. 2017;
Tanikawa et al. 2021b; Hijikawa et al. 2021; Kinugawa et al. 2021b). In this
scenario for close binaries in isolation, interactions such as mass transfer,
tidal effects, and common envelope evolution shrink the binary orbits to
facilitate merging within a Hubble time. The predictions from these studies
exhibit significant discrepancies due to different assumptions on the Pop III
star formation (SF) history, initial binary statistics and BSE parameters. For
instance, even if the Pop III SF history is constrained, the predicted local
($z\sim 0$) MRDs have large scatters of two orders of magnitude, such that it
is unclear which fraction (none or up to 50% ) of the LIGO events can be
attributed to Pop III remnants.

Unfortunately, it is difficult to reduce the uncertainties in the binary
statistics and BSE parameters without direct observations of Pop III stars.
Besides, all these studies rely on the assumption that Pop III binaries are
similar to their Pop I/II counterparts (with observational constraints), which
can be an over-idealization, given that Pop III SF occurs under markedly
different conditions (see e.g. Chon et al. 2021). Actually, recent high-
resolution (radiative) hydrodynamic simulations of Pop III SF and N-body
simulations of young Pop III star clusters (Susa, 2019; Sugimura et al., 2020;
Liu et al., 2021) indicate that close binaries of Pop III stars are likely
rare because of the expansion of the system by angular momentum conservation
during protostellar accretion.

The challenges for the BSE channel motivate us to explore alternative pathways
of forming Pop III binary remnant mergers that can also work even if Pop III
binaries are initially wide (with separations $\gtrsim 10\ \rm au$). Dynamical
hardening (DH) in nuclear star clusters (NSCs) has been shown to be an
efficient mechanism of shrinking binary orbits for Pop I/II remnants and
supermassive BHs (see e.g. Sesana & Khan 2015; Antonini & Rasio 2016; Choksi
et al. 2019; Sedda 2020; Mapelli et al. 2021). In this scenario, massive, hard
binaries of compact objects will be further hardened by three-body encounters
with surrounding stars, which can drive them to mergers within a Hubble time
under the high stellar densities of NSCs. In our case, Pop III (binary)
remnants need to fall into NSCs for this mechanism to work, which is a complex
process involving the dynamics of Pop III remnants during cosmic structure
formation. The exploratory study by Liu & Bromm (2020c) applies this NSC-DH
channel to Pop III remnants, based on a simple model for the dynamics of Pop
III remnants in their host galaxies, finding that it is possible to explain
the GW190521 event. In this work, we design a more advanced model to keep
track of the motions of Pop III remnants driven by dynamical friction (DF)
from field stars, coupled with a semi-analytical framework for early structure
formation based on halo merger trees. By taking into account the cosmological
context for the dynamics of Pop III remnants in their host galaxies and during
galaxy mergers, we can predict when and where Pop III remnants fall into NSCs.
The subsequent evolution of Pop III binary remnants inside NSCs is then
followed by post-processing to predict their mergers.

This paper is organised as follows. In Section 2, we describe our semi-
analytical framework for early structure formation (Sec. 2.1) and models for
the in-fall of Pop III binary remnants into NSCs (Sec. 2.2) and their
subsequent evolution (Sec. 2.3). In Section 3, we present our main results for
the properties of host galaxies/NSCs at the moments of NSC in-fall (Sec. 3.1)
and mergers of Pop III binary remnants in NSCs, compared with previous
theoretical predictions for the BSE channel and LIGO observations (Sec. 3.2).
Finally, in Section 4, we summarize the key features of the NSC-DH channel for
Pop III binary remnant mergers and discuss their implications for GW astronomy
in the early Universe.

##  2 Methodology

In this section, we present our method of modelling Pop III binary remnant
mergers in (high-$z$) NSCs. We start with a merger-tree based framework for
early structure formation that self-consistently models SF, stellar feedback,
and particularly, Pop III binary stars and the relevant binary remnants (Sec.
2.1). Then we keep track of the dynamics of such Pop III binary remnants in
their host galaxies, modelling the population of binaries falling into NSCs by
DF from field (low-mass Pop I/II) stars (Sec. 2.2). Finally, with an
observation-based NSC model we follow the evolution of such binaries in NSCs
(Sec. 2.3) via DH and (potential) ejection by three-body encounters, as well
as GW emission. The binary properties involved are derived from N-body
simulations of Pop III star clusters in Liu et al. (2021).

###  2.1 Early structure formation

Following earlier work (Hartwig et al., 2015; Hartwig et al., 2016), we model
early structure formation with halo merger trees generated by the galform code
(Parkinson et al., 2008). For simplicity, we simulate 500 trees starting from
an initial redshift of $z_{i}=42$, each targeted at a Milky-Way (MW)-like halo
of total mass $M_{\rm h}=1.26\times 10^{12}\ \rm M_{\odot}$ at $z_{\rm
vir}\simeq 2.5$. The mass resolution of the trees is $\Delta M_{\rm h}=5\times
10^{5}\ \rm M_{\odot}$, and the redshift resolution is $\Delta z\simeq 0.15$,
corresponding to timesteps of $\sim 0.3-180\ \rm Myr$. Since such MW-like
haloes with $M_{\rm h}\sim 10^{12}\ \rm M_{\odot}$, corresponding to $\sim
1-2\sigma$ peaks at $z\lesssim 2.5$, are massive enough to capture small-scale
perturbations that collapse at high-$z$, we expect their average assembly
histories to be cosmologically representative before turn-over ($z\gtrsim 5$),
when most Pop III stellar remnants are formed. Our results for the low-
redshift regime ($z\lesssim 5$) are derived from simple assumptions and
extrapolation (see Sec. 2.2), which may be incomplete. As the goal of this
study is to demonstrate the NSC-DH channel for Pop III remnant mergers, we
defer more accurate modelling at low-$z$ with complete halo populations to
future work.

The merger trees are coupled with a customized model for primordial SF that
self-consistently takes into account metal enrichment and Lyman-Werner (LW)
feedback (for details, see Hartwig et al. 2015, HT15 henceforth). We have
updated this model to achieve better agreement with observations (Campisi et
al., 2011; Madau & Dickinson, 2014; Aghanim et al., 2020), and implemented
recipes for Pop I/II SF, important for the dynamics of Pop III remnants in
their host galaxies. Below we describe our main refinements to the HT15 model.

Same as HT15, we use the cosmological parameters from Planck 2013 results (Ade
et al., 2014) for a flat $\Lambda$CDM universe: $H_{0}=67.77\ \rm km\ s^{-1}\
Mpc^{-1}$, $\Omega_{\rm m}=0.309$, $\Omega_{\rm b}=0.048$, $n_{\rm s}=0.9611$
and $\sigma_{8}=0.8288$. The dark matter power spectrum is calculated with the
camb code (Lewis et al., 2000) in the wave number ($k$) range $k/h\sim
10^{-6}-10^{6}\ \mathrm{Mpc^{-1}}$, where $h\equiv H_{0}/(100\ \rm km\ s^{-1}\
Mpc^{-1})$

####  2.1.1 Pop III star formation

The original Pop III SF model in HT15 does not consider the cooling/collapse
timescale of the star-forming cloud in the host halo, and therefore,
significantly overproduces Pop III stars at very high redshifts ($z\sim
15-35$) by instantaneous SF. This leads to early reionization ($z\gtrsim 8$)
with a high Thomson optical depth $\tau\gtrsim 0.09$, inconsistent with the
recent measurement by Planck, $\tau=0.0544\pm 0.0073$ (Aghanim et al., 2020),
which implies a late reionization completed at $z\sim 5.5$. In the refined
model, for each halo that meets the Pop III SF criteria, we delay Pop III SF
by a collapse timescale $t_{\rm col}=\max[10^{\sigma_{\rm col}}t_{\rm dyn},\rm
3\ Myr]$. Here $t_{\rm dyn}=t_{H}/10^{3/2}$, given the Hubble time $t_{H}$, is
an estimation of the free-fall timescale in the inner region ($r\lesssim
0.1R_{\rm vir}$, given the halo virial radius $R_{\rm vir}$) where the star-
forming cloud resides (Bromm et al., 2002). $\sigma_{\rm col}$ is a random
number generated from a uniform distribution in [-0.15,0.15], introduced to
capture the diversity in halo assembly histories and avoid spurious
fluctuations in the Pop III SF rate density (SFRD). The probability
distribution of $\sigma_{\rm col}$ is chosen (by trial-and-error) to reproduce
the SFRD from the cosmological simulation in Campisi et al. (2011, see below).
We also apply a lower bound of $10^{6}\ \rm M_{\odot}$ to the critical halo
mass threshold for Pop III SF, considering the effect of baryon-dark matter
streaming motion (Schauer et al., 2019).

To keep track of Pop III binary remnants in structure formation, we design a
scheme of sampling binary stars in Pop III star-forming haloes, based on
N-body simulations of Pop III star clusters from Liu et al. (2021). Following
the fiducial model in Liu et al. (2021), we assume the same Pop III IMF with a
log-flat form ($dN_{\star}/dm_{\star}\propto m_{\star}^{-1}$) in the range
between $m_{\min}=1\ \rm M_{\odot}$ and $m_{\max}=170\ \rm M_{\odot}$. In this
model, the average fraction of stellar mass in binaries is $f_{\rm B}=0.69$.
The mass range adopted here, $m_{\star}\sim 1-170\ \rm M_{\odot}$, is
generally consistent with that inferred from the MW stellar metallicity
distribution function in the low-metallicity regime, $m_{\star}\sim 2-180\ \rm
M_{\odot}$ (Tarumi et al., 2020). For each Pop III star-forming halo, the mass
$M_{\rm III}$ and number $N_{\star}$ of Pop III stars are given by the IMF
sampling scheme in HT15 for individual stars. If $N_{\star}>1$, we sample
binary stars until the total mass in binaries exceeds $f_{\rm os}f_{\rm
B}M_{\rm III}$, or the number of stars left is less than 2. Here, we introduce
$f_{\rm os}\simeq 0.3$ to suppress overshooting in our numerical sampling,
determined by trial-and-error, to make sure that the total mass of Pop III
binaries (ever formed) agrees with the total mass of Pop III stars multiplied
by $f_{\rm B}$. For each binary, we first generate the primary mass
$m_{\star,1}$ from a uniform distribution in [$m_{\min}$,$m_{\rm max}$], and
then generate the secondary mass from a uniform distribution in
[$m_{\min}$,$m_{\star,1}$]. The mass distributions adopted here are meant to
reproduce those in N-body simulations (see fig. 9 and 10 in Liu et al.
2021)111Here we have ignored the effect of stellar collisions, as they are
rare ($\lesssim 1\%$) in most cases (Liu et al., 2021)., which reflect the
fact that massive stars are favored for binaries.

Following Liu & Bromm (2020c, see their fig. 1), we map binary stars to binary
remnants with the remnant-initial mass relation $m_{\rm rem}(m_{\star})$ based
on the fitting formula from Tanikawa et al. (2020) for their $Z=10^{-6}\ \rm
Z_{\odot}$ models222In this case, NSs and BHs ($m_{\rm rem}\gtrsim 1.4\ \rm
M_{\odot}$) are formed at $m_{\star}\gtrsim 10\ \rm M_{\odot}$. Stars with
$m_{\star}\sim 25-85\ \rm M_{\odot}$ collapse directly into BHs. The gap for
(pulsational) pair-instability supernovae is $\sim 115-230\ (85-115)\ \rm
M_{\odot}$., taking into account PPISNe: $m_{\star,i}\rightarrow m_{i}\equiv
m_{\rm rem}(m_{\star,i})$ ($i=1,\ 2$)333We have ignored the effects of mass
transfer and supernova explosions on the survival and properties of binary
remnants. It is found in Liu et al. (2021) that mass transfer from Roche lobe
overflow is rare for Pop III binaries. It will be shown in Sec. 3.2.2 (Fig.
20) that most mergers of Pop III remnants involve BHs with $m_{\rm rem}\gtrsim
25\ \rm M_{\odot}$, which tend to form directly without supernova explosions..
As we are concerned with the GWs from binary compact object merges, we only
keep track of the remnant binaries with progenitor masses $m_{\star,i}>10\ \rm
M_{\odot}$ ($i=1,\ 2$), and total masses $m_{\bullet}\equiv m_{1}+m_{2}>3\ \rm
M_{\odot}$. On average, $\sim 10^{4}$ such binaries are formed in each tree.
It will be shown below that the distribution of total mass $m_{\bullet}$ in
remnant binaries obtained by our sampling scheme is generally consistent with
that derived from the binary statistics of N-body simulations in Liu et al.
(2021).

####  2.1.2 Pop I/II star formation

For simplicity, Pop I/II SF is modelled with an input SFRD in HT15 to
calculate the relevant contribution to the Thomson optical depth, while not
explicitly taking into account the chemical and radiative feedback from Pop
I/II stars. In our case, we need to derive the properties of individual
galaxies that potentially host NSCs in order to keep track of the dynamics of
Pop III (binary) remnants. Therefore, we now include a self-consistent model
for Pop I/II SF and feedback based on the stellar-halo mass relation (SHMR) in
Behroozi et al. (2019). We adopt their fitting formula for the true mass
values including both star-forming and quiescent galaxies and excluding
intrahalo light (see column 8 of their table J1), as shown in Fig. 1. We
further impose a lower limit for the SF efficiency (SFE)444Ratio of the
stellar to total baryon mass., $\eta_{\star,0}$, in (the extrapolation of)
this relation before reionization ($z>6$), which is particularly relevant for
high-$z$ haloes with relatively low masses ($M_{\rm h}\lesssim 10^{9}\ \rm
M_{\odot}$). We assume that all haloes enriched by metals with $M_{\rm
h}>10^{6}\ \rm M_{\odot}$ can host Pop I/II stars. The stellar masses of such
haloes are assigned and updated according to the SHMR: $M_{\star}\equiv
M_{\star}(M_{\rm h},z)$. For each new (existing) enriched halo, (the increase
of) $M_{\star}$ (with respect to the stellar mass carried by child haloes at
the previous timestep) is attributed to newly-formed Pop I/II stars at the
current timestep. To take into account the effect of reionization, we prevent
any halo with $M_{\rm h}<6.7\times 10^{8}\ \mathrm{M_{\odot}}\
[(1+z)/5]^{-3/2}$ from forming new stars at $z<6$ based on reionzation models
in Pawlik et al. (2015); Pawlik et al. (2017); Benitez-Llambay & Frenk (2020);
Hutter et al. (2021).

For chemical feedback from Pop I/II stars, we simply increase the volume
filling fraction of metals from Pop III enrichment by a factor of 5 to mimic
the additional enrichment from galactic outflows of Pop I/II galaxies formed
in regions initially enriched by Pop III stars. While for LW feedback, we add
the contribution from Pop I/II stars to the LW background (in units of
$10^{-21}\ \rm erg\ s^{-1}\ cm^{-2}\ Hz^{-1}\ sr^{-1}$) following Johnson et
al. (2013):

$\displaystyle J_{\rm
21,II,bg}=0.3\left(\frac{1+z}{16}\right)^{3}\left[\frac{{\rm SFRD_{\rm
I/II}}(z)}{10^{-3}\ \rm M_{\odot}\ yr^{-1}\ Mpc^{-3}}\right]\ ,$ (1)

where ${\rm SFRD_{\rm I/II}}(z)$ is the co-moving SF rate density of Pop I/II
stars self-consistently predicted by the above model, based on the SHMR.

Figure 1: SHMR from Behroozi et al. (2019) at $z=0$ (solid), 2 (dashed), 4
(dashed-dotted), 6 (dotted), 8 (long-dashed) and 10 (nest-dashed), given the
(pre-reionization) lower limit of SFE $\eta_{\star,0}=0.003$. The thick dashed
lines indicate constant SFEs in the range of $10^{-3}-1$. Figure 2: Co-moving
SFRDs for all stars (solid) and Pop III stars (dashed) with
$\eta_{\star,0}=0.003$ and $\eta_{\star,\rm III}=0.001$. $3\sigma$
uncertainties are illustrated with the shades. The reference total SFRD from
Campisi et al. (2011) is shown with the dashed-dotted curve. The observational
results in Madau & Dickinson (2014), inferred from UV and IR galaxy surveys
such as Finkelstein (2016, data points), are plotted as the long-dashed curve
(with scatters of 0.2 dex embodied by the shaded region). For comparison, we
also shown the total and Pop III SFRDs from Liu & Bromm (2020a, b) with the
nest-dashed and dotted curves, respectively. The thin vertical line denotes
the turn-over redshift $z\simeq 4.6$ of the target halo below which our
results may not be cosmologically representative. Figure 3: Thomson optical
depth as a function of redshift. The results with all contributions from Pop
III and I/II stars are shown with the solid and dotted curves, assuming a Pop
I/II escape fraction of $f_{\rm esc,I/II}=0.1$, and considering a range of Pop
III escape fractions of $f_{\rm esc,III}\sim 0.1-0.3$, where 0.125 is the best
fit value (solid). The results with only Pop I/II stars are shown with the
dashed-dotted curve (again for $f_{\rm esc,I/II}=0.1$). The recent measurement
by Planck (Aghanim et al., 2020) is shown with the horizontal dashed line,
whose $1\sigma$ scatters are denoted by the shaded region. For reference,
$z=5.5$ is shown with the thin vertical line.

Finally, we choose the (pre-reionization) lower limit for the Pop I/II SFE,
$\eta_{\star,0}$, and the Pop III SFE, $\eta_{\rm III}$, to reproduce the
total (Pop III+Pop I/II) SFRD inferred from observations and simulations.
Following HT15, we take the results from the cosmological simulation in
Campisi et al. (2011) as the reference, which are consistent with the
observational results in Madau & Dickinson (2014) at $z\lesssim 4$ where
observational constraints are most reliable. As shown in Fig. 2, the reference
SFRD can be well reproduced with $\eta_{\star,0}=0.003$ and $\eta_{\star,\rm
III}=0.001$, particularly for $z\gtrsim 5$ when the target halo has not
reached turn-over. We also compare our SFRDs with those derived from the
recent cosmological simulation FDbox_Lseed in Liu & Bromm (2020a, b). Our
results are generally consistent with simulation data within an order of
magnitude, especially around the peak of Pop III SF at $z\sim 6-12$. It will
be shown below that this is also the peak of NSC capture of Pop III remnants
(Sec. 3.1). Note that the early cosmic SF history is still highly uncertain,
particularly for Pop III stars with very limited constraints from observations
and significant difference in the predictions from theoretical studies with
different numerical models and assumptions for SF and stellar feedback (see
e.g. fig. 5 and 13 of Liu & Bromm 2020b). Hopefully more robust observational
constraints can be obtained in the JWST era. We defer more detailed
calibrations with SFRD to future work. Instead, when comparing our predicted
MRD of Pop III remnants with literature results, we constrain the Pop III SF
history by scaling their SFRDs (and MRDs) such that the scaled integrated Pop
III stellar mass densities ($\mathrm{ISMD}_{\rm
III}\equiv\int_{0}^{\infty}\mathrm{SFRD}_{\rm III}|dt/dz|dz$) are identical to
our value $\simeq 7\times 10^{4}\ \mathrm{M_{\odot}\ Mpc^{-3}}$.

Actually, the Thomson optical depth measured by Planck, $\tau=0.0544\pm
0.0073$ (Aghanim et al., 2020), can also be well reproduced by our SF model
under reasonable assumptions for the escape fractions of ionizing photons.
Fig. 3 shows the evolution of $\tau$ with redshift predicted by our merger
trees, which agrees well with the Planck result assuming escape fractions
$f_{\rm esc,II}=0.1$ for Pop I/II and $f_{\rm esc,III}\sim 0.1-0.3$ for Pop
III. With these escape fractions, the predicted volume filling fraction of
ionized gas reaches $\simeq 1\ (0.5)$ at $z\simeq 5.5\ (8)$, consistent with
the picture of late reionization. The escape fractions adopted here are within
the typical range $f_{\rm esc}\sim 0.1-0.7$ seen in simulations and semi-
analytical models (e.g. So et al. 2014; Paardekooper et al. 2015; Inayoshi et
al. 2016; Visbal et al. 2020; Katz et al. 2021).

###  2.2 In-spiral of Pop III remnants by dynamical friction

Once formed in small structures, Pop III (binary) remnants555We ignore the
growth of Pop III remnants through accretion, as previous studies based on 3D
cosmological simulations have found that stellar-mass Pop III remnants can
hardly grow via accretion at high-$z$ (e.g. Johnson & Bromm 2007; Alvarez et
al. 2009; Hirano et al. 2014; Smith et al. 2018). will fall into larger
structures during halo mergers and meanwhile spiral towards galaxy centres by
DF from field stars666For simplicity, we ignore the DF from gas, whose effect
is likely minor in the long-term (Chen et al., 2021). We also ignore the
effect of dark matter, as we are only concerned with the central region
($r\lesssim 300\ \rm pc$) dominated by baryons in most cases, where the DF
timescale is shorter than the Hubble time., where they can further fall into
NSCs. To model this process, we need to consider how remnants are (initially)
distributed in haloes during structure formation, as well as the host galaxy
properties.

For the former, we adopt the spatial distribution of remnants from the
cosmological simulation in Liu & Bromm (2020a, b), as shown in Fig. 4. The
distribution has little evolution at $z\sim 4-14$ in units of $r/R_{\rm vir}$,
where $r$ is the physical distance to the galaxy/halo centre and $R_{\rm vir}$
is the halo virial radius. All newly-formed Pop III remnants are distributed
in their initial hosts following this distribution. When two haloes merge, the
remnants at $r>3\ \rm pc$ in the smaller progenitor are distributed into the
merged halo with the same distribution, while the locations/orbits of the
remnants in the larger progenitor remain unaffected by the merger.

Figure 4: Enclosed fractions of Pop III remnants in terms of $r/R_{\rm vir}$
given the physical distance $r$ to the galaxy/halo centre and the halo virial
radius $R_{\rm vir}$, for 12 snapshots at $z\sim 4-14.2$ from the cosmological
simulation in Liu & Bromm (2020a, b). The combined distribution for all
snapshots is shown with the dotted curve, which is approximated with a triple
power-law fit (thick solid) in our numerical model. The inner power-law slope
is $\sim 2$.

For the latter, we adopt the galaxy model in Arca-Sedda & Capuzzo-Dolcetta
(2014, AC14, henceforth), which successfully reproduces several observed
scaling relations between NSCs and host galaxies for the globular cluster (GC)
in-spiral scenario. In this model, a galaxy is characterized by a Dehnen
sphere with density profile (Dehnen, 1993)

$\displaystyle\rho_{\gamma}(r|M,R,\gamma)=\frac{\rho_{\gamma
0}}{(r/R)^{\gamma}(1+r/R)^{4-\gamma}}\ ,$ (2)

where $\rho_{\gamma 0}=(3-\gamma)M/(4\pi R^{3})$, given $M$, $R$ and $\gamma$
as the total mass, scale length and inner slope. For each galaxy777Each halo
labelled as enriched by the stochastic enrichment model in HT15 hosts a
galaxy. To take into account self enrichment, we further assign a galaxy to
each halo that has hosted Pop III SF even if it is not enriched according to
the stochastic model. The (Pop I/II) galaxy formation in such self-enriched
haloes is delayed by $100\ \rm Myr$ with respect to the initial Pop III SF
event, which is the typical re-collapse timescale after supernovae., the
galaxy (stellar) mass $M_{\star}$ is given by our Pop I/II SF model. We then
estimate the galaxy size (scale length) $R_{\rm g}$ with the size-mass
relation (AC14)

$\displaystyle R_{\rm g}/\mathrm{kpc}=2.37\ [2^{1/(3-\gamma_{\rm
g})}-1]M_{\rm\star,11}^{0.14}\ ,$ (3)

where $M_{\star,11}\equiv M_{\star}/(10^{11}\ \rm M_{\odot})$, and the galaxy
inner slope $\gamma_{\rm g}$ is generated from a uniform distribution in the
ranges of $0-0.5$ for $M_{\star,11}<0.1$ and $0.5-1$ for $M_{\star,11}\geq
0.1$ (for newly-formed galaxies). During halo mergers, $\gamma_{\rm g}$ is
inherited from the larger progenitor (i.e. along the main branch). Since we
are concerned with high-$z$ low-mass galaxies, we always have $\gamma_{\rm
g}\sim 0-0.5$. Considering a broader range $\gamma_{\rm g}\sim 0-2$ hardly
changes our results.

Now that the initial conditions and galaxy backgrounds are set up, we can
estimate the DF timescale, $\tau_{\rm DF}$, for remnants to experience their
in-spiral. For an object of mass $m$ with initial apocentric distance $r$ and
velocity $v$ with respect to the centre in a system of total mass $M$ and size
$R$, the DF timescale given by Chandrasekhar’s formula is (Binney & Tremaine,
2011)

$\displaystyle\frac{\tau_{\rm
DF}}{\mathrm{Myr}}=\frac{34.2}{\ln\Lambda}\left(\frac{r}{\rm 3\
pc}\right)^{2}\left(\frac{v}{\rm 10\ km\ s^{-1}}\right)\left(\frac{m}{100\ \rm
M_{\odot}}\right)^{-1}\ ,$ (4)

where $\ln\Lambda\sim\ln[Mr/(0.8mR)]$ is the Coulomb logarithm, and the
initial velocity is estimated with the velocity dispersion
$v\sim\sigma\sim\sqrt{GM/R}$.

The work of Chandrasekhar has been generalized to describe DF in both cusped
and cored density profiles by Arca-Sedda et al. (2015); Arca-Sedda (2016)
which give an updated DF timescale formula for the Dehnen profile (Equ. 2):

$\displaystyle\frac{\tau_{\rm DF}}{\mathrm{Myr}}=0.3g\sqrt{\frac{r^{3}}{\rm 1\
kpc^{3}}M_{11}^{-1}}\left(\frac{m}{M}\right)^{\alpha}\left(\frac{r}{R}\right)^{\beta}\
,$ (5)

where $M_{11}=M/(10^{11}\ \rm M_{\odot})$, $\alpha=-0.67$, $\beta=1.76$, and

$\displaystyle g$ $\displaystyle\equiv g(e_{\rm if},\gamma)$
$\displaystyle=(2-\gamma)\left\\{a_{1}\left[\frac{1}{(2-\gamma)^{a_{2}}}+a_{3}\right](1-e_{\rm
if})+e_{\rm if}\right\\}\ .$ (6)

Here, $e_{\rm if}$ is the eccentricity of the in-fall orbit, $a_{1}=2.63\pm
0.17$, $a_{2}=2.26\pm 0.08$, and $a_{3}=0.9\pm 0.1$. Note that this formula is
only valid for $0<\gamma<2$. In our implementation, we take the minimum of the
above two DF timescales, which are coupled to our galaxy model via
$m=m_{\bullet}$, $M=M_{\star}$, $R=R_{\rm g}$ and $\gamma=\gamma_{\rm g}$.
Actually, in most cases, the recent formula (5) gives a shorter DF timescale.

Figure 5: DF timescale as a function of galaxy mass, for $m_{\bullet}=20\ \rm
M_{\odot}$, $\gamma_{\rm g}=0$ and $e_{\rm if}=0$, at $r\sim 1,\ 3,\ 10\
\text{and}\ 30\ \rm pc$ (from bottom to top), given by the Chandrasekhar’s
formula (Binney & Tremaine, 2011, BT08 solid, see Equ. 4), and more recent
work in Arca-Sedda (2016, AS16, dashed, see Equ. 5).

Fig. 5 demonstrates the dependence of $\tau_{\rm DF}$ on $M_{\star}$ for
$m_{\bullet}=20\ \rm M_{\odot}$, $\gamma_{\rm g}=0$ and $e_{\rm if}=0$. It
turns out that smaller galaxies are more efficient at DF. For most binaries
with $m_{\bullet}\gtrsim 20\ \rm M_{\odot}$ in our case, we have $\tau_{\rm
DF}\lesssim 100\ \rm Myr$ for $r\lesssim 3\ \rm pc$, shorter than the typical
timescale for galaxy evolution, and $\tau_{\rm DF}\gtrsim 10\ \rm Gyr$ for
$r\gtrsim 300\ \rm pc$, comparable/longer than the Hubble time. Therefore, we
only compute the insprial of remnants in the distance range $r\sim 3-300\ \rm
pc$. The evolution of $r$ is described by the differential equation

$\displaystyle\frac{dr}{dt}\simeq-\frac{r}{\tau_{\rm DF}(r|m,M,R,\gamma,e_{\rm
if})}\ ,$ (7)

with $m=m_{\bullet}$, $M=M_{\star}$, $R=R_{\rm g}$ and $\gamma=\gamma_{\rm
g}$, where $M_{\star}$ and $R_{\rm g}$ are updated on-the-fly. As a
conservative choice, we set $e_{\rm if}=0$. We have verified that considering
higher $e_{\rm if}$ up to 0.9 only slightly increases the fraction of remnants
with $r<3\ \rm pc$. We record all remnants that fall into galaxy centres
($r<3\ \rm pc$), together with their host properties. For remnants within the
target halo at $z\simeq 2.5$, we further evolve them to $z=0$ by post-
processing with galaxy properties fixed at $z\simeq 2.5$, and record the ones
that will fall into the centre at $z>0$. This offers a rough estimation for
the low-redshift regime, which will be explored in future work.

Figure 6: Mass distributions of Pop III binary remnants at galaxy centres
($r<3\ \rm pc$, solid) as well as formation sites, from our binary star
sampling scheme (dashed) and the N-body simulations in Liu et al. (2021,
histograms). Figure 7: Enclosed fractions of Pop III binary remnants in terms
of the physical distance $r$, at formation sites (solid) and the target halo
with (dashed) and without (dashed-dotted) DF. All distributions are normalized
to the total number of binary remnants formed in the merger trees. Without DF,
the inner distribution follows $f_{\bullet}(<r)\propto r^{2}$ approximately
(dotted), at both the target halo and formations sites, consistent with the
input simulation data (Fig. 4).

On average, 7.5% of Pop III remnants fall into the centres ($r<3\ \rm pc$) of
galaxies with $M_{\rm NSC}>10^{5}\ \rm M_{\odot}$ in our merger trees. These
remnants can further fall into NSCs, depending on the occupation frequency of
NSCs (see below). As DF is more efficient for more massive in-falling objects,
the mass distribution of remnants with $r<300\ \rm pc$ is more top-heavy
compared with the initial mass distribution at formation sites, as shown in
Fig. 6 in arbitrary units (a.u.). It is also shown that the mass distribution
obtained by our sampling scheme is generally consistent with that derived from
the N-body simulations in Liu et al. (2021), although the most massive
($m_{\bullet}\gtrsim 100\ \rm M_{\odot}$) and low-mass ($m_{\bullet}\lesssim
30\ \rm M_{\odot}$) objects are slightly overproduced.

To further illustrate the dynamics of Pop III remnants in their host galaxies
during structure formation, Fig. 7 shows the spatial distributions of Pop III
binary remnants at formation sites and the target halo. In general, due to
cosmic expansion and the nature of DF (see Fig. 5), it is more difficult for
Pop III remnants to enter inner regions of galaxies in more massive haloes at
later times. In the formation sites with $M_{\rm h}\sim 10^{6-7}\ \rm
M_{\odot}$ at $z\sim 5-20$, $\sim 90$% of remnants are initially within $300\
\rm pc$, while only $\sim 1-2$% of remnants can fall into the inner region
($r\lesssim 300\ \rm pc$) of the target halo ($M_{\rm h}\sim 10^{12}\ \rm
M_{\odot}$ at $z\simeq 2.5$). DF slightly increases the fraction in the target
halo but the effect is rather minor. This implies that NSCs formed in low-mass
galaxies at high-$z$ will be the dominant hosts of Pop III binary remnants,
which is confirmed by our results (see Sec. 3.1).

###  2.3 Evolution of compact object binaries in NSCs

For the $7.5\%$ of Pop III binary remnants at galaxy centres ($r<3\ \rm pc$),
a fraction of them will further fall into NSCs888We ignore the binary
evolution outside NSCs since most of our binary remnants are initially wide
(with separations larger than a few au) with very weak GW emission, and the
density of stars in the galaxy field is too low for efficient DH. in which
their evolution is driven by DF and three-body encounters from low-mass stars,
as well as GW emission. We follow this evolution with a NSC model based on
(local) observations (Neumayer et al., 2020) by post-processing. First, we
need to estimate the occupation frequency of NSCs, $f_{\rm NSC}$, which is
still uncertain for high-$z$ galaxies. We only consider NSCs with masses
$M_{\rm NSC}>10^{5}\ \rm M_{\odot}$. The lower bound here is consistent with
the typical mass of NSCs in the first galaxies predicted by theoretical models
(Devecchi & Volonteri, 2009; Devecchi et al., 2010, 2012), and the smallest
NSCs in observations. Note that the mass of Pop I/II star clusters in the
first galaxies is still uncertain, for which lower values ($\sim 3\times
10^{3}\ \rm M_{\odot}$) are predicted by high-resolution hydrodynamic
simulations (e.g. Safranek-Shrader et al. 2016). In the optimistic case, we
assume full occupation for galaxies with expected NSC masses $M_{\rm
NSC}>10^{5}\ \rm M_{\odot}$ (see below). In the conservative case, we adopt
the empirical occupation frequency as a function of galaxy mass, shown in Fig.
8. In our implementation of $f_{\rm NSC}$, we have $f_{\rm NSC}=0.15$ for
$M_{\star}\sim 10^{5}-3\times 10^{6}\ \rm M_{\odot}$, which is also consistent
with theoretical predictions for typical host haloes for these low-mass
galaxies with $M_{\rm h}\sim 10^{8}\ \rm M_{\odot}$ at $z\sim 10-20$ (Devecchi
et al., 2012, see their fig. 5) and observations of dwarf galaxies.

Figure 8: Occupation frequency of NSCs as a function of galaxy mass for early-
type (dashed) and late-type (dashed-dotted) galaxies in observations (Neumayer
et al. 2020, see their fig 3). The dotted curve shows the mean, combining the
two types of galaxies assuming a 1:1 mixture ratio, as observed in the general
field (Calvi et al., 2012). The thick solid curve is our approximation to the
observed mean, in which $f_{\rm NSC}=0.15$ at the low-mass end.

Next, we build a stochastic model for NSC properties from empirical NSC-host
scaling relations and the corresponding scatters. We again assume Dehnen
profile (Equ. 2) for NSCs. The NSC mass $M_{\rm NSC}$ is estimated with the
NSC-galaxy mass relation (Neumayer et al., 2020)

$\displaystyle\log(M_{\rm NSC}/\mathrm{M_{\odot}})=0.48\log[M_{\star}/(10^{9}\
\mathrm{M_{\odot}})]+6.51\ ,$ (8)

with scatters of $\sim 0.6$ dex. The (expected) NSC mass is assigned and
updated on-the-fly for each galaxy with $M_{\star}>10^{5}\ \rm M_{\odot}$. The
scatter is generated from an uniform distribution in logarithmic space, which
is also inherited along the main branch. We also require $M_{\rm NSC}\leq
M_{\star}$. Once $M_{\rm NSC}$ is known, we derive the NSC size following the
NSC size-mass relation (see fig. 7 of Neumayer et al. 2020). We fit the
observational data for effective radius $R_{\rm eff}$ and $M_{\rm NSC}$ with

$\displaystyle\log\left(\frac{R_{\rm eff}}{\rm
pc}\right)\simeq\begin{cases}0.54\ ,\ &M_{\rm NSC}<M_{\rm crit}\ ,\\\
0.34\log\left(\frac{M_{\rm NSC}}{\mathrm{M_{\odot}}}\right)-1.59\ ,\ &M_{\rm
NSC}\geq M_{\rm crit}\ ,\end{cases}$

where $M_{\rm crit}\simeq 2\times 10^{6}\ \mathrm{M_{\odot}}$. The standard
deviation around this fit is $\sim 0.32$ dex, which implies scatters of $\sim
0.55$ dex, assuming a uniform distribution (of scatters) in logarithmic space.
For a Dehnen sphere, the relation between scale length and effective radius is
(AC14)

$\displaystyle R_{\rm NSC}=(4/3)R_{\rm eff}\left[2^{1/(3-\gamma_{\rm
NSC})}-1\right]\ ,$ (9)

where the inner slope $\gamma_{\rm NSC}\sim 0.65-2.55$ tends to decrease with
higher $M_{\star}$ and $M_{\rm NSC}$ (Neumayer et al., 2020; Pechetti et al.,
2020). To capture this trend, we generate $\gamma_{\rm NSC}$ from uniform
distributions, whose upper bounds depend on $M_{\rm NSC}$. To be specific, we
have $\gamma_{\rm NSC}\sim 0.65-2.55$ for $M_{\rm NSC}<M_{\rm crit}$,
$\gamma_{\rm NSC}\sim 0.65-1.35$ for $M_{\rm NSC}>10^{8}\ \rm M_{\odot}$, and
the upper bound evolves linearly with $\log M_{\rm NSC}$ in-between. Fig. 9
shows the predictions of our stochastic model of NSC size-mass relations for
both $R_{\rm eff}$ (top) and $R_{\rm NSC}$ (bottom). For $R_{\rm eff}$, we
show the observational data complied by Neumayer et al. (2020). For $R_{\rm
NSC}$, the results for observed NSCs are derived from Equation (9) with
randomly generated $\gamma_{\rm NSC}$. We further illustrate the size
distribution by a mock sample of NSCs with a log-flat mass function in the
range of $M_{\rm NSC}\sim 10^{5}-10^{9}\ \rm M_{\odot}$. Similar to the sample
of host NSCs produced by our merger trees, for each NSC in the mock sample,
both $R_{\rm eff}$ and $\gamma_{\rm NSC}$ are drawn randomly to determine
$R_{\rm NSC}$. It turns out that $R_{\rm NSC}$ slowly increases with $M_{\rm
NSC}$, and $R_{\rm NSC}\sim 1-20\ \rm pc$ for most NSCs.

Figure 9: Size-mass relations for NSCs in terms of the effective radius (top)
and scale length (bottom). Observational data compiled by Neumayer et al.
(2020, see their fig. 7) are shown with circles and triangles for early-type
and late-type galaxies, respectively. The MW is labelled with the green star
for $M_{\rm NSC}=2.5\pm 1.8\times 10^{7}\ \rm M_{\odot}$, $R_{\rm eff}=6.5\pm
2.7\ \rm pc$ (Neumayer et al., 2020) and $R_{\rm NSC}\sim 4.4\ \rm pc$. In the
top panel, scatters of 0.55 dex are denoted by the shaded region. In the
bottom panel, typical scatters of $R_{\rm NSC}$ are shown with gray dots for a
mock sample of NSCs with a log-flat mass function in the range of $M_{\rm
NSC}\sim 10^{5}-10^{9}\ \rm M_{\odot}$.

Finally, once the NSC properties ($M_{\rm NSC}$, $R_{\rm NSC}$ and
$\gamma_{\rm NSC}$) are given, we can calculate the evolution of an in-falling
compact object binary. Similar to the galaxy model, the remnant will sink into
the NSC centre by DF. We fix the eccentricity $e_{\rm if}$ of the in-fall
orbit 0 because the orbit will be mostly circularized by DF in the galaxy when
the object reaches the NSC. For simplicity, we assume that NSCs do not evolve
with time999Before compact object binaries merge, their host galaxy may fall
into larger structures and even be destroyed/stripped in a galaxy merger.
Nevertheless, NSCs can survive such mergers due to their compact nature and
appear as compact GCs and ultra compact dwarfs (UCDs). In this work, we assume
a 100% survive rate of NSCs in galaxy mergers for simplicity., such that
Equation (7) has an analytical solution:

$\displaystyle r(t)\simeq r_{0}\left(1-\beta\frac{t}{\tau_{\rm
DF,0}}\right)^{1/\beta}\ ,\quad t<\hat{t}_{\rm DF}\equiv\tau_{\rm DF,0}/\beta\
,$ (10)

where the sink time $\hat{t}_{\rm DF}$ is the expected time taken for the
binary to reach the cluster centre101010A small fraction of highly eccentric
binaries can merge before reaching the centre/core., $\beta=1.76$ and
$\tau_{\rm DF,0}=\tau_{\rm DF}(r_{0}|m_{\bullet},M_{\rm NSC},R_{\rm
NSC},\gamma_{\rm NSC},e_{\rm if}=0)$, given the initial distance $r_{0}=3\ \rm
pc$. As the binary moves inside the NSC, the stellar environment around it
evolves with time according to $r(t)$, with the stellar density
$\rho_{\star}(r)$ given by the Dehnen profile (Equ. 2), and the velocity
dispersion following (Dehnen, 1993; Sedda, 2020)

$\displaystyle\sigma_{\star}(r)\simeq\sqrt{\frac{GM_{\rm NSC}}{R_{\rm
NSC}}}(r/R_{\rm NSC})^{\delta/2}\ ,\quad r\lesssim R\ ,$ (11)

where $\delta=\gamma_{\rm NSC}$ for $\gamma_{\rm NSC}<1$ and
$\delta=2-\gamma_{\rm NSC}$ for $\gamma_{\rm NSC}\geq 1$. Since the Dehnen
profile can be divergent at the centre, we assume that each NSC has a core
with a radius $r_{\rm c}=R_{\rm NSC}/\Delta_{\rm c}^{1/\gamma_{\rm NSC}}$,
defined by the core overdensity parameter $\Delta_{\rm c}$. Once the binary
enters the core ($r<r_{\rm c}$), we fix the environment to that evaluated at
$r_{\rm c}$. The overdensity parameter $\Delta_{\rm c}$ determines the rate of
DH, as most binary evolution happens in the core with high density. We
consider two cases with $\Delta_{\rm c}=100$ and 20.

Within certain stellar environments, hard binaries will be further hardened by
three-body encounters and GW emission until the final mergers111111For
simplicity, we ignore higher-order effects such as exchanges and repeated
mergers. Note that in our case the dominant ($\gtrsim 90$%) hosts of compact
object mergers are low-mass NSCs ($M_{\rm NSC}\lesssim 10^{6}\ \rm M_{\odot}$,
see Sec. 3.1) which may not be able to retain second-generation BHs
efficiently. We also defer the inclusion of central massive BHs to future work
(see e.g. Sedda 2020)., while soft binaries will be disrupted/destroyed. For
each binary at galaxy centres, we draw the initial separation $a_{0}$ and
eccentricity $e_{0}$ from the catalog of binaries produced by the N-body
simulations in Liu et al. (2021), matching the total mass $m_{\bullet}$. We
only keep track of hard binaries whose initial separations are below the
critical separation $a_{\rm
HDB}=Gm_{1}m_{2}/\left(m_{\star}\sigma_{\star}^{2}\right)$ (Mapelli et al.,
2021), where we evaluate the stellar environment at $r=R_{\rm NSC}$, since
$\sigma_{\star}$ peaks around $R_{\rm NSC}$ and $\rho_{\star}$ drops rapidly
at $r\gtrsim R_{\rm NSC}$. We also require $a_{0}<r_{\rm c}$ for efficient DH.
The evolution of binary parameters driven by three-body encounters
($|_{\star}$) and emission of GWs ($|_{\rm GW}$) is described by the following
differential equations:

$\displaystyle\begin{split}\frac{da}{dt}=\left.\frac{da}{dt}\right|_{\rm\star}+\left.\frac{da}{dt}\right|_{\rm
GW}\
,\quad\frac{de}{dt}=\left.\frac{de}{dt}\right|_{\star}+\left.\frac{de}{dt}\right|_{\rm
GW}\ ,\end{split}$ (12)

in which the three-body terms are (Sesana & Khan, 2015; Mapelli et al., 2021;
Sedda, 2020)

$\displaystyle\begin{split}\left.\frac{da}{dt}\right|_{\rm\star}=-\frac{GH\rho_{\star}(r)}{\sigma_{\star}(r)}a^{2}\
,\quad\left.\frac{de}{dt}\right|_{\star}=\kappa\frac{GH\rho_{\star}(r)}{\sigma_{\star}(r)}a\
,\end{split}$ (13)

where $H\sim 1-20$ and $\kappa\sim 0.01-0.1$ (Sesana et al., 2006; Sesana &
Khan, 2015; Mapelli et al., 2021) are two dimensionless factors. We adopt
$H=20$ and explore two cases with $\kappa=0.01$ and 0.1. The GW terms can be
written as (Peters, 1964; Mapelli et al., 2021)

$\displaystyle\begin{split}\left.\frac{da}{dt}\right|_{\rm
GW}&=-\frac{64}{5}\frac{A}{a^{3}(1-e^{2})^{7/2}}f_{1}(e)\ ,\\\
\left.\frac{de}{dt}\right|_{\rm
GW}&=-\frac{304}{15}e\frac{A}{a^{4}(1-e^{2})^{5/2}}f_{2}(e)\ ,\end{split}$
(14)

where $A=G^{3}m_{1}m_{2}(m_{1}+m_{2})/c^{5}$, and

$\displaystyle\begin{split}f_{1}(e)=1+\frac{73}{24}e^{2}+\frac{37}{96}e^{4}\
,\quad f_{2}(e)=1+\frac{121}{304}e^{2}\ .\end{split}$ (15)

For each binary, we integrate Equations (12) with the 4th-order Runge-Kutta
method for at most 15 Gyr and stop the integration if a merger happens with
$a<6G(m_{1}+m_{2})/c^{2}$. Note that in this calculation, we ignore the
possibility of ejection by three-body encounters. To take this into account,
following Mapelli et al. (2021), we compare the separation below which the
evolution is dominated by GW emission (Mapelli et al., 2021)

$\displaystyle a_{\rm
GW}=\left[\frac{64}{5H}\frac{A\sigma_{\star}f_{1}(e)}{G\rho_{\star}(1-e^{2})^{7/2}}\right]^{1/5}$
(16)

with the critical separation below which the binary, if not dominated by GW
emission, will be ejected by three-body encounters (Mapelli et al.,
2021)121212By our definition, the numerical factor $H$ is related to the $\xi$
in the formalism of Mapelli et al. (2021, see their equ. 11 and 17) by
$H=2\pi\xi$.

$\displaystyle a_{\rm
ej}=\frac{Hm_{\star}^{2}}{\pi(m_{1}+m_{2})^{3}}\frac{Gm_{1}m_{2}}{v_{\rm
esc}^{2}}\ ,$ (17)

where the escape velocity is given by (Fragione & Silk, 2020)

$\displaystyle\frac{v_{\rm esc}}{\mathrm{km\ s^{-1}}}=40\left(\frac{M_{\rm
NSC}}{10^{5}\ \rm M_{\odot}}\right)^{1/3}\left(\frac{\rho_{\star}}{10^{5}\ \rm
M_{\odot}\ pc^{-3}}\right)^{1/6}\ .$ (18)

Here we evaluate $\rho_{\star}$ at $r=R_{\rm NSC}$ as the typical stellar
density. If $a_{\rm gw}>a_{\rm ej}$ always holds, as in massive NSCs ($M_{\rm
NSC}\gtrsim 10^{6}\ \rm M_{\odot}$), the binary will not be ejected before
merger. If $a<a_{\rm ej}<a_{\rm GW}$ occurs, the binary will be ejected at
that moment ($t_{\rm ej}$). The subsequent evolution is no longer affected by
three-body encounters and the time to merger for such isolated binaries driven
by GW emission is

$\displaystyle\tau_{\rm
GW}=f_{3}(e)\frac{5}{256}\frac{a^{4}(1-e^{2})^{7/2}}{Af_{1}(e)}\ ,$ (19)

where $f_{3}(e)$ is a factor that captures the circularization of binary
orbits by GW emission, and $a=a(t_{\rm ej})$, $e=e(t_{\rm ej})$ for the
ejected binary. We adopt $f_{3}(e)=(1-0.8e)^{-1}$ as an approximation to the
result of numerical integration for simplicity131313See equ. 24 in Liu & Bromm
(2020a) for the exact formula for $\tau_{\rm GW}$.. The delay time $t_{\rm
delay}$ is defined as the time between NSC in-fall, $t_{\rm if}$, and the
final merger (or end of integration), $t_{\rm mg}$. With this formalism, the
population of Pop III binary remnants in NSCs is mapped to the population of
GW events of compact object mergers via $t_{\rm mg}=t_{\rm if}+t_{\rm delay}$.

Figure 10: Evolution of binary separation for a NS-BH binary with $m_{1}=20\
\rm M_{\odot}$, $m_{2}=2\ \rm M_{\odot}$, $a_{0}=100\ \rm au$ and $e_{0}=0.1$,
in a low-mass NSC with $M_{\rm NSC}=10^{5}\ \rm M_{\odot}$, $R_{\rm NSC}=3\
\rm pc$ and $\gamma_{\rm NSC}=1.5$ (for $\Delta_{\rm c}=100$), assuming no
ejection. The $a(t)$ curve is color coded by $\tau_{\rm GW}$. The critical
separations for GW-dominance and ejection by three-body encounters are shown
with the dashed and dashed-dotted curves, respectively. The critical
separation for hard binaries, reduced by a factor of 100, is shown with the
dotted curve. The sink time is $\hat{t}_{\rm DF}\simeq 60\ \rm Myr$. According
to the criterion, $a<a_{\rm ej}<a_{\rm GW}$, the binary is expected to be
ejected at $t\sim 3~{}\ \rm Gyr$ and merge at $t\sim 4.7\ \rm Gyr$. If not
ejected, it will merge at $t\sim 4.1\ \rm Gyr$. Figure 11: Same as Fig. 10 but
for a BH-BH binary similar to that observed in the event GW190521 ($m_{1}=85\
\rm M_{\odot}$, $m_{2}=66\ \rm M_{\odot}$, $a_{0}=200\ \rm au$, $e_{0}=0$) in
a MW-like NSC ($M_{\rm NSC}=2.5\times 10^{7}\ \rm M_{\odot}$, $R_{\rm
NSC}=4.4\ \rm pc$, $\gamma_{\rm NSC}=1.8$) with a core density
$\rho_{\star}\sim 2.5\times 10^{7}\ \rm M_{\odot}\ pc^{-3}$ at $r\lesssim
0.01\ \rm pc$ (Schödel et al., 2018). Here $a_{\rm ej}$ is much smaller than
$a_{\rm GW}$ and thus not shown. The sink time is $\hat{t}_{\rm DF}\simeq 100\
\rm Myr$, and the merge time is $t_{\rm mg}\simeq 170\ \rm Myr$.

We end this section with two (extreme) examples for the evolution of Pop III
binary remnants in NSCs under $\kappa=0.1$. In the first case, we consider a
NS-BH binary with $m_{1}=20\ \rm M_{\odot}$, $m_{2}=2\ \rm M_{\odot}$,
$a_{0}=100\ \rm au$ and $e_{0}=0.1$, in a low-mass NSC with $M_{\rm
NSC}=10^{5}\ \rm M_{\odot}$, $R_{\rm NSC}=3\ \rm pc$ and $\gamma_{\rm
NSC}=1.5$ (for $\Delta_{\rm c}=100$). As shown in Fig. 10 (assuming no
ejection), the binary sinks to the cluster core after $\hat{t}_{\rm DF}\simeq
60\ \rm Myr$, and is ejected at $t_{\rm ej}\sim 3\ \rm Gyr$. The delay time
with(out) ejection is $t_{\rm delay}\simeq 4.7\ (4.1)\ \rm Gyr$. In the second
case, we consider a BH-BH binary similar to that observed in the event
GW190521_030229 (GW190521 for short) with $m_{1}=85\ \rm M_{\odot}$,
$m_{2}=66\ \rm M_{\odot}$ (Abbott et al., 2020a), $a_{0}=200\ \rm au$ and
$e_{0}=0$, in a MW-like NSC ($M_{\rm NSC}=2.5\times 10^{7}\ \rm M_{\odot}$,
$R_{\rm NSC}=4.4\ \rm pc$, $\gamma_{\rm NSC}=1.8$) with $\Delta_{\rm c}=1000$.
Such a massive binary in a massive dense NSC will not be ejected before
merger, for which we have $\hat{t}_{\rm DF}\simeq 100\ \rm Myr$ and $t_{\rm
delay}\sim 170\ \rm Myr$, as shown in Fig. 11. Here $\Delta_{\rm c}=1000$ is
chosen to be consistent with observations of the realistic MW NSC by Schödel
et al. (2018), which show $\rho_{\star}\sim 2.6\pm 0.3\times 10^{7}\ \rm
M_{\odot}\ pc^{-3}$ at $r\lesssim 0.01\ \rm pc$.

##  3 Model predictions

In this section, we present the demography, rates and host properties of Pop
III binary remnants falling into NSCs and the subsequent mergers. We explore
two sets of parameters for three-body encounters: (i) the optimistic (OP) case
for efficient dynamical and GW hardening with $\Delta_{\rm c}=100$,
$\kappa=0.1$, and (ii) the conservative (CS) case for moderate hardening with
$\Delta_{\rm c}=20$, $\kappa=0.01$. We also consider two assumptions for the
occupation frequency of NSCs (see Sec. 2.3): (i) full occupation (F), i.e.
$f_{\rm NSC}=1$, in galaxies with $M_{\rm NSC}>10^{5}\ \rm M_{\odot}$, and
(ii) partial occupation (P) with $f_{\rm NSC}(M_{\star})<1$, based on local
observations (see Fig. 8). Combining these parameter choices leads to 4 NSC
models: OP_F, OP_P, CS_F, and CS_P.

Fig. 12 shows the in-fall rate density (IFRD) of Pop III binary remnants at
galaxy centres ($r<3\ \rm pc$) for the entire sample and those merged at
$z>0$, on top of the formation rate density of Pop III binary remnants. It
turns out that the in-fall of Pop III binary remnants into galaxy centres
starts at $z\sim 15-17$, and the IFRD peaks around $z\sim 6-10$, with a shape
that is almost independent of NSC models. Note that most Pop III SF (as well
as remnant formation) happens at $z\sim 6-20$ with 2 peaks at $z\sim 8$ and
14. The delay between remnant formation and in-fall reflects the time taken
for (Pop I/II) galaxies/NSCs to form and grow, and DF to drive the in-spiral
of remnants.

In Sec. 3.1, we discuss the host galaxy and NSC properties for the entire
sample of Pop III binary remnants at galaxy centres as well as those merged at
$z>0$. In Sec. 3.2, we focus on the mergers of Pop III remnants in NSCs,
predicting their statistics and cosmic MRD.

Figure 12: IFRDs of Pop III binary remnants at galaxy centres ($r<3\ \rm pc$)
for the entire sample (shaded region) and those merged at $z>0$ (solid
curves). The upper and lower bounds of the former denote the models with full
occupation (OP/CS_F) and observation-based occupation frequency (OP/CS_P) of
NSCs. While for the latter, the upper and lower bounds correspond to OP_F and
CS_P. To illustrate the delay between remnant formation and in-fall into
galaxy centres, we also show the formation rate density of Pop III binary
remnants with the dashed curve, which can be expressed in terms of the SFRD as
$\dot{\rho}_{\bullet,\rm form}=f_{\rm B}f_{\rm rem}\rm SFRD_{\rm III}$, where
$f_{\rm B}=0.69$ is the mass fraction of binary stars, and $f_{\rm rem}=0.31$
is the overall ratio of remnant mass to initial stellar mass (in binaries).

###  3.1 Host systems of Pop III remnants

Figure 13: Distributions of NSC masses hosting Pop III binary remnants for the
OP models, for the entire sample of NSCs and those merged at $z>0$, with full
occupation (histograms and the solid contour) and observation-based occupation
frequency (dashed and dashed-dotted contours). It is evident that low-mass
NSCs are the dominant hosts, such that $90$% of remnants merged at $z>0$ are
from NSCs with $M_{\rm NSC}\lesssim 3.9\ (8.1)\times 10^{5}\ \rm M_{\odot}$ in
the OP_F (P) model. Figure 14: Distributions of galaxy (stellar) masses
containing NSCs that host Pop III binary remnants from the OP models, for the
entire sample and those merged at $z>0$ with full occupation (histograms and
the solid contour) and observation-based occupation frequency (dashed and
dashed-dotted contours). As can be seen, Pop III remnants mostly fall into the
centres of low-mass dwarf galaxies, such that $90$% of remnants merged at
$z>0$ are from galaxies with $M_{\star}\lesssim 7.5\ (34)\times 10^{6}\ \rm
M_{\odot}$ in the OP_F (P) model. Figure 15: Sink time distributions of Pop
III binary remnants in NSCs from the OP models, for the entire sample in NSCs
and those merged at $z>0$ with full occupation (histograms and the solid
contour) and observation-based occupation frequency (dashed and dashed-dotted
contours) of NSCs.

For conciseness, here we only show the results from the OP models for
different NSC occupation assumptions, as the trends seen in the CS models are
similar. Fig. 13 and 14 show the distributions (weighted by the number of
binary remnants) of NSC masses hosting Pop III binary remnants, and of the
corresponding galaxy stellar masses (at the moment when the binary reaches the
centre), respectively. In all cases, the distributions are significantly
biased towards low-mass systems, especially for the binaries merged at $z>0$.
The reason is that DF is more efficient in low-mass galaxies (see Fig. 5), and
Pop III binary remnants are more likely to be destroyed in more massive NSCs
with higher $\sigma_{\star}$. Actually, assuming full occupation of NSCs, $50\
(90)\%$ of Pop III binary remnant mergers are hosted by NSCs with $M_{\rm
NSC}\lesssim 1.2\ (3.9)\times 10^{5}\ \rm M_{\odot}$, (initially) in galaxies
with $M_{\star}\lesssim 4.5\ (75.0)\times 10^{5}\ \rm M_{\odot}$. With the
observation-based occupation frequency, the mass distributions are less
dominated by small systems, such that we have $M_{\rm NSC}\lesssim 1.4\
(8.1)\times 10^{5}\ \rm M_{\odot}$ and $M_{\star}\lesssim 7.0\ (340)\times
10^{5}\ \rm M_{\odot}$ for $50\ (90)\%$ of mergers.

Note that the dominant (initial) host galaxies of our Pop III binary remnant
mergers are high-redshift ($z\gtrsim 4$) dwarf galaxies ($M_{\star}\lesssim
10^{7}\ \rm M_{\odot}$) which will further grow and merge into larger
galaxies. Some of them will even be destroyed/stripped in galaxy mergers where
the NSCs can survive as compact GCs and ultra compact dwarfs (UCDs). As the
final merger can be significantly delayed with respect to the in-fall (see
below), we expect a more diverse population of host systems at the moments of
merger, including NSCs of more massive central galaxies and satellites, GCs
and UCDs, compared with the initial hosts at the moments of in-fall. We defer
a detailed investigation into the host systems at the moments of merger to
future work.

Fig. 15 shows the distributions of sink time $\hat{t}_{\rm DF}$ (see Equ. 10
for definition) of Pop III binary remnants in NSCs. In general, we have
$\hat{t}_{\rm DF}\sim 5-10^{3}\ \rm Myr$, and the distribution peaks at
$\hat{t}_{\rm DF}\sim 20\ \rm Myr$. The results from different NSC occupation
models are similar, although the distribution is narrower in the observation-
based model (with $\hat{t}_{\rm DF}\lesssim 500\ \rm Myr$) than that with full
occupation. When restricted to the binaries merged at $z>0$, cases with
$\hat{t}_{\rm DF}\gtrsim 200\ \rm Myr$ are reduced, especially for the full
occupation model. We expect these cases to be extreme low-mass, low-density
NSCs where three-body hardening is inefficient. This reduction is more
efficient in the CS models.

###  3.2 Mergers of Pop III remnants in NSCs

Table 1: Key characteristics of Pop III binary remnant mergers for 4 NSC models, where $\Delta_{\rm c}$ is the core overdensity parameter that measures how dense the NSC core is compared with the average, $\kappa$ is the dimensionless efficiency at which DH increases the binary eccentricity (see Equ. 13), $f_{\rm NSC}$ is the NSC occupation frequency (see Fig. 8 for the case of $f_{\rm NSC}<1$), and $f_{\rm mg}$ is the fraction of Pop III binary remnants in galaxy centres ($r<3\ \rm pc$) that eventually merge at $z>0$. The peak MRD and the corresponding peak redshift $z_{\rm peak}$ are also shown (see Fig. 19 for details). Note that on average $760$ Pop III binary remnants will fall into galaxy centres out of the $2.1\times 10^{5}\ \rm M_{\odot}$ of Pop III stars formed in a merger tree, the efficiency of GW events $\epsilon_{\rm GW}$, defined as the number of compact object mergers (at $z>0$) per unit stellar mass, can be related to the merge fraction $f_{\rm mg}$ with $\epsilon_{\rm GW}\simeq 0.0036f_{\rm mg}\ \rm M_{\odot}^{-1}$. Model | $\Delta_{\rm c}$ | $\kappa$ | $f_{\rm NSC}$ | $f_{\rm mg}$ | Peak MRD | $z_{\rm peak}$
---|---|---|---|---|---|---
|  |  |  |  | $[\rm yr^{-1}\ Gpc^{-3}]$ |
OP_F | 100 | 0.1 | 1 | 0.50 | 10.4 | 7.2
OP_P | 100 | 0.1 | $<1$ | 0.090 | 1.78 | 7.2
CS_F | 20 | 0.01 | 1 | 0.30 | 2.28 | 5.2
CS_P | 20 | 0.01 | $<1$ | 0.055 | 0.378 | 5.2

In Table 1, we summarize the key characteristics of Pop III binary remnant
mergers (in NSCs) for different NSC models. As mentioned in Sec. 2.2, $7.5\%$
of Pop III remnants will fall into the centres ($r<3\ \rm pc$) of galaxies
with $M_{\rm NSC}>10^{5}\ \rm M_{\odot}$, of which $\sim 5.5-50\%$ can merge
at $z>0$. The corresponding efficiency of GW events $\epsilon_{\rm GW}$,
defined as the number of compact object mergers (at $z>0$) per unit stellar
mass, takes the range of $0.2-2\times 10^{-4}\ \rm M_{\odot}$, which is
comparable to the values in the classical binary stellar evolution (BSE)
channel (Kinugawa et al., 2014; Hartwig et al., 2016; Belczynski et al., 2017;
Hijikawa et al., 2021). We also find that ejection of binaries by three-body
encounters is rare ($\sim 1\%$), independent of NSC models.

####  3.2.1 Delay time and binary statistics

Figure 16: Delay time distributions of Pop III hard binary remnants in NSCs
for the OP (histograms) and CS (solid contour) scenarios, with full occupation
of NSCs. Power-law approximations for the distributions at $t_{\rm
delay}\gtrsim 20\ \rm Myr$ are shown with the dotted lines. The tail at $t<20\
\rm Myr$ makes up $\simeq 7.6\ (0.6)\%$ of all mergers in the OP (CS)_F model.
Here, the binaries with $t_{\rm delay}>15\ \rm Gyr$ are pilled up at the last
bin of $t_{\rm delay}\sim 15\ \rm Gyr$.  Figure 17: Distributions of initial
separation of Pop III binary remnants in NSCs from the OP_F model, for the
entire sample (histograms) and those merged at $z>0$ (solid contour). The sub-
samples of NS-BH and NS-NS mergers are also shown with the dashed and dashed-
dotted contours, accounting for 6.7% and 0.014% of all mergers, respectively.
The thin vertical line denotes the typical maximum radius of massive Pop III
stars, $3000\ \rm R_{\odot}$, which is approximately the critical separation
for close binary interactions (such as common envelope evolution). Figure 18:
Cumulative distributions of initial eccentricity and mass ratio of Pop III
binary remnants in NSCs from the OP_F model, for the entire sample (thick
solid and dashed) and those merged at $z>0$ (solid and dashed-dotted). The
dotted line denotes the case of a uniform distribution in [0,1]. The thin
horizontal line is plotted to point out the medians. The mass ratio
distribution for binaries merged at $z>0$ is more biased to the high end than
that for the entire sample since low-$q$ binaries are less massive and thus
merge less efficiently in NSCs.

Now we consider the distributions of delay times141414In our case, the delay
time is defined as the time between NSC in-fall and the final merger, instead
of the time between the initial remnant/star formation and merger. Therefore,
the MRD is a convolution of the delay time distribution with the IFRD (see
Fig. 12) rather than the SFRD. and binary orbital parameters for our Pop III
binary remnants. Fig. 16 shows the delay time distributions of hard151515On
average 60% of our Pop III binary remnants in NSCs are hard binaries. Pop III
binary remnants for the OP and CS models with full occupation of NSCs. The
results with the observation-based occupation frequency of NSCs are similar
and not shown. Among these hard binaries, $\sim 90\ (50)\%$ can merge at $z>0$
under the OP (CS) hardening. In both cases, the delay time distribution is
almost log-flat at $t_{\rm delay}\gtrsim 20\ \rm Myr$, while $t_{\rm delay}$
is generally shorter in OP models than in CS models. When approximated with a
power law, $dN/d\log t_{\rm delay}\propto t_{\rm delay}^{\zeta}$, the
distribution has a power-law slope $\zeta\sim-0.2$ with the OP hardening and
$\zeta\sim 0.3$ for the CS case. Moreover, there is a tail of lower delay
times extending down to $t_{\rm delay}\sim 0.1\ \rm Myr$, which is more
important with more efficient hardening. For instance, the fraction of hard
binaries with $t_{\rm delay}<20\ \rm Myr$ is $\simeq 0.076$ (0.006) for the OP
(CS) models. We expect this tail to be made of binaries with high initial
eccentricities ($e_{0}\gtrsim 0.9$) that can merge even before reaching the
NSC core due to hardening by GW emission (see Equs. 14).

For the statistics of binary orbital parameters of our Pop III binary remnants
in NSCs, we only show the results of the OP_F model since similar trends are
found in other NSC models. Fig. 17 shows the distributions of initial
separation $a_{0}$ of Pop III binary remnants in NSCs (including sub-samples
of NS-NS161616We assume that remnants with masses less than $3\ \rm M_{\odot}$
are NSs. and NS-BH mergers). For the entire sample in NSCs, the distribution
is very broad, ranging from a few au to a few pc, which is dominated ($\gtrsim
99.7\%$) by wide binaries ($a_{0}\gtrsim 10\ \rm au$)171717As discussed in Liu
et al. (2021), Pop III star clusters expand during the protostellar accretion
phase due to angular momentum conservation, making it difficult to form close
binaries of Pop III stars/remnants.. For such wide-binaries, close binary
interactions are likely inactive such that they cannot merge within a Hubble
time if isolated. However, with DH in NSCs, binaries with initial separations
up to $2\times 10^{5}\ \rm au$ can still merge at $z>0$, making up about half
of the entire sample in NSCs. These mergers are dominated by BH-BH binaries
($\gtrsim 90\%$), while NS-BH and NS-NS mergers only count for $\sim 4-7\%$
and $\lesssim 0.014\%$ of all GW events. NS-BH and NS-NS mergers also have
lower initial separations ($a_{0}\lesssim 6000$ and 40 au, respectively),
compared with the dominant BH-BH mergers. This is due to the fact that in
young Pop III star clusters, low-mass binaries tend to be disrupted by
encounters with massive stars, leading to lower initial fractions of binaries
involving NSs. Besides, such binaries with relatively lower masses can be more
easily destroyed in massive NSCs, and the hardening from GW emission is also
less efficient.

Fig. 18 displays the cumulative distributions of initial eccentricity $e_{0}$
and mass ratio $q$ of Pop III binary remnants in NSCs. The distribution of
$e_{0}$ is almost uniform for the entire sample, consistent with the results
from N-body simulations (Liu et al., 2021). For binaries merged at $z>0$, on
the other hand, the distribution is slightly biased towards higher $e_{0}$
since hardening by GW emission is more efficient for eccentric binaries (see
Equs. 14 and 19). This bias is slightly more pronounced for the CS models,
such that the median initial eccentricity of mergers is 0.57 (0.60) with OP
(CS) hardening. Interestingly, we find that the initial eccentricity does not
correlate with $a_{0}$ and $q$.

The mass ratio distribution is more complex, with three groups of binaries at
$q\sim$ (a) $0.02-0.07$, (b) $0.07-0.5$ and (c) $0.5-1$. Each regime can be
approximated with a uniform distribution. For Group (a) binaries, the
primaries are mostly direct-collapse BHs with $m_{1}\gtrsim 25\ \rm
M_{\odot}$, and the secondaries are NSs with $m_{2}\sim 1.4-3\ \rm M_{\odot}$.
In Group (b), most binaries are made of two BHs, one with $m_{2}\sim 5-40\ \rm
M_{\odot}$ including those (with $m_{2}\sim 5-25\ \rm M_{\odot}$) born in
core-collapse supernovae, and another with $m_{1}\sim 40-85\ \rm M_{\odot}$,
including those (with $m_{1}\sim 40-57\ \rm M_{\odot}$) from PPISNe. Finally,
Group (c) binaries are mostly made of massive BHs with $m_{1}\sim m_{2}\sim
40-85\ \rm M_{\odot}$. When restricted to the binaries merged at $z>0$, the
median mass ratio shifts from $0.6$ (for the entire sample in NSCs) to 0.7.
The latter is actually consistent with the median mass ratio of GW events
detected by the LIGO collaboration (Abbott et al., 2020c; Kinugawa et al.,
2021b). Meanwhile, the fraction of Group (a) binaries is significantly reduced
(from $\sim 20\%$ to $\lesssim 10\%$), consistent with the trend seen in the
distribution of $a_{0}$, since it is more difficult for lower-mass binaries to
merge in NSCs.

####  3.2.2 Comparison with the BSE channel and LIGO events

Next, we compare our predictions for Pop III binary remnant mergers in NSCs
with the results for the BSE channel (Kinugawa et al., 2014; Hartwig et al.,
2016; Belczynski et al., 2017; Hijikawa et al., 2021) and observations by LIGO
(Abbott et al., 2020c, a; Nitz et al., 2021). In Fig. 19, we show the MRD
redshift evolution for our NSC models. Since our merger trees are only
representative at $z\gtrsim 5$, the low-$z$ rate can be underestimated,
especially for the OP models with shorter delay times. We extrapolate the MRDs
from $z\sim 5$ to $z=0$ in the OP models, to approximate the evolution in the
low-$z$ regime. For each work considered for comparison (Kinugawa et al.,
2014; Hartwig et al., 2016; Belczynski et al., 2017; Hijikawa et al., 2021),
we plot the lower and upper bounds for their different assumptions. As
mentioned in Sec. 2.1.1, we here constrain the Pop III SF history by scaling
the SFRDs (and MRDs) of these studies such that the scaled integrated Pop III
stellar mass densities ($\mathrm{ISMD}_{\rm
III}\equiv\int_{0}^{\infty}\mathrm{SFRD}_{\rm III}|dt/dz|dz$) are identical to
our value $\simeq 7\times 10^{4}\ \mathrm{M_{\odot}\ Mpc^{-3}}$. In general,
the MRD of Pop III binary remnants in NSCs is comparable to those for the BSE
channel, especially in the post-reionization era ($z\lesssim 6$). For our DH
channel with NSCs, the local ($z\sim 0$) MRD is $\sim 0.02-0.6\ \rm yr^{-1}\
Gpc^{-3}$, within the range $\sim 0.01-1\ \rm yr^{-1}\ Gpc^{-3}$ predicted by
most BSE models. However, our MRD only starts at $z\sim 15$ and peaks at a
later time ($z\sim 6$) compared with most BSE MRDs that extend to very high
redshifts ($z\sim 20$) with non-negligible rates ($\gtrsim 0.01\ \rm yr^{-1}\
Gpc^{-3}$). The reason is that galaxy/NSC formation/growth and in-spiral of
remnants by DF take time. Measuring the MRD to very high redshifts ($z\gtrsim
15$) can reveal the relative importance of the two channels.

Figure 19: MRDs of Pop III binary remnants in NSCs for our OP (shaded region)
and CS (solid) models with upper and lower bounds corresponding to full
occupation and observation-based frequency of NSCs. The extrapolations of the
MRDs from the OP models are shown with the dotted lines. For comparison, we
show the results of 4 studies for the BSE channel, which are scaled according
to our $\rm ISMD_{\rm III}$: Kinugawa et al. (2014, KT14, dashed), Hartwig et
al. (2016, HT16, dashed-dotted), Belczynski et al. (2017, BK17, sparsely-
dashed) and Hijikawa et al. (2021, HK21, long-dashed). The local rate density
inferred from the LIGO events, $\dot{n}_{\rm GW}(z\sim
0.33)=23.9_{-8.6}^{+14.3}\ \rm yr^{-1}\ Gpc^{-3}$ (Abbott et al., 2020c), is
shown with the (red-crossed) symbol (for 90% confidence interval). The thin
vertical line marks the turn-over redshift of the target halo, $z\simeq 4.6$,
below which our results may not be cosmologically representative.

In all the models for Pop III mergers considered here (see Fig. 19), with the
Pop III SF history constrained by the Planck measurement of the Thomson
optical depth, the local MRD is lower than that of LIGO events, $\dot{n}_{\rm
GW}(z\sim 0.33)=23.9_{-8.6}^{+14.3}\ \rm yr^{-1}\ Gpc^{-3}$ (Abbott et al.,
2020c), by at least one order of magnitude. Given our predicted local MRD
$\sim 0.02-0.6\ \rm yr^{-1}\ Gpc^{-3}$, up to $\sim 4\%$ (2) of the 57 events
detected by LIGO so far are expected to originate from Pop III binary remnants
in NSCs. Note that our merger trees are only cosmologically representative at
$z\gtrsim 5$, such that the local MRD can be underestimated since we do not
fully take into account low-mass ($M_{\rm h}\lesssim 10^{12}\ \rm M_{\odot}$)
haloes at lower redshifts, especially for $z<2.5$. The actual fraction of LIGO
events from Pop III binary remnants in NSCs could thus be higher. Besides,
most of our predicted GW events from Pop III binary remnant mergers in NSCs
can be detected at promising rates by the 3rd generation of GW instruments
capable of reaching $z\gtrsim 10$, planned for the next decades, such as the
Einstein Telescope (Hild et al., 2009) and Decihertz Observatory (Arca Sedda
et al., 2020). For instance, we predict a full-sky detection rate of $\sim
170-2720\ \rm yr^{-1}$ for a horizon redshift $z_{\max}=10$.

Figure 20: Mass distributions of Pop III binary remnant mergers in NSCs from
the OP_F model. Top panel: joint distribution of primary and secondary masses,
where the 57 events from the 3-OGC (Nitz et al., 2021) are plotted as data
points with 90% confidence intervals. We highlight with triangles the events
involving massive BHs that are most likely to have Pop III progenitors, as
well as the special event GW190814 that may be a NS-BH merger (see the main
text for details). Bottom pandel: distributions of chirp mass (histograms) and
total mass (solid contour), where the corresponding LIGO observations (Nitz et
al., 2021) are shown with the dashed and dashed-dotted contours, respectively.
All distributions are normalized by the total number of events. Here we use
the results from the NRSur7dq4 waveform model for the event GW190521 (Abbott
et al., 2020a). Note that the horizontal axis for the total mass is
$m_{\bullet}/2$.

Fig. 20 shows the joint distribution of primary and secondary masses $m_{1}$
and $m_{2}$ (top) as well as distributions of chirp mass $m_{\rm
chirp}\equiv(m_{1}m_{2})^{3/5}/(m_{1}+m_{2})^{1/5}$ and total mass
$m_{\bullet}\equiv m_{1}+m_{2}$ (bottom) for the Pop III binary remnant
mergers in the OP_F model, where the statistics of the 57 events from the the
third Open Gravitational-wave Catalog (Nitz et al., 2021, 3-OGC) are also
shown for comparison. We again use OP_F as the example since other models
predict similar distributions. Here, for the special event GW190521 that may
involve two BHs in the PPISN mass gap for Pop I/II stars, we show the results
from the NRSur7dq4 waveform model with $m_{1}=85_{-14}^{+21}\rm M_{\odot}$ and
$m_{2}=66_{-18}^{+17}\ \rm M_{\odot}$ (Abbott et al., 2020a). Most of the
mergers in our models are made of massive BHs with $m_{1}\sim m_{2}\sim 40-85\
\rm M_{\odot}$. Particularly, two regions in the $m_{2}$-$m_{1}$ space are
most densely populated, one with $m_{1}\sim m_{2}\sim 40-57\ \rm M_{\odot}$,
and another with $m_{1}\sim 65-85\ \rm M_{\odot}$ and $m_{2}\sim 40-57\ \rm
M_{\odot}$. The reason is that Pop III stars with initial masses
$m_{\star}\sim 40-57\ \rm M_{\odot}$ will collapse into BHs directly with no
supernovae, while stars with $m_{\star}\sim 85-115\ \rm M_{\odot}$ will also
form BHs in the same mass range from PPISNe, such that this mass range is
favored, especially for the secondary. We also predict a few percent of NS-BH
mergers with $m_{1}\sim 20-85\ \rm M_{\odot}$ and $m_{2}\sim 1.4-3\ \rm
M_{\odot}$. Our models predict a broad distribution of $m_{\rm chirp}$ with
two peaks. The dominant peak is at $m_{\rm chirp}\sim 50\ \rm M_{\odot}$, and
there is a smaller peak at $m_{\rm chirp}\sim 5\ \rm M_{\odot}$ corresponding
to NS-BH mergers. The majority of mergers in our models have $m_{\bullet}\sim
20-170\ \rm M_{\odot}$ with a peak in the total mass distribution at
$m_{\bullet}\sim 125\ \rm M_{\odot}$.

Figure 21: Primary mass distribution of Pop III binary remnant mergers from
the OP_F model (histograms), compared with the observed (solid) and reweighted
(dashed) primary mass distributions inferred from the 3-OGC (Nitz et al.,
2021, see their fig. 5).

Different from our Pop III binary remnant mergers, most LIGO events involve
relatively less massive BH-BH mergers with $m_{1}\lesssim 50\ \rm M_{\odot}$
and $m_{2}\lesssim 40\ \rm M_{\odot}$. The chirp (total) mass distribution
peaks at $\sim 30\ (60)\ \rm M_{\odot}$. This low mass regime is also
populated by a few percent of our mergers. Particularly, we have a small peak
at $m_{1}\sim 40\ \rm M_{\odot}$, $m_{2}\sim 25\ \rm M_{\odot}$ in the
$m_{2}$-$m_{1}$ map. Nevertheless, considering the low local MRD of Pop III
mergers, in terms of overall probability, if there are any LIGO events from
Pop III stars, they are likely to be the most massive ones, such as GW170403,
GW190519, GW190602, GW190701, GW190706 and GW190521 (highlighted with
triangles in the top panel of Fig. 20). Actually, the distribution of $m_{1}$
for LIGO events cannot be described by a single truncated power-law, given the
presence of features implying a high-mass component at $m_{1}\gtrsim 40\ \rm
M_{\odot}$ (Abbott et al., 2020c; Nitz et al., 2021; Roulet et al., 2021). Our
results indicate that Pop III remnants can make important contributions to
this high-mass regime, as shown in Fig. 21. Moreover, the special event
GW190814 with $m_{1}=23_{-2.4}^{+1.6}\ \rm M_{\odot}$ and
$m_{2}=2.6_{-0.1}^{+0.2}\ \rm M_{\odot}$, which can be a NS-BH merger, is also
covered by our population of Pop III NS-BH mergers.

Figure 22: MRDs of GW190521-like mergers from the OP (histograms) and CS
models (solid curves), compared with the MRD inferred from GW190521,
$\dot{n}_{\rm GW}(z\sim 1)=0.08_{-0.07}^{+0.19}\ \rm yr^{-1}\ Gpc^{-3}$ (The
LIGO Scientific Collaboration et al., 2021), denoted by the triangle with the
90% confidence interval. Similar to Fig. 19, extrapolations of the MRDs from
the OP models are shown with the dotted lines. The thin vertical line denotes
the turn-over redshift $z\simeq 4.6$ of the target halo below which our
results may not be cosmologically representative.

Finally, we focus on the special event GW190521 that involves unusually large
BH masses $85_{-14}^{+21}\ \rm M_{\odot}$ and $66_{-18}^{+17}\ \rm M_{\odot}$
(Abbott et al., 2020a, b; Mehta et al., 2021; O’Brien et al., 2021), within
the standard PPISN mass gap $\sim 55-130\ \rm M_{\odot}$ for Pop I/II stars
(see e.g. Heger et al. 2003; Belczynski et al. 2016; Woosley 2017; Marchant et
al. 2019). This event likely also falls within the high-mass ($m_{i}\gtrsim
50\ \rm M_{\odot}$), low-spin ($\chi_{i}\lesssim 0.2$) regime ($i=1$ or 2)
that cannot be easily populated by hierarchical mergers (Gerosa et al., 2021).
Several studies have shown that GW190521-like mergers can originate from Pop
III stars (Farrell et al., 2021; Tanikawa et al., 2021a; Kinugawa et al.,
2021a; Liu & Bromm, 2020c), since Pop III stars can retain most of their
hydrogen envelopes (due to their compactness and little mass loss), avoid the
PPISN regime, and form BHs in this mass range181818A similar mechanism can
also form BHs in the standard mass gap from metal-enriched stars with
metallicities up to $0.1\ \rm Z_{\odot}$ with enhanced core overshooting and
reduced mass loss (Vink et al., 2021)..

In light of this, we select the events similar to GW190521 from our population
of Pop III mergers in NSCs and calculate their MRDs for different NSC models.
Here we define GW190521-like mergers as those with $m_{1}\sim 71-106\ \rm
M_{\odot}$ and $m_{\bullet}\sim 133-179\ \rm M_{\odot}$, based on the 90%
confidence intervals from the NRSur7dq4 waveform model (Abbott et al., 2020a).
The results are shown in Fig. 22. We find a MRD range of $\sim 0.01-0.09\ \rm
yr^{-1}\ Gpc^{-3}$ at $z\sim 1$, overlapping with that inferred from GW190521,
$\dot{n}_{\rm GW}(z\sim 1)=0.08_{-0.07}^{+0.19}\ \rm yr^{-1}\ Gpc^{-3}$ (The
LIGO Scientific Collaboration et al., 2021), within the 90% confidence
interval. This confirms the conclusion in Liu & Bromm (2020c) that DH in dense
star clusters is an important channel of forming massive BH binary mergers
like GW190521 with BHs in the PPISN mass gap for Pop I/II stars.

##  4 Summary and discussions

We study the in-fall and evolution of Pop III binary (compact object) remnants
in NSCs in their full cosmological context, based on halo merger trees
equipped with self-consistent modelling for SF, stellar feedback, galaxy and
NSC properties, DF, hardening by three-body encounters, and GW emission. We
find that on average 7.5% of Pop III binary remnants will fall into galaxy
centres ($<3\ \rm pc$) by DF from field stars. Considering different
assumptions on the NSC occupation frequency and hardening efficiency, we then
construct the population of Pop III binary remnant mergers in NSCs and derive
their host system statistics, MRDs and mass spectra. We find that $5-50\%$ of
the remnants at galaxy centres can merge in NSCs at $z>0$. The dominant
($>90\%$) hosts of these mergers are low-mass ($M_{\rm NSC}\lesssim 10^{6}\
\rm M_{\odot}$) NSCs where the in-fall happens in high-redshift ($z\gtrsim
4.5$) dwarf galaxies with $M_{\star}\lesssim 3\times 10^{7}\ \rm M_{\odot}$.
Below we summarize the key features of Pop III binary remnant mergers formed
in this NSC-DH channel and discuss their implications, in comparison with the
results for the BSE channel (Kinugawa et al., 2014; Hartwig et al., 2016;
Belczynski et al., 2017; Hijikawa et al., 2021), and GW events detected by
LIGO (Abbott et al., 2020c, a; Nitz et al., 2021):

  * •

The MRDs from our NSC channel start to rise at $z\sim 13-15$ and peak at lower
redshifts $z\sim 5-7$ with comparable peak values $\sim 0.4-10\ \rm yr^{-1}\
Gpc^{-3}$ and local values $\sim 0.02-0.6\ \rm yr^{-1}\ Gpc^{-3}$ compared to
the results for the BSE channel, most of which peak at $z\sim 8-13$ and extend
to very high redshifts ($z\sim 20$). This difference in time evolution is
caused by the delay between initial formation of Pop III remnants and the
emergence of the galaxy/NSC systems required for in-spiral by DF to take
place. In this way, the relative importance of the NSC-DH and BSE channels can
be evaluated by measuring the MRD at high redshifts ($z\gtrsim 10$) with next-
generation instruments.

  * •

Unlike the BSE channel, which is only relevant for close binaries
($a_{0}\lesssim 10\ \rm au$), in our NSC-DH channel, the initial binary
separations are no longer important such that binaries with very large initial
separations (up to $1\rm\ pc$) can still merge in NSCs driven by DH. The NSC-
DH channel offers a fundamentally different way of enabling Pop III mergers
even if most Pop III binaries are initially wide as indicated by N-body
simulations of young Pop III star clusters (Liu et al., 2021).

  * •

Although insensitive to initial binary separations, the NSC-DH channel is
sensitive to the remnant masses and thus the Pop III IMF. Note that both DF
and GW emission are more efficient for more massive binaries such that the
fraction of remnants merged at $z>0$ in NSCs decreases from $\sim 7\%$ for
$m_{\bullet}\sim 100\ \rm M_{\odot}$ to $\sim 0.2\%$ for $m_{\bullet}\sim 10\
\rm M_{\odot}$, and the MRD of Pop III binary remnant mergers in NSCs is lower
when the IMF is less top-heavy. For instance, in this work, we consider a log-
flat Pop III IMF with an upper bound of $170\ \rm M_{\odot}$. Given a fixed
total stellar mass, if Pop III stars were less massive with an upper mass
limit of $20\ \rm M_{\odot}$, the MRD in the NSC-DH channel would be reduced
by a factor of 15, while that of the BSE channel would be increased by a
factor of 5.

  * •

Lacking close binary interactions, Pop III binary remnant mergers formed in
the NSC-DH channel will have different binary properties with respect to those
formed via BSE.

    * –

First, we expect more binaries with small mass ratios without the regulation
of binary mass transfer.

    * –

Second, the NSC-DH channel can produce eccentric mergers efficiently since
(initially) wide remnant binaries will not be pre-circularized by BSE, and
three-body encounters enhance eccentricity. For instance, under $\kappa=0.1$
and $\Delta_{\rm c}=100$, $\sim 0.1\%$ of our Pop III remnant binaries with
high initial eccentricities ($e_{0}\gtrsim 0.999$) have $e\gtrsim 0.1(f_{\rm
GW}/0.1\ \rm Hz)^{-1}$ when observed at a rest-frame frequency $f_{\rm GW}$,
given typical binary parameters $m_{1}\sim 70\ \rm M_{\odot}$, $m_{2}\sim 40\
\rm M_{\odot}$ and $a_{0}\sim 300\ \rm au$ in NSCs of $M_{\rm NSC}\sim 10^{6}\
\rm M_{\odot}$, $R_{\rm NSC}\sim 3\ \rm pc$ and $\gamma_{\rm NSC}\sim 1.5$.
Such eccentric sources can be identified with $e\gtrsim 10^{-4}-0.1$ at
$z\lesssim 10$ by 3rd generation GW detectors that are most sensitive at
$f_{\rm GW}\sim 0.1-10\ \rm Hz$.

    * –

Finally, the spin distributions can also be different. In this work we do not
consider the spins of remnants for simplicity. Here we briefly discuss the key
processes that determine the spin distributions for the two channels. It is
found in small-scale hydrodynamic simulations (e.g. Stacy et al. 2010, 2012,
2016; Greif et al. 2012; Stacy & Bromm 2013; Susa et al. 2014; Machida &
Nakamura 2015; Hirano & Bromm 2017; Susa 2019; Sugimura et al. 2020; Chiaki &
Yoshida 2020) that Pop III stars form in small groups from disk fragmentation.
In this case, the initial spins of individual stars and binary orbit vectors
will be mostly aligned imperfectly as the angular momenta are inherited from
the same disk and meanwhile affected by local flows during fragmentation and
mergers as well as N-body dynamics. In the NSC-BH channel, the binary orbital
plane can be tilted by dynamical interactions (see e.g. Trani et al. 2021),
while for close binaries in the BSE channel, the stellar spins are regulated
by tidal effects and mass transfer, whose effects are rather complex (see e.g.
Kinugawa et al. 2020; Tanikawa et al. 2021b; Stegmann & Antonini 2021). Natal
kicks can also alter the binary orbit vector with respect to the spins
(Tanikawa et al., 2021b).

We defer a detailed analysis of the detectable eccentricities and spins of Pop
III binary remnant mergers for the NSC-DH channel to future work.

  * •

Compared with LIGO events mostly involving BHs with $m_{\rm rem}\lesssim 50\
\rm M_{\odot}$, our Pop III mergers are generally more massive with $m_{\rm
rem}\sim 40-85\ \rm M_{\odot}$. Considering the local MRD inferred from LIGO
events, $\dot{n}_{\rm GW}(z\sim 0.33)=23.9_{-8.6}^{+14.3}\ \rm yr^{-1}\
Gpc^{-3}$ (Abbott et al., 2020c), Pop III remnants mergers from our NSC-DH
channel can count for up to a few percent of the 57 events detected so far
(Nitz et al., 2021), particularly the most massive ones. In general we expect
the low-redshift ($z\lesssim 2$) regime to be dominated by mergers from Pop
I//II progenitors, while detection of Pop III mergers with be most promising
at high redshifts ($z\gtrsim 5$). Actually, most of our Pop III binary remnant
mergers in NSCs can be detected by the 3rd generation of future GW
instruments, capable of reaching $z\gtrsim 10$, such as the Einstein Telescope
(Hild et al., 2009) and Decihertz Observatory (Arca Sedda et al., 2020), for
which we predict a full-sky detection rate of $\gtrsim 170-2720\ \rm yr^{-1}$.

  * •

Our models can produce BH-BH mergers similar to the special event GW190521
(Abbott et al., 2020a, b; Mehta et al., 2021; O’Brien et al., 2021) that
involves two BHs in the standard PPISN mass gap $\sim 55-130\ \rm M_{\odot}$
for Pop I/II stars (see e.g. Heger et al. 2003; Belczynski et al. 2016;
Woosley 2017; Marchant et al. 2019). The predicted MRD for such mergers agrees
with that inferred from GW190521 within the 90% confidence interval. This
further supports the proposal that mergers of BHs in the standard PPISN mass
gap can originate from Pop III stars (Farrell et al., 2021; Tanikawa et al.,
2021a; Kinugawa et al., 2021a; Liu & Bromm, 2020c). Moreover, the NSC-DH
channel does not suffer from the possible reduction of BH mass due to close
binary interactions in the BSE channel (Farrell et al., 2021). Actually, at
low redshifts when mergers from different formation channels are mixed
together, the mass gaps for (P)PISNe can be useful windows to distinguish the
mergers from different stellar populations since the location of mass gap
depends on metallicity191919The exact dependence of mass gap location on
metallicity is still uncertain in current stellar evolution models due to
uncertainties in e.g. mass loss and convective overshooting. (Spera & Mapelli,
2017; Heger & Woosley, 2010; Abbott et al., 2020a, b; Vink et al., 2021). This
is particularly relevant for mergers involving high-mass ($\gtrsim 50\ \rm
M_{\odot}$), low-spin ($\lesssim 0.2$) BHs that are not likely to be
hierarchical mergers (Gerosa et al., 2021). In our case, Pop III stars can be
the unique progenitors for mergers with low-spin BHs in the mass range of
$\sim 55-85\ \rm M_{\odot}$. For mergers involving at least one BH in this
mass range, our models predict a full-sky detection rate of $\sim 110-1640\
\rm yr^{-1}$ for $z_{\max}=10$.

In general, our results show that the NSC-DH channel is efficient at producing
Pop III binary remnant mergers, which not only complements the well-known BSE
channel, but also highlights the importance of environment for the evolution
of compact object binaries. The GW events from this channel will be detected
at promising rates in the next decades. Once distinguished from other
populations by their unique features with a large enough sample, they will
become a useful probe to the first stars and also high-$z$ NSCs formed in the
first galaxies. In concluding, we would like to discuss the caveats of our
model and suggest avenues for future work:

  * •

Our merger trees are targeted towards a MW-like halo virialized at $z\simeq
2.5$, such that they are only cosmologically representative at high redshifts
($z\gtrsim 5$, before the turn-over of the target halo). As a result, the
population of Pop III binary remnant mergers in NSCs at $z\lesssim 5$ may be
incomplete in our models. And we have used simple assumptions and
extrapolation to estimate the MRD in this low-redshift regime. In future work,
we will generate a broader range of target haloes at $z=0$ and combine their
results weighted by the halo mass function to construct a cosmologically
representative population of Pop III binary remnant mergers in NSCs down to
$z=0$ (see e.g. Magg et al. 2016).

  * •

As a common shortcoming of merger trees generated by the extended Press-
Schechter formalism (Parkinson et al., 2008; Mo et al., 2010),
positional/clustering information of haloes is lacking in our stellar feedback
model, such that radiative feedback can only be described by a uniform
background, and metal enrichment is implemented in a stochastic manner.
Actually, the stochastic metal enrichment model in HT15 does not achieve
numerical convergence with redshift/time resolution such that we have further
introduced a scaling factor for the metal-enriched volume to calibrate our
predicted SFRDs with those in simulations and observations. Note that stellar
feedback regulates not only the overall Pop III SF history, but also the host
haloes and basic units of Pop III SF that are important for binary statistics
(Liu & Bromm, 2020b; Liu et al., 2021). Future studies should adopt more
advanced feedback models taking into account local effects with merger trees
constructed from cosmological simulations in which the positions of haloes are
available (see e.g. Magg et al. 2018; Hartwig et al. 2018).

  * •

Crucially, as DF and DH are complex processes in the context of galaxy
formation and evolution, our semi-analytical approach relies on idealized
assumptions for galaxy and NSC properties based on empirical constraints in
the local Universe. In this way, we are only able to predict a broad MRD range
for the Pop III NSC-DH channel by considering limits for idealized (optimistic
and pessimistic) scenarios. The detailed complex physics involved in the in-
fall and evolution of compact object binaries in NSCs, such as galaxy dynamics
and mergers, NSC formation, growth and structure, as well as central massive
BHs, can only be modelled self-consistently with cosmological hydrodynamic
simulations. Simulating structure formation with these physical processes
included and enough resolution for Pop III SF in a representative cosmological
volume (down to $z=0$) is a challenging but promising task. As an exploratory
step in this direction, we plan to bulid sub-grid models for NSCs and central
massive BHs based on existing analytical models (e.g. Devecchi & Volonteri
2009; Devecchi et al. 2010, 2012), and sub-grid models designed for
cosmological simulations with limited resolution (e.g. the E-MOSAICS project,
Pfeffer et al. 2018). These models will be implemented and tested in meso-
scale simulations (in a co-moving volume of a few hundred $\rm Mpc^{3}$) that
can resolve minihaloes (Jaacks et al., 2019), coupled with existing schemes
for Pop III and Pop II SF and stellar feedback (see e.g. Johnson et al. 2013;
Xu et al. 2016; Sarmento et al. 2018; Liu & Bromm 2020a, b).

Identifying observational probes of the elusive first generation of stars has
been a long-term challenge, given that any such signatures are typically
hidden in the dominant foreground from subsequent SF. Actually, Pop III stars,
formed mostly before reionization ($z\gtrsim 6$), only make up $\sim 10^{-5}$
of all stars (ever formed) in the Universe. However, as shown by our results
for the NSC-DH channel and previous studies for the BSE channel, a much higher
fraction ($\sim 10^{-3}-0.1$) of massive ($m_{\bullet}\gtrsim 20\ \rm
M_{\odot}$) compact object mergers can originate from Pop III stars, including
those observed by LIGO in the local Universe. This is a natural
consequence/extrapolation from the strong metallicity dependence of the
efficiency of producing massive compact object binaries found in Pop I/II
stars (e.g. Dominik et al. 2015; Giacobbo & Mapelli 2018; Mapelli et al. 2018,
2019). In the next decades when 3rd-generation GW detectors map the demography
of BHs and NSs with thousands of compact object mergers, ‘gravitational-wave
archaeology’, targeting the massive BH remnants from Pop III, will provide an
unusually efficient channel to test theories of the first stars.

## Acknowledgements

The authors would like to thank Tilman Hartwig for helpful comments. This work
was supported by National Science Foundation (NSF) grant AST-1413501.

## Data availability

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

## References

  * Abbott et al. (2020a) Abbott R., et al., 2020a, Phys. Rev. Lett., 125, 101102
  * Abbott et al. (2020b) Abbott R., et al., 2020b, ApJ, 900, L13
  * Abbott et al. (2020c) Abbott R., et al., 2020c, ApJ, 913, L7
  * Ade et al. (2014) Ade P. A., et al., 2014, A&A, 571, A16
  * Adhikari et al. (2020) Adhikari S., Fishbach M., Holz D. E., Wechsler R. H., Fang Z., 2020, ApJ, 905, 21
  * Aghanim et al. (2020) Aghanim N., et al., 2020, A&A, 641, A1
  * Alvarez et al. (2009) Alvarez M. A., Wise J. H., Abel T., 2009, ApJ, 701, L133
  * Antonini & Rasio (2016) Antonini F., Rasio F. A., 2016, ApJ, 831, 187
  * Arca-Sedda (2016) Arca-Sedda M., 2016, MNRAS, 455, 35
  * Arca-Sedda & Capuzzo-Dolcetta (2014) Arca-Sedda M., Capuzzo-Dolcetta R., 2014, MNRAS, 444, 3738
  * Arca-Sedda et al. (2015) Arca-Sedda M., Capuzzo-Dolcetta R., Antonini F., Seth A., 2015, ApJ, 806, 220
  * Arca Sedda et al. (2020) Arca Sedda M., et al., 2020, Classical and Quantum Gravity, 37, 215011
  * Barack et al. (2019) Barack L., et al., 2019, Classical and quantum gravity, 36, 143001
  * Behroozi et al. (2019) Behroozi P., Wechsler R. H., Hearin A. P., Conroy C., 2019, MNRAS, 488, 3143
  * Belczynski et al. (2016) Belczynski K., et al., 2016, A&A, 594, A97
  * Belczynski et al. (2017) Belczynski K., Ryu T., Perna R., Berti E., Tanaka T. L., Bulik T., 2017, MNRAS, 471, 4702
  * Benitez-Llambay & Frenk (2020) Benitez-Llambay A., Frenk C., 2020, MNRAS, 498, 4887
  * Binney & Tremaine (2011) Binney J., Tremaine S., 2011, Galactic dynamics. Princeton university press
  * Bouffanais et al. (2021) Bouffanais Y., Mapelli M., Santoliquido F., Giacobbo N., Di Carlo U. N., Rastello S., Artale M. C., Iorio G., 2021, arXiv e-prints, p. arXiv:2102.12495
  * Bromm (2013) Bromm V., 2013, Reports on Progress in Physics, 76, 112901
  * Bromm et al. (2002) Bromm V., Coppi P. S., Larson R. B., 2002, ApJ, 564, 23
  * Bromm et al. (2009) Bromm V., Yoshida N., Hernquist L., McKee C. F., 2009, Nature, 459, 49
  * Calvi et al. (2012) Calvi R., Poggianti B. M., Fasano G., Vulcani B., 2012, MNRAS, 419, L14
  * Campisi et al. (2011) Campisi M., Maio U., Salvaterra R., Ciardi B., 2011, MNRAS, 416, 2760
  * Chen et al. (2021) Chen N., Ni Y., Tremmel M., Di Matteo T., Bird S., DeGraf C., Feng Y., 2021, arXiv e-prints, p. arXiv:2104.00021
  * Chiaki & Yoshida (2020) Chiaki G., Yoshida N., 2020, arXiv e-prints, p. arXiv:2008.06107
  * Choksi et al. (2019) Choksi N., Volonteri M., Colpi M., Gnedin O. Y., Li H., 2019, ApJ, 873, 100
  * Chon et al. (2021) Chon S., Omukai K., Schneider R., 2021, arXiv e-prints, p. arXiv:2103.04997
  * De Luca et al. (2021) De Luca V., Franciolini G., Pani P., Riotto A., 2021, J. Cosmology Astropart. Phys., 2021, 003
  * Dehnen (1993) Dehnen W., 1993, MNRAS, 265, 250
  * Devecchi & Volonteri (2009) Devecchi B., Volonteri M., 2009, ApJ, 694, 302
  * Devecchi et al. (2010) Devecchi B., Volonteri M., Colpi M., Haardt F., 2010, MNRAS, 409, 1057
  * Devecchi et al. (2012) Devecchi B., Volonteri M., Rossi E., Colpi M., Portegies Zwart S., 2012, MNRAS, 421, 1465
  * Dominik et al. (2015) Dominik M., et al., 2015, ApJ, 806, 263
  * Dvorkin et al. (2016) Dvorkin I., Vangioni E., Silk J., Uzan J.-P., Olive K. A., 2016, MNRAS, 461, 3877
  * Farr et al. (2019) Farr W. M., Fishbach M., Ye J., Holz D. E., 2019, ApJ, 883, L42
  * Farrell et al. (2021) Farrell E., Groh J. H., Hirschi R., Murphy L., Kaiser E., Ekström S., Georgy C., Meynet G., 2021, MNRAS, 502, L40
  * Finkelstein (2016) Finkelstein S. L., 2016, Publ. Astron. Soc. Australia, 33, e037
  * Fishbach et al. (2018) Fishbach M., Holz D. E., Farr W. M., 2018, ApJ, 863, L41
  * Fishbach et al. (2021) Fishbach M., et al., 2021, ApJ, 912, 98
  * Fragione & Silk (2020) Fragione G., Silk J., 2020, MNRAS, 498, 4591
  * Gardner et al. (2006) Gardner J. P., et al., 2006, Space Sci. Rev., 123, 485
  * Gerosa et al. (2021) Gerosa D., Giacobbo N., Vecchio A., 2021, arXiv e-prints, p. arXiv:2104.11247
  * Giacobbo & Mapelli (2018) Giacobbo N., Mapelli M., 2018, MNRAS, 480, 2011
  * Greif et al. (2012) Greif T. H., Bromm V., Clark P. C., Glover S. C., Smith R. J., Klessen R. S., Yoshida N., Springel V., 2012, MNRAS, 424, 399
  * Haemmerlé et al. (2020) Haemmerlé L., Mayer L., Klessen R. S., Hosokawa T., Madau P., Bromm V., 2020, Space Sci. Rev., 216, 1
  * Hartwig et al. (2015) Hartwig T., Bromm V., Klessen R. S., Glover S. C., 2015, MNRAS, 447, 3892
  * Hartwig et al. (2016) Hartwig T., Volonteri M., Bromm V., Klessen R. S., Barausse E., Magg M., Stacy A., 2016, MNRAS, 460, L74
  * Hartwig et al. (2018) Hartwig T., et al., 2018, MNRAS, 478, 1795
  * Heger & Woosley (2010) Heger A., Woosley S. E., 2010, ApJ, 724, 341
  * Heger et al. (2003) Heger A., Fryer C., Woosley S., Langer N., Hartmann D. H., 2003, ApJ, 591, 288
  * Hijikawa et al. (2021) Hijikawa K., Tanikawa A., Kinugawa T., Yoshida T., Umeda H., 2021, MNRAS,
  * Hild et al. (2009) Hild S., Chelkowski S., Freise A., Franc J., Morgado N., Flaminio R., DeSalvo R., 2009, Classical and Quantum Gravity, 27, 015003
  * Hirano & Bromm (2017) Hirano S., Bromm V., 2017, MNRAS, 470, 898
  * Hirano et al. (2014) Hirano S., Hosokawa T., Yoshida N., Umeda H., Omukai K., Chiaki G., Yorke H. W., 2014, ApJ, 781, 60
  * Hutter et al. (2021) Hutter A., Dayal P., Yepes G., Gottlöber S., Legrand L., Ucci G., 2021, MNRAS, 503, 3698
  * Inayoshi et al. (2016) Inayoshi K., Kashiyama K., Visbal E., Haiman Z., 2016, MNRAS, 461, 2722
  * Jaacks et al. (2019) Jaacks J., Finkelstein S. L., Bromm V., 2019, MNRAS, 488, 2202
  * Johnson & Bromm (2007) Johnson J. L., Bromm V., 2007, MNRAS, 374, 1557
  * Johnson et al. (2013) Johnson J. L., Dalla V. C., Khochfar S., 2013, MNRAS, 428, 1857
  * Katz et al. (2021) Katz H., et al., 2021, arXiv e-prints, p. arXiv:2101.11624
  * Kinugawa et al. (2014) Kinugawa T., Inayoshi K., Hotokezaka K., Nakauchi D., Nakamura T., 2014, MNRAS, 442, 2963
  * Kinugawa et al. (2015) Kinugawa T., Miyamoto A., Kanda N., Nakamura T., 2015, MNRAS, 456, 1093
  * Kinugawa et al. (2020) Kinugawa T., Nakamura T., Nakano H., 2020, MNRAS, 498, 3946
  * Kinugawa et al. (2021a) Kinugawa T., Nakamura T., Nakano H., 2021a, MNRAS, 501, L49
  * Kinugawa et al. (2021b) Kinugawa T., Nakamura T., Nakano H., 2021b, MNRAS, 504, L28
  * Lewis et al. (2000) Lewis A., Challinor A., Lasenby A., 2000, ApJ, 538, 473
  * Liu & Bromm (2020a) Liu B., Bromm V., 2020a, MNRAS, 495, 2475
  * Liu & Bromm (2020b) Liu B., Bromm V., 2020b, MNRAS, 497, 2839
  * Liu & Bromm (2020c) Liu B., Bromm V., 2020c, ApJ, 903, L40
  * Liu et al. (2021) Liu B., Meynet G., Bromm V., 2021, MNRAS, 501, 643
  * Machida & Nakamura (2015) Machida M. N., Nakamura T., 2015, MNRAS, 448, 1405
  * Madau & Dickinson (2014) Madau P., Dickinson M., 2014, ARA&A, 52, 415
  * Magg et al. (2016) Magg M., Hartwig T., Glover S. C., Klessen R. S., Whalen D. J., 2016, MNRAS, 462, 3591
  * Magg et al. (2018) Magg M., Hartwig T., Agarwal B., Frebel A., Glover S. C., Griffen B. F., Klessen R. S., 2018, MNRAS, 473, 5308
  * Mapelli et al. (2018) Mapelli M., Giacobbo N., Toffano M., Ripamonti E., Bressan A., Spera M., Branchesi M., 2018, Monthly Notices of the Royal Astronomical Society, 481, 5324
  * Mapelli et al. (2019) Mapelli M., Giacobbo N., Santoliquido F., Artale M. C., 2019, MNRAS, 487, 2
  * Mapelli et al. (2021) Mapelli M., et al., 2021, MNRAS,
  * Marchant et al. (2019) Marchant P., Renzo M., Farmer R., Pappas K. M., Taam R. E., De Mink S. E., Kalogera V., 2019, ApJ, 882, 36
  * Mastrogiovanni et al. (2021) Mastrogiovanni S., et al., 2021, arXiv e-prints, p. arXiv:2103.14663
  * Mehta et al. (2021) Mehta A. K., Buonanno A., Gair J., Miller M. C., Farag E., deBoer R. J., Wiescher M., Timmes F. X., 2021, arXiv e-prints, p. arXiv:2105.06366
  * Mo et al. (2010) Mo H., Van den Bosch F., White S., 2010, Galaxy formation and evolution. Cambridge University Press
  * Neumayer et al. (2020) Neumayer N., Seth A., Böker T., 2020, A&ARv, 28, 4
  * Nitz et al. (2021) Nitz A. H., Capano C. D., Kumar S., Wang Y.-F., Kastha S., Schäfer M., Dhurkunde R., Cabero M., 2021, arXiv e-prints, p. arXiv:2105.09151
  * O’Brien et al. (2021) O’Brien B., Szczepanczyk M., Gayathri V., Bartos I., Vedovato G., Prodi G., Mitselmakher G., Klimenko S., 2021, arXiv e-prints, p. arXiv:2106.00605
  * Paardekooper et al. (2015) Paardekooper J.-P., Khochfar S., Dalla Vecchia C., 2015, MNRAS, 451, 2544
  * Parkinson et al. (2008) Parkinson H., Cole S., Helly J., 2008, MNRAS, 383, 557
  * Pawlik et al. (2015) Pawlik A. H., Schaye J., Vecchia C. D., 2015, MNRAS, 451, 1586
  * Pawlik et al. (2017) Pawlik A. H., Rahmati A., Schaye J., Jeon M., Dalla Vecchia C., 2017, MNRAS, 466, 960
  * Pechetti et al. (2020) Pechetti R., Seth A., Neumayer N., Georgiev I., Kacharov N., den Brok M., 2020, ApJ, 900, 32
  * Perna et al. (2019) Perna R., Wang Y.-H., Farr W. M., Leigh N., Cantiello M., 2019, ApJ, 878, L1
  * Peters (1964) Peters P. C., 1964, Phys. Rev., 136, 1224
  * Pfeffer et al. (2018) Pfeffer J., Kruijssen J. D., Crain R. A., Bastian N., 2018, MNRAS, 475, 4309
  * Roulet et al. (2021) Roulet J., Chia H. S., Olsen S., Dai L., Venumadhav T., Zackay B., Zaldarriaga M., 2021, arXiv e-prints, p. arXiv:2105.10580
  * Safarzadeh (2020) Safarzadeh M., 2020, ApJ, 892, L8
  * Safarzadeh & Berger (2019) Safarzadeh M., Berger E., 2019, ApJ, 878, L12
  * Safranek-Shrader et al. (2016) Safranek-Shrader C., Montgomery M. H., Milosavljević M., Bromm V., 2016, MNRAS, 455, 3288
  * Sarmento et al. (2018) Sarmento R., Scannapieco E., Cohen S., 2018, ApJ, 854, 75
  * Schauer et al. (2019) Schauer A. T., Glover S. C., Klessen R. S., Ceverino D., 2019, MNRAS, 484, 3510
  * Schauer et al. (2020) Schauer A. T. P., Drory N., Bromm V., 2020, ApJ, 904, 145
  * Schödel et al. (2018) Schödel R., Gallego-Cano E., Dong H., Nogueras-Lara F., Gallego-Calvente A., Amaro-Seoane P., Baumgardt H., 2018, A&A, 609, A27
  * Sedda (2020) Sedda M. A., 2020, ApJ, 891, 47
  * Sesana & Khan (2015) Sesana A., Khan F. M., 2015, MNRAS, 454, L66
  * Sesana et al. (2006) Sesana A., Haardt F., Madau P., 2006, ApJ, 651, 392
  * Smith et al. (2018) Smith B. D., Regan J. A., Downes T. P., Norman M. L., O’Shea B. W., Wise J. H., 2018, MNRAS, 480, 3762
  * So et al. (2014) So G. C., Norman M. L., Reynolds D. R., Wise J. H., 2014, ApJ, 789, 149
  * Spera & Mapelli (2017) Spera M., Mapelli M., 2017, MNRAS, 470, 4739
  * Stacy & Bromm (2013) Stacy A., Bromm V., 2013, MNRAS, 433, 1094
  * Stacy et al. (2010) Stacy A., Greif T. H., Bromm V., 2010, MNRAS, 403, 45
  * Stacy et al. (2012) Stacy A., Greif T. H., Bromm V., 2012, MNRAS, 422, 290
  * Stacy et al. (2016) Stacy A., Bromm V., Lee A. T., 2016, MNRAS, 462, 1307
  * Stegmann & Antonini (2021) Stegmann J., Antonini F., 2021, Phys. Rev. D, 103, 063007
  * Sugimura et al. (2020) Sugimura K., Matsumoto T., Hosokawa T., Hirano S., Omukai K., 2020, ApJ, 892, L14
  * Susa (2019) Susa H., 2019, ApJ, 877, 99
  * Susa et al. (2014) Susa H., Hasegawa K., Tominaga N., 2014, ApJ, 792, 32
  * Tang et al. (2020) Tang P. N., Eldridge J., Stanway E. R., Bray J. C., 2020, MNRAS, 493, L6
  * Tanikawa et al. (2020) Tanikawa A., Yoshida T., Kinugawa T., Takahashi K., Umeda H., 2020, MNRAS, 495, 4170
  * Tanikawa et al. (2021a) Tanikawa A., Kinugawa T., Yoshida T., Hijikawa K., Umeda H., 2021a, MNRAS,
  * Tanikawa et al. (2021b) Tanikawa A., Susa H., Yoshida T., Trani A. A., Kinugawa T., 2021b, ApJ, 910, 30
  * Tarumi et al. (2020) Tarumi Y., Hartwig T., Magg M., 2020, ApJ, 897, 58
  * The LIGO Scientific Collaboration et al. (2021) The LIGO Scientific Collaboration the Virgo Collaboration the KAGRA Collaboration 2021, arXiv e-prints, p. arXiv:2105.15120
  * Trani et al. (2021) Trani A. A., Tanikawa A., Fujii M. S., Leigh N. W. C., Kumamoto J., 2021, MNRAS, 504, 910
  * Vink et al. (2021) Vink J. S., Higgins E. R., Sander A. A. C., Sabhahit G. N., 2021, MNRAS, 504, 146
  * Visbal et al. (2020) Visbal E., Bryan G. L., Haiman Z., 2020, ApJ, 897, 95
  * Vitale et al. (2019) Vitale S., Farr W. M., Ng K. K., Rodriguez C. L., 2019, ApJ, 886, L1
  * Wang et al. (2021) Wang Y.-Z., Tang S.-P., Liang Y.-F., Han M.-Z., Li X., Jin Z.-P., Fan Y.-Z., Wei D.-M., 2021, ApJ, 913, 42
  * Woosley (2017) Woosley S., 2017, ApJ, 836, 244
  * Xu et al. (2016) Xu H., Norman M. L., O’Shea B. W., Wise J. H., 2016, ApJ, 823, 140