Core

class gwGenealogy.core.hierarchical.HierarchicalMergersInCluster(Mcl=None, rh=None, v_esc=None, Z=None, stellar_model=None, n_pool=5000, m_pool=None, m_min=None, m_max=None, imf=None, imf_gamma=-2.5, chi_max=0.2, spin_dist='uniform', beta_a=1.4, beta_b=3.6, evolve_v_esc=False, pairing='random', pairing_beta=None, max_gen=5, n_samples=None, precessing=True, mass_spin_model=None, kick_model=None, seed=None)

Bases: object

Monte Carlo simulation of hierarchical BH mergers in a single cluster.

Supports two modes for the escape velocity and two modes for the 1G BH mass pool:

Escape velocity (provide exactly one):

• Physical: Mcl + rh → v_esc computed via virial theorem.

• Direct: v_esc → used as-is.

If both Mcl/rh and v_esc are given, Mcl/rh takes precedence and a warning is printed.

1G BH masses (priority order):

1. m_pool (pre-built array) — used directly.

2. m_min/m_max/imf — simple IMF sampling (no stellar evolution). If Z is also provided it is ignored with a warning.

3. Z (+ optional stellar_model) — Kroupa IMF → stellar evolution → BH remnant masses. If stellar_model is not specified it defaults to 'SEVN_delayed' and a note is printed.

Generation convention: remnant_gen = max(parent_gen1, parent_gen2) + 1.

Parameters:

• Mcl (float or None) – Total cluster mass [Msun].

• rh (float or None) – Half-mass radius [pc].

• v_esc (float or None) – Cluster escape velocity [km/s]. Ignored when Mcl and rh are given.

• Z (float or None) – Metallicity for stellar evolution. Ignored when m_min/m_max/imf or m_pool are provided.

• stellar_model (str or None) – Stellar evolution model: ‘Fryer12_delayed’ or ‘SEVN_delayed’. Defaults to ‘SEVN_delayed’ when Z is used.

• n_pool (int) – Number of ZAMS stars to sample for the 1G BH pool (default: 5000).

• m_pool (array or None) – Pre-built 1G BH mass pool. If provided, skips all mass generation.

• m_min (float or None) – 1G BH mass range [Msun] for simple IMF sampling.

• m_max (float or None) – 1G BH mass range [Msun] for simple IMF sampling.

• imf (str or None) – Initial mass function when using m_min/m_max: ‘uniform’ or ‘powerlaw’ (p(m) ~ m^imf_gamma).

• imf_gamma (float) – Power-law index for imf=’powerlaw’ (default: -2.5).

• chi_max (float) – Maximum natal spin magnitude for 1G BHs (default: 0.2).

• spin_dist (str) – 1G spin magnitude distribution: ‘uniform’ (default) or ‘beta’.

• beta_a (float) – Beta distribution parameters for spin_dist=’beta’ (default: 1.4, 3.6).

• beta_b (float) – Beta distribution parameters for spin_dist=’beta’ (default: 1.4, 3.6).

• evolve_v_esc (bool) – If True, escape velocity decreases as BHs are ejected: v_esc(g) = v_esc_0 * (1 - cumulative_ejected_mass / Mcl)^(1/3). Default: False. Requires Mcl to be set.

• pairing (str) – Binary pairing model: ‘random’ (default), ‘secondary_mass_power_law’, or ‘total_mass_power_law’.

• pairing_beta (float or None) – Power-law index for pairing. None uses defaults.

• max_gen (int) – Maximum generation to simulate (default: 5).

• precessing (bool) – If True (default), use precessing remnant models.

• mass_spin_model (str or None) – Remnant mass/spin model for BBHRemnant (default: ‘hbr’).

• kick_model (str or None) – Kick model for BBHRemnant. Default: ‘gwmodel’ (precessing) or ‘gwmodel_kick_q200’ (nonprecessing).

• seed (int or None) – Master random seed for reproducibility.

simulate(verbose=False)

Run the hierarchical merger simulation across generations.

Parameters:

verbose (bool) – Print generation counts (default: False).

Returns:

dict

Return type:

{generation: {‘m’: array, ‘spin’: array, ‘q’: array}}

plot_generations(data, compare_data=None, label=None, compare_label=None, figsize=None)

Paper-style figure: mass-spin scatter(s) + mass histogram.

Parameters:

• data (dict) – Output from simulate().

• compare_data (dict or None) – Second dataset for comparison.

• label (str or None) – Label for the primary dataset (default: self.kick_model).

