21cmfast · steven-murray · Mar 10, 2026 · Jul 4, 2025 · Jul 8, 2025 · Jul 10, 2025
diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml
@@ -9,6 +9,7 @@ repos:
   hooks:
       - id: trailing-whitespace
       - id: check-added-large-files
+        args: ['--maxkb=2000']
       - id: check-ast
       - id: check-json
       - id: check-merge-conflict

diff --git a/docs/acknowledge.rst b/docs/acknowledge.rst
@@ -15,11 +15,24 @@ If you use ``21cmFAST v3+`` in your research please cite both of:
 In addition, the following papers introduce various features into ``21cmFAST``. If you use
 these features, please cite the relevant papers.
 
+Lyman alpha multiple scattering
+
+    Flitter, J., Munoz, J. B., Mesinger, A.,
+    "Semi-analytical approach to Lyman-alpha multiple-scattering in 21-cm signal simulations",
+    eprint arXiv:2601.14360, 2026. https://doi.org/10.48550/arXiv.2601.14360.
+
 Discrete Halo Sampler / version 4:
 
     Davies, J. E., Mesinger, A., Murray, S. G.,
     "Efficient simulation of discrete galaxy populations and associated radiation fields during the first billion years",
-    eprint arXiv:2504.17254, 2025. https://doi.org/10.48550/arXiv.2504.17254
+    Astronomy and Astrophysics, vol. 701, A. 236, 2025. https://doi.org/10.1051/0004-6361/202554951.
+
+Photon non-conservation correction:
+
+    Park, J., Greig, B., Mesinger, A.,
+    "Calibrating excursion set reionization models to approximately conserve ionizing photons",
+    Monthly Notices of the Royal Astronomical Society, vol. 517, no. 1,
+    pp 192-200, 2022 https://doi.org/10.1093/mnras/stac2756.
 
 Mini-halos:
 
@@ -42,6 +55,20 @@ Mass-dependent ionizing efficiency:
     functions and the 21-cm signal”, Monthly Notices of the Royal Astronomical Society,
     vol. 484, no. 1, pp. 933–949, 2019. https://doi.org/10.1093/mnras/stz032.
 
+Lightcone and redshift space distortions:
+
+    Greig, B., Mesinger, A.,
+    "21CMMC with a 3D light-cone: the impact of the co-evolution approximation on the astrophysics of reionisation and cosmic dawn.",
+    Monthly Notices of the Royal Astronomical Society, vol. 477, no. 3,
+    pp 3217-3229, 2018. https://doi.org/10.1093/mnras/sty796.
+
+Inhomogeneous recombination:
+
+    Sobacchi, E., Mesinger, A.,
+    "Inhomogeneous recombinations during cosmic reionization",
+    Monthly Notices of the Royal Astronomical Society, vol. 440, no. 2,
+    pp 1662-1673, 2014. https://doi.org/10.1093/mnras/stu377.
+
 If you are unsure which modules are used within your simulations, we provide a handy function
 to print out which works to refer ``py21cmfast.utils.show_references``, which accepts a single instance of
 the ``InputParameters`` class and shows which papers are relevant for your simulation.
diff --git a/docs/conf.py b/docs/conf.py
@@ -84,7 +84,7 @@ def check_module(self):
     'undoc-members', 
     'show-inheritance', 
     'show-module-summary', 
-    'special-members', 
+    #'special-members', 
     'imported-members',
     'inherited-members',
 ]
@@ -94,7 +94,8 @@ def check_module(self):
 autoapi_member_order = 'groupwise'
 autoapi_own_page_level = 'class'
 autoapi_keep_files = True
-autodoc_typehints = 'description'
+autodoc_typehints = 'both'
+autoapi_template_dir = 'templates/autoapi'
 
 autosummary_generate = False
 numpydoc_show_class_members = False
@@ -159,11 +160,3 @@ def check_module(self):
     "templates",
     "**.ipynb_checkpoints",
 ]
-
-def setup(app):
-    """Setup function for Sphinx."""
-    print("SETTING UP THE CLASS BASED DECORATOR DOCUMENTER")
-    app.add_autodocumenter(ClassDecoratedDocumenter)
-    return {
-        'parallel_read_safe': True
-    }
diff --git a/docs/images/papers/multiple_scattering/coevals_z11.png b/docs/images/papers/multiple_scattering/coevals_z11.png
diff --git a/docs/images/papers/multiple_scattering/coevals_z20.png b/docs/images/papers/multiple_scattering/coevals_z20.png
diff --git a/docs/index.rst b/docs/index.rst
@@ -8,13 +8,16 @@
    quickstart
    tutorials
    performance
