© The Authors 2025

1 Introduction

The ever-growing list of gravitational wave signals detected by the LIGO-Virgo-KAGRA Collaboration (LVK)¹ has now definitively proven the existence of double compact objects and, in particular, of binary black holes (BBHs), with initial orbital separations small enough for the two objects to merge in less than a Hubble time (see the third Gravitational Wave Transient Catalog, GWTC-3; Abbott et al. 2023b). However, the astrophysical processes leading to the formation of these binaries remain challenging to constrain. Recent analyses of the physical properties of these detected BBHs show growing evidence that several formation channels are likely to be involved (see e.g. Zevin et al. 2021; Arca Sedda et al. 2023; Ray et al. 2024).

As massive stars are expected to end their lives as black holes and are observed to exist primarily in binaries or multiples (see e.g. Sana et al. 2012; Moe & Stefano 2017), a natural scenario to explain the formation of these BBHs seems to be the conjoint evolution of two massive stars. This formation channel has been extensively studied over the past decades (e.g. Bethe & Brown 1998; Belczynski et al. 2002; Dominik et al. 2012; Belczynski et al. 2016; Eldridge & Stanway 2016; Stevenson et al. 2017; Giacobbo & Mapelli 2018; Neijssel et al. 2019; Santoliquido et al. 2021; Tauris & van den Heuvel 2023) and significant progress has been made in our understanding of binary stellar evolution. Most studies rely on population synthesis techniques that make it possible to numerically build large populations of double compact objects. Using models of metallicity-dependent star formation rate to describe the evolution of the Universe, to distribute their mergers over cosmic time. However, a number of uncertainties still persist in this framework, both in the description of the Universe and in the prescriptions used to evolve binary stars (see e.g. Broekgaarden et al. 2022; Santoliquido et al. 2022; Srinivasan et al. 2023, for studies on the relative importance of each of these two aspects).

A second possible pathway that could lead to the formation of two BHs orbiting each other with very close separation is to consider dynamical interactions in dense environments. In star clusters, where the presence of stellar BHs has been confirmed both through X-ray and radio observations (see e.g. Maccarone et al. 2007; Strader et al. 2012; Chomiuk et al. 2013; Miller-Jones et al. 2015; Shishkovsky et al. 2018) and radial velocity measurements of detached star-BH binaries (see e.g. Giesers et al. 2018, 2019), it is expected that mass segregation and three-body interactions would enhance the pairing of BHs in dynamically-hard binaries (see e.g. Morscher et al. 2015). These dense stellar structures, such as young clusters, globular clusters, or even nuclear clusters, constitute environments favourable to the formation of GW sources in the form of tight BBHs (see e.g. Sigurdsson & Hernquist 1993; Portegies Zwart & McMillan 1999; Rodriguez et al. 2015, 2016; Banerjee 2017; Di Carlo et al. 2020). As both direct and Monte Carlo simulations of dynamics in star clusters can be numerically quite expensive, most studies consider a rather modest number of simulated star clusters in a predetermined parameter space. The recent development of semi-analytical codes for the rapid evolution of star clusters, such as FASTCLUSTER (Mapelli et al. 2021) or RAPSTER (Kritos et al. 2024), now makes it possible to operate BBH population synthesis for much larger numbers of star clusters with a wide range of physical properties. The trends observed in the results of their evolution can then be used to predict cosmologically statistically significant quantities, such as the evolution of the BBH merger rate in star clusters over cosmic time (Rodriguez & Loeb 2018; Kremer et al. 2020; Mapelli et al. 2022; Ye & Fishbach 2024).

Although stellar populations are in most cases described with the combination of well-known models, such as a star formation rate (SFR, Madau & Dickinson 2014) and an initial mass function (Kroupa 2001), it is much more complex to describe the evolution of star cluster physical properties over cosmic time (see e.g. Krumholz et al. 2019, for a review). In particular, several quantities used to describe the formation of star clusters have been found to vary depending on the galactic environment, including the cluster formation efficiency (Ginsburg & Kruijssen 2018), the cluster initial mass function with a potential high-mass truncation (Wainer et al. 2022), and the mass-radius relation (Brown & Gnedin 2021). In addition, the evolution of a single star cluster over time is virtually impossible to predict using simple models alone, and shows major changes between its initial state and what we can observe as its present-day aspect. Both the formation and the evolution of star clusters are also highly dependent on a number of environmental effects (see e.g. Rossi et al. 2016; Suin et al. 2022). This implies that a realistic model of the formation and evolution of star clusters over cosmic time should take into account the evolution of the environmental conditions specific to each galaxy, and this with a sufficiently large number of different galaxies to be able to represent their diversity in the Universe. Here again, a number of uncertainties remain as to the physical properties of star clusters used to describe a complete cosmological population, and as to the diversity of results predicted from cluster dynamics for different initial conditions.

A second approach to study the cosmological population of dynamically formed BBHs is to use a large volume cosmological simulation as the natural environment in which stars and clusters can form. Given the high computational cost of this type of simulation, a compromise must always be found between the spatial resolution of the simulation, its volume, and the redshift to which it evolves. To date, the majority of numerical simulations with a sub-parsec scale resolution necessary to resolve the formation of massive star clusters have only been made in small volumes and/or at large redshifts (see e.g. Boley et al. 2009; Kimm et al. 2016; Lahén et al. 2019; Calura et al. 2022, 2025, the latter resolving the formation of star clusters in zoom-in simulations of cosmological volume of 5cMpc³ from z = 100 down to z = 10.5). To predict the formation rate and physical properties of the star clusters that would populate any given simulated galaxy up to the present-day, semi-analytical models are then required. Studies of this kind have already been carried out, with favourable results when compared with the population of globular clusters (GCs) observed in local galaxies (e.g. Li & Gnedin 2014; Pfeffer et al. 2018; Grudić et al. 2021).

This paper follows the Great Balls of FIRE series (Grudić et al. 2022; Rodriguez et al. 2023; Bruel et al. 2024, hereafter Paper I, Paper II, and Paper III respectively). Here, we elaborate upon these analyses of the formation and evolution of massive star clusters inside cosmological zoom-in simulations of galaxies to develop models of cluster formation that can be applied to a larger cosmological simulation containing a wide variety of simulated galaxies. We consider here only the production of BBHs in massive star clusters, defined as the star clusters with initial masses larger than 6 × 10⁴ M_⊙ (value selected to maintain consistency with our previous studies Paper II and Paper III, see in Sect. 2.2 below for more details). Making use of the large catalogue of massive clusters already integrated forward in time in Paper II; Paper III with the Cluster Monte Carlo code (CMC Pattabiraman et al. 2013; Rodriguez et al. 2022), we are able to make predictions about the contribution of different galaxies and environments to the dynamical population of merging BBHs. Throughout this study and in all that follows, we only consider BBHs that have already merged by z = 0 and every reference to the term ‘BBHs’ should be interpreted as ‘merging BBHs’.

After describing the cosmological volume simulation that we use as a large sample of realistic galaxies and the three models of cluster formation that we consider (Section 2), we present our predictions on massive star clusters in galaxies at a cosmological scale and on the BBH mergers that these clusters produce (Section 3). We discuss implications of our findings and outline future extensions of this work (Section 4) and conclude on the significance of our results (Section 5).

2 Method

In this work, we aim to estimate the physical properties and merger rate of dynamically formed BBHs at a cosmological scale. We use a large sample of galaxies identified in a cosmological volume simulation (Sect. 2.1) and propose three different semi-analytic models for massive star cluster formation (Sect. 2.2) to model the population of star clusters in each galaxy. We build upon the large collection of massive star clusters (1500) already integrated with the code CMC in Paper II; Paper III to develop a grid-matching algorithm that makes it possible to estimate the BBH mergers formed dynamically in any population of massive star clusters (Sect. 2.3).

2.1 Realistic star formation histories from individual galaxies in FIREbox

FIREbox (Feldmann et al. 2023) is a large scale cosmological simulation created as part of the FIRE project (Hopkins et al. 2014). It represents a cosmological volume of (22.1cMpc)³ evolved from redshift 120 to the present-day, with 1200 snapshots spaced in time. It uses the same FIRE-2 physics model as the zoom-in simulations used in Paper III. FIRE-2 models radiative cooling and heating across the range 10–10¹⁰ K, including free-free, photo-ionisation and recombination, Compton, photo-electric and dust collisional, cosmic ray, molecular, metal-line, and fine-structure processes. Photo-ionisation and heating from a redshift-dependent, spatially uniform ultraviolet background (Faucher-Giguère et al. 2009) is also taken into account. Star formation occurs in self-gravitating, self-shielding, and Jeans unstable dense molecular gas (n > 300 cm⁻³). All feedback event rates, including energy, momentum, mass and metal injection from type Ia and type II supernovae, and stellar winds luminosities, are tabulated from stellar evolution models (STARBUST99; Leitherer et al. 1999) assuming a Kroupa (2001) initial stellar mass function (IMF). A sub-grid model is applied for the turbulent diffusion of metals in gas (Escala et al. 2017).

We chose FIREbox specifically in the present work to maintain a certain consistency with our previous analysis, but also because it offers a spatial dynamic range of ~10⁶, which is about an order of magnitude higher than most large-volume simulations. This corresponds to a spatial resolution of ~20 pc, well suited to explore the internal structure of galaxies, but still too large to resolve individual star clusters. The baryonic mass resolution is m_b ~ 6.3 × 10⁴ M_⊙. Although this simulation is considered a large volume and contains thousands of galaxies, the total volume of the simulated cube is not yet sufficient to be representative of our Universe. Indeed, most analyses suggest that the typical homogeneity scale of our Universe is of the order of ℛ_H ≳ 100 h⁻¹Mpc (see e.g. Scrimgeour et al. 2012; Ntelis et al. 2018). This is an obvious limitation of the analysis we carry out here and should be borne in mind in all that follows. The effect of the cosmic variance inherent in this restricted volume will be discussed in more detail below.

