ler.gw_source_population.prior_functions

Module for JIT-compiled functions used in GW source population simulations.

This module provides Numba JIT-compiled functions for efficient sampling and computation of merger rate densities, star formation rates, and mass distributions for various gravitational wave source populations including BBH, BNS, and NSBH.

Key Features:

• Merger rate density functions for PopI/II, PopIII, and Primordial BBH

• Star formation rate models (Madau & Dickinson 2014, Madau & Fragos 2017)

• Mass distribution samplers (power-law Gaussian, broken power-law, bimodal)

• JIT-compiled for high performance with Numba

Copyright (C) 2026 Hemanta Ph. Distributed under MIT License.

Module Contents

Functions

merger_rate_density_bbh_oguri2018_function(zs[, R0, ...])	Compute the merger rate density for PopI/II BBH.
merger_rate_density_bbh_popIII_ken2022_function(zs[, ...])	Compute the unnormalized merger rate density for PopIII BBH.
merger_rate_density_madau_dickinson2014_function(zs[, ...])	Compute the merger rate density for BBH using Madau & Dickinson (2014) model.
merger_rate_density_madau_dickinson_belczynski_ng_function(zs)	Compute BBH merger rate density following Ng et al. (2021).
merger_rate_density_bbh_primordial_ken2022_function(zs)	Compute the merger rate density for Primordial BBH.
sfr_madau_fragos2017_with_bbh_td(zs[, R0])	Compute the merger rate density for BBH. This is computed from star formation rate, Madau & Fragos (2017), with an additional time delay. This function relies on pre-generated data points.
sfr_madau_dickinson2014_with_bbh_td(zs[, R0])	Compute the merger rate density for BBH. This is computed from star formation rate, Madau & Dickinson (2014), with an additional time delay. This function relies on pre-generated data points.
sfr_madau_fragos2017_with_bns_td(zs[, R0])	Compute the merger rate density for BNS. This is computed from star formation rate, Madau & Fragos (2017), with an additional time delay. This function relies on pre-generated data points.
sfr_madau_dickinson2014_with_bns_td(zs[, R0])	Compute the merger rate density for BNS. This is computed from star formation rate, Madau & Dickinson (2014), with an additional time delay. This function relies on pre-generated data points.
sfr_madau_fragos2017(zs[, a, b, c, d])	Compute star formation rate using Madau & Fragos (2017) model.
sfr_madau_dickinson2014(zs[, a, b, c, d])	Compute star formation rate using Madau & Dickinson (2014) model.
ng2022_lognormal_joint_pdf(m1, m2, Mc, sigma)	Compute joint probability density for lognormal 2D mass distribution.
binary_masses_BBH_popIII_lognormal_rvs(size[, m_min, ...])	Generate random samples of binary masses for PopIII BBH (lognormal model).
binary_masses_BBH_primordial_lognormal_rvs(size[, ...])	Generate random samples of binary masses for Primordial BBH (lognormal model).
bimodal_pdf(m[, w, muL, sigmaL, muR, sigmaR, mmin, mmax])	Evaluate fully normalized bimodal Gaussian PDF for BNS mass distribution.
bimodal_cdf(size[, w, muL, sigmaL, muR, sigmaR, mmin, ...])	Compute cumulative distribution function for BNS bimodal mass distribution.
binary_masses_BNS_bimodal_rvs(size[, w, muL, sigmaL, ...])	Generate random samples of binary masses for BNS (bimodal Gaussian model).
broken_powerlaw_pdf(m[, mminbh, mmaxbh, alpha_1, ...])	Compute normalized PDF for broken power-law distribution.
gaussian_plus_isotropic_pdf(x[, mu_t, sigma_t, zeta])	1D marginal PDF for a single cosine tilt x = cos(theta).
gaussian_plus_isotropic_joint_pdf(x1, x2[, mu_t, ...])	2D joint PDF for (x1, x2) = (cos(theta1), cos(theta2)).
powerlaw_pdf(x[, alpha, x_min, x_max])	Compute normalized power-law distribution.
powerlaw_rvs(size, alpha, x_min, x_max)	Inverse transform sampling for a power-law distribution.
truncated_normal_pdf(x, mu, sigma, x_min[, x_max])	Compute left-truncated or left-and-right-truncated normal probability density function.
truncated_normal_rvs(size, mu, sigma, x_min[, x_max])	Draw samples from a truncated normal distribution using analytical inverse CDF.
powerlaw_with_smoothing(q, m, mmin, beta, delta_m)	Compute power-law distribution with low-mass smoothing.
create_gw_parameters_sampler(zs_rvs, m1_rvs, q_rvs, ...)

