BEoRN: Modelling the 21-cm Signal from the Epoch of Reionization

BEoRN: Modelling the 21-cm Signal from the Epoch of Reionization#

This notebook demonstrates how to use BEoRN to model the 21-cm brightness temperature signal from the Epoch of Reionization (EoR). BEoRN paints astrophysical radiation fields — X-ray heating, Lyman-alpha coupling, and ionization — onto a cosmological density field by placing pre-computed 1D profiles around dark matter halos.

The density fields and halo catalogs are generated here using py21cmfast, a semi-numerical simulation code for the EoR. BEoRN then uses these as input to produce full 3D maps of the 21-cm signal.

Requirements: Install BEoRN with the optional dependencies (py21cmfast and hmf):

pip install "git+https://github.com/cosmic-reionization/beorn.git[extra]"

Output: All results are saved under ./output/ and intermediate cached computations under ./cache/. Both directories are created automatically.

import numpy as np
from pathlib import Path
import matplotlib.pyplot as plt
import matplotlib.gridspec as gridspec
%matplotlib inline

import logging
logging.basicConfig(level=logging.INFO)
logger = logging.getLogger(__name__)

FILE_ROOT = Path(".")
CACHE_ROOT = FILE_ROOT / "cache"
OUTPUT_ROOT = FILE_ROOT / "output"

import beorn

import matplotlib as mpl
mpl.rcParams.update({
    'font.size': 14,
    'axes.labelsize': 14,
    'xtick.labelsize': 12,
    'ytick.labelsize': 12,
    'legend.fontsize': 12,
    'axes.titlesize': 14,
})

INFO:numexpr.utils:NumExpr defaulting to 10 threads.

1. Set simulation and astrophysical parameters#

Here we define the cosmology, the simulation volume, and the astrophysical source model.

  • Simulation grid (simulation): box size and grid resolution determine the spatial dynamic range.

  • Cosmo-sim inputs (cosmo_sim): settings specific to py21cmfast — internal resolution factor, random seed, and the snapshot redshifts at which py21cmfast is run and BEoRN paints the signal.

  • Solver (solver): the high-resolution redshift grid over which the 1D RT profiles are integrated, and the halo mass and accretion-rate bins that define the profile lookup table. More profile redshift steps → more accurate bubble growth integration. The snapshot grid can be coarser — BEoRN matches each snapshot to the nearest profile redshift.

  • Cosmology (cosmology): standard ΛCDM parameters. Must be consistent with the py21cmfast run.

  • Source model (source): star formation efficiency, ionizing photon production, X-ray emission, and Lyman-alpha flux. These are the key parameters to vary when exploring different reionization histories.

Tip: Start with a coarse grid (Ncell=64) for quick exploration, then increase resolution for production runs.

parameters = beorn.structs.Parameters()

## Simulation grid
parameters.simulation.cores = 4
parameters.simulation.Lbox  = 100       # comoving Mpc/h
parameters.simulation.Ncell = 128

## Cosmo-sim inputs (py21cmfast-specific)
py21cmfast_high_res_factor = 3
parameters.cosmo_sim.random_seed = 12345
parameters.cosmo_sim.snapshot_redshifts = np.arange(20, 5, -0.50)    # coarser snapshot grid

## Solver: profile redshift grid and halo mass/accretion bins
parameters.solver.redshifts                = np.arange(20, 5, -0.25) # fine profile grid
parameters.solver.halo_mass_bin_min        = 1e7
parameters.solver.halo_mass_bin_max        = 1e15
parameters.solver.halo_mass_nbin           = 40
parameters.solver.halo_mass_accretion_alpha = np.array([0.785, 0.795])

# Source age: how far back to integrate the emission history.
# z_source_start: no sources before this redshift (replaces hardcoded z=35).
# t_source_age: cap the lookback to this many Myr (None = use z_source_start only).
parameters.solver.z_source_start  = 35.0  # default
parameters.source.t_source_age    = None  # default (no age cap)

## Cosmology
parameters.cosmology.Om = 0.315
parameters.cosmology.Ob = 0.049
parameters.cosmology.Ol = 1-0.315
parameters.cosmology.h0  = 0.673

## Source parameters
# X-ray emission
parameters.source.energy_cutoff_min_xray = 500
parameters.source.energy_cutoff_max_xray = 10000
parameters.source.energy_min_sed_xray = 200
parameters.source.energy_max_sed_xray = 10000
parameters.source.alS_xray = 1.5
parameters.source.xray_normalisation = 0.3 * 3.4e40

# Lyman-alpha
parameters.source.n_lyman_alpha_photons = 9690
parameters.source.lyman_alpha_power_law = 0.0

# Ionization
parameters.source.Nion = 2000

# Escape fraction
parameters.source.f0_esc = 0.2
parameters.source.pl_esc = 0.5

# Star formation efficiency
parameters.source.f_st = 0.2
parameters.source.g1 = 0.49
parameters.source.g2 = -0.61
parameters.source.g3 = 4
parameters.source.g4 = -4
parameters.source.Mp = 1.6e11 * parameters.cosmology.h0
parameters.source.Mt = 1e9

# Minimum star-forming halo mass
parameters.source.halo_mass_min = 1e8

Source model visualization#

Before running anything, it is useful to inspect the astrophysical source model implied by the parameters above. The two key efficiency functions are:

  • Star-formation efficiency $f_*(M_h)$: fraction of accreted baryons converted into stars (double power-law with a low-mass suppression below $M_t$).

  • Escape fraction $f_\mathrm{esc}(M_h)$: fraction of ionising photons that escape into the IGM (power-law in halo mass).

These two functions, multiplied together and by $N_\mathrm{ion}$, set the ionising emissivity of each halo.

fig, axes = plt.subplots(1, 2, figsize=(10, 4))

beorn.plotting.draw_star_formation_rate(axes[0], parameters, color='C0')
beorn.plotting.draw_f_esc(axes[1], parameters, color='C1')