We used the halo structures already identified inside FIREbox with the AMIGA Halo Finder (AHF, Knollmann & Knebe 2009). In particular, in the last snapshot (corresponding to redshift 0), we identified 979 halos with present-day virial masses M_vir ≥ 10¹⁰ M_⊙ (which corresponds to a lower limit on galaxy stellar mass of roughly ℳ_⋆ ≥ 10⁸ M_⊙). This provides us with a large sample of realistic galaxies that can be used to study the formation of star clusters inside each of them across cosmic time. Given the empirical near-constant fraction of the total mass of GCs in a galaxy divided by its total mass (dark and baryonic matter), no massive star clusters are expected to form in less massive dwarf galaxies (Harris et al. 2017; Chen & Gnedin 2023; Bruel et al. 2024). Using high-resolution optical imaging of quiescent galaxies in the Virgo cluster, Sánchez-Janssen et al. (2019) found that the GC occupation fraction rapidly falls off for galaxy stellar masses below 10⁹ M_⊙. Thus, our galaxy stellar mass cut-off ℳ_⋆ ≥ 10⁸ M_⊙ was expected to have only a minor impact on the results presented here on massive cluster formation and dynamically formed BBHs. For each galaxy we extracted the following present-day parameters: the virial radius, R_vir, the halo mass, ℳ_h, the stellar radius, R_⋆,90, computed as the spherical radius that encloses 90% of the stellar mass within 20 kpc, see (as in Wetzel et al. 2023), and the stellar mass, M_⋆,90, computed as the mass enclosed in that radius. Throughout the rest of this paper, each reference to the stellar mass ℳ_⋆ of a galaxy identified in FIREbox corresponds to this calculation of its present-day stellar mass M_⋆,90.

We display the halo mass function (HMF) obtained from these halos in the left-hand panel of Figure 1. The HMF presented in Tinker et al. (2008), obtained from a large set of collisionless dark-matter only cosmological simulations with flat ΛCDM cosmology, is also shown as a comparison. There appears to be a mismatch between the HMF realised in FIREbox and this ‘universal’ form of the HMF. Feldmann et al. (2023) explain this discrepancy by arguing that the HMF that exists and evolves throughout the cosmological volume of the FIREbox simulation differs from a true local HMF because of a cosmic variance that particularly affects smaller volumes. This variance is due to the finite size of the box and its limited numerical resolution.

As this study is aimed at comparing the predictions of this numerical approach with the LVK detections of BBH mergers, this sample of simulated galaxies needed to be compare with the actual distribution of galaxies in the Universe. A reasonable match at redshift 0 would not be sufficient to guarantee that the properties of our galaxies are realistic at all epochs, but it is already a first step towards constraining the significance of our subsequent analysis. To account for the cosmic variance described above, we re-weighted all the identified halos to ensure that their HMF matches the one of Tinker et al. (2008) (as indicated by the arrows in the left-hand panel of Figure 1). In practice, we computed the present-day HMF in FIREbox using the logarithm of the halo masses at redshift z = 0 and 0.2 wide bins between 10 and 14, and then calculate the fraction HMF_FIREbox/HMF_Tinker for each bin (see also Feldmann et al. 2023, for a more complete and detailed version of such a reweighting process). These values were saved as the weights for all of the 979 galaxies. This approach is equivalent to assuming that the halos and galaxies in FIREbox may very well have a realistic evolution, but are not in the right number to represent a global realistic population and their abundance (or number density) needs to be modified accordingly.

The galaxy stellar mass function (GSMF) obtained from these galaxies is shown in the right-hand panel of Figure 1. As a comparison, the best-fit GSMF presented in Tomczak et al. (2014) at redshift z ~ 0.1 is plotted in the same Figure. Furthermore, we also provide an estimate of the GSMF at z = 0 built using the non-parametric modelling of Leja et al. (2020). With the exception of an excess of both the ≲10¹⁰ M_⊙ and the highest mass galaxies, the GSMF in FIREbox is in qualitative agreement with the latter estimate.

For each identified galaxy, we extracted the formation time and metallicity (computed here as the mass fraction of all metals tracked in the simulation) for all the star particles that lie (at the present-day) inside the sphere of stellar radius, R_⋆,90. In the right-hand panel of Figure 2, we show the evolution of the average metallicity of newly formed stars as a function of time. The mean metallicity follows an overall similar trend when considering the entire simulation and when grouping galaxies by their present-day stellar masses. We refer to Bassini et al. (2024) for a more detailed analysis of the mass-metallicity relation in the FIREbox simulation. Consistently with observations of the mass-metallicity relation (see e.g. Mannucci et al. 2010; Nakajima et al. 2023), the most massive galaxies are found to be also the most-metal rich.

For all star particles associated with identified galaxies, we used their formation time, metallicity, and mass at z = 0 to calculate their initial masses. Thus, we were able to gain access to the mass-weighted formation history of this population of stars. This calculation yields an ‘archaeological’ star formation history for all the FIREbox galaxies. The resulting cosmic SFRD is shown in the left-hand panel of Figure 2. It is in good agreement with the observed cosmic SFRD (Madau & Dickinson 2014) up to redshifts z ≥ 1, but over-predicts the SFRD at lower redshifts (see also Feldmann et al. 2023). Even with the re-weighting operation to account for the cosmic variance in this type of simulation, there is still an excess of ~0.6 dex in the total SFRD at present-day (z = 0). This excess of star formation is mostly driven by a few very massive galaxies (17 galaxies with stellar masses in the range [10¹¹, 10¹²] M_⊙).

It is important to note that in this large scale cosmological simulation, no prescriptions for supermassive BHs and AGN feedback are included. As a result, the fraction of quiescent massive galaxies is significantly under-predicted in FIREbox and the contribution of these massive galaxies to the cosmic SFRD at low redshifts is over-estimated (Feldmann et al. 2023). Despite their overall negative impact on star formation and galaxy growth, recent simulations have shown that quasar winds could potentially induce local positive feedback on both star formation (Mercedes-Feliz et al. 2023) and star clusters (Mercedes-Feliz et al. 2024).

Since the overall excess of star formation in FIREbox occurs mainly in very massive galaxies and at low redshifts, it involves the formation of metal-rich stars and star clusters. These stellar structures can be expected to make only a small contribution to the total population of merging BBHs (Di Carlo et al. 2020; Bruel et al. 2024). Indeed, massive stars formed in metal-rich environments experience higher mass loss through stellar winds (see e.g. Vink et al. 2001), which means that a smaller fraction of them will have a final mass high enough to collapse into a BH. Furthermore, for massive stars in binary systems (whether isolated or inside a star cluster), mass loss produces a widening of their orbits, which results in a higher fraction of BBHs that will merge in more than a Hubble time. In practice, we found that ~98% of all the star particles formed at redshifts z ≤ 1 in galaxies with present-day stellar masses ℳ_⋆ ≥ 10¹¹ M_⊙ have values of metallicity higher than Z_⊙ (see the evolution of the average metallicity of newly formed stars for different galaxy masses in the right-hand panel of Figure 2). Such high metallicity values tend to prevent the production of merging BBHs, even in the most massive star clusters (see Figure 3 further on). For this reason, it appears reasonable to assume that the lack of quiescent galaxies in FIREbox will not have a major impact on further analysis. We present a more detailed analysis of this question in Sect. 4.1.

To account for the quenching of star formation due to additional physical processes not implemented in FIREbox, Feldmann et al. (2023) propose to re-weight all galaxies by a quenching factor 1 − f_Q, using f_Q = 20, 45, and 90% for galaxies with stellar masses in the ranges [10⁹, 10¹⁰] M_⊙, [10¹⁰, 10¹¹] M_⊙, and ≥10¹¹ M_⊙, respectively (with values taken from Behroozi et al. 2019). This operation brings the predictions for the SFRD in the simulation at the lowest redshifts in much better agreement with observations (Feldmann et al. 2023). However, as the evolution of BBHs through the emission of GWs can be an extremely slow process, a number of their mergers observed in the local Universe can very well originate from stellar binaries or star clusters formed at much higher redshifts (see e.g. Fishbach & Kalogera 2021). The quenched fractions are used to describe the observed properties of local galaxies and they cannot be applied at all redshifts, which makes this method inapplicable to predictions of BBH mergers. For this reason, we chose not to model the quenching of star formation in massive galaxies with such considerations and, instead, to quantify the impact of this discrepancy in our analyses (see the discussion in Sect. 4.1).

Fig. 1 Left: halo mass function (HMF) of the 979 galaxies identified inside FIREbox within the mass range 10 ≤ log(Mvir/M⊙) ≤ 14 (red). Vertical dashed lines show the range of the virial masses used for halo identification. Grey arrows indicate the re-weighting operated to match the HMF in FIREbox with a ‘universal’ HMF extracted from dark matter cosmological simulation (in black line, from Tinker et al. 2008). Right: galaxy stellar mass function (GSMF) for this same set of galaxies identified inside FIREbox (red), and with each galaxy re-weighted with the operation presented in the left-hand panel (orange). Two examples of GSMF built on observations of galaxies in the local Universe are also shown. The black solid line is the best-fit GSMF at z ~ 0.1 from Tomczak et al. (2014) and the dashed line is the median GSMF generated using the continuity model parameters of Leja et al. (2020) at z = 0. The shaded grey area corresponds to the 16–84% uncertainties associated with this model.

Fig. 2 Left: star formation rate density (SFRD) for the 979 galaxies identified inside FIREbox(black line) and for galaxies grouped in different bins of stellar mass (coloured lines). The number of galaxies in each mass bin is indicated in brackets. As a comparison, the grey dashed line shows the functional form of the cosmic SFRD proposed by Madau & Dickinson (2014). Right: mass-weighted average of the metallicity of newly formed stars in the galaxies identified in FIREbox. The black line indicates the average metallicity is in the entire volume, and the different colours indicate galaxies grouped in different bins of stellar mass. In both panels, the contribution of each galaxy is re-weighted as described in Sect. 2.1 (and as pointed out in the left-hand panel of Figure 1). By way of comparison with the different Figures shown in the FIREbox paper (Feldmann et al. 2023), we do not apply a bootstrapping method here and we therefore do not indicate error bars for these estimates. We refer to that particular paper for more details on these uncertainties.

