zeus21.reionization

Models reionization using an analogy of a halo mass function to ionized bubbles.

Authors: zeus21 v2 collaboration - June 2026

Emily Bregou ; Hector Afonso G. Cruz ; Sarah Libanore ; Julian B. Muñoz ; Yonny Sklansky ; Emilie Thélie ; Alessandra Venditti

arXiv:2302.08506, arXiv:2306.09403, arXiv:2407.18294, Sklansky et al. (in prep)

Classes

reionization_global

Computes the global reionization history and bubble mass function (BMF).

Module Contents

class zeus21.reionization.reionization_global(CosmoParams, AstroParams, HMFintclass, z_Init, SFRD_Init, PRINT_SUCCESS=True)[source]

Computes the global reionization history and bubble mass function (BMF).

Parameters:
  • CosmoParams (CosmoParams class)

  • AstroParams (AstroParams class)

  • HMFintclass (HMFinterp class)

  • z_Init (Z_init class)

  • SFRD_Init (SFRD_class class)

  • PRINT_SUCCESS (bool, optional) – Whether to print convergence status for BMF iteration. Default is True.

PRINT_SUCCESS[source]

Whether to print convergence status messages.

Type:

bool

zlist[source]

Redshift array.

Type:

array

Rs[source]

Smoothing radii array matching the rest of the code, in cMpc.

Type:

array

Rs_BMF[source]

Bubble radii only used for the BMF, in cMpc.

Type:

array

ds_array[source]

Sample overdensity values used to compute the ionization barrier.

Type:

array

gamma[source]

Linear coefficient multiplying delta in niondot fit.

Type:

array

gamma2[source]

Quadratic coefficient multiplying delta in niondot fit.

Type:

array

sigma[source]

Matter fluctuation on (zlist, Rs) grid.

Type:

array

gamma_int[source]

Interpolator for gamma as a function of z and log R.

Type:

RegularGridInterpolator

gamma2_int[source]

Interpolator for gamma2 as a function of z and log R.

Type:

RegularGridInterpolator

sigma_BMF[source]

Matter fluctuation evaluated on (zlist, Rs_BMF) grid.

Type:

array

sigma_int[source]

Interpolator for sigma as a function of z and log R.

Type:

RegularGridInterpolator

trec0[source]

Recombination time normalization at z=0.

Type:

float

trec[source]

Recombination time evaluated on zlist, in years.

Type:

array

trec_int[source]

Interpolator for recombination time as a function of z.

Type:

interp1d

niondot_avg[source]

Average ionizing photon production rate as a function of z.

Type:

array

niondot_avg_int[source]

Interpolator for the average ionizing photon production rate.

Type:

interp1d

ion_frac[source]

Global ionized fraction as a function of z.

Type:

array

ion_frac_initial[source]

Initial guess of the global ionized fraction before BMF convergence, from Madau equation.

Type:

array

nion_norm[source]

Normalization factor for niondot fit to match the average niondot, evaluated on (zlist, Rs) grid.

Type:

array

nion_norm_int[source]

Interpolator for nion_norm as a function of z and log R.

Type:

RegularGridInterpolator

prebarrier_xHII[source]

Ionized fraction evaluated on delta, z, and R before computing the barrier.

Type:

array

barrier[source]

Density barrier threshold for ionization as a function of z and R.

Type:

array

barrier_initial[source]

Initial barrier before BMF convergence.

Type:

array

barrier_int[source]

Interpolator for the barrier as a function of z and log R.

Type:

RegularGridInterpolator

prebarrier_xHII_int[source]

Interpolator for prebarrier_xHII as a function of delta, z, and log R.

Type:

RegularGridInterpolator

R_linear_sigma_fit_idx[source]

Index of the smoothing radius closest to AstroParams.R_linear_sigma_fit_input.

Type:

int

R_linear_sigma_fit[source]

Radius used as the initial guess of the peak scale for the linear barrier fit, in cMpc.

Type:

float

BMF[source]

BMF evaluated on (zlist, Rs_BMF) grid.

Type:

array

BMF_initial[source]

Initial BMF based on first guess values before convergence.

Type:

array

peakRofz[source]

Peak bubble radius as a function of z, in cMpc.

Type:

array

peakRofz_int[source]