ler.gw_source_population.prior_functions.merger_rate_density_bbh_oguri2018_function(zs, R0=19 * 1e-09, b2=1.6, b3=2.1, b4=30)

Compute the merger rate density for PopI/II BBH.

Reference: Oguri et al. (2018). The source-frame rate density is

\[R(z) = R_0 \frac{(b_4 + 1) e^{b_2 z}}{b_4 + e^{b_3 z}}.\]

Parameters:

zsfloat or numpy.ndarray

Source redshifts.

R0float

Local merger rate density at low redshift (Mpc^-3 yr^-1).

default: 19e-9 (GWTC-4)

b2float

Fitting parameter.

default: 1.6

b3float

Fitting parameter.

default: 2.1

b4float

Fitting parameter.

default: 30

Returns:

rate_densityfloat or numpy.ndarray

Merger rate density.

Examples

>>> import numpy as np
>>> from ler.gw_source_population import merger_rate_density_bbh_oguri2018_function
>>> rate_density = merger_rate_density_bbh_oguri2018_function(zs=np.array([0.1]))

ler.gw_source_population.prior_functions.merger_rate_density_bbh_popIII_ken2022_function(zs, R0=19.2 * 1e-09, aIII=0.66, bIII=0.3, zIII=11.6)

Compute the unnormalized merger rate density for PopIII BBH.

Reference: Ng et al. (2022). The model is

\[R(z) = R_0 \frac{e^{a_{\rm III}(z-z_{\rm III})}} {b_{\rm III} + a_{\rm III} e^{(a_{\rm III}+b_{\rm III})(z-z_{\rm III})}}.\]

Parameters:

zsfloat or numpy.ndarray

Source redshifts.

R0float

Normalization constant.

default: 19.2e-9

aIIIfloat

Fitting parameter.

default: 0.66

bIIIfloat

Fitting parameter.

default: 0.3

zIIIfloat

Characteristic redshift.

default: 11.6

Returns:

rate_densityfloat or numpy.ndarray

Merger rate density.

Examples

>>> import numpy as np
>>> from ler.gw_source_population import merger_rate_density_bbh_popIII_ken2022_function
>>> rate_density = merger_rate_density_bbh_popIII_ken2022_function(zs=np.array([0.1]))

ler.gw_source_population.prior_functions.merger_rate_density_madau_dickinson2014_function(zs, R0=19 * 1e-09, a=0.015, b=2.7, c=2.9, d=5.6)

Compute the merger rate density for BBH using Madau & Dickinson (2014) model.

The shape follows

\[\psi(z) = a \frac{(1+z)^b}{1 + \left((1+z)/c\right)^d}, \qquad R(z) = R_0 \frac{\psi(z)}{\psi(0)}.\]

Reference: Eq. 15 of https://arxiv.org/pdf/1403.0007

Parameters:

zsfloat or numpy.ndarray

Source redshifts.

R0float

Local merger rate density (Mpc^-3 yr^-1).

default: 19e-9

afloat

Normalization parameter.

default: 0.015

bfloat

Low-redshift power-law slope.

default: 2.7

cfloat

Turnover redshift parameter.

default: 2.9

dfloat

High-redshift power-law slope.

default: 5.6

Returns:

rate_densityfloat or numpy.ndarray

Merger rate density (Mpc^-3 yr^-1).

Examples

>>> import numpy as np
>>> from ler.gw_source_population import merger_rate_density_madau_dickinson2014_function
>>> rate_density = merger_rate_density_madau_dickinson2014_function(zs=np.array([0.1]))

ler.gw_source_population.prior_functions.merger_rate_density_madau_dickinson_belczynski_ng_function(zs, R0=19 * 1e-09, alpha_F=2.57, beta_F=5.83, c_F=3.36)