2.2 Models of cluster formation in a large scale cosmological simulation

Once we have identified the galaxies within the FIREbox simulation and extracted their physical properties and their evolution over time, we sought to describe the star clusters that would populate them. As individual star clusters are not resolved in FIREbox, a cluster formation model is needed to get an estimate of the populations of massive star clusters that would exist within the galaxies of FIREbox. In the following, we present three different sampling procedures used throughout this study for massive cluster formation in various galaxies located in a larger volume.

As FIREbox uses a different value of gas density threshold to trigger star formation (n ≥ 300 cm⁻³ vs. n ≥ 1000 cm⁻³ in the usual FIRE-2 simulations) and has a lower spatial resolution than the zoom-in simulations, we decided to not directly use the same cloud-to-cluster formation model to sample star cluster in FIREbox. In order to take advantage of the unique dynamic range of FIREbox and to maintain the same FIRE-2 physics models as in our previous studies (Paper I; Paper II; Paper III), we develop new models of massive star cluster formation that can be applied to this cosmological simulation.

We created three different simulation-based empirical models of massive star cluster formation based on the star clusters already sampled in our set of six zoom-in simulations (m11q, m11i, m11e, m11h, m11d, and m12i, with respective stellar masses of ℳ_⋆ = 9.2 × 10⁸, 6.1 × 10⁸, 1.4 × 10⁹, 3.6 × 10⁹, 3.9 × 10⁹, and 6.7 × 10¹⁰ M_⊙ respectively, El-Badry et al. 2017, Paper I; Paper III). Here, we enriched this sample of galaxies with three additional simulations from the FIRE-2 project: m12r, m12c, and m12f (Samuel et al. 2019; Garrison-Kimmel et al. 2019). In particular, m12r is a medium-mass galaxy with ℳ_⋆ ≃ 1.7 × 10¹⁰ M_⊙, while m12c and m12f are both MW-like galaxies with ℳ_⋆ ≃ 5.8 × 10¹⁰ M_⊙ and ℳ_⋆ ≃ 7.9 × 10¹⁰ M_⊙, respectively. In all those simulated galaxies, we used the CloudPhinder algorithm (Guszejnov et al. 2019) to locate the giant molecular clouds (GMCs) in each snapshot and applied the cluster formation framework described in Grudić et al. (2021) to predict the population of star clusters that each GMC would produce.

The three models of massive star cluster formation proposed here each have a different level of complexity and are based on different assumptions. This makes it possible to highlight certain points of comparison when studying their predictions. In the following, we describe each of these models in detail.

2.2.1 ‘Gamma’ model: Constant cluster formation efficiency

In this first approach, we made the straightforward assumption that in a given galaxy, the formation rate of all star cluster is at any time a constant fraction of the star formation rate: $${\dot{M}}_{\text{clusters}}[t]\equiv \mathrm{\Gamma }\times \mathrm{\Psi }[t],$$(1)

where ${\dot{M}}_{\text{clusters}}[t]$ is the formation rate of all star clusters at time t in this specific galaxy, Γ is its cluster formation efficiency assumed constant over time, and Ψ[t] is its star formation rate at time t.

For our sample of zoom-in simulations of galaxies, we observed that the most massive galaxies are more efficient at forming star clusters (see Paper III, for more details). To account for this trend, and for the possible variability of cluster formation in different galaxies with similar stellar masses, we proceed as follows:

• based on the results obtained with individual simulations of galaxies, we fit the average cluster formation efficiency as a function of the galaxy stellar mass such that Γ ≃ 0.043 + 0.019 × log(ℳ_⋆/10⁹M_⊙) with a standard deviation σ ≃ 2.9 × 10⁻³. For a given FIREbox galaxy, we can use its stellar mass, ℳ_⋆, to randomly draw one value of Γ from the fit described above. Equation (1) provides us with an estimate of the cluster formation rate in each galaxy.

• similarly, we observe in Paper III and with the additional zoom-in cosmological simulations considered here that the distribution of cluster initial masses typically gets steeper with decreasing galaxy stellar mass. Consequently, we fit the average power-law index of the cluster initial mass function such that α ≃ −3.01 + 0.43 × log(ℳ_⋆/10⁹ M_⊙) with a standard deviation σ ≃ 0.11. We used exactly the same approach in all three models to describe the dependence of the cluster mass function on the galaxy present-day stellar.

• finally, we used the combination of the cluster formation rate and the power-law slope of the distribution of cluster initial masses to sample individual cluster masses in all the snapshots of FIREbox. In this sampling process, the minimum cluster mass possible is taken to be 10³ M_⊙ (same value used in Paper I and Paper III, and in the additional zoom-in simulations considered here). Since in this study, we are only interested in the production of BBHs in the most massive star clusters, while maintaining consistency with the cluster catalogues presented in Paper I; Paper II; Paper III, we applied a mass threshold M_cl ≥ 6 × 10⁴ M_⊙ to the collections of star clusters sampled in each galaxy. This value roughly corresponds to the an initial number of particles (~10⁵) which is the minimum required to ensure that the cluster relaxation timescale is always significantly longer than its dynamical timescale (a necessary condition in the evolution of star clusters with the code CMC that we used to integrate the clusters forward in time in Paper II and Paper III, see Sect. 2.3 below for a brief description of this code).

Fig. 3 Upper left panel: 2D projection on the plane (Mcl, rhm) of the 1500 clusters integrated with CMC. Clusters in which no merging BBH are formed are indicated as grey dots, while the coloured dots represent the number of merging BBHs formed in each cluster. The diagonal dashed grey lines represent lines of constant density ${\rho }_{\mathrm{c}\mathrm{l}}\equiv M_{\mathrm{c}\mathrm{l}}/(4/3)\pi r_{\mathrm{h}\mathrm{m}}^3$. Upper right panel: 2D projection on the plane (Mcl, Zcl). Colours indicate the number of merging BBHs with the same norm as in the left-hand panel. Lower panel: 3D grid in parameter space (Mcl, rhm, Zcl).

2.2.2 ‘EB18’ model: Gas surface density

In this second step, we present a sampling procedure similar to that described in El-Badry et al. (2018). In this model, the formation efficiency of massive clusters is an increasing function of the gas surface density, which is precisely one of the features of the cluster formation model of Grudić et al. (2021) that has been used in our set of zoom-in simulations. Both theoretical models of star formation (see e.g. Kruijssen 2012) and observation of star clusters in nearby galaxies (see e.g. Portegies Zwart et al. 2010) support this idea, with massive young star clusters being observed in high density environments. On the other hand, less dense environments have longer free-fall times and are less efficient at forming stars. This implies that these regions remain gas-rich and that the expulsion of this gas is more likely to unbound any cluster in formation. Using high resolution simulations of collapsing clouds of gas (also evolved with the FIRE-2 model), Grudić et al. (2018) found that both the fraction of gas converted to stars, ε, and the fraction of stars formed in bound clusters, f_bound, can be expressed as functions of the physical properties of the host GMC, and primarily its surface density ${\mathrm{\Sigma }}_{\mathrm{G}\mathrm{M}\mathrm{C}}=M_{\mathrm{G}\mathrm{M}\mathrm{C}}/R_{\mathrm{G}\mathrm{M}\mathrm{C}}^2$. El-Badry et al. (2018) propose a functional form of this cluster formation efficiency as: $${\mathrm{\Gamma }}_{\text{massive}}[t]\equiv \frac{{\dot{M}}_{\text{massive}}[t]}{\mathrm{\Psi }[t]}=\frac{{\alpha }_{\mathrm{\Gamma }}}{1+{({\mathrm{\Sigma }}_{\mathrm{G}\mathrm{M}\mathrm{C}\mathrm{s}}[t]/{\mathrm{\Sigma }}_{\mathrm{c}\mathrm{r}\mathrm{i}\mathrm{t}})}^{-{\beta }_{\mathrm{\Gamma }}}},$$(2)

where ${\dot{M}}_{\text{massive}}[t]$ is the formation rate of massive clusters at time t, Ψ[t] is the SFR at time t, and α_Γ, β_Γ and Σ_crit are constants that we estimate here from the massive star clusters already sampled in our zoom-in simulations of individual galaxies.

Insofar as our analyses indicated that the physical properties of the GMCs present in FIREbox did not correspond exactly to what was expected from the zoom-in simulations, we decided not to take them into account here and to estimate the GMC surface density Σ_GMC following the process described in El-Badry et al. (2018). Although this semi-analytical model was built and applied on a dark-model only merger tree, we decided to apply it to FIREbox as is to serve as a point of comparison. We describe here the different steps of this process.

Given a galaxy identified in FIREbox, we compute at each snapshot, corresponding to the galaxy evolution time t, its star formation rate surface density as Σ_Ψ[t] = Ψ[t]/πR_d[t]² where R_d[t] is the scale length of the gas disc at time t. This scale length is assumed to be a function of the halo specific angular momentum (Mo et al. 1998) such that $$R_{\mathrm{d}}[t]=\frac{\lambda }{\sqrt{2}}{R}_{\mathrm{v}\mathrm{i}\mathrm{r}}[t],$$(3)

where λ is the halo spin parameter taken to be fixed at a typical value of λ = 0.035 (Bullock et al. 2001). Following El-Badry et al. (2018) the surface density of GMCs is then taken to be $${\mathrm{\Sigma }}_{\mathrm{G}\mathrm{M}\mathrm{C}\mathrm{s}}=5\times 1.2{10}^3\left(\frac{{\mathrm{\Sigma }}_{\mathrm{\Psi }}}{10{\mathrm{M}}_{\odot }{\mathrm{k}\mathrm{p}\mathrm{c}}^{-2}{\mathrm{y}\mathrm{r}}^{-1}}\right){\mathrm{M}}_{\odot }{\mathrm{p}\mathrm{c}}^{-2}.$$(4)

