Open Access
Issue
A&A
Volume 698, June 2025
Article Number A144
Number of page(s) 14
Section Extragalactic astronomy
DOI https://doi.org/10.1051/0004-6361/202452757
Published online 11 June 2025

© The Authors 2025

Licence Creative CommonsOpen Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

This article is published in open access under the Subscribe to Open model. Subscribe to A&A to support open access publication.

1. Introduction

The advent of gravitational wave (GW) detections has significantly expanded our understanding of stellar-mass black holes (BHs) and neutron stars (NSs, Abbott et al. 2016a,b, 2017a,b, 2021, 2023a, 2024). LIGO–Virgo–KAGRA (LVK) data offers valuable insights into the mass function of binary BHs (BBHs), challenging existing stellar evolutionary models with observations of BHs within the lower and upper mass gaps (e.g., Abbott et al. 2020a,b,c, 2024). Moreover, as the dataset grows, a possible evolution of BBH properties with redshift z is beginning to surface, even with current observations limited to z ≲ 1.5 (Biscoveanu et al. 2022; Abbott et al. 2023b; Nitz et al. 2023; Ray et al. 2023; Callister & Farr 2024; Rinaldi et al. 2024). Next-generation GW detectors will further advance our ability to explore the Universe, allowing us to probe BBH mergers at redshifts up to z ≲ 100 (Hall & Evans 2019; Kalogera et al. 2021; Branchesi et al. 2023). Specifically, the study of the BBH merger rates as a function of redshift can help us to test and distinguish among different formation channels of such objects (e.g., Belczynski et al. 2002; Dominik et al. 2013; Neijssel et al. 2019; Baibhav et al. 2019; Mapelli 2020; Mapelli et al. 2021; Zevin et al. 2021; Broekgaarden et al. 2021, 2022; Mandel & Broekgaarden 2022; Belczynski et al. 2022; van Son et al. 2022). The merger rate density is the result of the interplay between binary star evolution and the history of cosmic star formation (Lamberts et al. 2018; Santoliquido et al. 2020, 2022; Broekgaarden et al. 2021; Boesky et al. 2024a,b; Chruślińska 2024).

Most previous work relies on averaged star formation rate (SFR) density and metallicity distribution relations (e.g., van Son et al. 2022; Santoliquido et al. 2020; Neijssel et al. 2019). This assumption makes it extremely straightforward to calculate the merger rate density evolution, but contains several approximations in the description of the star formation process and metallicity evolution in the Universe. Such assumptions affect both the shape and absolute values of the estimated merger rates, as shown by several authors (Boco et al. 2019, 2021; Santoliquido et al. 2022; Briel et al. 2022). Current predictions for BBH merger rates span over three orders of magnitude (see e.g., Mandel & Broekgaarden 2022). However, emerging trends indicate that state-of-the-art binary population synthesis models tend to predict higher rates than those inferred by LVK, particularly when incorporating the latest observational data to model the Universe (e.g., Neijssel et al. 2019; Broekgaarden et al. 2021, 2022; Santoliquido et al. 2022; Srinivasan et al. 2023; Boesky et al. 2024a).

This discrepancy is largely driven by the interplay of merger efficiency (i.e., the number of BHs that merge per total stellar mass) and metal-poor star formation (Belczynski et al. 2010; Mapelli et al. 2019; Chruslinska & Nelemans 2019; Neijssel et al. 2019; Tang et al. 2020; Santoliquido et al. 2021; Broekgaarden et al. 2021; Chruślińska 2024). In fact, the merger efficiency of BBHs is 1−4 orders of magnitude higher in metal-poor than in metal-rich progenitor systems (Klencki et al. 2018; Giacobbo & Mapelli 2018; Spera et al. 2019; Broekgaarden et al. 2021; Iorio et al. 2023; van Son et al. 2025). Hence, the merger rate density evolution should be largely boosted by star formation in metal-poor regions, often associated with low-mass galaxies. As discussed by Chruślińska et al. (2021), Boco et al. (2021), Chruślińska (2024) the large uncertainties concerning the metallicity evolution in the Universe and the role of metal-poor galaxies might explain a fraction of the overall discrepancy.

Additionally, there is growing evidence that current models provide a poor description of the common envelope (CE) phase and binary evolution via stable mass transfer episodes might be more important than previously assumed for the formation of BBHs (Neijssel et al. 2019; Bavera et al. 2021; Marchant et al. 2021; Gallegos-Garcia et al. 2021; Olejak et al. 2021). Specifically, the stability criteria commonly adopted in population synthesis codes appear to produce unstable mass transfer for a wider range of parameters compared to detailed stellar modeling (Ge et al. 2010, 2015, 2020; Marchant et al. 2021). This in turn may lead to an overestimate of the BBH merger rates (Gallegos-Garcia et al. 2021).

This work is motivated by the considerations above. Firstly, we aim to assess the impact of the metal-dependent SFR on the BBH merger rates. To this end, we adopt and upgrade a detailed Universe model built from observational scaling relations (Boco et al. 2019, 2021) with GALAXYATE (Santoliquido et al. 2022). We test several SFR – galaxy mass relations, fundamental metallicity relations, and we explore the impact of low-mass galaxies. Secondly, we investigate the relative contribution of different formation channels to the BBH merger rates. We show that our state-of-the-art data-driven approach produces BBH merger rates consistently above the LVK inferred range. Such a discrepancy is robust against variations in the metal-dependent SFR model. In particular, low-mass galaxies (M* < 108 M) do not contribute more than a factor of ∼2 to the overall BBH merger rate density and their relative contribution strongly depends on the assumed SFR – galaxy mass relation.

The paper is organized as follows: Section 2 describes the methodology and the main codes used in this work. Section 3 provides an overview and a comparison of the adopted Universe models. Sections 4.2 and 4.3 present the results of the cosmic merger rate density as a function of the cosmic model. Section 4.4 details the merger rate density results for the different BBH formation channels. Finally, we discuss the implications of our results in Section 5. We summarize our findings in Section 6.

2. Method

We investigated the BBH merger rate density varying five model parameters, for a total of 360 parameter combinations. We used the code GALAXYATE (Santoliquido et al. 2022, 2023) to couple binary compact object populations to a synthetic Universe built from observational scaling relations (Mapelli et al. 2017; Artale et al. 2020; Santoliquido et al. 2021, 2022, 2023). The methodology adopted is described below.

2.1. Population-synthesis code SEVN

The stellar evolution for N-body code (SEVN)1 evolves stellar properties through on-the-fly interpolation of pre-computed stellar tracks (Spera & Mapelli 2017; Spera et al. 2019; Mapelli 2020; Iorio et al. 2023). In this work, we used the stellar tracks evolved with PARSEC (Bressan et al. 2012; Costa et al. 2019a,b; Nguyen et al. 2022). SEVN models binary evolution with semi-analytical prescriptions. We refer to Iorio et al. (2023) for a detailed description of the code.

For the purpose of this work, we adopted the default SEVN settings as described by Iorio et al. (2023). We assumed the DelayedGauNS supernova (SN) prescription, with BH masses determined according to the delayed model by Fryer et al. (2012). In our fiducial model, natal kicks are sampled as Giacobbo & Mapelli (2020, hereafter, GM20). Specifically, kick magnitudes are drawn from a Maxwellian distribution with one-dimensional root mean square σ = 265 km s−1 (Hobbs et al. 2005) and then re-scaled by a factor ∝Mej/Mrem, where Mej and Mrem are the mass of the ejecta and the compact remnant, respectively. We varied the kick prescription testing also a model in which BHs get the same natal kicks as derived for Galactic neutron stars (H05, Hobbs et al. 2005).

One of the biggest uncertainties in population-synthesis codes are the mass transfer stability criteria. The most common approach relies on the evaluation of the mass ratio q = Md/Ma and its comparison to a critical value qcrit, depending on the evolutionary phases of the stars involved. Here, Md is the mass of the donor star and Ma of the accretor. If q > qcrit the mass transfer episode is considered unstable on a dynamical timescale, triggering a CE (Ivanova et al. 2013). Our fiducial model assumes the SEVN default qcrit prescription (QCRS), that is the same as described by Hurley et al. (2002), but for one difference: we assume that mass transfer with donor stars in the main sequence and Hertzsprung gap is always stable (Iorio et al. 2023; Sgalletta et al. 2023).

SEVN handles CE evolution with the standard (αλ)−formalism (Webbink 1984; Hurley et al. 2002). We adopted the λ parameters specified by Claeys et al. (2014). Our fiducial model assumes a CE efficiency parameter α = 1 (this corresponds to assuming that the change in the orbital energy of the core contributes to unbind the envelope with an efficiency of 100%). However, we also explored different α parameters within our models, as detailed in Section 2.6.

2.2. Initial conditions

We sampled the masses of the primary star (M1) from a Kroupa initial mass function (Kroupa 2001), in a range between 5 and 150 M:

F ( M 1 ) M 1 2.3 . $$ \begin{aligned} \mathcal{F} (M_1) \propto M_{1}^{-2.3}. \end{aligned} $$(1)

We drew the secondary star mass (M2) within the range [2.2, 150] M from the observational distribution on q = M2/M1 (Sana et al. 2012):

F ( q ) q 0.1 , $$ \begin{aligned} \mathcal{F} (q) \propto q^{-0.1}, \end{aligned} $$(2)

with max ( 2.2 M M 1 , 0.1 ) q 1 $ \max \left(\frac{2.2\,\mathrm{M}_\odot}{M_{1}}, 0.1\right)\leq q \leq 1 $. The orbital periods and eccentricities were also drawn from the distributions derived by Sana et al. (2012):

F ( P orb ) ( log P orb ) 0.55 , $$ \begin{aligned} \mathcal{F} (P_{\rm orb}) \propto (\log P_{\rm orb})^{-0.55}, \end{aligned} $$(3)

with 0.15 ≤ log(Porb/d) ≤ 5.5, and

F ( e ) e 0.42 , $$ \begin{aligned} \mathcal{F} (e) \propto e^{-0.42}, \end{aligned} $$(4)