Interpolator for peak bubble radius as a function of z.

Type:

interp1d

compute_prebarrier_xHII[source]

Computes the ionized fraction before solving for the barrier.

Type:

method

compute_barrier[source]

Computes the density barrier threshold for ionization.

Type:

method

nion_normalization[source]

Computes the normalization factor for the niondot fit.

Type:

method

nrec[source]

Computes the cumulative number of recombinations over delta and z.

Type:

method

niondot_delta_r[source]

Computes the delta and R-dependent ionizing photon production rate.

Type:

method

nion_delta_r_int[source]

Computes the cumulative number of ionizing photons produced since the maximum redshift.

Type:

method

Madau_Q[source]

Computes the global ionized fraction by solving the Madau equation.

Type:

method

B_1[source]

Computes the slope term of the linear ionization barrier.

Type:

method

B_0[source]

Computes the y-intercept term of the linear ionization barrier.

Type:

method

B[source]

Computes the linear ionization barrier as a function of z and R.

Type:

method

dlogsigma_dlogR[source]

Computes dlogsigma/dlogR.

Type:

method

VRdn_dR[source]

Computes the volume-weighted BMF. Integrating this quantity gives the global xHII.

Type:

method

Rdn_dR[source]

Computes the number-weighted BMF.

Type:

method

BMF_peak_R[source]

Finds the radius at which the BMF peaks.

Type:

method

monotonic_after_peak[source]

Enforces monotonic growth of the bubble size peak with time.

Type:

method

analytic_Q[source]

Analytically integrate the BMF to compute global xHII.

Type:

method

converge_BMF[source]

Iteratively updates the ionization barrier, BMF, and xHII until convergence.

Type:

method

interpR[source]

Evaluates any interpolator at fixed z and an array of R.

Type:

method

interpz[source]

Evaluates any interpolator at an array of z and fixed R.

Type:

method

sigmaR_int[source]

Interpolates sigma at fixed z and an array of R.

Type:

method

sigmaz_int[source]

Interpolates sigma at an array of z and fixed R.

Type:

method

barrierR_int[source]

Interpolates the ionization barrier at fixed z and an array of R.

Type:

method

barrierz_int[source]

Interpolates the ionization barrier at an array of z and fixed R.

Type:

method

gammaR_int[source]

Interpolates gamma at at fixed z and an array of R.

Type:

method

gammaz_int[source]

Interpolates gamma at an array of z and fixed R.

Type:

method

gamma2R_int[source]

Interpolates gamma2 at at fixed z and an array of R.

Type:

method

gamma2z_int[source]

Interpolates gamma2 at an array of z and fixed R.

Type:

method

nion_normR_int[source]

Interpolates nion_norm at at fixed z and an array of R.

Type:

method

nion_normz_int[source]

Interpolates nion_norm at an array of z and fixed R.

Type:

method

interp_zR[source]

Evaluates any RegularGridInterpolator on z and R arrays.

Type:

method

sigma_zR_int[source]

Interpolates sigma on z and R arrays.

Type:

method

barrier_zR_int[source]

Interpolates the barrier on z and R arrays.

Type:

method

gamma_zR_int[source]

Interpolates gamma on z and R arrays.

Type:

method

gamma2_zR_int[source]

Interpolates gamma2 on z and R arrays.

Type:

method

nion_norm_zR_int[source]

Interpolates nion_norm on z and R arrays.

Type:

method

prebarrier_xHII_int_grid[source]

Evaluates prebarrier xHII on a delta field at fixed z and R.

Type:

method

PRINT_SUCCESS = True[source]
zlist[source]
Rs[source]
Rs_BMF[source]
ds_array[source]
gamma[source]
gamma2[source]
sigma[source]
zr[source]
gamma_int[source]
gamma2_int[source]
sigma_BMF[source]
zr_BMF[source]
sigma_int[source]
trec0[source]
trec[source]
trec_int[source]
niondot_avg[source]
niondot_avg_int[source]
ion_frac[source]
ion_frac_initial[source]
nion_norm[source]
nion_norm_int[source]
prebarrier_xHII[source]
barrier[source]
barrier_initial[source]
barrier_int[source]
dzr[source]
prebarrier_xHII_int[source]
R_linear_sigma_fit_idx[source]
R_linear_sigma_fit[source]
BMF[source]
peakRofz[source]
peakRofz_int[source]
BMF_initial[source]
compute_prebarrier_xHII(CosmoParams, ion_frac, z, R)[source]