Compute BBH merger rate density following Ng et al. (2021).

This model uses a Madau-Dickinson-like functional form to fit the merger rate density of field BHs, accounting for time delays and metallicity effects. Coefficients from Madau & Dickinson (2014) are translated as: B-> alpha_F, D-> beta_F, C-> c_F.

\[\psi_F(z) = \frac{(1+z)^{\alpha_F}} {1 + \left((1+z)/c_F\right)^{\beta_F}}, \qquad R(z) = R_0 \frac{\psi_F(z)}{\psi_F(0)}.\]

Parameters:

zsfloat or numpy.ndarray

Source redshifts.

R0float

Local merger rate density (Mpc^-3 yr^-1).

default: 19e-9

alpha_Ffloat

Low-redshift power-law slope.

default: 2.57

beta_Ffloat

High-redshift power-law slope.

default: 5.83

c_Ffloat

Turnover redshift parameter.

default: 3.36

Returns:

rate_densityfloat or numpy.ndarray

Merger rate density (Mpc^-3 yr^-1).

Examples

>>> import numpy as np
>>> from ler.gw_source_population import merger_rate_density_madau_dickinson_belczynski_ng_function
>>> rate_density = merger_rate_density_madau_dickinson_belczynski_ng_function(zs=np.array([0.1]))

ler.gw_source_population.prior_functions.merger_rate_density_bbh_primordial_ken2022_function(zs, cosmology=None, R0=0.044 * 1e-09, t0=13.786885302009708)

Compute the merger rate density for Primordial BBH.

Reference: Ng et al. (2022). The rate is parameterized by cosmic age:

\[R(z) = R_0 \left(\frac{t(z)}{t_0}\right)^{-34/37}.\]

Parameters:

zsfloat or numpy.ndarray

Source redshifts.

cosmologyastropy.cosmology or None

Cosmology object for age calculations.

default: LambdaCDM(H0=70, Om0=0.3, Ode0=0.7, Tcmb0=0.0, Neff=3.04, m_nu=None, Ob0=0.0)

R0float

Normalization constant.

default: 0.044e-9

t0float

Present age of the Universe (Gyr).

default: 13.786885302009708

Returns:

rate_densityfloat or numpy.ndarray

Merger rate density.

Examples

>>> import numpy as np
>>> from ler.gw_source_population import merger_rate_density_bbh_primordial_ken2022_function
>>> rate_density = merger_rate_density_bbh_primordial_ken2022_function(zs=np.array([0.1]))

ler.gw_source_population.prior_functions.sfr_madau_fragos2017_with_bbh_td(zs, R0=19 * 1e-09)

Compute the merger rate density for BBH. This is computed from star formation rate, Madau & Fragos (2017), with an additional time delay. This function relies on pre-generated data points.

Parameters:

zsfloat or numpy.ndarray

Source redshifts.

R0float

Local merger rate density (Mpc^-3 yr^-1).

default: 19e-9

Returns:

SFRfloat or numpy.ndarray

Star formation rate (Mpc^-3 yr^-1).

ler.gw_source_population.prior_functions.sfr_madau_dickinson2014_with_bbh_td(zs, R0=19 * 1e-09)

Compute the merger rate density for BBH. This is computed from star formation rate, Madau & Dickinson (2014), with an additional time delay. This function relies on pre-generated data points.

Parameters:

zsfloat or numpy.ndarray

Source redshifts.

R0float

Local merger rate density (Mpc^-3 yr^-1).

default: 19e-9

Returns:

SFRfloat or numpy.ndarray

Star formation rate (Mpc^-3 yr^-1).

ler.gw_source_population.prior_functions.sfr_madau_fragos2017_with_bns_td(zs, R0=89 * 1e-09)

Compute the merger rate density for BNS. This is computed from star formation rate, Madau & Fragos (2017), with an additional time delay. This function relies on pre-generated data points.

Parameters:

zsfloat or numpy.ndarray

Source redshifts.

R0float

Local merger rate density (Mpc^-3 yr^-1).

default: 89e-9

Returns:

SFRfloat or numpy.ndarray

Star formation rate (Mpc^-3 yr^-1).