• compare_label (str or None) – Label for the comparison dataset.

• figsize (tuple or None) – Figure size.

Return type:

fig, axes

class gwGenealogy.core.hierarchical.HierarchicalMergersInClusterPopulation(n_samples=5000, chi_max=0.2, m_min=3.0, m_max=60.0, imf='uniform', imf_gamma=-2.5, spin_dist='uniform', beta_a=1.4, beta_b=3.6, v_esc_min=1.0, v_esc_max=300.0, v_esc_dist='uniform', v_esc_mean=150.0, v_esc_sigma=45.0, pairing='random', pairing_beta=None, evolve_v_esc=False, v_esc_decay_index=-0.35, max_gen=5, precessing=True, mass_spin_model=None, kick_model=None, seed=None)

Bases: object

Monte Carlo simulation of hierarchical BH mergers across generations.

Generation convention: remnant_gen = max(parent_gen1, parent_gen2) + 1. This means 2G+2G -> 3G, not 4G.

Parameters:

• n_samples (int) – Number of 1G BHs to generate (default: 5000).

• chi_max (float) – Maximum natal spin magnitude for 1G BHs (default: 0.2).

• m_min (float) – 1G BH mass range [Msun] (default: 3.0-60.0).

• m_max (float) – 1G BH mass range [Msun] (default: 3.0-60.0).

• imf (str) – Initial mass function for 1G BH masses: ‘uniform’ (default) or ‘powerlaw’ (p(m) ~ m^imf_gamma).

• imf_gamma (float) – Power-law index when imf=’powerlaw’ (default: -2.5).

• spin_dist (str) – 1G spin magnitude distribution: ‘uniform’ (default) or ‘beta’.

• beta_a (float) – Beta distribution parameters for spin_dist=’beta’ (default: 1.4, 3.6; from arXiv:2111.03634).

• beta_b (float) – Beta distribution parameters for spin_dist=’beta’ (default: 1.4, 3.6; from arXiv:2111.03634).

• v_esc_min (float) – Escape velocity range [km/s] for uniform sampling (default: 1-300).

• v_esc_max (float) – Escape velocity range [km/s] for uniform sampling (default: 1-300).

• v_esc_dist (str) – Escape velocity distribution: ‘uniform’ (default) or ‘gaussian’.

• v_esc_mean (float) – Mean and std for Gaussian v_esc (default: 150, 45 km/s).

• v_esc_sigma (float) – Mean and std for Gaussian v_esc (default: 150, 45 km/s).

• pairing (str) – Binary pairing model: ‘random’ (default), ‘secondary_mass_power_law’ (p(m2|m1) ~ m2^beta, m2 < m1), or ‘total_mass_power_law’ (p(m2) ~ (m1+m2)^beta).

• pairing_beta (float or None) – Power-law index for pairing. None uses defaults: 6.7 for secondary_mass_power_law, 4.0 for total_mass_power_law.

• evolve_v_esc (bool) – If True, escape velocity decays with generation as v_esc(g) = v_esc_sampled * g^(v_esc_decay_index). Default: False.

• v_esc_decay_index (float) – Power-law index for v_esc evolution (default: -0.35). Only used when evolve_v_esc=True.

• max_gen (int) – Maximum generation to simulate (default: 5).

• precessing (bool) – If True (default), use precessing remnant models.

• mass_spin_model (str or None) – Remnant mass/spin model for BBHRemnant (default: ‘hbr’).

• kick_model (str or None) – Kick model for BBHRemnant. Default: ‘gwmodel’ (precessing) or ‘gwmodel_kick_q200’ (nonprecessing).

• seed (int or None) – Master random seed for reproducibility.

simulate(verbose=False)

Run the hierarchical merger simulation across generations.

Parameters:

verbose (bool) – Print generation counts (default: False).

Returns:

dict – Mass, spin, and mass ratio arrays for retained BHs at each generation. Generation 1 contains the initial 1G population (q is empty for 1G).

Return type:

{generation: {‘m’: array, ‘spin’: array, ‘q’: array}}

plot_generations(data, compare_data=None, label=None, compare_label=None, figsize=None)

Paper-style figure: mass-spin scatter(s) + mass histogram.

Parameters:

• data (dict) – Output from simulate().

• compare_data (dict or None) – Second dataset for comparison.

• label (str or None) – Label for the primary dataset (default: self.kick_model).

• compare_label (str or None) – Label for the comparison dataset.

• figsize (tuple or None) – Figure size.

Return type:

fig, axes

