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

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 |$\sim 8\!-\!95{{\ \rm 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 ^{\rm peak}_{\rm GW}\sim 10^{-11}\!-\!8\times 10^{-11}$| around observer-frame frequencies |$\nu \sim 10\!-\!100\ \rm Hz$|⁠. The Pop III contribution can be a non-negligible (⁠|$\sim 2\!-\!32{{\ \rm per\ cent}}$|⁠) component in the total SGWB at |$\nu \lesssim 10\ \rm Hz$|⁠. The estimated detection rates of Pop III BBH mergers by the Einstein Telescope are |$\sim 6\!-\!230$| and |$\sim 30\!-\!1230\ \rm 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\ \rm M_\odot$| by the Einstein Telescope is |$\sim 0.5\!-\!200\ \rm yr^{-1}$| in NSCs but remains below |$0.1\ \rm yr^{-1}$| for IBSE.

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 & 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\ \rm m$|⁠) extraterrestrial telescopes (Angel et al. 2008; Rhodes et al. 2020; Schauer, Drory & Bromm 2020). Detection of Pop III features¹ 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² 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 electro-magnetic 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³ 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_{\star }\sim 20\!-\!40\ \rm 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 stellar-mass compact objects (mostly BHs with |$\sim 5\!-\! 100\ \rm M_\odot$|⁠) within |$z\lesssim 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\sim 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 (⁠|$\sim 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 (⁠|$\lesssim 10^{-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 (⁠|$\gtrsim 20\ \rm M_\odot$|⁠) 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_{-0.34}^{+0.28}$| with unusual BH masses |$85_{-14}^{+21}$| and |$66_{-18}^{+17}\ \rm M_\odot$| (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.⁴ GW190521 is special because the mass range |$50\!-\!130\ \rm M_{\odot }$| 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 |$\sim 100\ \rm M_\odot$| 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\lesssim 10\ \rm \mathrm{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\gtrsim 100\ \rm \mathrm{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 (⁠|$\sim 10^{4}\!-\!10^{7}\ \rm M_{\odot }$|⁠) 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, |$\gtrsim 100\ \rm M_{\odot }$|⁠) 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.⁵

To fully unleash the power of the GW probe, deeper and more systematic investigations of the relevant physics of BBH formation 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⁶ (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ęcioł, 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.⁷

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 (⁠|$z\lesssim 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.

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 (⁠|$\sim 1\!-\!2\times 10^{7}$|⁠) 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 _{\rm 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_{\rm BBH}$| and the average initial mass of BBH progenitor stellar binaries |$\bar{m}_{\rm 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 ,\rm III}$| of Pop III stars formed at the current time-step t. We assume that on average a fraction |$f_{\rm 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⁸ 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_{\rm B}$| is drawn from a Poisson distribution with parameter |$f_{\rm boost}f_{\rm B}f_{\rm BBH}\delta M_{\star ,\rm III}/\bar{m}_{\rm p, BBH}$|⁠, where |$f_{\rm boost}=50$| is a boost factor introduced to obtain better statistics. If |$N_{\rm B}\gt 0$|⁠, we randomly draw |$N_{\rm B}$| BBH progenitor binary stars from the input BPS catalogue, which will become BBHs at |$t+\tau _{\rm 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]\ \rm 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).

• A log-flat distribution |$p(m_{\star ,1})\equiv \mathrm{d}N/\mathrm{d}m_{\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.

• 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}\gt 0.5\ \rm M_\odot$|⁠) in local observations.

• A top-heavy distribution |$p(m_{\star ,1})\propto m_{\star ,1}^{-0.17}e^{-m^{2}_{\rm cut}/m_{\star ,1}^{2}}$| given |$m_{\rm cut}^{2}=20\ \rm 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_{\star }\equiv m_{\star ,2}/m_{\star ,1}$|⁠, orbital period P, which is characterized by the variable |$\pi \equiv \log (P/{\rm 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}\gt 2.2\ \rm M_\odot$| is also imposed in both cases.

• 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.

• 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 (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_{\rm h}\sim 10^{7}\ \rm 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_{\rm h}\sim 10^{5-6}\ \rm 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 ,\rm III}$|⁠, a-sloth also samples individual Pop III stars from an input power-law IMF |$p(m_{\star })\propto m_{\star }^{-\alpha _{\rm 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 _{\rm s}=\alpha =1$| for LOG (fiducial), 2.3 for KRO, and 0.17 for TOP,⁹ with |$m_{\star }\in [5,550]\ \rm 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ötberg 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{{\ \rm 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\ \rm 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_{\rm vir}$| of the (post-merger) host halo, its apocentre distance r is drawn randomly from a (cumulative) probability distribution of |$x\equiv r/R_{\rm vir}$|⁠:

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_{\star }$| 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 inspirial of Pop III BBHs in galaxies hosting NSCs, as such galaxies reside in relatively massive haloes with |$M_{\rm h}\gtrsim 10^{9}\ \rm 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_{\rm 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).

2.2.1.2 Inspiral by dynamical friction

Within the host galaxy, we follow the inspiral of each binary with

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 _{\rm DF}$| is the DF time-scale formula (see below), |$\vec{\lambda }_{\chi }$| captures the dark matter halo properties (⁠|$M_{\rm vir}$| and |$R_{\rm vir}$|⁠), and |$\vec{\lambda }_{\star }$| denotes parameters of the background stellar system, which are the mass |$M_{\star }$|⁠, size |$R_{\star }$| and inner slope of density profile |$\gamma _{\star }$| of the host galaxy in our case. Here, |$M_{\star }$| 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_{\star }$| 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}\ \rm M_\odot$|⁠), and |$\tau _{\rm 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.¹⁰

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

where |$\ln \Lambda \sim \ln [M_{\star }r/(0.8mR_{\star })]$| is the Coulomb logarithm, and |$v\sim \sigma _{\star }\sim \sqrt{GM_{\star }/(0.8R_{\star })}$|⁠. 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):

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\lt \gamma _{\star }\lt 2$| and equation (3) is used otherwise.¹¹

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_{\star }$| 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):

where |$M_{\star ,11}\equiv M_{\star }/(10^{11}\ \rm M_{\odot })$|⁠, and the inner slope |$\gamma _{\star }$| is generated from a uniform distribution in the ranges of |$0\!-\!2$|⁠. We also impose a lower limit of |$0.02R_{\rm 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_{\star }\sim 10^8\!-\!10^{11}\ \rm 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_{\star })$| predicted by equation (5) with a random number (drawn once¹² and inherited along the main branch¹³) 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_{\star }$| can be regarded as an upper limit. It will be shown below (equation 4) that the dependence of DF time-scale |$\tau _{\rm DF}$| on galaxy size is weak (⁠|$\tau _{\rm DF}\propto R_{\star }^{-0.26}$|⁠), which, combined with the weak mass and redshift dependence of |$R_{\star }$| in equation (5), implies that our results are insensitive to the modelling of |$R_{\star }(M_{\star })$|⁠. Actually, we find by numerical experiments that changing |$R_{\star }$| 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_{\rm vir}r/(0.8mR_{\rm vir})]$| and |$v\sim \sqrt{GM_{\rm vir}/(0.8R_{\rm vir})}$|⁠, implicitly assuming a power-law slope |$\gamma _{\chi }\gt 2$| for the dark matter density profile,¹⁴ 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.¹⁵

We further introduce a free parameter |$M_{\star ,\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 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_{\star }\sim 10^{10}\ \rm 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_{\star ,\min }=10^{6}\ \rm 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_{\star }\lesssim 10^{9}\ \rm 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 _{\rm B}$| after the initial formation of the progenitor binary stars. The earlier stage (⁠|$t\lt \tau _{\rm 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¹⁶ is used to evolve a and e in SCs, we further do sub-cycling with time-steps |$\delta t_{\rm binary}\le 0.01 a/(\mathrm{d}a/\mathrm{d}t)$|⁠. A merger happens when the expected time taken to merge has passed since BBH formation or |$a\lt 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¹⁷ (Peters 1964):

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

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

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_{\rm 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_{\rm GW,ej}(a_{\rm ej},e_{\rm ej})$| calculated from the orbital parameters (⁠|$a_{\rm ej}$| and |$e_{\rm 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\lt R_{\rm 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\lt R_{\rm SC}$| and |$a\lt a_{\rm HDB}=Gm_{1}m_{2}/\left[m_{\star ,\rm SC}\sigma _{\star }(r)^{2}\right]$| (Mapelli et al. 2021), given the SC size |$R_{\rm SC}$|⁠, local velocity dispersion |$\sigma _{\star }(r)$|⁠, and typical mass of objects in the cluster. Throughout this study, we adopt |$m_{\star ,\rm SC}=0.5\ \rm M_{\odot }$| for simplicity, although |$m_{\star ,\rm 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\gt a_{\rm 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):

where |$\rho _{\star }(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.

As in Arca Sedda (2020), here a SC of mass |$M_{\rm SC}$| is characterized by a cored Dehnen sphere (Dehnen 1993) with

where |$x\equiv r/R_{\rm SC}$|⁠, |$\delta =\gamma _{\rm SC}$| for |$\gamma _{\rm SC}\lt 1$| and |$\delta =2-\gamma _{\rm SC}$| for |$\gamma _{\rm SC}\ge 1$|⁠, given the scale length |$R_{\rm SC}$| and inner slope |$\gamma _{\rm SC}$| of the Dehnen sphere, which are derived from |$M_{\rm 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 _{\rm c}$| that determines the core size |$r_{\rm c}=\max (R_{\rm SC}/\Delta _{\rm c}^{1/\gamma _{\rm SC}}, r_{\bullet })$|⁠, which is no smaller than the influence radius of the BBH |$r_{\bullet }=R_{\rm SC}(2m/M_{\rm SC})^{1/(3-\gamma _{\rm SC})}$|⁠. We adopt an optimistic core overdensity |$\Delta _{\rm 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 _{\star }(r)$|⁠, |$\sigma _{\star }(r)$|⁠, and |$\tau _{\rm 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_{\rm GW}\lt a\lt a_{\rm ej}$|⁠, given the critical separations for GW-dominance |$a_{\rm GW}$| and ejection |$a_{\rm ej}$| (Miller & Lauburg 2009; Antonini & Rasio 2016; Fragione & Silk 2020; Mapelli et al. 2021):

where |$\rho _{\star }(R_{\rm SC})$| is the typical stellar density evaluated at the scale length |$R_{\rm 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\ \rm M_{\odot }$|⁠) in low-mass (⁠|$M_{\rm SC}\lesssim 10^{5}\ \rm M_{\odot }$|⁠) SCs where |$a_{\rm ej}\gt a_{\rm GW}$|⁠. In our case, the dynamical hardening time-scale in such low-mass (⁠|$M_{\rm SC}\lesssim 10^{5}\ \rm 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}\ \rm \mathrm{au}$|⁠) is much larger than that required for efficient evolution in fields by GW emission (⁠|$a\lesssim 0.5\ \rm \mathrm{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; Fahrion 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 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 }\gt M_{\star ,\rm min}$|⁠) is governed by the input occupation fraction of NSCs as a function of galaxy mass |$f_{\rm occ}(M_{\star })\equiv f_{\rm occ}(\tilde{M}_{\star })$|⁠, |$\tilde{M}_{\star }\equiv \log (M_{\star }/\rm M_{\odot })$|⁠. We adopt the fit formula in Liu & Bromm (2021, see their fig. 8)

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

In the GC accretion scenario of NSC formation relevant for low-mass (⁠|$M_{\star }\lesssim 10^{9}\ \rm M_{\odot }$|⁠) galaxies, |$f_{\rm 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}\ \rm M_{\odot }$| at |$z\gtrsim 5$| (see also Biernacki, Teyssier & Bleuler 2017; Pfister et al. 2019; Bortolas et al. 2020).¹⁹

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:

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

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

• In the conservative model (labelled by ‘low’), we reduce the efficiency of DF²⁰ (as well as NSC formation) with |$f_{\rm occ}=\hat{f}_{\rm occ}$| and |$M_{\star ,\min }=\min \lbrace 1,[(1+z)/6]^{5.14}\rbrace \times 10^{10}\ \rm M_{\odot }$|⁠, assuming that |$M_{\star ,\min }$| evolves from |$10^{10}\ \rm M_{\odot }$| at |$z=5$| to |$10^{6}\ \rm M_{\odot }$| at |$z=0$| following a power law of |$(1+z)$|⁠.

Given the input |$f_{\rm occ}(M_{\star })$|⁠, 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}\lt f_{\rm occ}(M_{\star ,i})$|⁠, which indicates the presence of an NSC, where |$M_{\star ,i}$| is the current galaxy mass of node i. Given |$p_{i}\lt f_{\rm occ}(M_{\star ,i})$|⁠, if this node contains no SCs, we create a new SC object as its NSC, whose properties (i.e. mass |$M_{\rm SC}$|⁠, size |$R_{\rm SC}$|⁠, and inner slope of density profile |$\gamma _{\rm SC}$|⁠, see equation 11) are derived from the empirical NSC–galaxy scaling relations,²¹ 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_{\star ,i}$| as the input.²² 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_{\star ,i}$|⁠, assuming adiabatic evolution of NSCs regulated solely by host galaxy mass. A NSC host galaxy will lose its NSC if |$p_{i}\gt f_{\rm occ}(M_{\star ,i})$|⁠, which will only happen for |$M_{\star }\gt 10^{9.5}\ \rm 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 |$\hat{f}_{\rm occ}$| with |$M_{\star }$| for |$M_{\star }\gt 10^{9.5}\ \rm 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_{\rm SC}$| also shrinks and its eccentricity |$e_{\star }$| decays under DF according to equation (2), now replacing m with |$M_{\rm 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:

where |$n_{0}$| is fixed by |$4\pi \int _{0}^{R_{\rm vir}}n_{\rm SC}(r)r^{2}\mathrm{d}r=N_{\rm SC}$|⁠, given the halo virial radius |$R_{\rm vir}$| and number of (satellite) SCs |$N_{\rm SC}$|⁠. Also, |$R_{\rm core}=\min (31.62R_{\star }M_{\star , 11}^{1/6}, R_{\rm 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\rightarrow \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_{\rm SC}(r)r^{2}/N_{\rm SC}$|⁠. The eccentricity |$e_{\star }$| 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\lt R_{\rm SC}$|⁠. Immediately after the merger, the Pop III BBHs contained in the normal SC are placed at |$r=R_{\rm SC}$| and their external orbital eccentricities |$e_{\star }$| 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²³

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

where |$\ln \Lambda ^{\prime }=0.11 N$|⁠, |$r_{\rm g}=\max (r, R_{\rm SC})$|⁠, |$v_{\rm g}=\sqrt{GM_{\star }/(0.8R_{\star })}$| is the typical circular velocity of the host galaxy of mass |$M_{\star }$| and scale length |$R_{\star }$| (see their equation 18 Arca-Sedda & Capuzzo-Dolcetta 2014), and |$N=M_{\rm SC}/m_{\star ,\rm 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_{\rm SC}\propto t^{2/3}$| at |$t\rightarrow \infty$| (Portegies Zwart et al. 2010), if ignoring mass-loss, given |$t_{\rm th}\propto R_{\rm SC}^{3/2}$|⁠. For simplicity, we assume that |$\gamma _{\rm SC}$| remains constant. A SC evaporates when |$M_{\rm SC}$| drops below |$100\ \rm 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_{\rm NSC}=1\ \rm 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}\rm Mpc)^{3}$| with a mass resolution of |$5000\ h^{-1}\rm M_{\odot }$| down to |$z_{\rm 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_{\rm 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_{\rm f}$| and grow smoothly at the cosmic average rate as a function of halo mass |$M_{\rm h}$| and redshift z, inferred from the millennium cosmological simulations (Boylan-Kolchin et al. 2009) for |$\Lambda \rm CDM$| (Fakhouri et al. 2010):

Once the halo mass |$M_{\rm h}(z)$| is known, we can derive the galaxy (stellar) mass with

where |$M_{\star ,\rm SHMR}(M_{\rm 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_{\star }\Omega _{m}/(\Omega _{b}M_{\rm 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²⁴ under the same halo mass. Such stochastic effects are captured by a-sloth and can cause an offset between |$M_{\star ,\rm SHMR}$| and |$M_{\star }(z_{\rm f})$| predicted by a-sloth at |$z_{\rm f}$|⁠. Assuming that this offset decays after cosmic noon at |$z_{\rm noon}=2$|⁠, for each halo we evolve B from the initial condition |$\left[(B+1)M_{\star ,\rm SHMR})\right]|_{z=z_{\rm f}}=M_{\star }(z_{\rm f})$| with

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\gt B_{\max }+1=25$|⁠, whose halo/galaxy mass remains constant.²⁵ We also consider that in the post-reionization epoch (⁠|$z\lt 6$|⁠), haloes with |$M_{\rm h}\lt 6.7\times 10^{8}\ \mathrm{ 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_{\star }\lt \eta _{\max }(\Omega _{b}/\Omega _{m})M_{\rm peak}$|⁠, given |$M_{\rm peak}$| the maximum mass the halo has ever reached and |$\eta _{\max }=0.5$| the maximum star formation efficiency.

We adopt fixed (star formation) time-steps of |$\delta t_{\rm post}=1\ \rm Myr$| for the low-z evolution of galaxies/haloes with the above formalism (equations 21–23) where the values of |$z_{\rm noon}$|⁠, |$\xi$|⁠, |$B_{\max }$|⁠, and |$\eta _{\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\rightarrow 0$| compared with observations. The reason is that we did not include explicitly 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\lt z_{\rm f}$| are likely rather optimistic.

To better characterize the host galaxies of Pop III mergers at |$z\lt z_{\rm f}$|⁠, we also estimate the average gas-phase metallicity Z of each galaxy by extrapolating from the value at |$z=z_{\rm f}$| predicted by a-sloth with the redshift-dependent mass–metallicity relation |$Z\propto M_{\star }^{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_{\rm com}=(8\ h^{-1}\rm Mpc)^{3}\simeq 1650\ Mpc^{3}$| with a dark matter mass resolution of |$5000\ h^{-1}\rm M_{\odot }$|⁠, which only cover |$z \ge z_{\rm 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)²⁶ 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_{\rm B}=0.7$|⁠. Since various stochastic processes are involved in our modelling and the simulation volume underlying the merger trees is not very large,²⁷ we boost the number of Pop III binaries to be sampled in a-sloth by a factor of |$f_{\rm boost}=50$| to achieve better statistics and reduce numerical noise in the merger history, so the effective simulation volume |$\tilde{V}_{\rm com}=50 V_{\rm 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 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 (⁠|$\lesssim 2{{\ \rm 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²⁸ and the LVK network during O4²⁹ using the python package gwtoolbox³⁰ (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 this model, most (⁠|$\sim 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 |$\sim 35\ \rm M_\odot$| and a width of |$\sigma =0.41$|⁠, and make up |$f_{\rm 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\gt 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_{\rm occ}=\hat{f}_{\rm occ}$| (see equation 15) with |$M_{\star ,\min }=10^{6}\ \rm M_{\odot }$|⁠. 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_{\rm BBH}\simeq 12.6{{\ \rm per\ cent}}$|⁠, which (given |$f_{\rm B}=0.7$| and |$\bar{m}_{\rm p,BBH}\simeq 193\ \rm M_\odot$|⁠) corresponding to a BBH formation efficiency (number of BBHs formed per unit stellar mass) |$\epsilon _{\rm BBH}\simeq 4.78\times 10^{-4}\ \rm 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\gt 0$|⁠. The overall merger efficiency (number of mergers at |$z\gt 0$| per unit stellar mass) is |$\epsilon _{\rm GW}^{\rm all}=6.49\times 10^{-5}\ \rm M_{\odot }^{-1}$|⁠. 17.9 per cent of Pop III mergers occur in NSCs across cosmic history, among which only |$\sim 2{{\ \rm 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 |$\sim 60$| per cent of Pop III BBH mergers in NSCs.

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}_{\rm obs}(z= 0)=19.3_{-9}^{+15.1}\ \rm yr^{-1}\ 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}/(\tilde{V}_{\rm 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}_{\rm peak}\sim 3.6\ \rm yr^{-1} 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_{\star }\gt 10^{6}\ \rm M_{\odot }$|⁠) 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_{\rm 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_{\rm delay}=\tau _{\rm B}+t_{\rm 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_{\rm delay}\sim 5\ \rm Myr$|⁠. However, it takes at least |$\sim 200\ \rm Myr$| for Pop III BBHs to fall into NSCs in the first place. For BBH mergers in NSCs, we have |$t_{\rm delay}\gtrsim 300\ \rm Myr$|⁠, and the delay time distribution is almost uniform (i.e. |$\mathrm{d}N/\mathrm{d}\log t_{\rm delay}\propto t_{\rm delay}$|⁠) at |$t_{\rm 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|$\sim 5$| Gyr.

The predicted local MRD of Pop III BBH mergers is |$\sim 0.06\ \rm yr^{-1}\ Gpc^{-3}$|⁠, which only counts for |$0.31_{-0.13}^{+0.27}{{\ \rm 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 |$\sim 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_{\rm horizon}$| can be written as (Hartwig et al. 2016)

where |$\mathrm{d}V/\mathrm{d}z=4\pi cd_{\rm c}^{2}(z)/H(z)$| is the co-moving volume element, given the co-moving distance |$d_{\rm c}(z)=\int _{0}^{z}\mathrm{d}z^{\prime }[c/H(z^{\prime })]$| 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 {\dot{N}}_{\rm 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}_{\rm X}(z_{j})=\sum _{i=1}^{N_{j}}p_{{\rm det},i}$| for a redshift bin |$z_j$| given the detection probabilities |$p_{{\rm det},i}$| of the |$N_{j}$| mergers within the bin. In the calculation of the all-sky merger rate |$\mathcal {\dot{N}}$|⁠, we simply use |$p_{{\rm 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_{{\rm det},i}$| from gwtoolbox (Yi et al. 2022a, b), ET can detect most (⁠|$\gtrsim 90{{\ \rm per\ cent}}$|⁠) Pop III BBH mergers with |$\rm SNR\gt 8$|⁠. Similarly, most Pop I/II and PBH mergers in the reference model from Franciolini et al. (2022b) are detectable by ET with |$\rm SNR\gt 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\lt 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 {\dot{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 |$\dot{\mathcal {N}}\sim 650\ (9)\ \rm yr^{-1}$| within |$z\lt 20\ (1)$|⁠, which counts for |$\sim 3\ (0.4){{\ \rm 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_{\rm horizon}=20$| and 1. The contribution of NSC mergers saturates at |$z_{\rm 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}\gt 100\ \rm M_\odot$|⁠) |$\mathcal {N}_{\rm IMBH}=11.6\ (1.34)\ \rm yr^{-1}$| within |$z\lt 7\ (1)$|⁠. The Pop III BBH all-sky merger rate remains below that of PBHs by a factor of |$\sim 10$| for |$z_{\rm horizon}\lesssim 20$|⁠, and becomes comparable to the rate for Pop I/II BBHs in GCs for |$z_{\rm horizon}\gtrsim 20$|⁠. Here, the PBH merger rate keeps increasing with |$z_{\rm horizon}$|⁠, while the GC merger rate saturates around |$z_{\rm 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\ \rm 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}_{\rm com}=50 V_{\rm 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):

in which |$\nu$| is the observer-frame frequency, |$\rho _{\rm c}=3H_{0}^{2}/(8\pi G)$| is the critical density, |$H(z)$| is the Hubble parameter at redshift z, |$\mathrm{d}E_{\rm GW}/\mathrm{d}\nu _{\rm s}$| is the source-frame GW energy spectrum of an individual source evaluated at the source-frame frequency |$\nu _{\rm 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 \mathrm{d}z\int \mathrm{d}\theta \dot{n}(z)p(\theta ,z)=\lim _{\Delta z_{j}\rightarrow 0} \sum _{j=1}^{N_{z}} \Delta z_{j}\sum _{i=1}^{N_{j}}\Delta t_{j}^{-1}\tilde{V}_{\rm 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 |$\mathrm{d}E_{{\rm GW},i}/\mathrm{d}\nu _{\rm s}$| are added up, where |$\nu _{\rm s}=(1+z_{i})\nu$|⁠, and |$\Delta t_{j}$| is the corresponding cosmic age evolution at redshift bin j. Finally, since |$\mathrm{d}z/\mathrm{d}t=(1+z)H(z)$|⁠, we arrive at the last line with |$\lim _{\Delta z_{j}\rightarrow 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 |$\mathrm{d}E_{{\rm GW},i}/\mathrm{d}\nu _{\rm s}$|⁠, assuming random inclinations and a fixed effective spin³¹|$\chi _{\rm eff}=0.06$| for all sources, which is around the peak of the |$\chi _{\rm 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 _{\rm GW}$| at low frequencies³² (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_{\rm c}=(m_{1}m_{2})^{3/5}/(m_{1}+m_{2})^{1/5}$| of a BBH merger at redshift z, its contribution to the SGWB |$\delta \Omega _{\rm GW}=(c^{2}\rho _{\rm c}\tilde{V}_{\rm com})^{-1}\nu (\mathrm{d}E_{\rm GW}/\mathrm{d}\nu _{\rm s})$| peaks at an observer-frame frequency |$\nu _{\rm p}\propto 1/[(1+z)m_{\rm c}]$| with |$\delta \Omega _{\rm GW}(\nu _{\rm p})\propto m_{\rm c}/(1+z)$| and follows |$\delta \Omega _{\rm GW}\propto m_{\rm c}^{5/3}\nu ^{2/3}(1+z)^{-1/3}$| in the low-frequency inspiral regime (⁠|$\nu \ll \nu _{\rm 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 for |$\sim 18{{\ \rm 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\ \rm Hz$| (contributing |$\sim 80$| per cent of the total energy spectrum at |$\nu \lesssim 10$| Hz). The Pop III SGWB has a peak of |$\Omega _{\rm GW}\sim 10^{-10}$| at |$\nu =14$| Hz and converges to the standard power-law |$\Omega _{\rm GW}\propto \nu ^{2/3}$| form for BBH inspiral at |$\nu \lesssim 10$| Hz. The energy spectrum is almost flat with |$\Omega _{\rm GW}\sim 5\times 10^{-11}$| at |$\nu \sim 30\!-\!100~\rm Hz$| with similar contributions from field and NSC mergers. The field mergers produce a flat energy spectrum for |$\nu \sim 10\!-\!200$| Hz, while the SGWB from NSC mergers shows a small secondary peak at |$\nu \sim 100$| 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_{\rm SGWB}^{\rm PopIII}=13.5{{\ \rm per\ cent}}$| achieved in the low-frequency regime (⁠|$\nu \lesssim 10$| 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$| 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 _{\rm GW}(f=25\ {\rm Hz})=6.9_{-2.1}^{+3.0}\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 |$\sim 5\!-\!8$| in the frequency range |$\nu \sim 1\!-\!200$| 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 |$\sim 3$| per cent. These features of the Pop III SWGB, in 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_{\rm c}\sim 20\ \rm M_\odot$| and has a sharp cutoff at |$m_{\rm c}\sim 40\ \rm M_\odot$|⁠. The number of mergers also declines rapidly for |$m_{\rm c}\lesssim 10\ \rm M_\odot$|⁠. In between (⁠|$m_{\rm c}\sim 10\!-\!40\ \rm 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_{\rm c}\lesssim 10\ \rm 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_{\rm c}\sim 30$|⁠, 100, and |$300\ \rm 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_{\rm c}\sim 100\ \rm M_{\odot }$| insignificant. However, there is a (narrow) gap at |$m_{\rm c}\sim 200\ \rm 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 _{\rm GW}(\nu _{\rm p})\propto m_{\rm c}/(1+z)$| with the peak location |$\nu _{\rm p}$| inversely proportional to |$m_{\rm c}(1+z)$| to the first order. The strong dependence of |$\delta \Omega _{\rm GW}(\nu _{\rm p})$| on |$m_{\rm c}$| indicates that it is the (low-z) massive NSC mergers with |$m_{\rm c}\gtrsim 60\ \rm M_\odot$| around the second and third peaks in the chirp mass distribution that produce the primary peak of the SGWB at |$\nu =14$| Hz in Fig. 7, while less massive NSC mergers shape the secondary SGWB peak at |$\nu \sim 100$| Hz. On the other hand, due to the redshift dependence, Pop III BBH mergers in galaxy fields (with |$m_{\rm c}\sim 10\!-\!40\ \rm M_{\odot }$|⁠) make their peak contributions to the SGWB in a frequency range |$\nu \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 |$\nu \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\gt 0$|⁠). The mean values of |$m_{\rm c}$|⁠, |$m_{1}$|⁠, and q are summarized in Table 4. In galaxy fields, almost all Pop III BBH mergers come from BHs below |$50\ \rm M_\odot$|⁠. This upper bound is even significantly lower than the lower edge of the Pop III single-star PISN mass gap (⁠|$86\ \rm M_\odot$|⁠). In fact, only |$\sim 0.31\ (0.14)$| per cent of field mergers have |$m_{1}\gt 50\ (100)\ \rm M_\odot$|⁠, among which |$51\ (41){{\ \rm per\ cent}}$| are ejected from NSCs, while |$59.2\ (35.5){{\ \rm per\ cent}}$| of NSC mergers have |$m_{1}\gt 50\ (100)\ \rm M_\odot$|⁠. The mass ratio distribution of field mergers 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 |$\sim 86 \!-\! 242\ \rm 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 |$\sim 130\!-\!242\ \rm M_\odot$| so another narrower gap around |$86\!-\!130\ \rm 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_{\rm c}\sim 30\ \rm M_\odot$| is made of two BHs below the gap. (2) The intermediate-mass peak at |$m_{\rm c}\sim 100\ \rm M_\odot$| involves one BH below the gap and another above the gap. (3) The massive peak at |$m_{\rm c}\sim 300\ \rm 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_{\rm p}$| distributions of the progenitors 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_{\rm p}\equiv a(1-e)\lesssim 0.16$| au, which are mainly produced by close stellar binaries with |$r_{\rm p}\lesssim 10$| au initially. However, NSC mergers mostly come from initial BBHs with |$r_{\rm p}\sim 0.1-10^{3}$| au, and the shape of the |$r_{\rm p}$| distribution for such NSC merger progenitors is similar to that of all BBHs for |$r_{\rm p}\gtrsim 0.1$| au. The corresponding progenitor stellar binaries also have a broad |$r_{\rm p}$| distribution extending to |$r_{\rm 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_{\rm 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 |$~\rm 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_{\rm GW}$|⁠, and (b) survive binary-single encounters and merge by DH in NSCs with a given velocity dispersion |$\sigma _{\star }$|⁠. On the left side of the dotted line for |$t_{\rm GW}=14$| Gyr, we have |$m_{1}\lesssim 50\ \rm M_\odot$| for most BBHs that can merge within a Hubble time in isolation. Such BBHs originate from close (⁠|$r_{\rm p}\lesssim 10$| au) binaries of stars with ZAMS masses below |$242\ \rm M_\odot$|⁠, which lose their envelopes through mass transfer and/or CE phases, so the resulting naked cores can only produce BHs below |$\sim 50\ \rm 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 on the MS also fills the Roche lobe due to significant expansion³³ (Costa et al. 2023). Initially wider (⁠|$r_{\rm 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)\ \rm M_\odot$| from stars with ZAMS masses below (above) |$242 \rm M_\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 _{\star }=30\ \rm 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\ \rm M_\odot$| originate from massive stellar binaries with |$m_{\star ,1}\gtrsim 242\ \rm M_\odot$| and form two groups in the |$\log m_{1}-\log a$| space. The first group with |$m_{1}\sim 130\!-\!320\ \rm 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_{\rm 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\ \rm M_\odot$| (above the PISN mass gap) and |$a\gtrsim 30$| AU comes from initially wider (⁠|$r_{\rm 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 _{\star }=300\ \rm 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{{\ \rm per\ cent}}$|⁠) field mergers occur in small haloes with |$M_{\rm h}\lesssim 5\times 10^{9}\ \rm M_\odot$| at |$z\gtrsim 0.8$|⁠, which host low-mass, metal-poor galaxies with |$M_{\star }\lesssim 1.3\times 10^{6}\ \rm M_\odot$| and |$Z\lesssim 0.22\ \rm Z_\odot$|⁠. The typical (median) host system has |$M_{\rm h}\sim 9.4\times 10^{7}\ \rm M_\odot$|⁠, |$M_{\star }\sim 3.5\times 10^{3}\ \rm M_\odot$|⁠, |$Z\sim 0.005\ \rm 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{{\ \rm per\ cent}}$|⁠) of NSC mergers happen in more massive structures with |$M_{\rm h}\sim 6.7\times 10^{9}-1.1\times 10^{11}\ \rm M_\odot$|⁠, |$M_{\star }\sim 1.3\times 10^{6}-8\times 10^{8}\ \rm M_\odot$|⁠, |$Z\sim 0.017-0.48\ \rm Z_\odot$|⁠, and |$M_{\rm SC}\sim 9\times 10^{4}\!-\!3.4\times 10^{9}\ \rm M_\odot$| at |$z\sim 0.1\!-\!3$|⁠. The typical host system has |$M_{\rm h}\sim 2\times 10^{10}\ \rm M_\odot$|⁠, |$M_{\star }\sim 1.5\times 10^{7}\ \rm M_\odot$|⁠, |$Z\sim 0.075\ \rm Z_\odot$|⁠, |$M_{\rm SC}\sim 4.5\times 10^{5}\ \rm 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}\ \rm 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_{\rm h}\sim 10^{10}\ \rm M_\odot$|⁠, |$M_{\star }\sim 3.6\times 10^{6}\ \rm M_\odot$|⁠, |$Z\sim 0.05\ \rm Z_\odot$|⁠, |$M_{\rm SC}\sim 2.8\times 10^{5}\ \rm M_\odot$|⁠, and |$z\sim 1.1$|⁠. In the opposite case of KRO1_low, we have |$M_{\rm h}\sim 3.7\times 10^{10}\ \rm M_\odot$|⁠, |$M_{\star }\sim 4.9\times 10^{7}\ \rm M_\odot$|⁠, |$Z\sim 0.17\ \rm Z_\odot$|⁠, |$M_{\rm SC}\sim 8\times 10^{5}\ \rm 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 wide IBS model in which close binary interactions are rare, |$\epsilon _{\rm BBH}$| is not very sensitive to the IMF, and the fiducial log-flat IMF shows the highest |$\epsilon _{\rm 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 _{\rm 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 _{\rm BBH}$| and more massive BBHs, as shown in Fig. 14.

For Pop III BBH mergers in galaxy fields, the merger efficiency |$\epsilon _{\rm GW}^{\rm field}$| is sensitive to the NSC parameters since only a small (⁠|$\sim 3\!-\!12{{\ \rm per\ cent}}$|⁠) fraction of Pop III BBHs fall into NSCs. Therefore, we focus on the results of |$\epsilon _{\rm GW}^{\rm field}$| for the fiducial NSC model (obs). We find that |$\epsilon _{\rm GW}^{\rm field}$| always increases with |$\alpha$|⁠, and the trend is stronger with the close IBS model. The evolution of |$\epsilon _{\rm GW}^{\rm field}$| with |$\alpha$| is also stronger compared with the case of |$\epsilon _{\rm BBH}$|⁠, as the fraction |$f_{\rm GW}^{\rm field}$| of BBHs in galaxy fields that merge at |$z\gt 0$| increases with |$\alpha$|⁠. These trends can be explained by the fact that field mergers mostly originate from interacting close (⁠|$r_{\rm p}\lesssim 10$| au) binaries of relatively small stars (⁠|$\lesssim 240\ \rm 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 _{\rm GW}^{\rm 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 _{\rm GW}^{\rm NSC}$| with IMF and IBSE is very similar to that of |$\epsilon _{\rm BBH}$|⁠. Indeed, if the probability for a BBH to merge in a NSC at |$z\gt 0$| is constant, |$\epsilon _{\rm GW}^{\rm NSC}$| should be proportional to |$\epsilon _{\rm 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 _{\rm BBH}$|⁠, |$\epsilon _{\rm GW}^{\rm 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_{\rm GW}^{\rm NSC}$| of BBHs in NSCs that merge at |$z\gt 0$| and the fraction |$f_{\rm infall}$| of BBHs that fall into NSCs (see the last two columns of Table 2) with IMF and IBS. As expected, |$f_{\rm infall}$| and thus |$\epsilon _{\rm GW}^{\rm 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_{\rm 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_{\rm 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_{\rm 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 LOG1_obs (i.e. the fiducial model) and LOG5_obs 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 _{\rm GW}$| (Fig. 13), the MRDs of field and NSC mergers are reduced and increased when there are more wide binaries in the initial condition (LOG5_obs), 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 _{\rm GW}^{\rm NSC}$| and |$\epsilon _{\rm BBH}$|⁠. For field mergers, the MRD evolution is 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_{\rm 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 LOG1_obs with the close IBS model, the contribution from field (NSC) mergers is lower (higher) for LOG5_obs 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 LOG5_obs, 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\ \rm Hz$| deviating from the canonical shape |$\Omega _{\rm 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 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\ \rm 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\ \rm 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 LOG1_obs, TOP_obs, and KRO1_obs 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 _{\star ,\rm III}\simeq 4.6$|⁠, 5.6 and |$5.7\times 10^{4}\ \rm 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 _{\star ,\rm 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 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}\ \rm 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 _{\rm 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_{\rm c}\lesssim 60\ \rm 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_{\rm c}\gtrsim 60\ \rm M_\odot$| 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_{\star }\gt 10^{6}\ \rm 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{{\ \rm 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 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 |$\nu \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{\mathcal {N}}_{\rm ET,SNR\gt 8}$| of Pop III BBH mergers by ET for all the 18 models (Table 2) given the detection probability of each source with |$\rm SNR\gt 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\lt 20$| will be detectable by ET thanks to its high sensitivity across a broad frequency range, so the all-sky merger rate |$\dot{\mathcal {N}}(z\lt 20)$| within |$z\lt 20$| is a good approximation of the ET detection rate (for |$\rm SNR\gt 8$|⁠). Therefore, we derive the fraction of |$\dot{\mathcal {N}}(z\lt 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{\mathcal {N}}_{\rm ET,SNR\gt 8}\sim 50\!-\!1370\ \rm yr^{-1}$|⁠, and we have |$\dot{\mathcal {N}}_{\rm ET,SNR\gt 8}\sim 30\!-\!1230\ \rm yr^{-1}$| and |$\dot{\mathcal {N}}_{\rm ET,SNR\gt 8}\sim 6\!-\!230\ \rm 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\lt 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_{\rm NSC}\sim 65\!-\!95{{\ \rm 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\ \rm M_\odot$|⁠) BHs, especially those above |$130\ \rm M_\odot$|⁠. In light of this, we further look into the ET detection rates |$\dot{\mathcal {N}}_{\rm ET,IMBH,SNR\gt 8}$| of Pop III mergers with at least one IMBH (⁠|$m_{1}\gt 100\ \rm M_\odot$|⁠). We find that such mergers are extremely rare in galaxy fields, with ET detection rates below |$0.1\ \rm yr^{-1}$| in all the 18 models considered here. However, we have |$\dot{\mathcal {N}}_{\rm ET,IMBH,SNR\gt 8}\sim 0.5\!-\!200\ \rm 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 sensitive to all model parameters on Pop III IMF, IBS, and occupation fraction of high-z NSCs, implying that BBH mergers with IMBHs³⁴ of |$\sim 100\!-\!1000\ \rm 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³⁵: |$\dot{\mathcal {N}}_{\rm CEO,SNR\gt 8y}\sim 30\!-\!1270\ \rm yr^{-1}$| for field mergers, |$\dot{\mathcal {N}}_{\rm CEO,SNR\gt 8}\sim 6\!-\!210\ \rm yr^{-1}$| for NSC mergers, and |$\dot{\mathcal {N}}_{\rm CEO,SNR\gt 8}\sim 0.5\!-\!120\ \rm 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\ \rm 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){{\ \rm per\ cent}}$| for CEO (in TOP5_full), while it remains |$\gtrsim 95{{\ \rm per\ cent}}$| for ET.

Finally, we estimate the detection rate of Pop III BBH mergers by the LVK network during the O4 run with |$\rm SNR\gt 8$| as |$\dot{\mathcal {N}}_{\rm LVK,SNR\gt 8}\sim 0.9\!-\!9\ \rm yr^{-1}$|⁠. This small detection rate is consistent with the small local MRD |$\dot{n}(z=0)\sim 0.034\!-\!0.17\ \rm yr^{-1}\ Gpc^{-3}$| that is much lower than the observed local MRD |$\dot{n}_{\rm obs}(z=0)=19.3_{-9}^{+15.1}\ \rm yr^{-1}\ 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{\mathcal {N}}_{\rm LVK,IMBH,SNR\gt 8}\sim 0.04\!-\!1.5\ \rm yr^{-1}$|⁠, and the expected number of events detected in 18 months³⁶ of observations is only |$\sim 0.07\!-\!2.2$|⁠. Considering the typical horizon redshift |$z_{\rm horizon}\sim 1$| of the LVK network during O3 for BBH mergers with total masses |$\sim 100\!-\!300\ \rm M_\odot$| (see fig. 8 of De Luca et al. 2021), we derive the MRD of Pop III mergers with IMBHs for |$z\lt 1$| as |$\dot{n}_{\rm IMBH}(z\lt 1)\sim 0.003\!-\!0.15\ \rm yr^{-1}\ 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 mass gap |$50\!-\!130\ \rm M_\odot$| 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}_{\rm gap}(z\lt 1)\sim 0.0013\!-\!0.034\ \rm yr^{-1}\ Gpc^{-3}$| for |$z\lt 1$|⁠, also consistent with the rate |$0.08_{-0.07}^{+0.19}\ \rm 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:

• Although only a small fraction (⁠|$3\!-\!12{{\ \rm per\ cent}}$|⁠) of Pop III BBHs fall into NSCs, a significant fraction (⁠|$\sim 45\!-\!64{{\ \rm per\ cent}}$|⁠) of the BBHs in NSCs merge at |$z\gt 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_{\rm NSC}\sim 8\!-\!95{{\ \rm per\ cent}}$| of Pop III BBH mergers across cosmic history, and is generally more important at lower redshifts, accounting for |$\mathcal {F}_{\rm NSC}\sim 34\!-\!98$| per cent (⁠|$0.7\!-\!88$| per cent) of the all-sky merger rates of Pop III BBH mergers within |$z\lt 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.

• 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_{\rm NSC}\sim 8\!-\!51$| per cent of all Pop III mergers across cosmic history and explains |$\mathcal {F}_{\rm NSC}\sim 0.7\!-\!22$| per cent of the all-sky merger rate within |$z\lt 20$|⁠. However, when the IBS is dominated by wide binaries, we have |$f_{\rm NSC}\sim 65\!-\!95$| per cent and |$\mathcal {F}_{\rm NSC}\sim 12\!-\!88$| per cent for |$z\lt 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).

• Pop III BBH mergers in NSCs are more massive than those from IBSE. The latter mostly involve BHs below |$\sim 50\ \rm M_\odot$| from the collapse of naked cores of stellar progenitors that have relatively low initial masses (⁠|$\lesssim 240\ \rm M_\odot$|⁠) 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\ \rm M_\odot$|⁠) 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{{\ \rm per\ cent}}$|⁠) of Pop III BBH mergers in NSCs involve at least one intermediate mass BH (IMBH) above |$100\ \rm M_\odot$|⁠, while such mergers are extremely rare from IBSE.

• 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 |$\dot{\mathcal {N}}_{\rm ET,SNR\gt 8}\sim 50\!-\!1370\ \rm yr^{-1}$|⁠, and we have |$\dot{\mathcal {N}}_{\rm ET,SNR\gt 8}\sim 30\!-\!1230\ \rm yr^{-1}$| and |$\dot{\mathcal {N}}_{\rm ET,SNR\gt 8}\sim 6\!-\!230\ \rm 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}\gt 100\ \rm M_\odot$|⁠) per year from the NSC-DH channel (mostly at |$z\lesssim 10$|⁠). The detection rate of such mergers from IBSE is below |$0.1\ \rm 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.

• 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\ \rm M_\odot$|⁠) 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{{\ \rm 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 (⁠|$\dot{\mathcal {N}}_{\rm LVK,IMBH,SNR\gt 8}\sim 0.04\!-\!1.5\ \rm 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).

• The SGWB produced by Pop III BBH mergers has a peak value of |$\Omega _{\rm GW}\sim 10^{-11}\!-\!8\times 10^{-11}$| around observer-frame frequencies |$\nu _{\rm 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.³⁷

  The SGWB from Pop III BBH mergers is equal to |$\sim 2\!-\!32{{\ \rm per\ cent}}$| of the total SGWB of compact object mergers inferred from observations (Abbott et al. 2023a) at |$\nu \lesssim 10\ \rm 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 |$\nu \lesssim 10\ \rm Hz$|⁠. Considering the cumulative Pop III stellar mass density |$\rho _{\star ,\rm III}\sim 5\times 10^4\ \rm M_\odot \ Mpc^{-3}$| found in our simulations, such enhancement of Pop III star formation is still allowed by the constraints |$\rho _{\star ,\rm III}\lesssim 10^{5\!-\!6}\ \rm 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 (⁠|$\gt 100\ \rm M_\odot$|⁠) and BHs in the (standard) PISN mass gap (⁠|$50\!-\!130\ \rm M_\odot$|⁠) efficiently. This is consistent with previous studies showing that such massive BBH mergers form efficiently via dynamical processes in SCs (Rodriguez 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\ \rm 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\ \rm M_\odot$|⁠) can still form via IBSE.

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 |$\sim 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_{\rm h}\sim 10^{5\!-\!6}\ \rm 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 miss-aligned 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 (⁠|$\sim 15^{\circ }$|⁠) 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\sim 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 50OO2206). 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.

Footnotes

REFERENCES