Semi-Analytic Model

This page documents the main public class that drives the package.

SemiAnalyticPopulation

SemiAnalyticPopulation(population_params=None, integration_limits=None, sampling_grids=None, PTA_params=None)

Semi-analytic population model for supermassive black hole binaries.

This is the main high-level interface of :mod:fastropop. It wraps the phenomenological merger-rate density, the cosmology helpers, the gravitational-wave background integrals, and the Monte Carlo sampling tools used to realize discrete SMBHB populations.

The underlying population ansatz is

\[ \frac{d^2 n}{dz\,d\mathcal{M}} \propto \frac{1}{\mathcal{M}\ln 10} \left(\frac{\mathcal{M}}{10^7 M_\odot}\right)^{-\alpha_M} e^{-\mathcal{M}/\mathcal{M}_*} (1+z)^{\beta_z} e^{-z/z_0} \frac{dt_r}{dz}, \]

with free parameters n0, alphaM, Mstar, betaz, and z0. Public methods then derive quantities such as \(h_c^2(f)\), the expected number of binaries, Monte Carlo source realizations, and skymaps.

Initialize one semi-analytic population model.

Parameters:

• population_params (dict, default: None ) –

  Population-law parameters. Supported keys are n0, alphaM, Mstar, betaz, and z0.

• integration_limits (dict, default: None ) –

  Integration bounds with optional keys Mbounds, zbounds, and fbounds.

• sampling_grids (dict, default: None ) –

  Sampling grids with optional keys Mgrid, zgrid, and fgrid. These control the inverse-CDF sampling resolution.

• PTA_params (dict, default: None ) –

  PTA frequency-grid configuration with optional keys Tobs, fmin, fmax, and Nfreqs.

d2ndzdM

d2ndzdM(z, M)

Compute the differential merger-rate density \(d^2n/(dz\,d\mathcal{M})\).

This is the core semi-analytic population law used throughout the package:

\[ \frac{d^2 n}{dz\,d\mathcal{M}} = \frac{n_0}{\mathcal{M}\ln 10} \left(\frac{\mathcal{M}}{10^7 M_\odot}\right)^{-\alpha_M} e^{-\mathcal{M}/\mathcal{M}_*} (1+z)^{\beta_z} e^{-z/z_0} \frac{dt_r}{dz}. \]

Parameters:

• z (float or array - like) –

  Redshift.

• M (float or array - like) –

  Chirp mass in kilograms.

Returns:

• float or Array –

  Differential number density in SI-based units consistent with the rest of the package.

dndlog10M

dndlog10M(M, zmin=0, zmax=5)

Integrate the population over redshift to obtain \(dn/d\log_{10}\mathcal{M}\).

Parameters:

• M (float) –

  Chirp mass in kilograms.

• zmin (float, default: 0 ) –

  Lower redshift integration bound.

• zmax (float, default: 5 ) –

  Upper redshift integration bound.

Returns:

• float –

  Mass function in \({\rm Mpc^{-3}}\).

d3ndzdMdlnf

d3ndzdMdlnf(M, f, z)

Compute \(d^3n/(dz\,d\mathcal{M}\,d\ln f)\).

The extra frequency dimension is obtained by converting the population into time spent per logarithmic frequency interval:

\[ \frac{d^3 n}{dz\,d\mathcal{M}\,d\ln f} = \frac{d^2 n}{dz\,d\mathcal{M}} \left(\frac{d\ln f_r}{dt_r}\right)^{-1} \left(\frac{dt_r}{dz}\right)^{-1} \frac{dV_c}{dz}. \]

hc2

hc2(ff)

Compute the ensemble characteristic strain squared \(h_c^2(f)\).

This method evaluates the smooth background integral implied by the semi-analytic population model. In the current implementation the integration is performed with SciPy using NumPy scalar helpers for efficiency.

Parameters:

• ff (float) –

  Observer-frame frequency in Hz.

Returns:

• float –

  Ensemble characteristic strain squared at frequency ff.

compute_Nbinaries

compute_Nbinaries(Mbounds=None, zbounds=None, fbounds=None, verbose=False)

Integrate the expected number of binaries in the chosen domain.

The integral is carried out over chirp mass, redshift, and logarithmic frequency:

\[ \langle N \rangle = \int d\mathcal{M}\,dz\,d\log_{10}f\; \ln(10)\,\mathcal{M}_{\rm unit}\, \frac{d^3 n}{dz\,d\mathcal{M}\,d\ln f}. \]

Parameters:

• Mbounds (sequence of float, default: None ) –

  Mass bounds [M_min, M_max] in the package sampling convention.

• zbounds (sequence of float, default: None ) –

  Redshift bounds [z_min, z_max].

• fbounds (sequence of float, default: None ) –

  Frequency bounds [f_min, f_max] in SI-scaled package units.

• verbose (bool, default: False ) –

  If True, print the resulting mean binary count.