ler.gw_source_population.prior_functions.sfr_madau_dickinson2014_with_bns_td(zs, R0=89 * 1e-09)

Compute the merger rate density for BNS. This is computed from star formation rate, Madau & Dickinson (2014), with an additional time delay. This function relies on pre-generated data points.

Parameters:

zsfloat or numpy.ndarray

Source redshifts.

R0float

Local merger rate density (Mpc^-3 yr^-1).

default: 89e-9

Returns:

SFRfloat or numpy.ndarray

Star formation rate (Mpc^-3 yr^-1).

ler.gw_source_population.prior_functions.sfr_madau_fragos2017(zs, a=0.01, b=2.6, c=3.2, d=6.2)

Compute star formation rate using Madau & Fragos (2017) model.

Reference: https://arxiv.org/pdf/1606.07887.pdf

Parameters:

zsfloat or numpy.ndarray

Source redshifts.

afloat

Normalization parameter.

default: 0.01

bfloat

Low-redshift power-law slope.

default: 2.6

cfloat

Turnover redshift parameter.

default: 3.2

dfloat

High-redshift power-law slope.

default: 6.2

Returns:

SFRfloat or numpy.ndarray

Star formation rate (Msun yr^-1 Mpc^-3).

ler.gw_source_population.prior_functions.sfr_madau_dickinson2014(zs, a=0.015, b=2.7, c=2.9, d=5.6)

Compute star formation rate using Madau & Dickinson (2014) model.

Reference: Eq. 15 of https://arxiv.org/pdf/1403.0007

Parameters:

zsfloat or numpy.ndarray

Source redshifts.

afloat

Normalization parameter.

default: 0.015

bfloat

Low-redshift power-law slope.

default: 2.7

cfloat

Turnover redshift parameter.

default: 2.9

dfloat

High-redshift power-law slope.

default: 5.6

Returns:

SFRfloat or numpy.ndarray

Star formation rate (Msun yr^-1 Mpc^-3).

Examples

>>> from ler.gw_source_population import sfr_madau_dickinson2014
>>> sfr = sfr_madau_dickinson2014(zs=np.array([0.1]))

ler.gw_source_population.prior_functions.ng2022_lognormal_joint_pdf(m1, m2, Mc, sigma)

Compute joint probability density for lognormal 2D mass distribution.

Parameters:

m1numpy.ndarray

Primary masses (solar masses).

m2numpy.ndarray

Secondary masses (solar masses).

Mcfloat

Characteristic mass scale (solar masses).

sigmafloat

Width parameter of the lognormal distribution.

Returns:

pdfnumpy.ndarray

Joint probability density values.

ler.gw_source_population.prior_functions.binary_masses_BBH_popIII_lognormal_rvs(size, m_min=1.0, m_max=100.0, Mc=20.0, sigma=0.3, chunk_size=10000)

Generate random samples of binary masses for PopIII BBH (lognormal model).

Draws samples from a 2D lognormal distribution in mass space using rejection sampling. Reference: Ng et al. (2022).

Parameters:

sizeint

Number of binary systems to sample.

m_minfloat

Minimum mass (solar masses). default: 1.0

m_maxfloat

Maximum mass (solar masses). default: 100.0

Mcfloat

Characteristic mass scale (solar masses). default: 20.0

sigmafloat

Width parameter of the lognormal distribution. default: 0.3

chunk_sizeint

Number of samples per rejection sampling chunk. default: 10000

Returns:

m1numpy.ndarray

Primary mass samples (solar masses).

m2numpy.ndarray

Secondary mass samples (solar masses).

Examples

>>> from ler.gw_source_population.prior_functions import binary_masses_BBH_popIII_lognormal_rvs
>>> m1, m2 = binary_masses_BBH_popIII_lognormal_rvs(size=1000)
>>> print(m1.shape, m2.shape)
(1000,) (1000,)

ler.gw_source_population.prior_functions.binary_masses_BBH_primordial_lognormal_rvs(size, m_min=1.0, m_max=100.0, Mc=20.0, sigma=0.3, chunk_size=10000)

Generate random samples of binary masses for Primordial BBH (lognormal model).

Draws samples from a 2D lognormal distribution in mass space using rejection sampling. Reference: Ng et al. (2022).