# Overlay f_star * f_esc on the left panel to show the combined ionising efficiency
Mh_grid = parameters.solver.halo_mass_bins
keep = (Mh_grid > parameters.source.halo_mass_min) & (Mh_grid <= parameters.source.halo_mass_max)
Mh_plot = Mh_grid[keep]
from beorn.astro import f_star_Halo, f_esc as f_esc_func
combined = f_star_Halo(parameters, Mh_plot) * f_esc_func(parameters, Mh_plot)
axes[0].loglog(Mh_plot, combined, color='C2', ls='--', label=r'$f_* \cdot f_\mathrm{esc}$')
axes[0].legend(fontsize=11)

axes[0].set_title(r'Star-formation efficiency $f_*(M_h)$')
axes[1].set_title(r'Escape fraction $f_\mathrm{esc}(M_h)$')

for ax in axes:
    ax.set_xlabel(r'$M_h \; [M_\odot]$')
    ax.grid(True, which='both', ls=':', alpha=0.4)
    ax.set_xlim(1e8,1e15)
axes[0].set_ylim(1e-3,3e-1)
axes[1].set_ylim(1e-3,1.1)

plt.tight_layout()
plt.show()
../_images/3495db58e2e01d24ee192813e658ed56b7514d65f820499d231b85c3ed9e253c.png

2. Precompute radiation profiles#

BEoRN pre-computes 1D radial profiles of X-ray heating, Lyman-alpha flux, and ionization around halos of different masses. These profiles capture how radiation from the first galaxies affects the surrounding IGM as a function of distance. Pre-computing them once and caching the result means they can be reused efficiently when painting different realisations or astrophysical models.

Note: This step depends only on the astrophysical parameters above — it does not require the cosmological density field or halo catalogs generated in the next section.

cache_handler = beorn.io.Handler(CACHE_ROOT)

solver = beorn.precomputation.RadiationProfileSolver(
    parameters, parameters.cosmo_sim.snapshot_redshifts
)
profiles = solver.get_or_compute_profiles(cache_handler)

INFO:beorn.io.handler:Using persistence directory at cache
INFO:beorn.precomputation.solver:Loaded radiation profiles from cache (z=20.00→5.50, 30 steps).

Inspecting the pre-computed 1D profiles#

The figure below mirrors Figure 2 of Schaeffer et al. 2023. Each panel corresponds to a different halo mass and shows how the radial profile evolves with redshift (blue → red):

  1. Halo mass accretion history — exponential MAR model (grey) vs. Behroozi et al. 2020 simulations (gold).

  2. Lyman-α flux profile $\rho_\alpha(r)$ — sets the WF coupling $x_\alpha$.

  3. Kinetic temperature profile $T_k(r)$ — driven by X-ray heating.

  4. Ionization fraction $x_\mathrm{HII}(r)$ — sharp bubble front expanding with time.

Adjust mass_index and the redshift/alpha lists to explore different halo masses or accretion histories.

# Select a mid-range mass bin and three representative redshifts
# (matching roughly the three panels of Fig. 2 in arXiv:2305.15466)
mass_index = 20   # index into the 40-bin halo mass grid (~10^10–10^11 Msol depending on z)
profile_redshifts = [13.0, 9.5, 6.5]
# Use only the first alpha value to match the number of alpha bins in the cached profiles
profile_alphas = [parameters.solver.halo_mass_accretion_alpha[0]]

beorn.plotting.plot_1D_profiles(
    parameters,
    profiles,
    mass_index=mass_index,
    redshifts=profile_redshifts,
    alphas=profile_alphas,
    label=f"1D radiation profiles — mass bin {mass_index}",
    fontsize=15,
)
../_images/a6e6fd287c7b8fa0fdb0188f7a7c859e4fa7dc32924fe7f5ae83d8f27bb923a3.png

3. Generate the cosmological density field and halo catalogs#

We use py21cmfast to produce the raw cosmological inputs: matter density fields and dark matter halo catalogs at each redshift. These encode where and how massive the first collapsed structures are — the seeds of the first galaxies.

This step only needs to run once for a given cosmology, box size, and random seed. If the data already exists from a previous run it will be reused automatically, so you can safely re-run this cell.

loader = beorn.load_input_data.Py21cmFastLoader(parameters, high_res_factor=py21cmfast_high_res_factor)
output_handler = beorn.io.Handler(OUTPUT_ROOT, input_tag=loader.input_tag)

loader.generate(output_handler)

INFO:beorn.io.handler:Using persistence directory at output
INFO:beorn.load_input_data.cosmo_sim_py21cmfast:py21cmfast setup:
  Output directory : output/py21cmfast_N128_D384_L100_seed12345_afc10c94
  Grid             : HII_DIM=128, DIM=384 (factor 3x)
  Box size         : 100.0 Mpc/h (148.6 Mpc)
  Threads          : 4
  Random seed      : 12345
  Cosmology        : Om=0.315, Ob=0.049, h=0.673, sigma_8=0.83, ns=0.96
INFO:beorn.load_input_data.cosmo_sim_py21cmfast:All 30 snapshots already cached. Skipping generation.

4. Inspect the halo catalog and density field#

Before running the full simulation, it is useful to check the input data. The halo mass function describes the abundance of dark matter halos as a function of mass and should match the analytical prediction (Sheth & Tormen 1999) for the given cosmology and redshift. We also show a 2D slice through the matter overdensity field, which reveals the large-scale cosmic web structure.

z_target = 7.0
z_index = np.argmin(np.abs(loader.redshifts - z_target))

halo_catalog = loader.load_halo_catalog(z_index)
print(f"Number of halos at z={z_target}: {len(halo_catalog.masses)}")
print(f"Mass range: {halo_catalog.masses.min():.2e} — {halo_catalog.masses.max():.2e} Msol")

