# Gravitational waves from mergers of Population III binary black holes: roles played by two evolution channels

Boyuan Liu, Tilman Hartwig, Nina S Sartorio, Irina Dvorkin, Guglielmo Costa, Filippo Santoliquido, Anastasia Fialkov, Ralf S Klessen, Volker Bromm

<figure><img src="image_1.png" /></figure>

## To cite this version:

Boyuan Liu, Tilman Hartwig, Nina S Sartorio, Irina Dvorkin, Guglielmo Costa, et al. Gravitational waves from mergers of Population III binary black holes: roles played by two evolution channels. Monthly Notices of the Royal Astronomical Society, 2024, 534 (3), pp.1634-1667. 10.1093/mnras/stae2120. hal-04650435

HAL Id: hal-04650435

https://hal.science/hal-04650435v1

Submitted on 18 Jul 2025

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>
# Gravitational waves from mergers of Population III binary black holes: roles played by two evolution channels

Boyuan Liu<sup>D, 1,2,3★</sup> Tilman Hartwig<sup>D, 4,5,6</sup> Nina S. Sartorio<sup>D, 1,7</sup> Irina Dvorkin<sup>D, 8</sup> Guglielmo Costa<sup>D, 9</sup> Filippo Santoliquido<sup>D, 10,11</sup> Anastasia Fialkov,<sup>1,12</sup> Ralf S. Klessen<sup>D, 3,13</sup> and Volker Bromm<sup>D, 2,14</sup>

Accepted 2024 September 9. Received 2024 August 26; in original form 2024 June 25

## ABSTRACT

The gravitational wave (GW) signal from binary black hole (BBH) mergers is a promising probe of Population III (Pop III) stars. To fully unleash the power of the GW probe, one important step is to understand the relative importance and features of different BBH evolution channels. We model two channels, isolated binary stellar evolution (IBSE) and nuclear star cluster-dynamical hardening (NSC-DH), in one theoretical framework based on the semi-analytical code A-SLOTH, under various assumptions on Pop III initial mass function (IMF), initial binary statistics and high-z nuclear star clusters (NSCs). The NSC-DH channel contributes ~ 8–95 per cent of Pop III BBH mergers across cosmic history, with higher contributions achieved by initially wider binary stars, more top-heavy IMFs, and more abundant high-z NSCs. The dimensionless stochastic GW background (SGWB) produced by Pop III BBH mergers has peak values $\Omega_{\text{GW}}^{\text{peak}} \sim 10^{-11}–8 \times 10^{-11}$ around observer-frame frequencies $\nu \sim 10–100$ Hz. The Pop III contribution can be a non-negligible (~ 2–32 per cent) component in the total SGWB at $\nu \lesssim 10$ Hz. The estimated detection rates of Pop III BBH mergers by the Einstein Telescope are ~ 6–230 and ~ 30–1230 $\text{yr}^{-1}$ for the NSC-DH and IBSE channels, respectively. Pop III BBH mergers in NSCs are more massive than those from IBSE, so they dominate the Pop III SGWB below 20 Hz in most cases. Besides, the detection rate of Pop III BBH mergers involving at least one intermediate-mass BH above 100 $M_{\odot}$ by the Einstein Telescope is ~ 0.5–200 $\text{yr}^{-1}$ in NSCs but remains below 0.1 $\text{yr}^{-1}$ for IBSE.

**Key words:** gravitational waves – stars: Population III – galaxies: star clusters: general – dark ages, reionization, first stars – black hole mergers.

## 1 INTRODUCTION

One of the major goals of modern astrophysics is to understand the onset of star and galaxy formation at Cosmic Dawn, during the first billion years after the big bang (Loeb & Furlanetto 2013). In particular, the first generation of so-called Population III (Pop III) stars, formed in extremely metal-poor primordial gas via cooling by molecular hydrogen, are believed to have distinct features compared with present-day, Population I/II (Pop I/II), stars (reviewed by e.g.

Bromm et al. 2009; Bromm 2013; Haemmerlé et al. 2020; Klessen & Glover 2023).

The first luminous objects play important roles in early cosmic history through their metal enrichment, radiation fields, and cosmic ray production (e.g. Pan, Scannapieco & Scalo 2013; Salvador-Solé et al. 2017; Jaacks et al. 2018; Ohira & Murase 2019; Gessey-Jones et al. 2023; Sartorio et al. 2023; Yamaguchi, Furlanetto & Trapp 2023), significantly impacting the chemical and thermal evolution of the intergalactic medium (IGM; see e.g. Karlsson, Bromm & Bland-Hawthorn 2013; Barkana 2016; Dayal & Ferrara 2018 for reviews). This early feedback also establishes the conditions for the formation and evolution of subsequent populations of stars, galaxies, and (supermassive) black hole (BH) seeds (see e.g. Bromm &

© 2024 The Author(s).

Published by Oxford University Press on behalf of Royal Astronomical Society. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted reuse, distribution, and reproduction in any medium, provided the original work is properly cited.

<footnote n="1"><sup>1</sup>Institute of Astronomy, University of Cambridge, Madingley Road, Cambridge CB3 0HA, UK</footnote>

<footnote n="2"><sup>2</sup>Department of Astronomy, University of Texas, Austin, TX 78712, USA</footnote>

<footnote n="3"><sup>3</sup>Universität Heidelberg, Zentrum für Astronomie, Institut für Theoretische Astrophysik, D-69120 Heidelberg, Germany</footnote>

<footnote n="4"><sup>4</sup>Institute for Physics of Intelligence, School of Science, The University of Tokyo, Bunkyo, Tokyo 113-0033, Japan</footnote>

<footnote n="5"><sup>5</sup>Department of Physics, School of Science, The University of Tokyo, Bunkyo, Tokyo 113-0033, Japan</footnote>

<footnote n="6"><sup>6</sup>AI-Lab, German Environment Agency, Alte Messe 6, D-04103 Leipzig, Germany</footnote>

<footnote n="7"><sup>7</sup>Department of Physics and Astronomy, University of Ghent, Technologiepark 903, Ghent, B-9052 Zwijnaarde, Belgium</footnote>

<footnote n="8"><sup>8</sup>Institut d'Astrophysique de Paris, Sorbonne Université & CNRS, UMR 7095, 98 bis bd Arago, F-75014 Paris, France</footnote>

<footnote n="9"><sup>9</sup>Univ Lyon, Univ Lyon1, ENS de Lyon, CNRS, Centre de Recherche Astrophysique de Lyon UMR5574, F-69230 Saint-Genis-Laval, France</footnote>

<footnote n="10"><sup>10</sup>Gran Sasso Science Institute (GSSI), E-67100 L'Aquila, Italy</footnote>

<footnote n="11"><sup>11</sup>INFN, Laboratori Nazionali del Gran Sasso, E-67100 Assergi, Italy</footnote>

<footnote n="12"><sup>12</sup>Kavli Institute for Cosmology, Madingley Road, Cambridge CB3 0HA, UK</footnote>

<footnote n="13"><sup>13</sup>Universität Heidelberg, Interdiszipläres Zentrum für Wissenschaftliches Rechnen, D-69120 Heidelberg, Germany</footnote>

<footnote n="14"><sup>14</sup>Weinberg Institute for Theoretical Physics, University of Texas, Austin, TX 78712, USA</footnote>

<footnote>*E-mail: treibeis1995@gmail.com</footnote>
Yoshida 2011; Johnson, Dalla Vecchia & Khochfar 2013; Pawlik, Milosavljević & Bromm 2013; Jeon et al. 2015; Smith et al. 2015; Regan et al. 2017; Sakurai et al. 2017; Haemmerlé et al. 2020; Inayoshi, Visbal & Haiman 2020; Schauer et al. 2021; Chon et al. 2022; Sarmento & Scannapieco 2022; Chiaki et al. 2023; Sanati et al. 2023; Regan & Volonteri 2024), providing powerful diagnostics for early structure formation and even fundamental physics, such as the nature of dark matter (e.g. Hirano & Bromm 2018; Hirano, Sullivan & Bromm 2018a; Sullivan, Hirano & Bromm 2018; Liu et al. 2019a; Liu, Schauer & Bromm 2019b; Nebrin, Ghara & Mellema 2019; Ilie et al. 2021; Cappelluti, Hasinger & Natarajan 2022; Driskell et al. 2022; Hibbard et al. 2022; Kulkarni et al. 2022; Liu & Bromm 2022, 2023; Liu, Zhang & Bromm 2022; Hirano & Yoshida 2024; Zhang, Bromm & Liu 2024a; Zhang, Liu & Bromm 2024b; Zhang, Ilie & Freese 2024c).

Although the importance of Pop III stars is evident in theoretical predictions, their detailed properties are still unclear in the absence of direct observations, which are still challenging (especially for pure Pop III systems in minihaloes at $z \gtrsim 15$) even with the JWST (Gardner et al. 2006; Nakajima & Maiolino 2022; Riaz, Hartwig & Latif 2022; Bovill et al. 2024) and calls for extremely large ($\gtrsim 100$ m) extraterrestrial telescopes (Angel et al. 2008; Rhodes et al. 2020; Schauer, Drory & Bromm 2020). Detection of Pop III features<fnref n="1" /> is possible in rare cases of magnification via extreme gravitational lensing (Schauer et al. 2022; Vanzella et al. 2023; Welch et al. 2022; Zackrisson et al. 2024) or trace Pop III populations<fnref n="2" /> co-existing with Pop I/II stars and active galactic nuclei (AGN) at $z \lesssim 10$ (e.g. Grisdale et al. 2021; Maiolino et al. 2024; Nanayakkara et al. 2023; Wang et al. 2024). Another complementary approach to constrain Pop III stars, in particular their initial mass function (IMF), with future telescopes, e.g. Roman Space Telescope, is to consider electromagnetic transients, such as gamma-ray bursts and supernova (SN) explosions (e.g. Fryer et al. 2022; Hartwig et al. 2023; Lazar & Bromm 2022; Venditti et al. 2024b; Wiggins et al. 2024). Currently, observational constraints on Pop III stars are mainly derived from two indirect probes: stellar archaeology and the 21-cm signal<fnref n="3" /> from neutral hydrogen at Cosmic Dawn.

