BBH merger populations using cosmic star formation histories of cosmological simulations

This code is used for the paper:

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

Submitted to ApJ and also available on arXiv 2601.20202

Model the binary black hole (BBH) merger population by obtaining the metallicity-dependent cosmic star formation rate density SFRD(Z, z) from an IllustrisTNG cosmological simulation. The "cosmic integration" step is done using the COMPAS cosmic integration post-processing tool. The SFRD(Z, z) can be modeled using a 2D array (obtained from a TNG simulation) or an analytical fit to this 2D array following the methods described in Neijssel et al. 2019 and van Son et al. 2023.

Required inputs:

• COMPAS output file, download here
• GWTC-4 BBH B-spline model file (BBHMassSpinRedshift_BSplineIID.h5), download here

To obtain the SFRD(Z, z) from a TNG simulation:

• Run get_TNGdata_withmetals.py
• Returns SFRMetallicityFromGasWithMetalsTNG[version].hdf5

To fit the analytical model to a TNG simulation SFRD(Z, z):

• Run Fit_model_TNG_SFRD.py
• Returns test_best_fit_parameters_TNG[version].txt

To model BBH merger population using TNG simulation or analytical fit SFRD(Z, z)

• See cosmic_int_guide.txt inside the CosmicIntegration directory for additional detail
• Run CosmicIntegration/FastCosmicIntegration-withdata.py or CosmicIntegration/FastCosmicIntegration.py
• If using slurm:
  • Set up init_values.py
  • TNG simulation SFRD(Z, z) (as a 2D array): CallCosmicIntegration_data.py
  • Fitted SFRD(Z, z): CallCosmicIntegration_full.py OR CallCosmicIntegration_data.py without specifying a file containing the SFRD (remove " --datafname " + data_file_name)
• Returns Rate_info_TNG[version].h5

To reproduce the figures:

• Fig 1-2, 11: compare_SFRDparams.py
  • compare_SFR (Fig 1), compare_Zdist (Fig 2)
  • Zdist_3panel_plots (Fig 11)
• Fig 3, 10: SFRD_2Dplots.py
  • SFRD_2Dplot_sidepanels with plotregions=False (Fig 3) and plotregions=True (Fig 10)
• Fig 4-9, 12: TNG_BBHpop_properties.py
  • compare_BBH_data_and_model_rates (Fig 4)
  • plot_BBH_mass_dist_formation_channels (Fig 5)
  • compare_BBH_data_and_model_mass_dist_over_z with plotdiff=False (Fig 6) and plotdiff=True (Fig 12)
  • residuals_BBH_data_and_model_mass_dist (Fig 7)
  • plot_BBH_mass_Z_z with fractionalerror=False (Fig 8) and fractionalerror=True (Fig 9)