class gwGenealogy.core.seed_growth.MonteCarloBHSeedGrowth(v_esc=100.0, Z=0.005, chi_max=0.2, m_seed=10.0, m_targets=None, beta=4.0, max_generations=20, precessing=True, mass_spin_model=None, kick_model=None, evolve_v_esc=False, v_esc_decay_index=-0.35, m_pool=None, n_pool=5000, seed=None)

Bases: object

Monte Carlo simulation of BH seed growth through hierarchical mergers.

Parameters:

• v_esc (float) – Cluster escape velocity [km/s].

• Z (float) – Metallicity (default: 0.005). Controls the 1G BH mass pool via stellar evolution. Ignored if m_pool is provided.

• chi_max (float) – Maximum natal spin magnitude for both seed and 1G BHs (default: 0.2).

• m_seed (float) – Initial seed BH mass [Msun] (default: 10.0).

• m_targets (list of float or None) – Target BBH masses [Msun] to track (default: [100, 250]).

• beta (float) – Pairing steepness: P(partner) ∝ (m_seed + m_partner)^beta (default: 4.0). beta=0 gives random (uniform) pairing.

• max_generations (int) – Maximum number of mergers per chain (default: 20).

• precessing (bool) – If True (default), use precessing remnant models. If False, use nonprecessing (aligned-spin) models.

• mass_spin_model (str or None) – Remnant mass/spin model passed to BBHRemnant. None uses the default (‘hbr’ for precessing, ‘uib’ for nonprecessing).

• evolve_v_esc (bool) – If True, escape velocity decays with generation as v_esc(g) = v_esc * g^(v_esc_decay_index). Default: False.

• v_esc_decay_index (float) – Power-law index for v_esc evolution (default: -0.35). Only used when evolve_v_esc=True.

• kick_model (str) – Kick model passed to BBHRemnant: ‘gwmodel’ (default) or ‘hlz’ for precessing; ‘gwmodel_kick_q200’ (default) or ‘hlz’ for nonprecessing.

• m_pool (array or None) – Pre-built 1G BH mass pool. If provided, skips stellar evolution and uses this array directly. Useful when sweeping over v_esc with the same metallicity.

• n_pool (int) – Number of ZAMS stars to sample for the 1G pool (default: 5000). Ignored if m_pool is provided.

• seed (int or None) – Master random seed for reproducibility.

simulate(n_experiments=10000, store_history=False, verbose=False)

Run the full MC simulation.

Parameters:

• n_experiments (int) – Number of independent seed growth chains (default: 10000).

• store_history (bool) – If True, store per-generation merger details for each chain (default: False).

• verbose (bool) – Print progress every 10% (default: False).

Returns:

final_masses, final_spins, final_generations : arrays retained : bool array reached_target : dict of {m_target: bool array} P_ret : float (retention probability) P_target : dict of {m_target: float} median_mass : float (median final mass of retained seeds) params : dict (simulation parameters) histories : list of list-of-dicts (only if store_history=True)

Return type:

dict with keys

simulate_grid(v_esc_values, n_experiments=10000, verbose=False)

Sweep over escape velocities, reusing the same 1G pool.

Parameters:

• v_esc_values (array-like) – Escape velocities [km/s] to sweep over.

• n_experiments (int) – Number of chains per v_esc value (default: 10000).

• verbose (bool) – Print progress per v_esc value (default: False).

Returns:

v_esc_values : array P_ret : array of retention probabilities P_target : dict of {m_target: array of probabilities} median_mass : array of median retained masses

Return type:

dict with keys

class gwGenealogy.core.seed_growth.IMBHFormationProbability(seed_mass_values, chi_max_values, v_esc_values, m_target=100.0, n_experiments=10000, **mc_kwargs)

Bases: object

IMBH formation probability on a (seed_mass, chi_max, v_esc) grid.

At each grid point, runs a full Monte Carlo seed growth simulation via MonteCarloBHSeedGrowth and records the probability that the seed BH reaches the target mass.

Implements the parameter sweep from Section III.C / Figure 4 of Islam, Wadekar & Kritos (2026, arXiv:2603.10170).

Parameters:

• seed_mass_values (array-like) – Seed BH masses to sweep over [Msun].

• chi_max_values (array-like) – Maximum natal spin values to sweep over.

• v_esc_values (array-like) – Escape velocities to sweep over [km/s].

• m_target (float) – Target IMBH mass [Msun] (default: 100.0).

• n_experiments (int) – Number of MC chains per grid point (default: 10000).