Stellar archaeology connects the statistics and chemical patterns of extremely metal-poor stars observed in the local Universe with the IMF and SN properties of Pop III stars (reviewed by, e.g. Frebel & Norris 2015). Such extremely metal-poor stars are expected to be the bonafide second-generation (Pop II) stars that preserve the chemical imprints of Pop III metal enrichment (such as carbon enhancement, see e.g. Yoon et al. 2016, 2018; Hansen et al. 2019; Dietz et al. 2021; Zepeda et al. 2023). Recent progress in this direction (e.g. Hartwig et al. 2015; Ji, Frebel & Bromm 2015; Salvadori, Skúladóttir & Tolstoy 2015; de Bennassuti et al. 2017; Jeon, Besla & Bromm 2017; Sarmento, Scannapieco & Pan 2017; Ishigaki et al. 2018; Magg et al. 2018, 2019; Sharma et al. 2018; Salvadori et al. 2019; Chiaki et al. 2020; Komiya et al. 2020; Tarumi, Hartwig & Magg 2020; Jeon et al. 2021; Liu et al. 2021b; Rossi, Salvadori & Skúladóttir 2021;

Lucey et al. 2022; Hartwig et al. 2023; Koutsouridou, Salvadori & Skúladóttir 2024; Rossi et al. 2024; Skúladóttir et al. 2024) finds that although there is evidence for diverse explosion mechanisms of Pop III SNe including scenarios of jet-induced hypernovae, 'faint' SNe with fallback and mixing, as well as pair-instability SNe (PISNe; e.g. Heger & Woosley 2002, 2010; Heger et al. 2003; Maeda & Nomoto 2003; Umeda & Nomoto 2003, 2005; Iwamoto et al. 2005; Kobayashi et al. 2006; Tominaga 2009), the chemical patterns of extremely metal-poor stars are mostly consistent with the enrichment by core-collapse SNe from Pop III stars with initial masses $m_* \sim 20-40$ M$_\odot$, and more massive stars are rather rare or collapse directly into BHs at the end of metal-free stellar evolution (e.g. Heger et al. 2003; Fryer et al. 2012; Tanikawa et al. 2020), which potentially contribute to metal enrichment via strong post-main-sequence (MS) winds (Liu et al. 2021b; Jeena et al. 2023; Nandal et al. 2024a, b; Nandal, Sibony & Tsiatsiou 2024c; Tsiatsiou et al. 2024). Similarly, the imprints of Pop III metal enrichment can also be inferred from observations of high-$z$ extremely metal-poor absorption systems in the spectra of distant quasars (e.g. Welsh, Cooke & Fumagalli 2019; Welsh et al. 2022, 2023; Yoshii et al. 2022; Saccardi et al. 2023; Salvadori et al. 2023; Sodini et al. 2024; Vanni et al. 2024; Zou et al. 2024) and metal-poor galaxies with strong metal emission lines (Cameron et al. 2023; D'Eugenio et al. 2023; Ji et al. 2024b; Schaerer et al. 2024; Senchyna et al. 2024; Topping et al. 2024). In general, stellar archaeology confirms the massive nature of Pop III stars (with a broad, top-heavy IMF) as well as their formation in small clusters predicted by semi-analytical models and (magneto)hydrodynamic simulations of primordial star formation (e.g. Greif et al. 2011, 2012; Stacy & Bromm 2013; Hirano et al. 2014, 2015, 2018b; Susa, Hasegawa & Tominaga 2014; Stacy, Bromm & Lee 2016; Hirano & Bromm 2017; Riaz et al. 2018, 2023; Susa 2019; McKee, Stacy & Li 2020; Sharda, Federrath & Krumholz 2020; Sugimura et al. 2020, 2023; Wollenberg et al. 2020; Chon, Omukai & Schneider 2021; Latif et al. 2021; Park, Ricotti & Sugimura 2021, 2023, 2024; Sharda et al. 2021; Sharda & Krumholz 2022; Latif, Whalen & Khochfar 2022; Prole et al. 2022a, b, 2023; Liu et al. 2024a; Sadanari et al. 2024; Sharda & Menon 2024).

The 21-cm signal from the high-$z$ IGM is shaped by the UV and X-ray photons as well as cosmic rays from the first stars, SNe, and galaxies (reviewed by e.g. Barkana 2016), from which the timing, efficiency, and mode of early star formation, as well as source properties (e.g. Pop III IMF), can be constrained (e.g. Fialkov et al. 2013, 2017; Fialkov & Barkana 2014; Ma et al. 2018; Madau 2018; Mirocha & Furlanetto 2019; Schauer, Liu & Bromm 2019; Chatterjee et al. 2020; Qin et al. 2020; Gessey-Jones et al. 2022, 2023; Kamran et al. 2022; Kaur et al. 2022; Kovlakas et al. 2022; Magg et al. 2022b; Muñoz et al. 2022; Bera, Samui & Datta 2023; Fialkov, Gessey-Jones & Dhandha 2023; Hassan et al. 2023; Ma et al. 2023; Mondal & Barkana 2023; Ventura et al. 2023; Bevins et al. 2024; Lewis et al. 2024; Pochinda et al. 2024). For instance, the average star formation efficiency of Pop III stars is constrained to be below 5.5 at 68 per cent confidence by the 21-cm data from the Hydrogen Epoch of Reionization Array and Shaped Antenna measurement of the background RAdio Spectrum (Pochinda et al. 2024).

Another promising probe of the first stars that has just become possible thanks to the LIGO-Virgo-KAGRA (LVK) network (Acernese et al. 2015; LIGO Scientific Collaboration 2015; Kagra Collaboration 2019) is provided by gravitational waves (GWs) from compact object mergers (see e.g. Mapelli 2021; Mandel & Farmer 2022; Spera, Trani & Mencagli 2022 for recent reviews). So far, the LVK network has reported about one hundred events of mergers between

<footnote n="1">¹See Windhorst et al. (2018), Nakajima & Maiolino (2022), Trussler et al. (2023), Katz et al. (2023), Larkin, Gerasimov & Burgasser (2023), Zackrisson et al. (2024), and Venditti et al. (2024a) for detailed discussions of the observability and characterization of Pop III stars with JWST.</footnote>

<footnote n="2">²For instance, recently Maiolino et al. (2024) discover a broad (~170 Å) He IIλ1640 emission feature from NIRSpec-IFU and NIRSpec-MSA observations of GN-z11 at $z = 10.6$ that can be explained by a Pop III starburst of ~$6 \times 10^5$ M$_\odot$ (for an interpretation, see e.g. Venditti et al. 2023, 2024a).</footnote>

<footnote n="3">³In general, intensity mapping of other lines such as He IIλ1640 and H α can also constrain early star formation (e.g. Parsons et al. 2022).</footnote>
stellar-mass compact objects (mostly BHs with ~5–100 M⊙) within
z ≲ 1 in the third gravitational wave transient catalogue (GWTC-3;
Abbott et al. 2023b). Moreover, the 3rd-generation GW detectors
such as the Einstein Telescope (ET; Punturo et al. 2010; Maggiore
et al. 2020) and the Cosmic Explorer Observatory (CEO; Reitze
et al. 2019; Evans et al. 2023) will observe thousands of compact
object mergers per year up to z ~ 30. The mass, spin, and redshift
distributions of these mergers, especially at high redshifts, can place
novel constraints on Pop III star formation and stellar evolution.
Indeed, compared with Pop I/II stars, Pop III stars are more massive,
compact, and experience little mass-loss in the absence of metal-line-
driven winds, so they are more likely to produce (massive) compact
object remnants. Besides, recent hydrodynamic simulations converge
on the picture that Pop III stars are typically born in small clusters
(Haemmerlé et al. 2020; Hartwig et al. 2023; Klessen & Glover 2023)
in which the formation of binary and multiple systems are common.
Therefore, Pop III stars are expected to be efficient progenitors of
compact object binaries, in particular binary black holes (BBHs).
Recent theoretical studies (e.g. Kinugawa, Nakamura & Nakano
2020, 2021b; Liu & Bromm 2020a, 2021; Hijikawa et al. 2021;
Tanikawa et al. 2022b) reveal the possibility that mergers of Pop III
BBH remnants are abundant and energetic GW sources, contributing
a significant fraction (~0.1–10 per cent) of BBH mergers across
cosmic history, especially for massive BHs at high-z, much higher
than the Pop III fraction in the total mass budget of stars ever
formed in the Universe (≈ 10⁻<fnref n="4" />). A large fraction of the currently
observed BBH merger events may originate from Pop III stars (e.g.
Kinugawa et al. 2020, 2021b; Iwaya, Kinugawa & Tagoshi 2023),
which can explain the high-mass regime (> 20 M⊙) of the chirp
mass distribution for the detected events (Abbott et al. 2021, 2023b).

For instance, the recently detected event GW190521 at z = 0.82<sup>+0.28</sup><sub>-0.34</sub> with unusual BH masses 85<sup>+21</sup><sub>-14</sub> and 66<sup>+17</sup><sub>-18</sub> M<sub>☉</sub> (inferred with the NRSur7dq4 waveform model by Abbott et al. 2020a, see also Abbott et al. 2020b; O'Brien et al. 2021; Mehta et al. 2022) is suspected to have a Pop III origin.<sup>4</sup> GW190521 is special because the mass range 50–130 M<sub>☉</sub> is mostly forbidden for BHs born by Pop I/II stars according to the standard PISN models unless under peculiar conditions (e.g. Heger et al. 2003; Belczynski et al. 2016; Yoshida et al. 2016; Spera & Mapelli 2017; Woosley 2017, 2019; Farmer et al. 2019, 2020; Leung, Nomoto & Blinnikov 2019; Marchant et al. 2019; Mapelli et al. 2020; Marchant & Moriya 2020; Renzo et al. 2020a; Vink et al. 2021; Costa et al. 2021; Winch et al. 2024), which is called the (standard) PISN mass gap. The lower edge of the mass gap can be shifted up to ~100 M<sub>☉</sub> for Pop III stars with their unique evolution tracks, so Pop III stars can produce BBH mergers like GW190521 at a similar rate as observed under proper conditions (Liu & Bromm 2020c, 2021; Farrell et al. 2021; Kinugawa et al. 2021a; Tanikawa et al. 2021a, 2022b; Volpato et al. 2023; Tanikawa 2024).

However, there are significant discrepancies in current theoretical
predictions for Pop III BBH mergers, making it difficult to derive

accurate constraints on Pop III stars from current GW data. The
reason is that different studies consider different BBH formation and
evolution channels and adopt different assumptions on a variety of
physical processes over a broad range of scales. The most intensively
studied channel to date is isolated binary stellar evolution (IBSE)
of close binary stars (with initial separations *a* ≲ 10 au, see e.g.
Kinugawa et al. 2014, 2016, 2017, 2020; Hartwig et al. 2016;
Belczynski et al. 2017; Inayoshi et al. 2017; Hijikawa et al. 2021;
Tanikawa et al. 2021b, 2022b; Costa et al. 2023; Santoliquido et al.
2023), which relies on binary interactions (e.g. mass transfer, tidal
effects, and common envelope evolution) and SN natal kicks to
shrink the binary orbits and facilitate merging. The predicted merger
efficiencies and properties are highly sensitive to the uncertain initial
binary properties of Pop III stars (see e.g. Stacy & Bromm 2013;
Liu, Meynet & Bromm 2021a) and parameters of (binary) stellar
evolution (for broad discussions including Pop I/II stars, see e.g.
de Mink & Belczynski 2015; Kinugawa et al. 2020; Bavera et al.
2021; Marchant et al. 2021; Olejak, Belczynski & Ivanova 2021;
Santoliquido et al. 2021, 2023; Belczynski et al. 2022; Stevenson &
Clarke 2022; Tanikawa et al. 2022b; Iorio et al. 2023; Willcox
et al. 2023; Dorozsmai & Toonen 2024). Besides, it is implied
by recent hydrodynamic simulations of Pop III star formation (e.g.
Sugimura et al. 2020, 2023; Park, Ricotti & Sugimura 2023, 2024)
that most binaries of Pop III stars have wide orbits (*a* ≳ 100 au, due
to outward migrations of Pop III protostars and their circumstellar
discs by accretion of gas with high angular momentum), potentially
hampering the efficiency of the IBSE channel (Liu et al. 2021a; Costa
et al. 2023; Santoliquido et al. 2023).

In light of this, Liu & Bromm (2020c, 2021) propose an alternative
channel based on dynamical hardening (DH) in nuclear star clusters
(NSCs). In this scenario, Pop III BBHs fall into NSCs (made of
Pop I/II stars) by dynamical friction, in which the orbits of hard
binaries are shrunk by binary-single encounters, so BBHs from
initially wide binary stars can also merge within a Hubble time.
The NSC-DH channel is expected to be less sensitive to initial
binary properties and binary stellar evolution processes and can
be as efficient as the IBSE channel under favourable conditions.
Nevertheless, Pop III BBH mergers from the NSC-DH channel are
highly sensitive to the unknown properties of high-z dwarf galaxies
and their NSCs. Another dynamical channel (in situ) is given by N-
body dynamics within massive (~10⁴–10⁷ M☉) clusters of Pop III
stars and remnants themselves (Wang, Tanikawa & Fujii 2022;
Mestichelli et al. 2024; Liu et al. 2024b), which can also produce
non-negligible merger rates. In particular, it is shown in Mestichelli
et al. (2024) that massive BBH mergers involving intermediate-mass
BHs (IMBHs, ≳ 100 M☉) are common in massive Pop III clusters,
while such mergers can be extremely rare from isolated evolution
(Costa et al. 2023; Santoliquido et al. 2023). However, massive
Pop III star clusters (SCs) are only expected to form under special
conditions (e.g. strong dynamical heating, streaming motion, and
Lyman–Werner radiation, Li, Inayoshi & Qiu 2021; Lupi, Haiman &
Volonteri 2021), and their abundance and detailed properties are still
uncertain.<fnref n="5" />

To fully unleash the power of the GW probe, deeper and more
systematic investigations of the relevant physics of BBH formation

<footnote n="4">⁴The origin and properties of GW190521 are still in debate (see e.g. Fishbach & Holz 2020; Gayathri et al. 2020; Romero-Shaw et al. 2020). Beyond the Pop III scenario, other scenarios, such as IBSE of Pop I/II stars with special stellar evolution models, hierarchical BH mergers and stellar collisions in star clusters and discs of AGNs, can also explain this event (e.g. Belczynski 2020; Di Carlo et al. 2020; Fragione, Loeb & Rasio 2020; Kremer et al. 2020; Renzo et al. 2020b; Arca-Sedda et al. 2021; Costa et al. 2021, 2022; Gayathri et al. 2021; Gerosa & Fishbach 2021; Kimball et al. 2021; Liu & Lai 2021; Mapelli et al. 2021; Tagawa et al. 2021; Anagnostou, Trenti & Melatos 2022; Antonini, Romero-Shaw & Callister 2024).</footnote>

<footnote n="5">⁵In our simulations, most Pop III stars form in small clusters with total stellar masses of ~10³ M☉, and the contribution of Pop III stars in massive (~10⁴–10⁵ M☉) clusters to the total number of Pop III stars ever formed is very small (≪ 3 per cent). Therefore, we do not consider the pure Pop III (in situ) dynamical channel in this work.</footnote>
and evolution are required to establish robust connections between properties of the first stars (e.g. IMF, binary statistics at birth, and star formation history) and GW observations. One important step towards this ultimate goal is to understand the relative importance of different BBH formation/evolution channels as well as the corresponding features of GW sources. To do so, we implement the IBSE and NSC-DH channels in one theoretical framework based on the public semi-analytical model A-SLOTH<fnref n="6" /> (Ancient Stars and Local Observables by Tracing haloes; Hartwig et al. 2022; Magg et al. 2022a). Given the input of binary population synthesis (BPS) data, our framework for the first time follows the (external) dynamics (in host NSCs/galaxies/haloes and large-scale structure formation) as well as the (internal) orbital evolution of Pop III BBHs (in NSCs) on-the-fly, together with the underlying star/galaxy/structure formation process. This enables us to model different BBH evolution channels self-consistently and characterize the resulting Pop III BBH mergers and their host systems. In this paper, we explore the key properties of Pop III BBH mergers from the NSC-DH and IBSE channels in 18 simulations combining different assumptions on the yet uncertain (initial) properties of Pop III binary stars and high-z NSCs, based on the BPS results produced by the SEVN code (Costa et al. 2023) and the merger trees from the cosmological simulation in (Ishiyama et al. 2016). We plan to investigate the properties of host haloes/galaxies/NSCs of Pop III BBH mergers, the rate and spatial distribution of mergers in the host galaxy, and the dependence of merger properties on star formation, (binary) stellar evolution, and feedback parameters in future work (see e.g. Pacucci, Loeb & Salvadori 2017; Artale et al. 2019; Santoliquido et al. 2021, 2022, 2023; Iorio et al. 2023; Rauf et al. 2023; Srinivasan et al. 2023).

A-SLOTH was calibrated based on six independent observational constraints. It has been refined and applied to a variety of topics, including GWs from mergers of Pop III remnants (but treating the IBSE and NSC-DH channels separately, Hartwig et al. 2016; Liu & Bromm 2021), stellar archaeology constraints on the Pop III IMF (Hartwig et al. 2015; Magg et al. 2018; Tarumi et al. 2020) and local baryonic streaming velocity (Uysal & Hartwig 2023), hypernova signatures of Pop III stars in nearby dwarf satellite galaxies (Lee, Jeon & Bromm 2024a), impact of the transition from Pop III to Pop II star formation on the global 21-cm signal (Magg et al. 2022b), rates of Pop III SNe (Magg et al. 2016), formation of Milky Way (MW) satellites (Chen et al. 2022b, c), Pop III signatures and evolution of dust in high-z galaxies (Riaz et al. 2022; Dzięciol, Hartwig & Yoshida 2024), and the imprint of Pop III stars on the faint end of the MW white dwarf luminosity function. This broad range of applications implies that our model can be easily used to explore the correlations between the GW signals from Pop III (and Pop I/II) BBH mergers and other observational properties of the first stars and galaxies (see e.g. Dvorkin et al. 2016; Inayoshi et al. 2016; Santoliquido et al. 2022; Tanikawa et al. 2022a; Veronesi, van Velzen & Rossi 2024), although this is not the focus of this paper. Considering such broader applications, we make our code public.<fnref n="7" />

This paper is organized as follows. In Sections 2–4, we describe our implementation of the relevant physics in A-SLOTH, covering the formation of Pop III BBHs (Section 2.1), their subsequent evolution in different environments (Section 2.2), the formation and evolution of NSCs (Section 3), and an extension of A-SLOTH to follow the evolution of Pop III BBHs and their host galaxies/haloes in late epoch

<table><caption>**Table 1.** List of key abbreviations.</caption><thead><tr><td>Abbreviation</td><td>Full name</td></tr></thead><tbody><tr><td>GW</td><td>Gravitational wave</td></tr><tr><td>BH</td><td>Black hole</td></tr><tr><td>BBH</td><td>Binary black hole</td></tr><tr><td>SN</td><td>Supernova</td></tr><tr><td>PISN</td><td>Pair-instability supernova</td></tr><tr><td>MS</td><td>Main sequence</td></tr><tr><td>ZAMS</td><td>Zero-age main sequence</td></tr><tr><td>IMF</td><td>Initial mass function</td></tr><tr><td>BPS</td><td>Binary population synthesis</td></tr><tr><td>IBS</td><td>Initial binary statistics</td></tr><tr><td>IBSE</td><td>Isolated binary stellar evolution</td></tr><tr><td>NSC</td><td>Nuclear star cluster</td></tr><tr><td>SC</td><td>Star cluster</td></tr><tr><td>GC</td><td>Globular cluster</td></tr><tr><td>DF</td><td>Dynamical friction</td></tr><tr><td>DH</td><td>Dynamical hardening</td></tr><tr><td>CE</td><td>Common envelope</td></tr><tr><td>SMT</td><td>Stable mass transfer</td></tr><tr><td>SFRD</td><td>Star formation rate density</td></tr><tr><td>MRD (n̄)</td><td>Merger rate density (co-moving)</td></tr><tr><td>SGWB (ΩGW)</td><td>Stochastic gravitational wave background</td></tr><tr><td>SNR</td><td>Signal-to-noise ratio</td></tr><tr><td>GWTC</td><td>Gravitational wave transient catalogue</td></tr><tr><td>ET</td><td>The Einstein Telescope</td></tr><tr><td>CEO</td><td>The Cosmic Explorer Observatory</td></tr><tr><td>LVK</td><td>LIGO-Virgo-KAGRA</td></tr><tr><td>JWST</td><td>James Webb Space Telescope</td></tr><tr><td>AGN</td><td>Active galactic nuclei/nucleus</td></tr></tbody></table>

(z ≲ 4.5) that is not covered by merger trees from cosmological simulations (Section 4). In Section 5, we present the predicted properties of Pop III BBH mergers (e.g. merger rate, contribution to the stochastic GW background, mass distribution, and simple demography of host systems) from the IBSE and NSC-DH channels, starting with the fiducial model (Section 5.1), and then exploring the dependence on Pop III IMF, initial binary statistics (IBS), and abundances of (high-z) NSCs (Section 5.2). Finally, in Section 6, we summarize our results and discuss their implications for the GW signatures of the first stars. A list of the key abbreviations used in this paper is given in Table 1.

## 2 POP III BINARY STARS AND BBHS

We implement the formation, dynamics, and internal orbital evolution of Pop III BBHs in galaxy fields and NSCs, as well as the formation and evolution of NSCs in the public semi-analytical model A-SLOTH (Hartwig et al. 2022) for cosmic star/structure formation along halo merger trees. The reader is referred to Hartwig et al. (2022) and Magg et al. (2022a) for a detailed introduction of A-SLOTH. Our implementation is based on the earlier work by Liu & Bromm (2021) with updated treatments for the sampling/formation of Pop III BBHs, formation and evolution of galaxies and NSCs. Below we explain our BBH, galaxy, and NSC models, focusing on the new features with respect to Liu & Bromm (2021). We start with the BBH model in this section covering their formation (Section 2.1) and evolution (Section 2.2) in the context of external (galactic) dynamics (Section 2.2.1) and internal orbital dynamics (Section 2.2.2). The flow chart of our BBH routine is shown in Fig. 1. The model of NSCs will be described in the next section.

<footnote n="6">⁶https://gitlab.com/thartwig/asloth</footnote>

<footnote n="7">⁷https://gitlab.com/Treibeis/a-sloth-cob</footnote>
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Overview of the Pop III BBH routine (Section 2) within A-SLOTH. Red rounded boxes: interfaces with the merger tree. Blue boxes: decision-making steps. Orange boxes: other (sub-)routines that contain or are called by the BBH routine. Here, the star formation routine of A-SLOTH predicts the expected mass $\delta M_{\star,\text{III}}$ of Pop III stars formed at each time-step and updates the galaxy mass $M_{\star}$ (total mass of Pop I/II stars) that is crucial for DF of BBHs in galaxy fields. Green boxes: key processes for BBH formation (Section 2.1) and evolution (Section 2.2) via the (1) IBSE and (2) NSC-DH channels.ª The external (galactic) orbits of the BBHs from the smaller halo during a halo merger event are randomized according to equation (1).ᵇ Newly born BBHs are sampled from the input SEVN catalogue (see Section 2.1). Their galactic orbits are initialized with equation (1).</figcaption></figure>

## 2.1 Formation of Pop III BBHs

We couple the star formation routine in A-SLOTH with input catalogues of Pop III BBHs from BPS simulations of SEVN (Costa et al. 2023) to sample Pop III BBHs from Pop III star formation events. Each simulation follows the evolution of a large number (~1–2 × 10⁷) of Pop III binary stars randomly sampled from the given IMF (for the primary star) and distributions of initial binary parameters (mass ratio, orbital period, and eccentricity), based on the stellar evolution tracks of effectively metal-free stars (with an absolute metallicity of $Z = 10^{-11}$) generated by the PARSEC code (Bressan et al. 2012; Costa et al. 2021, 2022; Nguyen et al. 2022). The masses and orbital parameters of the resulting Pop III BBHs at the moment when both stars become BHs (and the binary remains bound), together with the initial masses, orbital parameters, and lifetimes of their progenitor stellar binaries are recorded in a catalogue. Here, the lifetime $\tau_B$ of a BBH progenitor stellar binary is the time taken for both stars to evolve from zero-age main sequence (ZAMS) to a BH. In this way, the simulation also predicts the mass fraction of binaries that become BBHs $f_{\text{BBH}}$ and the average initial mass of BBH progenitor stellar binaries $\bar{m}_{\text{p,BBH}}$, which are crucial parameters for our sampling process (see below).

In each Pop III star-forming halo, A-SLOTH predicts the expected mass $\delta M_{\star,\text{III}}$ of Pop III stars formed at the current time-step $t$. We assume that on average a fraction $f_B = 0.7$ of the mass of newly born Pop III stars is made up by binaries based on the N-body simulations of Pop III SCs in Liu et al. (2021a). Treating Pop III star formation as a stochastic<fnref n="8" /> process at the scales resolved by A-SLOTH, the number of Pop III BBHs (as well as the corresponding progenitor stellar binaries) to be sampled at this time-step $N_B$ is drawn from a Poisson distribution with parameter $f_{\text{boost}} f_B f_{\text{BBH}} \delta M_{\star,\text{III}} / \bar{m}_{\text{p,BBH}}$, where $f_{\text{boost}} = 50$ is a boost factor introduced to obtain better statistics. If $N_B > 0$, we randomly draw $N_B$ BBH progenitor binary stars from the input BPS catalogue, which will become BBHs at $t + \tau_B$.

Since the (initial/ZAMS) properties of Pop III binary stars are largely unknown in the lack of direct observations, to explore the parameter space, we consider six SEVN simulations (LOG1, TOP1, KRO1, LOG5, TOP5, and KRO5) with diverse initial conditions from Costa et al. (2023, see their section 2 for details), combining three IMF models and two sets of initial binary parameter distributions. In particular, we fix the mass range of primary stars as $m_{\star,1} \in [5, 550] M_{\odot}$ and consider three models for (the shape of) the primary star IMF as follows (see also fig. 4 of Costa et al. 2023).

(i) A log-flat distribution $p(m_{\star,1}) \equiv dN/dm_{\star,1} \propto m_{\star,1}^{-1}$ (labelled by ‘LOG’), which is motivated by the mass distributions of Pop III protostars in (magneto)hydrodynamic simulations (see e.g. fig 6 in Klessen & Glover 2023, and references therein) and chosen as the fiducial model.

(ii) A Kroupa (2001) IMF $p(m_{\star,1}) \propto m_{\star,1}^{-2.3}$ (labelled by ‘KRO’), which is consistent with the IMF of Pop I/II stars (for $m_{\star,1} > 0.5 M_{\odot}$) in local observations.

(iii) A top-heavy distribution $p(m_{\star,1}) \propto m_{\star,1}^{-0.17} e^{-m_{\star,1}^2/m_{\star,1}^2}$ given $m_{\text{cut}}^2 = 20 M_{\odot}^2$ (labelled by ‘TOP’) as an extreme case, which is adopted in some cosmological simulations (Jaacks et al. 2018; Jaacks, Finkelstein & Bromm 2019; Liu & Bromm 2020b, c).

In all cases, the IMF follows a power-law at the high-mass end. Therefore, each IMF model can also be characterized by the (asymptotic) power-law slope $\alpha$, with $\alpha = 1$ for LOG, 2.3 for KRO and 0.17 for TOP.

For IBS, we consider two distinct sets of distributions of the mass ratio $q_s \equiv m_{\star,2}/m_{\star,1}$, orbital period $P$, which is characterized by the variable $\pi \equiv \log(P/\text{day})$ in practice, and eccentricity $e$ (see also fig. 3 of Costa et al. 2023). A lower bound on the secondary mass $m_{\star,2} > 2.2 M_{\odot}$ is also imposed in both cases.

(i) The first model is based on observations of present-day massive binary stars (Sana et al. 2012, henceforth S12) with $p(q) \propto p^{-0.1}$ for $q \in [0.1, 1]$, $p(\pi) \propto \pi^{-0.55}$ for $\pi \in [0.15, 5.5]$, and $p(e) \propto e^{-0.42}$ for $e \in [0, 1)$. This model favours initially close binary stars and is labelled by the number ‘1’ following the convention in Costa et al. (2023). We call it the close IBS model from now on.

(ii) The second model adopts the power-law and log-normal fits to the mass ratio and orbital period distributions of Pop III protostars in the hydrodynamic simulations by Stacy & Bromm

<footnote n="8">⁸It is likely that star formation is self-regulated and not fully stochastic (Lewis et al. 2023; Yan, Jerabkova & Kroupa 2023; Liu et al. 2024a), especially considering Pop III systems made of small numbers of massive stars (but see also Grudic et al. 2023a). Since we are only concerned with the global population of Pop III BBH mergers in this paper, a detailed investigation of IMF sampling schemes (see e.g. Kroupa et al. 2013) is deferred to future work.</footnote>
(2013, hereafter SB13): $p(q) \propto q^{-0.55}$ for $q \in [0.1, 1]$, and $p(\pi) \propto \exp[-(\pi - \mu)^2/(2\sigma^2)]$, given $\mu = 5.5$ and $\sigma = 0.85$. The eccentricity distribution is assumed to be thermal, i.e. $p(e) = 2e$ for $e \in [0, 1]$, following previous BPS studies of Pop III BBHs (e.g. Kinugawa et al. 2014; Hartwig et al. 2016; Tanikawa et al. 2021b). This model is dominated by initially wide binaries and labelled by the number '5'. We refer to it as the wide IBS model henceforth.

As it is still challenging to simulate the entire Pop III star formation process with well-resolved disc fragmentation, protostellar feedback, and N-body dynamics, the properties of Pop III SCs are still unclear in current simulations, and so is the statistics of binary stars (see e.g. Liu et al. 2021a). We hope to bracket the reality with these two models, which are supported by different simulations with different setups (numerical scheme of hydrodynamics, resolution, modelling of feedback, and cosmological context). For instance, the simulations by Hirano et al. (2018b) show that it is possible to form highly eccentric, close binaries via violent N-body dynamics in compact Pop III SCs formed by mergers of multiple star-forming gas clumps in relatively large haloes ($M_h \sim 10^7 M_\odot$) under strong streaming motions between gas and dark matter. The formation of close binaries by dynamical friction is also favoured in dense gas structures with collision-induced emission cooling (Riaz et al. 2023), at least in the initial stage (a few thousand years after the formation of the first protostar) when protostellar feedback is expected to be unimportant (Jaura et al. 2022). On the other hand, the model dominated by wide binaries is motivated by the trend that Pop III protostar clusters/binaries formed in typical minihaloes ($M_h \sim 10^{5-6} M_\odot$) tend to expand due to accretion/inflow of gas with high angular momentum in (magneto-)hydrodynamic simulations of hot, thick star-forming discs (e.g. Chon & Hosokawa 2019; Heath & Nixon 2020; Sugimura et al. 2020, 2023; Franchini et al. 2023; Mignon-Risse, González & Commerçon 2023; Park et al. 2023, 2024).

Given $\delta M_{\star,\text{III}}$, A-SLOTH also samples individual Pop III stars from an input power-law IMF $p(m_\star) \propto m_\star^{-\alpha_s}$, which are treated as single stars to model stellar feedback, i.e. photoheating, ionization, SN-driven outflows, and metal enrichment (see section 2 of Hartwig et al. 2023). In principle, one can pair a fraction of the sampled stars into binaries and follow their evolution on the fly to predict the formation of BBHs (and meanwhile model the feedback from binary stars). However, as described above, we instead adopt a simple approach by sampling BBHs directly from the input SEVN catalogue to explore the parameter space of the poorly constrained properties of Pop III binary stars. In this way, the progenitor stellar binaries of Pop III BBHs sampled separately by our scheme in a star formation event are not consistent with the (underlying) population of Pop III stars sampled by the stellar feedback routine of A-SLOTH, and the impact of binary interactions on stellar feedback is also ignored. To approximately capture the variations of stellar feedback and BBH formation with IMF spontaneously, we set the power-law slope and mass range of the single-star IMF $p(m_\star)$ for stellar feedback to those of the primary star IMF $p(m_{\star,1})$ adopted by the input BPS catalogue, i.e. $\alpha_s = \alpha = 1$ for LOG (fiducial), 2.3 for KRO, and 0.17 for TOP,<fnref n="9" /> with $m_\star \in [5, 550] M_\odot$. We defer a more complete modelling of binary stars, particularly their feedback features such as enhanced

production of (ionizing) UV photons and different metal yields (e.g. Sansom, Izzard & Ocvirk 2009; Chen et al. 2015; Göttberg et al. 2019; Secunda et al. 2020; Ma et al. 2022; Tsai et al. 2023; Lecroq et al. 2024; Liu et al. 2024c; Yates et al. 2024), to future work. We expect the effects of binary interactions to be comparable to those of varying the IMF (see Section 5.2.3), which will be overwhelmed by the uncertainties in other aspects. For instance, the star formation history of Pop III stars is still highly uncertain in theory, showing orders of magnitude of discrepancies between different studies with different feedback models, resolution, and simulation volumes (see e.g. fig. 13 in Hartwig et al. 2022). The poorly understood initial binary properties (captured by the six SEVN models considered here) also lead to order-of-magnitude uncertainties in the merger rate of Pop III BBHs (see Section 5.2 and e.g. Santoliquido et al. 2023).

## 2.2 Evolution of Pop III BBHs

Once a Pop III BH binary is sampled, we follow on-the-fly the evolution of its orbit in the host galaxy or NSC as well as the internal evolution of the binary orbit.

### 2.2.1 Galactic dynamics

The dynamics of Pop III BBHs in host galaxies is driven by galaxy mergers, which randomize the (galactic/external) orbits of Pop III BBHs coming from the smaller progenitor halo, and dynamical friction (DF) by field stars and dark matter, through which massive Pop III BBHs sink towards galaxy centres where they can further fall into NSCs. We use the same method as in Liu & Bromm (2021, see their section 2.2 for details) to model such dynamics but further include the effects of dark matter in addition to stars. Here, we only consider the DF of stars and dark matter, ignoring the effects of gas. It is shown in simulations (e.g. Chen et al. 2022a) that DF is dominated ($\gtrsim 99$ per cent) by collisionless particles. Besides, for Pop III BBHs, which are typically much lighter than supermassive BH seeds and SCs, the DF time-scale is only shorter than the Hubble time in the central region ($r \lesssim 300$ pc) dominated by stars in most cases.

#### 2.2.1.1 Initial distribution
When a Pop III BBH is born or falls into a larger halo during a halo merger event, given the virial radius $R_{\text{vir}}$ of the (post-merger) host halo, its apocentre distance $r$ is drawn randomly from a (cumulative) probability distribution of $x \equiv r/R_{\text{vir}}$:

$$F(< x) \approx \begin{cases} 10^{-4}(1000x)^{2.23}, & x < 0.03, \\ 0.2(100x/3)^{0.88}, & 0.03 \le x < 0.1, \\ 0.58(10x)^{0.24}, & 0.1 \le x < 0.95, \end{cases} \quad (1)$$

which is a broken-power-law fit of the spatial distribution of Pop III remnants in high-z haloes predicted by the cosmological simulation FDbox_Lseed in Liu & Bromm (2020a, b). The eccentricity $e_*$ of a (galactic) orbit is initially drawn from a uniform distribution in $[0,1]$, following Arca-Sedda & Capuzzo-Dolcetta (2014). Here, we do not consider the effects of SN natal kicks during the formation of Pop III BBHs on their initial galactic orbits. We expect the initial orbits at birth to be unimportant for the inspiral of Pop III BBHs in galaxies hosting NSCs, as such galaxies reside in relatively massive haloes with $M_h \gtrsim 10^9 M_\odot$ in which most Pop III BBHs come from cosmic accretion flows and halo mergers, so their dynamics are mainly determined by halo assembly processes captured by the distribution of $x \equiv r/R_{\text{vir}}$ in cosmological simulations (Liu & Bromm 2020a, b). Besides, we assume that during a halo merger event, all Pop III BBHs in the smaller halo will be immediately stripped and redistributed in the post-merger halo according to equation (1).

<footnote n="9">⁹For the TOP model with $p(m_{\star,1}) \propto m_{\star,1}^{-0.17} e^{-m_{\star,1}^2/m_{\star,1}^2}$, we do not include the exponential cutoff term in the IMF used to sample single stars for stellar feedback for simplicity. We have verified by numerical experiments that the exact shape of the IMF at the low-mass ($\gtrsim 20 M_\odot$) end has negligible impact on the star formation history in this case since stellar feedback is dominated by massive stars for such a top-heavy IMF.</footnote>
### 2.2.1.2 Inspiral by dynamical friction
Within the host galaxy, we follow the inspiral of each binary with

$$ \begin{aligned} \frac{dr}{dt} & \approx -r \left[ \frac{1}{\tau_{\text{DF}}(r, e_*, m|\vec{\lambda}_*)} + \frac{1}{\tau_{\text{DF}}(r, e_*, m|\vec{\lambda}_\chi)} \right], \\ \frac{de_*}{dt} & \approx -e_* \left[ \frac{1}{\tau_{\text{DF}}(r, e_*, m|\vec{\lambda}_*)} + \frac{1}{\tau_{\text{DF}}(r, e_*, m|\vec{\lambda}_\chi)} \right], \end{aligned} \quad (2) $$

where $m \equiv m_1 + m_2$ is the total mass of the binary given the primary BH mass $m_1$ and secondary BH mass $m_2$, $\tau_{\text{DF}}$ is the DF time-scale formula (see below), $\vec{\lambda}_\chi$ captures the dark matter halo properties ($M_{\text{vir}}$ and $R_{\text{vir}}$), and $\vec{\lambda}_*$ denotes parameters of the background stellar system, which are the mass $M_*$, size $R_*$ and inner slope of density profile $\gamma_*$ of the host galaxy in our case. Here, $M_*$ is set to the total mass of Pop I/II stars predicted by the baryon cycle in A-SLOTH or an extrapolation scheme for the low-$z$ regime not covered by the merger trees (see Section 4). We update $r$ and $e_*$ at each star formation time-step of the host halo. Note that Pop III stars typically form in small clusters, which, if they remain bound, can sink into the galaxy centre more efficiently. The reason is that in this case, $m$ in equation (2) should be the total mass of the cluster ($\sim 10^3 M_\odot$), and $\tau_{\text{DF}}$ is shorter for higher $m$ (see below). By using the masses of individual binaries in equation (2), we make the conservative assumption that Pop III SCs have been dissolved by internal relaxation, mass-loss from stellar evolution, and tidal disruption during halo/galaxy mergers before falling into galaxies hosting NSCs.<fnref n="10" />

We use the minimum of the values predicted by two DF time-scale formulae to estimate $\tau_{\text{DF}}$. The first is the Chandrasekhar formula (Binney & Tremaine 2011)

$$ \tau_{\text{DF,C}} = \frac{342}{\ln \Lambda} \left(\frac{r}{3 \text{ pc}}\right)^2 \left(\frac{v}{10 \text{ km s}^{-1}}\right) \left(\frac{m}{100 M_{\odot}}\right)^{-1}, \quad (3) $$

where $\ln \Lambda \sim \ln[M_*r/(0.8mR_*)]$ is the Coulomb logarithm, and $v \sim \sigma_* \sim \sqrt{GM_*/(0.8R_*)}$. The second formula is a generalization of equation (3) and can be applied to both cusped and cored density profiles (Arca-Sedda et al. 2015; Arca-Sedda 2016):

$$ \begin{aligned} \frac{\tau_{\text{DF,A}}}{\text{Myr}} &= 0.3g \sqrt{\left(\frac{R_*^3}{1 \text{ kpc}^3}\right) M_{*11}^{-1} \left(\frac{m}{M_*}\right)^\alpha \left(\frac{r}{R_*}\right)^\beta}, \\ g &= (2 - \gamma_*) \left\{ a_1 \left[ \frac{1}{(2 - \gamma_*)^{a_2}} + a_3 \right] (1 - e_*) + e_{\text{if}} \right\}, \end{aligned} \quad (4) $$

given $\alpha = -0.67$, $\beta = 1.76$, $a_1 = 2.63 \pm 0.17$, $a_2 = 2.26 \pm 0.08$, and $a_3 = 0.9 \pm 0.1$. This formula is only valid for $0 < \gamma_* < 2$ and equation (3) is used otherwise.<fnref n="11" />

### 2.2.1.3 Host galaxy/halo model
To evaluate the first dynamical friction term driven by stars, we assume that the galaxy (stellar) density profile is described by a Dehnen sphere (see also equation 11 of Dehnen 1993), and the galaxy size (scale length) $R_*$ is derived from the empirical relation between galaxy mass and half-mass radius $R_{1/2}$ in local observations (see their equation 18 and fig. 1, Arca-Sedda &

Capuzzo-Dolcetta 2014):

$$ \begin{aligned} R_* &= [2^{1/(3-\gamma_*)} - 1]R_{1/2}, \\ R_{1/2} &= \max(0.02R_{\text{vir}}, 2.37M_{*11}^{0.14} \text{ kpc}), \end{aligned} \quad (5) $$

where $M_{*11} \equiv M_*/(10^{11} M_\odot)$, and the inner slope $\gamma_*$ is generated from a uniform distribution in the ranges of 0–2. We also impose a lower limit of $0.02R_{\text{vir}}$ on $R_{1/2}$ according to the results (for median galaxy half-mass radii) from sub-halo abundance matching (Somerville et al. 2018) for galaxies with $M_* \sim 10^8-10^{11} M_\odot$ at $z \le 3$.

In reality, galaxy size/morphology is regulated by a variety of processes, such as gas accretion, galaxy mergers, stellar, and AGN feedback (see the discussions in e.g. Hopkins et al. 2023; Karmakar, Genel & Somerville 2023), which introduce dispersion to the scaling relation. To capture such dispersion, we perturb the value of $\log(R_*)$ predicted by equation (5) with a random number (drawn once<fnref n="12" /> and inherited along the main branch<fnref n="13" />) following a uniform distribution in $[-0.36, 0.36]$ dex, motivated by the scatter in observations (Leigh, Böker & Knigge 2012). Besides, additional uncertainties may arise when extrapolating the local empirical scaling relation to the high-$z$ regime ($z \sim 3-10$) that is important for our purposes. For instance, it is found in recent cosmological simulations (Roper et al. 2023) that galaxies at $z \gtrsim 5$ are expected to follow a different size–mass relation with a negative slope, and are more compact, embedded in central concentration of non-star forming gas (i.e. ‘inside-out’ galaxy formation), which is also supported by JWST observations (e.g. Baggen et al. 2023; Langeroodi & Hjorth 2023; Ito et al. 2024; Ji et al. 2024a; Ormerod et al. 2024). Our $z$-independent estimate of $R_*$ can be regarded as an upper limit. It will be shown below (equation 4) that the dependence of DF time-scale $\tau_{\text{DF}}$ on galaxy size is weak ($\tau_{\text{DF}} \propto R_*^{0.26}$), which, combined with the weak mass and redshift dependence of $R_*$ in equation (5), implies that our results are insensitive to the modelling of $R_*(M_*)$. Actually, we find by numerical experiments that changing $R_*$ by a factor of 10 can only alter the number of Pop III BBHs (and mergers) in NSCs by up to $\sim 40$ per cent.

For the second term in equation (2) driven by dark matter, which is only important for Pop III BBHs at galaxy outskirts ($r \gtrsim 300$ pc), we only use the Chandrasekhar formula (equation 3) with $\ln \Lambda \sim \ln[M_{\text{vir}}r/(0.8mR_{\text{vir}})]$ and $v \sim \sqrt{GM_{\text{vir}}/(0.8R_{\text{vir}})}$, implicitly assuming a power-law slope $\gamma_\chi > 2$ for the dark matter density profile,<fnref n="14" /> which describes the asymptotic behaviour of the NFW profile at large radii. During a merger between two haloes, the orbits of binaries from the smaller halo are randomized with the same distributions adopted for initialization.<fnref n="15" />

We further introduce a free parameter $M_{*\text{min}}$, as the (Pop I/II stellar) mass threshold above which galaxies have regular morphology and well-defined dynamical centres to which massive objects can sink. DF is turned off in galaxies below this threshold. The morphology and dynamics of high-$z$ galaxies are still poorly

<footnote n="10">¹⁰Most Pop III stars form in small clusters with total masses of $\sim 100-10^4 M_\odot$ in our simulations. Assuming a half-mass radius of $\sim 1$ pc and an average stellar mass of $\sim 10 M_\odot$, the dissolution time-scale of such clusters by internal relaxation is $\sim 100$ Myr, similar to the time-scale $\sim 100$ Myr in which minihaloes hosting Pop III stars/remnants experience major mergers that bring Pop III remnants to large haloes ($M_h \gtrsim 10^9 M_\odot$) with NSCs.</footnote>

<footnote n="11">¹¹We have $\tau_{\text{DF,A}} < \tau_{\text{DF,C}}$ in most cases.</footnote>

<footnote n="12">¹²It is found in cosmological simulations that most galaxies oscillate around the median size several times in their lifetimes (see fig. 8 and 9 of Karmakar et al. 2023). We ignore such oscillations for simplicity.</footnote>

<footnote n="13">¹³In a merger tree, starting with the root, the main branch goes through the most massive child halo in each halo merger event (branching).</footnote>

<footnote n="14">¹⁴Apparently, this simple assumption ($\gamma_\chi > 2$) does not hold in the central region, where the dynamical friction by dark matter can be important for the dynamics of (normal) SCs (see Section 3.2). Nevertheless, we find that our results are insensitive to the dynamics of SCs.</footnote>

<footnote n="15">¹⁵For simplicity, we ignore possible enhancement of inspiral by DF during galaxy mergers implied by recent observations (Román et al. 2023).</footnote>
constrained observationally. The recent numerical experiments by Hopkins et al. (2023) found that the formation of a centrally concentrated density profile that defines the dynamical centre drives the transition from irregular to disc morphology, which mostly occurs around $M_* \sim 10^{10} M_\odot$, although in principle the transition can happen at any mass scale and less massive disc galaxies are also seen in simulations. In this study, we adopt $M_{*,min} = 10^6 M_\odot$ by default as an optimistic choice, which is approximately the minimum mass of galaxies hosting NSCs observed in the local Universe (Neumayer, Seth & Böker 2020). The motivation is that the observed properties of NSCs in local low-mass galaxies ($M_* \lesssim 10^9 M_\odot$), including both early- and late-type, can be well explained by the globular cluster (GC) accretion scenario (see e.g. Fahrion et al. 2022b; Leaman & van de Ven 2022; Román et al. 2023), in which multiple GCs inspiral to the centre by DF and build up the NSC (e.g. Arca-Sedda & Capuzzo-Dolcetta 2014; Tsatsi et al. 2017; Arca-Sedda, Kocsis & Brandt 2018; Fahrion et al. 2022a; Wang & Lin 2023), implying that DF must work in such NSC host galaxies.

### 2.2.2 Binary orbital evolution

The internal evolution of a BBH orbit starts at $\tau_B$ after the initial formation of the progenitor binary stars. The earlier stage ($t < \tau_B$) is governed by binary stellar evolution whose outcome is given by the input BPS catalogue and serves as the initial condition for BBH evolution in A-SLOTH. We use the same set of differential equations detailed in Liu & Bromm (2021, see their section 2.3 and below) to follow the internal evolution of BBHs in the field and NSCs as well as the inspiral of BBHs in NSCs. Here, we integrate the equations for evolution inside SCs on the fly with the explicit Euler method in A-SLOTH and adopt an approximation based on the merger time-scale formula from Iorio et al. (2023) to facilitate the computation for the internal evolution of binaries in the field. When numerical integration<fnref n="16" /> is used to evolve a and e in SCs, we further do sub-cycling with time-steps $\delta t_{\text{binary}} \le 0.01 a/(da/dt)$. A merger happens when the expected time taken to merge has passed since BBH formation or $a < 6G(m_1 + m_2)/c^2$ according to numerical integration. Then the merger properties are recorded and the binary is removed from the data structure in A-SLOTH. This approach is more accurate and flexible than post-processing since the evolution of host SCs is self-consistently taken into account (see Section 3).

#### 2.2.2.1 Isolated evolution in galaxy fields
In the field, we assume that binary evolution is only driven by GW emission<fnref n="17" /> (Peters 1964):

$$ 
\begin{aligned}
\dot{a}(a, e) &\equiv \frac{da}{at} = \frac{da}{dt}\bigg|_{\text{GW}} = -\frac{64}{5} \frac{A}{a^3(1-e^2)^{7/2}} f_1(e), \\
\dot{e}(a, e) &\equiv \frac{de}{dt} = \frac{de}{dt}\bigg|_{\text{GW}} = -\frac{304}{15} e \frac{A}{a^4(1-e^2)^{5/2}} f_2(e),
\end{aligned}
\quad (6) $$

where $A = G^3 m_1 m_2 (m_1 + m_2)/c^5$, and

$$ f_1(e) = 1 + \frac{73}{24}e^2 + \frac{37}{96}e^4, \quad f_2(e) = 1 + \frac{121}{304}e^2. \quad (7) $$

An approximated solution of these equations is given by Iorio et al. (2023, see their appendix C) based on Zwick et al. (2020) that predicts the time taken to merge from the initial BBH condition ($a_i$ and $e_i$) as

$$ 
\begin{aligned}
t_{\text{GW}}(a_i, e_i) &= \begin{cases} t_{\text{GW},0}(a_i)/[1 + f_{\text{corr}}(e_i)], & e_i < 0.999, \\ t_{\text{peters}}(a_i)8^{1-\sqrt{1-e_i}}, & e_i \ge 0.999, \end{cases} \\
t_{\text{GW},0}(a_i) &= \frac{[5a_i^4(1-e_i^2)^{7/2}]}{(256A)}, \\
t_{\text{peters}}(a_i) &= t_{\text{GW},0}/f_1(e_i), \\
f_{\text{corr}} &= e_i^2 [-0.443 + 0.58(1-e_i^{3.074})^{f_3(e_i)}], \\
f_3(e_i) &= 1.105 - 0.807e_i + 0.193e_i^2,
\end{aligned}
\quad (8) $$

which is accurate enough for our purposes. Therefore, instead of integrating equations (6), we set a clock for each BBH at formation and trigger the merger after $t_{\text{GW}}$ has passed. If a BBH falls into a NSC and is ejected later (either by binary-single encounters or disruption of the SC, see below), we update the clock on an individual basis with a new merger time-scale $t_{\text{GW},ej}(a_{ej}, e_{ej})$ calculated from the orbital parameters ($a_{ej}$ and $e_{ej}$) at the moment of ejection, which have evolved inside the SC (see below).

#### 2.2.2.2 Dynamical evolution in NSCs
A binary falls into an NSC when the apocentre distance of its (galactic) orbit is smaller than the characteristic size (i.e. scale length) of the NSC in the host galaxy (i.e. $r < R_{\text{SC}}$). Once inside SCs, soft binaries will be destroyed and hard binaries will be hardened by binary-single encounters. The criterion for hard binaries is $a < R_{\text{SC}}$ and $a < a_{\text{HDB}} = Gm_1m_2/[m_{*,\text{SC}}\sigma_*(r)^2]$ (Mapelli et al. 2021), given the SC size $R_{\text{SC}}$, local velocity dispersion $\sigma_*(r)$, and typical mass of objects in the cluster. Throughout this study, we adopt $m_{*,\text{SC}} = 0.5 M_\odot$ for simplicity, although $m_{*,\text{SC}}$ may vary from cluster to cluster and evolve with redshift in reality. The hard binary criterion is checked at each (star formation) time-step and soft binaries are immediately removed. For simplicity, here we have ignored the hardening of BBHs by interactions with gas in NSCs, which can even shrink the orbits of initially soft binaries with $a > a_{\text{HDB}}$ to make them hard, shifting the soft-hard boundary to larger separations (Rozner & Perets 2024).

For the surviving hard binaries, we include additional binary-single encounter terms in the binary evolution equations (Sesana & Khan 2015; Arca Sedda 2020; Mapelli et al. 2021):

$$ \frac{da}{dt} = \frac{da}{dt}\bigg|_{\star} + \frac{da}{dt}\bigg|_{\text{GW}}, \quad \frac{de}{dt} = \frac{de}{dt}\bigg|_{\star} + \frac{de}{dt}\bigg|_{\text{GW}}, \quad (9) $$

$$ \frac{da}{dt}\bigg|_{\star} = -\frac{GH\rho_{\star}(r)}{\sigma_{\star}(r)}a^2, \quad \frac{de}{dt}\bigg|_{\star} = \kappa \frac{GH\rho_{\star}(r)}{\sigma_{\star}(r)}a, \quad (10) $$

where $\rho_*(r)$ is the local stellar density of the cluster, $\mathcal{H} \sim 1-20$ and $\kappa \sim 0.01-0.1$ (Sesana, Haardt & Madau 2006; Sesana & Khan 2015; Mapelli et al. 2021) are two dimensionless factors. We adopt $\mathcal{H} = 20$ following Liu & Bromm (2021) and fix $\kappa = 0.01$ as a conservative assumption. We also impose a lower limit of $10^{-8}$ on $1-e$ under the optimistic assumption that the binary will merge immediately when driven to this limit by binary-single interactions.

<footnote n="16">¹⁶At each time-step $\Delta t_i$, we first update a by $\Delta t_i/2$: $a_{i+1/2} = a_i + \dot{a}(a_i, e_i)\Delta t_i/2$. Then, we evolve e for the whole step: $e_{i+1} = e_i + \dot{e}(a_{i+1/2}, e_i)\Delta t_i$. Finally, we assign $a_{i+1} = a_{i+1/2} + \dot{a}(a_{i+1/2}, e_{i+1})\Delta t_i/2$.</footnote>

<footnote n="17">¹⁷BH binaries with highly eccentric ($e \gtrsim 0.99$), wide ($a \sim 10^4$ au) orbits (and wide BH triples) can also be driven to merge by perturbations from fly-by field stars (Michaely & Perets 2019, 2020). This channel can be important for Pop III BBH mergers if the Pop III IBS is dominated by wide binaries and Pop III stellar/BH triples are common (see e.g. fig. 8 of Sharda et al. 2020). We plan to consider this channel in the future.</footnote>
As in Arca Sedda (2020), here a SC of mass $M_{\text{SC}}$ is characterized by a cored Dehnen sphere (Dehnen 1993) with

$$ 
\begin{align} 
\rho_{*}(r) &= \begin{cases} \frac{(3-\gamma_{\text{SC}})M_{\text{SC}}/(4\pi R_{\text{SC}}^3)}{\chi^{\gamma_{\text{SC}}}(1+x)^{4-\gamma_{\text{SC}}}}, & r \ge r_{c}, \\ \rho_{*}(r_{c}), & r < r_{c}, \end{cases} \tag{11} \\ 
\sigma_{*}(r) &= \begin{cases} \sqrt{\frac{G M_{\text{SC}}}{R_{\text{SC}}}} \left(\frac{r}{R_{\text{SC}}}\right)^{\delta/2}, & r_{c} \le r < R_{\text{SC}}, \\ \sigma_{*}(r_{c}), & r < r_{c}, \end{cases} \nonumber 
\end{align} 
$$

where $x \equiv r/R_{\text{SC}}$, $\delta = \gamma_{\text{SC}}$ for $\gamma_{\text{SC}} < 1$ and $\delta = 2 - \gamma_{\text{SC}}$ for $\gamma_{\text{SC}} \ge 1$, given the scale length $R_{\text{SC}}$ and inner slope $\gamma_{\text{SC}}$ of the Dehnen sphere, which are derived from $M_{\text{SC}}$ using the empirical scaling relations in Liu & Bromm (2021) based on local observations (Neumayer et al. 2020; Pechetti et al. 2020). In this model, most hardening happens in the core, and the hardening efficiency is sensitive to the core overdensity parameter $\Delta_c$ that determines the core size $r_c = \max(R_{\text{SC}}/\Delta_c^{1/\gamma_{\text{SC}}}, r_*)$, which is no smaller than the influence radius of the BBH $r_* = R_{\text{SC}}(2m/M_{\text{SC}})^{1/(3-\gamma_{\text{SC}})}$. We adopt an optimistic core overdensity $\Delta_c = 100$. The external orbit of the binary also evolves by DF within the cluster as described by equation (2) in which the background stellar system now becomes the SC. In our case, $\rho_*(r)$, $\sigma_*(r)$, and $\tau_{\text{DF}}(r)$ are all updated at each star formation time-step to capture the evolution of the host SC (see Section 3).

Similar to Liu & Bromm (2021), we assume spherical symmetry and ignore higher-order processes involving tidal fields, general relativity effects, and interactions with binary Pop I/II BHs, central massive BHs, and gas, such as relativistic phase space diffusion, tides-driven eccentricity excitation, the Kozai-Lidov mechanism, hierarchical and repeated mergers, dynamical hardening/softening and accretion in AGN discs (see e.g. Hoang et al. 2018; Antonini, Gieles & Gualandris 2019; Yang et al. 2019; Arca Sedda 2020; McKernan et al. 2020; Tagawa, Haiman & Kocsis 2020; Mapelli et al. 2021, 2022; Zhang et al. 2021; Fragione et al. 2022; Atallah et al. 2023; Chattopadhyay et al. 2023; Chen & Lin 2023; DeLaurentiis, Epstein-Martin & Haiman 2023; Fragione & Rasio 2023; Li, Lin & Yuan 2023; Rowan et al. 2023; Rozner, Generozov & Perets 2023; Wang, Ma & Wu 2023a; Winter-Granić, Petrovich & Peña-Donaire 2024; Balberg 2024; Dall'Amico et al. 2024; Fabj & Samsing 2024; Gangardt et al. 2024; Hamilton & Rafikov 2024; Kritos et al. 2024; Purohit et al. 2024; Rozner & Perets 2024; Torniamenti et al. 2024; Trani, Quaini & Colpi 2024; Vaccaro et al. 2024). Here, interactions with existing Pop I/II BBHs are expected to be important in high-z metal-poor SCs (Barber et al. 2024), which we plan to take into account in future work. For completeness, we also consider ejections of Pop III binaries by binary-single encounters, which occur when $a_{\text{GW}} < a < a_{\text{ej}}$, given the critical separations for GW-dominance $a_{\text{GW}}$ and ejection $a_{\text{ej}}$ (Miller & Lauburg 2009; Antonini & Rasio 2016; Fragione & Silk 2020; Mapelli et al. 2021):

$$ a_{\text{GW}} = \left[ \frac{64}{5H} \frac{A\sigma_* f_1(e)}{G\rho_*(1 - e^2)^{7/2}} \right]^{1/5}, \quad (12) $$

$$ a_{\text{ej}} = \frac{Hm_{*}^{2} M_{\text{SC}}}{\pi (m_{1} + m_{2})^{3}} \frac{Gm_{1}m_{2}}{v_{\text{esc}}^{2}}, \quad (13) $$

$$ \frac{v_{\text{esc}}}{\text{km s}^{-1}} = 40 \left( \frac{M_{\text{SC}}}{10^5 M_{\odot}} \right)^{1/3} \left[ \frac{\rho_{*}(R_{\text{SC}})}{10^5 M_{\odot} \text{pc}^{-3}} \right]^{1/6}, \quad (14) $$

where $\rho_*(R_{\text{SC}})$ is the typical stellar density evaluated at the scale length $R_{\text{SC}}$ of the cluster. Once a binary is ejected from a cluster, its subsequent evolution is only driven by GWs. Besides, an ejected binary can no longer sink towards the galaxy centre by DF until its host halo merges into a larger halo and a new orbit is assigned to this binary in the post-merger system. The ejection criterion can

only be satisfied by low-mass BBHs ($m \lesssim 20 M_{\odot}$) in low-mass ($M_{\text{SC}} \lesssim 10^5 M_{\odot}$) SCs where $a_{\text{ej}} > a_{\text{GW}}$. In our case, the dynamical hardening time-scale in such low-mass ($M_{\text{SC}} \lesssim 10^5 M_{\odot}$) SCs is typically comparable to the Hubble time, and light BBHs in low-mass SCs are rare, so ejection by binary-single encounters is unimportant.

Since equation (6) is not integrated explicitly for BBHs in galaxy fields, integration of equations (9) starts from the initial conditions at the moment of BBH formation rather than infall. That is to say, we make an approximation for the initial condition of evolution inside SCs, which is justified by the fact that the progenitors of most mergers in SCs seldom evolve in fields and the evolution inside SCs is insensitive to the initial condition. This approximation does not affect the fate of binaries in SC, since the maximum separation for a BBH to be a hard binary inside SCs ($a \sim 10^2-10^4$ au) is much larger than that required for efficient evolution in fields by GW emission ($a \lesssim 0.5$ au). We have verified with numerical experiments that the difference between the results from explicit integration of equations (6) and the approximated merger time-scale solution (equations 8) is negligible.

## 3 NUCLEAR STAR CLUSTERS

For simplicity, we only consider NSCs and their descendants (e.g. from disrupted satellite galaxies during galaxy mergers) as the potential sites that can provide efficient dynamical hardening of Pop III binaries. All SCs considered in this work are initially NSCs although later on some of them become normal SCs (see Section 3.2) either as remnant NSCs from stripped satellite galaxies or NSCs disrupted by internal processes captured by the NSC occupation fraction (equation 15).

Since Pop III star formation typically peaks at $z \sim 10-20$ as predicted by cosmological simulations (e.g. Tornatore, Ferrara & Schneider 2007; Wise et al. 2012; Johnson et al. 2013; Muratov et al. 2013; Pallottini et al. 2014; Xu et al. 2016; Sarmento et al. 2017; Liu & Bromm 2020b), we have to consider NSCs across the entire cosmic history of galaxy formation (particularly at $z \gtrsim 1$) to model their interactions with Pop III BBHs. However, the properties of NSCs (and SCs in general) beyond the local volume and nearby galaxy clusters are poorly constrained in observations. Although the formation and evolution of (high-z) SCs (in dwarf galaxies) have been intensively studied with semi-analytical models, N-body and (cosmological) hydrodynamic (zoom-in) simulations (e.g. Devecchi & Volonteri 2009; Devecchi et al. 2010, 2012; Renaud, Bournaud & Duc 2015; Guillard, Emsellem & Renaud 2016; Safranek-Shrader et al. 2016; Li et al. 2017, 2022; Brown, Gnedin & Li 2018; Howard, Pudritz & Harris 2018; Li, Gnedin & Gnedin 2018; Pfeffer et al. 2018; Li & Gnedin 2019; Lahén et al. 2020, 2023; Ma et al. 2020; Chen, Li & Vogelsberger 2021; McKenzie & Bekki 2021; Fahriion et al. 2022b; Hislop et al. 2022; Grudić et al. 2023b; Livernois, Vesperini & Pavlík 2023; Sameie et al. 2023; van Donkelaar et al. 2023, 2024; Arca Sedda et al. 2024a, b; Chen, Mo & Wang 2024; Gao et al. 2024; Gray et al. 2024), our understanding of high-z SCs is still limited due to theoretical and numerical problems (see e.g. Chen et al. 2021; Hislop et al. 2022).

In the absence of a robust and universal theoretical model for NSCs in the broad redshift range ($z \sim 0-30$) involved in our simulations, we extrapolate the empirical scaling relations and occupation fraction of NSCs in local observations (Neumayer et al. 2020; Pechetti et al. 2020) to high redshifts, as detailed in section 2.3 of Liu & Bromm (2021). In this work, we incorporate their NSC model into the merger trees of A-SLOTH to keep track of NSCs on the fly (Section 3.1). We also model the dynamics, mergers, and evaporation of normal
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Same as Fig. 1 but for the NSC routine. Here, $M_{\text{SC}}$, $R_{\text{SC}}$, and $\gamma_{\text{SC}}$ denote the mass, size, and density profile inner slope of the SC characterized by a cored Dehnen sphere (equation 11). The orbit of a normal SC in the host halo is described by the apocentre distance $r$ and eccentricity $e_\star$. The green boxes denote the key input physics that evolve the SC properties (cyan boxes), through two pathways: ‘NSC formation/evolution’ and ‘normal SCs’ (white boxes), detailed in Sections 3.1 and 3.2, respectively. The orbits of the SCs from the smaller halo during a halo merger event are randomized according to equation (16).</figcaption></figure>

SCs as NSC descendants (Section 3.2). An illustration of our NSC routine is shown in Fig. 2. For simplicity, we have ignored the effects of Pop III BBHs on SC evolution (e.g. heating via binary-single encounters), which are expected to be small because Pop III BBHs are completely sub-dominant in the NSCs made of Pop I/II stars considered here (see Section 5.1.4). Building a physically motivated model for NSCs, other populations of SCs (e.g. GCs and young star clusters) and SMBHs that are important for NSC formation and evolution in A-SLOTH is an intriguing direction for future research (see e.g. El-Badry et al. 2019; Fahrion et al. 2022a, b; Park et al. 2022b; Polkas et al. 2024; Chen et al. 2024; De Lucia et al. 2024; Gao et al. 2024; Kaur, Rom & Sari 2024).

## 3.1 Formation and evolution of nuclear star clusters

In our model, formation of NSCs (in galaxies with $M_\star > M_{\star,\text{min}}$) is governed by the input occupation fraction of NSCs as a function of galaxy mass $f_{\text{occ}}(M_\star) \equiv f_{\text{occ}}(\tilde{M}_\star)$, $\tilde{M}_\star \equiv \log(M_\star/M_\odot)$. We adopt

the fit formula in Liu & Bromm (2021, see their fig. 8)

$$f_{\text{occ}} \approx \begin{cases} 0.15, & \tilde{M}_\star < 6.6, \\ 0.15 + 0.36(\tilde{M}_\star - 6.6), & 6.6 \le \tilde{M}_\star < 8.2, \\ 0.72, & 8.2 \le \tilde{M}_\star < 9.5, \\ 0.72 - 0.185(\tilde{M}_\star - 9.5), & 9.5 \le \tilde{M}_\star, \end{cases} \quad (15)$$

for our fiducial model (labelled by ‘-obs.’), based on observations of nearby galaxies in all environments<fnref n="18" /> (Neumayer et al. 2020). We further impose a low-mass truncation at $M_{\star,\text{min}} = 10^6 M_\odot$, i.e. $f_{\text{occ}} = 0$ for $M_\star \le M_{\star,\text{min}}$, consistent with our modeling of DF. Here, we also ignore the redshift evolution of $f_{\text{occ}}$ for simplicity.

In the GC accretion scenario of NSC formation relevant for low-mass ($M_\star \lesssim 10^9 M_\odot$) galaxies, $f_{\text{occ}}$ is sensitive to galaxy sizes, so the occupation fraction can be higher if the majority of NSC formation happens before the dominant epoch of size growth at $z \sim 2$ (see their fig. 6 Leaman & van de Ven 2022). In fact, galaxies tend to be more compact at higher $z$ according to recent JWST observations (e.g. Baggen et al. 2023; Langeroodi & Hjorth 2023; Ito et al. 2024; Ji et al. 2024a; Ormerod et al. 2024) and cosmological (zoom-in) simulations (e.g. Roper et al. 2023), which support the inside-out scenario of galaxy formation. Besides, it is found in the Romulus cosmological simulations that dwarf galaxies formed earlier are more likely to host massive BHs (Tremmel et al. 2024), which is consistent with the enhanced NSC occupation observed in galaxy clusters. These outcomes imply that the dense environments at high-$z$ favour BH and NSC formation. On the other hand, it is also likely that high-$z$ galaxies are typically irregular and clumpy where DF is ineffective, as predicted by e.g. the simulations in Ma et al. (2021) for galaxies with $M_\star \lesssim 10^{10} M_\odot$ at $z \gtrsim 5$ (see also Biernacki, Teyssier & Bleuler 2017; Pfister et al. 2019; Bortolas et al. 2020).<fnref n="19" />

In light of this, in addition to the fiducial model, we also consider two extreme cases with higher and lower NSC occupation fractions. The three NSC models are summarized as follows:

(i) In the fiducial model (labelled by ‘obs.’), we have $f_{\text{occ}} = f_{\text{occ}}$ and $M_{\star,\text{min}} = 10^6 M_\odot$ based on local observations.

(ii) In the optimistic model (labelled by ‘full’), we consider full occupation with $f_{\text{occ}}(M_\star) \equiv 1$ for $M_{\star,\text{min}} = 10^6 M_\odot$.

(iii) In the conservative model (labelled by ‘low’), we reduce the efficiency of DF<fnref n="20" /> (as well as NSC formation) with $f_{\text{occ}} = f_{\text{occ}}$ and $M_{\star,\text{min}} = \min\{1, [(1+z)/6]^{5.14}\} \times 10^{10} M_\odot$, assuming that $M_{\star,\text{min}}$ evolves from $10^{10} M_\odot$ at $z = 5$ to $10^6 M_\odot$ at $z = 0$ following a power law of $(1+z)$.

<footnote n="18">¹⁸ $f_{\text{occ}}$ is generally higher in denser environments (Hoyer et al. 2021; Neumayer et al. 2020; Zanatta et al. 2024). For instance, the $f_{\text{occ}}$ measured in the Virgo cluster is higher than that in the local ($\lesssim 11$ Mpc) volume by $\sim 20$ per cent at $M_\star \sim 10^7-10^9 M_\star$ (Hoyer et al. 2021). Our fiducial model $f_{\text{occ}}$ may slightly overestimate the cosmic average by including observations of galaxies in all environments.</footnote>

<footnote n="19">¹⁹ Recent JWST observations find that (purely) irregular galaxies do not dominate at $z \gtrsim 3$ (especially for $z \lesssim 5$) as previously thought (Ferreira et al. 2023; Kwartalte et al. 2023). In particular, disc galaxies dominate the galaxy populations with $M_\star \gtrsim 10^9 M_\odot$ at $z \sim 0.6-8$, while the fraction of irregular galaxies remains below $\lesssim 20$ per cent, consistent with cosmological simulations (Lee et al. 2021, 2024b; Park et al. 2022a). However, one should be cautious with the interpretations of observational data in comparison with galaxy formation simulations (see e.g. Vega-Ferrero et al. 2024).</footnote>

<footnote n="20">²⁰ In this scenario, NSCs can still form without DF in ‘low-mass’ galaxies ($M_\star \lesssim 10^{10} M_\odot$) at high $z$ via *in situ* star formation (e.g. Brown et al. 2018; Guillard et al. 2016; Howard et al. 2018; Gray et al. 2024). However, without DF, Pop III BBHs can hardly fall into such NSCs, so they are irrelevant to our work. Therefore, we do not track NSCs in galaxies with $M_\star < M_{\star,\text{min}}$.</footnote>
Given the input $f_{\text{occ}}(M_*)$, we first generate a random number from a uniform distribution in [0,1] for each leaf of the merger tree. This random number is then inherited along the main branch such that every node $i$ has a random number $p_i$. At each (star formation) time-step, we check the criterion $p_i < f_{\text{occ}}(M_{*,i})$, which indicates the presence of an NSC, where $M_{*,i}$ is the current galaxy mass of node $i$. Given $p_i < f_{\text{occ}}(M_{*,i})$, if this node contains no SCs, we create a new SC object as its NSC, whose properties (i.e. mass $M_{\text{SC}}$, size $R_{\text{SC}}$, and inner slope of density profile $\gamma_{\text{SC}}$, see equation 11) are derived from the empirical NSC–galaxy scaling relations,<fnref n="21" /> detailed in section 2.3 of Liu & Bromm (2021) based on the observational data compiled in Neumayer et al. (2020) and Pechetti et al. (2020), given the host galaxy stellar mass $M_{*,i}$ as the input.<fnref n="22" /> If the node already contains one or multiple SC(s) but does not host a NSC in the previous time-step, we convert the oldest SC (remnant) along the main branch into the NSC with properties also specified by the NSC–galaxy scaling relations, regardless of the status of this SC (see Section 3.2 below). If the node already hosts an NSC, we simply update the NSC properties according to the scaling relations that reflect the variation of $M_{*,i}$, assuming adiabatic evolution of NSCs regulated solely by host galaxy mass. A NSC host galaxy will lose its NSC if $p_i > f_{\text{occ}}(M_{*,i})$, which will only happen for $M_* > 10^{9.5} M_{\odot}$. In this case, the original NSC is not removed but turned into a normal SC that no longer follows the NSC–galaxy relations and will be evaporated and redistributed during galaxy mergers (see below). In reality, such internal disruption of NSCs can be caused by merging with supermassive black holes, which is more common in more massive (elliptical) galaxies (Neumayer et al. 2020), as reflected in the decrease of $f_{\text{occ}}$ with $M_*$ for $M_* > 10^{9.5} M_{\odot}$.

In our model, any snapshot of the merger tree(s) satisfies (statistically) the input occupation fraction and NSC–galaxy scaling relations. During a halo/galaxy merger, the parent node inherits the SCs from child nodes and only one of them can be the NSC of the post-merger galaxy. The rest are normal SCs, i.e. NSCs from stripped satellite galaxies (Wang et al. 2023b), which are subject to orbital redistribution, DF, tidal forces, and internal relaxation heating, so they can merge with NSCs or evaporate, as discussed below.

## 3.2 Dynamics and evolution of normal star clusters

Similar to Pop III BBHs, the orbit of a normal SC of mass $M_{\text{SC}}$ also shrinks and its eccentricity $e_*$ decays under DF according to equation (2), now replacing $m$ with $M_{\text{SC}}$. To integrate these equations in A-SLOTH, we must specify the conditions after galaxy mergers. We assume that the orbits of SCs coming from the larger (target) halo remain unchanged by the merger while the SCs from the smaller (satellite/infalling) halo are redistributed in the post-merger halo following a cored power-law profile:

$$n_{\text{SC}}(r) = \begin{cases} n_0, & r \le R_{\text{core}}, \\ n_0(r/R_{\text{core}})^{-3}, & r > R_{\text{core}}, \end{cases} \quad (16)$$

where $n_0$ is fixed by $4\pi \int_0^{R_{\text{vir}}} n_{\text{SC}}(r)r^2 dr = N_{\text{SC}}$, given the halo virial radius $R_{\text{vir}}$ and number of (satellite) SCs $N_{\text{SC}}$. Also, $R_{\text{core}} = \min(31.62 R_* M_{*,11}^{1/6}, R_{\text{vir}})$ is the core size, which is the minimum of the halo virial radius and total galaxy radius given by the size–mass relation in Arca-Sedda & Capuzzo-Dolcetta (2014, see their equations 18 and 19) with a scatter of 0.36 dex. An inner core of SCs with a uniform distribution is supported by the GC accretion model in Arca-Sedda & Capuzzo-Dolcetta (2014) that can reproduce the observed relations between NSC mass, galaxy mass, and velocity dispersion (Erwin & Gadotti 2012; Leigh et al. 2012; Scott & Graham 2013), while the power-law form of the outer region with a slope of –3 is motivated by the asymptotic feature of the NFW profile at $r \to \infty$, assuming that SCs and dark matter follow similar spatial distributions at halo outskirts.

Now for each SC from the smaller halo, we draw its apocentre distance $r$ in the post-merger halo from a probability distribution $p(r) = 4\pi n_{\text{SC}}(r)r^2/N_{\text{SC}}$. The eccentricity $e_*$ is again generated from a uniform distribution in [0,1) (Arca-Sedda & Capuzzo-Dolcetta 2014). This orbit randomization process for infalling SCs is supported by N-body simulations (e.g. Pagnini et al. 2023). Here, we assume instantaneous disruption of infalling galaxies so it is always the SC mass that is used to evaluate equation (2), even for the NSC of the infalling galaxy.

Since our modelling is highly simplified under spherical symmetry, we only consider mergers between NSCs and normal SCs, not among normal SCs themselves. A merger happens when a SC gets too close to an NSC with $r < R_{\text{SC}}$. Immediately after the merger, the Pop III BBHs contained in the normal SC are placed at $r = R_{\text{SC}}$ and their external orbital eccentricities $e_*$ are set to that of the infalling SC. For simplicity, we do not model the effects of individual merger events on NSC properties. We assume that the long-term effects of mergers are already captured by the empirical scaling relations (Neumayer et al. 2020) that govern NSC evolution in our model.

In addition to dynamics, we also calculate the mass-loss and expansion of normal SCs under external tidal fields and internal relaxation by<fnref n="23" />

$$\frac{dM_{\text{SC}}}{dt} \approx -\frac{M_{\text{SC}}}{t_{\text{dis}}}, \quad (17)$$

$$\frac{dR_{\text{SC}}}{dt} \approx \frac{2}{3} \frac{R_{\text{SC}}}{t_{\text{rh}}}, \quad (18)$$

given the (half-mass) two-body relaxation time-scale $t_{\text{th}}$, and the overall dissolution time-scale $t_{\text{dis}}$ which is estimated with the minimum of the two time-scales of internal and external effects (Portegies Zwart, McMillan & Gieles 2010):

$$t_{\text{rh}} = 200 \text{ Myr} \left( \frac{M_{\text{SC}}}{10^6 M_{\odot}} \right)^{1/2} \left( \frac{R_{\text{SC}}}{\text{pc}} \right)^{3/2} \left( \frac{m_{*,\text{SC}}}{M_{\odot}} \right)^{-1} \frac{10}{\ln \Lambda'}, \quad (19)$$

$$t_{\text{dis}} = \min(t_{\text{dis,relax}}, t_{\text{dis,tidal}}), \quad t_{\text{dis,relax}} \approx 30 t_{\text{rh}}, \quad t_{\text{dis,tidal}} \approx 2 \text{ Myr} \\ \times \left(\frac{N}{\ln \Lambda'}\right)^{3/4} \left(\frac{r_g}{\text{kpc}}\right) \left(\frac{v_g}{220 \text{ km s}^{-1}}\right)^{-1} (1 - e_*), \quad (20)$$

where $\ln \Lambda' = 0.11 N$, $r_g = \max(r, R_{\text{SC}})$, $v_g = \sqrt{GM_*/(0.8R_*)}$ is the typical circular velocity of the host galaxy of mass $M_*$ and scale

<footnote n="21">²¹More recent results for the (local) NSC–galaxy scaling relations are available in Hannah et al. (2024), which show higher NSC occupation fractions ($f_{\text{occ}} \simeq 1$ for $M_* \sim 10^{9.5}–10^{9.5} M_{\odot}$) than those adopted in our fiducial model based on Neumayer et al. (2020). Their impact on our final results is expected to be within the range covered by the three NSC models considered here. We plan to consider these new relations in future work.</footnote>

<footnote n="22">²²The scatter in NSC–galaxy scaling relations is also implemented with random variables inherited along main branches, similar to the NSC occupation and scatter of galaxy size.</footnote>

<footnote n="23">²³It is possible to calculate the mass-loss of SCs more accurately based on simulation results (see e.g. Madrid et al. 2017; Gieles & Gnedin 2023), which is particularly important for the formation and growth of NSCs by SC mergers. Since modelling such processes is beyond the scope of our study, we adopt this simple equation to capture the long-term decrease of dynamical hardening efficiency during SC relaxation.</footnote>
length $R_*$ (see their equation 18 Arca-Sedda & Capuzzo-Dolcetta 2014), and $N = M_{\text{SC}}/m_{*,\text{SC}}$ is the total number of stars in the cluster. In equation (18), we only consider long-term expansion driven by internal heating, which leads to an asymptotic solution $R_{\text{SC}} \propto t^{2/3}$ at $t \to \infty$ (Portegies Zwart et al. 2010), if ignoring mass-loss, given $t_{\text{th}} \propto R_{\text{SC}}^{3/2}$. For simplicity, we assume that $\gamma_{\text{SC}}$ remains constant. A SC evaporates when $M_{\text{SC}}$ drops below $100 M_{\odot}$. At this moment, the Pop III remnants contained in it inherit the orbit of the SC in the galaxy. All NSC properties are updated in time-steps no larger than $\delta t_{\text{NSC}} = 1$ Myr, and sub-cycling is introduced within star formation time-steps if needed.

## 4 LOW-REDSHIFT EXTRAPOLATION

Cosmological simulations that resolve Pop III star-forming minihaloes in large representative volumes down to $z=0$ are still prohibitively expensive. Therefore, cosmological simulations resolving minihaloes typically have small volumes and are only representative at high redshifts. For instance, the simulation data from Ishiyama et al. (2016) adopted in our A-SLOTH runs provides merger trees for every halo in a co-moving volume of $(8 h^{-1} \text{Mpc})^3$ with a mass resolution of $5000 h^{-1} M_{\odot}$ down to $z_f \simeq 4.5$ when the box is marginally representative. We design an extrapolation scheme to follow the evolution of Pop III BBHs and their host galaxies/haloes further down to $z=0$ from the last simulation snapshot ($z=z_f$), based on the halo growth model in Fakhouri, Ma & Boylan-Kolchin (2010) and the stellar-halo mass relation (SHMR) in Behroozi et al. (2019).

We assume that all haloes, except for significantly stripped (satellite) ones, remain quasi-isolated from $z_f$ and grow smoothly at the cosmic average rate as a function of halo mass $M_h$ and redshift $z$, inferred from the millennium cosmological simulations (Boylan-Kolchin et al. 2009) for $\Lambda$CDM (Fakhouri et al. 2010):

$$ \frac{dM_h}{dt} = 46.1 M_{\odot} \text{ yr}^{-1} \left( \frac{M_h}{10^{12} M_{\odot}} \right)^{1.1} \times (1 + 1.11z) \sqrt{\Omega_m (1+z)^3 + \Omega_{\Lambda}} \quad (21) $$

Once the halo mass $M_h(z)$ is known, we can derive the galaxy (stellar) mass with

$$ M_{*}(z) = (B + 1) M_{*\text{,SHMR}}(M_{h}(z), z), \quad (22) $$

where $M_{*\text{,SHMR}}(M_h, z)$ is given by the SHMR fitting formula for the true stellar mass values including both star-forming and quiescent galaxies and excluding intrahalo light in Behroozi et al. (2019), see column 8 of their table J1), and $B$ is a normalization factor. Here, we extrapolate the fitting formula to low-mass haloes with a minimum star formation efficiency of $\eta \equiv M_* \Omega_m / (\Omega_b M_h) = 10^{-3}$. In this way, we ignore the fluctuations of the star formation rate with respect to the cosmic average caused by different halo assembly histories and baryon cycles<fnref n="24" /> under the same halo mass. Such stochastic effects are captured by A-SLOTH and can cause an offset between $M_{*\text{,SHMR}}$ and $M_{*}(z_f)$ predicted by A-SLOTH at $z_f$. Assuming that this offset decays after cosmic noon at $z_{\text{noon}} = 2$, for each halo we evolve $B$ from the initial condition $[(B+1)M_{*\text{,SHMR}}]|_{z=z_f} = M_{*}(z_f)$ with

$$ \frac{dB}{dt} = \begin{cases} 0, & z > z_{\text{noon}}, \\ -\xi \frac{B}{t_H}, & z \le z_{\text{noon}}, \end{cases} \quad (23) $$

<figure><img src="image_5.png" /><figcaption>**Figure 3.** Co-moving SFRDs of Pop I/II (solid) and Pop III stars (dashed), predicted by the fiducial Pop III IMF ($\alpha = 1$) and the best-fitting star formation and stellar feedback parameters in A-SLOTH (see their table 3 Hartwig et al. 2023) coupled with our extrapolation scheme, applied on the merger trees from the cosmological simulation by Ishiyama et al. (2016). 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 a scatter of 0.2 dex embodied by the shaded region). For comparison, we also show the Pop I/II and Pop III SFRDs from Liu & Bromm (2020a) with the dashed-dotted and dotted curves, respectively. The thin vertical line denotes the final redshift $z_f \simeq 4.5$ above which merger trees are constructed from the cosmological simulation.</figcaption></figure>

where $t_H = 1/H(t)$ given the Hubble parameter at $H(t)$, and $\xi = 2$ is a parameter that governs the decay rate. A satellite halo is regarded as significantly stripped when $B+1 > B_{\text{max}}+1 = 25$, whose halo/galaxy mass remains constant.<fnref n="25" /> We also consider that in the post-reionization epoch ($z<6$), haloes with $M_h < 6.7 \times 10^8 M_{\odot}$ $[(1+z)/5]^{-3/2}$ cannot form new stars according to the reionization models in Pawlik, Schaye & Dalla Vecchia (2015), Pawlik et al. (2017), Benitez-Llambay & Frenk (2020), and Hutter et al. (2021). These haloes still grow with equation (21) but their stellar masses are no longer updated with equation (22). Besides, we impose an upper bound to galaxy mass (growth) $M_* < \eta_{\text{max}}(\Omega_b/\Omega_m)M_{\text{peak}}$, given $M_{\text{peak}}$ the maximum mass the halo has ever reached and $\eta_{\text{max}} = 0.5$ the maximum star formation efficiency.

We adopt fixed (star formation) time-steps of $\delta t_{\text{post}} = 1$ Myr for the low-$z$ evolution of galaxies/haloes with the above formalism (equations 21–23) where the values of $z_{\text{noon}}$, $\xi$, $B_{\text{max}}$, and $\eta_{\text{max}}$ are chosen to reproduce the cosmic star formation rate density (SFRD) at $z \lesssim 6$ in observations (e.g. Madau & Dickinson 2014; Finkelstein 2016). For instance, Fig. 3 shows the star formation histories of Pop III and Pop I/II stars predicted by the default star formation and stellar feedback parameters in A-SLOTH (see their table 3 Hartwig et al. 2022) coupled with our extrapolation scheme (equations 21 and 22), based on the merger trees from the cosmological simulation in Ishiyama et al. (2016), for the fiducial Pop III IMF ($\alpha = 1$, see Section 2.1). The Pop I/II SFRD, which dominates the total SFRD at $z \lesssim 15$ in our model, is slightly overestimated at $z \to 0$ compared with observations. The reason is that we did not include explicitly

<footnote n="24">²⁴There is evidence from recent observations that fluctuations of star formation rates are weak for intermediate-mass galaxies ($M_* \sim 10^9-10^{10} M_{\odot}$) at $z \sim 0.3-0.4$ (Patel et al. 2023).</footnote>

<footnote n="25">²⁵The parameter $B+1$ measures the mass fraction of stars in the halo, and $B \sim 25$ means that stars make up the entire halo mass under a cosmic average star formation efficiency $\eta \equiv M_{*\text{,SHMR}} \Omega_m / (\Omega_b M_h) \sim 0.2$, or at least the entire baryonic mass for $\eta \gtrsim 0.03$. Here, $\eta \sim 0.03-0.2$ is typical for galaxies with $M_* \gtrsim 10^9 M_{\odot}$ at $z \lesssim 4$ (Behroozi et al. 2019).</footnote>
<table><caption>**Table 2.** Summary of models. Column 1 gives the model name. Columns 2, 3, 4, and 5 define the initial condition model underlying the input SEVN catalogue, in terms of the initial distributions of primary stellar mass, mass ratio, orbital period, and eccentricity, respectively (see Section 2.1). Here, the primary star IMF is specified with the power-law slope at the high-mass end (given the fixed ZAMS mass range $m_{*,1} \in [5, 550]$ M$_{\odot}$), which is identical to that of the Pop III IMF adopted in A-SLOTH to model stellar feedback. Columns 6 and 7 define the NSC (and DF) model with the NSC occupation fraction as a function of galaxy (stellar) mass and the minimum galaxy mass for DF and NSC formation at $z > 5$. See Section 3.1 for a detailed description of how these parameters are used to determine whether Pop III BBHs can fall into the NSC by DF in a galaxy. Column 8 shows the BBH formation efficiency $\epsilon_{BBH}$, i.e. the number of BBHs formed per unit stellar mass. Similarly, Column 9 shows the BBH merger efficiency $\epsilon_{GW}^{\text{all}}$ (field/NSC), i.e. the number of BBH mergers at $z > 0$ per unit stellar mass, for all (field/NSC) mergers. Column 10 shows the fraction of mergers in NSCs $f_{\text{NSC}} = \epsilon_{\text{GW}}^{\text{NSC}}/\epsilon_{\text{GW}}^{\text{all}}$. Column 11 shows the fraction $f_{\text{GW}}^{\text{field/NSC}}$ of BBHs that merge at $z > 0$ in galaxy fields/NSCs. The last column shows the fraction $f_{\text{infall}}$ of Pop III BBHs that fall into NSCs.</caption><thead><tr><th>Model</th><th>α</th><th>p(q)</th><th>p(π)</th><th>p(e)</th><th>f<sub>occ</sub></th><th>M<sub>*</sub>,min<br>[M<sub>☉</sub>]</th><th>ε<sub>BBH</sub><br>[10<sup>-4</sup> M<sub>☉</sub><sup>-1</sup>]</th><th>ε<sub>GW</sub><sup>all</sup> (field/NSC)<br>[10<sup>-5</sup> M<sub>☉</sub><sup>-1</sup>]</th><th>f<sub>NSC</sub></th><th>f<sub>GW</sub><sup>field/NSC</sup></th><th>f<sub>infall</sub></th></tr></thead><tbody><tr><td>LOG1obs</td><td>1</td><td>S12</td><td>S12</td><td>S12</td><td>f<sub>occ</sub></td><td>10<sup>6</sup></td><td>4.78</td><td>6.49 (5.33/1.16)</td><td>17.9%</td><td>11.6%/62.8%</td><td>3.87%</td></tr><tr><td>TOP1obs</td><td>0.17</td><td>S12</td><td>S12</td><td>S12</td><td>f<sub>occ</sub></td><td>10<sup>6</sup></td><td>2.58</td><td>2.38 (1.66/0.722)</td><td>30.3%</td><td>6.73%/64.1%</td><td>4.36%</td></tr><tr><td>KRO1obs</td><td>2.3</td><td>S12</td><td>S12</td><td>S12</td><td>f<sub>occ</sub></td><td>10<sup>6</sup></td><td>7.05</td><td>11.8 (10.2/1.57)</td><td>13.3%</td><td>15.1%/60.6%</td><td>3.67%</td></tr><tr><td>LOG5obs</td><td>1</td><td>SB13</td><td>SB13</td><td>2e</td><td>f<sub>occ</sub></td><td>10<sup>6</sup></td><td>10</td><td>3.64 (0.685/2.95)</td><td>81.2%</td><td>0.716%/61%</td><td>4.82%</td></tr><tr><td>TOP5obs</td><td>0.17</td><td>SB13</td><td>SB13</td><td>2e</td><td>f<sub>occ</sub></td><td>10<sup>6</sup></td><td>7.94</td><td>2.76 (0.283/2.48)</td><td>89.8%</td><td>0.376%/62.5%</td><td>4.99%</td></tr><tr><td>KRO5obs</td><td>2.3</td><td>SB13</td><td>SB13</td><td>2e</td><td>f<sub>occ</sub></td><td>10<sup>6</sup></td><td>8.9</td><td>3.14 (0.755/2.38)</td><td>75.9%</td><td>0.888%/58.6%</td><td>4.56%</td></tr><tr><td>LOG1full</td><td>1</td><td>S12</td><td>S12</td><td>S12</td><td>1</td><td>10<sup>6</sup></td><td>4.78</td><td>7.98 (5.27/2.71)</td><td>34%</td><td>12.2%/58.9%</td><td>9.65%</td></tr><tr><td>TOP1full</td><td>0.17</td><td>S12</td><td>S12</td><td>S12</td><td>1</td><td>10<sup>6</sup></td><td>2.58</td><td>3.33 (1.64/1.68)</td><td>50.6%</td><td>7.14%/60.3%</td><td>10.8%</td></tr><tr><td>KRO1full</td><td>2.3</td><td>S12</td><td>S12</td><td>S12</td><td>1</td><td>10<sup>6</sup></td><td>7.05</td><td>13.7 (10.1/3.64)</td><td>26.5%</td><td>15.7%/56.4%</td><td>9.14%</td></tr><tr><td>LOG5full</td><td>1</td><td>SB13</td><td>SB13</td><td>2e</td><td>1</td><td>10<sup>6</sup></td><td>10</td><td>7.29 (0.674/6.62)</td><td>90.8%</td><td>0.761%/55.6%</td><td>11.9%</td></tr><tr><td>TOP5full</td><td>0.17</td><td>SB13</td><td>SB13</td><td>2e</td><td>1</td><td>10<sup>6</sup></td><td>7.94</td><td>5.8 (0.276/5.53)</td><td>95.2%</td><td>0.396%/56.7%</td><td>12.3%</td></tr><tr><td>KRO5full</td><td>2.3</td><td>SB13</td><td>SB13</td><td>2e</td><td>1</td><td>10<sup>6</sup></td><td>8.9</td><td>6.01 (0.741/5.27)</td><td>87.7%</td><td>0.937%/52.9%</td><td>11.2%</td></tr><tr><td>LOG1low</td><td>1</td><td>S12</td><td>S12</td><td>S12</td><td>f<sub>occ</sub></td><td>10<sup>10</sup></td><td>4.78</td><td>6.07 (5.34/0.73)</td><td>12%</td><td>11.5%/50.2%</td><td>3.04%</td></tr><tr><td>TOP1low</td><td>0.17</td><td>S12</td><td>S12</td><td>S12</td><td>f<sub>occ</sub></td><td>10<sup>10</sup></td><td>2.58</td><td>2.16 (1.66/0.495)</td><td>22.9%</td><td>6.69%/53.8%</td><td>3.57%</td></tr><tr><td>KRO1low</td><td>2.3</td><td>S12</td><td>S12</td><td>S12</td><td>f<sub>occ</sub></td><td>10<sup>10</sup></td><td>7.05</td><td>11.2 (10.3/0.894)</td><td>8%</td><td>15%/45.8%</td><td>2.77%</td></tr><tr><td>LOG5low</td><td>1</td><td>SB13</td><td>SB13</td><td>2e</td><td>f<sub>occ</sub></td><td>10<sup>10</sup></td><td>10</td><td>2.71 (0.686/2.02)</td><td>74.7%</td><td>0.711%/51.9%</td><td>3.89%</td></tr><tr><td>TOP5low</td><td>0.17</td><td>SB13</td><td>SB13</td><td>2e</td><td>f<sub>occ</sub></td><td>10<sup>10</sup></td><td>7.94</td><td>2.09 (0.28/1.81)</td><td>86.6%</td><td>0.369%/55.2%</td><td>4.13%</td></tr><tr><td>KRO5low</td><td>2.3</td><td>SB13</td><td>SB13</td><td>2e</td><td>f<sub>occ</sub></td><td>10<sup>10</sup></td><td>8.9</td><td>2.18 (0.755/1.42)</td><td>65.4%</td><td>0.879%/45.4%</td><td>3.52%</td></tr></tbody></table>

galaxy mergers and mechanisms of quenching (e.g. AGN feedback and environmental quenching) in the low-$z$ regime. Nevertheless, by experimenting with other star formation models, we find that varying the late-time ($z \lesssim 1$) SFRD by up to a factor of $\sim 10$ has little impact on Pop III BBH mergers.

Given the star formation histories $M_{\star}(t)$ of individual haloes, all galaxy and NSC properties can be derived, which set the background for the evolution of SCs and Pop III binaries as explained in the previous sections and Liu & Bromm (2021). In reality, haloes are not isolated and will merge into larger haloes. This process disrupts the inspiral of Pop III BBHs in satellite galaxies but may enhance DF during galaxy mergers (Román et al. 2023). Galaxy mergers also strip satellite galaxies, so their NSCs can no longer grow and become vulnerable to tidal forces from the central galaxy, reducing their ability to facilitate Pop III mergers. Our results for the NSC infall rate and merger rate in NSCs at $z < z_f$ are likely rather optimistic.

To better characterize the host galaxies of Pop III mergers at $z < z_f$, we also estimate the average gas-phase metallicity $Z$ of each galaxy by extrapolating from the value at $z = z_f$ predicted by A-SLOTH with the redshift-dependent mass–metallicity relation $Z \propto M_*^{0.3}(1+z)^{-0.9}$ based on Langeroodi et al. (2023).

## 5 RESULTS

We apply A-SLOTH to the merger trees constructed from the cosmological simulation by Ishiyama et al. (2016) in a co-moving volume of $V_{\text{com}} = (8 h^{-1} \text{Mpc})^3 \simeq 1650$ Mpc$^3$ with a dark matter mass resolution of $5000 h^{-1} \text{M}_{\odot}$, which only cover $z \ge z_f \simeq 4.5$

where the box size is large enough to be marginally representative. We use the extrapolation scheme described in Section 4 to follow the subsequent evolution down to $z = 0$. Throughout this work, we adopt the cosmological parameters from Planck Collaboration XIII (2016)$^{26}$ and the default A-SLOTH star formation and feedback parameters (see their table 3 Hartwig et al. 2022). Following Liu et al. (2021a), Liu & Bromm (2021), and Hartwig et al. (2023), the mass fraction of Pop III stars in binaries is fixed as $f_B = 0.7$. Since various stochastic processes are involved in our modelling and the simulation volume underlying the merger trees is not very large,$^{27}$ we boost the number of Pop III binaries to be sampled in A-SLOTH by a factor of $f_{\text{boost}} = 50$ to achieve better statistics and reduce numerical noise in the merger history, so the effective simulation volume $\tilde{V}_{\text{com}} = 50V_{\text{com}}$ is used to calculate the cosmic average quantities from our merger populations, assuming that the actual simulation volume in Ishiyama et al. (2016) is cosmologically representative.

We consider 18 models combining different choices of the Pop III IMF, IBS, and NSC parameters (defined in Sections 2.1 and 3.1), as summarized in Table 2 (see also Tables 3 and 4 for the key properties of Pop III mergers in these models). We divide Pop III mergers into two groups based on their sites: galaxy fields and NSCs. Here and henceforth, we use the term ‘NSC’ to denote all SCs including both NSCs and NSC descendants (i.e. normal

<footnote>$^{26}$ $\Omega_m = 0.3089$, $\Omega_b = 0.0486$, and $H_0 = 100h \text{ km s}^{-1} \text{ Mpc}^{-1}$ with $h = 0.6774$.</footnote>

<footnote>$^{27}$ The total mass of Pop III stars formed in the original simulation volume is $\sim 7.6–9.4 \times 10^7 \text{ M}_{\odot}$, which only corresponds to $\sim 2.4–8.4 \times 10^4$ BBHs.</footnote>
<table><caption>**Table 3.** Basic observational signatures of Pop III BBH mergers for the 18 models listed in Table 2. Column 1: model name. Columns 2: all-sky merger rate $\mathcal{N}$ (equation 24) for all mergers within $z < 1$. Column 3: fraction $\mathcal{F}_{NSC}$ of all-sky merger rate contributed by NSC mergers for $z < 1$. Column 4: fraction $\mathcal{F}_{PopIII}$ of all-sky merger rate of BBH mergers contributed by Pop III stars for $z < 1$. Column 5: all-sky merger rate $\mathcal{N}_{IMBH}$ of Pop III BBHs with at least one IMBH (i.e. $m_1 > 100 M_\odot$) for $z < 1$. Columns 6–9: same as Columns 2–5 but for $z < 20$. Columns 10 and 11: peak value $\Omega_{\text{GW}}^{\text{peak}}$ and location $v_{\text{peak}}$ of the SGWB energy spectrum from Pop III BBH mergers. Column 12: maximum ratio $f_{\text{SGWB}}^{\text{PopIII}}$ between the SGWB produced by Pop III BBH mergers and the conservative prediction of the SGWB from Pop I/II and PBH mergers by Bavera et al. (2022). Given the optimistic results for the total SGWB inferred from GWTC-3 events (Abbott et al. 2023a) as the reference, the maximum contribution of Pop III BBH mergers to the total SGWB is approximately $f_{\text{SGWB}}^{\text{PopIII}}/5$.</caption><thead><tr><th>Model</th><th>N [yr<sup>-1</sup>]<br>(z &lt; 1)</th><th>F<sub>NSC</sub><br>(z &lt; 1)</th><th>F<sub>PopIII</sub><br>(z &lt; 1)</th><th>N<sub>IMBH</sub><br>(z &lt; 1)</th><th>N [yr<sup>-1</sup>]<br>(z &lt; 20)</th><th>F<sub>NSC</sub><br>(z &lt; 20)</th><th>F<sub>PopIII</sub><br>(z &lt; 20)</th><th>N<sub>IMBH</sub><br>(z &lt; 20)</th><th>&Omega;<sub>GW</sub><sup>peak</sup><br>[10<sup>-12</sup>]</th><th>v<sub>peak</sub><br>[Hz]</th><th>f<sub>PopIII</sub><br><sub>SGWB</sub></th></tr></thead><tbody><tr><td>LOG1obs</td><td>8.95</td><td>41.6%</td><td>0.423%</td><td>1.34</td><td>650</td><td>4.75%</td><td>2.93%</td><td>11.6</td><td>10.4</td><td>14</td><td>13.5%</td></tr><tr><td>TOP1obs</td><td>4.87</td><td>63%</td><td>0.231%</td><td>2.18</td><td>264</td><td>8.41%</td><td>1.21%</td><td>16.2</td><td>13.8</td><td>12.6</td><td>18.7%</td></tr><tr><td>KRO1obs</td><td>21.2</td><td>30.9%</td><td>0.998%</td><td>0.304</td><td>1400</td><td>3.35%</td><td>6.12%</td><td>2.56</td><td>11.8</td><td>105</td><td>9.05%</td></tr><tr><td>LOG5obs</td><td>10.1</td><td>91.9%</td><td>0.476%</td><td>5.35</td><td>151</td><td>54%</td><td>0.697%</td><td>48.3</td><td>28.7</td><td>13.1</td><td>40.6%</td></tr><tr><td>TOP5obs</td><td>10.6</td><td>96%</td><td>0.499%</td><td>8.53</td><td>117</td><td>70.4%</td><td>0.539%</td><td>70.7</td><td>44.8</td><td>12.6</td><td>64.8%</td></tr><tr><td>KRO5obs</td><td>10.5</td><td>88.9%</td><td>0.494%</td><td>1.19</td><td>174</td><td>44.3%</td><td>0.799%</td><td>11.2</td><td>10.8</td><td>92.8</td><td>11.7%</td></tr><tr><td>LOG1full</td><td>13.4</td><td>62.4%</td><td>0.631%</td><td>2.8</td><td>707</td><td>12.7%</td><td>3.17%</td><td>32.5</td><td>17.6</td><td>13.1</td><td>26.4%</td></tr><tr><td>TOP1full</td><td>8.1</td><td>78.7%</td><td>0.383%</td><td>4.33</td><td>309</td><td>21.9%</td><td>1.41%</td><td>48.2</td><td>24.6</td><td>11.6</td><td>39.2%</td></tr><tr><td>KRO1full</td><td>28.2</td><td>50.2%</td><td>1.32%</td><td>0.615</td><td>1490</td><td>9.56%</td><td>6.48%</td><td>6.91</td><td>16.7</td><td>98.6</td><td>13%</td></tr><tr><td>LOG5full</td><td>20.7</td><td>96.2%</td><td>0.973%</td><td>10.9</td><td>310</td><td>77.7%</td><td>1.42%</td><td>140</td><td>50.4</td><td>11.6</td><td>83.6%</td></tr><tr><td>TOP5full</td><td>20.7</td><td>98.2%</td><td>0.974%</td><td>16.7</td><td>279</td><td>87.7%</td><td>1.28%</td><td>207</td><td>78.9</td><td>11.2</td><td>134%</td></tr><tr><td>KRO5full</td><td>21</td><td>94.7%</td><td>0.988%</td><td>2.38</td><td>320</td><td>69.9%</td><td>1.46%</td><td>30</td><td>19.3</td><td>83.9</td><td>23.3%</td></tr><tr><td>LOG1low</td><td>8</td><td>33.9%</td><td>0.379%</td><td>1.08</td><td>626</td><td>1.05%</td><td>2.82%</td><td>2.76</td><td>11</td><td>14.8</td><td>13.3%</td></tr><tr><td>TOP1low</td><td>4.02</td><td>54.9%</td><td>0.191%</td><td>1.61</td><td>248</td><td>2.3%</td><td>1.14%</td><td>4.32</td><td>14.7</td><td>14.5</td><td>17.9%</td></tr><tr><td>KRO1low</td><td>19</td><td>21.6%</td><td>0.893%</td><td>0.235</td><td>1370</td><td>0.66%</td><td>5.97%</td><td>0.56</td><td>10.7</td><td>116</td><td>8.37%</td></tr><tr><td>LOG5low</td><td>8</td><td>89.7%</td><td>0.378%</td><td>4.47</td><td>87.4</td><td>20.2%</td><td>0.404%</td><td>11.7</td><td>33.2</td><td>14.5</td><td>40.3%</td></tr><tr><td>TOP5low</td><td>8.22</td><td>95.1%</td><td>0.389%</td><td>6.81</td><td>54.8</td><td>36.9%</td><td>0.253%</td><td>17.9</td><td>52.3</td><td>14.2</td><td>64.1%</td></tr><tr><td>KRO5low</td><td>7.26</td><td>83.9%</td><td>0.343%</td><td>0.945</td><td>110</td><td>12.2%</td><td>0.508%</td><td>2.44</td><td>9.47</td><td>101</td><td>10.5%</td></tr></tbody></table>

<table><caption>**Table 4.** Summary statistics of BH masses in Pop III BBH mergers for the 18 models listed in Table 2. The second, third, and last columns show the average chirp mass $m_c$, primary mass $m_1$, and mass ratio $q$ for BBH mergers in galaxy fields/NSCs.</caption><thead><tr><th>Model</th><th>m&#773;<sub>c</sub><sup>field/NSC</sup><br>[M<sub>☉</sub>]</th><th>m&#773;<sub>1</sub><sup>field/NSC</sup><br>[M<sub>☉</sub>]</th><th>q&#773;<sup>field/NSC</sup></th></tr></thead><tbody><tr><td>LOG1obs</td><td>23.4/97.2</td><td>31.3/160</td><td>0.75/0.64</td></tr><tr><td>TOP1obs</td><td>25.3/169</td><td>33.8/282</td><td>0.76/0.58</td></tr><tr><td>KRO1obs</td><td>20.4/38.4</td><td>28.1/56.7</td><td>0.71/0.71</td></tr><tr><td>LOG5obs</td><td>25.8/131</td><td>37/242</td><td>0.69/0.55</td></tr><tr><td>TOP5obs</td><td>27.9/185</td><td>41.7/345</td><td>0.68/0.48</td></tr><tr><td>KRO5obs</td><td>24.4/58.3</td><td>34.6/94.3</td><td>0.7/0.68</td></tr><tr><td>LOG1full</td><td>23.4/94.2</td><td>31.3/154</td><td>0.75/0.64</td></tr><tr><td>TOP1full</td><td>25.2/165</td><td>33.5/276</td><td>0.77/0.58</td></tr><tr><td>KRO1full</td><td>20.4/37.5</td><td>28/55.4</td><td>0.71/0.7</td></tr><tr><td>LOG5full</td><td>25.6/129</td><td>36.5/238</td><td>0.69/0.55</td></tr><tr><td>TOP5full</td><td>26.2/182</td><td>38.1/341</td><td>0.68/0.48</td></tr><tr><td>KRO5full</td><td>24.3/57.3</td><td>34.4/92.2</td><td>0.7/0.69</td></tr><tr><td>LOG1low</td><td>23.4/107</td><td>31.3/174</td><td>0.75/0.64</td></tr><tr><td>TOP1low</td><td>25.2/180</td><td>33.6/295</td><td>0.76/0.58</td></tr><tr><td>KRO1low</td><td>20.4/40.5</td><td>28.1/59.9</td><td>0.71/0.71</td></tr><tr><td>LOG5low</td><td>25.6/143</td><td>36.6/261</td><td>0.69/0.54</td></tr><tr><td>TOP5low</td><td>26.3/193</td><td>38.3/357</td><td>0.68/0.48</td></tr><tr><td>KRO5low</td><td>24.2/63.4</td><td>34.4/104</td><td>0.7/0.68</td></tr></tbody></table>

SCs) defined in Section 3 if not specially clarified. Most mergers in galaxy fields come from BBHs that evolve in isolation throughout their lifetimes, representative for the IBSE channel, except for a tiny (≤ 2 per cent) fraction of mergers from BBHs returned to galaxy fields from NSCs (via ejection by binary-single encounters

or SC evaporation/disruption), which are influenced by the NSC-DH channel. We do not exclude these special field mergers when discussing the features of the IBSE channel. In this section, we first present the results for the fiducial model and discuss the general properties of Pop III BBH mergers (Section 5.1), in comparison with those of mergers of Pop I/II BHs and primordial back holes (PBHs) predicted by Franciolini et al. (2022b) and Bavera et al. (2022) based on low-z observations (Abbott et al. 2021, 2023a). Next, we explore the dependence of our results on the parameters of Pop III IMF, IBS, and NSCs (Section 5.2). Finally, in Section 5.3, we briefly discuss the observational perspective of Pop III BBH mergers, focusing on the detection rates of Pop III BBH mergers by ET<fnref n="28" /> and the LVK network during O4²⁹ using the python package GWTOOLBOX<fnref n="30" /> (Yi et al. 2022a, b), and current observations of BBH mergers involving massive BHs (like GW190521).

Following Bavera et al. (2022), we adopt the best-fitting model in Franciolini et al. (2022b) as the reference to evaluate the Pop III contributions to the GW signals of BBH mergers. This model combines PBH mergers with the Pop I/II BBH mergers formed in the IBSE channel via common envelope (CE) and stable mass transfer (SMT), as well as dynamically in GCs, with branching ratios inferred from the events in the second gravitational wave transient catalogue (GWTC-2; Abbott et al. 2021). According to

<footnote n="28">²⁸We use the ET-D-sum sensitivity model from https://www.et-gw.eu/index.php/etsensitivities (also see Hild, Chelkowski & Freise 2008; Hild et al. 2010, 2011; Branchesi et al. 2023).</footnote>

<footnote n="29">²⁹We use the sensitivity model for advanced LIGO from https://dcc.ligo.org/LIGO-T200012-v1/public (aligo_O4high.txt) to represent the maximum capacity of the LVK network during O4.</footnote>

<footnote n="30">³⁰https://bitbucket.org/radboudradiolab/gwtoolbox/src/master/</footnote>
this model, most (~94 per cent) events in GWTC-2 are explained by Pop I/II BBH mergers, and the remaining ones (especially the special event GW190521) are attributed to PBHs that follow a log-normal distribution with a characteristic mass of ~35 M⊙ and a width of σ = 0.41, and make up $f_{PBH} \sim 2 \times 10^{-4}$ of dark matter (see their fig. 4 Franciolini et al. 2022b). Here, the contribution of PBHs to the observed events is statistically insignificant and subject to uncertainties in the modelling of astrophysical populations. Considering every possible origin of BBH mergers for comparison is beyond the scope of this work, and we do not expect the reference model to be the most representative model for the populations of BBH mergers other than that from Pop III stars. In fact, large uncertainties remain in current theoretical predictions on the properties of BBH mergers, so inferences about the origins of BBH mergers and their branching ratios from the limited sample of detected events are highly model-dependent and tentative. The reference model considered here produces relatively conservative predictions for the merger rate and SGWB of BBH mergers involving Pop I/II BHs and PBHs compared with other models in the literature. Therefore, our results on the relative contributions of Pop III mergers should be regarded as heuristic optimistic estimates (given the specific Pop III star formation histories predicted by A-SLOTH) meant to illustrate the special features of Pop III BBH mergers, and the exact values can be off by a factor of a few. For simplicity, we focus on the intrinsic properties of entire BBH merger populations at $z > 0$ when comparing our results with those in the reference model (Bavera et al. 2022; Franciolini et al. 2022b). We plan to investigate the detailed (redshift-dependent) observational signatures of BBH mergers of different origins and evolution channels in future work fully taking into account source detectability and characterization (see, e.g. Santoliquido et al. 2023; Tanikawa et al. 2022b).

## 5.1 Fiducial model

The fiducial model LOG1-obs is defined with a log-flat IMF for primary stars ($\alpha = 1$), the IBS from local observations (S12), and the empirical NSC occupation fraction $f_{occ} = \hat{f}_{occ}$ (see equation 15) with $M_{*,min} = 10^6$ M⊙. The basic statistics of Pop III BBHs and their mergers in the fiducial model are summarized as follows: The mass fraction of binary stars that become BBHs is $f_{BBH} \simeq 12.6$ per cent, which (given $f_B = 0.7$ and $\tilde{m}_{p,BBH} \simeq 193$ M⊙) corresponding to a BBH formation efficiency (number of BBHs formed per unit stellar mass) $\epsilon_{BBH} \simeq 4.78 \times 10^{-4} M_{\odot}^{-1}$. 3.55 per cent of Pop III BBHs fall into NSCs, among which 21.3 per cent are disrupted as soft binaries, and 37 per cent have reached the galaxy centre by DF before the formation of the NSC. In galaxy fields (NSCs), 11.6 per cent (62.8 per cent) of Pop III BBHs have merged at $z > 0$. The overall merger efficiency (number of mergers at $z > 0$ per unit stellar mass) is $\epsilon_{GW}^{all} = 6.49 \times 10^{-5} M_{\odot}^{-1}$. 17.9 per cent of Pop III mergers occur in NSCs across cosmic history, among which only ~2 per cent happen in NSC descendants (i.e. normal SCs). The small fraction of NSC mergers in NSC descendants is an artificial outcome of our low-z extrapolation scheme that ignores halo mergers. Since mergers between galaxies hosting NSCs are rare at high-z, NSC descendants (produced by galaxy mergers) are expected to be important after cosmic noon ($z \lesssim 2$), which is not covered by the merger trees from Ishiyama et al. (2016). In fact, for merger trees targeting Milky Way-mass galaxies (reaching $z = 0$) from the Caterpillar project (Griffen et al. 2016), we find that NSC descendants can host up to ~60 per cent of Pop III BBH mergers in NSCs.

<figure><img src="image_6.png" /><figcaption>**Figure 4.** Co-moving MRD of Pop III BBHs as a function of redshift in the fiducial model LOG1-obs. The results for mergers in galaxy fields, NSCs, and combined are shown with the thin solid, dotted, and thick solid curves, respectively. The rate density of NSC infall is also shown with the dashed-dotted curve. For comparison, we plot the MRDs of all Pop I/II BBHs (long-dashed) in both galaxy fields (formed via CE and SMT) and GCs and only those in GCs (dashed) as well as that of PBH binaries (dashed-dotted-dotted), from Bavera et al. (2022, see their fig. 3), normalized to the model selection results in Franciolini et al. (2022b). The thick dotted curve shows the total MRD combining Pop I/II and PBH mergers, while the thick dashed curve shows the MRD of all stellar BBHs (Pop I/II + III). We find that Pop III mergers account for~0.3 per cent of the (redshift-dependent) local BBH MRD $\dot{n}_{obs}(z=0) = 19.3^{+15.1}_{-9} \text{ yr}^{-1} \text{ Gpc}^{-3}$ inferred from the GWTC-2 events (triangle data point, Abbott et al. 2021).</figcaption></figure>

### 5.1.1 Merger rate

Fig. 4 shows the (co-moving) merger rate density (MRD) $\dot{n}$ as a function of redshift for Pop III BBHs in comparison with the merger histories of Pop I/II BBHs and PBHs from Franciolini et al. (2022b) and the (redshift-dependent) local BBH MRD $\dot{n}_{obs}(z=0) = 19.3^{+15.1}_{-9} \text{ yr}^{-1} \text{ Gpc}^{-3}$ inferred from the GWTC-2 events (Abbott et al. 2021). Here, the MRD at a redshift bin $z_j$ is derived by counting the number of mergers $N_j$ in the redshift range $z \in [z_j - 0.5\Delta z_j, z_j + 0.5\Delta z_j)$, i.e. $\dot{n}(z_j) = N_j/(\bar{V}_{\text{com}}\Delta t_j)$, where $\Delta t_j$ is the cosmic age evolution across the redshift bin. Sixty evenly spaced bins are used for $z \sim 0-30$ with $\Delta z_j = 0.5$.

We find that the Pop III MRD is dominated by field mergers at $z \gtrsim 1$ and peaks at $z \sim 15$ with $\dot{n}_{\text{peak}} \sim 3.6 \text{ yr}^{-1} \text{ Gpc}^{-3}$, which is close to the peak of Pop III star formation (Fig. 3). The MRD of field mergers is consistent with that derived in Santoliquido et al. (2023) for the IBSE channel using the same Pop III star formation history predicted by the fiducial A-SLOTH model. The consistency also holds for the other 17 models considered in Table 2. The MRD from Pop III mergers in NSCs increases towards lower redshift and reaches a plateau at $z \lesssim 6$, and it slightly exceeds the MRD of field mergers at $z = 0$. It is also shown in Fig. 4 that NSC infall only starts at $z \sim 11$, because massive enough galaxies that host NSCs ($M_{*} > 10^6$ M⊙) only form at $z \lesssim 12$, and the DF time-scale for Pop III BBHs can be large. The different trends in the MRDs of field and NSC mergers can be understood through the delay time ($t_{\text{delay}}$) distribution, shown in Fig. 5. Here, the delay time is defined as the time taken from the initial formation of the stellar binary to the final BBH merger (or NSC infall), which is simply $t_{\text{delay}} = \tau_{\text{B}} + t_{\text{GW}}$ for BBHs that evolve in isolation (Section 2.2.2). The delay time distribution of field mergers is approximately log-flat down to $t_{\text{delay}} \sim 5$ Myr.
<figure><img src="image_7.png" /><figcaption>Figure 5. Delay time distribution for Pop III BBH mergers in galaxy fields (solid) and NSCs (dotted), as well as NSC infall (dash-dotted) and DH (long-dashed) of BBHs in the fiducial model LOG1obs.</figcaption></figure>

However, it takes at least ~200 Myr for Pop III BBHs to fall into NSCs in the first place. For BBH mergers in NSCs, we have $t_{delay} \gtrsim 300$ Myr, and the delay time distribution is almost uniform (i.e. $dN/d\log t_{delay} \propto t_{delay}$) at $t_{delay} \gtrsim 800$ Myr. In Fig. 5, we also show the distribution of the DH time-scale, i.e. the time taken from NSC infall to the final merger, which turns out to be nearly uniform below ~5 Gyr.

The predicted local MRD of Pop III BBH mergers is ~0.06 $yr^{-1}$ Gpc$^{-3}$, which only counts for $0.31_{-0.13}^{+0.27}$ per cent of the value inferred from observations (Abbott et al. 2021). The Pop III MRD is overwhelmed by those from Pop I/II and PBH mergers at low $z$ but its importance increases beyond cosmic noon from $z \sim 2$ to $z \sim 15$. The Pop III MRD exceeds that of Pop I/II mergers in GCs at $z \gtrsim 6$ and that of all Pop I/II mergers at $z \gtrsim 13$, although it remains below the PBH MRD by at least a factor of ~4, since the PBH rate increases monotonically towards higher $z$. This indicates that although we are more likely to detect Pop III mergers at higher $z$ for $z \sim 2-15$, distinguishing them from PBH mergers can be a challenge (see e.g. Franciolini et al. 2022a).

Given the (intrinsic) MRD $\dot{n}(z)$, the all-sky merger rate of Pop III BBHs observed at $z=0$ as a function of the horizon redshift $z_{horizon}$ can be written as (Hartwig et al. 2016)

$$ \mathcal{N}(z < z_{\text{horizon}}) = \int_{0}^{z_{\text{horizon}}} dz \frac{\dot{n}(z)}{(1+z)} \frac{dV}{dz}, \qquad (24) $$

where $dV/dz = 4\pi c d_c^2(z)/H(z)$ is the co-moving volume element, given the co-moving distance $d_c(z) = \int_0^z dz' [c/H(z')]$ and Hubble parameter $H(z)$. One can apply a weight to each source using its detection probability above a certain signal-to-noise ratio (SNR) to calculate the detection rate $\mathcal{N}_X$ by a given detector X with finite sensitivity using the same formalism. In this case, $\dot{n}(z)$ should be replaced with the MRD of detectable sources, i.e. $\dot{n}_X(z_j) = \sum_{i=1}^{N_j} p_{\text{det},i}$ for a redshift bin $z_j$ given the detection probabilities $p_{\text{det},i}$ of the $N_j$ mergers within the bin. In the calculation of the all-sky merger rate $\mathcal{N}$, we simply use $p_{\text{det},i} = 1$. Interestingly, this turns out to be a good approximation for our Pop III BBHs observed by ET. According to the results of $p_{\text{det},i}$ from GWTOOLBOX (Yi et al. 2022a, b), ET can detect most ($\gtrsim 90$ per cent) Pop III BBH mergers with SNR > 8. Similarly, most Pop I/II and PBH mergers in the reference model from Franciolini et al. (2022b) are detectable by ET with

<figure><img src="image_8.png" /><figcaption>Figure 6. All-sky merger rate (for an observer at $z=0$) of Pop III BBHs as a function of horizon redshift in the fiducial model LOG1obs. The notations of merger populations are the same as in Fig. 4. The vertical lines label the reference horizon redshift $z_{horizon} = 1$ and 20.</figcaption></figure>

SNR > 8 as well given the large horizon redshift of ET at their mass scales (see their fig. 8 De Luca et al. 2021). In both cases (Pop III and Pop I/II + PBH), the contribution of sources at very high redshifts is very small. Therefore, the Pop III contribution to the all-sky BBH merger rate within $z < 20$ is a good indicator of the chance of finding Pop III sources in observations of BBH mergers by ET.

Fig. 6 (see also Table 3) shows the results of $\mathcal{N}$ for our Pop III BBHs and the Pop I/II and PBH mergers from Franciolini et al. (2022b), where 300 bins are used for $z \sim 0-30$ with $\Delta z_j = 0.1$ to integrate equation (24) numerically. We predict an all-sky merger rate of Pop III BBHs $\mathcal{N} \sim 650 (9) yr^{-1}$ within $z < 20 (1)$, which counts for $\sim 3 (0.4)$ per cent of the total merger rate dominated by Pop I/II mergers. NSC mergers contribute $\sim 5$ and 42 per cent of the total Pop III BBH merger rates for $z_{horizon} = 20$ and 1. The contribution of NSC mergers saturates at $z_{horizon} \sim 7$. However, NSC mergers dominate ($\gtrsim 99.7$ per cent) the merger rate of massive BBHs containing at least one IMBH ($m_1 > 100 M_{\odot}$) $\mathcal{N}_{\text{IMBH}} = 11.6 (1.34) \text{yr}^{-1}$ within $z < 7 (1)$. The Pop III BBH all-sky merger rate remains below that of PBHs by a factor of $\sim 10$ for $z_{horizon} \lesssim 20$, and becomes comparable to the rate for Pop I/II BBHs in GCs for $z_{horizon} \gtrsim 20$. Here, the PBH merger rate keeps increasing with $z_{horizon}$, while the GC merger rate saturates around $z_{horizon} \sim 5$. Considering the facts that in our model universe, Pop III stars only make up $\sim 5 \times 10^{-5}$ of the total mass density of stars ever formed, and that Pop III stars are significantly outnumbered (by a factor of $\sim 1000$) by the PBHs (with a characteristic mass of $\sim 35 M_{\odot}$) that make up $\sim 2 \times 10^{-4}$ of dark matter (see their fig. 4 Franciolini et al. 2022b) to achieve merger rates $\sim 6-9$ times higher than those of Pop III BBHs within $z \lesssim 20$, the contribution by Pop III stars to the BBH merger rate is already quite large, which reflects the high efficiency of Pop III stars at forming BBH mergers.

### 5.1.2 Stochastic gravitational wave background

Next, we drive the contribution to SGWB by Pop III BBH mergers. Given $N$ mergers predicted by A-SLOTH in our effective simulation volume $\tilde{V}_{\text{com}} = 50V_{\text{com}}$, the SGWB (energy density spectrum) can be characterized by the dimensionless parameter (see e.g. Abbott et al. 2018; Braglia, García-Bellido & Kuroyanagi 2021; Périgois et al.
2021; Martinovic et al. 2022; Lehoucq et al. 2023):

$$ 
\begin{aligned}
\Omega_{\text{GW}}(\nu) &= \frac{\nu}{c^2 \rho_c} \int d\theta \int_0^\infty dz \frac{\dot{n}(z) p(\theta, z)}{(1+z)H(z)} \frac{dE_{\text{GW}}(\theta)}{dv_s} \\
&= \frac{\nu}{c^2 \rho_c} \lim_{\Delta z_j \to 0} \sum_{j=1}^{N_z} \sum_{i=1}^{N_j} \frac{(\Delta z_j / \Delta t_j) \tilde{V}_{\text{com}}^{-1}}{(1+z_j)H(z_j)} \frac{dE_{\text{GW},i}}{dv_s} \\
&= \frac{\nu}{c^2 \rho_c} \sum_{i=1}^{N} \frac{1}{\tilde{V}_{\text{com}}} \frac{dE_{\text{GW},i}}{dv_s},
\end{aligned}
\quad (25) $$

in which $\nu$ is the observer-frame frequency, $\rho_c = 3H_0^2/(8\pi G)$ is the critical density, $H(z)$ is the Hubble parameter at redshift $z$, $dE_{\text{GW}}/dv_s$ is the source-frame GW energy spectrum of an individual source evaluated at the source-frame frequency $v_s = (1+z)\nu$, and the factor $(1+z)$ in the denominator captures time dilation due to cosmic expansion. In the first line, $\dot{n}(z)$ denotes the source-frame co-moving merger rate density, and the source properties are captured by the parameter $\theta$ and its probability distribution $p(\theta, z)$. In the second line, we re-write the integration in its discrete form (with $\int dz \int d\theta \dot{n}(z) p(\theta, z) = \lim_{\Delta z_j \to 0} \sum_{j=1}^{N_z} \Delta z_j \sum_{i=1}^{N_j} \Delta t_j^{-1} \tilde{V}_{\text{com}}^{-1}$) as the summation over $N_z$ redshift bins $(z_j, \Delta z_j)$ each containing $N_j$ mergers at $z_i \in [z_j - 0.5\Delta z_j, z_j + 0.5\Delta z_j)$, whose GW spectra $dE_{\text{GW},i}/dv_s$ are added up, where $v_s = (1+z_i)\nu$, and $\Delta t_j$ is the corresponding cosmic age evolution at redshift bin $j$. Finally, since $dz/dt = (1+z)H(z)$, we arrive at the last line with $\lim_{\Delta z_j \to 0}(\Delta z_j / \Delta t_j) = (1+z_j)H(z_j)$. We use the phenomenological single-source energy spectrum model (for mergers with circular orbits) from Ajith et al. (2011) to calculate $dE_{\text{GW},i}/dv_s$, assuming random inclinations and a fixed effective spin$^{31}$ $\chi_{\text{eff}} = 0.06$ for all sources, which is around the peak of the $\chi_{\text{eff}}$ distribution, inferred from the BBH merger population detected by the LVK network (e.g. Callister & Farr 2024). The ignorance of eccentricities in our calculation is likely an oversimplification for some NSC mergers considering that large eccentricities can arise in dense star clusters and AGN discs (e.g. Hoang et al. 2018; Mapelli et al. 2021; Zhang et al. 2021; Arca Sedda et al. 2024a; Chattopadhyay et al. 2023; Dall’Amico et al. 2024; Fabj & Samsing 2024; Trani et al. 2024). We plan to take into account the effects of eccentricities on the GW signals in future work, which can reduce $\Omega_{\text{GW}}$ at low frequencies$^{32}$ (see e.g. Chen, Sesana & Del Pozzo 2017; Buskirk & Babiuc Hamilton 2023; Islam 2024; Raidal et al. 2024; Xuan et al. 2024). To the first order, given the intrinsic chirp mass $m_c = (m_1m_2)^{3/5}/(m_1+m_2)^{1/5}$ of a BBH merger at redshift $z$, its contribution to the SGWB $\delta\Omega_{\text{GW}} = (c^2\rho_c\tilde{V}_{\text{com}})^{-1}\nu(dE_{\text{GW}}/dv_s)$ peaks at an observer-frame frequency $v_p \propto 1/[(1+z)m_c]$ with $\delta\Omega_{\text{GW}}(v_p) \propto m_c/(1+z)$ and follows $\delta\Omega_{\text{GW}} \propto m_c^{5/3}v^{2/3}(1+z)^{-1/3}$ in the low-frequency inspiral regime ($v \ll v_p$). This indicates that the SGWB is mostly sensitive to massive mergers at low $z$.

Our results are shown in Fig. 7 (see also Table 3), compared with the SGWB energy spectra for Pop I/II and PBH mergers from Bavera et al. (2022). We find that although NSC mergers only account

<figure><img src="image_9.png" /><figcaption>Figure 7. SGWB (dimensionless energy density spectrum) contributed by Pop III BBH mergers in the fiducial model LOG1_0bs, derived with the phenomenological single-source energy spectrum model in Ajith et al. (2011) assuming a universal fixed effective spin $\chi_{\text{eff}} = 0.06$. Our Pop III results are compared and combined with those for Pop I/II and PBH mergers from Bavera et al. (2022, see their fig. 4) with the same notations of merger populations as in Fig. 4. The smaller and larger shaded regions denote the PI sensitivity curves of 2-year observations by the LVK network and ET for a signal-to-noise ratio SNR > 2, calculated by Bavera et al. (2022, see their appendix C) using the public code SCHNELL (Alonso et al. 2020). Here, for LVK, we consider the detector configuration including LIGO, Virgo, and KAGRA at design sensitivity, assuming no cross-correlations between detectors. The triangle shows the total SGWB energy density from compact object mergers $\Omega_{\text{GW}}(f = 25 \text{ Hz}) = 6.9^{+3.0}_{-2.1} \times 10^{-10}$ estimated by population analysis through GWTC-3 (Abbott et al. 2023a).</figcaption></figure>

for ~18 per cent of the number of Pop III mergers across cosmic history, they dominate the SGWB from Pop III BBH mergers for $\nu \lesssim 20 \text{ Hz}$ (contributing ~80 per cent of the total energy spectrum at $\nu \lesssim 10 \text{ Hz}$). The Pop III SGWB has a peak of $\Omega_{\text{GW}} \sim 10^{-10}$ at $\nu = 14 \text{ Hz}$ and converges to the standard power-law $\Omega_{\text{GW}} \propto v^{2/3}$ form for BBH inspiral at $\nu \lesssim 10 \text{ Hz}$. The energy spectrum is almost flat with $\Omega_{\text{GW}} \sim 5 \times 10^{-11}$ at $\nu \sim 30-100 \text{ Hz}$ with similar contributions from field and NSC mergers. The field mergers produce a flat energy spectrum for $\nu \sim 10-200 \text{ Hz}$, while the SGWB from NSC mergers shows a small secondary peak at $\nu \sim 100 \text{ Hz}$. The maximum ratio between the SGWB of our Pop III mergers and the total SGWB of Pop I/II and PBH mergers from Bavera et al. (2022) is $f_{\text{PopIII}} = 13.5$ per cent achieved in the low-frequency regime ($\nu \lesssim 10 \text{ Hz}$). Compared with the results for sub-components in Bavera et al. (2022), the SGWB from our Pop III BBH mergers exceeds that from Pop I/II BBH mergers in GCs at $\nu \lesssim 20 \text{ Hz}$, but remains below the SGWB from PBHs. In Fig. 7, we also show the total SGWB (including contributions from mergers involving neutron stars) $\Omega_{\text{GW}}(f = 25 \text{ Hz}) = 6.9^{+3.0}_{-2.1} \times 10^{-10}$ inferred from GWTC-3 events (Abbott et al. 2023a, see also Renzini & Golomb 2024). Actually, the total SGWB estimated by Abbott et al. (2023a) is higher than the total SGWB of Pop I/II and PBH mergers from Bavera et al. (2022) by a factor of ~5–8 in the frequency range $\nu \sim 1-200 \text{ Hz}$ that we are mostly concerned with. Therefore, if we take the results from Abbott et al. (2023a) as the reference, we arrive at a conservative estimate of the maximum contribution of Pop III BBH mergers to the total SGWB as ~3 per cent. These features of the Pop III SGWB, in

<footnote>$^{31}$The effect of BH spins is small for the total energy of GWs emitted in the inspiral phase (Lehoucq et al. 2023; Zhou et al. 2023).</footnote>

<footnote>$^{32}$In this paper, we are mostly interested in the frequency range $\nu \sim 1-200 \text{ Hz}$ accessible by ET, where most Pop III BBH mergers have their peak GW emission. Since we adopt a conservative coefficient $\kappa = 0.01$ for eccentricity enhancement by binary-single encounters, the fraction of Pop III BBH mergers in NSCs with large eccentricities ($e \gtrsim 0.1$) at $\nu \gtrsim 1 \text{ Hz}$ is expected to be small ($\lesssim 1$ per cent).</footnote>
<figure><img src="image_10.png" /><figcaption>**Figure 8.** Intrinsic chirp mass distributions of Pop III BBH mergers in the fiducial model LOG1obs, compared with those of Pop I/II mergers formed via CE (long-dashed curve), SMT (dashed-dotted curve), and in GCs (dashed curve), as well as PBH mergers (dashed-dotted curve) from the best-fitting model in Franciolini et al. (2022b) adopted by Bavera et al. (2022, see their fig. 1). The results for Pop III BBH mergers in galaxy fields and NSCs are shown with the solid and dotted contours, respectively. We also plot the intrinsic chirp mass distribution inferred from GWTC-3 events using the flexible mixture model framework (see their fig. 2 Abbott et al. 2023a) as the thick solid curve.</figcaption></figure>

particular the dominance of NSC mergers at low frequencies, can be interpreted with the different merger histories and mass distributions of field and NSC mergers (see below).

### 5.1.3 Mass distribution

In Fig. 8, we plot the intrinsic chirp distributions of Pop III BBH mergers in galaxy fields and NSCs, compared with the results for Pop I/II mergers formed by CE, SMT, and in GCs, as well as PBH mergers from the best-fitting model in Franciolini et al. (2022b) adopted by Bavera et al. (2022, see their fig. 1), and the distribution inferred from observations (see their fig. 2 Abbott et al. 2023a). The chirp mass distribution of Pop III BBH mergers in galaxy fields peaks around $m_c \sim 20 M_\odot$ and has a sharp cutoff at $m_c \sim 40 M_\odot$. The number of mergers also declines rapidly for $m_c \lesssim 10 M_\odot$. In between ($m_c \sim 10-40 M_\odot$), the distribution of Pop III field mergers is similar to that of Pop I/II BBH mergers in GCs. However, the latter is broader, shaped by repeated mergers of Pop I/II BHs (smaller than those from Pop III stars) via N-body dynamics in GCs. In galaxy fields, Pop III BBH mergers produced mostly by IBSE are generally more massive than their Pop I/II counterparts (from CE and SMT during IBSE) due to the massive and compact nature of Pop III stars. The Pop I/II mergers involve significant fractions of low-mass ($m_c \lesssim 10 M_\odot$) mergers that are rare in our case. Compared with field mergers, Pop III BBH mergers in NSCs are even more massive. Their chirp mass distribution has three peaks at $m_c \sim 30, 100, \text{ and } 300 M_\odot$. Interestingly, the location and shape of the first peak are similar to the (quasi-log-normal) chirp mass distribution of PBH mergers, but the peak of Pop III mergers is broader, which renders the second peak at $m_c \sim 100 M_\odot$ insignificant. However, there is a (narrow) gap at $m_c \sim 200 M_\odot$ between the second and third peaks.

As mentioned in Section 5.1.2, the peak contribution from a BBH merger to the SGWB satisfies $\delta\Omega_{GW}(v_p) \propto m_c/(1+z)$ with the peak location $v_p$ inversely proportional to $m_c(1+z)$ to the first order. The

<figure><img src="image_11.png" /><figcaption>**Figure 9.** Intrinsic primary (top) and secondary (bottom) mass distributions of Pop III BBH mergers in galaxy fields (solid) and NSCs (dotted) and initial BBHs (dashed) in the fiducial model LOG1obs. Here, initial BBHs denote all Pop III BBHs sampled in A-SLOTH (including those that do not merge at $z > 0$). In the top panel, we also plot the primary mass distribution inferred from the 69 confident BBH mergers in GWTC-3 (see their fig. 1 Edelman, Farr & Doctor 2023) as the thick solid curve.</figcaption></figure>

strong dependence of $\delta\Omega_{GW}(v_p)$ on $m_c$ indicates that it is the (low-z) massive NSC mergers with $m_c \gtrsim 60 M_\odot$ around the second and third peaks in the chirp mass distribution that produce the primary peak of the SGWB at $v = 14$ Hz in Fig. 7, while less massive NSC mergers shape the secondary SGWB peak at $v \sim 100$ Hz. On the other hand, due to the redshift dependence, Pop III BBH mergers in galaxy fields (with $m_c \sim 10-40 M_\odot$) make their peak contributions to the SGWB in a frequency range $v \sim 10-200$ Hz much broader than their chirp mass range given their extremely broad redshift distribution with the MRD rapidly rising from $z=0$ to $z \sim 15$ (Fig. 4). This explains the broad, flat peak of the SGWB spectrum from Pop III field mergers at $v \sim 10-200$ Hz.

To better understand the chirp mass distributions of Pop III BBH mergers from different channels, we plot the underlying primary and secondary mass distributions in Fig. 9 as well as the mass ratio distributions in Fig. 10, where we also show the results for initial Pop III BBHs (including those that do not merger at $z > 0$). The mean values of $m_c$, $m_1$, and $q$ are summarized in Table 4. In galaxy fields, almost all Pop III BBH mergers come from BHs below $50 M_\odot$. This upper bound is even significantly lower than the lower edge of the Pop III single-star PISN mass gap ($86 M_\odot$). In fact, only $\sim 0.31$ (0.14) per cent of field mergers have $m_1 > 50$ (100) $M_\odot$, among which 51 (41) per cent are ejected from NSCs, while 59.2 (35.5) per cent of NSC mergers have $m_1 > 50$ (100) $M_\odot$. The mass ratio distribution of field mergers
<figure><img src="image_12.png" /><figcaption>**Figure 10.** Intrinsic mass ratio distributions of Pop III BBH mergers in galaxy fields (solid) and NSCs (dotted) and initial BBHs (dashed) in the fiducial model LOG1obs. The mass ratio distribution inferred from the 69 confident BBH mergers in GWTC-3 (see their fig. 2 Edelman et al. 2023) is shown with the thick solid curve.</figcaption></figure>

is dominated by nearly equal-mass mergers with $q \gtrsim 0.8$, similar to the distribution inferred from GWTC-3 events (Edelman et al. 2023).

However, NSC mergers involve (both primary and secondary) BHs below, above, and also inside the PISN mass gap ~ 86–242 $M_\odot$ expected from the single-star evolution models adopted in SEVN (Costa et al. 2023). The BHs inside the mass gap are relatively rare and only cover the mass range ~ 130–242 $M_\odot$ so another narrower gap around 86–130 $M_\odot$ exists in the mass distribution of BHs in binaries, which can be regarded as a stricter definition of the Pop III PISN mass gap. This feature is also seen in other BPS studies (e.g. Tanikawa et al. 2021b, 2022b). It is the mass gap that results in the three peaks in the chirp mass distribution of NSC mergers (dotted contour in Fig. 8): (1) The low-mass peak around $m_c \sim 30 M_\odot$ is made of two BHs below the gap. (2) The intermediate-mass peak at $m_c \sim 100 M_\odot$ involves one BH below the gap and another above the gap. (3) The massive peak at $m_c \sim 300 M_\odot$ is produced by BHs above the gap. Here, the (2) mergers between one BH below and another above the mass gap produce a peak at $q \sim 0.15$ in the mass ratio distribution. The mass and mass ratio distributions of BHs in NSC mergers are similar to those for all initial BBHs. The difference is that for NSC mergers, the mass distribution is tilted to the massive end, and the peak at $q \sim 0.15$ in the mass ratio distribution is more obvious. The reason is that more massive BBHs are more likely to fall into NSCs by DF (see equations 3 and 4). Such distinct features of the NSC-DH channel with respect to IBSE, i.e. large fractions of mergers with low mass ratios and massive BHs (inside/above the mass gap), are also seen in the in situ dynamical channel for BBH mergers in massive Pop III clusters (Wang et al. 2022; Liu et al. 2024b; Mestichelli et al. 2024).

### 5.1.4 Progenitor and host system properties

The difference between the IBSE and NSC-DH channels in the mass distribution of BHs can be explained by the facts that (i) only close BBHs can merge in isolation within a Hubble time while the NSC-DH channel can also drive initially wide BBHs to merge, and (ii) the masses and initial separations of BHs are correlated as a result of binary stellar evolution. To illustrate the first point (i), we plot the initial pericenter separation $r_p$ distributions of the progenitors

<figure><img src="image_13.png" /><figcaption>**Figure 11.** Initial pericenter separation distribution for the progenitors of Pop III BBH mergers in the fiducial model LOG1obs. Shaded region: (ZAMS) stellar progenitors of field mergers. Dashed-dotted contour: stellar progenitors of NSC mergers. Solid contour: initial BBHs that become field mergers. Dotted contour: initial BBHs that become NSC mergers. Dashed contour: all initial BBHs sampled in ASLOTH. Hatched shaded region: initial BBHs that once fall into NSCs but are eventually ejected and merge in galaxy fields, which make up a tiny fraction (~0.04 per cent) of field mergers.</figcaption></figure>

of Pop III BBH mergers in Fig. 11, for both initial stellar binaries and BBHs. It can be seen that field mergers mostly originate from BBHs with initial pericenter separations $r_p \equiv a(1-e) \lesssim 0.16$ au, which are mainly produced by close stellar binaries with $r_p \lesssim 10$ au initially. However, NSC mergers mostly come from initial BBHs with $r_p \sim 0.1-10^3$ au, and the shape of the $r_p$ distribution for such NSC merger progenitors is similar to that of all BBHs for $r_p \gtrsim 0.1$ au. The corresponding progenitor stellar binaries also have a broad $r_p$ distribution extending to $r_p \sim 10^3$ au. Interestingly, there is a tiny fraction (0.04 per cent) of field mergers from initially wide BBHs. These binaries have once fallen into NSCs and been processed by DH, but are ejected later and eventually merge in galaxy fields. In general, the progenitor BBHs of field mergers and NSC mergers are approximately separated by a critical initial pericenter separation $r_{p,crit} \sim 0.1$ au, which is a common feature in all cases considered in this work (Table 2). Below this critical separation, Pop III BBHs typically merge within 1 Gyr, so there is usually not enough time for them to sink into NSCs by DF. Above the critical separation, due to the strong dependence of merger time-scale on initial separation for isolated evolution (equation 8), most BBHs cannot merge within a Hubble time in isolation.

For the second point (ii), Fig. 12 shows the distribution of all Pop III BBHs at formation in the $\log m_1 - \log a$ space from the LOG1 SEVN catalogue used in our fiducial model, where we also demonstrate the criteria for equal-mass BBHs on initially circular orbits to (a) merge in galaxy fields by GW emission for a given merger time-scale $t_{GW}$, and (b) survive binary-single encounters and merge by DH in NSCs with a given velocity dispersion $\sigma_*$ On the left side of the dotted line for $t_{GW} = 14$ Gyr, we have $m_1 \lesssim 50 M_\odot$ for most BBHs that can merge within a Hubble time in isolation. Such BBHs originate from close ($r_p \lesssim 10$ au) binaries of stars with ZAMS masses below 242 $M_\odot$, which lose their envelopes through mass transfer and/or CE phases, so the resulting naked cores can only produce BHs below ~50 $M_\odot$ (Iorio et al. 2023). Close binaries of more massive stars tend to merge during unstable mass transfer when the donor leaves the MS as a red super-giant while the accretor still
<figure><img src="image_14.png" /><figcaption>**Figure 12.** Distribution of Pop III BBHs at formation in the log $m_1 - \log a$ space for the LOG1 simulation by SEVN from Costa et al. (2023) with the close IBS model. To estimate which BBHs can merge by GW emission in galaxy fields, we plot the primary mass as a function of initial separation given the merger time-scales $t_{GW} = 14000$ (dotted), 100 (long-dashed), and 1 Myr (dashed-dotted-dotted) for equal-mass BBHs on initially circular orbits. Similarly, for the NSC-DH channel, we derive the minimum primary mass required for a BH binary to be hard in NSCs with $\sigma_* = 30$ (solid), 100 (dashed), and 300 km s$^{-1}$ (dashed-dotted) given the criterion $a < a_{HDB} = Gm_1^2 / [m_{\star,SC}\sigma_*^2]$.</figcaption></figure>

on the MS also fills the Roche lobe due to significant expansion<fnref n="33" /> (Costa et al. 2023). Initially wider ($r_p \gtrsim 10$ au) binaries undergo less and even negligible mass-loss during binary stellar evolution and produce more massive BHs up to 86 (550) M$_\odot$ from stars with ZAMS masses below (above) 242M$_\odot$.

However, the resulting BBHs are mostly too wide ($a \gtrsim 0.3$ au) to merge within a Hubble time in isolation. Nevertheless, the majority of such massive BBHs with wide orbits are still hard binaries in typical host NSCs with $\sigma_* = 30$ km s$^{-1}$ (solid line), so a significant fraction ($\sim 45–64$ per cent) of them can merge within a Hubble time by DH. In particular, BBHs with at least one BH inside or above the PISN mass gap $\sim 86–242$ M$_\odot$ originate from massive stellar binaries with $m_{\star,1} \gtrsim 242$ M$_\odot$ and form two groups in the log $m_1 - \log a$ space. The first group with $m_1 \sim 130–320$ M$_\odot$ (including all BHs inside the PISN mass gap) and $a \sim 0.1–30$ AU is produced by stellar binaries with initial pericenter separations of $r_p \sim 5–30$ au, in which the primary star loses the H envelope by CE evolution and stable mass transfer during the red super-giant phase before collapsing into a BH directly. The second group with $m_1 \sim 242–550$ M$_\odot$ (above the PISN mass gap) and $a \gtrsim 30$ AU comes from initially wider ($r_p \gtrsim 30$ AU) stellar binaries in which mass transfer/loss via Roche lobe overflow never happens, so the primary star collapses entirely into a BH (as mass-loss by stellar winds is negligible in the Pop III tracks adopted by SEVN). Almost all BBHs with BHs inside or above the PISN mass gap are massive enough to merge by DH even in massive NSCs with $\sigma_* = 300$ km s$^{-1}$ (dashed-dotted line).

Beyond distinct progenitors, the two evolution channels also produce Pop III BBH mergers in very different environments. For conciseness, here we briefly summarize the host halo/galaxy/NSC properties of Pop III BBH mergers in our fiducial model and comment on how they vary with model parameters. Most ($\gtrsim 90$ per cent) field mergers occur in small haloes with $M_h \lesssim 5 \times 10^9$ M$_\odot$ at $z \gtrsim 0.8$, which host low-mass, metal-poor galaxies with $M_\star \lesssim 1.3 \times 10^6$ M$_\odot$ and $Z \lesssim 0.22$ Z$_\odot$. The typical (median) host system has $M_h \sim 9.4 \times 10^7$ M$_\odot$, $M_\star \sim 3.5 \times 10^3$ M$_\odot$, $Z \sim 0.005$ Z$_\odot$, and $z \sim 6.8$. In fact, $\sim 10$ per cent of Pop III BBHs merge before the formation of the second generation of stars (i.e, $M_\star = 0$). Therefore, direct observations of the host galaxies of Pop III BBH mergers from IBSE are very challenging. On the contrary, the majority ($\sim 80$ per cent) of NSC mergers happen in more massive structures with $M_h \sim 6.7 \times 10^9 – 1.1 \times 10^{11}$ M$_\odot$, $M_\star \sim 1.3 \times 10^6 – 8 \times 10^8$ M$_\odot$, $Z \sim 0.017 – 0.48$ Z$_\odot$, and $M_{SC} \sim 9 \times 10^4–3.4 \times 10^9$ M$_\odot$ at $z \sim 0.1–3$. The typical host system has $M_h \sim 2 \times 10^{10}$ M$_\odot$, $M_\star \sim 1.5 \times 10^7$ M$_\odot$, $Z \sim 0.075$ Z$_\odot$, $M_{SC} \sim 4.5 \times 10^5$ M$_\odot$, and $z \sim 0.7$. Nevertheless, the host galaxies of our Pop III BBH mergers in NSCs are generally less massive than the galaxies hosting NSCs in local observations (see their fig. 12 Neumayer et al. 2020), especially for Late-Type NSC-host galaxies that mainly have $M_\star \gtrsim 10^8$ M$_\odot$. This indicates that the NSC-DH channel favours low-mass, metal-poor, compact (Early-Type) galaxies in small haloes where the inspiral of Pop III BBHs into NSCs by DF is efficient. We also find that among the host NSCs of Pop III BBHs, each NSC typically only swallows a few Pop III BBHs throughout its lifetime, so the total mass (as well as binding energy) carried by the Pop III BBHs is minor (up to a few per cent) compared with that of the NSC itself. This justifies our ignorance of the heating effect by DH of Pop III BBHs on NSC evolution.

In general, the host system properties of field mergers are insensitive to model parameters, although more massive systems at lower redshifts are slightly favoured by initially wider binaries with more bottom-heavy IMFs. However, increasing the abundances of Pop III BBHs and NSCs will shift NSC mergers to less massive systems with lower metallicities at higher redshifts. For instance, in the model LOG5_full with the highest abundances of Pop III BBHs and NSCs, the typical host system has $M_h \sim 10^{10}$ M$_\odot$, $M_\star \sim 3.6 \times 10^6$ M$_\odot$, $Z \sim 0.05$ Z$_\odot$, $M_{SC} \sim 2.8 \times 10^5$ M$_\odot$, and $z \sim 1.1$. In the opposite case of KRO1_low, we have $M_h \sim 3.7 \times 10^{10}$ M$_\odot$, $M_\star \sim 4.9 \times 10^7$ M$_\odot$, $Z \sim 0.17$ Z$_\odot$, $M_{SC} \sim 8 \times 10^5$ M$_\odot$, and $z \sim 0.3$. These trends can be explained by the distinct merger histories in different models, which are discussed below.

## 5.2 Exploration of parameter space

The general features of Pop III BBH mergers from the IBSE and NSC-DH channels described above hold in all the 18 models considered in this paper. However, the detailed properties of BBH mergers and the relative importance of the two channels do vary greatly with the underlying assumptions on (1) Pop III IMF, (2) IBS, and (3) occupation fraction of high-$z$ NSCs. In this section, we first discuss the general trends in the BBH formation and merger efficiencies (Section 5.2.1). Then we illustrate the effects of each aspect on the merger history, SGWB, and mass distribution (Section 5.2.2–5.2.4).

### 5.2.1 BBH formation and merger efficiencies

The BBH formation and merger efficiencies of the 18 simulations are summarized in Table 2 and visualized in Fig. 13. Under the

<footnote n="33">³³If massive Pop III stars remain compact during MS (by, e.g. suppressed convective overshooting or chemically homogeneous evolution), some of these close binaries can undergo stable mass transfer without stellar mergers and produce massive BBHs (with BHs above 130 M$_\odot$) merging within a Hubble time (Hijikawa et al. 2021; Tanikawa et al. 2021a, b, 2022b; Santoliquido et al. 2023; Mestichelli et al. 2024). This indicates that the fates of massive binary stars and the efficiency of forming massive BBH mergers via IBSE are sensitive to the modelling of single star evolution.</footnote>
<figure><img src="image_15.png" /><figcaption>**Figure 13.** Pop III BBH formation and merger efficiencies from the 18 models listed in Table 2. The BBH formation efficiency (crosses) is independent of the NSC parameters, so we only show 6 data points for the 6 input SEVN catalogues. The efficiency of field mergers (diamonds) is insensitive to the NSC parameters, and we only show the results for the fiducial NSC model (obs). For NSC mergers, we show the results for the obs, full, and low NSC models (Section 3.1) with triangles, circles, and squares, respectively. The left (right) section of the plot shows the results for the close (wide) IBS model (Section 2.1). Within each section, the IMF becomes more bottom-heavy (with higher $\alpha$) from left to right. We have $\epsilon_{GW}^{\text{field}} \simeq \epsilon_{GW}^{\text{NSC}}$ in the model TOP1_full, causing the two data points to overlap.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>**Figure 14.** Same as Fig. 12 but for the LOG5 simulation by SEVN from Costa et al. (2023) with the wide IBS model.</figcaption></figure>

wide IBS model in which close binary interactions are rare, $\epsilon_{BBH}$ is not very sensitive to the IMF, and the fiducial log-flat IMF shows the highest $\epsilon_{BBH}$. This is produced by two competing effects: When individual stars become more massive (with increasing $\alpha$), the number of stars/binaries is simply smaller for a fixed total stellar mass, but meanwhile, the fraction of stars massive enough to form BHs is larger. However, under the close IBS model where binary interactions are important, $\epsilon_{BBH}$ becomes higher when the IMF is more bottom-heavy (with higher $\alpha$). The reason is that more massive stars are more likely to merge during close binary interactions due to their larger radii, so a larger fraction of stellar binaries that would form BBHs without interactions will be lost in mergers for a larger $\alpha$. Similarly, when there are more wide binaries in the IBS model, mergers of massive stars are suppressed, resulting in higher $\epsilon_{BBH}$ and more massive BBHs, as shown in Fig. 14.

For Pop III BBH mergers in galaxy fields, the merger efficiency $\epsilon_{GW}^{\text{field}}$ is sensitive to the NSC parameters since only a small ($\sim 3–12$ per cent) fraction of Pop III BBHs fall into NSCs. Therefore, we focus on the results of $\epsilon_{GW}^{\text{field}}$ for the fiducial NSC model (obs). We find that $\epsilon_{GW}^{\text{field}}$ always increases with $\alpha$, and the trend is stronger with the close IBS model. The evolution of $\epsilon_{GW}^{\text{field}}$ with $\alpha$ is also stronger compared with the case of $\epsilon_{BBH}$, as the fraction $f_{GW}^{\text{field}}$ of BBHs in galaxy fields that merge at $z > 0$ increases with $\alpha$. These trends can be explained by the fact that field mergers mostly originate from interacting close ($r_p \lesssim 10$ au) binaries of relatively small stars ($\lesssim 240$ M$_\odot$) while close binaries of massive stars will merge, as discussed in Section 5.1.4 (see Fig. 12). For the same reason, with a fixed IMF, $\epsilon_{GW}^{\text{field}}$ is significantly reduced (by a factor of $\sim 6–14$) when we switch from the close IBS model to that dominated by wide binaries. The reduction is stronger when the IMF is more bottom-heavy because more massive stars are more likely to interact given a fixed separation.

For NSC mergers, the variation of merger efficiency $\epsilon_{GW}^{\text{NSC}}$ with IMF and IBSE is very similar to that of $\epsilon_{BBH}$. Indeed, if the probability for a BBH to merge in a NSC at $z > 0$ is constant, $\epsilon_{GW}^{\text{NSC}}$ should be proportional to $\epsilon_{BBH}$. In fact, the probability is slightly higher for more massive BBHs that sink into NSCs more efficiently by DF and require less hardening to merge rapidly via GW emission. As a result, on top of the trend driven by $\epsilon_{BBH}$, $\epsilon_{GW}^{\text{NSC}}$ is further enhanced by a more top-heavy IMF or IBS with more wide binaries, which tend to produce more massive BBHs. This point is reflected in the weak evolution of the fraction $f_{GW}^{\text{NSC}}$ of BBHs in NSCs that merge at $z > 0$ and the fraction $f_{\text{infall}}$ of BBHs that fall into NSCs (see the last two columns of Table 2) with IMF and IBS. As expected, $f_{\text{infall}}$ and thus $\epsilon_{GW}^{\text{NSC}}$ increase when the occupation fraction of NSCs is higher. Combing these distinct features of field and NSC mergers, we find that the fraction $f_{\text{NSC}}$ of mergers in NSCs across cosmic history is higher when the IMF is more top-heavy, the initial condition includes more wide stellar binaries, and high-$z$ NSCs are more abundant. The most important factor here is IBS: For the close IBS model, field mergers are more abundant than NSC mergers in most cases with $f_{\text{NSC}} \sim 8–51$ per cent except for the TOP1_full model in which the two populations of mergers have almost identical numbers. However, for the wide IBS model, NSC mergers always dominate with $f_{\text{NSC}} \sim 65–95$ per cent. Besides, the BBH merger efficiency of the NSC-DH channel is less sensitive to the IMF and IBS compared with that of the IBSE channel.

### 5.2.2 Effects of initial binary statistics

To illustrate the effects of IBS, as an example, we compare the results from LOG1obs (i.e. the fiducial model) and LOG5obs with the close and wide IBS models (see Section 2.1), respectively, where the IMF and NSC models are fixed to the fiducial choices. Consistent with the trends seen in merger efficiencies $\epsilon_{GW}$ (Fig. 13), the MRDs of field and NSC mergers are reduced and increased when there are more wide binaries in the initial condition (LOG5obs), respectively, as shown in Fig. 15. As a result, NSC mergers dominate the MRD at $z \lesssim 6$ with the wide IBS model. Similarly, the reduction (enhancement) of the MRD of field (NSC) mergers is weaker (stronger) when the IMF is more bottom-heavy. The MRD redshift evolution for NSC mergers is not significantly affected by IBS, which holds for other choices of IMFs and NSC models. This indicates that the NSC-DH channel is mainly sensitive to the abundances of Pop III BBHs rather than their initial properties, which is also reflected in the correlation between $\epsilon_{GW}^{\text{NSC}}$ and $\epsilon_{BBH}$. For field mergers, the MRD evolution is
<figure><img src="image_17.png" /><figcaption>**Figure 15.** Same as Fig. 4 but comparing the MRDs for two IBS models dominated by close (blue) and wide (orange) binaries assuming the log-flat IMF and fiducial NSC occupation fraction.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>**Figure 16.** Same as Fig. 7 but comparing the energy spectra of SGWB for the two IBS models dominated by close (blue) and wide (orange) binaries assuming the log-flat IMF and fiducial NSC occupation fraction.</figcaption></figure>

insensitive to IBS at $z \gtrsim 5$, which always shows a peak at $z \sim 15$ closely following the peak of the Pop III SFRD (Fig. 3), since this regime is dominated by mergers with short delay times. However, at lower redshifts when Pop III star formation has terminated, and the MRD is dominated by mergers of long delay times, the decrease of MRD with decreasing redshift is slower for the wide IBS model compared with the case of the close IBS model. The reason is that the wide IBS model produces a larger fraction of BBHs with long merger time-scales ($t_{GW} \gtrsim 300$ Myr) (see fig. 11 in Costa et al. 2023). This effect is stronger when the IMF is more top-heavy. In fact, the close and wide IBS models prefer different binary evolution pathways to BBH mergers in isolation, which produce different delay time distributions (see fig. 12 and 17 in Costa et al. 2023). The former favours stable mass transfer during early evolutionary stages, while the latter is dominated by CE evolution of red giants with large radii ($\sim 1–10$ au). The reader is referred to Costa et al. (2023) and Santoliquido et al. (2023) for in-depth discussions on the roles played by different binary evolution pathways in shaping the properties of BBH mergers via IBSE.

Fig. 16 shows the contributions of Pop III BBH mergers to the SGWB for the two IBS models. Compared with LOG1obs with

<figure><img src="image_19.png" /><figcaption>**Figure 17.** Same as Fig. 8, but comparing the chirp mass distributions of Pop III BBH mergers in galaxy fields (middle) and NSCs (bottom) for the two IBS models dominated by close (blue shaded region) and wide (orange contour) binaries assuming the log-flat IMF and fiducial NSC occupation fraction. Here, we only show the results for Pop I/II BBH mergers in GCs (dashed curve) and PBH mergers (dashed-dotted-dotted curve) from the best-fitting model in Franciolini et al. (2022b) adopted by Bavera et al. (2022, see their fig. 1), and the intrinsic chirp mass distribution (thick solid curve) inferred from GWTC-3 events from (see their fig. 2 Abbott et al. 2023a).</figcaption></figure>

the close IBS model, the contribution from field (NSC) mergers is lower (higher) for LOG5obs with the wide IBS model, where NSC mergers completely dominate the Pop III SGWB at $\nu \lesssim 200$ Hz. This dominance of NSC mergers is a common feature of the wide IBS model regardless of the chosen IMF and NSC parameters. This can be understood by the fact that the SGWB (observed at $z=0$) captures the summation of the GW energies from all mergers across cosmic history and more massive mergers at lower $z$ contribute higher energies to the SGWB (Section 5.1.2). As mentioned in the last section and shown in Fig. 15, the merger efficiency and local MRD of NSC mergers are always larger than those of field mergers under the wide IBS model. Moreover, NSC mergers are more massive than field mergers, and thus release more energies in GWs (Section 5.1.3). Therefore, in LOG5obs, the contribution from NSC mergers is higher than that from field mergers by about 100 times in the Pop III SGWB, being as strong as the SGWB from PBH mergers predicted by Bavera et al. (2022) at $\nu \lesssim 10$ Hz. In this case, Pop III NSC mergers produce a significant flattening feature around $\nu \sim 10–20$ Hz deviating from the canonical shape $\Omega_{GW} \propto f^{2/3}$ in the total SGWB from stellar BBH mergers (Pop I/II + III).

Finally, we compare the chirp mass distributions between the two models in Fig. 17 (see also Table 4). We find that both field and NSC mergers become more massive with the wide IBS model, and the
<figure><img src="image_20.png" /><figcaption>**Figure 18.** Same as Fig. 4, but comparing the MRDs for the three IMFs with power-law slopes of $\alpha = 1$ (blue), 0.17 (orange), and 2.3 (green), assuming the close IBS and fiducial NSC models.</figcaption></figure>

effect is stronger for NSC mergers. The average chirp mass increases by 10 and 35 per cent for field and NSC mergers given the wide IBS model compared with the case of the close IBS model. On the one hand, mass-loss and stellar mergers during close binary interactions are suppressed in the wide IBS model, which enhances the formation of massive BBHs (see Figs 12 and 14) and thus greatly increases the masses of NSC mergers (whose mass distribution closely correlates with the mass distribution of all BBHs at formation, see Fig. 9). Besides, the fraction of NSC mergers with small mass ratios $q \lesssim 0.4$ is higher in the wide IBS model. On the other hand, since close binary interactions are required to produce BBHs that can merge within a Hubble time in isolation, field mergers mostly involve low-mass BHs ($\lesssim 50 M_\odot$) regardless of the IBS model. Therefore, changing from the close IBS model to the wide IBS model only shifts the primary mass distribution of field mergers to the massive end by a few solar masses (but still limited by $m_1 \lesssim 50 M_\odot$), leading to slightly larger chirp masses.

### 5.2.3 Effects of IMF

Next, we look into the effects of IMF variations using the results from LOG1obs, TOPobs, and KRO1obs with IMF slopes of $\alpha = 1, 0.17$, and 2.3 (see Section 2.1), respectively, where we adopt the close IBS model and fiducial NSC occupation fraction. As shown in Fig. 18, when the IMF is more bottom-heavy (with a larger $\alpha$), the MRD of field (NSC) mergers is significantly (mildly) increased. For both field and NSC mergers, the MRD redshift evolution is insensitive to the IMF, while the normalization of the MRD curve is approximately proportional to the merger efficiency (Fig. 13), whose dependence on the IMF is explained in Section 5.2.1. These trends reflect that the IMF has minor effects on the delay time distributions for both field and NSC mergers, and that the Pop III formation history does not vary greatly with $\alpha$, as shown in Fig. 19. In fact, the (co-moving) stellar mass density of Pop III stars ever formed in our simulation box is $\rho_{*,III} \simeq 4.6, 5.6$ and $5.7 \times 10^4 M_\odot Mpc^{-3}$ for $\alpha = 1, 0.17$, and 2.3, respectively. The small variation of the Pop III star formation history with $\alpha$ justifies our simple modelling of Pop III stellar feedback (Section 2.1). Our $\rho_{*,III}$ values are generally consistent with the predictions from cosmological simulations (e.g. Johnson et al. 2013; Smith et al. 2015; Xu et al. 2016; Sarmento et al. 2017; Liu & Bromm 2020a, b; Incatasciato, Khochfar & Oñorbe

<figure><img src="image_21.png" /><figcaption>**Figure 19.** Same as Fig. 3 but comparing the SFRDs of Pop III (dashed) and Pop I/II (solid) stars from LOG1obs (blue), TOPobs (orange), and KRO1obs (green) with Pop III IMF slopes of $\alpha = 1, 0.17$, and 2.3 under the close IBS model and the fiducial NSC model.</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>**Figure 20.** Same as Fig. 7, but comparing the energy spectra of SGWB for the three IMFs with power-law slopes of $\alpha = 1$ (blue), 0.17, (orange), and 2.3 (green), under the close IBS model and the fiducial NSC model.</figcaption></figure>

2023; Venditti et al. 2023) and other semi-analytical models (e.g. Dayal et al. 2020; Visbal, Bryan & Haiman 2020; Hegde & Furlanetto 2023; Feathers et al. 2024; Ventura et al. 2024) and the upper limit $\sim 10^{5-6} M_\odot Mpc^{-3}$ placed by reionization (Visbal, Haiman & Bryan 2015; Inayoshi et al. 2016).

The SGWB contributed by Pop III BBH mergers in galaxy fields also increases with $\alpha$, as shown in Fig. 20. The shape of the energy spectrum at $\nu \sim 15-200$ Hz varies slightly with the IMF. Compared with the flat spectrum for $\alpha = 1$, $\Omega_{GW}$ slowly increases (decreases) with $\nu$ for $\alpha = 2.3$ (0.17). The reason is that BBH mergers become slightly more massive when the IMF is more top-heavy (see Fig. 21 and Table 4), which enhances GW emission at lower frequencies. However, for NSC mergers, the SGWB increases with $\alpha$ at $\nu \gtrsim 30$ Hz, but decreases with $\alpha$ at $\nu \lesssim 30$ Hz. Here, the high-frequency regime is dominated by low-mass mergers with $m_c \lesssim 60 M_\odot$ whose fraction increases with $\alpha$ as shown in Fig. 21. Since the merger efficiency and local MRD of NSC mergers also increase with $\alpha$ (Figs 13 and 18), it is reasonable that the SGWB from such low-mass mergers is stronger. On the other hand, the low-frequency part of the SGWB is produced by massive mergers with $m_c \gtrsim 60 M_\odot$
<figure><img src="image_23.png" /><figcaption>**Figure 21.** Same as Figs 8 and 17, but comparing the chirp mass distributions of Pop III BBH mergers in galaxy fields (middle) and NSCs (bottom) for the three IMFs with power-law slopes of $\alpha = 1$ (blue shaded region), 0.17 (orange contour), and 2.3 (green contour), again under the close IBS and fiducial NSC models.</figcaption></figure>

whose fraction is reduced when the IMF is more bottom-heavy (see Fig. 21 and Table 4). This reduction overcomes the increase of merger number with $\alpha$ and drives the decrease of SGWB with $\alpha$.

### 5.2.4 Effects of NSC occupation fraction

Finally, we consider the effects of the NSC occupation fraction in Figs 22 and 23 with the wide IBS model and log-flat IMF. We find that the properties of field mergers are hardly affected by the NSC parameters, which is reasonable as only a small ($\lesssim 12$ per cent) fraction of Pop III BBHs ever fall into NSCs. Therefore, we focus on NSC mergers in this section. The chirp mass distributions of NSC mergers are very similar in the three NSC models and not shown here for conciseness (see Table 4 for the average chirp masses). As expected, the MRD of NSC mergers is enhanced when more galaxies host NSCs. Unlike the fiducial NSC model where the MRD of NSC mergers flattens at $z \lesssim 4$, the MRD of NSC mergers shows a peak at $z \sim 4$ in the optimistic NSC model (where every galaxy with $M_* > 10^6 M_\odot$ hosts a NSC). The MRD keeps increasing rapidly towards lower redshifts in the conservative NSC model, producing a local MRD that is slightly higher than that in the fiducial NSC model. The difference between the fiducial and conservative NSC models in the SGWB is very small, which can be understood with their similar merger efficiencies and local MRDs, although the peak value of the SGWB in the conservative NSC model is slightly ($\sim 16$ per cent) higher than that in the fiducial NSC model due to the slightly higher local MRD. The SGWB is stronger in the optimistic NSC model

<figure><img src="image_24.png" /><figcaption>**Figure 22.** Same as Fig. 4 but comparing the MRDs for the fiducial (blue), optimistic (orange), and conservative (green) models of NSC occupation fraction under the wide IBS model and log-flat IMF.</figcaption></figure>

<figure><img src="image_25.png" /><figcaption>**Figure 23.** Same as Fig. 7 but comparing the energy spectra of SGWB for the fiducial (blue), optimistic (orange), and conservative (green) models of NSC occupation fraction under the wide IBS model and log-flat IMF.</figcaption></figure>

as expected from its higher merger efficiency and local MRD. In this case (LOG5_full), the Pop III SGWB becomes comparable ($\sim 84$ per cent) to the SGWB from Pop I/II and PBH mergers predicted by Bavera et al. (2021) at $v \lesssim 10$ Hz, which may leave a detectable feature in the total SGWB (Périgois et al. 2021; Martinovic et al. 2022).

## 5.3 Observational perspective

We estimate the detection rates $\dot{N}_{ET,SNR>8}$ of Pop III BBH mergers by ET for all the 18 models (Table 2) given the detection probability of each source with SNR > 8 derived from the python package GWTOOLBOX (Yi et al. 2022a, b), using the formalism of equation (24). For simplicity, we do not distinguish high-$z$ and low-$z$ sources but sum up all events detectable by ET. We also do not consider eccentricities in this calculation as in the case of SGWB (Section 5.1.2). The results are summarized in Fig. 24, where we only show the detection rates of field mergers with the fiducial NSC model since the results for the other two NSC models are almost identical. As discussed in Section 5.1.1, most ($\gtrsim 90$ per cent) BBH mergers with $z < 20$ will be detectable by ET thanks to its high
<figure><img src="image_26.png" /><figcaption>**Figure 24.** Detection rates of Pop III BBH mergers by ET with SNR > 8 for the 18 models listed in Table 2. The detection rates of field mergers (diamonds) are insensitive to the NSC parameters, and we only show the results for the fiducial NSC model (obs). For NSC mergers, we show the detection rates for the obs, full, and low NSC models (Section 3.1) with triangles, circles, and squares, respectively. The left (right) section of the plot shows the results for the close (wide) IBS model (Section 2.1). Within each section, the IMF becomes more bottom-heavy (with higher $\alpha$) from left to right.</figcaption></figure>

sensitivity across a broad frequency range, so the all-sky merger rate $\dot{N}(z < 20)$ within $z < 20$ is a good approximation of the ET detection rate (for SNR > 8). Therefore, we derive the fraction of $\dot{N}(z < 20)$ contributed by Pop III BBHs based on the reference rates for Pop I/II and PBH mergers from Franciolini et al. (2022b) to estimate the chance of finding Pop III BBHs by ET. The results are summarized in Table 3.

In general, the ET detection rate of all Pop III BBH mergers is $\dot{N}_{ET,SNR>8} \sim 50-1370 \text{ yr}^{-1}$, and we have $\dot{N}_{ET,SNR>8} \sim 30-1230 \text{ yr}^{-1}$ and $\dot{N}_{ET,SNR>8} \sim 6-230 \text{ yr}^{-1}$ for field and NSC mergers, respectively. The highest rate is achieved by KRO1_full, in which case Pop III BBHs contribute $\sim 6.5$ per cent of the all-sky merger rate within $z < 20$. The dependence of detection rates on model parameters follows similar trends as those seen in merger efficiencies (Fig. 13), which are explained in Section 5.2.1. However, here the effects of NSC parameters are stronger, and field mergers become more important. The reason is that the detection rate is more sensitive to the high-$z$ MRD than the merger efficiency. IBS still plays the most important role in determining the relative importance of the NSC-DH and IBSE channels. Given the close IBS model, field mergers always dominate the total detection rate, as the contribution of NSC mergers remains low ($\sim 0.7-22$ per cent). With the wide IBS model, the NSC mergers account for $\sim 13-88$ per cent of the total detection rate. Their contributions are only higher than those of field mergers in 5 of the 9 models with relatively high NSC occupation fractions and large $\alpha$, although they always dominate the total number of mergers across cosmic history with $f_{NSC} \sim 65-95$ per cent.

As discussed in Section 5.1.3, an important feature of the NSC-DH channel is the ability to produce Pop III BBH mergers with massive ($\gtrsim 50 \text{ M}_\odot$) BHs, especially those above 130 $\text{M}_\odot$. In light of this, we further look into the ET detection rates $\dot{N}_{ET,IMBH,SNR>8}$ of Pop III mergers with at least one IMBH ($m_1 > 100 \text{ M}_\odot$). We find that such mergers are extremely rare in galaxy fields, with ET detection rates below $0.1 \text{ yr}^{-1}$ in all the 18 models considered here. However, we have $\dot{N}_{ET,IMBH,SNR>8} \sim 0.5-200 \text{ yr}^{-1}$ for NSC mergers as shown in Fig. 25. The completeness of detection of such mergers by ET is also very high ($\gtrsim 95$ per cent). Here, the detection rate is

<figure><img src="image_27.png" /><figcaption>**Figure 25.** Same as Fig. 24 but for the detection rates of Pop III BBH mergers including at least one IMBH ($m_1 > 100 \text{ M}_\odot$) in NSCs. The detection rates of such mergers in galaxy fields remain below $0.1 \text{ yr}^{-1}$ and are therefore not shown here.</figcaption></figure>

sensitive to all model parameters on Pop III IMF, IBS, and occupation fraction of high-$z$ NSCs, implying that BBH mergers with IMBHs$^{34}$ of $\sim 100-1000 \text{ M}_\odot$ are valuable for constraining the properties of Pop III binaries and high-$z$ SCs with GW observations.

We obtain similar detection rates for CEO$^{35}$: $\dot{N}_{CEO,SNR>8y} \sim 30-1270 \text{ yr}^{-1}$ for field mergers, $\dot{N}_{CEO,SNR>8} \sim 6-210 \text{ yr}^{-1}$ for NSC mergers, and $\dot{N}_{CEO,SNR>8} \sim 0.5-120 \text{ yr}^{-1}$ for mergers involving IMBHs. The main difference here is that the detection probability of massive sources by CEO is lower at $z \gtrsim 2$, because CEO is less sensitive than ET at low frequencies ($\nu \lesssim 10$ Hz, see fig. 8 of De Luca et al. 2021). The completeness of detection of NSC mergers (with IMBHs) can be as low as $\sim 70 (60)$ per cent for CEO (in TOP5_full), while it remains $\gtrsim 95$ per cent for ET.

Finally, we estimate the detection rate of Pop III BBH mergers by the LVK network during the O4 run with SNR > 8 as $\dot{N}_{LVK,SNR>8} \sim 0.9-9 \text{ yr}^{-1}$. This small detection rate is consistent with the small local MRD $\dot{n}(z=0) \sim 0.034-0.17 \text{ yr}^{-1} \text{ Gpc}^{-3}$ that is much lower than the observed local MRD $\dot{n}_{obs}(z=0) = 19.3_{-9}^{+15.1} \text{ yr}^{-1} \text{ Gpc}^{-3}$ (Abbott et al. 2021). This indicates that Pop III mergers are overwhelmed by other populations of BBH mergers in the low-$z$ Universe (see also Table 3). Although Pop III stars can still be important for the most massive mergers, such events are rare in LVK observations. For instance, the LVK O4 detection rate of Pop III mergers with IMBHs (in NSCs) is $\dot{N}_{LVK,IMBH,SNR>8} \sim 0.04-1.5 \text{ yr}^{-1}$, and the expected number of events detected in 18 months$^{36}$ of observations is only $\sim 0.07-2.2$. Considering the typical horizon redshift $z_{horizon} \sim 1$ of the LVK network during O3 for BBH mergers with total masses $\sim 100-300 \text{ M}_\odot$ (see fig. 8 of De Luca et al. 2021), we derive the MRD of Pop III mergers with IMBHs for $z < 1$ as $\dot{n}_{IMBH}(z < 1) \sim 0.003-0.15 \text{ yr}^{-1} \text{ Gpc}^{-3}$, which is lower than the upper limits derived from O3 data under most conditions (see their fig. 2 Abbott et al. 2022). For Pop III BBH mergers with both BHs in the standard PISN

<footnote>$^{34}$Because of PISNe, all IMBHs produced by Pop III stars in our models are above 130 $\text{M}_\odot$ (see Fig. 9), which is a common feature of the adopted BPS catalogues from SEVN (Costa et al. 2023), and also found in other BPS studies (Hijikawa et al. 2021; Tanikawa et al. 2021a, 2022b).</footnote>

<footnote>$^{35}$We use the CEO second stage (wideband) sensitivity curve (see fig. 2 of Abbott et al. 2017) from https://dcc.ligo.org/LIGO-P1600143/public (curve_data.txt).</footnote>

<footnote>$^{36}$https://observing/docs.ligo.org/plan/</footnote>
mass gap 50–130 M⊙ like GW190521 (detected at $z = 0.82_{-0.34}^{+0.28}$ $\lesssim$ 1, Abbott et al. 2020b), we predict the MRD to be $\dot{n}_{gap}(z < 1) \sim 0.0013–0.034$ yr$^{-1}$ Gpc$^{-3}$ for $z < 1$, also consistent with the rate $0.08_{-0.07}^{+0.19}$ yr$^{-1}$ Gpc$^{-3}$ estimated for GW190521-like mergers from observations (Abbott et al. 2022). In fact, we find that GW190521 can be marginally explained by our Pop III BBH mergers (in NSCs) assuming the log-flat or bottom-heavy IMF and the wide IBS model. In general, the results from all the 18 models are consistent with current observations of BBH mergers, but it will be challenging to derive strong constrains on Pop III BBH mergers from current and upcoming observations by the LVK network that is only sensitive to low-$z$ events among which Pop III mergers are expected to be sub-dominant (unless Pop III star formation is much more efficient than that modelled by A-SLOTH in our case). This highlights the importance of the 3rd-generation GW detectors like ET and CEO in observing Pop III BBH mergers at high $z$.

# 6 SUMMARY AND DISCUSSION

We study the roles played by two evolution channels in producing BBH mergers from Pop III) stars using the public semi-analytical model A-SLOTH (Hartwig et al. 2022; Magg et al. 2022a) applied to the halo merger trees produced by the cosmological simulation from Ishiyama et al. (2016). The first channel considers IBSE in which interactions of close ($\lesssim$ 10 au) binary stars shrink binary orbits to form tight ($\lesssim$ 0.1 au) BBHs that can merge within a Hubble time by gravitational wave (GW) emission in galaxy fields (i.e. effectively in isolation). In the second channel, Pop III BBHs fall into NSCs by DF and are subsequently driven to merge via DH by binary-single encounters, which also works for initially wide (up to $10^5$ au) Pop III BBHs given their massive nature. We combine the binary population synthesis (BPS) results of SEVN (Costa et al. 2023) with customized routines for the formation, dynamics, and internal orbital evolution of Pop III BBHs in galaxy fields and NSCs (Fig. 1), the formation and evolution of NSCs (Fig. 2), as well as halo growth and galaxy evolution at lower redshifts ($z \lesssim 5$, Section 4) in A-SLOTH based on the earlier work by Liu & Bromm (2021). We explore 18 models (Table 2) under different assumptions on Pop III IMF, IBS (see Section 2.1), and occupation fraction of high-$z$ NSCs (see Section 3.1). Our main findings are summarized as follows:

(i) Although only a small fraction (3–12 per cent) of Pop III BBHs fall into NSCs, a significant fraction ($\sim$ 45–64 per cent) of the BBHs in NSCs merge at $z > 0$. The NSC-DH channel can be as efficient as the classical IBSE channel for producing Pop III BBH mergers at $z \lesssim 8$, while the IBSE channel always dominates the merger rate density at $z \gtrsim 10$, producing a peak at $z \sim 15$ that closely follows the peak of Pop III star formation. The NSC-DH channel contributes $f_{NSC} \sim 8–95$ per cent of Pop III BBH mergers across cosmic history, and is generally more important at lower redshifts, accounting for $\mathcal{F}_{NSC} \sim 34–98$ per cent (0.7–88 per cent) of the all-sky merger rates of Pop III BBH mergers within $z < 1$ (20). Higher contributions from the NSC-DH channel are achieved by initially wider binary stars, more top-heavy IMFs, and higher occupation fractions of NSCs.

(ii) The most important factor that determines the relative importance of the two channels is IBS. When the IBS is dominated by close binaries, the NSC-DH channel produces $f_{NSC} \sim 8–51$ per cent of all Pop III mergers across cosmic history and explains $\mathcal{F}_{NSC} \sim 0.7–22$ per cent of the all-sky merger rate within $z < 20$. However, when the IBS is dominated by wide binaries, we have $f_{NSC} \sim 65–95$ per cent and $\mathcal{F}_{NSC} \sim 12–88$ per cent for $z < 20$. These outcomes are mainly driven by the strong dependence of the

IBSE channel on IBS. The NSC-DH channel is relatively insensitive to the initial properties of Pop III binary stars (i.e. IMF and IBS).

(iii) Pop III BBH mergers in NSCs are more massive than those from IBSE. The latter mostly involve BHs below $\sim 50$ M⊙ from the collapse of naked cores of stellar progenitors that have relatively low initial masses ($\lesssim 240$ M⊙) and lose their hydrogen envelopes during binary interactions. As pointed out in Mestichelli et al. (2024), binaries of more massive stars that can potentially form more massive BHs either merge during unstable mass transfer when they are initially too close or do not produce tight enough BBHs that can merge within a Hubble time in isolation when they are initially too far apart. However, the NSC-DH channel can produce BBH mergers from wide ($\gtrsim 0.1$ au) binaries of more massive BHs (in the mass ranges 50–86 and 130–550 M⊙) whose stellar progenitors start on wide ($\gtrsim 5$ au) orbits and experience less and even negligible mass-loss during binary stellar evolution. For instance, a significant fraction ($\sim$ 4–84 per cent) of Pop III BBH mergers in NSCs involve at least one intermediate mass BH (IMBH) above 100 M⊙, while such mergers are extremely rare from IBSE.

(iv) We find that the ET (Punturo et al. 2010; Maggiore et al. 2020) is able to detect most ($\gtrsim$ 90 per cent) Pop III BBH mergers in our models with SNR above 8. The total detection rate is estimated as $\mathcal{N}_{ET,SNR>8} \sim 50–1370$ yr$^{-1}$, and we have $\mathcal{N}_{ET,SNR>8} \sim 30–1230$ yr$^{-1}$ and $\mathcal{N}_{ET,SNR>8} \sim 6–230$ yr$^{-1}$ for the IBSE and NSC-DH channels, respectively. Moreover, ET will detect $\sim 0.5–200$ Pop III BBH mergers involving at least one IMBH ($m_1 > 100$ M⊙) per year from the NSC-DH channel (mostly at $z \lesssim 10$). The detection rate of such mergers from IBSE is below 0.1 yr$^{-1}$. In our simulations, Pop III stars only make up a tiny fraction ($\lesssim 10^{-4}$) of the total mass budget of stars ever formed in the universe. However, we estimate that Pop III BBH mergers will make up a much larger fraction ($\sim$ 0.3–6.5 per cent) of the BBH mergers that the 3rd-generation GW detectors will observe, given the conservative predictions on the merger rates of BBHs from metal-enriched Pop I/II stars and PBHs by Franciolini et al. (2022b) as the reference.

(v) Our results are consistent with current observations of BBH mergers by the LVK network, especially for the merger rates of IMBHs and BHs with masses in the standard PISN mass gap (50–130 M⊙) like those in GW190521 (Abbott et al. 2020b, 2022). However, it will be challenging to identify and characterize Pop III mergers with observations of the LVK network that are mostly sensitive to the low-$z$ regime ($z \lesssim 1$) where Pop III mergers in our models are completely sub-dominant (accounting for $\lesssim 1$ per cent of the total merger rate). Even for massive BBH mergers with IMBHs for which Pop III stars are expected to make higher contributions, the estimated detection rate of such mergers from Pop III remnants (mostly in NSCs) during the O4 run of the LVK network is very low ($\mathcal{N}_{LVK,IMBH,SNR>8} \sim 0.04–1.5$ yr$^{-1}$). This indicates that the 3rd-generation detectors like ET are required to efficiently constrain Pop III star formation (and high-$z$ SCs) with GW observations (Iwaya et al. 2023; Franciolini et al. 2024; Santoliquido et al. 2024).

(vi) The SGWB produced by Pop III BBH mergers has a peak value of $\Omega_{GW} \sim 10^{-11}–8 \times 10^{-11}$ around observer-frame frequencies $\nu_{peak} \sim 10–100$ Hz. Mergers in NSCs dominate the Pop III SGWB at $\nu \lesssim 20$ Hz in most cases due to their massive nature.<fnref n="37" />

<footnote n="37">³⁷Interestingly, our predictions on the SGWB from Pop III BBH mergers at $\nu \lesssim 20$ Hz (dominated by the NSC-DH channel in most cases) are similar to previous results (Périgois et al. 2021; Martinovic et al. 2022) that only consider the IBSE channel and assume much higher Pop III star formation</footnote>
The SGWB from Pop III BBH mergers is equal to ~2–32 per cent of the total SGWB of compact object mergers inferred from observations (Abbott et al. 2023a) at $v \lesssim 10$ Hz. Therefore, it is very difficult to observe such a sub-dominant signal from Pop III BBHs in the total SGWB. However, the contribution of Pop III mergers can be larger in the residual background observed by the 3rd-generation detectors like ET where individual detected sources can be properly subtracted (Zhong et al. 2024), under certain conditions for the other BBH merger populations (Périgois et al. 2021; Martinovic et al. 2022; Kouvatsos & Sakellariadou 2024). Moreover, if the cosmic star formation rate density of Pop III stars is higher than that predicted by A-SLOTH in our case by a factor of a few, Pop III BBHs in NSCs can dominate the SGWB at $v \lesssim 10$ Hz. Considering the cumulative Pop III stellar mass density $\rho_{*,III} \sim 5 \times 10^4 M_{\odot} Mpc^{-3}$ found in our simulations, such enhancement of Pop III star formation is still allowed by the constraints $\rho_{*,III} \lesssim 10^{5-6} M_{\odot} Mpc^{-3}$ from reionization under conservative assumptions for the production/escape of ionizing photons from Pop III stars (Visbal et al. 2015; Inayoshi et al. 2016).

In general, the IBSE and NSC-DH channels produce Pop III BBH mergers of distinct properties and shape the Pop III contributions to the GW signals from compact object mergers in different ways. In particular, an important feature of the NSC-DH channel is that it can produce mergers with IMBHs ($> 100 M_{\odot}$) and BHs in the (standard) PISN mass gap (50–130 $M_{\odot}$) efficiently. This is consistent with previous studies showing that such massive BBH mergers form efficiently via dynamical processes in SCs (Rodríguez et al. 2015, 2019; Arca Sedda et al. 2020, 2021; Di Carlo et al. 2020; Fragione & Silk 2020; Fragione et al. 2020, 2022; Gerosa & Fishbach 2021; Khan & Holley-Bockelmann 2021; Kimball et al. 2021; Liu & Lai 2021; Mapelli et al. 2021, 2022; Anagnostou et al. 2022; Chattopadhyay et al. 2023; Fragione & Rasio 2023; Bruel et al. 2024; Dall’Amico et al. 2024; Liu et al. 2024b; Mestichelli et al. 2024; Wang et al. 2022), which can also provides seeds for supermassive BHs (e.g. Askar, Davies & Church 2022; Rose et al. 2022; Kritos, Berti & Silk 2023; Kritos et al. 2024; Rantala, Naab & Lahén 2024).

Our results show strong, complex dependence on the largely unknown properties of Pop III binary stars and high-$z$ NSC clusters. In addition to the uncertain aspects explored in the current paper, there are other uncertainties/caveats in our modelling that may have significant effects on our results. Below we discuss the ones known to date in the hope that they provide directions for future work.

(1) The properties of Pop III BBHs at birth are crucial inputs in our model, which not only depend on initial conditions but are also highly sensitive to the parameters/algorithms governing binary interactions and single stellar evolution (see e.g. Iorio et al. 2023). Our conclusions are based on the BPS results from Costa et al. (2023) using SEVN with a specific set of parameters, and thus can be altered if other BPS results are considered. For instance, the finding that the IBSE channel hardly produces mergers with BHs above 100 $M_{\odot}$ results from stellar mergers caused by significant expansion of massive Pop III stars during the MS phase in the stellar evolution tracks adopted by SEVN (see Section 5.1.3). However, as shown in Tanikawa et al. (2021b), if massive Pop III stars remain compact during MS evolution, massive BBH mergers (involving BHs above 130 $M_{\odot}$) can still form via IBSE.

This indicates that significant degeneracy exists in the large parameter space of Pop III BBH mergers to explain a single observable.

(2) Our modelling of the dynamics of Pop III BBHs in host galaxies and during halo mergers under spherical symmetry is highly idealized. The adopted assumptions on the galaxy size and structure (Arca-Sedda & Capuzzo-Dolcetta 2014) are based on local observations and may not hold at high redshifts. Besides, we assume that Pop III BBHs are always on their own, i.e. not embedded by larger structures (e.g. satellite galaxies, bound clusters of Pop III stars and remnants) that can sink into galaxy centres more efficiently by dynamical friction, so the efficiency of the NSC-DH channel may be underestimated. Similarly, we assume that the smaller galaxies hosting NSCs during galaxy mergers are immediately destroyed, leaving behind naked NSCs.

(3) We use a phenomenological model of NSCs with adiabatic evolution based on empirical scaling relations from local observations (Neumayer et al. 2020). These relations may not be valid at high redshifts and cannot capture dynamical NSC assembly by mergers of globular clusters and young star clusters and in-situ star formation, which may affect the evolution of BBHs in NSCs. Moreover, we only consider dynamical hardening, softening/disruption and ejection by binary-single encounters to follow the internal orbital evolution of BBHs in NSCs, ignoring higher-order processes (e.g. relativistic phase space diffusion, tides-driven eccentricity excitation, and Kozai–Lidov mechanism) involving tidal fields, general relativity effects, and interactions with other components (e.g. single/binary Pop I/II BHs, central massive BHs and discs of active galactic nuclei) in NSCs (see Section 2.2.2).

(4) The merger trees adopted in our simulations only cover the high-$z$ regime ($z \gtrsim 4.5$) where the limited simulation volume is cosmologically representative. At lower redshifts, we use a simple galaxy evolution model to keep track of Pop III BBHs, calibrated to reproduce the observed star formation rate density within a factor of 2. This model only considers smooth halo growth and star formation without taking into account halo/galaxy mergers which are important for the infall of BBHs into NSCs. How exactly this simplification affects our results is unclear and related to the assumptions on the dynamics of BBHs and NSCs during galaxy mergers.

(5) We use a stochastic model for the formation of Pop III BBHs that is not self-consistently connected to the stellar feedback routine, since we do not use the same population of binary and single stars to model stellar feedback, which is instead derived from a separate population of single stars. Besides, we use the fiducial star formation and stellar feedback parameters of A-SLOTH (calibrated against observations) in all our simulations. These parameters may need re-calibrations to be consistent with observations when we change the IMF and include binary stars. In general, the cosmic star formation rate density of Pop III stars is poorly constrained without direct observations, and current theoretical predictions have a scatter of ~1.5 dex (see e.g. Liu & Bromm 2020b; Hartwig et al. 2022; Klessen & Glover 2023). Considering other Pop III star formation histories will change the merger histories of Pop III BBHs, especially for the IBSE model (Santoliquido et al. 2023).

In this paper, we focus on the merger rate, SGWB, and mass distribution of Pop III BBH mergers, whereas the spins of merging BHs (with respect to the orbital vector) can also provide clues for the formation/evolution channels of BBHs (e.g. Biscoveanu et al. 2022; Callister et al. 2022; Mould et al. 2022; Adamcewicz, Lasky & Thrane 2023; Guo et al. 2024; Pierra, Mastrogiovanni & Perriès 2024). We plan to investigate the spin distributions of Pop III BBH mergers from different channels in future work. Here, we briefly comment on the physics involved. In the standard pathway of Pop III star cluster formation (in minihaloes of $M_h \sim 10^{5-6} M_{\odot}$) via disc
fragmentation (Liu et al. 2021a; Klessen & Glover 2023; ), the initial stellar spins and binary orbital vectors are expected to be aligned as they both inherit the angular momentum of the same star-forming disc. On the other hand, Pop III binary stars with misaligned spins can form through non-standard pathways with turbulent fragmentation and mergers of sub-clusters in more massive haloes (under the influence of peculiar environmental effects such as baryon-dark matter streaming motion, Hirano et al. 2018b, 2023). The stellar spins are further regulated by stellar winds, tidal effects, and mass transfer during (binary) stellar evolution leading to diverse outcomes (e.g. alignment, tilt, and flipping, Kinugawa et al. 2020; Stegmann & Antonini 2021; Tanikawa et al. 2021b). Next, the angular momentum transport in stellar interiors (before and during the collapse) and SN natal kicks determine the BH spin magnitudes and orientations, which are still highly uncertain in current theoretical models (e.g. Jiménez-Forteza et al. 2017; Qin et al. 2018, 2019; Fuller & Ma 2019; Fuller, Piro & Jermyn 2019; Bavera et al. 2020; Belczynski et al. 2020; Olejak & Belczynski 2021; Ghodla & Eldridge 2023). Finally, if the BBHs reside in dense star clusters, the BH spins and orbit plane can also be altered by dynamical interactions, among which (fly-by) binary-single encounters typically cause small tilts (~15°) of the orbit plane while exchanges tend to result in isotropic orientations (e.g. Bouffanais et al. 2019; Trani et al. 2021; Banerjee, Olejak & Belczynski 2023; Dall'Amico et al.).

In conclusion, our results indicate that valuable information on the initial properties of Pop III binary stars and the evolution channels/environments of BBHs is encoded in the GWs from Pop III BBH mergers. The former can be fundamentally different from the observed properties of Pop I/II binaries in the local Universe, and capture the elusive interplay between fragmentation, (proto)stellar feedback, competitive accretion, migration, gravitational scatters, and mergers of protostars during Pop III star (cluster) formation, as well as the dynamics of Pop III star clusters (Sakurai et al. 2017; Hirano & Bromm 2018; Chon & Hosokawa 2019; Sugimura et al. 2020, 2023; Liu et al. 2021a, 2024b; Riaz et al. 2022; Wang et al. 2022; Franchini et al. 2023; Park et al. 2023, 2024; Mestichelli et al. 2024). The latter are governed by key processes in galaxy evolution at Cosmic Dawn, such as star cluster formation and galactic dynamics which have been increasingly revealed by JWST observations (e.g. Baggen et al. 2023; Ferreira et al. 2023; Kartaltepe et al. 2023; Langeroodi & Hjorth 2023; Adamo et al. 2024; Ito et al. 2024; Ji et al. 2024a; Ormerod et al. 2024). Although distinguishing Pop III BBHs from other populations of mergers on an individual basis is non-trivial (Franciolini et al. 2022a; Costa et al. 2023), identifying a subgroup of BBH mergers of Pop III origins from a large enough sample of mergers is feasible and promising given the unique features of Pop III BBHs and the non-proportionally high contributions to GW signals by Pop III stars (Tanikawa et al. 2022b; Iwaya et al. 2023; Franciolini et al. 2024; Santoliquido et al. 2024). In the next decades, the 3rd-generation GW detectors like ET and CEO (Reitze et al. 2019; Evans et al. 2023) will observe thousands of events per year reaching z ~ 30. The 'gravitational-wave archaeology' in synergy with other observations of Cosmic Dawn will provide us new insights into the first stars, star clusters, and galaxies in the early Universe.

## ACKNOWLEDGEMENTS

The authors thank the referee, Tomoya Kinugawa, for insightful comments that helped improve our paper. We thank Michela Mapelli for providing the input SEVN catalogues underlying Costa et al. (2023) and valuable discussions. We thank Giuliano Iorio for developing and making SEVN available. We thank Benedetta Mestichelli for useful

discussions on Pop III binary stellar evolution. We thank Mattis Magg and Li-Hsin Chen for their help in code development with A-SLOTH. We thank Simone S. Bavera for providing the data underlying Bavera et al. (2022). We thank Sebastian Ljung for feedback on the BBH and NSC modules in A-SLOTH. The authors acknowledge the Texas Advanced Computing Center (TACC) for providing HPC resources under FRONTERA allocation AST22003. BL is supported by the Royal Society University Research Fellowship. NS gratefully acknowledges the support of the Research Foundation - Flanders (FWO Vlaanderen) grant 1290123N. GC acknowledges support from the Agence Nationale de la Recherche grant POPSYCLE number ANR-19-CE31-0022. FS acknowledges financial support from the AHEAD2020 project (grant agreement no. 871158). GC and FS acknowledge financial support from the European Research Council for the ERC Consolidator grant DEMOBLACK, under contract no. 770017. RSK acknowledges financial support from the European Research Council via the ERC Synergy Grant 'ECOGAL' (project ID 855130), from the German Excellence Strategy via the Heidelberg Cluster of Excellence (EXC 2181 - 390900948) 'STRUCTURES', and from the German Ministry for Economic Affairs and Climate Action in project 'MAINN' (funding ID 50002206). RSK thanks for computing resources provided by the Ministry of Science, Research and the Arts (MWK) of the State of Baden-Württemberg through bwHPC and the German Science Foundation (DFG) through grants INST 35/1134-1 FUGG and 35/1597-1 FUGG, and also for data storage at SDS@hd funded through grants INST 35/1314-1 FUGG and INST 35/1503-1 FUGG. This work excessively used the public packages NUMPY (van der Walt, Colbert & Varoquaux 2011), MATPLOTLIB (Hunter 2007), and SCIPY (Virtanen et al. 2020). The authors wish to express their gratitude to the developers of these packages and to those who maintain them.

## DATA AVAILABILITY

A-SLOTH is publicly available at https://gitlab.com/thartwig/asloth. The version used in this paper (including the new modules for compact object mergers and NSCs) is public at https://gitlab.com/Treibeis/a-sloth-cob. SEVN is publicly available at https://gitlab.com/sevncodes/sevn. The data underlying this paper will be shared on reasonable request to the corresponding author.

## REFERENCES

* Abbott B. P. et al., 2017, Class. Quantum Gravity, 34, 044001

* Abbott B. P. et al., 2018, Phys. Rev. Lett., 120, 091101

* Abbott R. et al., 2020a, Phys. Rev. Lett., 125, 101102

* Abbott R. et al., 2020b, ApJ, 900, L13

* Abbott R. et al., 2021, ApJ, 913, L7

* Abbott R. et al., 2022, A&A, 659, A84

* Abbott R. et al., 2023a, Phys. Rev. X, 13, 011048

* Abbott R. et al., 2023b, Phys. Rev. X, 13, 041039

* Acernese F. et al., 2015, Class. Quantum Gravity, 32, 024001

* Adamcewicz C., Lasky P. D., Thrane E., 2023, ApJ, 958, 13

* Adamo A. et al., 2024, Nature, 632, 513

* Ajith P. et al., 2011, Phys. Rev. Lett., 106, 241101

* Alonso D., Contaldi C. R., Cusin G., Ferreira P. G., Renzini A. I., 2020, Phys. Rev. D, 101, 124048

* Anagnostou O., Trenti M., Melatos A., 2022, ApJ, 941, 4

* Angel R. et al., 2008, ApJ, 680, 1582

* Antonini F., Rasio F. A., 2016, ApJ, 831, 187

* Antonini F., Gieles M., Gualandris A., 2019, MNRAS, 486, 5008

* Antonini F., Romero-Shaw I. M., Callister T., 2024, preprint (arXiv:2406.19044)