From cosmological simulations to binary black hole mergers: The impact of using analytical star formation history models on gravitational-wave source populations

TL;DR

This work quantifies how well analytical, two-dimensional fits to the metallicity-dependent star formation history $\mathcal{S}(Z,z)$ reproduce the full IllustrisTNG simulation histories and investigates the consequences for binary black hole (BBH) merger populations modeled with COMPAS. By comparing simulation-based $\mathcal{S}(Z,z)$ to a nine-parameter analytic fit across three TNG resolutions, the authors show that the analytic form can overestimate high-redshift BBH rates by up to $\sim10^4$ and introduce artificial features in the BBH mass distribution, such as a spurious $\sim8\,M_\odot$ bump. The discrepancies arise from missing high-metallicity bumps and a flattened low-metallicity tail in the analytic $\mathcal{S}(Z,z)$, highlighting the sensitivity of BBH predictions to the detailed shape of the cosmic star formation history. The results underscore the need to couple cosmological simulations with population synthesis carefully and to treat high-$z$ and high-$Z$ regimes with flexible, simulation-informed models to interpret gravitational-wave observations accurately.

Abstract

Observations of binary black hole (BBH) mergers provide a unique window into the lives of massive stars across cosmic time. Connecting redshift-dependent merger properties to massive star progenitors requires accurate models of cosmic star formation and chemical enrichment histories. Analytical fits for the metallicity-specific cosmic star formation rate density S(Z, z) are commonly used as proxies for the complex underlying star formation history, yet they remain unconstrained. Using the IllustrisTNG cosmological simulations, we evaluate the accuracy of these analytical S(Z, z) prescriptions and assess how simulation resolution and volume affect the inferred S(Z, z). By coupling the simulated and analytical S(Z, z) to the population synthesis code COMPAS, we investigate the resulting BBH merger rates and mass distributions. We find that analytical S(Z, z) prescriptions can overestimate BBH merger rates at high redshift ($z \gtrsim 6$) by up to a factor of $10$-$10^4$, depending on cosmological simulation resolution, and can introduce spurious features in the BBH mass distribution. For example, they can produce an artificial feature near $8\,M_\odot$ in the primary mass distribution at $z \lesssim 2$, which is absent when using the full simulation-based S(Z, z), while simultaneously suppressing high-mass features. These discrepancies arise because simple analytical models fail to capture a high-metallicity bump and a more flattened low-metallicity tail in the simulated S(Z, z) metallicity distribution. Our results highlight the importance of accurate star formation histories for modeling BBH populations, demonstrate the limitation of widely used analytical S(Z, z) fits, and underscore the need for careful integration of cosmological simulations, analytical fits, and population synthesis when interpreting gravitational-wave observations.

Figures (12)

• Figure 1: Bottom panel: Star formation rate density, $\mathcal{S}(z)$, as a function of redshift, $z$, for the IllustrisTNG simulations TNG50-1 (blue), TNG100-1 (magenta), and TNG300-1 (green). Solid lines show $\mathcal{S}(z)$ obtained by summing the full TNG $\mathcal{S}(Z,z)$ from each simulation over metallicity, while dashed lines are analytical fits to the TNG $\mathcal{S}(Z,z)$. For comparison, $\mathcal{S}(z)$ model from MadauFragos2017, which has the parameters $a=0.01$, $b=2.6$, $c=3.2$, and $d=6.2$ as defined in Equation \ref{['eq:sfr']}, is shown as a gray line. Top panel: The fractional error for $\mathcal{S}(z)$ as a function of redshift, as described in Section \ref{['sec:comparison methods']}.
• Figure 2: Bottom panel: Star formation rate density, $\mathcal{S}(Z)$, as a function of $Z/Z_\odot$. The colors and line styles are analogous to Figure \ref{['fig:sfr_z']}. Top panel: The percent error for $\mathcal{S}(Z)$.
• Figure 3: TNG100-1 $\mathcal{S}(Z,z)$ as a function of lookback time (Gyr) and $Z/Z_\odot$. The color bar shows the simulation $\mathcal{S}(Z,z)$ and the black dashed contours show $\mathcal{S}(Z,z)$ from the analytical fit. The top panel shows $\mathcal{S}(Z,z)$ marginalized over metallicity, and the right panel over lookback time. Note that the side panels are plotted using a linear scale for $\mathcal{S}(z)$ and $\mathcal{S}(Z)$. Black regions near the upper right reflect numerical noise from the interpolation edge.
• Figure 4: Bottom panel:BBH merger-rate density as a function of redshift. The observed BBH merger rate based on GWTC-4 data from LVK2025 is shown in gray with the 95% credible regions. Top panel: the ratio between the rate from the analytical fit and TNG model (see Section \ref{['sec:comparison methods']}). The noisy features in the ratio around $z\gtrsim13$ arise from division by very small numbers. Colors and line styles match those in Figure \ref{['fig:sfr_z']}.
• Figure 5: Redshift evolution of the BBH primary mass distribution, $M_{\rm{BH,1}}$ [${\rm Gpc}^{-3}{\rm yr}^{-1}\,\rm{M}_{\odot}\xspace^{-1}$], for the TNG100-1, computed using the full simulation $\mathcal{S}(Z,z)$ (left column) and the analytical fit (right column). The darkest color represents the local mass distribution $z_\mathrm{merger} = 0.2$ and the lightest is at $z_\mathrm{merger}=8$. The B-spline BBH primary mass distribution from GWTC-4 LVK2025 is shown in gray. The rows represent different formation channels: the top row includes mergers from both the stable mass transfer and common-envelope (CE) channels, the middle row shows only stable mass transfer channel mergers, and the bottom row includes only CE channel mergers.