• **mc_kwargs – Additional keyword arguments passed to MonteCarloBHSeedGrowth (e.g., Z, beta, max_generations, kick_model, precessing, n_pool, evolve_v_esc, v_esc_decay_index, seed).

compute(verbose=False)

Run the IMBH probability grid computation.

Returns:

‘p_imbh’ : 3D array (n_seed, n_chi, n_vesc) ‘p_retention’ : 3D array (n_seed, n_chi, n_vesc) ‘seed_mass_values’ : array ‘chi_max_values’ : array ‘v_esc_values’ : array

Return type:

dict with keys

plot_heatmap(v_esc, axes='seed_chi', ax=None, cmap='inferno', vmin=None, vmax=None)

Plot a single P_IMBH heatmap slice at a given v_esc.

Parameters:

• v_esc (float) – Escape velocity value (must be in v_esc_values).

• axes (str) – Which 2D slice: ‘seed_chi’ (default) plots seed_mass vs chi_max.

• ax (matplotlib.axes.Axes or None) – Axes to plot on. If None, creates a new figure.

• cmap (str) – Colormap (default: ‘inferno’).

• vmin (float or None) – Colorbar limits.

• vmax (float or None) – Colorbar limits.

Returns:

fig, ax, im

Return type:

Figure, Axes, and ScalarMappable

plot_heatmap_all_vesc(cmap='inferno', vmin=None, vmax=None, figsize=None)

Multi-panel P_IMBH heatmap across v_esc values.

Parameters:

• cmap (str) – Colormap (default: ‘inferno’).

• vmin (float or None) – Colorbar limits.

• vmax (float or None) – Colorbar limits.

• figsize (tuple or None) – Figure size. Default scales with number of panels.

Returns:

fig, axes

Return type:

Figure and array of Axes

class gwGenealogy.core.retention.BBHRetentionProbabilityOverChiq(q_values, chi_max_values, v_esc_values, kick_models=None, n_samples=10000, spin_dist='uniform', beta_a=1.4, beta_b=3.6, precessing=True, seed=None)

Bases: object

Retention probability grid for BBH mergers over the (chi, q) plane.

Pre-computes shared isotropic spin angles, then sweeps over a (q, chi_max) grid computing P_ret = fraction with kick < v_esc for one or more kick models and v_esc values simultaneously.

Parameters:

• q_values (array-like) – Mass ratio values (q >= 1).

• chi_max_values (array-like) – Maximum spin magnitude values. For spin_dist=’uniform’, spins are drawn from U(0, chi_max). For spin_dist=’beta’, chi_max scales the Beta distribution so spins lie in [0, chi_max].

• v_esc_values (array-like) – Escape velocity values [km/s].

• kick_models (list of str) – Kick model names supported by BBHRemnant: ‘hlz’, ‘gwmodel’, ‘gwmodel_kick_q200’, ‘sur7dq4remnant’, etc. Default: [‘hlz’, ‘gwmodel’].

• n_samples (int) – Number of spin samples per (q, chi_max) grid point (default: 10000).

• spin_dist (str) – Spin magnitude distribution: ‘uniform’ (default) or ‘beta’.

• beta_a (float) – Beta distribution shape parameters (default: 1.4, 3.6). Only used when spin_dist=’beta’.

• beta_b (float) – Beta distribution shape parameters (default: 1.4, 3.6). Only used when spin_dist=’beta’.

• precessing (bool) – If True (default), use precessing kick models.

• seed (int or None) – Random seed for reproducibility.

compute(verbose=False)

Run the retention probability grid computation.

Returns:

‘q_values’ : array ‘chi_max_values’ : array ‘v_esc_values’ : array ‘p_ret’ : dict of {model_name: 3D array (n_q, n_chi, n_vesc)}

Return type:

dict with keys

plot_heatmap(model, v_esc, ax=None, cmap='plasma', vmin=0, vmax=1)

Plot a single P_ret heatmap for one model and v_esc value.

Parameters:

• model (str) – Kick model name (must be in self.kick_models).

• v_esc (float) – Escape velocity value (must be in self.v_esc_values).

• ax (matplotlib.axes.Axes or None) – Axes to plot on. If None, creates a new figure.

• cmap (str) – Colormap (default: ‘plasma’).

• vmin (float) – Colorbar limits (default: 0-1).

• vmax (float) – Colorbar limits (default: 0-1).

Returns:

fig, ax, im

Return type:

Figure, Axes, and ScalarMappable

plot_heatmap_all_vesc(cmap='plasma', vmin=0, vmax=1, figsize=None)