# Plot halo mass function from the catalog with analytical Sheth & Tormen (1999) reference
fig, ax = plt.subplots()
beorn.plotting.plot_halo_mass_function(ax, halo_catalog, bin_count=15,
                                       label=f'py21cmfast z={z_target}',
                                       analytical_model='ST')
ax.set_title(f'Halo Mass Function at z={z_target}', fontsize=14)
ax.legend()
plt.show()

INFO:beorn.structs.halo_catalog:No alpha values provided, using default value of 0.79 for all halos.

Number of halos at z=7.0: 320931
Mass range: 1.78e+09 — 9.63e+11 Msol
../_images/ce15afdc947bb51fa43bbc9ef937b91eb744be17fae9dad1aa6c0cf562b782cb.png 
delta_m = loader.load_density_field(z_index)

Lbox = parameters.simulation.Lbox  # cMpc/h
extent = [0, Lbox, 0, Lbox]

plt.figure()
plt.imshow(delta_m[:, delta_m.shape[1] // 2, :], origin='lower', cmap='RdBu_r', extent=extent)
plt.colorbar(label=r'$\delta_\mathrm{m}$')
plt.xlabel('[cMpc/h]')
plt.ylabel('[cMpc/h]')
plt.title(f'Matter overdensity at z={z_target}')
plt.show()
../_images/43c17365ecd9d2d4729f2610c4448810c5bc27d30a926218000631f8b12266a2.png

5. Paint the 3D signal maps#

BEoRN places the pre-computed profiles around each halo in the catalog and accumulates their contributions on the grid. This produces 3D maps of the ionization fraction, gas temperature, and Lyman-alpha coupling at each redshift — the ingredients needed to compute the 21-cm brightness temperature.

Snapshots that have already been computed are skipped automatically, so interrupted runs resume from where they left off.

p = beorn.painting.PaintingCoordinator(
    parameters,
    loader=loader,
    output_handler=output_handler,
)

# paint_full() returns immediately — grid fields are lazy.
# Each field (Grid_xHII, Grid_Temp, …) is backed by per-z HDF5 files;
# a slice is read from disk only when first accessed and then cached.
multi_z_quantities = p.paint_full(profiles)

INFO:beorn.painting.coordinator:============================================================
BEoRN model summary
============================================================
  Cosmology   : Om=0.315, Ob=0.049, h0=0.673, sigma_8=0.83
  Grid        : Ncell=128, Lbox=100 Mpc/h
  Profile z   : z=20.0 -> 5.2 (60 steps)
  Snapshot z  : z=20.0 -> 5.5 (30 snapshots)
  1D RT bins  : 1.0e+07 - 1.0e+15 Msun at z=5.2 (40 bins, traced back via exp. accretion)
  Source      : f_st=0.2, Nion=2000, f0_esc=0.2
  X-ray       : norm=1.02e+40, E=[500, 10000] eV
  Lyman-alpha : n_phot=9690, star-forming above 1.0e+08 Msun
  Beorn hash  : f954c658
============================================================
INFO:beorn.painting.coordinator:Painting profiles onto grid for 30 of 30 redshift snapshots. Using 4 processes on a single node.

INFO:beorn.structs.halo_catalog:No alpha values provided, using default value of 0.79 for all halos.
INFO:beorn.painting.coordinator:No halos at z=20.00. Returning empty grids.
INFO:beorn.structs.base_struct:Data written to output/igm_data_py21cmfast_N128_D384_L100_seed12345_afc10c94_f954c658/CoevalCube_z20.000.h5 (48.02 MB)
INFO:beorn.structs.halo_catalog:No alpha values provided, using default value of 0.79 for all halos.
INFO:beorn.painting.coordinator:No halos at z=19.50. Returning empty grids.
INFO:beorn.structs.base_struct:Data written to output/igm_data_py21cmfast_N128_D384_L100_seed12345_afc10c94_f954c658/CoevalCube_z19.500.h5 (48.02 MB)
INFO:beorn.structs.halo_catalog:No alpha values provided, using default value of 0.79 for all halos.
INFO:beorn.painting.coordinator:No halos at z=19.00. Returning empty grids.
INFO:beorn.structs.base_struct:Data written to output/igm_data_py21cmfast_N128_D384_L100_seed12345_afc10c94_f954c658/CoevalCube_z19.000.h5 (48.02 MB)
INFO:beorn.structs.halo_catalog:No alpha values provided, using default value of 0.79 for all halos.
INFO:beorn.painting.coordinator:No halos at z=18.50. Returning empty grids.
INFO:beorn.structs.base_struct:Data written to output/igm_data_py21cmfast_N128_D384_L100_seed12345_afc10c94_f954c658/CoevalCube_z18.500.h5 (48.02 MB)
INFO:beorn.structs.halo_catalog:No alpha values provided, using default value of 0.79 for all halos.
INFO:beorn.painting.coordinator:Painting 2 halos at zgrid=18.00 (z_index=4).
INFO:beorn.painting.coordinator:Profile painting took 0:00:03.467739.
INFO:beorn.painting.spread:Universe not fully ionized : xHII is 0.0.
INFO:beorn.painting.coordinator:Redistributing excess photons from the overlapping regions took 0:00:00.027449.
INFO:beorn.painting.coordinator:Postprocessing of the grids took 0:00:00.041735.
INFO:beorn.painting.coordinator:Current snapshot took 0:00:03.557075.
INFO:beorn.structs.base_struct:Data written to output/igm_data_py21cmfast_N128_D384_L100_seed12345_afc10c94_f954c658/CoevalCube_z18.000.h5 (56.02 MB)
INFO:beorn.structs.halo_catalog:No alpha values provided, using default value of 0.79 for all halos.
INFO:beorn.painting.coordinator:Painting 4 halos at zgrid=17.50 (z_index=5).
INFO:beorn.painting.coordinator:Profile painting took 0:00:03.587410.
INFO:beorn.painting.spread:Universe not fully ionized : xHII is 0.0.
INFO:beorn.painting.coordinator:Redistributing excess photons from the overlapping regions took 0:00:00.028846.
INFO:beorn.painting.coordinator:Postprocessing of the grids took 0:00:00.042925.
INFO:beorn.painting.coordinator:Current snapshot took 0:00:03.676345.
INFO:beorn.structs.base_struct:Data written to output/igm_data_py21cmfast_N128_D384_L100_seed12345_afc10c94_f954c658/CoevalCube_z17.500.h5 (56.02 MB)
INFO:beorn.structs.halo_catalog:No alpha values provided, using default value of 0.79 for all halos.
INFO:beorn.painting.coordinator:Painting 6 halos at zgrid=17.00 (z_index=6).
INFO:beorn.painting.coordinator:Profile painting took 0:00:03.647132.
INFO:beorn.painting.spread:Universe not fully ionized : xHII is 0.0.
INFO:beorn.painting.coordinator:Redistributing excess photons from the overlapping regions took 0:00:00.340325.
INFO:beorn.painting.coordinator:Postprocessing of the grids took 0:00:00.040360.
INFO:beorn.painting.coordinator:Current snapshot took 0:00:04.047246.
INFO:beorn.structs.base_struct:Data written to output/igm_data_py21cmfast_N128_D384_L100_seed12345_afc10c94_f954c658/CoevalCube_z17.000.h5 (56.02 MB)
INFO:beorn.structs.halo_catalog:No alpha values provided, using default value of 0.79 for all halos.
INFO:beorn.painting.coordinator:Painting 10 halos at zgrid=16.50 (z_index=7).
INFO:beorn.painting.coordinator:Profile painting took 0:00:03.644197.
INFO:beorn.painting.spread:Universe not fully ionized : xHII is 0.0.
INFO:beorn.painting.coordinator:Redistributing excess photons from the overlapping regions took 0:00:00.344778.
INFO:beorn.painting.coordinator:Postprocessing of the grids took 0:00:00.040706.
INFO:beorn.painting.coordinator:Current snapshot took 0:00:04.048908.
INFO:beorn.structs.base_struct:Data written to output/igm_data_py21cmfast_N128_D384_L100_seed12345_afc10c94_f954c658/CoevalCube_z16.500.h5 (56.02 MB)
INFO:beorn.structs.halo_catalog:No alpha values provided, using default value of 0.79 for all halos.
INFO:beorn.painting.coordinator:Painting 20 halos at zgrid=16.00 (z_index=8).
INFO:beorn.painting.coordinator:Profile painting took 0:00:04.373759.
INFO:beorn.painting.spread:Universe not fully ionized : xHII is 0.0.
INFO:beorn.painting.coordinator:Redistributing excess photons from the overlapping regions took 0:00:00.335773.
INFO:beorn.painting.coordinator:Postprocessing of the grids took 0:00:00.047285.
INFO:beorn.painting.coordinator:Current snapshot took 0:00:04.778049.
INFO:beorn.structs.base_struct:Data written to output/igm_data_py21cmfast_N128_D384_L100_seed12345_afc10c94_f954c658/CoevalCube_z16.000.h5 (56.02 MB)
INFO:beorn.structs.halo_catalog:No alpha values provided, using default value of 0.79 for all halos.
INFO:beorn.painting.coordinator:Painting 38 halos at zgrid=15.50 (z_index=9).
INFO:beorn.painting.coordinator:Profile painting took 0:00:04.386065.
INFO:beorn.painting.spread:Universe not fully ionized : xHII is 0.0.
INFO:beorn.painting.coordinator:Redistributing excess photons from the overlapping regions took 0:00:00.338932.
INFO:beorn.painting.coordinator:Postprocessing of the grids took 0:00:00.045836.
INFO:beorn.painting.coordinator:Current snapshot took 0:00:04.791760.
INFO:beorn.structs.base_struct:Data written to output/igm_data_py21cmfast_N128_D384_L100_seed12345_afc10c94_f954c658/CoevalCube_z15.500.h5 (56.02 MB)
INFO:beorn.structs.halo_catalog:No alpha values provided, using default value of 0.79 for all halos.
INFO:beorn.painting.coordinator:Painting 72 halos at zgrid=15.00 (z_index=10).
INFO:beorn.painting.coordinator:Profile painting took 0:00:04.227437.
INFO:beorn.painting.spread:Universe not fully ionized : xHII is 0.0.
INFO:beorn.painting.coordinator:Redistributing excess photons from the overlapping regions took 0:00:00.364767.
INFO:beorn.painting.coordinator:Postprocessing of the grids took 0:00:00.046345.
INFO:beorn.painting.coordinator:Current snapshot took 0:00:04.661819.
INFO:beorn.structs.base_struct:Data written to output/igm_data_py21cmfast_N128_D384_L100_seed12345_afc10c94_f954c658/CoevalCube_z15.000.h5 (56.02 MB)
INFO:beorn.structs.halo_catalog:No alpha values provided, using default value of 0.79 for all halos.
INFO:beorn.painting.coordinator:Painting 214 halos at zgrid=14.50 (z_index=11).
INFO:beorn.painting.coordinator:Profile painting took 0:00:04.349535.
INFO:beorn.painting.spread:Universe not fully ionized : xHII is 0.0001.
INFO:beorn.painting.coordinator:Redistributing excess photons from the overlapping regions took 0:00:00.369475.
INFO:beorn.painting.coordinator:Postprocessing of the grids took 0:00:00.046251.
INFO:beorn.painting.coordinator:Current snapshot took 0:00:04.786655.
INFO:beorn.structs.base_struct:Data written to output/igm_data_py21cmfast_N128_D384_L100_seed12345_afc10c94_f954c658/CoevalCube_z14.500.h5 (56.02 MB)
INFO:beorn.structs.halo_catalog:No alpha values provided, using default value of 0.79 for all halos.
INFO:beorn.painting.coordinator:Painting 408 halos at zgrid=14.00 (z_index=12).
INFO:beorn.painting.coordinator:Profile painting took 0:00:04.451275.
INFO:beorn.painting.spread:Universe not fully ionized : xHII is 0.0002.
INFO:beorn.painting.coordinator:Redistributing excess photons from the overlapping regions took 0:00:00.384283.
INFO:beorn.painting.coordinator:Postprocessing of the grids took 0:00:00.046090.
INFO:beorn.painting.coordinator:Current snapshot took 0:00:04.904892.
INFO:beorn.structs.base_struct:Data written to output/igm_data_py21cmfast_N128_D384_L100_seed12345_afc10c94_f954c658/CoevalCube_z14.000.h5 (56.02 MB)
INFO:beorn.structs.halo_catalog:No alpha values provided, using default value of 0.79 for all halos.
INFO:beorn.painting.coordinator:Painting 736 halos at zgrid=13.50 (z_index=13).
INFO:beorn.painting.coordinator:Profile painting took 0:00:04.554303.
INFO:beorn.painting.spread:Universe not fully ionized : xHII is 0.0004.
INFO:beorn.painting.coordinator:Redistributing excess photons from the overlapping regions took 0:00:00.393743.
INFO:beorn.painting.coordinator:Postprocessing of the grids took 0:00:00.046472.
INFO:beorn.painting.coordinator:Current snapshot took 0:00:05.030433.
INFO:beorn.structs.base_struct:Data written to output/igm_data_py21cmfast_N128_D384_L100_seed12345_afc10c94_f954c658/CoevalCube_z13.500.h5 (56.02 MB)
INFO:beorn.structs.halo_catalog:No alpha values provided, using default value of 0.79 for all halos.
INFO:beorn.painting.coordinator:Painting 1286 halos at zgrid=13.00 (z_index=14).
INFO:beorn.painting.coordinator:Profile painting took 0:00:04.833470.
INFO:beorn.painting.spread:Universe not fully ionized : xHII is 0.0007.
INFO:beorn.painting.coordinator:Redistributing excess photons from the overlapping regions took 0:00:00.447453.
INFO:beorn.painting.coordinator:Postprocessing of the grids took 0:00:00.047950.
INFO:beorn.painting.coordinator:Current snapshot took 0:00:05.353263.
INFO:beorn.structs.base_struct:Data written to output/igm_data_py21cmfast_N128_D384_L100_seed12345_afc10c94_f954c658/CoevalCube_z13.000.h5 (56.02 MB)
INFO:beorn.structs.halo_catalog:No alpha values provided, using default value of 0.79 for all halos.
INFO:beorn.painting.coordinator:Painting 2093 halos at zgrid=12.50 (z_index=15).
INFO:beorn.painting.coordinator:Profile painting took 0:00:04.536998.
INFO:beorn.painting.spread:Universe not fully ionized : xHII is 0.0013.
INFO:beorn.painting.coordinator:Redistributing excess photons from the overlapping regions took 0:00:00.464412.
INFO:beorn.painting.coordinator:Postprocessing of the grids took 0:00:00.047941.
INFO:beorn.painting.coordinator:Current snapshot took 0:00:05.070357.
INFO:beorn.structs.base_struct:Data written to output/igm_data_py21cmfast_N128_D384_L100_seed12345_afc10c94_f954c658/CoevalCube_z12.500.h5 (56.02 MB)
INFO:beorn.structs.halo_catalog:No alpha values provided, using default value of 0.79 for all halos.
INFO:beorn.painting.coordinator:Painting 3442 halos at zgrid=12.00 (z_index=16).
INFO:beorn.painting.coordinator:Profile painting took 0:00:04.834804.
INFO:beorn.painting.spread:Universe not fully ionized : xHII is 0.0025.
INFO:beorn.painting.coordinator:Redistributing excess photons from the overlapping regions took 0:00:00.553181.
INFO:beorn.painting.coordinator:Postprocessing of the grids took 0:00:00.046880.
INFO:beorn.painting.coordinator:Current snapshot took 0:00:05.454581.
INFO:beorn.structs.base_struct:Data written to output/igm_data_py21cmfast_N128_D384_L100_seed12345_afc10c94_f954c658/CoevalCube_z12.000.h5 (56.02 MB)
INFO:beorn.structs.halo_catalog:No alpha values provided, using default value of 0.79 for all halos.
INFO:beorn.painting.coordinator:Painting 7204 halos at zgrid=11.50 (z_index=17).
INFO:beorn.painting.coordinator:Profile painting took 0:00:04.731901.
INFO:beorn.painting.spread:Universe not fully ionized : xHII is 0.0043.
INFO:beorn.painting.coordinator:Redistributing excess photons from the overlapping regions took 0:00:00.689226.
INFO:beorn.painting.coordinator:Postprocessing of the grids took 0:00:00.049123.
INFO:beorn.painting.coordinator:Current snapshot took 0:00:05.493695.
INFO:beorn.structs.base_struct:Data written to output/igm_data_py21cmfast_N128_D384_L100_seed12345_afc10c94_f954c658/CoevalCube_z11.500.h5 (56.02 MB)
INFO:beorn.structs.halo_catalog:No alpha values provided, using default value of 0.79 for all halos.
INFO:beorn.painting.coordinator:Painting 11407 halos at zgrid=11.00 (z_index=18).
INFO:beorn.painting.coordinator:Profile painting took 0:00:04.729517.
INFO:beorn.painting.spread:Universe not fully ionized : xHII is 0.0084.
INFO:beorn.painting.coordinator:Redistributing excess photons from the overlapping regions took 0:00:01.086798.
INFO:beorn.painting.coordinator:Postprocessing of the grids took 0:00:00.045408.
INFO:beorn.painting.coordinator:Current snapshot took 0:00:05.882654.
INFO:beorn.structs.base_struct:Data written to output/igm_data_py21cmfast_N128_D384_L100_seed12345_afc10c94_f954c658/CoevalCube_z11.000.h5 (56.02 MB)
INFO:beorn.structs.halo_catalog:No alpha values provided, using default value of 0.79 for all halos.
INFO:beorn.painting.coordinator:Painting 17795 halos at zgrid=10.50 (z_index=19).
INFO:beorn.painting.coordinator:Profile painting took 0:00:04.890882.
INFO:beorn.painting.spread:Universe not fully ionized : xHII is 0.0141.
INFO:beorn.painting.coordinator:Redistributing excess photons from the overlapping regions took 0:00:01.475744.
INFO:beorn.painting.coordinator:Postprocessing of the grids took 0:00:00.047131.
INFO:beorn.painting.coordinator:Current snapshot took 0:00:06.440837.
INFO:beorn.structs.base_struct:Data written to output/igm_data_py21cmfast_N128_D384_L100_seed12345_afc10c94_f954c658/CoevalCube_z10.500.h5 (56.02 MB)
INFO:beorn.structs.halo_catalog:No alpha values provided, using default value of 0.79 for all halos.
INFO:beorn.painting.coordinator:Painting 26895 halos at zgrid=10.00 (z_index=20).
INFO:beorn.painting.coordinator:Profile painting took 0:00:05.002832.
INFO:beorn.painting.spread:Universe not fully ionized : xHII is 0.0256.
INFO:beorn.painting.coordinator:Redistributing excess photons from the overlapping regions took 0:00:02.433860.
INFO:beorn.painting.coordinator:Postprocessing of the grids took 0:00:00.047754.
INFO:beorn.painting.coordinator:Current snapshot took 0:00:07.507069.
INFO:beorn.structs.base_struct:Data written to output/igm_data_py21cmfast_N128_D384_L100_seed12345_afc10c94_f954c658/CoevalCube_z10.000.h5 (56.02 MB)
INFO:beorn.structs.halo_catalog:No alpha values provided, using default value of 0.79 for all halos.
INFO:beorn.painting.coordinator:Painting 48926 halos at zgrid=9.50 (z_index=21).
INFO:beorn.painting.coordinator:Profile painting took 0:00:04.787054.
INFO:beorn.painting.spread:Universe not fully ionized : xHII is 0.0368.
INFO:beorn.painting.coordinator:Redistributing excess photons from the overlapping regions took 0:00:03.708379.
INFO:beorn.painting.coordinator:Postprocessing of the grids took 0:00:00.046448.
INFO:beorn.painting.coordinator:Current snapshot took 0:00:08.568313.
INFO:beorn.structs.base_struct:Data written to output/igm_data_py21cmfast_N128_D384_L100_seed12345_afc10c94_f954c658/CoevalCube_z9.500.h5 (56.02 MB)
INFO:beorn.structs.halo_catalog:No alpha values provided, using default value of 0.79 for all halos.
INFO:beorn.painting.coordinator:Painting 72344 halos at zgrid=9.00 (z_index=22).
INFO:beorn.painting.coordinator:Profile painting took 0:00:04.770899.
INFO:beorn.painting.spread:Universe not fully ionized : xHII is 0.0704.
INFO:beorn.painting.coordinator:Redistributing excess photons from the overlapping regions took 0:00:06.769397.
INFO:beorn.painting.coordinator:Postprocessing of the grids took 0:00:00.046608.
INFO:beorn.painting.coordinator:Current snapshot took 0:00:11.612958.
INFO:beorn.structs.base_struct:Data written to output/igm_data_py21cmfast_N128_D384_L100_seed12345_afc10c94_f954c658/CoevalCube_z9.000.h5 (56.02 MB)
INFO:beorn.structs.halo_catalog:No alpha values provided, using default value of 0.79 for all halos.
INFO:beorn.painting.coordinator:Painting 105154 halos at zgrid=8.50 (z_index=23).
INFO:beorn.painting.coordinator:Profile painting took 0:00:04.893232.
INFO:beorn.painting.spread:Universe not fully ionized : xHII is 0.116.
INFO:beorn.painting.coordinator:Redistributing excess photons from the overlapping regions took 0:00:10.761899.
INFO:beorn.painting.coordinator:Postprocessing of the grids took 0:00:00.046014.
INFO:beorn.painting.coordinator:Current snapshot took 0:00:15.726640.
INFO:beorn.structs.base_struct:Data written to output/igm_data_py21cmfast_N128_D384_L100_seed12345_afc10c94_f954c658/CoevalCube_z8.500.h5 (56.02 MB)
INFO:beorn.structs.halo_catalog:No alpha values provided, using default value of 0.79 for all halos.
INFO:beorn.painting.coordinator:Painting 146176 halos at zgrid=8.00 (z_index=24).
INFO:beorn.painting.coordinator:Profile painting took 0:00:05.034174.
INFO:beorn.painting.spread:Universe not fully ionized : xHII is 0.1763.
INFO:beorn.painting.coordinator:Redistributing excess photons from the overlapping regions took 0:00:13.271135.
INFO:beorn.painting.coordinator:Postprocessing of the grids took 0:00:00.045330.
INFO:beorn.painting.coordinator:Current snapshot took 0:00:18.377290.
INFO:beorn.structs.base_struct:Data written to output/igm_data_py21cmfast_N128_D384_L100_seed12345_afc10c94_f954c658/CoevalCube_z8.000.h5 (56.02 MB)
INFO:beorn.structs.halo_catalog:No alpha values provided, using default value of 0.79 for all halos.
INFO:beorn.painting.coordinator:Painting 235134 halos at zgrid=7.50 (z_index=25).
INFO:beorn.painting.coordinator:Profile painting took 0:00:05.070835.
INFO:beorn.painting.spread:Universe not fully ionized : xHII is 0.2657.
INFO:beorn.painting.coordinator:Redistributing excess photons from the overlapping regions took 0:00:13.043184.
INFO:beorn.painting.coordinator:Postprocessing of the grids took 0:00:00.044719.
INFO:beorn.painting.coordinator:Current snapshot took 0:00:18.187451.
INFO:beorn.structs.base_struct:Data written to output/igm_data_py21cmfast_N128_D384_L100_seed12345_afc10c94_f954c658/CoevalCube_z7.500.h5 (56.02 MB)
INFO:beorn.structs.halo_catalog:No alpha values provided, using default value of 0.79 for all halos.
INFO:beorn.painting.coordinator:Painting 320931 halos at zgrid=7.00 (z_index=26).
INFO:beorn.painting.coordinator:Profile painting took 0:00:05.265423.
INFO:beorn.painting.spread:Universe not fully ionized : xHII is 0.4227.
INFO:beorn.painting.coordinator:Redistributing excess photons from the overlapping regions took 0:00:08.699281.
INFO:beorn.painting.coordinator:Postprocessing of the grids took 0:00:00.045201.
INFO:beorn.painting.coordinator:Current snapshot took 0:00:14.035967.
INFO:beorn.structs.base_struct:Data written to output/igm_data_py21cmfast_N128_D384_L100_seed12345_afc10c94_f954c658/CoevalCube_z7.000.h5 (56.02 MB)
INFO:beorn.structs.halo_catalog:No alpha values provided, using default value of 0.79 for all halos.
INFO:beorn.painting.coordinator:Painting 419899 halos at zgrid=6.50 (z_index=27).
INFO:beorn.painting.coordinator:Profile painting took 0:00:05.201960.
INFO:beorn.painting.spread:Universe not fully ionized : xHII is 0.6667.
INFO:beorn.painting.coordinator:Redistributing excess photons from the overlapping regions took 0:00:06.623138.
INFO:beorn.painting.coordinator:Postprocessing of the grids took 0:00:00.040906.
INFO:beorn.painting.coordinator:Current snapshot took 0:00:11.902025.
INFO:beorn.structs.base_struct:Data written to output/igm_data_py21cmfast_N128_D384_L100_seed12345_afc10c94_f954c658/CoevalCube_z6.500.h5 (56.02 MB)
INFO:beorn.structs.halo_catalog:No alpha values provided, using default value of 0.79 for all halos.
INFO:beorn.painting.coordinator:Painting 616951 halos at zgrid=6.00 (z_index=28).
INFO:beorn.painting.coordinator:Profile painting took 0:00:05.323217.
INFO:beorn.painting.spread:Universe not fully ionized : xHII is 0.9955.
INFO:beorn.painting.coordinator:Redistributing excess photons from the overlapping regions took 0:00:27.274100.
INFO:beorn.painting.coordinator:Postprocessing of the grids took 0:00:00.036055.
INFO:beorn.painting.coordinator:Current snapshot took 0:00:32.671692.
INFO:beorn.structs.base_struct:Data written to output/igm_data_py21cmfast_N128_D384_L100_seed12345_afc10c94_f954c658/CoevalCube_z6.000.h5 (56.02 MB)
INFO:beorn.structs.halo_catalog:No alpha values provided, using default value of 0.79 for all halos.
INFO:beorn.painting.coordinator:Painting 787230 halos at zgrid=5.50 (z_index=29).
INFO:beorn.painting.coordinator:Profile painting took 0:00:05.896070.
INFO:beorn.painting.coordinator:Redistributing excess photons from the overlapping regions took 0:00:00.078394.
INFO:beorn.painting.coordinator:Postprocessing of the grids took 0:00:00.035820.
INFO:beorn.painting.coordinator:Current snapshot took 0:00:06.050421.
INFO:beorn.structs.base_struct:Data written to output/igm_data_py21cmfast_N128_D384_L100_seed12345_afc10c94_f954c658/CoevalCube_z5.500.h5 (56.02 MB)
INFO:beorn.painting.coordinator:Painting of 30 snapshots done.
INFO:beorn.structs.temporal_cube:Opened 30 snapshots from output/igm_data_py21cmfast_N128_D384_L100_seed12345_afc10c94_f954c658 (grid data is lazy — slices load on demand).

6. Visualise the 3D maps at a single redshift#

We show a 2D slice through the simulation box for the three key fields at the chosen redshift: the ionization fraction $x_\mathrm{HII}$, the gas kinetic temperature $T_k$, and the 21-cm brightness temperature $\delta T_b$. The ionized bubbles visible in the $x_\mathrm{HII}$ map are driven by the same halos that appear in the mass function above.

# Load only the single snapshot we want to visualise.
# snapshot() accepts a redshift value (nearest match) or an integer index,
# and opens only that one HDF5 file — minimal memory use.
snap = multi_z_quantities.snapshot(z_target)

xHII_grid = snap.Grid_xHII[:]  # 3D array: (Ncell, Ncell, Ncell)
dTb_grid  = snap.Grid_dTb      # derived from this snapshot only
Tk_grid   = snap.Grid_Temp[:]

mid = xHII_grid.shape[1] // 2
Lbox = parameters.simulation.Lbox  # cMpc/h
extent = [0, Lbox, 0, Lbox]

fig, axs = plt.subplots(1, 3, figsize=(14, 4))
for ax, grid, label in zip(axs,
                            [xHII_grid, Tk_grid, dTb_grid],
                            [r'$x_\mathrm{HII}$', r'$T_k$ [K]', r'$\delta T_b$ [mK]']):
    im = ax.imshow(grid[:, mid, :], origin='lower', cmap='viridis', extent=extent)
    ax.set_xlabel('[cMpc/h]', fontsize=15)
    fig.colorbar(im, ax=ax, label=label)
axs[0].set_ylabel('[cMpc/h]', fontsize=15)
fig.suptitle(f'$z={snap.z:.3f}$', fontsize=15)
plt.tight_layout()
plt.show()
../_images/6f9a2ed99da0b0138935ced2e85d18e26bfec3cdc47bc306936079971e4817ee.png

7. Global history and power spectra#

The global (volume-averaged) signal summarises the overall evolution of reionization. We show the redshift evolution of the mean brightness temperature, gas temperature, and ionization fraction, as well as the 21-cm power spectrum amplitude at a fixed wavenumber $k$. These are the primary observational targets of radio telescopes such as SKA, HERA, and MWA.

Memory note: The draw_* plotting functions use global_mean() internally, which iterates snapshots one at a time so peak memory is one 3D grid. The power spectrum panel requires Grid_dTb across all $z$ and will load all base fields. For a memory-efficient global signal without the power spectrum, call each draw_* function individually.

fig = plt.figure(constrained_layout=True, figsize=(18, 6))
gs = gridspec.GridSpec(2, 3, figure=fig)
ax1 = fig.add_subplot(gs[:, 0])
ax2 = fig.add_subplot(gs[0, 1])
ax3 = fig.add_subplot(gs[1, 1], sharex=ax2)
ax4 = fig.add_subplot(gs[:, 2])

beorn.plotting.draw_dTb_signal(ax1, multi_z_quantities, color='k', lw=3)
beorn.plotting.draw_Temp_signal(ax2, multi_z_quantities, color='k', lw=3)
beorn.plotting.draw_xHII_signal(ax3, multi_z_quantities, color='k', lw=3)
kplot = beorn.plotting.draw_dTb_power_spectrum_of_z(
                    ax4, multi_z_quantities, parameters, k_value=0.1, color='k', lw=3)
ax4.plot([], [], color='k', label=f'k={kplot:.2f} Mpc$^{{-1}}$')
ax4.legend(loc=0)
ax2.axes.get_xaxis().set_visible(False)
plt.show()

k=0.12 Mpc$^{-1}$
../_images/1914c14e239e88413f294552d9969ae0fb648208bf721ba85002db352786eadf.png

8. Lightcone#

A lightcone stitches together slices from different redshifts to mimic an actual observation: as we observe further into the past, the signal evolves. This is what a radio telescope effectively measures along the line of sight. The horizontal axis shows increasing redshift (earlier cosmic times) from right to left.

lightcone = beorn.structs.Lightcone.build(
    parameters,
    multi_z_quantities,
    quantity="Grid_dTb",
)

fig, ax = plt.subplots(figsize=(18, 6))
beorn.plotting.plot_lightcone(lightcone, ax, "", slice_number=30)
fig.show()

Making lightcone between 0.047619 < z < 0.153845

100%|██████████| 559/559 [00:00<00:00, 4261.44it/s]

...done
../_images/5395607e828197fc0998a0303c5804e60fae2b60eb35a12ae714293d56e03c23.png

9. Compare two runs#

To compare different astrophysical models side-by-side, point IGM_DATA_1 and IGM_DATA_2 at two igm_data_* output directories. BEoRN reads the parameters and grid data directly from those folders — no need to re-specify any parameters manually.

Each TemporalCube.read() call recovers the full Parameters object from the embedded HDF5 metadata, so param1 and param2 automatically reflect the settings used to produce each run.

# Edit these two paths to point at your igm_data_* output directories
IGM_DATA_1 = Path("output/igm_data_py21cmfast_N128_D384_L100_seed12345_4660308f_4024dbb5")
IGM_DATA_2 = Path("output/igm_data_py21cmfast_N128_D384_L100_seed12345_4660308f_457d9d8c")

multi_z_quantities_1 = beorn.structs.TemporalCube.read(IGM_DATA_1)
multi_z_quantities_2 = beorn.structs.TemporalCube.read(IGM_DATA_2)

param1 = multi_z_quantities_1.parameters
param2 = multi_z_quantities_2.parameters

INFO:beorn.structs.temporal_cube:Opened 30 snapshots from output/igm_data_py21cmfast_N128_D384_L100_seed12345_4660308f_4024dbb5 (grid data is lazy — slices load on demand).
INFO:beorn.structs.temporal_cube:Opened 30 snapshots from output/igm_data_py21cmfast_N128_D384_L100_seed12345_4660308f_457d9d8c (grid data is lazy — slices load on demand).

label1 = "Model 1"
label2 = "Model 2"

fig = plt.figure(constrained_layout=True, figsize=(18, 6))
gs = gridspec.GridSpec(2, 3, figure=fig)
ax1 = fig.add_subplot(gs[:, 0])
ax2 = fig.add_subplot(gs[0, 1])
ax3 = fig.add_subplot(gs[1, 1], sharex=ax2)
ax4 = fig.add_subplot(gs[:, 2])

for grid, params, label, color, lstyle in [
    (multi_z_quantities_1, param1, label1, 'C0', '-'),
    (multi_z_quantities_2, param2, label2, 'C1', '--'),
]:
    beorn.plotting.draw_dTb_signal(ax1, grid, label=label, color=color, ls=lstyle, lw=3)
    beorn.plotting.draw_Temp_signal(ax2, grid, label=label, color=color, ls=lstyle, lw=3)
    beorn.plotting.draw_xHII_signal(ax3, grid, label=label, color=color, ls=lstyle, lw=3)
    kplot = beorn.plotting.draw_dTb_power_spectrum_of_z(
                    ax4, grid, params, label=None, k_value=0.1, color=color, ls=lstyle, lw=3)

ax1.legend()
ax4.plot([], [], color='k', label=f'k={kplot:.2f} Mpc$^{{-1}}$')
ax4.legend(loc=0)
ax2.axes.get_xaxis().set_visible(False)
plt.show()

k=0.12 Mpc$^{-1}$

k=0.12 Mpc$^{-1}$
../_images/4d7622137ba485627d421eb9df483c6f1bd8851506fb122be9009dc3b64b3809.png