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

Technical Summary: From Cosmological Simulations to Binary Black Hole Mergers

Problem Statement
Observations of binary black hole (BBH) mergers provide critical insights into massive star evolution across cosmic time. However, connecting observed merger properties to their progenitor stars requires accurate models of the cosmic star formation history (SFRH) and chemical enrichment. Currently, the field frequently relies on analytical fits for the metallicity-dependent cosmic star formation rate density, $S(Z, z)$, as proxies for complex underlying histories derived from cosmological simulations or observations. These analytical prescriptions are often unconstrained and their impact on inferred BBH populations remains insufficiently understood. Furthermore, the effects of cosmological simulation resolution and box size on the derived $S(Z, z)$ and subsequent BBH merger rates have not been systematically quantified.

Methodology
The authors utilize the IllustrisTNG suite of cosmological simulations (specifically TNG50-1, TNG100-1, and TNG300-1), which vary in resolution and box size, to derive simulation-based $S(Z, z)$ directly from gas cell data. They compare these direct simulation results against a widely used analytical fit model (based on L. A. C. van Son et al. 2023), which combines a Madau & Dickinson (2014) functional form for total SFR density, $S(z)$, with a skewed log-normal metallicity distribution.

To assess the impact of these different $S(Z, z)$ inputs on BBH populations, the authors couple both the simulation-derived and analytical $S(Z, z)$ models to the COMPAS (Compact Object Mergers: Population Astrophysics and Statistics) rapid binary population synthesis code. They employ a "cosmic integration" framework to calculate merger rates and mass distributions across redshift. The study specifically investigates:

1. The fidelity of the analytical fit in reproducing the shape and normalization of the $S(Z, z)$ from different TNG resolutions.
2. The resulting deviations in BBH merger rates and primary mass distributions ($M_{BH,1}$) when using the analytical fit versus the full simulation data.
3. The influence of simulation resolution (ranging from high-resolution TNG50-1 to lower-resolution TNG300-1) on these outcomes.

Key Results
The study finds that analytical prescriptions for $S(Z, z)$ introduce significant systematic errors when modeling BBH populations:

• Discrepancies in $S(Z, z)$: The analytical fit fails to capture key structural features of the simulation-based metallicity distribution. Specifically, it misses a secondary "bump" in the high-metallicity regime ($Z \gtrsim 3 Z_\odot$) and a flattened low-metallicity tail present in the simulations. The fit tends to overestimate the star formation rate at high redshifts ($z \gtrsim 6$) and high metallicities, while underestimating contributions at low metallicities and very high redshifts, depending on the simulation resolution.
• Impact on Merger Rates: At high redshifts ($z \gtrsim 6$), the analytical fit overestimates BBH merger rates by factors ranging from 10 to $10^4$ compared to the simulation-based rates, with the magnitude of the error increasing as simulation resolution decreases. At low redshift ($z \approx 0$), the local merger rates derived from the fit generally agree with simulation-based rates within a factor of $\sim 1.5$, though TNG50-1 shows a $\sim 15$ Gpc$^{-3}$ yr$^{-1}$ overestimation.
• Impact on Mass Distributions: The use of analytical fits introduces spurious features in the BBH primary mass distribution. Most notably, at $z \lesssim 2$, the fit produces an artificial excess (bump) near $8 M_\odot$ in the primary mass distribution, which is absent when using the full simulation-based $S(Z, z)$. This feature is attributed to the stable mass transfer channel and is caused by the fit's inability to correctly model the high-metallicity star formation peak. Conversely, the fit suppresses high-mass features (near $33 M_\odot$ and $45 M_\odot$) by factors of 3–5 at $M_{BH,1} \gtrsim 35 M_\odot$.
• Resolution Dependence: Simulation resolution significantly affects the normalization and peak location of $S(z)$. Higher-resolution simulations (TNG50-1) yield higher star formation rates and peak at slightly lower redshifts ($z \approx 2$) compared to lower-resolution runs ($z \approx 2.5$). The analytical fit's accuracy degrades with lower simulation resolution, particularly in reproducing the high-redshift slope.

Significance and Claims
The authors conclude that while analytical fits offer computational convenience, they are unreliable proxies for the complex, underlying $S(Z, z)$ derived from cosmological simulations, particularly when interpreting gravitational-wave observations. The study demonstrates that simple analytical models fail to capture critical physical features—such as high-metallicity bumps and specific low-metallicity tails—that directly influence the predicted BBH merger rates and mass distributions.

The paper emphasizes that the discrepancies arise because the functional forms of current analytical fits cannot accommodate the multi-modal nature of metallicity distributions found in simulations. Consequently, relying on these fits may lead to incorrect inferences about the formation channels and evolutionary history of BBHs. The authors argue for the careful integration of cosmological simulations, analytical fits, and population synthesis to avoid introducing spurious features into the interpretation of gravitational-wave data. They do not claim that analytical fits are entirely useless but highlight their limitations and the necessity of validating them against high-fidelity simulation data, especially as next-generation detectors probe higher redshifts.