The expression of the massive cluster formation efficiency in Equation (2) implies that very few massive clusters can be formed when the GMC surface density is low (Σ_GMC ≪ Σ_crit), and it reaches the plateau value α_Γ when the GMC surface density is high (Σ_GMC ≫ Σ_crit). Following El-Badry et al. (2018), and supported by the values inferred from zoom-in simulations of collapsing clouds in Grudić et al. (2018), we set Σ_crit = 3000 M_⊙pc⁻² and β_Γ = −1. We finally set α_Γ such that it is consistent with the total mass of massive clusters sampled in our set of zoom-in simulations used in Paper III. This gives us a value α_Γ = 0.06. The strong difference between this value and the one proposed in El-Badry et al. (2018) of α_Γ = 2.1 × 10⁻³ comes from the fact that they only consider massive clusters that survive to the present-day, while we aim to estimate the formation of all massive clusters, including those that get disrupted on short timescales. These clusters, which are also included in the zoom-in simulations, can contribute to the production of merging BBHs before their disruption, and even afterwards, if we consider the long time delays that BBHs can take to reach coalescence.

Once we have obtained the massive cluster formation rate from the galaxy SFR Ψ[t] and the estimated formation efficiency Γ_massive[t], individual clusters masses are sampled from a power-law with a slope randomly drawn from the fit presented in Sect. 2.2.1 and with a minimum mass, M_min = 6 × 10⁴ M_⊙.

2.2.3 ‘SFRpeak’ model: Extreme episodes of star formation

The third approach presented here relies on the idea that massive star clusters are not ordinarily formed in the most typical inter-stellar medium, but are a natural by-product of intense episodes of star formation. These episodes could themselves be related to environments with high gas density such as disks in high redshift galaxies (see e.g. Kruijssen 2015) or to the outcome of galaxy interactions and mergers (see e.g. Renaud et al. 2016, 2017; Li et al. 2017). To determine the epochs of cluster formation in each of the FIREbox galaxies, we identify extreme episodes of star formation using the galaxy star formation history described in Sect. 2.1. In practice, we compute at each snapshot the ratio $$\gamma [t]=\frac{\mathrm{\Psi }[t]}{{\overline{\mathrm{\Psi }}}_{1\mathrm{G}\mathrm{y}\mathrm{r}}[t]},$$(5)

where ${\overline{\mathrm{\Psi }}}_{1\text{Gyr}}[t]$ is the mean SFRD around time t. We use an arbitrary value of one Gyr to compute this mean SFR, as it shows reasonable results when applied to our set of zoom-in simulations (see appendix). This ratio can then be used as a proxy for the formation efficiency of massive star clusters Γ_massive[t].

To predict the total mass of massive star clusters that should be sampled in each galaxy, we used the populations of star clusters sampled in the zoom-in FIRE-2 simulations of individual galaxies and built a linear fit estimator. Similarly to the relation observed in a wide range of galaxies between the total mass of their GCs and their virial mass (see e.g. Blakeslee et al. 1997; Harris et al. 2017), we relate the total mass of massive star clusters (M_cl ≥ 6 × 10⁴ M_⊙) that have ever formed in a given galaxy with its present-day stellar mass. With our set of zoom-in simulations, we obtain the relation log(ΣM_cl,massive) ≃ 5.97 + 2.01 × log(ℳ_⋆/10⁹ M_⊙) with a standard deviation σ ≃ 0.32.

The combination of the estimated total mass of massive clusters to sample M_cl,massive with the locations of massive cluster formation indicated by γ[t] provides us with the massive cluster formation rate of this third model. Here again, individual clusters masses are sampled from a power-law with a slope randomly drawn for each galaxy from the fit presented in Sect. 2.2.1 and with a minimum mass M_min = 6 × 10⁴ M_⊙.

2.2.4 Sampling of other cluster properties

Regardless of the sampling procedure for the number and masses of individual star clusters, we use the same following prescriptions to predict their other properties. To associate a value of metallicity to each cluster, we first determine the mean metallicity of its host galaxy at all snapshots following the same method used for the computation of the star formation histories (see Sect. 2.1). Cluster metallicities are then interpolated from their formation times, with a scatter of 0.1 dex to account for a non-homogeneous metallicity in the considered galaxy at a given time. Although the dispersion of metallicities has been observed to be greater than this value for distant galaxies (see e.g. Troncoso et al. 2014) or lower in nearby galaxies (see e.g. Williams et al. 2021), we find that this value gives is a good match to the metallicity dispersion for massive star clusters sampled in the individual zoom-in simulations.

Finally, the cluster half-mass radius is determined by sampling from a size-mass relation similar to that of Grudić et al. (2021) $$r_{\mathrm{h}\mathrm{m}}=3\mathrm{p}\mathrm{c}{\left(\frac{M_{\mathrm{c}\mathrm{l}}}{{10}^4{\mathrm{M}}_{\odot }}\right)}^{1/3}{\left(\frac{Z}{{\mathrm{Z}}_{\odot }}\right)}^{1/10},$$(6)

with a log-normal scatter of ±0.4 dex. Compared to the original formula presented in Grudić et al. (2021), here, we have removed the dependence of the clusters half-mass radii on the physical properties of their host GMCs $(\propto M_{\mathrm{G}\mathrm{M}\mathrm{C}}^{1/5}{\mathrm{\Sigma }}_{\mathrm{G}\mathrm{M}\mathrm{C}}^{-1})$, as we do not identify these structures individually.

2.3 BBH mergers in large populations of massive star clusters

With almost a thousand galaxies identified in the simulation, and up to thousands of massive clusters sampled in each of them, it is certainly impractical to integrate them all with any code for cluster evolution (with the possible exception of some recently developed rapid population synthesis codes for cluster evolution, such as the RAPSTER code Kritos et al. 2024). We develop a predictive model to estimate the BBH mergers formed in large populations of massive star clusters, by taking advantage of the 1500 massive clusters already integrated with the code CMC in our previous studies (Paper II; Paper III). These clusters have been sampled in six different zoom-in cosmological simulations of individual galaxies from the FIRE-2 project, taking into account the impact of the galactic environment on their evolution through tidal fields and dynamical friction.

CMC models collisional stellar dynamics in star clusters with an orbit averaging technique using a Hénon Monte Carlo approach (Hénon 1971). It is based on the hypothesis that clusters with a sufficiently large number of particles (≳10⁵) evolve mainly by two-body encounters that can be modelled as a statistical process. This condition translates into the initial mass threshold M_cl ≥ 6 × 10⁴ M_⊙ that we applied to our clusters. Various processes relevant to the formation of BBHs are taken into account, including two-body relaxation (Joshi et al. 2000), three-body binary formation (Morscher et al. 2015), and direct integration of small-N resonant encounters with post-Newtonian corrections (Rodriguez et al. 2018). The evolution of stars and stellar binaries in each cluster is modelled using the rapid population synthesis code COSMIC (Breivik et al. 2020).

From our collection of 1500 clusters evolved with CMC we build a 3-dimensional grid in the parameter space (M_cl, r_hm, Z_cl). We show this 3D grid, and its projections on the planes (M_cl, r_hm) and (M_cl, Z_cl), in Figure 3, with the number of BBH mergers predicted by CMC represented with different colours. The general trends that emerge from this collection of star clusters are that the most massive and densest clusters are the most likely to produce a large number of BBH mergers (see the upper left-hand panel of Figure 3), and that the most metal-rich clusters are also the most inefficient (upper right-hand panel). The latter result is due both to stellar evolution, which is particularly affected by metallicity through stellar winds and mass loss (see e.g. Vink et al. 2001; Mapelli et al. 2010; Spera et al. 2015; Di Carlo et al. 2020), and to the fact that the majority of these clusters appear at very late times in their respective galaxies and have therefore not had the opportunity to evolve long enough for their BBHs to merge.

We note that the peculiar distribution of cluster metallicities, with several star clusters located on lines of constant metallicity, is a natural result of the cluster formation model used in Paper Bruel et al. (2024). In this model, the metallicity of each star cluster is directly inherited from its parent GMC. As the more massive and densest GMCs can form more than one massive star cluster, this results in several clusters sharing the exact same metallicity in our 3D grid.

For every new cluster sampled in a given galaxy, we compute the Euclidean distance in log space with all the clusters in the grid as: $$d_{\mathrm{i}}^2={(\log \left(\frac{M}{M_{\mathrm{c}\mathrm{l}}}\right))}^2+{(\log \left(\frac{r}{r_{\mathrm{h}\mathrm{m}}}\right))}^2+{(\log \left(\frac{Z}{Z_{\mathrm{c}\mathrm{l}}}\right))}^2,$$(7)

where M, r, and Z are respectively the mass, radius and metallicity of the new cluster considered. To estimate the number of BBH mergers that this new cluster would produce, we first identify the 10 closest clusters in the grid according to the aforementioned distance (Equation (7)). We then compute the average of the number of BBH mergers predicted by CMC for these 10 neighbouring grid clusters, weighted by the inverse of their distance to the new cluster. Finally, we sample the properties of these merging BBHs, namely the two component masses and the time elapsed between cluster formation and BBH merger, from the collection of merging BBHs produced in the 10 neighbouring clusters (again including inverse distance weighting).

In this form, the grid-matching method only takes into account the intrinsic properties of each star cluster sampled and does not consider the impact of the galactic environment on its evolution. Although it has been shown that this effect can play an important role in determining the lifetime of star clusters (see e.g. Rodriguez et al. 2023), in this study we do not take into account the influence of the galactic environment on cluster evolution. We make here the assumption that the clusters already evolved with CMC in our set of 6 zoom-in simulations (Paper III) capture a range of tidal fields wide enough to be representative of a variety of distinct galactic environments, and that the selection of clusters in the grid is sufficiently random to ensure that the impact of these mechanisms on a large population of clusters is statistically respected. We leave it for future work to develop a more realistic and robust tool to predict the number and properties of BBH mergers formed in massive, and low-mass star clusters, embedded in different galactic environments.

Although we do not expect an exact match for every cluster sampled, especially since the effects of galactic environment are not taken into account in this grid-matching process, this method shows promising results when considering large populations of massive star clusters (see Appendix for a comparison with previous analyses from Paper III).

