Simulating the source population¶
GWForge can simulate a wide range of binary source populations (at the moment). This page tells you how you can use them.
The first step in generating a source population is to determine the distribution of sources in luminosity distance \(D_L\) (or equivalently redshift \(z\)) and the expected number of signals in the data. So, we will start by setting up the [Redshift] section.
Redshift¶
For this, you need to specify:
Redshift distribution model
Local merger rate density in \(Gpc^{-3}yr^{-1}\)
Maximum redshift of the source
Cosmological parameters such as \(H_0,~O_{m0}, O_{de0}\) and \(T^{cmb}_0\), assuming LambdaCDM cosmology
A reference start time when you switch on your detector.
The last is optional. If not provided, GWForge assumes Planck18 cosmology. $$ H_0 = 67.66~km/s/Mpc,~O_{m0} = 0.30966,~T^{cmb}_0 = 2.7255 K $$
Alternatively, you can select any of the cosmological realisations implemented in astropy.
If you choose to define a custom Universe, you can do so as follows:
[Redshift]
redshift-model = MadauDickinson
redshift-parameters = {'gamma': 2.7, 'kappa': 5.6, 'z_peak': 1.9}
local-merger-rate-density = 22
maximum-redshift = 30
; custom cosmology
cosmology = custom
H0 = 70
Om0 = 0.3
Ode0 = 0.7
Tcmb0 = 2.735
; analysis start time
gps-start-time = 1893024018
For reference, the Madau-Dickinson distribution function is: $$ p(z | \gamma, \kappa, z_\mathrm{peak}) \propto \frac{1}{1+z} \frac{dV_c}{dz} \psi(z | \gamma, \kappa, z_\mathrm{peak}) $$ where $$ \psi(z | \gamma, \kappa, z_\mathrm{peak}) \propto \frac{(1+z)^\gamma}{1+\big(\frac{1+z}{1+z_\mathrm{peak}}\big)^\kappa} M_\odot/\mathrm{year}/\mathrm{Mpc}^3 $$ and
parameter |
descsription |
|---|---|
\(\gamma\) |
Slope of the distribution at low redshift |
\(\kappa\) |
Slope of the distribution at high redshift |
\(z_\textrm{peak}\) |
Redshift at which the distribution peaks. |
\(z_\textrm{max}\) |
The maximum redshift allowed. |
It is important to note that \(\psi(z)\) describes the progenitor formation rate distribution. To obtain the compact binary merger rate distribution, GWForge convolves it with a time-delay distribution \(p(\tau)\) between the formation of the binary and its eventual merger: $$ \mathcal{R}(z_m) \propto \int_{z_m}^{\infty} \psi(z_f), p\big(t(z_f) - t(z_m)\big), \frac{dt}{dz_f}, dz_f, $$ where \(t(z)\) is the lookback time and the lower limit \(z_f \geq z_m\) enforces that a binary must form before it merges.
By default GWForge uses the inverse time-delay distribution, \(p(\tau) \propto 1/\tau\) (i.e. flat-in-log), between a minimum delay of \(20~\mathrm{Myr}\) and the age of the Universe. This is the canonical choice for compact binaries formed through isolated binary evolution via the common-envelope phase. If you want a different slope, use the powerlaw model, \(p(\tau) \propto \tau^{-\mathrm{slope}}\), and set time-delay-slope:
[Redshift]
redshift-model = MadauDickinson
redshift-parameters = {'gamma': 2.7, 'kappa': 5.6, 'z_peak': 1.9}
local-merger-rate-density = 22
maximum-redshift = 10
gps-start-time = 1893024018
; power-law time delay instead of the default 1/tau
time-delay-model = powerlaw
time-delay-slope = 0.8
parameter |
description |
|---|---|
|
|
|
Slope of the power-law time delay. Ignored (fixed to 1) for |
If, instead of Madau-Dickinson, you would rather scale the merger rate as a simple power law of redshift, \(\psi(z) = (1+z)^{\lambda}\), use the PowerLaw model:
[Redshift]
redshift-model = PowerLaw
redshift-parameters = {'lamb': 2.7}
local-merger-rate-density = 22
maximum-redshift = 10
gps-start-time = 1893024018
The reader can refer to the following for further details:
Note
The redshift, time-delay, mass and spin models are all implemented in-house on top of bilby — GWForge no longer depends on gwpopulation or pycbc to build the population. The numbers are validated against those original implementations, so nothing changes for you except a lighter set of dependencies.
Mass¶
The [Mass] section helps define the mass distribution of the binary population. Similar to the [Redshift] section, a model name and a dictionary of parameters must be provided. For example:
[Mass]
mass-model = PowerLaw+Peak
mass-parameters = {'alpha':3.37, 'beta': 0.76, 'delta_m':5.23, 'mmin':4.89, 'mmax':88.81, 'lam':0.04, 'mpp': 33.60, 'sigpp':4.59}
The currently available mass distribution models and their parameters are:
List of mass distribution models
Model Name |
Parameters |
Description |
|---|---|---|
|
Powerlaw + peak model for two-dimensional mass distribution with low mass smoothing. |
|
|
Powerlaw + two peak model for two-dimensional mass distribution with low mass smoothing. |
|
|
Broken power law for two-dimensional mass distribution with low mass smoothing. |
|
|
|
Broken power law + two Gaussian peaks — the fiducial BBH mass model from GWTC-4.0/5.0. See below. |
|
|
PowerLaw + Peak for primary mass and uniform for secondary |
|
|
Truncated Gaussian distribution for primary and secondary |
|
|
Log-normal distribution with mean mu and width sigma for primary and secondary |
|
|
Extension of power law break model |
|
|
Power law with bounds and alpha, spectral index for primary and secondary |
|
|
Both component masses drawn uniformly in |
|
|
Uniform in total mass and mass ratio |
|
|
Full “Power Law + Dip + Break” compact-binary mass function with a pairing function (GWTC-4.0). See below. |
|
|
Tabulated distributions supplied by you, e.g. from a population-synthesis run. See below. |
The parameter definitions follow gwpopulation and the LVK population papers, even though the models are now implemented in-house. It is therefore still worth keeping an eye on how those references define things.
For more details, refer to the following publications:
Note
GWForge overlooks special characters and converts everything to lower cases. So PowerLaw+Peak is equivalent to powerlawpeak.
The BGP model¶
BGP (Broken power law + Gaussian Peaks) is the fiducial BBH mass model used in the GWTC-4.0 and GWTC-5.0 population analyses. The primary mass is a mixture of a broken power law and two left-truncated Gaussian peaks, with a low-mass Planck taper \(S\) applied to the whole distribution:
$$
\pi(m_1) \propto \Big[\lambda_0, p_\mathrm{BP}(m_1) + \lambda_1, N_\mathrm{lt}(m_1 | \mu_1, \sigma_1) + (1 - \lambda_0 - \lambda_1), N_\mathrm{lt}(m_1 | \mu_2, \sigma_2)\Big], S(m_1 | m_\mathrm{min}, \delta_m),
$$
where the broken power law switches from slope \(-\alpha_1\) to \(-\alpha_2\) at \(m_\mathrm{break}\) over \([m_\mathrm{min}, m_\mathrm{high}]\), and \(N_\mathrm{lt}\) is a normal distribution truncated below at \(m_\mathrm{min}\). The mass ratio is a power law with its own low-mass taper on the secondary,
with \(m_\mathrm{min,2}\) and \(\delta_{m,2}\) independent of the primary’s pair (GWTC-5.0 Tab. 5, prior \(m_\mathrm{min,2} \sim U(3, m_\mathrm{min})\) — the secondary generally turns on below the primary). Both default to their primary counterparts if omitted. These are the shipped defaults, so mass-parameters may be left out entirely and you get exactly the block below:
[Mass]
mass-model = BGP
mass-parameters = {'alpha_1': 1.456442737, 'alpha_2': 5.100400428, 'm_break': 37.9806382, 'mmin': 4.489078237, 'delta_m': 3.123416382, 'mmin_2': 3.487749145, 'delta_m_2': 5.60056759, 'm_high': 300.0, 'lam_0': 0.4206178692, 'lam_1': 0.5228817702, 'mpp_1': 9.989384122, 'sigpp_1': 0.6601839103, 'mpp_2': 33.2656028, 'sigpp_2': 4.582221363, 'beta': 0.8049660633, 'maximum_mass': 300}
parameter |
description |
|---|---|
|
Power-law slopes below and above the break |
|
Mass at which the power law breaks |
|
Lower and upper edges of the power-law component (\(m_\mathrm{1,low}\), \(m_\mathrm{high}\)) |
|
Width of the low-mass taper on the primary (\(\delta_{m,1}\)) |
|
Edge and width of the low-mass taper on the secondary. Default to |
|
Mixing fractions of the power law and the first peak (the second peak gets \(1 - \lambda_0 - \lambda_1\)) |
|
Location and width of the first Gaussian peak |
|
Location and width of the second Gaussian peak |
|
Power-law slope of the mass ratio |
Note
Which peak is which is a matter of the values you supply, not of the name: at the O4b medians peak 1 is the narrow low-mass one (9.99 ± 0.66 \(M_\odot\)), carrying 52% of the mixture, and peak 2 is the broad 33.3 \(M_\odot\) one carrying 5.7%. An earlier version of this page used an illustrative example with the two the other way round and the weights almost reversed.
These values are written down once, as
GWForge.population.mass.BGP_PARAMETERS, and read from there by
everything that needs them — the Fisher model and the shipped configs included.
The Default BBH spin values below live the same way, in
GWForge.population.spin.DEFAULT_BBH_SPIN_PARAMETERS.
| maximum_mass | (optional) upper bound of the evaluation grid, default 200 \(M_\odot\). Must exceed m_high |
The exact definitions are Eqs. (B10)–(B14) of the GWTC-5.0 population paper.
User-defined populations¶
Sometimes you already have a population — say, the output of a population-synthesis code — and you simply want GWForge to draw from it. The UserDefined model lets you do exactly that: you hand it a JSON file that tabulates the support xx and the (un-normalised) probability density yy of each parameter, and GWForge turns each into a prior it can sample from. The file looks like this:
{
"mass_1_source": {"xx": [...], "yy": [...]},
"mass_ratio": {"xx": [...], "yy": [...]}
}
You then point the [Mass] section at the file and tell it which key is the primary mass, and either a mass-ratio key or a secondary-mass key:
[Mass]
mass-model = UserDefined
mass-parameters = {'file': '/path/to/priors.json', 'primary_parameter': 'mass_1_source', 'mass_ratio_parameter': 'mass_ratio'}
parameter |
description |
|---|---|
|
Path to the JSON population-priors file |
|
Key to use for the primary mass (default |
|
Key to use for the mass ratio |
|
Key to use for the secondary mass (use this instead of |
Note
If you provide the primary and secondary masses as two independent distributions, GWForge samples them independently and then orders each pair so that \(m_1 \geq m_2\); the pairing of the original population is therefore not preserved. If you care about the pairing, provide the mass ratio instead.
The FullPop_GWTC4 model¶
FullPop_GWTC4 is the strongly-parameterised “Power Law + Dip + Break” model that describes the whole compact-binary mass spectrum — from neutron stars, through the lower mass gap, to black holes — together with a pairing function that sets how strongly binaries favour equal masses. It carries a large number of hyperparameters (the notch edges, filter sharpnesses n0..n5, peak locations and pairing slopes listed in the table above). Because there is no closed-form inverse, the masses are drawn with an importance sampler; the default is importance_m1_m2. The model is described in App. B.2 of the GWTC-5.0 population paper (following Farah et al. 2022 and Mali & Essick 2025).
Spin¶
The [Spin] section determines the spin distribution of the population. For example:
[Spin]
spin-model = Beta-Aligned
spin-parameters = {'minimum_primary_spin' : 0, 'maximum_primary_spin': 0.99, 'minimum_secondary_spin' : 0, 'maximum_secondary_spin' : 0.5, 'mu_chi' : 0.26, 'sigma_squared_chi' : 0.02}
defines a quasi-circular (non-precessing) binary population whose spin magnitude is sampled from a beta distribution.
Here is the list of currently available spin distribution
List of spin distribution models
Model |
Parameters |
Description |
|---|---|---|
|
|
Non-spinning |
|
|
Primary aligned spin drawn from a Truncated Gaussian; secondary non-spinning |
|
|
|
|
|
|
|
|
Aligned component of spins are sampled from uniform distribution |
|
|
Bilby style aligned spin distribution Bilby-style with spin magnitudes obeying Beta distribution |
|
|
Aligned component of primary is sampled from Truncated Gaussian and secondary from uniform |
|
|
Spin Magnitudes sampled from Uniform distribution + Isotropic distribution of spin angles |
|
|
Spin Magnitudes sampled from Beta distribution. Isotropic distribution of spin angles |
|
|
Spin magnitudes sampled from Beta distribution. Truncated Gaussian distribution for cosine tilt angles. |
|
|
Spin magnitudes sampled from Beta distribution. A fraction of the binaries have cosine tilt angles from Truncated Gaussian distribution and the rest from a uniform distribution between (-1,1) |
|
|
The GWTC-4.0/5.0 Default BBH model, Eqs. B15–B16 of arXiv:2605.27226. See below. |
Warning
Default used to be an alias for Isotropic-Beta_Gaussian_Uniform. It is not
the same model, and a config written for the old meaning will now raise rather
than silently sample something else: sigma_squared_chi (a variance, for the
Beta) has become sigma_chi (a standard deviation, for a truncated
Gaussian), and the tilt Gaussian’s mean mu_t is now a parameter instead of
being pinned at \(\cos\theta = 1\). The old model is unchanged and still
available under its own name.
The Default BBH spin model¶
Spin magnitudes are independent and identically distributed truncated Gaussians on \([0, a_{\max}]\) (Eq. B15),
$$\pi(\chi_1, \chi_2 \mid \mu_\chi, \sigma_\chi) = N_{[0, a_{\max}]}(\chi_1 \mid \mu_\chi, \sigma_\chi), N_{[0, a_{\max}]}(\chi_2 \mid \mu_\chi, \sigma_\chi),$$
and the cosine tilts are identically but not independently distributed (Eq. B16),
$$\pi(\cos\theta_1, \cos\theta_2 \mid \mu_t, \sigma_t, \xi) = \xi, N_{[t_{\min}, 1]}(\cos\theta_1 \mid \mu_t, \sigma_t), N_{[t_{\min}, 1]}(\cos\theta_2 \mid \mu_t, \sigma_t)
\frac{1 - \xi}{(1 - t_{\min})^2}.$$
The mixture is over the binary: one draw decides whether both tilts come
from the Gaussian, so the two are correlated. Drawing the split independently
for each component gives identical marginals and the wrong joint, which is why
tests/test_population_spin.py checks the correlation between the two tilts
rather than their histograms.
parameter |
description |
|---|---|
|
Location and width of the spin-magnitude Gaussian, truncated to \([0, a_{\max}]\) |
|
Location and width of the cos-tilt Gaussian. The posterior files call these |
|
Fraction of binaries in the Gaussian component (\(\zeta\) in the paper) |
|
(optional) upper edge of the spin magnitude, default 1 |
|
(optional) lower edge of \(\cos\theta\), default −1 |
These too are the shipped defaults, so a [Spin] section naming the model is enough:
[Spin]
spin-model = Default
spin-parameters = {'mu_chi': 0.0751318233, 'sigma_chi': 0.3667469243, 'mu_t': 0.2788662234, 'sigma_t': 0.9661688072, 'xi_spin': 0.6650168221}
For more details, refer to the following publications:
Extrinsic¶
The [Extrinsic] section handles the sky location and binary orientation. By default GWForge assumes an isotropic sky (right ascension uniform, declination cosine-distributed), isotropic orientation (inclination sine-distributed), and a uniform polarization angle. An empty section is enough to get these defaults:
[Extrinsic]
If you want something other than the defaults, you can point to your own bilby prior file:
[Extrinsic]
extrinsic-prior-file = /path/to/my_extrinsic.prior
You can also override just the inclination without writing a full prior file. Setting inclination-distribution = schutz draws inclinations from the Schutz (2011) distribution, which weights orientations by their detectability — handy if you want a detection-like sample rather than a strictly isotropic one:
[Extrinsic]
inclination-distribution = schutz
parameter |
description |
|---|---|
|
(optional) bilby prior file overriding the default isotropic priors |
|
(optional) set to |
EOS (Equation of State)¶
The [EOS] section allows specifying an eos-file that provides the mass and tidal parameters for a neutron star equation of state. By default, the SLy EOS (Skyrme-Lyon) is used, but you can override this by specifying a different eos-file in the following way:
[EOS]
eos-file = /ligo/home/ligo.org/koustav.chandra/projects/Cosmic-Explorer-MDC/gwforge/GWForge/inject/eos_tables/TOVSeq_SLy.dat
provided it is consistent with how Rahul likes to define them. You can find examples of eos-tables in the GWForge repository.
Structure of eos-file The EOS tables are structured in columns, where each column corresponds to different physical parameters. Below is a breakdown of few of these columns:
C: Compactness of neutron star.
Mb: Baryonic mass in solar masses
M: Mass in solar mass units
R: Radius in solar mass units
kl : Second love number
Generating the population.¶
To generate the binary parameters for the population, execute the following:
gwforge_population --config-file bbh.ini --output-file bbh.h5
It should take at most a minute to generate the output file. By default gwforge_population assumes your source type is BBH. For other options, please check gwforge_population --help. Please note that the waveform approximant that you use for your waveform generation supports tidal parameters if the source-type is bns or nsbh.
If you want the population to be reproducible, pass a --seed:
gwforge_population --config-file bbh.ini --output-file bbh.h5 --seed 42
Running with the same seed and the same configuration gives you byte-for-byte the same population every time. Without a seed, each run gives you a fresh realisation.
Note
By default gwforge_population generates a year’s worth of population. If you want some other value, add the duration option (in seconds) to the [Redshift] section — for example duration = 4096. Please note that the population generated should be greater than the number of signals injected.
A few more example configuration files ship with the package under GWForge/population/population_configuration_files/ (inside your environment’s site-packages, or in the source tree). Feel free to modify them and see what you get.
The one to start from is bgp-gwtc5.ini: the GWTC-5.0 Default BBH population, with every mass and spin value a posterior median of the O4b analysis and a header recording where each came from. It generates about 33,700 BBH mergers a year out to \(z = 10\):
gwforge_population --config-file bgp-gwtc5.ini --output-file bbh.h5 \
--source-type bbh --seed 250114 --save-config
Note
Its [Redshift] block is deliberately not the paper’s. GWForge reads Madau-Dickinson as the star-formation rate and convolves it with a formation-to-merger time delay, so what comes out is a delayed population rather than a rate evolving as \(\psi(z)\) directly.
Naive way to check the population¶
You can check the binary parameters of the population by doing the following:
from GWForge.utils import cornerplot
cornerplot(file='bbh.h5', parameters=['mass_1_source', 'mass_2_source', 'spin_1z','spin_2z', 'redshift'], save='pop.png')
This will create a plot called pop.png in the current working directory with the parameters. The list of parameters can be found by doing h5ls -r bbh.h5. It list all the keys of an HDF5 file.
Checking the population against the model¶
That corner shows you what was drawn. To check it is what the model says, run
python validation/gwtc5_population_check.py --population bbh.h5 \
--config GWForge/population/population_configuration_files/bgp-gwtc5.ini \
--output-directory population_check
It draws an eleven-parameter corner with the analytic density overlaid on every diagonal, and behind each panel runs a numerical comparison: the histogram against the density integrated over each bin, in units of that bin’s counting noise. Two of the checks cannot be seen in any one-dimensional histogram, and both are places this package has had bugs:
\(p(q \mid m_1)\) is a conditional, so it is checked in bins of \(m_1\). Drawing \(q\) from the marginal instead reproduces the marginal perfectly and the conditional not at all.
the cos-tilt mixture is joint over the binary, so it is checked through the correlation between the two tilts. Factorising it leaves both marginals untouched and sets that correlation to zero.