Multi-panel P_ret heatmap: v_esc rows x model columns.

Reproduces Figure 2 layout from Islam, Wadekar & Kritos (2026).

Parameters:

• cmap (str) – Colormap (default: ‘plasma’).

• vmin (float) – Colorbar limits (default: 0-1).

• vmax (float) – Colorbar limits (default: 0-1).

• figsize (tuple or None) – Figure size. Default scales with number of panels.

Returns:

fig, axes

Return type:

Figure and 2D array of Axes

class gwGenealogy.core.retention.BBHRetentionProbabilities(v_esc_values, kick_models=None, n_samples=5000, q_min=1.0, q_max=10.0, q_dist='uniform', q_power=-1.0, q_beta_a=2.0, q_beta_b=5.0, a_min=0.0, a_max=1.0, spin_dist='uniform', spin_beta_a=1.4, spin_beta_b=3.6, spin_angles='isotropic', tilt_beta_a=None, tilt_beta_b=None, precessing=True, seed=None)

Bases: object

Population-level BBH retention probabilities.

Samples a BBH population from specified distributions of mass ratio, spin magnitudes, and spin orientations, then computes kick velocities and retention probabilities as a function of escape velocity.

Unlike BBHRetentionProbabilityOverChiq which sweeps a (q, chi) grid, this class draws from distributions and returns P_ret(v_esc) curves plus the raw kick velocity samples.

Parameters:

• v_esc_values (array-like) – Escape velocity values [km/s] at which to evaluate P_ret.

• kick_models (list of str or None) – Kick model names supported by BBHRemnant. Default: [‘hlz’, ‘gwmodel’].

• n_samples (int) – Number of BBH samples to draw (default: 5000).

• q_min (float) – Mass ratio range (default: 1.0-10.0).

• q_max (float) – Mass ratio range (default: 1.0-10.0).

• q_dist (str) – Mass ratio distribution: ‘uniform’ (default), ‘powerlaw’, or ‘beta’.

• q_power (float) – Power-law index for q_dist=’powerlaw’ (default: -1.0).

• q_beta_a (float) – Beta distribution parameters for q_dist=’beta’ (default: 2.0, 5.0). The raw Beta sample in [0,1] is rescaled to [q_min, q_max].

• q_beta_b (float) – Beta distribution parameters for q_dist=’beta’ (default: 2.0, 5.0). The raw Beta sample in [0,1] is rescaled to [q_min, q_max].

• a_min (float) – Spin magnitude range (default: 0.0-1.0).

• a_max (float) – Spin magnitude range (default: 0.0-1.0).

• spin_dist (str) – Spin magnitude distribution: ‘uniform’ (default) or ‘beta’.

• spin_beta_a (float) – Beta distribution parameters for spin_dist=’beta’ (default: 1.4, 3.6). The raw Beta sample in [0,1] is rescaled to [a_min, a_max].

• spin_beta_b (float) – Beta distribution parameters for spin_dist=’beta’ (default: 1.4, 3.6). The raw Beta sample in [0,1] is rescaled to [a_min, a_max].

• spin_angles (str) – Spin orientation distribution: ‘isotropic’ (default), ‘uniform’, or ‘beta’.

• tilt_beta_a (float or None) – Beta parameters for spin_angles=’beta’ (preferentially aligned).

• tilt_beta_b (float or None) – Beta parameters for spin_angles=’beta’ (preferentially aligned).

• precessing (bool) – If True (default), use precessing kick models.

• seed (int or None) – Random seed for reproducibility.

compute(verbose=False)

Sample the BBH population and compute kick velocities.

Returns:

‘q’ : array of sampled mass ratios ‘a1’, ‘a2’ : arrays of sampled spin magnitudes ‘v_esc_values’ : array ‘kicks’ : dict of {model: array of kick velocities} ‘p_ret’ : dict of {model: array of P_ret per v_esc}

Return type:

dict with keys

plot_kicks(bins=50, ax=None, log=True)

Plot kick velocity distributions for all models.

Parameters:

• bins (int) – Number of histogram bins (default: 50).

• ax (matplotlib.axes.Axes or None) – Axes to plot on. If None, creates a new figure.

• log (bool) – Log scale on x-axis (default: True).

Return type:

fig, ax

plot_retention(ax=None)

Plot P_ret(v_esc) curves for all models.

Parameters:

ax (matplotlib.axes.Axes or None) – Axes to plot on. If None, creates a new figure.

Return type:

fig, ax