Fig. 4 Formation rate of massive clusters (Mcl ≥ 6 104 M⊙) in FIREbox galaxies using the cluster formation model ‘Gamma’ (left panel), the ‘EB18’ model (middle panel) and the ‘SFRpeak’ model (right panel). A re-weighting factor corresponding to the halo mass of the host galaxy (see Sect. 2.1) is applied to each BBH obtained. Different colours indicate the contribution of galaxies grouped in different bins of stellar mass, with the same bins as those already used in Figure 2. The black solid line shows the mean total formation rate density obtained from 20 different random realisations, and the grey shaded area is the 90% credible interval.

3 Results

3.1 Massive star clusters sampled in FIREbox

The formation rates of massive clusters in all the galaxies identified in FIREbox with the three different methods presented here are displayed in Figure 4. All models predict a total formation rate density of massive star clusters in the range 10⁵−10⁶M_⊙Mpc⁻³Gyr⁻¹ at z = 0 and a peak of massive cluster formation at around redshift z ~ 2–3, but with different values for this peak. The ‘SFRpeak’, in particular, finds a formation rate density of massive clusters at high redshifts one order of magnitude higher than the two other models. It results in a decrease of this formation rate with a steeper, but also more stochastic, slope. The ‘sawtooth’ shape is a direct consequence of the small number of massive galaxies in FIREbox. Indeed, the ‘SFRpeak’ model places the epochs of massive cluster formation precisely in the most extreme episodes of star formation, and these massive galaxies each have important episodes of star formation at different epochs. This results in a formation rate density of massive star clusters that appears highly stochastic in such a reduced cosmological volume. In reality, we do not expect such strong stochasticity in the formation rate density of massive star clusters when considering a much larger cosmological volume. The ‘Gamma’ and ‘EB’ models, however, find a much smoother formation rate over time. These features are expected from the very construction of these models, as they both follow the SFRD in each galaxy, and therefore the global SFRD in FIREbox.

In all three models the massive cluster formation rate is dominated at all redshifts by the more massive galaxies (blue and orange lines), even though there are far fewer of them (see Figure 2). It is clear that these predictions depend heavily on the physical processes that govern the growth of galaxies, and the omission of AGN feedback here plays a major role in the relative importance of the very massive galaxies (ℳ_⋆ ≥ 10¹¹ M_⊙). However, we do not expect that the inclusion of such feedback mechanisms would change the global picture of massive star clusters preferentially forming in massive galaxies (ℳ_⋆ ≥ 10¹⁰ M_⊙). On the other hand, it is also clearly apparent in the three panels that the low-mass galaxies (ℳ_⋆ ≤ 10⁹ M_⊙) have a negligible contribution to the cosmological population of massive star clusters. Predictions for the intermediate-mass galaxies are fairly consistent across the three models.

This feature can be directly related to the excess of star formation that originates from the massive non-quiescent galaxies in FIREbox (discussed above, see Sect. 2.1). We note here that, due to the number of massive non-quiescent galaxies in FIREbox, the excess of star formation at low redshifts translates into an excess of galaxies with very high stellar masses (see Figure 1). For the ‘Gamma’ and ‘SFRpeak’ models, which both use fits with galaxy stellar mass to estimate the populations of massive star clusters, this implies a local formation rate density of massive star clusters somewhat higher than the one predicted with the ‘EB18’ model. However, this excess is associated with clusters with typically very high metallicities, which will therefore contribute very little to the overall production of merging BBHs (see the 3D grid in the left-hand panel of Figure 3 for the relation between cluster metallicity and BBH formation).

3.2 Cosmological population of BBH mergers

We combine all the star clusters sampled with the 3D grid-matching technique described in Sect. 2.3 to create populations of dynamically formed merging BBHs across cosmic time. The BBH merger rate densities obtained with the three models presented in Sect. 2.2 are shown in Figure 5. The three populations of BBH mergers give comparable estimates of the local BBH merger rate density from massive clusters in the range ℛ_GCs ~ 1–10 Gpc⁻³yr⁻¹. These estimated values are in qualitative agreement with predictions from previous studies for the contribution of GCs to the astrophysical population of merging BBHs (see e.g. Rodriguez et al. 2018; Mapelli et al. 2022). Fishbach & Fragione (2023) infer a local BBH merger rate from GCs of ℛ_GCs(z = 0) = ${10.9}_{-9.3}^{+16.8}$ Gpc⁻³yr⁻¹ and that increases to ${58.9}_{-46}^{+149.4}$ Gpc⁻³yr⁻¹ at z = 1. These predictions are more in line with those of the ‘SFRpeak’ model.

By way of comparison, the LVK Collaboration infers an intrinsic overall BBH merger rate density of ℛ_BBHs ~ 17.9–44 Gpc⁻³yr⁻¹ at a fiducial redshift z = 0.2 (Abbott et al. 2023a). This value is inferred from 69 GW observations corresponding to confident BBH events with a false alarm rate (FAR) <1 yr⁻¹, and takes into accounts the particular selection effects of the GW observatory network. In particular, the blue shaded region in Figure 5 shows the best-fit model of the merger rate density in the form ℛ(z) ∝ (1 + z)^κ when using their POWER LAW+PEAK mass model (taken from Figure 13 in Abbott et al. 2023a). The combination of the three cluster formation models ‘Gamma’, ‘EB18’, and ‘SFRpeak’ with the grid-matching algorithm predicts values of the BBH merger rate density at z = 0.2 of ${4.8}_{-2.0}^{+3.1}$ Gpc⁻³yr⁻¹, ${2.8}_{-0.6}^{+0.7}$ Gpc⁻³yr⁻¹ and ${6.8}_{-5}^{17}$ Gpc⁻³yr⁻¹ respectively. We emphasize here that the errors reported here do not represent the absolute uncertainties of our model, but rather the fluctuation inherent the stochasticity associated with the random processes taking place inside each of the three cluster formation models. The differences between these values and the one inferred by the LVK indicate that although it is not dominant, the formation of merging BBHs in massive star clusters could very well represent a non-negligible fraction of the astrophysical population.

It is particularly interesting to note that the three models have quite different slopes to describe the evolution of ℛ_GCs with redshift. In the hypothesis that we could accurately measure the evolution of the merger rate density with redshift and that we could separate BBH mergers coming from different formation channels (using, for example, the value of the chirp mass ℳ_c, or values of the spin projections, see e.g. Arca Sedda et al. 2020; Antonelli et al. 2023), constraining in particular the redshift evolution of the dynamical BBH merger rate density would be particularly instructive for better understanding the formation of star clusters through cosmic time.

All three models of massive star cluster formation predict that the merger rate density of BBHs formed in massive star clusters is highly dominated by the most massive galaxies. This result seems to direct the identification of the host galaxies for the ‘dynamical’ BBHs towards massive galaxies (ℳ_⋆ ≥ 10¹⁰ M_⊙). Through a more detailed analysis, such a feature could potentially be used in cosmological studies to assign astrophysically motivated probabilities to potential host galaxies for the most massive BBH mergers observed (see e.g. Mastrogiovanni et al. 2023, for different methods on how to use GW events as dark standard sirens for cosmology measurements).

Fig. 5 BBH merger rate density in FIREbox estimated from the results of our cluster sampling algorithms (model ‘Gamma’, ‘EB18’, and ‘SFRpeak’ presented from left to right respectively) combined with predictions of the 3D grid of massive clusters integrated with CMC. A re-weighting factor corresponding to the halo mass of the host galaxy (see Sect. 2.1) is applied to each BBH obtained. Different colours indicate the contribution of galaxies grouped in different bins of stellar mass. The black solid line shows the mean total merger rate density obtained from 20 different random realisations, and the grey shaded area is the 90% credible interval. The blue shape shows the evolution of the BBH merger rate density with redshift as reported by the LVK Collaboration in Abbott et al. (2023a).

3.3 Properties of dynamically formed BBHs

We show in Figure 6 the distribution of BBH primary masses for the BBH mergers occurring at z ≤ 1 estimated from the results of our three cluster formation models (‘Gamma’, ‘EB18’, and ‘SFRpeak’) combined with predictions of the 3D grid of massive clusters already integrated forward in time with CMC. We compare these results with the astrophysical primary mass distribution inferred by the LVK Collaboration from the same 69 GW observations with FAR < 1 yr⁻¹ using their POWER LAW+PEAK mass model (here in blue colour, taken from Figure 10 in Abbott et al. 2023a). We would like to stress out here that this mass distribution inferred from GW observations, although supposed to represent the intrinsic population of BBH mergers, corresponds really to a best-fit parametric model. In particular, it relies on the assumption that the distribution of primary BH masses does not evolve with redshift. Care must therefore be taken when comparing a simulated astrophysical population with this distribution, and our aim here is not to reproduce it, but rather to identify common characteristics.

The three differential merger rates obtained in this study fall below the differential merger rate inferred by the LVK Collaboration, as is naturally expected from the lower merger rate densities predicted (see Sect. 3.2 and Figure 5). With all three cluster formation models, we obtain a distinctive feature at around m_BH1 ~ 35 M_⊙ that is notably consistent with the location of the peak reported at ${34}_{-4.0}^{2.6}$ M_⊙ in Abbott et al. (2023a) using the POWER LAW+PEAK model (indicated as the vertical dashed line in Figure 6). Comparing these predicted astrophysical primary mass distributions with the mass distribution of all the BBH mergers directly obtained from our CMC runs (see Figure 8 in Bruel et al. 2024), it is interesting to note that this peak does not arise from the global distribution of all BBH mergers simulated across a few massive star clusters, but really from predictions over a cosmological population of them. We did not find any particular sampling effect responsible for its presence in our predictions presented here. Further in-depth analyses of stellar evolution and cluster dynamics in CMC would be necessary to understand this feature in more detail.

It is also interesting to note that all three models predict a distribution of BBH primary masses that is clearly different from predictions of the isolated evolution channel obtained with rapid population synthesis codes (see e.g. van Son et al. 2023). The merging BBHs formed through dynamical interactions in massive star clusters can have masses that extend well above the pair-instability mass limit (see e.g. Marchant et al. 2019), with an almost continuous distribution of primary masses from ~50 to ~80 M_⊙. The predictions of our astrophysical BBH primary masses in this mass range are remarkably consistent with the high-mass end of the POWER LAW+PEAK model, indicating a clear preference for the dynamical formation channel for such massive events.