with 0 ≤ e ≤ 1 − (P/2 days)−2/3, following the correction by Moe & Di Stefano (2017). These distributions are fits to observations of a sample of massive binary stars in open clusters.

We sampled 107 binaries and used them as initial conditions for 12 metallicities: Z = 0.0002, 0.0004, 0.0008, 0.0012, 0.0016, 0.002, 0.004, 0.006, 0.008, 0.012, 0.016, and 0.02. The total simulated mass for each metallicity is thus Msim = 2.2 × 108 M.

2.3. Formation channels

We distinguish four main formation channels of BBHs, following the classification by Broekgaarden et al. (2022) and Iorio et al. (2023).

  • Channel I: The systems experience a stable mass transfer episode before the first compact object formation and only after they evolve through at least one CE phase.

  • Channel II: The systems interact only via stable mass transfer episodes.

  • Channel III: At the time of the first compact remnant formation the system has already experienced at least one CE phase and the system is composed of one H-rich star and one without the H envelope.

  • Channel IV: Similarly to Channel III, the systems experience at least one CE before the first compact remnant formation, but in this case, at the time of compact remnant formation both stars have lost their H envelopes.

2.4. Observational scaling relations

In order to model the Universe and populate it with BBHs, we adopted the code GALAXYATE2 (Santoliquido et al. 2022). We generated a set of star-forming galaxies from observational scaling relations for an array of formation redshifts, sampling their masses M*, SFRs, and average metallicities. Here below, we summarize such scaling relations.

We adopted the star-forming galaxy stellar mass function derived by Chruslinska & Nelemans (2019). For each formation redshift, we kept the comoving volume fixed V ∼ (100 cMpc)3 and sampled a number of star-forming galaxies 𝒩gal(zform) proportional to the galaxy number density. We sampled galaxy masses in the range [Mmin,1012 M]. The minimum galaxy mass is a free parameter in our models, we tested different values: Mmin = 106, 107, and 108 M. We assigned a SFR to every sampled galaxy following a double lognormal distribution, with a first peak centered on the galaxy main sequence, and a second one, on the starbust sequence (Boco et al. 2021). The parameters adopted are described in Appendix A.

The slope of the galaxy main sequence is highly debated, especially moving toward low-mass galaxies (≲ 108.5 M) and high redshifts. For this reason, we tested several galaxy main sequence definitions. Specifically, we considered the definitions given by Speagle et al. (2014, hereafter S14), Boogaard et al. (2018, hereafter B18), and Popesso et al. (2023, hereafter P23). In the following equations, the SFR is expressed in M yr−1, and all masses are given in solar units. S14 define the galaxy main sequence as:

log SFR ( M , t ) = ( 0.84 0.026 t ) log M 6.51 0.11 t , $$ \begin{aligned} \log \mathrm{SFR} (M_*, t) = (0.84 - 0.026 t)\log M_*- 6.51 - 0.11 t, \end{aligned} $$(5)

where t is the age of the Universe in Gyr. The definition by B18 is as follows:

log SFR ( M , z ) = 0.83 log ( M M 0 ) 0.83 + 1.74 log ( 1 + z 1 + z 0 ) , $$ \begin{aligned} \log \mathrm{SFR} (M_*, z) = 0.83 \log \left(\frac{M_*}{M_0}\right) - 0.83 + 1.74 \log \left(\frac{1+z}{1+ z_0}\right), \end{aligned} $$(6)

where M0 = 108.5 M and z0 = 0.55 are the median values of the data. The fit by P23 is:

log SFR ( M , t ) = a 0 + a 1 t log [ 1 + ( M / 10 a 2 + a 3 t ) a 4 ] , $$ \begin{aligned} \log \mathrm{SFR} (M_*, t) = a_0 + a_1 t - \log \left[1 + \left(M_*/10^{a_2 + a_3 t}\right)^{-a_4}\right], \end{aligned} $$(7)

with a0 = 2.693 ± 0.012, a1 = −0.186 ± 0.009, a2 = 10.85 ± 0.05, a3 = −0.0729 ± 0.0024 and a4 = 0.99 ± 0.01.

Moreover, we built two additional parametric models that explore different slopes of the main sequence relation at low-masses (M* < 109 M). In this range, in fact, the data is scarce (see e.g., Popesso et al. 2023) and very uncertain. With this approach, we were able to bracket all possible behaviors of the galaxy main sequence at low-masses. For these models, we cut the fit in P23 at 109 M, and adopted a constant slope a across cosmic times for lower galaxy masses:

log SFR ( M , t ) = { a log M + b for M < 10 9 M Eq . 7 for M 10 9 M . $$ \begin{aligned} \log \mathrm{SFR}(M_*, t) = {\left\{ \begin{array}{ll} a \log M_*+ b&\mathrm{for} \; M_{*} < 10^9\,\mathrm{M}_{\odot } \\ \mathrm{Eq.}~7&\mathrm{for} \; M_{*} \ge 10^9\,\mathrm{M}_{\odot }. \end{array}\right.} \end{aligned} $$(8)

We assumed a = 0.5, for the SHALLOW model, and a = 1 for the STEEP model. We decided such values for the SHALLOW (STEEP) slope to always keep it as the shallowest (steepest) main sequence slope at all redshifts in the low-mass galaxies’ tail. The continuity assumption between the two equations sets the value of the second parameter b. See Appendix A.2 for more details on the main sequence definitions adopted.

By integrating the galaxies SFR for each redshift bin, we get the SFR density as a function of redshift, shown in Fig. 1. We compare the predicted SFR density for all the adopted galaxy models. The model by B18 predicts a lower SFR density with respect to the observational data. This is likely a result of the data sample utilized in B18 reaching up to z ∼ 0.9, wheareas here we extrapolate their relation to much higher redshifts. All the other galaxy main sequence models predict SFR densities that are consistent with UV and far-IR data (see Madau & Dickinson 2014, and references therein). The fit by S14 adopts galaxy main sequence data, with galaxy stellar masses within 109.5 < M* < 1011.5 M. P23 expand the work by S14, by incorporating the most updated and recent galaxy main sequence data, encompassing the widest range of galaxy stellar mass (108.5 < M* < 1011.5 M) and redshift (0 < z < 6) available to date. Thus, we adopted Equation (7) P23 as our fiducial SFR–M* relation. Moreover, varying the lower mass limit of the galaxy main sequence (represented by different line styles) does not significantly impact the SFR density. This is not surprising, as the SFR from low-mass galaxies is suppressed by the galaxy main sequence relation. Therefore, even if these small galaxies are more numerous, the contribution to the total SFR density is negligible compared to the contribution from higher mass galaxies (M* ≥ 108M).

thumbnail Fig. 1.

SFR density as a function of redshift for different SFR–M* relations (different colors). The different linestyles show the results by cutting the galaxy main sequence at different minimum galaxy masses, Mmin = 106 M (dotted), 107 M (solid) and 108 M (dashed). The wide blue line shows the fit by Madau & Fragos (2017). The markers show observed data points: the blue dots are are taken from Madau & Dickinson (2014), the orange diamonds are from Casey et al. (2018), and the purple squares are taken from Gruppioni et al. (2020).

For each galaxy, we sampled an average metallicity based on the fundamental metallicity relation, that links the average metallicity of a galaxy with its stellar mass and SFR (Mannucci et al. 2010; Maiolino & Mannucci 2019). In this work, we compared three different definitions of the fundamental metallicity relation, corresponding to the fits by Mannucci et al. (2011), Curti et al. (2020), and Andrews & Martini (2013). The recent work by Nakajima et al. (2023) compares known models of mass–metallicity relations with a comprehensive, up-to-date sample of James Webb Space Telescope (JWST) observations of galaxies at redshift 4 to 10. Their data (see their Figure 12) show the fundamental metallicity relation by Andrews & Martini (2013) to best fit these new data up to z ∼ 8. For this reason, we adopted the fundamental metallicity relation by Andrews & Martini (2013) as our fiducial model. We assumed that the metallicity is distributed following a lognormal function within each galaxy, with a scatter σZ = 0.14. Such an observed scatter (amounting to ≲0.1 dex in the local Universe and ≲0.2 dex at high redshift) is small since it refers to galaxies at a given redshift with the same stellar mass and SFR. However, the metallicity at a given z is substantially more dispersed over the whole galaxy populations; for instance, averaging the fundamental metallicity relation over stellar mass (via the stellar mass function) yields a mass-metallicity relation characterized by a large and asymmetric dispersion of ≳0.35 dex, in agreement with observational findings (Chruślińska et al. 2021; Boco et al. 2021). We describe the adopted observational scaling relations in more detail in Appendix A and review all the models’ parameters in Section 2.6.

We also tested a different approach based on the average evolution of the cosmic SFR and metallicity with COSMOATE (Santoliquido et al. 2020). In this case, the cosmic SFR density ψ(z) and the average metallicity evolution μ(z) = ⟨Z/Z⟩ are modeled using the fitting formulae provided by Madau & Fragos (2017):

ψ ( z ) = 0.01 ( 1 + z ) 2.6 1 + [ ( 1 + z ) / 3.2 ] 6.2 M Mpc 3 yr 1 , $$ \begin{aligned} \psi (z)&= 0.01 \frac{(1+z)^{2.6}}{1+ [(1+z)/3.2]^{6.2}}\,\mathrm{M}_{\odot }\,\mathrm{Mpc}^{-3}\,\mathrm{yr}^{-1}, \end{aligned} $$(9)

log μ ( z ) = 0.153 0.074 z 1.34 . $$ \begin{aligned} \log \mu (z)&= 0.153 - 0.074 z^{1.34}. \end{aligned} $$(10)

We spread the metallicities around this average value following a normal distribution with standard deviation σ Z $ \tilde{\sigma}_{\mathrm{Z}} $:

p ( Z | z ) = 1 2 π σ Z 2 exp [ [ log ( Z / Z ) log ( Z / Z ) ] 2 2 σ Z 2 ] , $$ \begin{aligned} p(Z \vert z) = \frac{1}{\sqrt{2 \pi \tilde{\sigma }_{\rm Z}^2}} \exp \left[-\frac{[\log (Z/Z_{\odot }) - \langle \log (Z/Z_{\odot }) \rangle ]^{2}}{2 \tilde{\sigma }_{\rm Z}^{2}}\right], \end{aligned} $$(11)

where log Z / Z = μ ( z ) ln 10 σ Z 2 / 2 $ \langle \log{Z/Z_{\odot}}\rangle = \mu(z) - \ln{10} \tilde{\sigma}_{\mathrm{Z}}^{2}/2 $. With this simplified approach, we can change the spread of metallicity at a given redshift just by tuning the parameter σ Z $ \tilde{\sigma}_Z $.

2.5. Merger rate density

We evaluate the merger rate density with GALAXYATE as Santoliquido et al. (2022):

R ( z ) = 1 V 3 z max z [ i = 1 N gal ( z ) Z min Z max S i ( z , Z ) F ( z , z , Z ) d Z ] d t ( z ) d z d z , $$ \begin{aligned} \mathcal{R} (z) = \frac{1}{V^3} \int _{z_{\rm max}}^{z} \left[\sum _{i = 1}^{\mathcal{N} _{\rm gal}(z\prime )} \int _{Z_{\rm min}}^{Z_{\rm max}} \mathcal{S} _{i} (z\prime , Z) \mathcal{F} (z\prime , z, Z) \mathrm{d}Z\right] \frac{\mathrm{d}t(z\prime )}{\mathrm{d}z\prime } \mathrm{d}z\prime , \end{aligned} $$(12)

where we sum over all of the formation galaxies, 𝒩gal(z′), at redshift z′ (i.e., the galaxies where the binary systems form at z′≥z), and then integrate over the redshift between the BBH merger redshift z (i.e., the redshift at which the BBH merges) and the maximum considered formation redshift zmax = 8. For the i-th galaxy 𝒮i(z′,Z) = ψi(z′)pi(Z|z′), where ψi(z′) is its SFR and pi(Z|z′) is its metallicity distribution at a given redshift z′. The second term in the integral accounts for the distribution of binaries in our catalogs:

F ( z , z , Z ) = 1 M sim N BBH ( z , z , Z ) d t f bin f corr , $$ \begin{aligned} \mathcal{F} (z\prime , z, Z) = \frac{1}{M_{\rm sim}} \frac{\mathcal{N} _{\rm BBH} (z\prime , z, Z)}{\mathrm{d}t} f_{\rm bin} f_{\rm corr}, \end{aligned} $$(13)

where 𝒩BBH(z′,z, Z)/dt is the rate of BBHs progenitors that form at redshift z′ and merge at redshift z, for a given metallicity Z. The factor fbin = 0.5 corrects for the fraction of binaries, and fcorr = 0.251 takes into account the cut at lower masses from our initial sampling conditions (see Section 2.2).

2.6. Summary of model parameters

We ran our SEVN and GALAXYATE simulations over an extensive grid of model parameters, exploring both binary evolution and galaxy modeling uncertainties. We summarize the grid of models in Table 1. Specifically, we ran sets of simulations with SEVN varying the CE efficiency parameter α = 0.5, 1, 3, and 5. We consider two natal-kick models: GM20, and H05, as defined in Section 2.1.

Table 1.

Summary of the grid of parameters varied in our study.

As for the synthetic galaxy models, we vary the minimum galaxy mass assuming Mmin = 106, 107 and 108 M. We explore the galaxy main-sequence definitions from S14, B18, P23, as well as our SHALLOW and STEEP models (Equation 8 with a = 0.5 and a = 1, respectively). For the fundamental metallicity relations, we compare the fits by Mannucci et al. (2011), Andrews & Martini (2013), and Curti et al. (2020). Our fiducial model incorporates the most updated main sequence definition using P23, the fundamental metallicity relation by Andrews & Martini (2013), and a minimum galaxy stellar mass Mmin = 107 M.

3. A comparative analysis of Universe models

Figure 2 provides a schematic overview and comparison of the two approaches to modeling a synthetic Universe. The detailed method offers several improvements over the conventional averaged approach. We note that σ Z $ \tilde{\sigma}_{\mathrm{Z}} $, used in the averaged model, and σZ, used in the detailed model, have two different physical meanings. σZ represents the metallicity dispersion around the FMR (i.e., at fixed galaxy stellar mass and SFR), that is then convoluted with the distribution of galaxy masses and SFR to obtain the spread in metallicity at a given redshift. Therefore, even if we assume the distribution around the FMR to be log-normal and with a small scatter ∼0.14, the overall distribution of metallicities is asymmetric and can have a large scatter, driven by the observed galaxy distribution at that specific redshift.

thumbnail Fig. 2.

Schematic overview of the different metal-dependent SFR models adopted here. Upper panel: overview of the detailed Universe model, based on observational scaling relations. Lower panel: overview of the model based on the average SFR and metallicity evolution across cosmic time. The plots on the right illustrate the metallicity distributions as a function of redshift for our different prescriptions. For the average metallicity evolution model, we show the distributions assuming both σ Z = 0.1 $ \tilde{\sigma}_{\mathrm{Z}} = 0.1 $ and 0.6.

In contrast, models relying on averaged metallicity evolution require the total scatter ( σ Z $ \tilde{\sigma}_{\mathrm{Z}} $) to artificially encompass the full metallicity distribution at a given redshift. For this reason typical σ Z $ \tilde{\sigma}_{\mathrm{Z}} $ values of ∼0.5 are required. In this work, we adopt σ Z = 0.6 $ \tilde{\sigma}_{\mathrm{Z}}=0.6 $ as fiducial value and σ Z = 0.1 $ \tilde{\sigma}_{\mathrm{Z}}=0.1 $, even if not realistic, to demonstrate the effect of reducing the low metallicity tail. Moreover, the distribution adopted is symmetric around the average metallicity. Therefore, even if the spread in the averaged model is large enough to mimic the amount of SFR at low-metallicities of the detailed approach, it inevitably produces an unrealistic high-metallicity tail, that is not supported by observations (see e.g., Chruślińska 2024).

We note that the distribution of metallicities of the averaged models could be improved by adopting a skewed distribution around the mean value (see e.g., van Son et al. 2023; Fishbach 2025). Such models allow for a more realistic representation of metallicities, when properly fitted to the data. We do not consider them here as they would ultimately lead to results consistent with our detailed approach.

4. Results

4.1. Merger efficiency and metallicity

The BBH merger efficiency (i.e., the number of BHs that merge per total initial stellar mass) depends on metallicity, being higher for metal-poor stars (Dominik et al. 2012; Giacobbo & Mapelli 2018; Klencki et al. 2018). With SEVN, we find an abrupt drop of three orders of magnitude in the BBH merger efficiency, above ∼1/3 Z (Fig. 3). The drop mildly depends on the α parameter: it is close to 1/3 Z for α ≤ 1 and shifts to lower metallicities (∼1/10 Z) for α = 3 or 5. For this reason, the fraction of SFR at low metallicity is crucial for determining the BBH merger rate.

thumbnail Fig. 3.

Merger efficiency η for BBHs as a function of metallicity Z. From top to bottom α = 0.5, 1, 3, 5. The different colors refer to formation channels I–IV (Section 2.3). The vertical red lines highlight metallicity values of Z = 1/100, 1/10, and 1/3 Z.

4.2. The impact of low-mass galaxies

Figure 4 shows the relative contribution of different galaxy masses to the total SFR density and to the SFR density below a certain Z threshold. We consider two reference metallicity thresholds: Z ≤ 1/3 Z and ≤1/10 Z. We compare the results for different main-sequence definitions. In Figure 4, we emphasize the SHALLOW and STEEP models (wider lines), because they represent the extreme cases: the former has the shallowest trend, whereas the latter has the steepest decrease in the SFR density at low galaxy masses. This figure shows the results for galaxies at z = 0; however, the following considerations hold for every considered redshift.

thumbnail Fig. 4.

Contribution of galaxies of mass M* to the SFR density evaluated for different metallicity thresholds at z = 0. The x axis represents the galaxy stellar mass M*, while the lines show the total contribution to the SFR density from all galaxies in a given stellar mass bin, with a metallicity Z ≤ Zth. The different line styles refer to different metallicity thresholds Zth; solid lines: Zth = 1/3 Z, dashed: Zth = 1/10 Z, and dotted: Zth = 1/100 Z. The different colors represent models with different SFR–M* relations, blue: B18, violet: S14, pink: P23, orange: SHALLOW, yellow: STEEP. Each panel shows the results for a different FMR, from left to right: Mannucci et al. (2011), Andrews & Martini (2013), Curti et al. (2020).

We see from Fig. 4 that only with the shallowest SFR–M* relations the contribution from low-mass galaxies is important, for metallicities Z ≤ 1/3 Z. For the steepest galaxy main-sequence models (e.g., the STEEP model and P23), most of the SFR for metallicities below 1/3 Z originates from galaxies with masses in the range of (109 − 1010) M for the fundamental metallicity relations by Andrews & Martini (2013) and Curti et al. (2020), while the contribution from 106 M galaxies is about half. The peak in SFR shifts toward galaxy masses of (108 − 109) M adopting the models by Mannucci et al. (2011). On the other end, shallower main sequence slopes (e.g., the SHALLOW model, S14) show a similar SFR contribution from low mass galaxies up to ∼108 M. Model B18 falls in between these two extreme cases.

If we restrict our attention to even lower metallicities (Z < 1/10 Z), the contribution from low-mass galaxies becomes important for all the SFR–M* relations. As a result, for higher α parameters (for which the merging efficiency peaks at smaller metallicities, see Fig. 3) the contribution from low-mass galaxies becomes important, regardless of the chosen SFR–M* relation. While the trend described above is true for all the fundamental metallicity relations, we see from the different panels that the prescriptions by Curti et al. (2020) and Andrews & Martini (2013) predict on average higher metallicities (Fig. A.1).

Several studies have shown how the contribution of starburst galaxies might increase at high redshifts compared to our current treatment. Both Bisigello et al. (2018) and Rinaldi et al. (2022) find that starburst galaxies might contribute to ≳50% of the cosmic SFR at redshifts z ∼ 4 − 5. Bisigello et al. (2018) finds that starburst galaxies might constitute the ∼16% of total galaxies at redshifts 2 < z < 3. A higher percentage of starburst galaxies at high redshift would have two main effects. First, we would see an increase in the global SFR, as starburst galaxies have an higher SFR at fixed mass. Second, starburst galaxies specifically increase the SFR at low metallicity. In fact, because starburst galaxies have a higher SFR at fixed mass, they typically exhibit lower metallicity due to the FMR. Both factors contribute to an increase in the BBH merger rates, although the impact at z = 0 is smaller, as this effect primarily affects redshifts z ≳ 2 − 3.

4.3. Cosmic merger rate density

Figure 5 shows the BBH cosmic merger rate density we obtain assuming different SFR–M* relations, adopting the fundamental metallicity plane by Andrews & Martini (2013). The results for the fundamental metallicity relations by Mannucci et al. (2011) and Curti et al. (2020) are shown in Appendix B. We compare our merger rate densities with the inferred fit ℛ(z)∝(1 + z)k by the LVK collaboration up to the third observing run (Abbott et al. 2023b).

thumbnail Fig. 5.

BBH merger rate density as a function of redshift for different observational scaling relations (different colors) and CE α parameters (different panels). Each color is associated with a different SFR–M* relation, as summarized in the legend in the top left-hand panel. The solid lines show the model with Mmin = 107 M. The shaded area encompasses the difference between Mmin = 106 M (upper boundary), and Mmin = 108 M (lower boundary). All the models assume the fiducial fundamental metallicity relation by Andrews & Martini (2013). The thick dashed lines show the merger rate density obtained with COSMOATE (simplified model), assuming the cosmic SFR density by Madau & Dickinson (2014) and a lognormal distribution for the metallicity with a spread σ Z = 0.1 $ \tilde{\sigma}_{\mathrm{Z}}=0.1 $ (dark gray) and 0.6 (light gray). The gray-shaded regions show the MRD inferred by the LVK collaboration: ℛ(z)∝(1 + z)k (Abbott et al. 2023b).

All predicted BBH merger rate densities exceed the LVK estimates across all our galaxy-derived models. The B18 model yields the lowest merger rate density. This is expected as B18 results in the lowest SFR density (see Figure 1). As discussed in Section 2.4, this is likely due to our extrapolation of the relation in B18 to high redshifts, whereas the data sample considered there is limited to z ≲ 0.9.

The merger rate densities evaluated by setting different minimum galaxy masses differ at most by a factor of ∼2. Even removing all the galaxies with mass below 108 M (lower bound of the shaded areas in Fig. 5), the merger rate density is still well above the 90% credible interval from LVK. The magnitude of the difference between assuming Mmin = 106, 107 and 108 M depends on the main sequence definition adopted. As expected, steeper SFR–M* relations result in smaller differences with the minimum galaxy stellar mass adopted. The differences between Mmin = 106, 107, and 108 M are bigger for α = 3 and α = 5, compared to α = 1 or α = 0.5, as expected from the results described in Section 4.2.

Figure 5 also compares the merger rate density we obtain through our host-galaxy models (i.e., with GALAXYATE) with the merger-rate density calculated adopting an average SFR density and metallicity evolution of the Universe (i.e., with COSMOATE). The merger rate density we obtain modeling the host galaxies with observational scaling relations and the one we estimate by assuming an average SFR density agree only when we assume a large metallicity spread ( σ Z = 0.6 $ \tilde{\sigma}_{\mathrm{Z}}=0.6 $, see Sect. 2.4) for the latter. Moreover, even assuming a large value of σ Z $ \tilde{\sigma}_{\mathrm{Z}} $, the predicted shape and steepness of the merger rate density differ in the two models, especially at high redshift and for models with α > 1.

4.4. Formation channels

Figure 6 shows the merger rate density for each BBH formation channel, as defined in Sect. 2.3, through Equation (12). Channel I (traditional CE scenario, with at least one CE phase after the formation of the first BH) gives the main contribution to the full BBH merger rate density at all redshifts, if α ≤ 1. Channel I is subdominant only at high redshift and when α > 1. This is the result of longer delay time distributions for channel I, compared to channels II and IV. As a consequence, even if the merger efficiency of channel I is higher at all metallicities (Figure 3), channels II and IV dominate the merger rate density at high redshift. Moreover, the longer delay times result in a much steeper shape of the merger rate density for channel I, especially if we adopt a large value of α.

thumbnail Fig. 6.

BBH merger rate density as a function of redshift for different BBH formation channels. The solid black lines show the merger rate density of the whole BBH population. The shaded area encompasses the difference between Mmin = 106 M (upper boundary), and Mmin = 108 M (lower boundary). The different panels show the results obtained assuming different CE α parameters. All the models assume the fiducial SFR–M* relation P23 and the fundamental metallicity relation by Andrews & Martini (2013). The gray-shaded regions show the MRD inferred by the LVK collaboration (Abbott et al. 2023b).

At low redshift, the merger-rate density associated with channel I is always largely above the LVK 90% credible interval. In contrast, the merger rate densities associated with channels II (stable mass transfer), III and IV lie within the observed range reported by LVK for all the considered α parameters. Channel III is by far the least common formation pathway.

5. Discussion

Predictions of BBH merger rate are highly uncertain, as they strongly depend on binary evolution physics and the metal-dependent SFR. Mandel & Broekgaarden (2022) show that predicted rates can vary by more than three orders of magnitude depending on the adopted assumptions. Binary population synthesis codes generally tend to overpredict BBH merger rates, suggesting that careful calibration of binary parameters is necessary to bridge the gap with observations. This trend is evident across multiple studies using different population synthesis codes, including STARTRACK (Klencki et al. 2018; Chruslinska et al. 2019; Olejak et al. 2021, 2022; Romagnolo et al. 2023), COMPAS (Neijssel et al. 2019; Broekgaarden et al. 2021, 2022; Riley et al. 2021; Stevenson & Clarke 2022; Rauf et al. 2023; Boesky et al. 2024a), COSMIC (Zevin et al. 2020; Srinivasan et al. 2023), BPASS (Eldridge et al. 2018; Tang et al. 2020), and POSYDON (Bavera et al. 2021; Román-Garza et al. 2021). For similar variations in binary parameters, Boesky et al. (2024a) find rates that are consistent with our results, adopting the COMPAS binary population synthesis code (Team COMPAS 2022). Several studies have explored the contribution of different formation channels to the BBH population. Bavera et al. (2021) and Román-Garza et al. (2021) identify the evolution through at least one CE as the dominant channel in most models. The model reported by Olejak et al. (2021) predicts a higher contribution to the BBH merger rates from channel I; however, their study highlights that different treatments of mass transfer stability can significantly alter the relative contributions of each channel. Similarly, Dorozsmai & Toonen (2024) find the stability of mass transfer to be a key factor in enhancing the contribution from the stable mass transfer channel. Notably, in the simulations run with COMPAS (Neijssel et al. 2019; Broekgaarden et al. 2021) channel II is dominant compared to channel I because of the different stability criteria. Finally, the choice of the metal-dependent SFR can have a large impact on the BBH merger rates, as shown by Neijssel et al. (2019), Broekgaarden et al. (2021, 2022), Srinivasan et al. (2023). Specifically, Tang et al. (2020) show that a narrower metallicity distribution leads to lower merger rates, in agreement with our findings.

In this work, we have modeled the host galaxies of BBHs, their SFR and metallicity with the most up-to-date observational scaling relations. Nevertheless, as we improve the description of the metal-dependent SFR, the BBH merger rate density we obtain at low redshift (102 − 103 Gpc−3 yr−1) is always well above the 90% credible interval inferred by LVK after the third observing run (Fig. 5). This is true for all of our models, considering different SFR–M* relations, different fundamental metallicity relations, and different values of the efficiency parameter for CE ejection.

As discussed above, we can exclude that the overestimate of the BBH merger rate density is a feature of our population-synthesis code SEVN (Iorio et al. 2023). Rather, it is a problem common to current population synthesis codes.

In the following, we discuss several potential solutions to reconcile theoretical predictions and GW data. One possible explanation is that the CE phase occurs less frequently than previously assumed, or that BBH progenitors do not survive the CE phase as often as currently believed. Conversely, the stable mass transfer (channel II) may play a considerably more significant role in BBH formation than previously expected, as already proposed by several authors (Pavlovskii et al. 2017; Neijssel et al. 2019; Misra et al. 2020; Klencki et al. 2020, 2021; Olejak et al. 2021; Shao & Li 2021; Bavera et al. 2021; Marchant et al. 2021; van Son et al. 2022). In binary population synthesis codes, the CE channel tends to dominate over the others; in contrast, stable mass transfer is prevalent in detailed stellar-structure calculations (e.g., Gallegos-Garcia et al. 2021). One of the main factors of uncertainty lies in the criteria for the onset of a CE phase. As proposed by Ge et al. (2020), mass transfer might be stable for a wider range of mass ratios than previously believed, at least under the assumption of adiabatic, conservative mass transfer. Another slightly different interpretation is that the traditional description of semi-major axis evolution during a CE phase (based on the α parameter, Hurley et al. 2002) is too simplified, and does not capture the actual evolution of the system (Hirai & Mandel 2022; Trani et al. 2022; Di Stefano et al. 2023; Everson et al. 2025). If the orbital separation during a CE phase shrinks less than predicted by the traditional α formalism, then the resulting merger rate density is likely going to be lower (e.g., Hirai & Mandel 2022; Lau et al. 2022).

Natal kicks also significantly influence the merger rate density (Mapelli & Giacobbo 2018; Iorio et al. 2023; Boesky et al. 2024a). In Fig. 7, we compare our fiducial model with a model in which natal kicks are drawn from a Maxwellian distribution with root-mean square velocity σ = 265 km s−1, without accounting for the fallback (Hobbs et al. 2005). This model yields substantially higher kicks, leading to lower BBH merger rate densities across all considered α parameters. Thus, another possible way to reconcile observed and predicted merger-rate density is to assume that natal kicks are higher than usually expected for BBHs. However, this also has an impact on the BH mass function, favoring the merger of more massive binary systems (Iorio et al. 2023). The recent work by Nagarajan & El-Badry (2025) reveal the existence of both strong (≳100 km s−1) and weak (≲50 km s−1) black hole kicks, potentially indicating a bimodal distribution (see also Verbunt et al. 2017). Although the uncertainties remain too large to conclusively confirm this scenario, these preliminary results suggest that the black hole kick distribution may lie between the two extreme models proposed in our work.

thumbnail Fig. 7.

Merger rate density as a function of redshift, for the natal-kick models GM20 (solid lines) and H05 (dashed lines). The line colors show the results obtained assuming different CE efficiency parameters (α). The upper boundaries of the shaded areas show the models obtained with minimum galaxy mass M* = 106 M, whereas for the lower boundary the minimum mass is 108 M. All the models assume the fiducial observational scaling relations. The gray-shaded region shows the MRD inferred by the LVK collaboration (Abbott et al. 2023b).

Finally, recent work by Schneider et al. (2021) shows that envelope stripping can affect the final fate and compact remnant mass of massive stars, increasing the minimum zero-age main sequence mass of BH progenitors from ∼30 M to ∼70 M. While population calculations still need to be done for this scenario, this might result in a decrease in the BBH merger rate density as well.

A few caveats should be added. Recent work by Rinaldi et al. (2022) suggests that starburst galaxies contribute to 60−90% of the total SFR. Rinaldi et al. (2022) find a clear bimodality between the galaxy main sequence and the starburst sequence for higher massive galaxies, M* ≥ 109 M. Our current approach likely underestimates the impact of starburst galaxies (see also Caputi et al. 2017; Bisigello et al. 2018). Improved modeling of starburst galaxies is beyond the scope of this work. Nevertheless, an higher contribution of starburst galaxies would further increase the predicted BBH merger rates as they tend to be more metal-poor than main sequence galaxies. We would like to remark that models based on the evolution of the average SFR and metallicity over cosmic time (e.g., Santoliquido et al. 2021) broadly capture our fiducial BBH merger rate densities (Fig. 5), when accounting for a wide spread in metallicity for a given redshift bin ( σ Z = 0.6 $ \tilde{\sigma}_{\mathrm{Z}} = 0.6 $). This assumption is needed to mimic the complex distribution of metallicities that in our detailed approach naturally results from the spread in galaxy masses and SFRs. We do not vary the qcrit prescriptions in this study. However, Figure 22 in Iorio et al. (2023) demonstrates that different qcrit prescriptions affect the BBH merger rate density by a factor of ∼2. Therefore, given the prescriptions available to date, we do not expect this parameter to significantly influence our results. Nevertheless, different conditions for mass-transfer stability may have an impact on the relative contribution from different formation channels (see e.g., Olejak et al. 2021, 2022; Dorozsmai & Toonen 2024). The parameter fbin, takes into account the fraction of binary systems in the Universe. This value strongly varies with stellar mass, being close to 1 for massive stars (M1 ≥ 16 M, see e.g., Sana et al. 2012; Moe & Di Stefano 2017) and decreasing for lower primary masses, reaching ∼0.4 for M1 ∼ 2 M (Moe & Di Stefano 2017). We adopt a constant fbin = 0.5 as a reference value, to account for the whole population of binary systems. It should be noted that fbin only acts as a renormalization factor, rigidly shifting the merger rates upward or downward. Adopting a different fbin value in the range [0.4, 1] would lead to at most a factor of 2 difference in the merger rates.

As a final remark, differences between the BBH mass function obtained from population-synthesis calculations and the one assumed to infer the LVK rates, might bias the comparison between the two merger rate densities. However, the discrepancy between predictions and LVK rates (at least a factor of ∼10) is too large to be fully accounted for by such BBH mass function differences.

Finally, here we have only considered isolated binary evolution. Including the dynamical formation of BBHs in star clusters and galactic nuclei would only lead to a further increase of the merger rate density because dynamics is known to boost the merger rate density even in a metal-rich environment (Rodriguez & Loeb 2018; Di Carlo et al. 2020; Mapelli et al. 2022; Barber & Antonini 2025).

6. Summary

We have studied the merger rate density of BBHs across cosmic time over 360 combinations of parameters, exploring the effects of different galaxy observational scaling relations and binary evolution models. Furthermore, we have studied the relative impact of different isolated BBH formation channels on the merger rate density. We can summarize our main results as follows.

  • Current models of binary population synthesis predict values of the BBH merger rate density at redshift z ≤ 1 that are well above the 90% credible interval inferred from LVK data (Abbott et al. 2023b).

  • The impact of low-mass galaxies (M* < 108 M) on the BBH merger rate density is highly dependent on the steepness of the SFR–M* relation. For steep SFR–M* relations P23, galaxies with mass M* ≲ 107 M do not contribute significantly to the merger rate density.

  • Even removing all the galaxies with mass below 108 M (hence, with lower metallicity), the merger rate density is still well above the 90% credible interval from LVK. Thus, uncertainties in the number of low-mass galaxies do not dramatically affect the merger rate density and do not solve the tension between predicted and observed BBH merger rate density.

  • Models based on the average SFR and metallicity evolution need to account for the whole metallicity spread at a given redshift σ Z 0.6 $ \tilde{\sigma}_{\mathrm{Z}} \sim 0.6 $ in order to mirror detailed Universe prescriptions. In this case, the merger rate densities predicted by the two approaches broadly agree, with differences that become more noticeable toward high redshifts.

  • The merger rate density of BBHs is dominated by Channel I, where the systems evolve through at least one CE phase only after the first compact object formation.

Overall, our results clearly indicate that there is a tension between current models of BBH merger rates and values inferred from GW data. We cannot explain this tension through uncertainties on the cosmic SFR density, as we have adopted state-of-the-art models, robustly grounded in observations. Future studies are requested to understand the origin of this tension, addressing –e.g.,– the nature of core collapse, the stability of mass transfer, and the magnitude of natal kicks.


1

In this work, we used the SEVN version V 2.10.1 (commit 22c9236). SEVN is publicly available at the gitlab repository https://gitlab.com/sevncodes/sevn

2

The code GALAXYATE can be found at https://gitlab.com/Filippo.santoliquido/galaxy_rate_open

Acknowledgments

The authors are grateful to Dylan Nelson and Annalisa Pillepich for their helpful insights. The authors thank the anonymous referee for the constructive report. MM acknowledges financial support from the European Research Council for the ERC Consolidator grant DEMOBLACK, under contract no. 770017. MM also acknowledges financial support from the German Excellence Strategy via the Heidelberg Cluster of Excellence (EXC 2181 – 390900948) STRUCTURES. GI also received the support of a fellowship from the “la Caixa Foundation (ID 100010434)”. The fellowship code is LCF/BQ/PI24/12040020. GI also acknowledges financial support under the National Recovery and Resilience Plan (NRRP), Mission 4, Component 2, Investment 1.4, – Call for tender No. 3138 of 18/12/2021 of Italian Ministry of University and Research funded by the European Union – NextGenerationEU. MCA acknowledges financial support from FONDECYT Iniciación 11240540. and ANID BASAL project FB210003. The authors acknowledge support by the state of Baden-Württemberg through bwHPC and the German Research Foundation (DFG) through grants INST 35/1597-1 FUGG and INST 35/1503-1 FUGG. Software. We use the SEVN version V 2.10.1 (commit 22c9236) to generate our BBHs catalogs. SEVN is publicly available at the gitlab repository https://gitlab.com/sevncodes/sevn. We use TRACKCRUNCHER (https://gitlab.com/sevncodes/trackcruncher) (Iorio et al. 2023) to produce the tables needed for the interpolation in SEVN. We use the code GALAXYATE to evolve our BBHs in a Universe model. GALAXYATE can be found at https://gitlab.com/Filippo.santoliquido/galaxy_rate_open. We use the code COSMOATE to model the averaged Universe models. COSMOATE can be found at https://gitlab.com/Filippo.santoliquido/cosmo_rate_public. This research made use of NUMPY (Harris et al. 2020), SCIPY (Virtanen et al. 2020), PANDAS (The Pandas Development Team 2024) and ASTROPY (Astropy Collaboration 2013, 2018, 2022). For the plots we used MATPLOTLIB (Hunter 2007).

References

  1. Abac, A. G., Abbott, R., Abouelfettouh, I., et al. 2024, ApJ, 970, L34 [CrossRef] [Google Scholar]
  2. Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2016a, Phys. Rev. Lett., 116, 241103 [Google Scholar]
  3. Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2016b, Phys. Rev. Lett., 116, 061102 [Google Scholar]
  4. Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2017a, Phys. Rev. Lett., 118, 221101 [CrossRef] [Google Scholar]
  5. Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2017b, Phys. Rev. Lett., 119, 161101 [Google Scholar]
  6. Abbott, R., Abbott, T. D., Abraham, S., et al. 2020a, ApJ, 896, L44 [Google Scholar]
  7. Abbott, R., Abbott, T. D., Abraham, S., et al. 2020b, Phys. Rev. Lett., 125, 101102 [Google Scholar]
  8. Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2020c, ApJ, 892, L3 [NASA ADS] [CrossRef] [Google Scholar]
  9. Abbott, R., Abbott, T. D., Abraham, S., et al. 2021, Phys. Rev. X, 11, 021053 [Google Scholar]
  10. Abbott, R., Abbott, T. D., Acernese, F., et al. 2023a, Phys. Rev. X, 13, 041039 [Google Scholar]
  11. Abbott, R., Abbott, T. D., Acernese, F., et al. 2023b, Phys. Rev. X, 13, 011048 [NASA ADS] [Google Scholar]
  12. Abbott, R., Abbott, T. D., Acernese, F., et al. 2024, Phys. Rev. D, 109, 022001 [NASA ADS] [CrossRef] [Google Scholar]
  13. Andrews, B. H., & Martini, P. 2013, ApJ, 765, 140 [NASA ADS] [CrossRef] [Google Scholar]
  14. Artale, M. C., Mapelli, M., Bouffanais, Y., et al. 2020, MNRAS, 491, 3419 [Google Scholar]
  15. Astropy Collaboration (Robitaille, T. P., et al.) 2013, A&A, 558, A33 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  16. Astropy Collaboration (Price-Whelan, A. M., et al.) 2018, AJ, 156, 123 [Google Scholar]
  17. Astropy Collaboration (Price-Whelan, A. M., et al.) 2022, ApJ, 935, 167 [NASA ADS] [CrossRef] [Google Scholar]
  18. Baibhav, V., Berti, E., Gerosa, D., et al. 2019, Phys. Rev. D, 100, 064060 [NASA ADS] [CrossRef] [Google Scholar]
  19. Barber, J., & Antonini, F. 2025, MNRAS, 538, 639 [Google Scholar]
  20. Bavera, S. S., Fragos, T., Zevin, M., et al. 2021, A&A, 647, A153 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  21. Belczynski, K., Kalogera, V., & Bulik, T. 2002, ApJ, 572, 407 [NASA ADS] [CrossRef] [Google Scholar]
  22. Belczynski, K., Dominik, M., Bulik, T., et al. 2010, ApJ, 715, L138 [Google Scholar]
  23. Belczynski, K., Romagnolo, A., Olejak, A., et al. 2022, ApJ, 925, 69 [NASA ADS] [CrossRef] [Google Scholar]
  24. Biscoveanu, S., Callister, T. A., Haster, C.-J., et al. 2022, ApJ, 932, L19 [NASA ADS] [CrossRef] [Google Scholar]
  25. Bisigello, L., Caputi, K. I., Grogin, N., & Koekemoer, A. 2018, A&A, 609, A82 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  26. Boco, L., Lapi, A., Goswami, S., et al. 2019, ApJ, 881, 157 [Google Scholar]
  27. Boco, L., Lapi, A., Chruslinska, M., et al. 2021, ApJ, 907, 110 [NASA ADS] [CrossRef] [Google Scholar]
  28. Boesky, A. P., Broekgaarden, F. S., & Berger, E. 2024a, ApJ, 976, 24 [NASA ADS] [Google Scholar]
  29. Boesky, A. P., Broekgaarden, F. S., & Berger, E. 2024b, ApJ, 976, 23 [NASA ADS] [Google Scholar]
  30. Boogaard, L. A., Brinchmann, J., Bouché, N., et al. 2018, A&A, 619, A27 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  31. Branchesi, M., Maggiore, M., Alonso, D., et al. 2023, JCAP, 2023, 068 [CrossRef] [Google Scholar]
  32. Bressan, A., Marigo, P., Girardi, L., et al. 2012, MNRAS, 427, 127 [NASA ADS] [CrossRef] [Google Scholar]
  33. Briel, M. M., Eldridge, J. J., Stanway, E. R., Stevance, H. F., & Chrimes, A. A. 2022, MNRAS, 514, 1315 [NASA ADS] [CrossRef] [Google Scholar]
  34. Broekgaarden, F. S., Berger, E., Neijssel, C. J., et al. 2021, MNRAS, 508, 5028 [NASA ADS] [CrossRef] [Google Scholar]
  35. Broekgaarden, F. S., Berger, E., Stevenson, S., et al. 2022, MNRAS, 516, 5737 [NASA ADS] [CrossRef] [Google Scholar]
  36. Caffau, E., Ludwig, H. G., Steffen, M., Freytag, B., & Bonifacio, P. 2011, Sol. Phys., 268, 255 [Google Scholar]
  37. Callister, T. A., & Farr, W. M. 2024, Phys. Rev. X, 14, 021005 [NASA ADS] [Google Scholar]
  38. Caputi, K. I., Deshmukh, S., Ashby, M. L. N., et al. 2017, ApJ, 849, 45 [Google Scholar]
  39. Casey, C. M., Zavala, J. A., Spilker, J., et al. 2018, ApJ, 862, 77 [Google Scholar]
  40. Chruślińska, M. 2024, Ann. Phys., 536, 2200170 [CrossRef] [Google Scholar]
  41. Chruslinska, M., & Nelemans, G. 2019, MNRAS, 488, 5300 [NASA ADS] [Google Scholar]
  42. Chruslinska, M., Nelemans, G., & Belczynski, K. 2019, MNRAS, 482, 5012 [Google Scholar]
  43. Chruślińska, M., Nelemans, G., Boco, L., & Lapi, A. 2021, MNRAS, 508, 4994 [CrossRef] [Google Scholar]
  44. Claeys, J. S. W., Pols, O. R., Izzard, R. G., Vink, J., & Verbunt, F. W. M. 2014, A&A, 563, A83 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  45. Costa, G., Girardi, L., Bressan, A., et al. 2019a, A&A, 631, A128 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  46. Costa, G., Girardi, L., Bressan, A., et al. 2019b, MNRAS, 485, 4641 [NASA ADS] [CrossRef] [Google Scholar]
  47. Curti, M., Mannucci, F., Cresci, G., & Maiolino, R. 2020, MNRAS, 491, 944 [Google Scholar]
  48. Di Carlo, U. N., Mapelli, M., Giacobbo, N., et al. 2020, MNRAS, 498, 495 [NASA ADS] [CrossRef] [Google Scholar]
  49. Di Stefano, R., Kruckow, M. U., Gao, Y., Neunteufel, P. G., & Kobayashi, C. 2023, ApJ, 944, 87 [NASA ADS] [CrossRef] [Google Scholar]
  50. Dominik, M., Belczynski, K., Fryer, C., et al. 2012, ApJ, 759, 52 [Google Scholar]
  51. Dominik, M., Belczynski, K., Fryer, C., et al. 2013, ApJ, 779, 72 [Google Scholar]
  52. Dorozsmai, A., & Toonen, S. 2024, MNRAS, 530, 3706 [NASA ADS] [CrossRef] [Google Scholar]
  53. Eldridge, J. J., Stanway, E. R., & Tang, P. N. 2018, MNRAS, 482, 870 [Google Scholar]
  54. Everson, R. W., MacLeod, M., & Ramirez-Ruiz, E. 2025, ApJ, 979, L11 [Google Scholar]
  55. Fishbach, M. 2025, Class. Quant. Grav., 42, 055009 [Google Scholar]
  56. Fryer, C. L., Belczynski, K., Wiktorowicz, G., et al. 2012, ApJ, 749, 91 [Google Scholar]
  57. Gallegos-Garcia, M., Berry, C. P. L., Marchant, P., & Kalogera, V. 2021, ApJ, 922, 110 [NASA ADS] [CrossRef] [Google Scholar]
  58. Ge, H., Hjellming, M. S., Webbink, R. F., Chen, X., & Han, Z. 2010, ApJ, 717, 724 [NASA ADS] [CrossRef] [Google Scholar]
  59. Ge, H., Webbink, R. F., Chen, X., & Han, Z. 2015, ApJ, 812, 40 [Google Scholar]
  60. Ge, H., Webbink, R. F., Chen, X., & Han, Z. 2020, ApJ, 899, 132 [NASA ADS] [CrossRef] [Google Scholar]
  61. Giacobbo, N., & Mapelli, M. 2018, MNRAS, 480, 2011 [Google Scholar]
  62. Giacobbo, N., & Mapelli, M. 2020, ApJ, 891, 141 [NASA ADS] [CrossRef] [Google Scholar]
  63. Gruppioni, C., Béthermin, M., Loiacono, F., et al. 2020, A&A, 643, A8 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  64. Hall, E. D., & Evans, M. 2019, Class. Quant. Grav., 36, 225002 [NASA ADS] [CrossRef] [Google Scholar]
  65. Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357 [NASA ADS] [CrossRef] [Google Scholar]
  66. Hirai, R., & Mandel, I. 2022, ApJ, 937, L42 [NASA ADS] [CrossRef] [Google Scholar]
  67. Hobbs, G., Lorimer, D. R., Lyne, A. G., & Kramer, M. 2005, MNRAS, 360, 974 [Google Scholar]
  68. Hunter, J. D. 2007, Comput. Sci. Eng., 9, 90 [NASA ADS] [CrossRef] [Google Scholar]
  69. Hurley, J. R., Tout, C. A., & Pols, O. R. 2002, MNRAS, 329, 897 [Google Scholar]
  70. Iorio, G., Mapelli, M., Costa, G., et al. 2023, MNRAS, 524, 426 [NASA ADS] [CrossRef] [Google Scholar]
  71. Ivanova, N., Justham, S., Chen, X., et al. 2013, A&ARv, 21, 59 [Google Scholar]
  72. Kalogera, V., Sathyaprakash, B. S., Bailes, M., et al. 2021, ArXiv e-prints [arXiv:2111.06990] [Google Scholar]
  73. Klencki, J., Moe, M., Gladysz, W., et al. 2018, A&A, 619, A77 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  74. Klencki, J., Nelemans, G., Istrate, A. G., & Pols, O. 2020, A&A, 638, A55 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  75. Klencki, J., Nelemans, G., Istrate, A. G., & Chruslinska, M. 2021, A&A, 645, A54 [EDP Sciences] [Google Scholar]
  76. Kroupa, P. 2001, MNRAS, 322, 231 [NASA ADS] [CrossRef] [Google Scholar]
  77. Lamberts, A., Garrison-Kimmel, S., Hopkins, P. F., et al. 2018, MNRAS, 480, 2704 [NASA ADS] [CrossRef] [Google Scholar]
  78. Lau, M. Y. M., Hirai, R., González-Bolívar, M., et al. 2022, MNRAS, 512, 5462 [NASA ADS] [CrossRef] [Google Scholar]
  79. Madau, P., & Dickinson, M. 2014, ARA&A, 52, 415 [Google Scholar]
  80. Madau, P., & Fragos, T. 2017, ApJ, 840, 39 [Google Scholar]
  81. Maiolino, R., & Mannucci, F. 2019, A&ARv, 27, 3 [Google Scholar]
  82. Mandel, I., & Broekgaarden, F. S. 2022, Liv. Rev. Relat., 25, 1 [Google Scholar]
  83. Mannucci, F., Cresci, G., Maiolino, R., Marconi, A., & Gnerucci, A. 2010, MNRAS, 408, 2115 [NASA ADS] [CrossRef] [Google Scholar]
  84. Mannucci, F., Salvaterra, R., & Campisi, M. A. 2011, MNRAS, 414, 1263 [NASA ADS] [CrossRef] [Google Scholar]
  85. Mapelli, M. 2020, Front. Astron. Space Sci., 7, 38 [NASA ADS] [CrossRef] [Google Scholar]
  86. Mapelli, M. 2021, in Handbook of Gravitational Wave Astronomy, eds. C. Bambi, S. Katsanevas, & K. D. Kokkotas, 16 [Google Scholar]
  87. Mapelli, M., & Giacobbo, N. 2018, MNRAS, 479, 4391 [NASA ADS] [CrossRef] [Google Scholar]
  88. Mapelli, M., Giacobbo, N., Ripamonti, E., & Spera, M. 2017, MNRAS, 472, 2422 [NASA ADS] [CrossRef] [Google Scholar]
  89. Mapelli, M., Giacobbo, N., Santoliquido, F., & Artale, M. C. 2019, MNRAS, 487, 2 [NASA ADS] [CrossRef] [Google Scholar]
  90. Mapelli, M., Bouffanais, Y., Santoliquido, F., Arca Sedda, M., & Artale, M. C. 2022, MNRAS, 511, 5797 [NASA ADS] [CrossRef] [Google Scholar]
  91. Marchant, P., Pappas, K. M. W., Gallegos-Garcia, M., et al. 2021, A&A, 650, A107 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  92. Misra, D., Fragos, T., Tauris, T. M., Zapartas, E., & Aguilera-Dena, D. R. 2020, A&A, 642, A174 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  93. Moe, M., & Di Stefano, R. 2017, ApJS, 230, 15 [Google Scholar]
  94. Nagarajan, P., & El-Badry, K. 2025, PASP, 137, 034203 [Google Scholar]
  95. Nakajima, K., Ouchi, M., Isobe, Y., et al. 2023, ApJS, 269, 33 [NASA ADS] [CrossRef] [Google Scholar]
  96. Neijssel, C. J., Vigna-Gómez, A., Stevenson, S., et al. 2019, MNRAS, 490, 3740 [Google Scholar]
  97. Nguyen, C. T., Costa, G., Girardi, L., et al. 2022, A&A, 665, A126 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  98. Nitz, A. H., Kumar, S., Wang, Y.-F., et al. 2023, ApJ, 946, 59 [NASA ADS] [CrossRef] [Google Scholar]
  99. Olejak, A., Belczynski, K., & Ivanova, N. 2021, A&A, 651, A100 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  100. Olejak, A., Fryer, C. L., Belczynski, K., & Baibhav, V. 2022, MNRAS, 516, 2252 [NASA ADS] [CrossRef] [Google Scholar]
  101. Pavlovskii, K., Ivanova, N., Belczynski, K., & Van, K. X. 2017, MNRAS, 465, 2092 [Google Scholar]
  102. Popesso, P., Concas, A., Cresci, G., et al. 2023, MNRAS, 519, 1526 [Google Scholar]
  103. Rauf, L., Howlett, C., Davis, T. M., & Lagos, C. D. P. 2023, MNRAS, 523, 5719 [Google Scholar]
  104. Ray, A., Hernandez, I. M., Mohite, S., Creighton, J., & Kapadia, S. 2023, ApJ, 957, 37 [NASA ADS] [CrossRef] [Google Scholar]
  105. Riley, J., Mandel, I., Marchant, P., et al. 2021, MNRAS, 505, 663 [Google Scholar]
  106. Rinaldi, P., Caputi, K. I., van Mierlo, S. E., et al. 2022, ApJ, 930, 128 [NASA ADS] [CrossRef] [Google Scholar]
  107. Rinaldi, S., Del Pozzo, W., Mapelli, M., Lorenzo-Medina, A., & Dent, T. 2024, A&A, 684, A204 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  108. Rodighiero, G., Brusa, M., Daddi, E., et al. 2015, ApJ, 800, L10 [Google Scholar]
  109. Rodriguez, C. L., & Loeb, A. 2018, ApJ, 866, L5 [NASA ADS] [CrossRef] [Google Scholar]
  110. Romagnolo, A., Belczynski, K., Klencki, J., et al. 2023, MNRAS, 525, 706 [NASA ADS] [CrossRef] [Google Scholar]
  111. Román-Garza, J., Bavera, S. S., Fragos, T., et al. 2021, ApJ, 912, L23 [CrossRef] [Google Scholar]
  112. Sana, H., de Mink, S. E., de Koter, A., et al. 2012, Science, 337, 444 [Google Scholar]
  113. Santoliquido, F., Mapelli, M., Bouffanais, Y., et al. 2020, ApJ, 898, 152 [Google Scholar]
  114. Santoliquido, F., Mapelli, M., Giacobbo, N., Bouffanais, Y., & Artale, M. C. 2021, MNRAS, 502, 4877 [NASA ADS] [CrossRef] [Google Scholar]
  115. Santoliquido, F., Mapelli, M., Artale, M. C., & Boco, L. 2022, MNRAS, 516, 3297 [NASA ADS] [CrossRef] [Google Scholar]
  116. Santoliquido, F., Mapelli, M., Iorio, G., et al. 2023, MNRAS, 524, 307 [NASA ADS] [CrossRef] [Google Scholar]
  117. Sargent, M. T., Béthermin, M., Daddi, E., & Elbaz, D. 2012, ApJ, 747, L31 [Google Scholar]
  118. Schneider, F. R. N., Podsiadlowski, P., & Müller, B. 2021, A&A, 645, A5 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  119. Schreiber, C., Pannella, M., Elbaz, D., et al. 2015, A&A, 575, A74 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  120. Sgalletta, C., Iorio, G., Mapelli, M., et al. 2023, MNRAS, 526, 2210 [NASA ADS] [CrossRef] [Google Scholar]
  121. Shao, Y., & Li, X.-D. 2021, ApJ, 920, 81 [NASA ADS] [CrossRef] [Google Scholar]
  122. Speagle, J. S., Steinhardt, C. L., Capak, P. L., & Silverman, J. D. 2014, ApJS, 214, 15 [Google Scholar]
  123. Spera, M., & Mapelli, M. 2017, MNRAS, 470, 4739 [Google Scholar]
  124. Spera, M., Mapelli, M., Giacobbo, N., et al. 2019, MNRAS, 485, 889 [Google Scholar]
  125. Srinivasan, R., Lamberts, A., Bizouard, M. A., Bruel, T., & Mastrogiovanni, S. 2023, MNRAS, 524, 60 [CrossRef] [Google Scholar]
  126. Stevenson, S., & Clarke, T. A. 2022, MNRAS, 517, 4034 [Google Scholar]
  127. Tang, P. N., Eldridge, J. J., Stanway, E. R., & Bray, J. C. 2020, MNRAS, 493, L6 [Google Scholar]
  128. Team COMPAS (Riley, J., et al.) 2022, ApJS, 258, 34 [NASA ADS] [CrossRef] [Google Scholar]
  129. The Pandas Development Team 2024, https://doi.org/10.5281/zenodo.3509134 [Google Scholar]
  130. Trani, A. A., Rieder, S., Tanikawa, A., et al. 2022, Phys. Rev. D, 106, 043014 [NASA ADS] [CrossRef] [Google Scholar]
  131. van Son, L. A. C., de Mink, S. E., Callister, T., et al. 2022, ApJ, 931, 17 [NASA ADS] [CrossRef] [Google Scholar]
  132. van Son, L. A. C., de Mink, S. E., Chruślińska, M., et al. 2023, ApJ, 948, 105 [NASA ADS] [CrossRef] [Google Scholar]
  133. van Son, L. A. C., Roy, S. K., Mandel, I., et al. 2025, ApJ, 979, 209 [Google Scholar]
  134. Verbunt, F., Igoshev, A., & Cator, E. 2017, A&A, 608, A57 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  135. Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nat. Meth., 17, 261 [Google Scholar]
  136. Webbink, R. F. 1984, ApJ, 277, 355 [NASA ADS] [CrossRef] [Google Scholar]
  137. Zevin, M., Spera, M., Berry, C. P. L., & Kalogera, V. 2020, ApJ, 899, L1 [Google Scholar]
  138. Zevin, M., Bavera, S. S., Berry, C. P. L., et al. 2021, ApJ, 910, 152 [Google Scholar]

Appendix A: Observational scaling relations

Below, we provide the details of the observational scaling relations adopted in this work. We follow the approach presented by Santoliquido et al. (2022).

A.1. Galaxy stellar mass function

Our adopted galaxy stellar mass function takes the form of a Schechter function:

ϕ ( M , z ) d M = ϕ N ( z ) e M / M cut ( z ) ( M M cut ( z ) ) α GSMF d M , $$ \begin{aligned} \phi (M_*, z)\,dM_*= \phi _{\rm N}(z)\,e^{-M_*/M_{\rm cut}(z)} \left(\frac{M_*}{M_{\rm cut}(z)}\right)^{-\alpha _{\rm GSMF}} dM_*, \end{aligned} $$(A.1)

where ϕN(z) is a normalization factor and Mcut(z) is the galaxy stellar mass at which the function transitions from a power law at low masses to an exponential decay at high masses. In our models, we use the fit derived by Chruslinska & Nelemans (2019) and based on a wide catalog of galaxy stellar mass function relations over separate redshift bins.

A.2. Star formation rate (SFR)

We draw the SFR of galaxies from the following double lognormal distribution (Sargent et al. 2012; Rodighiero et al. 2015; Schreiber et al. 2015):

P ( log SFR | M , z ) = A MS exp [ ( log SFR log SFR MS ) 2 2 σ MS 2 ] + A rmSB exp [ ( log SFR log SFR SB ) 2 2 σ SB 2 ] , $$ \begin{aligned} \mathcal{P} (\log \mathrm{SFR} | M_*, z) = A_{\rm MS} \exp \left[- \frac{\left(\log \mathrm{SFR} - \langle \log \mathrm{SFR} \rangle _{\rm MS} \right)^{2}}{ 2 \sigma _{\rm MS}^{2}}\right] + A_{rm SB} \exp \left[- \frac{\left(\log \mathrm{SFR} - \langle \log \mathrm{SFR} \rangle _{\rm SB} \right)^{2}}{2 \sigma _{\rm SB}^{2}}\right], \end{aligned} $$(A.2)

where AMS = 0.97 and ASB = 0.243 constants (Sargent et al. 2012). ⟨logSFR⟩MS and σMS are the average SFR of the main sequence and its dispersion. We adopt several definitions of the galaxy main sequence, as described in section 2.4. For all of our models we use the fiducial value σMS = 0.188 dex (Santoliquido et al. 2022). In a similar way, ⟨logSFR⟩SB and σSB are the average and the standard deviation of the galaxy starburst sequence. The latter is defined as in Sargent et al. (2012),

log SFR SB = log SFR MS + 0.59 $$ \begin{aligned} \langle \log \mathrm{SFR} \rangle _{\rm SB} = \langle \log \mathrm{SFR} \rangle _{\rm MS} + 0.59 \end{aligned} $$(A.3)

Moreover, we assume σSB = 0.243 dex. We couple the galaxy stellar mass function with the distribution of SFR as described above.

A.3. Fundamental metallicity relations

We consider three different fundamental metallicity relations. We adopt the fit from Mannucci et al. (2011):

12 + log ( O / H ) = 8.90 + 0.37 m 0.14 s 0.19 m 2 + 0.12 m s 0.054 s 2 for μ 0.32 9.5 = 8.93 + 0.51 ( μ 0.32 10 ) for μ 0.32 < 9.5 $$ \begin{aligned} \begin{aligned} 12 + \log \left(\mathrm{O/H} \right)&= 8.90+0.37 m - 0.14 s - 0.19 m^2 + 0.12 m s - 0.054 s^2 \quad \; \mathrm{for} \; \mu _{0.32} \ge 9.5 \\&= 8.93 + 0.51(\mu _{0.32} - 10) \quad \; \mathrm{for} \; \mu _{0.32} < 9.5 \end{aligned} \end{aligned} $$(A.4)

where m = log M* − 10 and s = logSFR and μ0.32 = log M* − 0.32logSFR, all the quantities are in solar units.

We derive an additional fit for the fundamental metallicity relation based on figure 12 by Andrews & Martini (2013)

12 + log ( O / H ) = 0.43 ( log M 0.66 log SFR ) + 4.59 $$ \begin{aligned} 12 + \log (\mathrm{O/H}) = 0.43 (\log M_*- 0.66 \log \mathrm{SFR}) + 4.59 \end{aligned} $$(A.5)

Lastly, we consider the metallicity relation calculated in Curti et al. (2020):

12 + log ( O / H ) = Z 0 γ / β log ( 1 + ( M / M 0 ( SFR ) ) β ) $$ \begin{aligned} 12 + \log (\mathrm{O/H}) = Z_{0} - \gamma /\beta \log \left(1 + (M_*/M_0(\mathrm{SFR}))^{-\beta }\right) \end{aligned} $$(A.6)

where M0(SFR) = m0 + m1logSFR and Z0 = 8.779 ± 0.005, m0 = 10.11 ± 0.03, m1 = 0.56 ± 0.01, γ = 0.31 ± 0.01 and β = 2.1 ± 0.4. For the purposes of our simulations, we convert these relations into absolute metallicities. We adopt the solar metallicity values from Caffau et al. (2011): Z = 0.0153 and 12 + log(O/H) = 8.76.

Figure A.1 compares the distributions of SFR density as a function of galaxy stellar mass and metallicity, for different fundamental metallicity relations at z = 0. The differences are particularly noticeable for low–mass galaxies, where the model by Mannucci et al. (2011) clearly predicts lower metallicities compared to the other two prescriptions. Curti et al. (2020) predicts the flattest relation, with less than an order of magnitude difference in metallicty between the low and high–mass galaxies.

thumbnail Fig. A.1.

SFR density per galaxy mass (x axes) and metallicity bins (y axes) at z = 0. Each panel shows the result for a different FMR relation, from left to right: Mannucci et al. (2011), Andrews & Martini (2013) and Curti et al. (2020).

Appendix B: Merger rate densities

Figures B.1a and B.1b show the BBH cosmic merger rate densities obtained for different SFR–M* relations, adopting the fundamental metallicity relations by Curti et al. (2020) and Mannucci et al. (2011), respectively.

thumbnail Fig. B.1.

BBH cosmic merger rate densities for different SFR–M* relations, using two different fundamental metallicity relations. (a) Same as Fig. 5 but for the FMR by Curti et al. (2020). (b) Same as Fig. 5 but for the FMR by Mannucci et al. (2011).

All Tables

Table 1.

Summary of the grid of parameters varied in our study.

All Figures

thumbnail Fig. 1.

SFR density as a function of redshift for different SFR–M* relations (different colors). The different linestyles show the results by cutting the galaxy main sequence at different minimum galaxy masses, Mmin = 106 M (dotted), 107 M (solid) and 108 M (dashed). The wide blue line shows the fit by Madau & Fragos (2017). The markers show observed data points: the blue dots are are taken from Madau & Dickinson (2014), the orange diamonds are from Casey et al. (2018), and the purple squares are taken from Gruppioni et al. (2020).

In the text
thumbnail Fig. 2.

Schematic overview of the different metal-dependent SFR models adopted here. Upper panel: overview of the detailed Universe model, based on observational scaling relations. Lower panel: overview of the model based on the average SFR and metallicity evolution across cosmic time. The plots on the right illustrate the metallicity distributions as a function of redshift for our different prescriptions. For the average metallicity evolution model, we show the distributions assuming both σ Z = 0.1 $ \tilde{\sigma}_{\mathrm{Z}} = 0.1 $ and 0.6.

In the text
thumbnail Fig. 3.

Merger efficiency η for BBHs as a function of metallicity Z. From top to bottom α = 0.5, 1, 3, 5. The different colors refer to formation channels I–IV (Section 2.3). The vertical red lines highlight metallicity values of Z = 1/100, 1/10, and 1/3 Z.

In the text
thumbnail Fig. 4.

Contribution of galaxies of mass M* to the SFR density evaluated for different metallicity thresholds at z = 0. The x axis represents the galaxy stellar mass M*, while the lines show the total contribution to the SFR density from all galaxies in a given stellar mass bin, with a metallicity Z ≤ Zth. The different line styles refer to different metallicity thresholds Zth; solid lines: Zth = 1/3 Z, dashed: Zth = 1/10 Z, and dotted: Zth = 1/100 Z. The different colors represent models with different SFR–M* relations, blue: B18, violet: S14, pink: P23, orange: SHALLOW, yellow: STEEP. Each panel shows the results for a different FMR, from left to right: Mannucci et al. (2011), Andrews & Martini (2013), Curti et al. (2020).

In the text
thumbnail Fig. 5.

BBH merger rate density as a function of redshift for different observational scaling relations (different colors) and CE α parameters (different panels). Each color is associated with a different SFR–M* relation, as summarized in the legend in the top left-hand panel. The solid lines show the model with Mmin = 107 M. The shaded area encompasses the difference between Mmin = 106 M (upper boundary), and Mmin = 108 M (lower boundary). All the models assume the fiducial fundamental metallicity relation by Andrews & Martini (2013). The thick dashed lines show the merger rate density obtained with COSMOATE (simplified model), assuming the cosmic SFR density by Madau & Dickinson (2014) and a lognormal distribution for the metallicity with a spread σ Z = 0.1 $ \tilde{\sigma}_{\mathrm{Z}}=0.1 $ (dark gray) and 0.6 (light gray). The gray-shaded regions show the MRD inferred by the LVK collaboration: ℛ(z)∝(1 + z)k (Abbott et al. 2023b).

In the text
thumbnail Fig. 6.

BBH merger rate density as a function of redshift for different BBH formation channels. The solid black lines show the merger rate density of the whole BBH population. The shaded area encompasses the difference between Mmin = 106 M (upper boundary), and Mmin = 108 M (lower boundary). The different panels show the results obtained assuming different CE α parameters. All the models assume the fiducial SFR–M* relation P23 and the fundamental metallicity relation by Andrews & Martini (2013). The gray-shaded regions show the MRD inferred by the LVK collaboration (Abbott et al. 2023b).

In the text
thumbnail Fig. 7.

Merger rate density as a function of redshift, for the natal-kick models GM20 (solid lines) and H05 (dashed lines). The line colors show the results obtained assuming different CE efficiency parameters (α). The upper boundaries of the shaded areas show the models obtained with minimum galaxy mass M* = 106 M, whereas for the lower boundary the minimum mass is 108 M. All the models assume the fiducial observational scaling relations. The gray-shaded region shows the MRD inferred by the LVK collaboration (Abbott et al. 2023b).

In the text
thumbnail Fig. A.1.

SFR density per galaxy mass (x axes) and metallicity bins (y axes) at z = 0. Each panel shows the result for a different FMR relation, from left to right: Mannucci et al. (2011), Andrews & Martini (2013) and Curti et al. (2020).

In the text
thumbnail Fig. B.1.

BBH cosmic merger rate densities for different SFR–M* relations, using two different fundamental metallicity relations. (a) Same as Fig. 5 but for the FMR by Curti et al. (2020). (b) Same as Fig. 5 but for the FMR by Mannucci et al. (2011).

In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.