+   models
    faqs/index
    updates_from_v3
    acknowledge
 
 .. toctree::
    :caption: API Reference
    :hidden:
+   :maxdepth: 1
+
    autoapi/py21cmfast/index
 
 .. toctree::

diff --git a/docs/models.rst b/docs/models.rst
@@ -0,0 +1,84 @@
+Description of Physical Models
+==============================
+
+``21cmFAST`` supports many physical models and effects that can often be toggled on or
+off interchangeably (though some depend on others). The
+:class:`~py21cmfast.MatterOptions` and
+:class:`~py21cmfast.AstroOptions` classes contain all of the flags
+that control which models to include in the simulation.
+
+
+Below we provide a brief explanation on how some of the flags modify the output of
+``21cmFAST``.
+
+Models of the Matter Field
+--------------------------
+
+The parameters and flags of :class:`~py21cmfast.MatterOptions`
+are used to control how cosmological matter fields (e.g. densities, velocities and halo
+properties) are evaluated in the simulation.
+
+Source Model
+~~~~~~~~~~~~
+
+.. note:: Set with the ``SOURCE_MODEL`` parameter of :class:`~py21cmfast.MatterOptions`.
+
+To be filled.
+
+Astrophysical Models
+--------------------
+
+The parameters and flags of :class:`~py21cmfast.AstroOptions`
+are used to control how astrophysical quantities (e.g. star formation rate, UV and
+ionizing radiation) are evaluated in the simulation.
+
+It is important to stress that the generation of
+:class:`~py21cmfast.PerturbedField` and
+:class:`~py21cmfast.HaloCatalog` objects do not depend on these
+parameters (nor on :class:`~py21cmfast.AstroParams`). Therefore, if
+the cache contains :class:`~py21cmfast.PerturbedField` and
+:class:`~py21cmfast.HaloCatalog` objects that had been previously
+generated with a different set of :class:`~py21cmfast.AstroOptions`, these objects will
+be loaded from the cache, instead of being re-evaluated.
+This architecture allows one to quickly simulate different astrophysical models, given a
+cosmological model.
+
+Spin Temperature Fluctuations
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+To be filled.
+
+Minihalos
+~~~~~~~~~
+
+To be filled.
+
+Multiple Scattering of Lyman Alpha Photons
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+.. note:: Toggled with the ``LYA_MULTIPLE_SCATTERING`` flag of
+   :class:`~py21cmfast.AstroOptions` (default: ``False``).
+
+The physical effect that enables the absorption feature in the 21-cm signal during
+cosmic dawn is the strong coupling between the spin temperature and the gas kinetic
+temperature. This coupling is obtained through Lyman alpha radiation that comes from the
+first stars and is absorbed by the IGM. Since the cross section for the interaction
+between photons near the Lyman alpha resonance frequency and HI atoms in the IGM has a
+non-negligible width, the Lyman alpha photons that are absorbed by the IGM had not
+traveled in straight lines, but rather had scattered along their path. This means that
+the effective Lyman alpha emissivity that the IGM "sees" becomes more local, thereby
+increasing the contrast in ``J_alpha`` (Lyman alpha flux) maps in the simulation, as can
+be seen below. In the context of the 21-cm signal, this effect becomes negligible at
+sufficiently low redshifts (normally below :math:`z<15`), once the spin temperature
+completely follows gas kinetic temperature. For more information on this effect in the
+simulation, see `Flitter, Munoz and Mesinger 2026 <https://arxiv.org/pdf/2601.14360>`_.
+
+.. image:: ./images/papers/multiple_scattering/coevals_z20.png
+    :width: 800px
+    :align: center
+    :alt: the effect of the Lyman alpha multiple scattering effect in the simulation, at z=20. Figure is taken from Flitter, Munoz and Mesinger 2026 (arxiv: 2601.14360).
+
+.. image:: ./images/papers/multiple_scattering/coevals_z11.png
+    :width: 800px
+    :align: center
+    :alt: the effect of the Lyman alpha multiple scattering effect in the simulation, at z=20. Figure is taken from Flitter, Munoz and Mesinger 2026 (arxiv: 2601.14360).
diff --git a/docs/updates_from_v3.rst b/docs/updates_from_v3.rst
@@ -1,6 +1,6 @@
-==============================
-Changes from 21cmFAST v3 to v4
-==============================
+=====================
+Changes from v3 to v4
+=====================
 
 This page outlines the main differences between the usage of ``21cmFAST`` versions 3 and
 4. Use this documentation as a reference for understanding and adapting to the updates in