Although our predictions on the physical properties of merging BBHs formed in massive star clusters do not provide an exact match to the astrophysical population inferred from GW observations, these features support the idea that different formation channels may populate different regions of the BBH mass distribution. Combined with BBHs formed through the isolated evolution of massive stellar binaries in the stable mass transfer channel (see e.g. van Son et al. 2022b, 2023), and with BBHs formed in low and intermediate mass clusters (see e.g. Torniamenti et al. 2022), both the predicted merger rate densities and the mass distribution could very well align with the astrophysical population inferred from GW observations.

We find no evolution of this BBH primary mass distribution with redshift for m_BH1 ∈ [10–50] M_⊙, but a clear preference for the production of massive BBHs m_BH1 > 65 M_⊙ in star clusters formed at z ≥ 1.5. This trend naturally arises from the fact that massive stars at low metallicity lose less mass in stellar winds, thus forming more massive compact objects, and low metallicity environments are preferably found at high redshifts.

Looking specifically at the galaxies in which these type of extreme mergers take place, we find that the majority of these events take place in the most massive galaxies (with present-day stellar masses ℳ_⋆ ≥ 10¹⁰ M_⊙). Although a significant fraction of these massive galaxies should in reality be quiescent now, due to several feedback mechanisms not modelled in FIREbox, the quiescent fraction is observed to be less important at higher redshifts (see e.g. Fontana et al. 2009), and the absence of these feedback mechanisms in FIREbox therefore has less of an impact on these predictions. This result indicates that the detection of such an extreme BBH merger could most certainly translate into some strong constraints on the physical properties of its host galaxy. With the deployment of the next generation of GW interferometers, we may be able to identify the exact host galaxy of some of the nearest BBH mergers (see e.g. Mo & Chen 2024). The potential correlations observed between the physical properties of the host galaxies and those of the BBH mergers would make it possible to go even further in analysing the astrophysical origin of merging BBHs.

We note here that the three distributions are quite similar in shape because all the sampled BBHs come from the same 1500 clusters, which produce in total 7631 BBH mergers. With around a million massive star clusters sampled in total, each cluster in our 3D grid can be expected to be selected around 1000 times in average. In practice, we find that the average number of occurrences of each massive star cluster as nearest neighbour in the grid is indeed around 1000, with a handful of clusters selected more than 10 000 times. We did not find any particular physical properties of these star clusters that are sampled most frequently. This redundancy is then naturally reflected in the properties of the sampled BBHs. In practice, we find that the median number of occurrences of each BBH merger is around 100 for all three models of cluster formation. However, a small number of BBHs can be repeated several thousand times. Here again, we find no particular preference for the physical properties of the most frequently sampled BBHs.

Finally, we point out that a non-negligible number of BBH mergers sampled with all three cluster formation models in FIREbox have primary masses well above 100 M_⊙. With the three models of cluster formation ‘Gamma’, ‘EB18’, and ‘SFRpeak’, we find respectively 2315, 966, and 3820 of such BBH mergers taking place at redshifts z < 1 in the volume of FIREbox. We find that these massive BBH mergers are sampled from 27 different clusters in our grid of massive star clusters simulated with CMC. Around half of these clusters have well-identifiable initial conditions and are among the most massive and dense in our catalogue. These environments are precisely the most favourable hosts for stellar collisions and repeated BBH mergers (see e.g. Fragione & Rasio 2023). However, the other half of these clusters have initial properties that lie in the bulk of the total catalogue, and the formation of massive BBH mergers in these environments appears to be a much more stochastic process. These predictions of the differential merger rate at the highest masses are therefore also the most uncertain, and we advise taking these results with caution.

The lack of detection of BHs more massive than 100 M_⊙ in LVK data seems to indicate that the models of stellar evolution and cluster dynamics presented here predict too massive objects (see however Wadekar et al. 2023, where higher harmonics are employed in the GW templates and, using a detection threshold p_astro > 0.5, they find 14 additional BBH mergers in LVK data, with the most massive one having M_BH1 = ${300}_{-120}^{+60}$ M_⊙). Stronger observational constraints on the existence and detectability of such massive mergers will certainly emerge thanks to a greater number of GW observations and the contribution of new-generation detectors (see e.g. Franciolini et al. 2024).

Fig. 6 Astrophysical BBH primary mass distribution of all the local (z ≤ 1) BBH mergers sampled in FIREbox. Different colours correspond to different massive star cluster formation models. We apply a re-weighting factor corresponding to the halo mass of the host galaxy (see Sect. 2.1) to each BBH merger. All the distributions are normalised to the predicted value for the BBH merger rate density at z = 0. The blue curve shows the differential merger rate as a function of primary mass from GWTC-3 analyses (Abbott et al. 2023a), with the shaded region showing the 90% credible interval.

3.4 Host galaxies of BBHs formed dynamically in massive star clusters

Next, we were interested in exploring the properties of the galaxies that contain the massive star clusters in which merging BBHs form dynamically. In Figure 7, we show the distribution of galaxy present-day stellar masses and merger times of the BBH mergers sampled with the cluster formation ‘Gamma’ in FIREbox. These mergers preferentially originate from massive star clusters formed at high redshifts (z ≳ 2) and in massive galaxies (with present-day stellar masses ℳ_⋆ ≳ 10¹¹M_⊙). We compare, for the three cluster formation models, the marginalised distributions of BBH merger over time (subplot on the left) and over galaxy masses (subplot on top). All three models qualitatively agree on the characteristic epochs and galaxy stellar masses for the production of BBH mergers in massive star clusters.

We also compare these distributions with all star formation in FIREbox. Marginalised over galaxy masses this gives the total SFRD in the cosmological volume, already presented in Figure 2, and marginalised over time this gives the SFR-weighted density of galaxies at z = 0. As mentioned before (Sect. 2.1), most of star formation takes place around 10 billion years ago (z ~ 2) in galaxies that are presently massive (11 ≤ log(ℳ_⋆/M_⊙) ≤ 12).

Inspecting the marginalised distributions in more detail, we find an even more prominent production of BBH mergers at around z ≳ 3 when compared with the normalised SFRD. The fact that relatively fewer massive clusters produce BBH mergers at low redshifts is the combined effect of higher metallicities and the necessary long delay times from cluster formation to the production of BBHs and all the way to their mergers due to GW emission alone. Comparing the galaxy mass functions, we find a similar preference for the highest present-day stellar masses, but also a clear relative dearth of low mass galaxies (ℳ_⋆ ≤ 10¹⁰ M_⊙) contributing to the production of massive star clusters in which merging BBHs form. This difference arises from the fact that, although a number of stars form within them, there are practically no massive clusters in these low-mass galaxies.

The increasing sensitivity of current ground-based interferometers, as well as the future contribution of next-generation detectors in the coming decades, will certainly put strong constraints on the redshift evolution of the BBH merger rate density (see e.g. van Son et al. 2022a). This will be invaluable in helping us to better understand the mechanisms involved in the formation of these BBHs. Furthermore, these next-generation detectors will likely make it possible to localise some GW events in volumes small enough to allow for the identification of the host galaxies (see e.g. Mo & Chen 2024). We predict here that the majority of merging BBHs formed in massive star clusters originate from massive galaxies. This trend towards the formation of merging BBHs in massive galaxies is in agreement with previous predictions exploring the isolated evolutionary channel (see e.g. Artale et al. 2019b,a), although it has been shown to be strongly dependent on assumptions about star formation and metallicity distribution (see e.g. Santoliquido et al. 2022; Srinivasan et al. 2023). On the other hand, Srinivasan et al. (2023) find that the stellar binary progenitors of detectable BBH mergers tend to form preferentially at lower redshifts (z ≲ 1) and in dwarf galaxies. Under the assumption that different formation channels of detectable BBH mergers typically take place in different environments and at different epochs, the identification of host galaxies for certain future GW events could shed new light on the question of the astrophysical origin of double compact objects.

Fig. 7 Merger times and galaxy present-day stellar masses of the BBH mergers in FIREbox estimated from the ‘Gamma’ cluster formation model (colormap). The left and top subplots show the marginalised distributions obtained with the three different models of cluster formation (coloured lines), as well as the corresponding distributions for all star formation in the 979 galaxies identified inside FIREbox (black solid line, with arbitrary normalisation).

4 Discussion

4.1 Caveats

One caveat to this study lies in the limits imposed by the cosmological volume simulation FIREbox itself. Indeed, as described above (Sect. 2.1) and in more detail in Feldmann et al. (2023), there is a notable excess of star formation at z < 1 in the simulation compared to observations. For the most part, this excess of star formation takes place in galaxies that evolve to become the most massive at present-day (ℳ_⋆ ≥ 10¹⁰ M_⊙). However, observations show that the majority of these massive galaxies should be quiescent at such redshifts (see e.g. Leja et al. 2022). This suggests that stellar feedback alone, as it is implemented in the FIRE-2 physics model, is not sufficient to explain the observed quenching of galaxies. Additional mechanisms such as feedback from cosmic rays or AGN, which are not implemented in FIREbox, could certainly help to reduce this inconsistency. We refer to Feldmann et al. (2023) for a more in-depth analysis of this excess of star formation in FIREbox.

With an excess of star formation at z = 0 of a factor ~3 in FIREbox (see Figure 2), we can put an upper limit of the same factor on the over-prediction of the local BBH merger rate density originating from this overabundance of stars. However, since a certain fraction of local BBH mergers can be produced in massive star clusters formed at much higher redshifts where FIREbox predicts a SFRD in good agreement with observations, the excess of star formation at low redshifts has in reality a more moderate impact. In practice, we find that, for the three models of cluster formation ‘Gamma’, ‘EB 18’, and ‘SFRpeak’ respectively, 21, 25 and 34% of the BBH mergers taking place at redshifts z ≤ 0.5 originate from massive star clusters formed at z > 1. Considering that, for the local BBH merger rate density, one fourth of the BBHs are correctly produced in clusters formed at z > 1 and the other three fourths are overestimated by a factor 3 for clusters formed at z < 1, this translates into a worst-case over-prediction of a factor of about 2–3.

