zeus21.reionization =================== .. py:module:: zeus21.reionization .. autoapi-nested-parse:: 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 ------- .. autoapisummary:: zeus21.reionization.reionization_global Module Contents --------------- .. py:class:: reionization_global(CosmoParams, AstroParams, HMFintclass, z_Init, SFRD_Init, PRINT_SUCCESS=True) Computes the global reionization history and bubble mass function (BMF). :param CosmoParams: :type CosmoParams: CosmoParams class :param AstroParams: :type AstroParams: AstroParams class :param HMFintclass: :type HMFintclass: HMFinterp class :param z_Init: :type z_Init: Z_init class :param SFRD_Init: :type SFRD_Init: SFRD_class class :param PRINT_SUCCESS: Whether to print convergence status for BMF iteration. Default is True. :type PRINT_SUCCESS: bool, optional .. attribute:: PRINT_SUCCESS Whether to print convergence status messages. :type: bool .. attribute:: zlist Redshift array. :type: array .. attribute:: Rs Smoothing radii array matching the rest of the code, in cMpc. :type: array .. attribute:: Rs_BMF Bubble radii only used for the BMF, in cMpc. :type: array .. attribute:: ds_array Sample overdensity values used to compute the ionization barrier. :type: array .. attribute:: gamma Linear coefficient multiplying delta in niondot fit. :type: array .. attribute:: gamma2 Quadratic coefficient multiplying delta in niondot fit. :type: array .. attribute:: sigma Matter fluctuation on (zlist, Rs) grid. :type: array .. attribute:: gamma_int Interpolator for gamma as a function of z and log R. :type: RegularGridInterpolator .. attribute:: gamma2_int Interpolator for gamma2 as a function of z and log R. :type: RegularGridInterpolator .. attribute:: sigma_BMF Matter fluctuation evaluated on (zlist, Rs_BMF) grid. :type: array .. attribute:: sigma_int Interpolator for sigma as a function of z and log R. :type: RegularGridInterpolator .. attribute:: trec0 Recombination time normalization at z=0. :type: float .. attribute:: trec Recombination time evaluated on zlist, in years. :type: array .. attribute:: trec_int Interpolator for recombination time as a function of z. :type: interp1d .. attribute:: niondot_avg Average ionizing photon production rate as a function of z. :type: array .. attribute:: niondot_avg_int Interpolator for the average ionizing photon production rate. :type: interp1d .. attribute:: ion_frac Global ionized fraction as a function of z. :type: array .. attribute:: ion_frac_initial Initial guess of the global ionized fraction before BMF convergence, from Madau equation. :type: array .. attribute:: nion_norm Normalization factor for niondot fit to match the average niondot, evaluated on (zlist, Rs) grid. :type: array .. attribute:: nion_norm_int Interpolator for nion_norm as a function of z and log R. :type: RegularGridInterpolator .. attribute:: prebarrier_xHII Ionized fraction evaluated on delta, z, and R before computing the barrier. :type: array .. attribute:: barrier Density barrier threshold for ionization as a function of z and R. :type: array .. attribute:: barrier_initial Initial barrier before BMF convergence. :type: array .. attribute:: barrier_int Interpolator for the barrier as a function of z and log R. :type: RegularGridInterpolator .. attribute:: prebarrier_xHII_int Interpolator for prebarrier_xHII as a function of delta, z, and log R. :type: RegularGridInterpolator .. attribute:: R_linear_sigma_fit_idx Index of the smoothing radius closest to AstroParams.R_linear_sigma_fit_input. :type: int .. attribute:: R_linear_sigma_fit Radius used as the initial guess of the peak scale for the linear barrier fit, in cMpc. :type: float .. attribute:: BMF BMF evaluated on (zlist, Rs_BMF) grid. :type: array .. attribute:: BMF_initial Initial BMF based on first guess values before convergence. :type: array .. attribute:: peakRofz Peak bubble radius as a function of z, in cMpc. :type: array .. attribute:: peakRofz_int Interpolator for peak bubble radius as a function of z. :type: interp1d .. attribute:: compute_prebarrier_xHII Computes the ionized fraction before solving for the barrier. :type: method .. attribute:: compute_barrier Computes the density barrier threshold for ionization. :type: method .. attribute:: nion_normalization Computes the normalization factor for the niondot fit. :type: method .. attribute:: nrec Computes the cumulative number of recombinations over delta and z. :type: method .. attribute:: niondot_delta_r Computes the delta and R-dependent ionizing photon production rate. :type: method .. attribute:: nion_delta_r_int Computes the cumulative number of ionizing photons produced since the maximum redshift. :type: method .. attribute:: Madau_Q Computes the global ionized fraction by solving the Madau equation. :type: method .. attribute:: B_1 Computes the slope term of the linear ionization barrier. :type: method .. attribute:: B_0 Computes the y-intercept term of the linear ionization barrier. :type: method .. attribute:: B Computes the linear ionization barrier as a function of z and R. :type: method .. attribute:: dlogsigma_dlogR Computes dlogsigma/dlogR. :type: method .. attribute:: VRdn_dR Computes the volume-weighted BMF. Integrating this quantity gives the global xHII. :type: method .. attribute:: Rdn_dR Computes the number-weighted BMF. :type: method .. attribute:: BMF_peak_R Finds the radius at which the BMF peaks. :type: method .. attribute:: monotonic_after_peak Enforces monotonic growth of the bubble size peak with time. :type: method .. attribute:: analytic_Q Analytically integrate the BMF to compute global xHII. :type: method .. attribute:: converge_BMF Iteratively updates the ionization barrier, BMF, and xHII until convergence. :type: method .. attribute:: interpR Evaluates any interpolator at fixed z and an array of R. :type: method .. attribute:: interpz Evaluates any interpolator at an array of z and fixed R. :type: method .. attribute:: sigmaR_int Interpolates sigma at fixed z and an array of R. :type: method .. attribute:: sigmaz_int Interpolates sigma at an array of z and fixed R. :type: method .. attribute:: barrierR_int Interpolates the ionization barrier at fixed z and an array of R. :type: method .. attribute:: barrierz_int Interpolates the ionization barrier at an array of z and fixed R. :type: method .. attribute:: gammaR_int Interpolates gamma at at fixed z and an array of R. :type: method .. attribute:: gammaz_int Interpolates gamma at an array of z and fixed R. :type: method .. attribute:: gamma2R_int Interpolates gamma2 at at fixed z and an array of R. :type: method .. attribute:: gamma2z_int Interpolates gamma2 at an array of z and fixed R. :type: method .. attribute:: nion_normR_int Interpolates nion_norm at at fixed z and an array of R. :type: method .. attribute:: nion_normz_int Interpolates nion_norm at an array of z and fixed R. :type: method .. attribute:: interp_zR Evaluates any RegularGridInterpolator on z and R arrays. :type: method .. attribute:: sigma_zR_int Interpolates sigma on z and R arrays. :type: method .. attribute:: barrier_zR_int Interpolates the barrier on z and R arrays. :type: method .. attribute:: gamma_zR_int Interpolates gamma on z and R arrays. :type: method .. attribute:: gamma2_zR_int Interpolates gamma2 on z and R arrays. :type: method .. attribute:: nion_norm_zR_int Interpolates nion_norm on z and R arrays. :type: method .. attribute:: prebarrier_xHII_int_grid Evaluates prebarrier xHII on a delta field at fixed z and R. :type: method .. py:attribute:: PRINT_SUCCESS :value: True .. py:attribute:: zlist .. py:attribute:: Rs .. py:attribute:: Rs_BMF .. py:attribute:: ds_array .. py:attribute:: gamma .. py:attribute:: gamma2 .. py:attribute:: sigma .. py:attribute:: zr .. py:attribute:: gamma_int .. py:attribute:: gamma2_int .. py:attribute:: sigma_BMF .. py:attribute:: zr_BMF .. py:attribute:: sigma_int .. py:attribute:: trec0 .. py:attribute:: trec .. py:attribute:: trec_int .. py:attribute:: niondot_avg .. py:attribute:: niondot_avg_int .. py:attribute:: ion_frac .. py:attribute:: ion_frac_initial .. py:attribute:: nion_norm .. py:attribute:: nion_norm_int .. py:attribute:: prebarrier_xHII .. py:attribute:: barrier .. py:attribute:: barrier_initial .. py:attribute:: barrier_int .. py:attribute:: dzr .. py:attribute:: prebarrier_xHII_int .. py:attribute:: R_linear_sigma_fit_idx .. py:attribute:: R_linear_sigma_fit .. py:attribute:: BMF .. py:attribute:: peakRofz .. py:attribute:: peakRofz_int .. py:attribute:: BMF_initial .. py:method:: compute_prebarrier_xHII(CosmoParams, ion_frac, z, R) Computes the ionized fraction before solving for the barrier. :param CosmoParams: :type CosmoParams: CosmoParams class :param ion_frac: Global xHII at z. :type ion_frac: array :param z: Redshifts. :type z: array :param R: Smoothing radii in cMpc. :type R: array :returns: **prebarrier_xHII** -- Ionized fraction evaluated over sample delta, z, and R before solving for the barrier. :rtype: array .. py:method:: compute_barrier(CosmoParams, AstroParams, ion_frac, z, R) Computes the density barrier threshold for ionization by finding when nion surpasses nH+nrec. :param CosmoParams: :type CosmoParams: CosmoParams class :param AstroParams: :type AstroParams: AstroParams class :param ion_frac: Global xHII at z. :type ion_frac: array :param z: Redshifts. :type z: array :param R: Smoothing radii in cMpc. :type R: array :returns: **barrier** -- Density threshold for ionization as a function of z and R. :rtype: array .. py:method:: nion_normalization(z, R) Computes the normalization factor for the niondot fit. :param z: Redshift-grid index or indices. :type z: int or array :param R: Radius-grid index or indices. :type R: int or array :returns: **nion_norm** -- Normalization for niondot. :rtype: float or array .. py:method:: nrec(CosmoParams, ion_frac, z, d_array=None) Vectorized computation of nrec over an array of overdensities d_array. :param CosmoParams: Stores cosmology. :type CosmoParams: CosmoParams class :param ion_frac: The ionized fraction over all redshifts. :type ion_frac: 1D np.array :param z: Redshifts :type z: array :param d_array: Sample delta values. Default is None, in which case self.ds_array is used. :type d_array: array, optional :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. :rtype: 2D np.array .. py:method:: niondot_delta_r(CosmoParams, z, R, d_array=None) Compute niondot over an array of overdensities d_array for a given R. :param CosmoParams: Stores cosmology. :type CosmoParams: CosmoParams class :param z: Redshifts :type z: array :param R: Radius value (cMpc) :type R: float :param d_array: Sample delta values. Default is None, in which case self.ds_array is used. :type d_array: array, optional :param Output: :param ----------: :param niondot: The rates of ionizing photon production. The first dimension is densities, the second dimension is redshifts. :type niondot: 2D np.array .. py:method:: nion_delta_r_int(CosmoParams, z, R, d_array=None) Vectorized computation of nion over an array of overdensities d_array for a given R. :param CosmoParams: Stores cosmology. :type CosmoParams: CosmoParams class :param z: Redshifts :type z: array :param R: Radius value (cMpc) :type R: float :param d_array: Sample delta values. Default is None, in which case self.ds_array is used. :type d_array: array, optional :param Output: :param ----------: :param nion: The total number of ionizing photons produced since z=zmax. The first dimension is densities, the second dimension is redshifts. :type nion: 2D np.array .. py:method:: Madau_Q(CosmoParams, z) Computes the global ionized fraction by solving the Madau equation. :param CosmoParams: :type CosmoParams: CosmoParams class :param z: Redshifts :type z: float or array :returns: **Q** -- Global xHII evaluated at z. :rtype: float or array .. py:method:: B_1(z) Computes the slope term of the linear ionization barrier. :param z: Redshifts :type z: float or array :returns: **B1** -- Slope term of the linear barrier as a function of z. :rtype: array .. py:method:: B_0(z) Computes the y-intercept term of the linear ionization barrier. :param z: Redshifts :type z: float or array :returns: **B0** -- Intercept term of the linear barrier as a function of z. :rtype: array .. py:method:: B(z, R, sig=None) Computes the linear ionization barrier as a function of z and R. :param z: Redshifts :type z: float or array :param R: Radii in cMpc :type R: float or array :param sig: Matter fluctuation sigma(z, R). Default is None, in which case sigma is interpolated internally. :type sig: array, optional :returns: **B** -- Linear barrier as a function of z and R. :rtype: array .. py:method:: dlogsigma_dlogR(z, R, sig) Computes dlogsigma/dlogR. :param z: Redshifts :type z: float or array :param R: Bubble radii in cMpc. :type R: array :param sig: Matter fluctuation sigma(z, R). :type sig: array :returns: **dlogsigma_dlogR** -- Logarithmic derivative dlog(sigma)/dlog(R). :rtype: array .. py:method:: VRdn_dR(z, R) Computes the volume-weighted BMF. Integrating this quantity gives the global xHII. :param z: Redshifts :type z: float or array :param R: Bubble radii in cMpc. :type R: array :returns: **VRdn_dR** -- Volume-weighted BMF evaluated over redshift and radius. :rtype: array .. py:method:: Rdn_dR(z, R) Computes the number-weighted BMF. :param z: Redshifts :type z: float or array :param R: Bubble radii in cMpc. :type R: array :returns: **Rdn_dR** -- Number-weighted bubble mass function evaluated over redshift and radius. :rtype: array .. py:method:: BMF_peak_R(z, fit_window=5, max_bubble=100, min_bubble=0.5) Finds the radius at which the BMF peaks. :param z: Redshift :type z: float :param fit_window: Number of points on each side of the coarse peak to use for the spline fit. Default is 5. :type fit_window: int, optional :param max_bubble: Maximum allowed bubble radius in cMpc. Default is 100. :type max_bubble: float, optional :param min_bubble: Minimum allowed bubble radius in cMpc. Default is 0.5. :type min_bubble: float, optional :returns: **peak_R** -- Peak bubble radius in cMpc. :rtype: float .. py:method:: monotonic_after_peak(x) Enforces monotonic growth of the bubble size peak with time. :param x: Input array. :type x: array :returns: **x** -- Array with values left of the peak forced to be non-decreasing. :rtype: array .. py:method:: analytic_Q(CosmoParams, z) 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. :param CosmoParams: :type CosmoParams: CosmoParams class :param z: Redshifts :type z: float or array :returns: **Q** -- Ionized fraction obtained from the analytic integral of the BMF. :rtype: array .. py:method:: converge_BMF(CosmoParams, AstroParams, ion_frac_input) Iteratively updates the ionization barrier, BMF, and xHII until convergence. :param CosmoParams: :type CosmoParams: CosmoParams class :param AstroParams: :type AstroParams: AstroParams class :param ion_frac_input: Initial ionized fraction used to begin the BMF convergence loop. :type ion_frac_input: array :rtype: None .. py:method:: interpR(z, R, func) Interpolator to find func(z, R), designed to take a single z but an array of R in cMpc .. py:method:: interpz(z, R, func) Interpolator to find func(z, R), designed to take a single R in cMpc but an array of z .. py:method:: sigmaR_int(z, R) .. py:method:: sigmaz_int(z, R) .. py:method:: barrierR_int(z, R) .. py:method:: barrierz_int(z, R) .. py:method:: gammaR_int(z, R) .. py:method:: gammaz_int(z, R) .. py:method:: gamma2R_int(z, R) .. py:method:: gamma2z_int(z, R) .. py:method:: nion_normR_int(z, R) .. py:method:: nion_normz_int(z, R) .. py:method:: interp_zR(z, R, func) Evaluate a RegularGridInterpolator defined on (z, logR). Accepts scalar, 1D, 2D, or ND z and R. z and R are broadcast against each other. .. rubric:: 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) .. py:method:: sigma_zR_int(z, R) .. py:method:: barrier_zR_int(z, R) .. py:method:: gamma_zR_int(z, R) .. py:method:: gamma2_zR_int(z, R) .. py:method:: nion_norm_zR_int(z, R) .. py:method:: prebarrier_xHII_int_grid(d, z, R) Evaluate prebarrier xHII on a density field d(x), at fixed redshift z and smoothing radius R. :param d: Density/overdensity field. Can be any shape (...). :type d: np.ndarray :param z: Redshift. :type z: float :param R: Smoothing radius (cMpc). :type R: float :param Output: :param ----------: :param values: xHII field with the same shape as d. :type values: np.ndarray