Computes the ionized fraction before solving for the barrier.

Parameters:
  • CosmoParams (CosmoParams class)

  • ion_frac (array) – Global xHII at z.

  • z (array) – Redshifts.

  • R (array) – Smoothing radii in cMpc.

Returns:

prebarrier_xHII – Ionized fraction evaluated over sample delta, z, and R before solving for the barrier.

Return type:

array

compute_barrier(CosmoParams, AstroParams, ion_frac, z, R)[source]

Computes the density barrier threshold for ionization by finding when nion surpasses nH+nrec.

Parameters:
  • CosmoParams (CosmoParams class)

  • AstroParams (AstroParams class)

  • ion_frac (array) – Global xHII at z.

  • z (array) – Redshifts.

  • R (array) – Smoothing radii in cMpc.

Returns:

barrier – Density threshold for ionization as a function of z and R.

Return type:

array

nion_normalization(z, R)[source]

Computes the normalization factor for the niondot fit.

Parameters:
  • z (int or array) – Redshift-grid index or indices.

  • R (int or array) – Radius-grid index or indices.

Returns:

nion_norm – Normalization for niondot.

Return type:

float or array

nrec(CosmoParams, ion_frac, z, d_array=None)[source]

Vectorized computation of nrec over an array of overdensities d_array.

Parameters:
  • CosmoParams (CosmoParams class) – Stores cosmology.

  • ion_frac (1D np.array) – The ionized fraction over all redshifts.

  • z (array) – Redshifts

  • d_array (array, optional) – Sample delta values. Default is None, in which case self.ds_array is used.

Returns:

nrecs – The total number of recombinations at each overdensity for a certain ionized fraction history at each redshift. The first dimension is densities, the second dimension is redshifts.

Return type:

2D np.array

niondot_delta_r(CosmoParams, z, R, d_array=None)[source]

Compute niondot over an array of overdensities d_array for a given R.

Parameters:
  • CosmoParams (CosmoParams class) – Stores cosmology.

  • z (array) – Redshifts

  • R (float) – Radius value (cMpc)

  • d_array (array, optional) – Sample delta values. Default is None, in which case self.ds_array is used.

  • Output

  • ----------

  • niondot (2D np.array) – The rates of ionizing photon production. The first dimension is densities, the second dimension is redshifts.

nion_delta_r_int(CosmoParams, z, R, d_array=None)[source]

Vectorized computation of nion over an array of overdensities d_array for a given R.

Parameters:
  • CosmoParams (CosmoParams class) – Stores cosmology.

  • z (array) – Redshifts

  • R (float) – Radius value (cMpc)

  • d_array (array, optional) – Sample delta values. Default is None, in which case self.ds_array is used.

  • Output

  • ----------

  • nion (2D np.array) – The total number of ionizing photons produced since z=zmax. The first dimension is densities, the second dimension is redshifts.

Madau_Q(CosmoParams, z)[source]

Computes the global ionized fraction by solving the Madau equation.

Parameters:
  • CosmoParams (CosmoParams class)

  • z (float or array) – Redshifts

Returns:

Q – Global xHII evaluated at z.

Return type:

float or array

B_1(z)[source]

Computes the slope term of the linear ionization barrier.

Parameters:

z (float or array) – Redshifts

Returns:

B1 – Slope term of the linear barrier as a function of z.

Return type:

array

B_0(z)[source]

Computes the y-intercept term of the linear ionization barrier.

Parameters:

z (float or array) – Redshifts

Returns:

B0 – Intercept term of the linear barrier as a function of z.

Return type:

array

B(z, R, sig=None)[source]

Computes the linear ionization barrier as a function of z and R.

Parameters:
  • z (float or array) – Redshifts

  • R (float or array) – Radii in cMpc

  • sig (array, optional) – Matter fluctuation sigma(z, R). Default is None, in which case sigma is interpolated internally.

Returns:

B – Linear barrier as a function of z and R.

Return type:

array

dlogsigma_dlogR(z, R, sig)[source]