This over-prediction contributes among other uncertainties associated with our cluster formation models and our grid-matching algorithm for predicting populations of dynamical BBH mergers. Further analyses need to be performed to quantify in more detail each of these uncertainties.

A second caveat lies precisely in the grid-matching algorithm, and in particular in the incomplete parameter space covered by the collection of our star clusters already integrated with the code CMC. While, for most of the massive star clusters sampled, the closest neighbours in the grid have relatively similar physical properties (see e.g. Figure A.1), a certain fraction of them can actually fall in a region of the parameter space that is not populated. This is notably the case for clusters with low metallicities, where our grid is only scarcely filled (see the 3D grid in Figure 3) This lack of low-metallicity clusters is largely due to the fact that FIRE cosmological simulations, including m12i which was used in Paper I; Paper II; Paper III, typically have too little star formation at high redshifts (see e.g. Wetzel et al. 2023).

This poses two problems. As FIREbox uses the same physical model, we can expect that there will also be a lack of star formation at high redshifts, and therefore a lack of low-metallicity star clusters sampled in this cosmological volume simulation. Furthermore, the scarcity of low-metallicity clusters in our 3D grid means that our interpolation algorithm is unlikely to be very robust in this region of the parameter space.

To account for this second issue, we looked at the impact of adding massive star clusters at low metallicities to the already existing 3D grid presented in Sect. 2.3. We select a set of 14 metal-poor clusters (Z = 0.01 Z_⊙) from the total set of 148 systems presented in Kremer et al. (2020). These clusters have all been evolved with the same code CMC although with slightly different physical prescriptions (including the initial fraction of stellar binaries, which is taken as a constant value f_b = 5% in Kremer et al. (2020) while we considered a mass-dependent binary fraction following van Haaften et al. (2013)). After using the updated 3D grid to run the grid-matching algorithm on the populations of massive star clusters sampled with our three formation models, we find no major difference in the predictions already presented in Section 3.

This does not come as a surprise, as the addition of only 14 star clusters to a grid containing already 1500 would not be expected to have a significant impact. To quantify the impact of this empty region of the parameter space, more clusters could potentially be run with CMC. Different codes for star cluster evolution could also be used for a more global comparison. We leave the exploration of these aspects to future studies.

4.2 Comparison with other studies

Our three models of massive cluster formation predict a total formation rate density of massive star clusters in the range 10⁵−10⁶ M_⊙Mpc⁻³Gyr⁻¹ at z = 0 and a peak of massive cluster formation at around z ~ 2–3 (see Figure 4). These high values of the massive cluster formation rate at z > 2 are in qualitative agreement with observations of star-forming galaxies at high redshifts (see e.g. Shapiro et al. 2010). Integrated over time, we find a value of the cluster mass formed per unit volume of ρ = 3.2, 1.7, and 7.3 × 10⁷ M_⊙Mpc⁻³ with our three cluster formation models ‘Gamma’, ‘EB18’, and ‘SFRpeak’ respectively. As a comparison, by constructing their own model of GC formation with a present-day mass density of GCs in the Universe consistent with its empirical value and a current mass function consistent with the observed distribution of galactic GCs, Antonini & Gieles (2020) find ρ = 2.4 × 10⁷ M_⊙Mpc⁻³ with an estimated uncertainty of a factor of 2. As a comparison, Antonini & Gieles (2020) construct their own model of GC formation using a present-day mass density of GCs in the Universe consistent with its empirical value, and a present-day mass function consistent with the distribution of observed galactic GCs. They find a value of the cluster mass formed per unit volume ρ = 2.4 × 10⁷ M_⊙Mpc⁻³ with an estimated uncertainty of a factor of 2. This falls reasonably within the predictions of our models.

Using the catalogue of GW sources GWTC-3 (Abbott et al. 2023a), and assuming that the GC formation rate number density follows a Madau-like (Madau & Dickinson 2014) functional form, Fishbach & Fragione (2023) infer a value for the GC formation rate at redshift 2 of ~2 × 10⁶ M_⊙Mpc⁻³Gyr⁻¹. This value is in good agreement with the predictions from the two models ‘Gamma’ and ‘EB18’. The cluster formation model ‘SFRpeak’ predicts overall a more important formation rate of massive star clusters. This feature could very well be explained by the fact that it uses scaling relations with galaxy present-day stellar masses, and the most massive galaxies are found to be in excess inside FIREbox. A more realistic cosmological volume population or a post-processing treatment for AGN feedback in FIREbox could potentially resolve this issue and produce more consistent predictions between our different cluster formation models.

5 Conclusions

In this work, we modelled the formation of massive star clusters in the cosmological volume simulations FIREbox (part of the FIRE-2 project). We used three different models of massive cluster formation empirically built from results of previous studies (Grudić et al. 2022; Rodriguez et al. 2023; Bruel et al. 2024) based on a cloud-to-cluster formation framework applied to the GMCs identified in a set of cosmological zoom-in simulations of individual galaxies. Extending on the work presented in these studies, we assembled all the massive star clusters already integrated with the Monte Carlo code CMC in a 3D grid in the parameter space (M_cl, r_hm, Z_cl). We then developed an algorithm to ‘predict the population of merging BBHs that would be produced by any large population of massive clusters. Combined with our different models of cluster formation applied to close to a thousand galaxies identified in FIREbox, this enabled us to create populations of merging BBHs that would have formed in these realistic environments.

We studied the properties of these BBH mergers and of their host galaxies. We found distinct features that could prove decisive in improving our understanding the formation history of massive star clusters and their contribution to the astrophysical population of BBH mergers in the local Universe. Our main findings are detailed below.

1. Different assumptions on the mechanisms driving the production of massive star clusters in various galaxies result in different predictions for their formation rate over cosmic time. In particular, assuming that the formation of massive star clusters occurs preferentially during episodes of intense star formation results in a formation rate density of massive star clusters with a larger peak at early times (z < 2) and with a steeper decrease up to the present day. In contrast, assuming that their formation follows the global star formation rate results in smoother evolutions over time.

2. Over different models of massive cluster formation, we predicted a consistent value of the local BBH merger rate density in the range ℛ_GCs ~ 1–10 Gpc⁻³yr⁻¹. Although this is below the value inferred by the LVK Collaboration from GW signals, it is a clear indication that this formation channel may contribute to a significant fraction of the list of detected BBH mergers. Different models of cluster formation model predict different vales for the local BBH merger rate density, but also for its evolution with redshift.

3. In agreement with previous studies, we found that the physical properties of BBHs formed in massive star clusters hold very distinct features. The possibility of stellar mergers and hierarchical BBH mergers in such dense environments allows for the formation of very massive BHs. These extreme systems, which are challenging to explain via different formation channels, bear the imprint of their dynamical formation in their unique physical properties.

4. The massive clusters that produce BBHs merging in the local Universe are found to exist in galaxies with similar properties as those driving star formation at a cosmological scale: they typically form at z ≳ 2 in galaxies with present-day stellar masses in the range 11 ≤ log(ℳ_⋆/M_⊙) ≤ 12. While the peaks of these two distributions lie in the same region, we also found a clear relative dearth of massive clusters producing BBH mergers in low mass galaxies (ℳ_⋆ ≤ 10⁹ M_⊙). It is worth noting that these results run counter to some predictions of the isolated evolution channel. This indicates that identifying the host galaxies of BBH mergers could provide valuable constraints on their formation channels and on the branching fraction between these channels.

This study represents the first attempt to estimate the populations of merging BBHs formed in massive star clusters in the cosmological volume simulation FIREbox. These results highlight the importance of modelling the formation of star clusters for predicting their contribution to the astrophysical population of local BBH mergers, but also to the redshift evolution of the BBH merger rate density. The ever-growing list of observed BBH mergers, the increasing sensitivities of the current GW interferometers, and the future advent of the third generation detectors will provide valuable information on the origin of the detected BBH mergers, the redshift evolution of the merger rate density, and, thus, the properties of the massive star clusters that produce them.

Acknowledgements

T.B. is supported by ERC Starting Grant No. 945155-GWmining, Cariplo Foundation Grant No. 2021-0555, MUR PRIN Grant No. 2022-Z9X4XS, MUR Grant “Progetto Dipartimenti di Eccellenza 2023-2027” (BiCoQ), and the ICSC National Research Centre funded by NextGenerationEU. T.B. and A.L. are supported by the ANR COSMERGE project, grant ANR-20-CE31-001 of the French Agence Nationale de la Recherche. C.L.R. acknowledges support from NSF Grant AST-2310362, NASA ATP Grant 80NSSC22K0722, a Alfred P. Sloan Research Fellowship, and a David and Lucile Packard Foundation Fellowship. This work was supported by the ‘Programme National des Hautes Energies’ (PNHE) of CNRS/INSU co-funded by CEA and CNES’ and the authors acknowledge HPC ressources from ‘Mesocentre SIGAMM’ hosted by Observatoire de la Côte d’Azur. This work made use of infrastructure services provided by the Science IT team of the University of Zurich (www.s3it.uzh.ch). The FIREbox simulation was supported in part by computing allocations at the Swiss National Supercomputing Centre (CSCS) under project IDs s697, s698, and uzh18. Support for M.Y.G. was provided by NASA through the NASA Hubble Fellowship grant #HST-HF2-51479 awarded by the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., for NASA, under contract NAS5-26555.

Appendix A Quantitative testing of the three sampling methods of massive star clusters combined with the grid-matching algorithm for predicting BBH mergers