Parameters:

sizeint

Number of binary systems to sample.

m_minfloat

Minimum mass (solar masses). default: 1.0

m_maxfloat

Maximum mass (solar masses). default: 100.0

Mcfloat

Characteristic mass scale (solar masses). default: 20.0

sigmafloat

Width parameter of the lognormal distribution. default: 0.3

chunk_sizeint

Number of samples per rejection sampling chunk. default: 10000

Returns:

m1numpy.ndarray

Primary mass samples (solar masses).

m2numpy.ndarray

Secondary mass samples (solar masses).

Examples

>>> from ler.gw_source_population.prior_functions import binary_masses_BBH_primordial_lognormal_rvs
>>> m1, m2 = binary_masses_BBH_primordial_lognormal_rvs(size=1000)
>>> print(m1.shape, m2.shape)
(1000,) (1000,)

ler.gw_source_population.prior_functions.bimodal_pdf(m, w=0.643, muL=1.352, sigmaL=0.08, muR=1.88, sigmaR=0.3, mmin=1.0, mmax=2.3)

Evaluate fully normalized bimodal Gaussian PDF for BNS mass distribution.

Computes the normalized probability density function combining two truncated Gaussian peaks. Reference: Will M. Farr et al. 2020.

Parameters:

mnumpy.ndarray

Mass values (solar masses).

wfloat

Weight of the left (low-mass) peak. default: 0.643

muLfloat

Mean of the left peak (solar masses). default: 1.352

sigmaLfloat

Standard deviation of the left peak (solar masses). default: 0.08

muRfloat

Mean of the right peak (solar masses). default: 1.88

sigmaRfloat

Standard deviation of the right peak (solar masses). default: 0.3

mminfloat

Minimum mass (solar masses). default: 1.0

mmaxfloat

Maximum mass (solar masses). default: 2.3

Returns:

pdfnumpy.ndarray

Normalized probability density values.

Examples

>>> from ler.gw_source_population.prior_functions import bimodal_pdf
>>> import numpy as np
>>> m = np.array([1.2, 1.4, 1.8])
>>> pdf_values = bimodal_pdf(m)
>>> print(pdf_values)

ler.gw_source_population.prior_functions.bimodal_cdf(size, w=0.643, muL=1.352, sigmaL=0.08, muR=1.88, sigmaR=0.3, mmin=1.0, mmax=2.3)

Compute cumulative distribution function for BNS bimodal mass distribution.

Parameters:

sizeint

Resolution of the mass grid.

wfloat

Weight of the left (low-mass) peak. default: 0.643

muLfloat

Mean of the left peak (solar masses). default: 1.352

sigmaLfloat

Standard deviation of the left peak (solar masses). default: 0.08

muRfloat

Mean of the right peak (solar masses). default: 1.88

sigmaRfloat

Standard deviation of the right peak (solar masses). default: 0.3

mminfloat

Minimum mass (solar masses). default: 1.0

mmaxfloat

Maximum mass (solar masses). default: 2.3

Returns:

cdf_valuesnumpy.ndarray

Cumulative distribution function values.

mass_gridnumpy.ndarray

Mass grid corresponding to CDF values.

ler.gw_source_population.prior_functions.binary_masses_BNS_bimodal_rvs(size, w=0.643, muL=1.352, sigmaL=0.08, muR=1.88, sigmaR=0.3, mmin=1.0, mmax=2.3, resolution=500)

Generate random samples of binary masses for BNS (bimodal Gaussian model).

Uses inverse transform sampling to draw samples from a bimodal Gaussian distribution of neutron star masses.

Parameters:

sizeint

Number of binary systems to sample.

wfloat

Weight of the left (low-mass) peak. default: 0.643

muLfloat

Mean of the left peak (solar masses). default: 1.352

sigmaLfloat

Standard deviation of the left peak (solar masses). default: 0.08

muRfloat

Mean of the right peak (solar masses). default: 1.88

sigmaRfloat

Standard deviation of the right peak (solar masses). default: 0.3

mminfloat

Minimum mass (solar masses). default: 1.0

mmaxfloat

Maximum mass (solar masses). default: 2.3