Returns:

• float –

  Expected number of binaries in the specified domain.

generate_poisson_realization

generate_poisson_realization(Nrealizations=1, Nbinaries_mean=None, key=None, verbose=False)

Draw Poisson realizations of the binary count.

Parameters:

• Nrealizations (int, default: 1 ) –

  Number of Poisson draws to generate.

• Nbinaries_mean (int or float, default: None ) –

  Mean count used as the Poisson expectation value. If omitted, self.Nbinaries_mean is used.

• key (PRNGKey or int, default: None ) –

  Random key or seed for reproducibility. If None, a fresh nondeterministic key is created internally.

• verbose (bool, default: False ) –

  If True, print the realized counts.

Returns:

• Array –

  Integer-valued array of shape (Nrealizations,).

sample_dist

sample_dist(Nbinaries=None, Mgrid=None, zgrid=None, fgrid=None, do_plot=False, key=None)

Sample masses, redshifts, and frequencies from the semi-analytic population.

The sampler builds one-dimensional inverse-CDF approximations for the mass, redshift, and logarithmic-frequency distributions and then draws Nbinaries independent sources.

Parameters:

• Nbinaries (int, default: None ) –

  Number of binaries to sample. If omitted, the method uses int(self.Nbinaries_mean).

• Mgrid (array - like, default: None ) –

  Custom mass grid used to build the inverse-CDF sampler.

• zgrid (array - like, default: None ) –

  Custom redshift grid used to build the inverse-CDF sampler.

• fgrid (array - like, default: None ) –

  Custom frequency grid used to build the inverse-CDF sampler.

• do_plot (bool, default: False ) –

  If True, produce diagnostic sampling plots.

• key (PRNGKey or int, default: None ) –

  JAX random key or integer seed for reproducibility. If None, a fresh nondeterministic key is created internally.

Returns:

• tuple of jax.Array –

  (distM, distz, distlog10f) with shapes (Nbinaries,). distM is returned in the package sampling mass convention, distz in redshift, and distlog10f in \(\log_{10} f\).

generate_skymaps

generate_skymaps(Nbinaries=None, PTA_frequencies=None, Mgrid=None, zgrid=None, fgrid=None, Nside=16, batch_size=int(10000000.0), verbose=False, key=None)

Generate HEALPix skymaps from a sampled realization.

Sources are sampled in batches, assigned isotropic sky positions and binary orientations, converted into plus and cross polarizations, and accumulated into HEALPix maps frequency bin by frequency bin.

Parameters:

• Nbinaries (int, default: None ) –

  Total number of binaries to realize. If omitted, the method uses int(self.Nbinaries_mean).

• PTA_frequencies (array - like, default: None ) –

  Observer-frame PTA frequency-bin centers in Hz.

• Mgrid (array - like, default: None ) –

  Custom mass grid used to build the inverse-CDF sampler.

• zgrid (array - like, default: None ) –

  Custom redshift grid used to build the inverse-CDF sampler.

• fgrid (array - like, default: None ) –

  Custom frequency grid used to build the inverse-CDF sampler.

• Nside (int, default: 16 ) –

  HEALPix resolution parameter.

• batch_size (int, default: int(10000000.0) ) –

  Number of sources processed per batch.

• verbose (bool, default: False ) –

  If True, print batch progress information.

• key (PRNGKey or int, default: None ) –

  JAX random key or integer seed for reproducibility. If None, a fresh nondeterministic key is created internally.

Returns:

• tuple of jax.Array –

  (skymaps_tot, skymaps_plus, skymaps_cross) where skymaps_tot has shape (2, npix, n_frequencies) and the polarization-resolved maps have shape (npix, n_frequencies).

compute_many_realizations

compute_many_realizations(Nbinaries_mean, nrealizations=100, freqs=None, hc2_values=None, do_plot=True, key=None)

Generate many Monte Carlo realizations of the binned GW spectrum.

Each realization first draws a Poisson-distributed binary count with mean Nbinaries_mean, then samples a discrete population and bins it into PTA-style frequency bins.

Parameters:

• Nbinaries_mean (int or float) –

  Poisson mean used for the source-count draw in each realization.

• nrealizations (int, default: 100 ) –

  Number of independent realizations to generate.

• freqs (array - like, default: None ) –

  Optional smooth comparison frequency grid in Hz.

• hc2_values (array - like, default: None ) –

  Optional smooth comparison strain curve passed to the plotting helper.

• do_plot (bool, default: True ) –

  If True, visualize the ensemble of realizations together with the median and central interval.

• key (PRNGKey or int, default: None ) –

  Random key or seed for reproducibility. If None, a fresh nondeterministic key is created internally.

Returns:

• tuple –

  (tabreal, log10f, median, q_low, q_high) where tabreal has shape (nrealizations, n_bins, 2) and stores the transformed spectra used for plotting.