Computes dlogsigma/dlogR.

Parameters:
  • z (float or array) – Redshifts

  • R (array) – Bubble radii in cMpc.

  • sig (array) – Matter fluctuation sigma(z, R).

Returns:

dlogsigma_dlogR – Logarithmic derivative dlog(sigma)/dlog(R).

Return type:

array

VRdn_dR(z, R)[source]

Computes the volume-weighted BMF. Integrating this quantity gives the global xHII.

Parameters:
  • z (float or array) – Redshifts

  • R (array) – Bubble radii in cMpc.

Returns:

VRdn_dR – Volume-weighted BMF evaluated over redshift and radius.

Return type:

array

Rdn_dR(z, R)[source]

Computes the number-weighted BMF.

Parameters:
  • z (float or array) – Redshifts

  • R (array) – Bubble radii in cMpc.

Returns:

Rdn_dR – Number-weighted bubble mass function evaluated over redshift and radius.

Return type:

array

BMF_peak_R(z, fit_window=5, max_bubble=100, min_bubble=0.5)[source]

Finds the radius at which the BMF peaks.

Parameters:
  • z (float) – Redshift

  • fit_window (int, optional) – Number of points on each side of the coarse peak to use for the spline fit. Default is 5.

  • max_bubble (float, optional) – Maximum allowed bubble radius in cMpc. Default is 100.

  • min_bubble (float, optional) – Minimum allowed bubble radius in cMpc. Default is 0.5.

Returns:

peak_R – Peak bubble radius in cMpc.

Return type:

float

monotonic_after_peak(x)[source]

Enforces monotonic growth of the bubble size peak with time.

Parameters:

x (array) – Input array.

Returns:

x – Array with values left of the peak forced to be non-decreasing.

Return type:

array

analytic_Q(CosmoParams, z)[source]

Analytically integrate the BMF to compute global xHII. Integrates down to some minimum sigma corresponding to the smallest relevant scale for bubbles. Any smaller makes the integral discrepant with the BMF, and the linear barrier becomes a poor fit.

Parameters:
  • CosmoParams (CosmoParams class)

  • z (float or array) – Redshifts

Returns:

Q – Ionized fraction obtained from the analytic integral of the BMF.

Return type:

array

converge_BMF(CosmoParams, AstroParams, ion_frac_input)[source]

Iteratively updates the ionization barrier, BMF, and xHII until convergence.

Parameters:
  • CosmoParams (CosmoParams class)

  • AstroParams (AstroParams class)

  • ion_frac_input (array) – Initial ionized fraction used to begin the BMF convergence loop.

Return type:

None

interpR(z, R, func)[source]

Interpolator to find func(z, R), designed to take a single z but an array of R in cMpc

interpz(z, R, func)[source]

Interpolator to find func(z, R), designed to take a single R in cMpc but an array of z

sigmaR_int(z, R)[source]
sigmaz_int(z, R)[source]
barrierR_int(z, R)[source]
barrierz_int(z, R)[source]
gammaR_int(z, R)[source]
gammaz_int(z, R)[source]
gamma2R_int(z, R)[source]
gamma2z_int(z, R)[source]
nion_normR_int(z, R)[source]
nion_normz_int(z, R)[source]
interp_zR(z, R, func)[source]

Evaluate a RegularGridInterpolator defined on (z, logR).

Accepts scalar, 1D, 2D, or ND z and R. z and R are broadcast against each other.

Examples

scalar z, vector R:

out.shape == R.shape

vector z, scalar R:

out.shape == z.shape

z[:, None], R[None, :]:

out.shape == (nz, nR)

sigma_zR_int(z, R)[source]
barrier_zR_int(z, R)[source]
gamma_zR_int(z, R)[source]
gamma2_zR_int(z, R)[source]
nion_norm_zR_int(z, R)[source]
prebarrier_xHII_int_grid(d, z, R)[source]

Evaluate prebarrier xHII on a density field d(x), at fixed redshift z and smoothing radius R.

Parameters:
  • d (np.ndarray) – Density/overdensity field. Can be any shape (…).

  • z (float) – Redshift.

  • R (float) – Smoothing radius (cMpc).

  • Output

  • ----------

  • values (np.ndarray) – xHII field with the same shape as d.