resolutionint

Resolution of the mass grid for CDF computation. default: 500

Returns:

mass_samplesnumpy.ndarray

Array of mass samples (solar masses).

ler.gw_source_population.prior_functions.broken_powerlaw_pdf(m, mminbh=26.0, mmaxbh=125.0, alpha_1=6.75, alpha_2=0.0, b=0.5, delta_m=5.0)

Compute normalized PDF for broken power-law distribution.

Parameters:

mnumpy.ndarray

Mass values.

mminbhfloat

Minimum BH mass.

mmaxbhfloat

Maximum BH mass.

alpha_1float

Power-law index below break.

alpha_2float

Power-law index above break.

bfloat

Break location parameter.

delta_mfloat

Smoothing width.

Returns:

pdfnumpy.ndarray

Normalized PDF values.

ler.gw_source_population.prior_functions.gaussian_plus_isotropic_pdf(x, mu_t=0.426, sigma_t=1.222, zeta=0.652)

1D marginal PDF for a single cosine tilt x = cos(theta).

p(x) = zeta * TruncNorm[-1,1](x | mu_t, sigma_t) + (1-zeta)/2

ler.gw_source_population.prior_functions.gaussian_plus_isotropic_joint_pdf(x1, x2, mu_t=0.426, sigma_t=1.222, zeta=0.652)

2D joint PDF for (x1, x2) = (cos(theta1), cos(theta2)).

p(x1, x2) =

zeta * TN(x1) * TN(x2) + (1-zeta) / 4

on [-1,1]^2.

ler.gw_source_population.prior_functions.powerlaw_pdf(x, alpha=-7.7, x_min=1.0, x_max=2.5)

Compute normalized power-law distribution.

p(x) is proportional to x**(-alpha) for x in [x_min, x_max].

Parameters:

xnumpy.ndarray

Input values.

alphafloat

Power-law spectral index.

x_minfloat

Minimum value.

x_maxfloat

Maximum value.

Returns:

pdfnumpy.ndarray

Normalized power-law PDF.

ler.gw_source_population.prior_functions.powerlaw_rvs(size, alpha, x_min, x_max)

Inverse transform sampling for a power-law distribution.

p(x) ∝ x^{-alpha}, x in [x_min, x_max]

Parameters:

sizeint

Number of samples to generate.

alphafloat

Power-law index (alpha).

x_minfloat

Minimum value (lower bound).

x_maxfloat

Maximum value (upper bound).

Returns:

xnumpy.ndarray

Array of sampled values.

ler.gw_source_population.prior_functions.truncated_normal_pdf(x, mu, sigma, x_min, x_max=np.inf)

Compute left-truncated or left-and-right-truncated normal probability density function.

Evaluates the truncated normal distribution $N_{[x_min, x_max]}(x | μ, σ)$, which is a Gaussian distribution with support only between a minimum and maximum value. If x_max is not provided (or set to np.inf), it defaults to left-only truncation.

Parameters:

xnumpy.ndarray

Input values.

mufloat

Mean of the Gaussian distribution.

sigmafloat

Standard deviation of the Gaussian distribution.

x_minfloat

Minimum value (left truncation point).

x_maxfloat, optional

Maximum value (right truncation point). Default is np.inf (no right truncation).

Returns:

pdfnumpy.ndarray

Probability density values, 0 for $x < x_min$ or $x > x_max$.

ler.gw_source_population.prior_functions.truncated_normal_rvs(size, mu, sigma, x_min, x_max=np.inf)

Draw samples from a truncated normal distribution using analytical inverse CDF.

ler.gw_source_population.prior_functions.powerlaw_with_smoothing(q, m, mmin, beta, delta_m)

Compute power-law distribution with low-mass smoothing.

ler.gw_source_population.prior_functions.create_gw_parameters_sampler(zs_rvs, m1_rvs, q_rvs, m2_rvs, tc_rvs, ra_rvs, dec_rvs, phase_rvs, psi_rvs, theta_jn_rvs, a_1_rvs, a_2_rvs, tilt_1_rvs, tilt_2_rvs, phi_12_rvs, phi_jl_rvs, use_njit_sampler=True, spin_zero=False, spin_precession=False)