How do analytical star formation history models affect binary black hole merger predictions?

Testing fits against simulated star formation

The paper asks whether the simple analytical formulas commonly used to describe cosmic star formation are good enough for predicting binary black hole (BBH) mergers. The star formation history in question is the metallicity-specific star formation rate density S(Z, z): how much star formation occurs per unit volume, per unit time, at a given metallicity Z and redshift z. Metallicity matters because it shapes the evolution of massive stars and therefore the black holes they leave behind.

The authors take S(Z, z) directly from three IllustrisTNG cosmological simulations of different mass resolution and box size — TNG50-1, TNG100-1 and TNG300-1 — and fit each one with an analytical form combining a Madau–Dickinson-style S(z) with a skewed log-normal metallicity distribution. They then feed both the full simulated S(Z, z) and the fitted version into COMPAS, a binary population synthesis code, using the same COMPAS simulations throughout so that differences between the two come from the star formation input rather than from binary physics.

This design isolates one modelling choice. The comparison is not between simulations and observations, but between a detailed simulated star formation history and the simplified fit meant to stand in for it.

The fit misses the metallicity distribution's shape

The analytical fit reproduces the simulation best around 10⁻¹ ≲ Z/Z☉ ≲ 10⁰, but it fails at both ends. The simulated S(Z) has one peak near Z/Z☉ ∼ 1 and a second bump near Z/Z☉ ≈ 2–5; the fit cannot capture that second bump and instead places a single peak at Z/Z☉ ≈ 1–2, with the exact position shifting with resolution. Above Z/Z☉ ≳ 6 the fit overestimates S(Z) by up to a factor of roughly 10–1000.

The low-metallicity tail is also misrepresented, in a way that depends on the simulation. TNG50-1 shows a change in slope between 10⁻³ ≲ Z/Z☉ ≲ 10⁻², while TNG100-1 and TNG300-1 flatten below Z/Z☉ < 10⁻². The authors note that a better fit would need a steeper high-metallicity cutoff and a low-metallicity tail with at least one additional change in slope.

In redshift, the fit systematically overestimates S(z) at z ≲ 0.5 by up to a factor of about 1.5, and it cannot reproduce the turning point near z ≈ 6–8 seen in the simulations. At z ≈ 14 it overestimates S(z) for all three simulations, by about a factor of 10 for TNG300-1 and about 2 for TNG50-1 and TNG100-1. The authors attribute this partly to an optimization that prioritizes the high star formation rate near the peak and gives less weight to low-rate regions such as high redshift.

High-redshift merger rates inflate by orders of magnitude

Feeding the fits into COMPAS changes the predicted BBH merger rate. At z ≈ 0 the local rate ranges from about 29 to 59 Gpc⁻³ yr⁻¹ across the three simulations when the full simulated S(Z, z) is used, and from about 28 to 74 Gpc⁻³ yr⁻¹ when the analytical fit is used. For TNG50-1 the fit overestimates the local rate by about 15 Gpc⁻³ yr⁻¹; for TNG100-1 and TNG300-1 the simulated and fitted local rates agree within ±2 Gpc⁻³ yr⁻¹.

The larger discrepancy appears at high redshift. The paper reports that analytical S(Z, z) prescriptions can overestimate BBH merger rates at z ≳ 6 by up to a factor of 10–10⁴, depending on simulation resolution, with the effect growing at higher redshift and at lower resolution or larger box size. The fit also shifts when the merger rate peaks: with decreasing resolution and increasing box size, the peak moves to higher redshift.

The local-rate comparison and the high-redshift comparison answer different questions and should not be conflated. A fit that matches the z ≈ 0 normalization for two of the three simulations can still diverge strongly at z ≳ 6.

An artificial 8 M☉ feature appears in the mass distribution

The fit also changes the shape of the primary mass distribution, not just its normalization. Using the analytical fit produces an artificial excess near 8 M☉ in the primary mass distribution at z_merger = 0.2, by up to a factor of 7 at z_merger ≲ 2, a feature absent when the full simulation-based S(Z, z) is used. The paper suggests this feature may originate primarily from the common-envelope channel. At the same time, the fit suppresses high-mass features: the high-mass rate above M₁ ≳ 20 M☉ is increasingly suppressed with increasing primary mass at z_merger = 0.2, including flattening of the peaks near 33 M☉ and 45 M☉.

The authors trace this back to the metallicity distribution. Because the fit smooths over the secondary bump near 3 Z☉, it places the global maximum of the metallicity distribution at a higher metallicity than the simulation does, which propagates into the predicted merger population.

One thing stays stable: the locations of features in the primary mass distribution remain consistent across redshift and between the fit and simulation-based S(Z, z) for all three simulations. The paper concludes that binary modelling assumptions more strongly set where features sit, while S(Z, z) affects the overall rate and the relative prominence of features.

Resolution changes normalization more than shape

The choice of IllustrisTNG run mainly shifts the normalization of the results. TNG50-1, with the highest mass resolution and smallest box, yields the largest star formation rate at all redshifts, while TNG300-1 is suppressed by roughly a factor of about 2. The local merger rate follows the same ordering, from about 29 Gpc⁻³ yr⁻¹ for TNG300-1 to about 59 Gpc⁻³ yr⁻¹ for TNG50-1 with the simulated S(Z, z). Near the merger-rate peak, TNG50-1 exceeds TNG300-1 by roughly a factor of three, and at z ≳ 6 by up to two orders of magnitude.

The paper is careful about what resolution means here. Higher resolution does not automatically mean a more reliable star formation history, because subgrid prescriptions remain in place and smaller boxes increase cosmic variance and limit statistics. The authors also note that resolution and box size vary together across these three runs, so separating their effects requires a larger simulation set, which they leave to future work.

The comparison to the observed GWTC-4 primary mass distribution is also limited. The models underpredict the rate at M_BH,1 ≳ 35 M☉, which the authors suggest may require additional formation channels such as hierarchical mergers in dynamical environments, beyond the isolated binary evolution modelled here.

Where the fit is usable and where it is not

The paper's practical split is between shape and normalization. The analytical fit for TNG50-1 most accurately represents the shape of its S(Z, z), while at lower resolution and larger volume the rates are lower and the disagreement between fit and simulation is greatest. The fit reproduces the merger-rate evolution most accurately between z ≈ 1–4 for TNG50-1, at z ≲ 2 for TNG100-1, and at z ≲ 1.5 for TNG300-1.

That makes the choice of star formation model a redshift-dependent decision rather than a single yes-or-no judgement. A fit calibrated to match the local rate can still distort the high-redshift rate and the prominence of mass features, and the size of that distortion depends on which simulation supplies the underlying S(Z, z).

The authors' stated priority is to improve S(Z, z) modelling rather than to abandon analytical fits outright. They call for checking whether the artificial features seen with IllustrisTNG also appear in other cosmological simulations and in observationally derived star formation histories; if they do, the case for moving toward detailed star formation histories, or for updating the fits, becomes stronger. Until then, the paper's clearest warning is that a fit matching the z ≈ 0 merger rate is not evidence that its high-redshift predictions or its mass-distribution shape can be trusted.

Sources used in this answer