As a consistency test, we apply the three different models of massive star cluster formation presented in §2.2 to the simulated galaxies m11e, m11h, and m12i (see I; II; III, for further analysis on the GMCs, massive star clusters, and BBH mergers in this simulated galaxy). The formation rates of massive star clusters obtained are shown in Figure A.2. As there is some uncertainty due to random draws in the cluster sampling models, in particular for the ‘Gamma’ and ‘SFRpeak’ models, we repeat the whole process 100 times to obtain different estimates of the star cluster populations in each simulated galaxy for each of the cluster models. There is clearly no exact match between the clusters sampled for the three cluster formation models and the ones obtained from identified GMCs, particularly at low redshifts. This is due to the fact that, while star formation still occurs at low redshifts in the simulation m12i, few massive and dense GMCs are found there, which translates into few massive star clusters. The ‘SFRpeak’ model, locating the epochs of massive cluster formation during extreme episodes of star formation only, provides a better match to the cluster formation rate at these low redshifts, but also over-predicts their formation at higher redshifts.

To measure the precision of the grid-matching method, we store the list of distances to the closest neighbour in the grid computed for all the clusters sampled for each of the three cluster formation models (one iteration of cluster sampling in the simulated galaxy m12i). We show the cumulative distributions of the obtained distances in Figure A.1. In this example, all the clusters sampled using the three models have a nearest neighbour in the 3D grid at a distance of less than 1, and more than 81% of them with d_i ≥ 0.15. This last value corresponds to a cluster whose physical properties are typically 1.2 times higher (or lower) than those of its nearest neighbour.

We then applied the grid-matching algorithm presented in §2.3 to estimate populations of dynamically formed merging BBHs. Figure A.3 shows the results of this test in the form of a BBH merger rate and compares with the merger rate obtained with massive star clusters sampled from individual GMCs and integrated forward in time with CMC(III). In this example, all models struggle to match the BBH merger rate at high redshifts, but they recover the local (z ≤ 1) BBH merger rate with a fairly good accuracy. This result is quite encouraging, as the three cluster formation models were built using the star clusters from our set of zoom-in simulations, and without taking their BBH mergers into account. We thus find that the combination of approximate models of cluster formation and BBH production developed here remains consistent with the more detailed study carried out on simulations of individual galaxies.

A more complete analysis shows that the discrepancy observed at redshifts z > 1 is not mainly due to the grid matching method, but rather to differences in the number and properties of massive clusters predicted to exist at these epochs. Indeed, examining in detail the predictions of the three proposed cluster formation models, it appears that they all tend to overestimate the formation of massive clusters at early times compared with clusters sampled directly from the GMCs. This is linked to the fact that our three models use the SFRD of a galaxy to estimate the formation of massive clusters within it, whereas the appearance of massive and dense GMCs, which were precisely used as the birthplaces of clusters in the previous study, is not necessarily correlated with this SFRD. This eventually results in a BBH merger rate always overestimated at these high redshifts. Since we aim to compare our predictions with LVK observations of BBH mergers which, with current detector sensitivities, are found at redshifts z ≤ 1, these discrepancies at high redshifts can be set aside for the moment.

Fig. A.1 Cumulative distribution of distances di obtained after sampling massive clusters in the simulation m12i with the three models ‘Gamma’, ‘EB18’, and ‘SFRpeak’, and applying the grid-matching scheme.

Fig. A.2 Massive star cluster formation rate estimated in the cosmological zoom-in simulations m11e, m11d, and m12i from top to bottom by sampling massive star clusters from individual GMCs (red, see Paper III), and from the results of our cluster sampling algorithms (models ‘Gamma’, ‘EB18’, and ‘SFRpeak’ respectively presented from left to right). For each model, the process is repeated 100 times and we show the mean formation rate (black solid line) as well as the 90% credible interval (grey shaded area).

Fig. A.3 BBH merger rate computed in the cosmological zoom-in simulations m11e, m11d, and m12i from top to bottom by sampling massive star clusters from individual GMCs and integrating them forward in time with the code CMC (red, see Paper III), and from the results of our cluster sampling algorithms (models ‘Gamma’, ‘EB18’, and ‘SFRpeak’ respectively presented from left to right) combined with predictions of the 3D grid. For each model, the process is repeated 100 times and we show the mean BBH merger rate (black solid line) as well as the 90% credible interval (grey shaded area).

References

All Figures

Fig. 1 Left: halo mass function (HMF) of the 979 galaxies identified inside FIREbox within the mass range 10 ≤ log(Mvir/M⊙) ≤ 14 (red). Vertical dashed lines show the range of the virial masses used for halo identification. Grey arrows indicate the re-weighting operated to match the HMF in FIREbox with a ‘universal’ HMF extracted from dark matter cosmological simulation (in black line, from Tinker et al. 2008). Right: galaxy stellar mass function (GSMF) for this same set of galaxies identified inside FIREbox (red), and with each galaxy re-weighted with the operation presented in the left-hand panel (orange). Two examples of GSMF built on observations of galaxies in the local Universe are also shown. The black solid line is the best-fit GSMF at z ~ 0.1 from Tomczak et al. (2014) and the dashed line is the median GSMF generated using the continuity model parameters of Leja et al. (2020) at z = 0. The shaded grey area corresponds to the 16–84% uncertainties associated with this model.

In the text

Fig. 2 Left: star formation rate density (SFRD) for the 979 galaxies identified inside FIREbox(black line) and for galaxies grouped in different bins of stellar mass (coloured lines). The number of galaxies in each mass bin is indicated in brackets. As a comparison, the grey dashed line shows the functional form of the cosmic SFRD proposed by Madau & Dickinson (2014). Right: mass-weighted average of the metallicity of newly formed stars in the galaxies identified in FIREbox. The black line indicates the average metallicity is in the entire volume, and the different colours indicate galaxies grouped in different bins of stellar mass. In both panels, the contribution of each galaxy is re-weighted as described in Sect. 2.1 (and as pointed out in the left-hand panel of Figure 1). By way of comparison with the different Figures shown in the FIREbox paper (Feldmann et al. 2023), we do not apply a bootstrapping method here and we therefore do not indicate error bars for these estimates. We refer to that particular paper for more details on these uncertainties.

In the text

Fig. 3 Upper left panel: 2D projection on the plane (Mcl, rhm) of the 1500 clusters integrated with CMC. Clusters in which no merging BBH are formed are indicated as grey dots, while the coloured dots represent the number of merging BBHs formed in each cluster. The diagonal dashed grey lines represent lines of constant density ${\rho }_{\mathrm{c}\mathrm{l}}\equiv M_{\mathrm{c}\mathrm{l}}/(4/3)\pi r_{\mathrm{h}\mathrm{m}}^3$. Upper right panel: 2D projection on the plane (Mcl, Zcl). Colours indicate the number of merging BBHs with the same norm as in the left-hand panel. Lower panel: 3D grid in parameter space (Mcl, rhm, Zcl).

In the text

Fig. 4 Formation rate of massive clusters (Mcl ≥ 6 104 M⊙) in FIREbox galaxies using the cluster formation model ‘Gamma’ (left panel), the ‘EB18’ model (middle panel) and the ‘SFRpeak’ model (right panel). A re-weighting factor corresponding to the halo mass of the host galaxy (see Sect. 2.1) is applied to each BBH obtained. Different colours indicate the contribution of galaxies grouped in different bins of stellar mass, with the same bins as those already used in Figure 2. The black solid line shows the mean total formation rate density obtained from 20 different random realisations, and the grey shaded area is the 90% credible interval.

In the text

Fig. 5 BBH merger rate density in FIREbox estimated from the results of our cluster sampling algorithms (model ‘Gamma’, ‘EB18’, and ‘SFRpeak’ presented from left to right respectively) combined with predictions of the 3D grid of massive clusters integrated with CMC. A re-weighting factor corresponding to the halo mass of the host galaxy (see Sect. 2.1) is applied to each BBH obtained. Different colours indicate the contribution of galaxies grouped in different bins of stellar mass. The black solid line shows the mean total merger rate density obtained from 20 different random realisations, and the grey shaded area is the 90% credible interval. The blue shape shows the evolution of the BBH merger rate density with redshift as reported by the LVK Collaboration in Abbott et al. (2023a).

In the text

Fig. 6 Astrophysical BBH primary mass distribution of all the local (z ≤ 1) BBH mergers sampled in FIREbox. Different colours correspond to different massive star cluster formation models. We apply a re-weighting factor corresponding to the halo mass of the host galaxy (see Sect. 2.1) to each BBH merger. All the distributions are normalised to the predicted value for the BBH merger rate density at z = 0. The blue curve shows the differential merger rate as a function of primary mass from GWTC-3 analyses (Abbott et al. 2023a), with the shaded region showing the 90% credible interval.

In the text

Fig. 7 Merger times and galaxy present-day stellar masses of the BBH mergers in FIREbox estimated from the ‘Gamma’ cluster formation model (colormap). The left and top subplots show the marginalised distributions obtained with the three different models of cluster formation (coloured lines), as well as the corresponding distributions for all star formation in the 979 galaxies identified inside FIREbox (black solid line, with arbitrary normalisation).

In the text

	Fig. A.1 Cumulative distribution of distances di obtained after sampling massive clusters in the simulation m12i with the three models ‘Gamma’, ‘EB18’, and ‘SFRpeak’, and applying the grid-matching scheme.
In the text

Fig. A.2 Massive star cluster formation rate estimated in the cosmological zoom-in simulations m11e, m11d, and m12i from top to bottom by sampling massive star clusters from individual GMCs (red, see Paper III), and from the results of our cluster sampling algorithms (models ‘Gamma’, ‘EB18’, and ‘SFRpeak’ respectively presented from left to right). For each model, the process is repeated 100 times and we show the mean formation rate (black solid line) as well as the 90% credible interval (grey shaded area).

In the text

Fig. A.3 BBH merger rate computed in the cosmological zoom-in simulations m11e, m11d, and m12i from top to bottom by sampling massive star clusters from individual GMCs and integrating them forward in time with the code CMC (red, see Paper III), and from the results of our cluster sampling algorithms (models ‘Gamma’, ‘EB18’, and ‘SFRpeak’ respectively presented from left to right) combined with predictions of the 3D grid. For each model, the process is repeated 100 times and we show the mean BBH merger rate (black solid line) as well as the 90% credible interval (grey shaded area).

In the text