diff --git a/setup.py b/setup.py
@@ -24,7 +24,7 @@ def _read(name: str):
 test_req = [
     "clang-format",
     "clang-tidy",
-    "hmf",
+    "mpmath",
     "pre-commit",
     "pytest>=5.0",
     "pytest-cov",

diff --git a/src/py21cmfast/__init__.py b/src/py21cmfast/__init__.py
@@ -12,7 +12,7 @@
         __version__ = "unknown"
 
 __all__ = [
-    "DATA_PATH",
+    "_DATA_PATH",
     "AngularLightconer",
     "AstroOptions",
     "AstroParams",
@@ -73,7 +73,7 @@
 
 from . import lightconers, plotting, wrapper
 from ._cfg import config
-from ._data import DATA_PATH
+from ._data import _DATA_PATH
 from ._logging import configure_logging
 from ._templates import create_params_from_template, list_templates, write_template
 from .drivers.coeval import Coeval, generate_coeval, run_coeval

diff --git a/src/py21cmfast/_cfg.py b/src/py21cmfast/_cfg.py
@@ -7,7 +7,7 @@
 from typing import ClassVar
 
 from . import yaml
-from ._data import DATA_PATH
+from ._data import _DATA_PATH
 from .c_21cmfast import ffi, lib
 from .wrapper.structs import StructInstanceWrapper
 
@@ -23,7 +23,7 @@ class Config(dict):
     _defaults: ClassVar = {
         "direc": "~/21cmFAST-cache",
         "ignore_R_BUBBLE_MAX_error": False,
-        "external_table_path": DATA_PATH,
+        "external_table_path": _DATA_PATH,
         "HALO_CATALOG_MEM_FACTOR": 1.2,
         "EXTRA_HALOBOX_FIELDS": False,
         "safe_read": True,
@@ -106,3 +106,4 @@ def load(cls, file_name: str | Path):
 
 # On import, load the default config
 config = Config()
+"""A singleton Config object used to manage global configuration options."""
diff --git a/src/py21cmfast/_data/__init__.py b/src/py21cmfast/_data/__init__.py
@@ -1,3 +1,3 @@
 from pathlib import Path
 
-DATA_PATH = Path(__file__).parent
+_DATA_PATH = Path(__file__).parent
diff --git a/src/py21cmfast/drivers/coeval.py b/src/py21cmfast/drivers/coeval.py
@@ -756,6 +756,7 @@ def _redshift_loop_generator(
                     this_xraysource = sf.compute_xray_source_field(
                         redshift=z,
                         hboxes=[*hbox_arr, this_halobox],
+                        previous_ionize_box=getattr(prev_coeval, "ionized_box", None),
                         write=write.xray_source_box,
                         **kw,
                     )

diff --git a/src/py21cmfast/drivers/single_field.py b/src/py21cmfast/drivers/single_field.py
@@ -9,6 +9,7 @@
 import warnings
 
 import numpy as np
+from astropy import constants
 from astropy import units as un
 from astropy.cosmology import z_at_value
 
@@ -462,6 +463,7 @@ def compute_xray_source_field(
     initial_conditions: InitialConditions,
     hboxes: list[HaloBox],
     redshift: float,
+    previous_ionize_box: IonizedBox | None = None,
 ) -> XraySourceBox:
     r"""
     Compute filtered grid of SFR for use in spin temperature calculation.
@@ -478,6 +480,9 @@ def compute_xray_source_field(
         The initial conditions of the run. The user and cosmo params
     hboxes: Sequence of :class:`~HaloBox` instances
         This contains the list of Halobox instances which are used to create this source field
+    previous_ionize_box: :class:`IonizedBox` or None
+        An ionized box at higher redshift. This is only used if `LYA_MULTIPLE_SCATTERING` is true.
+
 
     Returns
     -------
@@ -528,6 +533,34 @@ def compute_xray_source_field(
     # inner and outer redshifts (following the C code)
     zpp_avg = zpp_edges - np.diff(np.insert(zpp_edges, 0, redshift)) / 2
 
+    # Compute the comoving diffusion scale in the case of Lyman alpha multiple scattering
+    if inputs.astro_options.LYA_MULTIPLE_SCATTERING:
+        # TODO: In principle, the diffusion scale varies locally but for simplicty, we consider the global ionization value.
+        # See https://github.com/21cmfast/21cmFAST/issues/606.
+        if previous_ionize_box is None:
+            x_HI = 1.0
+        else:
+            x_HI = previous_ionize_box.neutral_fraction.value.mean()
+        A_alpha = 6.25e8 * un.Hz
+        nu_Lya = 2.46606727e15 * un.Hz
+        n_H_z0 = (
+            (1.0 - inputs.cosmo_params.Y_He)
+            * inputs.cosmo_params.cosmo.critical_density(0)
+            * inputs.cosmo_params.OMb
+            / constants.m_p
+        )
+        # Eq. (24) in arxiv: 2601.14360
+        R_star = 3.0 * constants.c**4 * A_alpha**2 * n_H_z0 * x_HI * (1.0 + redshift)
+        R_star /= (
+            32.0
+            * np.pi**3
+            * nu_Lya**4
+            * inputs.cosmo_params.cosmo.H0**2
+            * inputs.cosmo_params.OMm
+        )
+    else:
+        R_star = 0.0 * un.Mpc
+
     interp_fields = ["halo_sfr", "halo_xray"]
     if inputs.astro_options.USE_MINI_HALOS:
         interp_fields += ["halo_sfr_mini", "log10_Mcrit_MCG_ave"]
@@ -544,6 +577,9 @@ def compute_xray_source_field(
             if inputs.astro_options.USE_MINI_HALOS:
                 box.filtered_sfr_mini.value[i] = 0
                 box.mean_log10_Mcrit_LW.value[i] = inputs.astro_params.M_TURN  # minimum
+                if inputs.astro_options.LYA_MULTIPLE_SCATTERING:
+                    box.filtered_sfr_lw.value[i] = 0
+                    box.filtered_sfr_mini_lw.value[i] = 0
             logger.debug(f"ignoring Radius {i} which is above Z_HEAT_MAX")
             continue
 
@@ -563,6 +599,9 @@ def compute_xray_source_field(
             if inputs.astro_options.USE_MINI_HALOS:
                 box.filtered_sfr_mini.value[i] = 0
                 box.mean_log10_Mcrit_LW.value[i] = hbox_interp.log10_Mcrit_MCG_ave
+                if inputs.astro_options.LYA_MULTIPLE_SCATTERING:
+                    box.filtered_sfr_lw.value[i] = 0
+                    box.filtered_sfr_mini_lw.value[i] = 0
             logger.debug(f"ignoring Radius {i} due to no stars")
             continue
 
@@ -571,6 +610,7 @@ def compute_xray_source_field(
             R_inner=R_inner,
             R_outer=R_outer,
             R_ct=i,
+            R_star=R_star.to("Mpc").value,
             allow_already_computed=True,
         )
 

diff --git a/src/py21cmfast/io/h5.py b/src/py21cmfast/io/h5.py
@@ -205,7 +205,7 @@ def write_outputs_to_group(
         new = array.written_to_disk(H5Backend(group.file.filename, f"{group.name}/{k}"))
         setattr(output, k, new)
 
-    for k in output.struct.primitive_fields:
+    for k in output._struct.primitive_fields:
         try:
             group.attrs[k] = getattr(output, k)
         except TypeError as e:

diff --git a/src/py21cmfast/plotting.py b/src/py21cmfast/plotting.py
@@ -170,7 +170,7 @@ def coeval_sliceplot(
     """
     if kind is None:
         if isinstance(struct, outputs.OutputStruct):
-            kind = struct.struct.fieldnames[0]
+            kind = struct._struct.fieldnames[0]
         elif isinstance(struct, Coeval):
             kind = "brightness_temp"
 

diff --git a/src/py21cmfast/src/HaloCatalog.c b/src/py21cmfast/src/HaloCatalog.c
@@ -207,7 +207,7 @@ int ComputeHaloCatalog(float redshift_desc, float redshift, InitialConditions *b
 
             // now filter the box on scale R
             // 0 = top hat in real space, 1 = top hat in k space
-            filter_box(density_field, box_dim, matter_options_global->HALO_FILTER, R, 0.);
+            filter_box(density_field, box_dim, matter_options_global->HALO_FILTER, R, 0., 0.);
 
             // do the FFT to get delta_m box
             dft_c2r_cube(matter_options_global->USE_FFTW_WISDOM, grid_dim, z_dim,
@@ -369,7 +369,7 @@ int ComputeHaloCatalog(float redshift_desc, float redshift, InitialConditions *b
                 filter_box(density_field, box_dim, 0,
                            physconst.l_factor * simulation_options_global->BOX_LEN /
                                (simulation_options_global->HII_DIM + 0.0),
-                           0.);
+                           0., 0.);
             }
             dft_c2r_cube(matter_options_global->USE_FFTW_WISDOM, simulation_options_global->DIM,
                          D_PARA, simulation_options_global->N_THREADS, density_field);