# Exploring the effect of different cosmologies on the Epoch of Reionization 21-cm signal with POLAR

Anshuman Acharya, Qing-Bo Ma, Sambit K Giri, Benedetta Ciardi, Raghunath Ghara, Garrelt Mellema, Saleem Zaroubi, Ian Hothi, Ilian T Iliev, Léon V.E Koopmans, et al.

<figure><img src="image_1.png" /></figure>

## To cite this version:

Anshuman Acharya, Qing-Bo Ma, Sambit K Giri, Benedetta Ciardi, Raghunath Ghara, et al., Exploring the effect of different cosmologies on the Epoch of Reionization 21-cm signal with POLAR. Monthly Notices of the Royal Astronomical Society, 2025, 1058, pp.1078. 10.1093/mnras/staf1412. hal-04761590

HAL Id: hal-04761590

https://hal.science/hal-04761590v1

Submitted on 31 Oct 2025

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>
# Exploring the effect of different cosmologies on the Epoch of Reionization 21-cm signal with POLAR

Anshuman Acharya<sup>1,★</sup> Qing-bo Ma,<sup>2,3</sup> Sambit K. Giri,<sup>4</sup> Benedetta Ciardi,<sup>1</sup> Raghunath Ghara,<sup>5</sup> Garrelt Mellema,<sup>6</sup> Saleem Zaroubi,<sup>7,8</sup> Ian Hothi,<sup>9,10</sup> Ilian T. Iliev,<sup>11</sup> Léon V. E. Koopmans,<sup>12</sup> and Michele Bianco<sup>12</sup>

Accepted 2025 August 21. Received 2025 August 21; in original form 2024 October 15

## ABSTRACT

A detection of the 21-cm signal power spectrum from the Epoch of Reionization is imminent, thanks to consistent advancements from telescopes such as LOFAR, MWA, and HERA, along with the development of SKA. In light of this progress, it is crucial to expand the parameter space of simulations used to infer astrophysical properties from this signal. In this work, we explore the role of cosmological parameters such as the Hubble constant $H_0$ and the matter clustering amplitude $\sigma_8$, whose values as provided by measurements at different redshifts are in tension. We run N-body simulations using GADGET-4, and post-process them with the reionization simulation code POLAR, that uses L-GALAXIES to include galaxy formation and evolution properties and GRIZZLY to execute 1D radiative transfer of ionizing photons in the intergalactic medium (IGM). We compare our results with the latest James Webb Space Telescope (JWST) observations and explore which astrophysical properties for different cosmologies are necessary to match the observed UV luminosity functions at redshifts $z = 10$ and 9. Additionally, we explore the impact of these parameters on the observed 21-cm signal power spectrum upper limits, focusing on the redshifts within the range of LOFAR 21-cm signal observations ($z \approx 8.5–10$). Despite differences in cosmological and astrophysical parameters, our models cannot be ruled out by the current upper limits. This suggests the need for broader physical parameter spaces for inference modeling to account for all models that agree with observations. However, we also propose stronger constraining power by using a combination of galactic and IGM observables.

**Key words:** galaxies: formation – dark ages, reionization, first stars – cosmology: theory.

## 1 INTRODUCTION

The Epoch of Reionization (EoR) refers to the period of the Universe when sufficient astrophysical objects formed and emitted UV photons to (re)ionize the neutral hydrogen of the intergalactic medium (IGM). One of the most important probes for studying this epoch is the brightness temperature fluctuations of the 21-cm line of neutral hydrogen (H1) as observed in emission or absorption against the cosmic microwave background radiation (CMB; Field 1959; Hogan & Rees 1979; Madau, Meiksin & Rees 1997; Shaver et al. 1999; Tozzi et al. 2000; Ciardi & Madau 2003; Zaroubi 2013).

The 21-cm signal is produced by the forbidden hyperfine spin-flip transition of the ground state of neutral hydrogen. As the excited state has an extremely long lifetime and neutral hydrogen is abundant in the IGM at the start of the EoR, this signal can be used as a tracer of the IGM ionization state. A statistical detection of the strength of the fluctuations of its brightness temperature allows us to constrain models of the Universe during this early formation stage.

For this purpose, multiple interferometric low-frequency radio telescopes have been designed, such as PAPER<sup>1</sup>, MWA<sup>2</sup>, HERA<sup>3</sup>,

© The Author(s) 2025.

Published by Oxford University Press on behalf of Royal Astronomical Society. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0/), which permits unrestricted reuse, distribution, and reproduction in any medium, provided the original work is properly cited.

<footnote n="1"><sup>1</sup>Max-Planck-Institut für Astrophysik, Garching D-85748, Germany</footnote>

<footnote n="2"><sup>2</sup>School of Physics and Electronic Science, Guizhou Normal University, Guiyang 550001, PR China</footnote>

<footnote n="3"><sup>3</sup>Guizhou Provincial Key Laboratory of Radio Astronomy and Data Processing, Guizhou Normal University, Guiyang 550001, PR China</footnote>

<footnote n="4"><sup>4</sup>Nordita, KTH Royal Institute of Technology and Stockholm University, Hannes Alfvén's väg 12, SE-106 91 Stockholm, Sweden</footnote>

<footnote n="5"><sup>5</sup>Department of Physical Sciences, Indian Institute of Science Education and Research Kolkata, Mohanpur, WB 741 246, India</footnote>

<footnote n="6"><sup>6</sup>The Oskar Klein Centre, Department of Astronomy, Stockholm University, AlbaNova, SE-10691 Stockholm, Sweden</footnote>

<footnote n="7"><sup>7</sup>Kapteyn Astronomical Institute, University of Groningen, P.O. Box 800, NL-9700AV Groningen, the Netherlands</footnote>

<footnote n="8"><sup>8</sup>Astrophysics Research Centre of the Open University of Israel, Ra'anana 4353701, Israel</footnote>

<footnote n="9"><sup>9</sup>LERMA, Observatoire de Paris, PSL Research University, CNRS, Sorbonne Université, F-75014 Paris, France</footnote>

<footnote n="10"><sup>10</sup>Laboratoire de Physique de l'ENS, ENS, Université PSL, CNRS, Sorbonne Université, Universitée Paris Cité, F-75005 Paris, France</footnote>

<footnote n="11"><sup>11</sup>Astronomy Centre, Department of Physics & Astronomy, Pevensey III Building, University of Sussex, Falmer, Brighton, BN1 9QH, UK</footnote>

<footnote n="12"><sup>12</sup>Institute for Particle Physics and Astrophysics, ETH Zurich, Wolfgang-Pauli-Str 27, CH-8093 Zurich, Switzerland</footnote>

<footnote>*E-mail: anshuman@mpa-garching.mpg.de</footnote>

<footnote n="1"><sup>1</sup>Precision Array to Probe EoR, http://eor.berkeley.edu</footnote>

<footnote n="2"><sup>2</sup>Murchison Widefield Array, http://www.mwatelescope.org</footnote>

<footnote n="3"><sup>3</sup>Hydrogen Epoch of Reionization Array, https://reionization.org/</footnote>
LOFAR<fnref n="4" />, and the upcoming SKA<fnref n="5" /> While the signal is yet to be detected, upper limits of the 21-cm power spectrum have been obtained and are becoming increasingly tighter (e.g. Mertens et al. 2020; Trott et al. 2020; HERA Collaboration 2023; Acharya et al. 2024c). These have allowed us to rule out some extreme astrophysical models through comparisons with simulations (e.g. Ghara et al. 2020; Mondal et al. 2020; Greig et al. 2021a, b; Abdurashidova et al. 2022). New methods that utilize multiredshift power spectrum observations (Choudhury et al. 2024; Ghara et al. 2024), wavelet statistics (Hothi et al. 2024), and simulation-based inference (Saxena et al. 2023; Greig et al. 2024) have also been developed to improve the understanding of parameters governing the properties of the IGM.

At the same time, galaxies that formed during the EoR that are responsible for H1 ionization have also been studied using hydrodynamical and/or radiative transfer simulations with different prescriptions, such as CROC (Gnedin 2014; Gnedin & Kaurov 2014; Esmerian & Gnedin 2021, 2022, 2024), FIRE (Ma et al. 2018), TECHNICOLOR DAWN (Finlator et al. 2018), SPHINX (Rosdahl et al. 2018), CRASH simulations (Eide et al. 2018, 2020; Ma et al. 2021; Kostyuk et al. 2023; Basu, Garaldi & Ciardi 2024), C²RAY simulations (Mellema et al. 2006; Dixon et al. 2016; Giri & Mellema 2021; Giri et al. 2024; Hirling et al. 2024), CODA I,II, III (Ocvirk et al. 2016, 2020; Lewis et al. 2022), ASTRID (Bird et al. 2022), COLDSIM (Maio, Péroux & Ciardi 2022; Maio & Viel 2023; Casavecchia et al. 2024), THESAN (Garaldi et al. 2022; Kannan et al. 2022; Smith et al. 2022; Garaldi et al. 2024), and SPICE (Bhagwat et al. 2024). Although these simulations are ideal tools to investigate the IGM and galactic properties, they are computationally expensive and thus limited in box size and/or in the number of simulations that can be run. Thus, they cannot be used for building a wide range of models for comparison with the 21-cm power spectrum or their upper limits.

To this aim, semi-numerical approaches such as SIMFAST21 (Santos et al. 2010a, b), BEARS (Thomas et al. 2009), 21CMFAST (Mesinger & Furlanetto 2007; Mesinger, Furlanetto & Cen 2011), GRIZZLY (Ghara, Choudhury & Datta 2015; Ghara et al. 2018), REIONYUGA (Mondal, Bharadwaj & Majumdar 2017), ARTIST (Molaro et al. 2019), AMBER (Trac et al. 2022), and BEORN (Schaeffer, Giri & Schneider 2023) are typically employed. To incorporate a more realistic modelling of galactic properties into these methods, codes such as ASTRAEUS (Hutter et al. 2021, 2024), MERAXES (Mutch et al. 2016; Balu et al. 2023), and POLAR (Ma et al. 2023, hereafter M23) have been developed.

The LOFAR EoR Key Science Project (KSP) team has designed POLAR, a semi-numeric approach that strikes a balance between speed and complexity of relevant physical processes, by post-processing N-body Dark Matter (DM) simulations with Semi-Analytic Models (SAMs) of galaxy formation and evolution, and subsequently applying the 1D radiative transfer (RT) code GRIZZLY (Ghara et al. 2015, 2018), an updated version of BEARS (Thomas et al. 2009). This allows a fast modelling of the 21-cm signal power spectrum, while at the same time accounting for important galactic properties and RT effects. The 21-cm signal can be used to derive physical properties of the Universe (McQuinn et al. 2006), but as it primarily comes from the IGM, it provides only a partial picture. Thus, instead, we combine observations of the 21-cm line with data obtained in other frequency bands (e.g. from JWST) to jointly constrain IGM and galactic properties. This has been explored in the past, for example, by using cross-power spectra of galaxy clustering observations (Lidz

et al. 2009; Vrbanec et al. 2016; Berti, Spinelli & Viel 2024), or intensity mapping (Mao et al. 2008; La Plante et al. 2023; Berti et al. 2024) at different redshifts. Eventually, the goal of such studies would be to jointly constrain the physical properties of the Universe from multiple observables of galaxies and the IGM.

In this work, we investigate the impact of varying cosmological parameters on the 21-cm signal at $z \approx 10.11, 9.16$, and 8.3, i.e. the three redshifts targeted by LOFAR (Mertens et al. 2025). For this, we develop models with the most extreme choices of cosmological parameters observed so far. Thus, we focus on those parameters whose low and high redshift measurements show maximum tension (here we define 'tension' as having more than $2\sigma$ difference in the measurements of the parameters), i.e. the Hubble constant $H_0 = 100h \text{ km s}^{-1} \text{ Mpc}^{-1}$ (where $h$ is a dimensionless parametrization of the Hubble constant), and the matter clustering amplitude $\sigma_8$<fnref n="6" />. In future work, though, we plan to run models with various combinations of observed cosmological parameters to build a suite of N-body simulations that span a large parameter space. Specifically, we will build a range of models whose corresponding astrophysical parameters are derived by fitting multiwavelength observables with an MCMC-based algorithm, which in turn can be used for parameter forecasting (see M23 and Ma et al. 2025 for more details).

Our fiducial model is the one of Planck Collaboration VI (2020), with $h = 0.6766$ and $\sigma_8 = 0.8102$, while a higher $h$ value is adopted to explore the low redshift measurement of $h = 0.733$ from studies of Cepheid variables in the host galaxies of 42 Type Ia supernovae (Riess et al. 2022). Similarly, a lower $\sigma_8$ value of 0.702 is adopted from the low-redshift measurement of anisotropic galaxy clustering measurement analysis (Tröster et al. 2020). Lastly, we also adopt a $\sigma_8 = 0.88$ from recent eROSITA results (Ghirardini et al. 2024) to consider both extremes of $\sigma_8$ measurements. Motivated by results obtained in the context of the 21-cm signal shown by Giri et al. (2023) and Acharya et al. (2024b), we additionally implement the Fixed & Paired (F&P) approach (Angulo & Pontzen 2016) to boost effective volumes of the simulations and suppress cosmic variance.

The paper is structured as follows. In Section 2, we discuss the setup of the POLAR simulations with different cosmologies. In Section 3, we present the resulting galactic and IGM properties, and in Section 4, we discuss the implications of varying astrophysical and cosmological parameters for inference modeling from 21-cm signal observations, jointly constrained with other observables of the EoR. Finally, we summarize our results in Section 5.

## 2 METHODOLOGY

Building on M23, we utilize a similar setup by running N-body DM simulations, and then post-processing them with the L-GALAXIES sam (Henriques et al. 2015, 2020; Barrera et al. 2023) to model the formation and evolution of galaxies, and with the 1D radiative transfer code GRIZZLY (Ghara et al. 2015; Ghara et al. 2018) to model the gas ionization and the 21-cm signal from neutral hydrogen. In

<footnote n="6">⁶While the $S_8$ parameter is also occasionally used, it is a derived quantity given by $S_8 = \sigma_8 \sqrt{\Omega_m}/0.3$. Thus, when a tension is present, it is driven by the underlying tensions of $\sigma_8$ and/or $\Omega_m$. For example, the $3.4\sigma$ tension between the Tröster et al. (2020) and Planck Collaboration VI (2020) measurements is completely driven by the $3.7\sigma$ tension between the $\sigma_8$ parameters. We note that in some cases $S_8$ is measured directly (and not derived from $\sigma_8$ and $\Omega_m$), and a tension up to $\approx 2\sigma$ can be found (for example Amon et al. 2023; Dark Energy Survey and Kilo-Degree Survey Collaboration, 2023; Longley et al. 2023). Most simulations, though, adopt $\sigma_8$ and $\Omega_m$ as independent parameters, and we follow this typical convention.</footnote>

<footnote n="4">⁴Low-Frequency Array, http://www.lofар.org</footnote>

<footnote n="5">⁵Square Kilometre Array, https://www.skao.int/en</footnote>
Sections 2.1 and 2.2, we highlight the key parameters used for setting up the simulations and the analysis done in this work. Further, in Section 2.3, we propose two cases: one where we keep the astrophysical parameters the same for all four cosmological models, and one in which we tune them to match UV luminosity functions (UVLFs) observed with JWST and HST at z = 10 and 9. We refer to the first and second case as ‘unconstrained’ and ‘constrained’, respectively.

## 2.1 Dark-matter simulations

We use the GADGET-4 code (Springel et al. 2021) for running N-body DM simulations with a volume of $(150 h^{-1} \text{ cMpc})^3$ and $2048^3$ DM particles. This translates into a DM particle mass of $5 \times 10^7 \text{ M}_\odot$, matching the particle resolution of larger simulations used by the LOFAR EoR KSP team (see for example Giri et al. 2019a, b). The box size is chosen to probe wave-modes of $k \le 1.0h \text{ cMpc}^{-1}$, where the best results from observations with LOFAR, HERA, MWA (and eventually SKA) are expected (Koopmans et al. 2015). This condition requires simulations with box sizes $> 100 h^{-1} \text{ cMpc}$ (Iliev et al. 2014), although Kaur, Gillet & Mesinger (2020) suggests box sizes $> 175 h^{-1} \text{ cMpc}$ being necessary when accounting for X-ray heating of the IGM. To sidestep this issue, we assume that at our redshifts of interest, the IGM has already been heated above the CMB temperature (see Section 2.2.2). We additionally employ the F&P approach (Angulo & Pontzen 2016) to mitigate sample variance, and thus run two realizations of each DM simulation. Each pair of realizations has the mode amplitudes fixed to the square root of the initial matter power spectrum, and the phases of the second realizations (B-series) are obtained by mirroring those of the first realizations (A-series). This is in contrast to the traditional method of generating initial conditions, where the mode amplitudes are sampled from a Gaussian distribution centred on the expectation value, which is given by the initial matter power spectrum. Taking the F&P approach suppresses the variance-induced fluctuations (which in traditional initial conditions mainly affect the large-scale modes, where the sampling of the power spectrum amplitude is scarce and therefore more susceptible to deviations from its expectation value) up to the fourth perturbative order (Angulo & Pontzen 2016). This means that while small-scale properties, such as individual halos, may not show any noticeable difference, averaged large-scale statistics of an F&P pair reproduce very closely the ensemble averages of traditional simulations. The F&P approach of averaging observables of these two realizations has been shown to boost the statistical precision of the matter power spectrum, bispectrum and halo mass function (Angulo & Pontzen 2016; Villaescusa-Navarro et al. 2018; Chartier et al. 2021; Maion, Angulo & Zennaro 2022), Lyman-$\alpha$ power spectra (Anderson et al. 2019), 21-cm signal power spectrum and bispectrum (Giri et al. 2023; Acharya et al. 2024b) and several other quantities derived from simulations (Villaescusa-Navarro et al. 2018; Klypin, Prada & Byun 2020). Beyond this, we use the same random seed to minimize the effect of randomized initial conditions for both the A and B-series of simulations.

For our reference simulation, we assume a ‘fiducial’ $\Lambda$-Cold Dark Matter ($\Lambda$CDM) cosmological model based on Planck Collaboration VI (2020; specifically the TT, TE, EE + lowE + lensing + BAO case), setting the dark energy density parameter as $\Omega_\Lambda = 0.6889$, the total matter density parameter, including both baryons and dark matter as $\Omega_m = \Omega_b + \Omega_{dm} = 0.3111$, the baryonic matter density parameter as $\Omega_b = 0.04897$, the Hubble constant as $H_0 = 100 h \text{ km s}^{-1} \text{ Mpc}^{-1}$ with $h = 0.6766$, the amplitude of matter density fluctuations on scales of $8 h^{-1} \text{ Mpc}$ as $\sigma_8 = 0.8102$, and

<table><caption>**Table 1.** The four N-body dark matter simulation models considered in this work. From left to right, the model name, the adopted value of $h$ and $\sigma_8$, and the reference for the values. All other cosmological parameters are fixed to the Planck Collaboration (2020) values.</caption><thead><tr><th>Model</th><th><i>h</i></th><th>σ<sub>8</sub></th><th>Reference</th></tr></thead><tbody><tr><td>fiducial</td><td>0.6766</td><td>0.8102</td><td>Planck Collaboration VI (2020)</td></tr><tr><td><i>h</i> high</td><td>0.7330</td><td>0.8102</td><td>Riess et al. (2022)</td></tr><tr><td>σ<sub>8</sub> low</td><td>0.6766</td><td>0.7020</td><td>Tröster et al. (2020)</td></tr><tr><td>σ<sub>8</sub> high</td><td>0.6766</td><td>0.8800</td><td>Ghirardini et al. (2024)</td></tr></tbody></table>

the tilt of the primordial power spectrum of scalar perturbations as $n_s = 0.9665$. Further, we consider three additional cosmologies, where we vary $h$ and $\sigma_8$ while keeping the other parameters fixed to the above-mentioned values. This is necessary to consider the maximum impact of changing cosmological parameters across the range of observed values. First, we adopt $h = 0.7330$ (from Riess et al. 2022), and refer to this as the ‘$h$ high’ model. In addition, we explore two extreme values of $\sigma_8$, namely the ‘$\sigma_8$ low’ case with $\sigma_8 = 0.702$ (Tröster et al. 2020), and the ‘$\sigma_8$ high’ case with $\sigma_8 = 0.88$ (Ghirardini et al. 2024). The details of the four simulation models are listed in Table 1. As the $\sigma_8$ parameter represents the amplitude of matter density fluctuations, a higher (lower) $\sigma_8$ value leads to larger (smaller) dark matter haloes, which in turn can lead to more (less) numerous and more (less) massive galaxies as compared to the fiducial case. If all astrophysical parameters are kept the same, this would result in a faster (slower) reionization. Thus, if we want, e.g. to produce a similar overall star formation history, the astrophysical parameters which affect the amount of stars formed would have to be reduced (boosted) in their impact. A higher (lower) Hubble constant leads to faster (slower) expansion of space-time, i.e. to larger (smaller) voids. This translates into more (less) matter being pushed together, which leads to the formation of larger (smaller) dark matter halos. However, as the Hubble constant does not directly control the matter density, halo sizes are less sensitive to changes in $h$ than in $\sigma_8$. As a consequence, the impact on galaxy formation and evolution of the former is expected to be similar to the latter, albeit smaller. As the interplay between $\sigma_8$ and/or $h$ and the astrophysical parameters governing galaxy formation and evolution is complex, understanding their impact on the reionization process is not straightforward and should be explored carefully.

The initial conditions for all simulations are generated at $z = 199$ with a second-order Lagrangian perturbation theory based on the NGENIC algorithm implemented into GADGET-4, using the same linear theory power spectrum as Hernández-Aguayo et al. (2023). We generate 90 snapshots between $z = 20$ and $z = 5$, for which we save the full snapshot information, i.e. the positions and velocities of all the dark matter particles. However, in order to ensure a maximum step size of 10 Myrs between snapshots used to build halo merger trees with the Friend-Of-Friends (FOF) group finding algorithm (Springel et al. 2001), we perform a finer gridding to generate a total of 156 output time steps for FOF groups. This is necessary because the dynamical time of a typical galactic disc at high redshift is of the order of tens of Myrs (Poole et al. 2016; Balu et al. 2023), which also is the lifetime of massive stars (Barkana & Loeb 2001; Heger & Woosley 2002; Schaerer 2002; Bromm & Larson 2004). Thus, a longer time-step would miss how reionization progresses when post-processed with SAMs like L-GALAXIES and radiative transfer codes like GRIZZLY. We use a standard linking length of 0.2 times the mean particle spacing and a minimum group size of 64 DM particles (corresponding to a minimum halo size of
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Halo mass functions (HMFs) for our four models: fiducial (cyan), *h* high (red), *σ*<sub>8</sub> low (purple), and *σ*<sub>8</sub> high (green) at *z* = 10 (top panel) and 9 (bottom). The A- and B-series are shown with dashed and dotted lines, respectively. The corresponding fit for each simulation using the Tinker et al. (2008) model is shown with solid lines of the corresponding colours.</figcaption></figure>

≈ 3 × 10⁹ M☉). Substructures within halos are identified using the SUBFIND-HBT algorithm (Han et al. 2018; Springel et al. 2021). Finally, the gravitational softening length is set to 0.025 of the mean particle spacing, i.e. ≈ 1.83 h⁻¹ ckpc.

To evaluate the performance of the simulations, in Fig. 1, we present the halo mass function (HMF) as $M_h \frac{d\ln h}{dM_h}$, where $M_h$ and $n_h$ are the halo mass and number density respectively, and compare them with fits by Tinker et al. (2008), computed using the PYTHON package HMF (Murray 2014) at $z=9$ and 10. As in this study, Tinker et al. (2008) take into account the non-universality of the HMF by considering the impact of varying the cosmological parameters in the ΛCDM model. We note that the fiducial and *h* high models at these redshifts produce very similar HMFs, despite some noticeable differences in their corresponding fits. A deeper exploration indicates that our simulation set-up is less sensitive to differences in the Hubble parameter, given the choice of the starting redshift, and that the initial conditions were generated using second-order Lagrangian perturbation theory. A significantly higher starting redshift or third/fourth-order Lagrangian perturbation theory may thus be required to better match the fits. Nevertheless, some differences do exist, especially at the extremities of the HMF, which lead to the differences noted in the subsequent sections.

Further, we also find a mild under-prediction of the HMF at the low mass end and an over-prediction at the high mass end. This may be due to the choice of the linking length of 0.2, which at significantly high redshifts may lead to overlinking. A lower value of 0.17, as suggested by Watson et al. (2013), may resolve this issue. However, we still use 0.2 for minimizing differences between the set-ups used

for past applications of L-GALAXIES (e.g. Henriques et al. 2015, 2020; Barrera et al. 2023, and M23) and our set-up. Lastly, we note that the B-series has slightly more massive halos in all models. This is purely because of the choice of the random seed; the same seed that leads to voids in the A-series produces regions of massive halo formation in the B-series due to the F&P method described above.

## 2.2 Galactic and IGM properties with POLAR

POLAR is a semi-numerical model designed to obtain high-z galaxy properties and the 21-cm signal from the IGM in a fast and robust manner. It combines the semi-analytical galaxy formation and evolution model of L-GALAXIES (Henriques et al. 2015, 2020; Barrera et al. 2023) with the one-dimensional radiative transfer code GRIZZLY (Ghara et al. 2015, 2018).

### 2.2.1 Semi-Analytic modeling of galaxies

We use L-GALAXIES as described by Barrera et al. (2023) as a post-processing module of GADGET-4. This is an updated version of the publicly available L-GALAXIES 2020 that was used by M23. L-GALAXIES implements most major physical processes of gas cooling, star formation, galaxy mergers, supernovae feedback, black hole growth, AGN feedback, and dust attenuation. While Barrera et al. (2023) largely builds on Henriques et al. (2015), we also consider the parameters used by Henriques et al. (2020) before adapting them for our purposes. In particular, we focus on those that scale the star formation efficiency ($\alpha_{SF}$), star formation efficiency during galaxy mergers ($\alpha_{SF, burst}$), AGN accretion rate ($k_{AGN}$), reheating of cold gas by star formation ($\epsilon_{reheat}$, $V_{reheat}$), and the energy released by each supernova ($E_{SN}$). L-GALAXIES assumes that stars form from the cold gas within the disc of each galaxy. Thus, these parameters control the amount of cold gas available in the disc and the efficiency of converting it into stars. For more details, refer to the supplementary material of Henriques et al. (2015).

We tune the parameters to match photometric observations of the UVLFs from JWST and HST at $z=10$ and 9, as these redshifts are observationally relevant for the LOFAR EoR KSP. We note that L-GALAXIES assumes a 100 per cent escape fraction for the UV photons, and thus we also implement the dust attenuation approach of Henriques et al. (2015). However, this is a simplified model and more complex ones may lead to greater suppression of the UV luminosity function. In Table 2, we list all the possible parameters available in L-GALAXIES, and their values set to match UVLF observations as discussed in Section 2.3.

### 2.2.2 Radiative transfer and the 21-cm signal

Modelling the EoR requires the inclusion of radiative transfer to describe the hydrogen ionization and heating. For this, we take the results of the *N*-body simulations from Section 2.1 and the semi-analytic modeling of galaxies from Section 2.2.1, and post-process them with the 1D radiative transfer code GRIZZLY, as done by M23. GRIZZLY uses pre-computed ionization and temperature profiles of gas for different source and density properties at various redshifts to model the ionization and heating processes and the differential brightness temperature of the 21-cm signal ($\delta T_b$).

More specifically, GRIZZLY requires as input the gridded matter density field and dark matter halo masses from the *N*-body simulation, as well as the corresponding galactic stellar masses and stellar ages obtained with L-GALAXIES. Our reference simulation has a grid
<table><caption>**Table 2.** List of astrophysical parameters in L-GALAXIES (see Henriques et al. 2015, 2020, for more details) that can be tuned to control galaxy formation and evolution. We show the names of the parameters in column 1, and Henriques et al. (2015, 2020) values tuned to $z < 3$ observations in columns 2 and 3. Our cosmological models (fiducial, $h$ high, $\sigma_8$ low, and $\sigma_8$ high) tuned to UVLF observations at $z = 10$ and 9 are shown in columns 4, 5, 6, and 7, respectively. We highlight the parameters that differ from Henriques et al. (2015, 2020) values in bold and fix the values of the unchanged parameters to the latter. For the unconstrained case (see the text for more details), we use the fiducial model parameters for all our N-body simulations.</caption><thead><tr><th>Parameter</th><th>Henriques et al. (2015)</th><th>Henriques et al. (2020)</th><th>fiducial</th><th><i>h</i> high</th><th>σ<sub>8</sub> low</th><th>σ<sub>8</sub> high</th></tr></thead><tbody><tr><td>α<sub>SF</sub></td><td>0.025</td><td>0.060</td><td><b>0.20</b></td><td><b>0.20</b></td><td><b>0.50</b></td><td><b>0.20</b></td></tr><tr><td>α<sub>SF,burst</sub></td><td>0.60</td><td>0.50</td><td><b>0.80</b></td><td><b>0.80</b></td><td><b>0.90</b></td><td><b>0.80</b></td></tr><tr><td>β<sub>SF,burst</sub></td><td>1.90</td><td>0.38</td><td>0.38</td><td>0.38</td><td>0.38</td><td>0.38</td></tr><tr><td>k<sub>AGN</sub> [10<sup>-3</sup>M<sub>☉</sub> yr<sup>-1</sup>]</td><td>5.3</td><td>2.5</td><td><b>2.0</b></td><td><b>2.0</b></td><td><b>2.0</b></td><td><b>2.0</b></td></tr><tr><td>f<sub>BH</sub></td><td>0.041</td><td>0.066</td><td>0.066</td><td>0.066</td><td>0.066</td><td>0.066</td></tr><tr><td>V<sub>BH</sub> [km s<sup>-1</sup>]</td><td>750</td><td>700</td><td>700</td><td>700</td><td>700</td><td>700</td></tr><tr><td>M<sub>r.p.</sub> [10<sup>14</sup>M<sub>☉</sub>]</td><td>1.2</td><td>5.1</td><td>5.1</td><td>5.1</td><td>5.1</td><td>5.1</td></tr><tr><td>α<sub>dyn.fric.</sub></td><td>2.5</td><td>1.8</td><td>1.8</td><td>1.8</td><td>1.8</td><td>1.8</td></tr><tr><td>ε<sub>reheat</sub></td><td>2.6</td><td>5.6</td><td><b>8.0</b></td><td><b>8.0</b></td><td><b>8.0</b></td><td><b>8.0</b></td></tr><tr><td>V<sub>reheat</sub> [km s<sup>-1</sup>]</td><td>480</td><td>110</td><td><b>250</b></td><td><b>250</b></td><td><b>250</b></td><td><b>250</b></td></tr><tr><td>β<sub>reheat</sub></td><td>0.72</td><td>2.90</td><td>2.90</td><td>2.90</td><td>2.90</td><td>2.90</td></tr><tr><td>η<sub>eject</sub></td><td>0.62</td><td>5.50</td><td>5.50</td><td>5.50</td><td>5.50</td><td>5.50</td></tr><tr><td>V<sub>eject</sub> [km s<sup>-1</sup>]</td><td>100</td><td>220</td><td>220</td><td>220</td><td>220</td><td>220</td></tr><tr><td>β<sub>eject</sub></td><td>0.8</td><td>2.0</td><td>2.0</td><td>2.0</td><td>2.0</td><td>2.0</td></tr><tr><td>γ<sub>reinc</sub> [10<sup>10</sup> yr<sup>-1</sup>]</td><td>3.0</td><td>1.2</td><td>1.2</td><td>1.2</td><td>1.2</td><td>1.2</td></tr><tr><td>E<sub>SN</sub> [10<sup>51</sup> erg]</td><td>1.0</td><td>1.0</td><td><b>0.80</b></td><td>1.0</td><td><b>0.15</b></td><td><b>2.00</b></td></tr><tr><td>Z<sub>yield</sub></td><td>0.046</td><td>0.030</td><td>0.030</td><td>0.030</td><td>0.030</td><td>0.030</td></tr></tbody></table>

of 256³ cells, resulting in a cell size of ≈ 600 h⁻¹ ckpc⁷. Next, we assume that the gas density scales by a factor of Ω<sub>b</sub>/Ω<sub>dm</sub> with the dark matter density. Lastly, as done in M23, we use the Binary Population and Spectral Synthesis (BPASS; Stanway & Eldridge 2018) code to model the Spectral Energy Distributions (SEDs) of stellar sources. In future work, we will also explore the impact of other source types such as X-ray binaries, shock-heated interstellar medium, and accreting black holes (as done in Eide et al. 2018, 2020; Ma et al. 2021).

GRIZZLY computes δT<sub>b</sub> as follows (see Furlanetto, Oh & Pierpaoli 2006):

$$ \delta T_b = 27x_{HI}(1 + \delta_B) \left(1 - \frac{T_{CMB}}{T_S}\right) \times \left[ \left(\frac{\Omega_b h^2}{0.023}\right) \left(\frac{0.15}{\Omega_m h^2} + \frac{1+z}{10}\right)^{1/2} \right] \text{mK}, \quad (1) $$

where $x_{HI}$ is the fraction of neutral hydrogen, $\delta_B$ is the fractional overdensity of baryons, $T_S$ is the hydrogen spin temperature, $T_{CMB}$ is the temperature of the CMB photons at redshift $z$, $h$ is the dimensionless Hubble constant, and $\Omega_b$ and $\Omega_m$ are the baryonic matter density and total matter density parameters respectively. As done by M23, we assume $T_S \gg T_{CMB}$ which is valid when the IGM has been sufficiently heated by X-ray sources and expected to be the case in the range of redshift of interest here. Additionally, we also ignore the impact of redshift space distortions. Thus equation (1) is reduced to

$$ \delta T_b = 27x_{HI}(1 + \delta_B) \left[ \left( \frac{\Omega_b h^2}{0.023} \right) \left( \frac{0.15}{\Omega_m h^2} + \frac{1+z}{10} \right)^{1/2} \right] \text{mK}. \quad (2) $$

⁷The grid size is chosen after checking the convergence of the volume-averaged neutral hydrogen fraction and the 21-cm power spectrum in the fiducial model for a range of grid sizes (from 64³ to 512³ cells in steps of 2³) across redshifts 12 > z > 5. We note that the values are within 5 per cent, and thus opt for 256³ to ensure we do not miss out on the impact of low-density regions in any redshift regime.

We define the power spectrum of δT<sub>b</sub> as (Mertens et al. 2020; Metha & Berger 2024):

$$ P_{21\text{cm}}(\mathbf{k}) = \delta_D(\mathbf{k} + \mathbf{k}') \langle \delta T_b(\mathbf{k}) \delta T_b(\mathbf{k}') \rangle, \quad (3) $$

where $\delta_D$ is the Dirac delta function, δT<sub>b</sub>(k) is the differential brightness temperature in Fourier space, and (...) is the ensemble average. In the following, we will report our results in terms of the normalized form of the power spectrum, given by (see Peacock 1999, for details):

$$ \Delta_{21\text{cm}}^2(k) = \frac{k^3}{2\pi^2} \times P_{21\text{cm}}(k). \quad (4) $$

## 2.3 UV luminosity functions

Because Henriques et al. (2015, 2020) tuned the L-GALAXIES parameters to the low-redshift Universe ($z < 3$), it is necessary to adapt them to match the redshift regime that we are interested in (see Vani et al. 2024, on discrepancies from observations at higher redshifts). In particular, for the LOFAR EoR KSP team, focussing on $10 > z > 8.5$ is crucial, as this is the regime in which the LOFAR telescope is most sensitive. Thus, we change the values of the parameters listed in Section 2.2.1 to match the observations of the UVLFs from HST legacy fields and JWST programmes (Bouwens et al. 2015, 2021, 2023a, b; Finkelstein et al. 2015; Harikane et al. 2022, 2023; Leung et al. 2023; Adams et al. 2024; McLeod et al. 2024) at $z = 10$ and 9 by employing the least square minimisation⁸. We iteratively vary parameters one at a time, with the range going from 1/100th of the lower value to 100 times the higher value between the fitted values obtained by Henriques et al. (2015, 2020) for each parameter. While typically the average energy released by each supernova is expected to be ≈ 10⁵¹ ergs (see Jones et al. 1998), we find that a mild reduction

⁸In the least square minimization the square of the residuals between UVLF observations and models is summed for both redshifts over a range of different values for the parameters. The values that lead to the minimum sum of the squares of the residuals are accepted.
to $0.8 \times 10^{51}$ ergs provides a better fit, while still being consistent with observations (Rubin et al. 2016; Kozyreva et al. 2025). Further, the resulting values for $\alpha_{SF}$ and $\alpha_{SF,burst}$ suggest higher star formation efficiency, in agreement with past observations (Stefanon et al. 2021), abundance matching estimates (Tacchella et al. 2018) and models (Sun & Furlanetto 2016; Sipple & Lidz 2024). The mild reduction of $k_{AGN}$ with respect to the values in Henriques et al. (2020) leads to lower gas accretion by AGNs, allowing for more star formation and thus more galaxies at the brightest magnitudes ($M_{1600,AB} < -20$), which is necessary to match recent observations. Similarly, the small differences in the supernova feedback-based heating efficiency parameters $\epsilon_{reheat}$ and $V_{reheat}$ are required to improve the agreement for $M_{1600,AB} > -19$. Indeed, as supernovae can effectively disrupt smaller galaxies, reducing the impact of their heating efficiency allows for more cold gas to form, and thus more star formation, further increasing the UVLF at the low brightness end and providing a closer match with observations.

We have also varied other parameters affecting star formation ($\beta_{SF,burst}$), the heating of cold gas by stars and supernova feedback (such as $\beta_{reheat}$, $\eta_{eject}$, $V_{eject}$, $\beta_{eject}$, $\gamma_{reinc}$), black hole properties ($f_{BH}$, $V_{BH}$), and galaxy mergers ($M_{r,p}$, $\alpha_{dyn.fric}$), but their impact on the UVLFs is minimal (see M23, for a detailed analysis of the impact of different parameters). For such parameters, we fix their value to that obtained by Henriques et al. (2020). In future work, we will broaden the range of astrophysical parameters and analyse their effect on other observables. Additionally, we also plan to employ more rigorous MCMC-based techniques for fitting parameters (as shown in Ma et al. 2025).

So far, we only considered the fiducial model. While the values of the parameters obtained may be sensitive to the method used to derive them (in this case, the least squares minimization), the qualitative difference in their values for the different cosmological models introduced in Section 2.1 can still be explored. Thus, we consider two cases. In the first one, the same L-GALAXIES parameters of the fiducial model above are adopted for all N-body simulations to investigate the impact of different cosmologies on the 21-cm signal independently from observations of galaxies at high redshifts. We refer to this model as ‘unconstrained’. This case enables us to gain a clearer understanding of the impact of changing cosmologies on the 21-cm signal, without astrophysical processes potentially cancelling out their effects on observables. However, not changing the astrophysical parameters for the different cosmological models may produce UVLFs that are inconsistent with observations. Thus, we also consider another scenario, called the ‘constrained’ case. Here, the L-GALAXIES parameters (specifically $\alpha_{SF}$, $\alpha_{SF,burst}$, and $E_{SN}$) are changed for each cosmological model to match the observed UVLFs. As the average energy released per supernova is the most sensitive parameter in scaling the UVLF while maintaining its overall shape, we first focus on tuning it. For the h high model, we find that fixing it at $10^{51}$ ergs and using the fiducial model’s values for $\alpha_{SF}$ and $\alpha_{SF,burst}$ is enough to get the best fit via the least squares minimization described above. The increase in $E_{SN}$ as compared to the fiducial model is required to reduce the amount of cold gas in the galaxy available for star formation, balancing the increased clustering of matter due to the higher Hubble constant. Similarly, due to even stronger matter clustering due to a larger $\sigma_8$ in the $\sigma_8$ high model, we note that a value of $2 \times 10^{51}$ ergs for $E_{SN}$ provides the best fit. Because of this, we expect that just reducing $E_{SN}$ for the $\sigma_8$ low model to improve star formation efficiency is enough despite the lower matter clustering. However, we note that in addition to reducing the average supernova energy to $0.15 \times 10^{51}$ ergs, it is also necessary to increase the parameters regulating the star formation

<figure><img src="image_4.png" /><figcaption>**Figure 2.** A and B-series averaged UVLFs for the unconstrained (top) and constrained (bottom) cases at $z = 10$ and 9 for the fiducial (cyan), h high (red), $\sigma_8$ low (purple), and $\sigma_8$ high (green). The averaging of the F&P pairs allows to compensate for stochasticity in underfilled bins of the UVLFs, by suppressing the impact of cosmic variance. We also show the dust attenuated (solid) and unattenuated (dashed) UVLFs, and JWST and HST observations (Bouwens et al. 2015, 2021, 2023a, b; Finkelstein et al. 2015; Harikane et al. 2022, 2023; Leung et al. 2023; Adams et al. 2024; McLeod et al. 2024, grey circles).</figcaption></figure>

efficiency, $\alpha_{SF}$, and the burstiness of star formation efficiency during mergers, $\alpha_{SF,burst}$. Further raising the values beyond 0.50 for $\alpha_{SF}$ and 0.90 for $\alpha_{SF,burst}$ does not induce any appreciable difference in the sum of the squares of the residuals.

These differences in the astrophysical parameters for various cosmologies are potentially more interesting, as they allow for a joint constraint of astrophysical and cosmological parameters based on multifrequency observations. The list of all the L-GALAXIES parameters adopted for the four models is provided in Table 2.

We show the UVLFs for the unconstrained (top) and constrained (bottom) case in Fig. 2 for all cosmological models, along with observations at $z = 10$ and 9. In the unconstrained case, we note that the fiducial and h high models have similar UVLFs, while the $\sigma_8$ low model underpredicts the UVLF as compared to observations. On the other hand, the UVLF in the $\sigma_8$ high model is mildly overpredicted.
These results are in agreement with the impact that the parameters are expected to have on matter clustering, as discussed above. In the constrained case, the boost in star formation due to higher $\alpha_{SF}$ and $\alpha_{SF,burst}$ values and a lower $E_{SN}$ value, allows the $\sigma_8$ low model to match the bright end of the UVLF. However, due to lower matter clustering, it does not produce as many faint galaxies as the other cases, and sees a steep drop at $M_{1600,AB} > -18$. On the other hand, the high $E_{SN}$ in the $\sigma_8$ high model strongly suppresses star formation and thus matches other cosmological models, as well as JWST and HST observations. Additionally, it produces more faint galaxies even at $M_{1600,AB} > -16$. Between the $h$ high and fiducial models, the difference is minimal at $z = 10$ due to the mildly higher $E_{SN}$ in the $h$ high model. However, at $z = 9$ the impact of higher matter clustering in the $h$ high model shows up in the form of more galaxies at the low luminosity end.

Lastly, we note that in the fiducial model the brightest galaxy at $z = 14$ has $M_{1600,AB} = -20.36$, while at $z = 13$ it has $M_{1600,AB} = -20.60$, which is in agreement with the recent spectroscopically confirmed galaxies at $z = 14.32$ and 13.90 (Carniani et al. 2024). For completeness, in Appendix A, we also show the UVLFs in the range $12 > z > 5$ for the constrained case. We find that our models reasonably agree with UVLF observations across this redshift range.

## 2.4 Reionization history

Before discussing the reionization history in more detail, we note that the escape fraction of UV photons had been set to $f_{esc} = 12.5$ per cent in order for the fiducial model to reionize completely (i.e. the average neutral fraction $\langle x_{HI} \rangle$ becomes zero) by $z = 5$. For the other models, though, full reionization is reached at different redshifts. Further, we note that this is a simplified assumption, as $f_{esc}$ is expected to vary with galaxy and halo properties (as shown by Ma et al. 2020; Kostyuk et al. 2023) and maybe with redshift (while e.g. Ma et al. 2020 and Kostyuk, Ciardi & Ferrara 2025 find redshift variation, Gorce et al. 2018 favours $f_{esc}$ to be constant with redshift), and thus can be strongly affected by modelling choices and approximations. Nevertheless, our choice aligns with recent estimates from observations (Mascia et al. 2023; Bosman & Davies 2024; Jaskot et al. 2024, which report median values of 10–12 per cent). In the future, we plan to consider $f_{esc}$ as a free parameter as well, and, as a reference, in Appendix B, we show the reionization history of the fiducial model with $f_{esc} = 25$ per cent. This is of particular interest, as observations like Jaskot et al. (2024) have reported higher escape fractions ($f_{esc} > 20$ per cent) for $\ge 33$ per cent of observed galaxies.

For a qualitative view of the reionization history, in Fig. 3, we show the A-series (i.e. with ‘fixed’ initial conditions in the F&P pair of simulations) middle slices of the $\delta T_{\text{h}}$ cubes for the four cosmological models in the constrained case at $z = 12, 10, 8$ and 6. For completeness, in Appendix C1, we also show the corresponding B-series slices, and in Appendix C2 both the A and B-series middle slices for the unconstrained case. We note that there are no major differences until $z = 8$, when more and larger ionized regions appear in the $\sigma_8$ low model, while fewer and smaller ones are present in the $\sigma_8$ high model. This leads to the $\sigma_8$ low model being almost completely ionized by $z = 6$ (note that while in Fig. 3, the shown slice at $z = 6$ happens to be fully ionized, there are still some neutral regions in the full simulation box. The same is shown below by computing the volume-averaged neutral hydrogen fraction), while the $\sigma_8$ high model still has large neutral regions. This is because in these cases the choice of the parameters controlling the astrophysical processes at high redshift is balanced by the impact on the matter clustering

due to the choice of $\sigma_8$, but at lower redshift, once enough haloes have formed, the astrophysical processes dominate. In particular, a higher (lower) average energy per supernova in the $\sigma_8$ high (low) case significantly impedes (bolsters) star formation as too much cold gas is blown away (collapsed). However, the differences between the astrophysical parameters of the $h$ high and the fiducial model are smaller, so the differences in the size of the ionized regions from $z = 10$ are mainly driven by matter clustering. Thus, overall, the choice of both astrophysical and cosmological parameters leads to significantly different speeds at which reionization happens across all four models. Indeed, despite the relatively similar maps observed at $z = 8$, by $z = 6$ the differences are much more evident.

In Fig. 4, we present the redshift evolution of the average of the A and B-series volume-averaged neutral hydrogen fraction, $\langle x_{HI} \rangle$, in the unconstrained (top) and constrained (bottom) cases, together with observational constraints. As in some models, the redshift of reionization is not properly captured due to the coarse redshift resolution<fnref n="9" />, we additionally provide a fit to these curves (shown as a dotted line) to estimate more accurately the end of reionization. We note that in the unconstrained case, all models reionize at different redshifts. Reionization is the fastest (at $z = 6.4$) in the $\sigma_8$ high model, as structure formation happens earlier. The $\sigma_8$ low model is still only 45 per cent reionized by $z = 5$ because of the lack of sources, and shows the poorest agreement with observational constraints. In the constrained case, the differences are reduced, but still significant. In particular, we note that the $\sigma_8$ low model reionizes the fastest at $z = 5.75$, while the $\sigma_8$ high model is only $\approx 65$ per cent ionized by $z = 5$. This reversal of reionization histories with respect to the unconstrained case is due to the impact of different astrophysical parameters. Specifically, the significantly higher $E_{SN}$ in the $\sigma_8$ high model blows away gas and thus suppresses star formation, which consequently reduces the production of ionizing photons. On the other hand, the higher $\alpha_{SF}$ and $\alpha_{SF,burst}$ in the $\sigma_8$ low model along with the lower $E_{SN}$ significantly boosts star formation.

As the four models reionize at different redshifts and the end of the EoR can be observationally constrained (e.g. detections of Gunn–Peterson troughs in the Lyman $\alpha$ forest as shown in Becker et al. 2015, Qin et al. 2021, and Bosman et al. 2022), this in turn limits the possible choice of astrophysical and cosmological parameters. We note that to model more realistically the final phases of reionization, the unresolved Lyman limit systems which govern absorption during the final stages of the reionization process should be accounted for (see Georgiev, Mellema & Giri 2024; Giri et al. 2024). These, though, are not included here.

Finally, the four models produce a Thomson scattering optical depth in the range (0.041–0.067) and (0.048–0.060) for the unconstrained and constrained case respectively, assuming an Helium (He) II fraction equal to the H II fraction, and an instantaneous He ii reionization at $z = 3$. These numbers are in agreement with Planck observations (Planck Collaboration VI 2020; de Belsunce et al. 2021; Giarè, Di Valentino & Melchiorri 2024).

## 3 RESULTS

In this section, we present the results from our simulations with respect to various galactic and IGM properties for the unconstrained and constrained case. To reduce the effect of cosmic variance, we

<footnote n="9">⁹While we have a large number of snapshots as discussed in Section 2.1 for the sake of simplicity, to run POLAR we adopt a uniform step size of $\Delta z = 0.5$.</footnote>
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Maps of $\delta T_b$ of the A-series middle slices of single cell thickness (i.e. $\approx 600 \ h^{-1} \ \text{ckpc}$) for the four cosmological models in the constrained case at $z = 12, 10, 8$, and 6 (from top to bottom). Here, the dark areas represent the ionized regions with $\delta T_b = 0$. Note that $\delta T_b$ cannot have negative values due to the assumption of $T_S \gg T_{\text{CMB}}$.</figcaption></figure>

always consider the A and B-series averages when comparing to
observations.

## 3.1 Unconstrained case

In the unconstrained case, we discuss the galactic and IGM observables for models with the UVLFs reported in the top panel of Fig. 2.

### 3.1.1 Galactic properties

In Fig. 5, we present the mass-binned star formation rate (SFR) of galaxies in the four models at $z = 10$ and 9, with the top panel showing the unconstrained case. The solid lines are the median SFR, with the shaded regions referring to the 16th to 84th percentiles. The star formation main sequence (SFMS) in each case is additionally compared with results from various recent JWST programmes repre-

sented by grey circles with a binning of $\pm 0.25$ around the redshift. As all parameters affecting the star formation rate were kept the same across all four models, we note that the overall trends are similar, and agree with observations. The impact of matter clustering is seen at the low-mass end, where a larger clustering allows for the formation of a statistically significant sample of smaller mass galaxies in the $\sigma_8$ high and $h$ high models, while the smallest mass galaxies formed in the $\sigma_8$ low case are an order of magnitude more massive. Further, at the high mass end, we note that the $\sigma_8$ high model produces the most massive galaxies, which also have the highest SFR.

To confirm the agreement in star formation rates, in Fig. 6 we additionally look at the stellar fraction, which serves as a proxy for the global star formation efficiency (SFE) at the same redshifts with the top panel showing the unconstrained case. We define the stellar fraction as the ratio of the stellar and halo mass scaled by $f_b = \Omega_b / \Omega_m$, i.e. $M_{\text{star}} / M_{\text{halo}} f_b$, as a function of the halo mass $M_{\text{halo}}$. We also compare this to *Spitzer* observations at $z = 7$ (Stefanon
<figure><img src="image_6.png" /></figure>

**Figure 4.** Redshift evolution of the average of the A and B-series volume-averaged neutral hydrogen fraction (x<sub>HI</sub>) for the fiducial (cyan), *h* high (red), σ<sub>8</sub> low (purple) and σ<sub>8</sub> high (green) models for the unconstrained (top) and constrained (bottom) cases. Dotted lines refer to a fit to the curves, which is used for a better estimate of the redshift of reionization when the redshift resolution is too coarse. The vertical grey dashed lines indicate the redshifts observationally relevant for LOFAR (z = 10.11, 9.16 and 8.3), and the black solid line is drawn at (x<sub>HI</sub>) = 0.5 to guide the eye. Grey circles are a collection of observational constraints (Fan et al. 2006; Totani et al. 2006; Ota et al. 2008; Ouchi et al. 2010; Bolton et al. 2011; Dijkstra, Mesinger & Wyithe 2011; McGreer, Mesinger & Fan 2011; Mortlock et al. 2011; Ono et al. 2012; Chornock et al. 2013; Jensen et al. 2013; Robertson et al. 2013; Schroeder, Mesinger & Haiman 2013; Pentericci et al. 2014; Schenker et al. 2014; Choudhury et al. 2015; McGreer, Mesinger & D'Odorico 2015; Mesinger et al. 2015; Sobacchi & Mesinger 2015; Greig et al. 2017; Davies et al. 2018; Mason et al. 2018; Greig, Mesinger & Bañados 2019; Hoag et al. 2019; Jones et al. 2024, collected in the CORECON module; Garaldi 2023).

et al. 2021) and to abundance matching estimates from Tacchella
et al. (2018). As mentioned above, the parameters controlling the
star formation efficiency are the same across all models, and thus the
corresponding stellar fraction is the same as well. The only difference
is in the largest galaxies formed, which is directly proportional to
σ<sub>8</sub> as it controls the clustering of mass. We note that the stellar
fraction is marginally higher than the abundance matching estimates
for the lower mass halos, but the difference is ≲ 0.3 dex. Further, the
weak mass dependence of their trend is in agreement with zoom-in
hydrodynamical simulations like FIREBOX<sup>HR</sup> (Feldmann et al. 2024
and COLDSIM (Maio & Viel 2023).

3.1.2 21-cm signal from the IGM

We note that despite the disagreement in UVLFs for the four
cosmological models shown in the top panel of Fig. 2, the models
still reproduce fairly well observable quantities such as the SFR
and the stellar fraction, as shown in Section 3.1.1. Further, they
are broadly in agreement with available data. We now look at

<figure><img src="image_7.png" /></figure>

**Figure 5.** A and B-series averaged star formation rate (SFR) as a function of stellar mass for the unconstrained (top) and constrained (bottom) cases, for the fiducial (cyan), *h* high (red), σ<sub>8</sub> low (purple) and σ<sub>8</sub> high (green) models at z = 10 and 9. Solid lines refer to the median star formation rate in each mass bin, with 16th and 84th percentiles shown as shaded regions. Grey circles are observations from various JWST programmes (Arrabal Haro et al. 2023a, b; Atek et al. 2023; Fujimoto et al. 2023; Heintz et al. 2023a, b; Jin et al. 2023; Leethochawalit et al. 2023; Long et al. 2023; Looser et al. 2023; Papovich et al. 2023; Robertson et al. 2023; Treu et al. 2023; Bouwens et al. 2023b; Helton et al. 2024; Jung et al. 2024).
<figure><img src="image_8.png" /><figcaption>**Figure 6.** A and B-series averaged stellar fraction (which serves as a proxy for the star formation efficiency) as a function of halo mass for the unconstrained (top) and constrained (bottom) cases, for the fiducial (cyan), *h* high (red), *σ*<sub>8</sub> low (purple) and *σ*<sub>8</sub> high (green) models at *z* = 10 and 9. Solid lines refer to the median stellar fraction in each mass bin, with 16<sup>th</sup> and 84<sup>th</sup> percentiles shown as shaded regions. Grey circles are *Spitzer* observations from Stefanon et al. (2021), and grey triangles are abundance matching estimates from Tacchella et al. (2018).</figcaption></figure>

IGM observables, and more specifically at the 21-cm signal power spectrum (calculated according to equation 4), shown in the top panel of Fig. 7 for the unconstrained case at *z* = 9. In particular, we focus on the 0.05 < *k* < 0.5*hcMpc*<sup>-1</sup> regime, which is where the LOFAR telescope is the most sensitive. At redshift *z* = 9, we note from the top panel of Fig. 4 that the matter clustering primarily governs the reionization history. Thus, the *σ*<sub>8</sub> low model having the highest ⟨*x*<sub>HI</sub>⟩ suggests that the power amplitude up to *k* = 0.40 *hcMpc*<sup>-1</sup> for it is higher than the other three models, as it is primarily driven by the matter power spectrum. Additionally, the power at these scales is further boosted, because while ⟨*x*<sub>HI</sub>⟩ ≈ 1, there are still voxels where it is ≪ 1 as they are included in ionized regions. These ionized regions in the *σ*<sub>8</sub> low model are small and few enough to lead to higher contrast on average in the differential brightness temperature across all scales greater than the size of the ionized bubbles. In other words, as matter clustering translates into the neutral hydrogen fraction term

<figure><img src="image_9.png" /><figcaption>**Figure 7.** A and B-series averaged normalized 21-cm signal power spectrum (Δ<sup>2</sup><sub>21cm</sub>) for the unconstrained (top) and constrained (bottom) cases, for the fiducial (cyan), *h* high (red), *σ*<sub>8</sub> low (purple) and *σ*<sub>8</sub> high (green) models at *z* = 9.</figcaption></figure>

of equation (2) as well as the overdensity term, in a scenario where the ionized regions are not large and numerous enough to wash out contrast, it leads to a marginal boost in power.

The lower ⟨*x*<sub>HI</sub>⟩ in the other three models suggests suppression of the power spectrum for the 21-cm signal as compared to the matter power spectrum due to more numerous and larger ionized regions, especially in the *σ*<sub>8</sub> high and *h* high models. This is due to the washing out of the contrast driven purely by the large-scale structure. However, the more numerous ionized regions present at smaller scales induce on average a greater contrast, and thus the *h* high model reaches values similar to those of the fiducial model, and the *σ*<sub>8</sub> high model has the highest power overall at *k* > 0.40 *hcMpc*<sup>-1</sup>.

For a better understanding of the 21-cm signal power spectrum evolution, in Fig. 8, we additionally present its redshift dependence at large scales, i.e. *k* = 0.15 *hcMpc*<sup>-1</sup>, with the top panel showing the unconstrained case. This allows us to compare how the 21-cm signal evolves in the four cosmological models across redshifts relevant to LOFAR observations. As suggested by the reionization histories, we note that the *σ*<sub>8</sub> low model has a larger contrast due to slower reionization, and thus has a higher amplitude of Δ<sup>2</sup><sub>21cm</sub> (*k* = 0.15 *hcMpc*<sup>-1</sup>) at the beginning of reionization. However, the *σ*<sub>8</sub> high and *h* high models reach a higher amplitude by *z* ≈ 8 due to the formation of larger ionized regions.

## 3.2 Constrained case

As in this case we have also changed parameters that control star formation efficiency and supernova feedback, it is crucial to
<figure><img src="image_10.png" /><figcaption>**Figure 8.** Redshift evolution of the A and B-series averaged normalized 21-cm signal power spectrum ($\Delta_{21\text{cm}}^2$) at $k = 0.15 \text{ hCpc}^{-1}$ for the unconstrained (top) and constrained (bottom) cases, for the fiducial (cyan), h high (red), $\sigma_8$ low (purple), and $\sigma_8$ high (green) models. The vertical dashed lines indicate the redshifts relevant for LOFAR, i.e. $z = 10.11, 9.16$, and 8.3.</figcaption></figure>

understand their impact on observables other than the UVLFs shown in the bottom panel of Fig. 2.

### 3.2.1 Galactic properties

As in Section 3.1.1, in the bottom panel of Fig. 5 we look at the A and B-series averaged SFMS and compare them with observations from JWST at $z = 10$ and 9 for the constrained case. We note that despite differences in the energy of supernovae ($\alpha_{\text{SF}}$ and $\alpha_{\text{SF, burst}}$ also differ for the $\sigma_8$ low model), the high mass ends of the SFMS of the four models still agree with each other and with observations. This is because the most massive galaxies are also the most luminous, and matching their UVLFs in Section 2.3 consequently matches their star formation properties. However, at the low-mass end, the SFMS of the four models shows some deviation at both redshifts. This is because at lower masses, the impact of astrophysical parameters starts to dominate, and while the $\sigma_8$ high model forms galaxies with masses smaller than the $\sigma_8$ low model, the median SFR of these low mass galaxies is an order of magnitude smaller than the lowest mass galaxies in the $\sigma_8$ low model, because the impact of the energy released per supernova is sufficient to clear out gas and suppress star formation in the low-mass galaxies.

To further explore the impact of the energy released by supernovae, in the bottom panel of Fig. 6, we also look at the global SFE for the constrained case by again considering the stellar fraction. Here, we notice that, despite similar SFRs at the high-mass end, the efficiency of star formation in each case differs. The stellar fraction for the $\sigma_8$

low model is significantly boosted compared to observations, while the fiducial and h high models only mildly overpredict the observed stellar fraction. The $\sigma_8$ high case, on the other hand, is the best match. We note that the stellar fraction is inversely proportional to $E_{\text{SN}}$ for all halo masses. The weak mass dependence trend of the stellar fraction is the same compared to the unconstrained case, and in agreement with the trends of the FIREbox^HR simulations (Feldmann et al. 2024) and COLDSIM (Maio & Viel 2023) for the fiducial, $\sigma_8$ high and h high models, while the $\sigma_8$ low model is an order of magnitude higher but with the same slope.

### 3.2.2 21-cm signal from the IGM

While the UVLFs at $M_{1600,AB} < -18$ in the constrained case are similar for all four cosmological models, the SFR and stellar fraction show differences, and in some cases do not match observations. In such a scenario, it becomes interesting to analyze the impact on the 21-cm signal power spectrum, which is presented in the bottom panel of Fig. 7 at $z = 9$ for the constrained case. Unlike the unconstrained case where the differences in matter clustering contributed to differences not only in the overdensity term but also in the neutral hydrogen fraction term in equation (2), here the changes in astrophysical parameters also contribute to the neutral hydrogen fraction and reduce the impact of matter clustering, leading to a drastic reduction in the differences between the value of $\langle x_{HI} \rangle$ across the four models. Thus, the differences in the 21-cm signal power spectra of the four models are governed primarily by the overdensity term. This leads to a higher power across all wave modes for the $\sigma_8$ high model due to greater matter clustering. While one would expect the h high model to have slightly higher power than the fiducial one due to the dependence on overdensity, the neutral-fraction term in equation (2) reduces the contrast due to a marginally faster reionization process, so that the power in the two cases becomes similar. On the other hand, due to lower matter clustering, the $\sigma_8$ low model has a consistently lower power amplitude as compared to the fiducial model.

In the bottom panel of Fig. 8, we present the redshift evolution of the 21-cm signal power spectrum at $k = 0.15 \text{ hCpc}^{-1}$ for the constrained case and notice a significant reduction in the differences of $\Delta_{21\text{cm}}^2(k = 0.15 \text{ hCpc}^{-1})$ values between the various cosmological models at $12 > z > 5$ compared to the extent of differences in the unconstrained case (top panel of Fig. 8). Indeed, at $10 > z > 8.5$, where LOFAR observations are focused, variations in the $\sigma_8$ parameter lead to differences in the power spectrum of only a few $\text{mK}^2$ (in contrast, even with improvements in signal extraction, with current LOFAR data, Acharya et al. 2024c reported the lowest $2\sigma$ measurement error of $\approx 25^2 \text{ mK}^2$ on the recovered power spectrum). While the exact magnitude of the difference is subject to the position of the peak of $\Delta_{21\text{cm}}^2(k = 0.15 \text{ hCpc}^{-1})$, which in turn is governed by modelling assumptions like the choice of $f_{\text{esc}}$ (see Appendix B for a comparison between $f_{\text{esc}} = 12.5$ per cent and 25 per cent), the key takeaway is that significantly different choices of cosmological and astrophysical parameters can still lead to similar 21-cm signal observables, even when constrained by the UVLFs from JWST and HST observations. Thus, to explore the parameter space of models that agree with the upper limits of 21-cm signal power spectrum observations at the redshifts of interest to LOFAR, one should still consider various cosmological and astrophysical parameters as free, unlike what has been done so far for parameter inference not just for observations from LOFAR, but also from HERA and MWA. However, jointly constraining these models with a more diverse set of galactic observables from JWST
and upcoming surveys with Euclid and SPHEREx can allow us to
narrow down the choice of viable ones. Additionally, tomographic
imaging of the distribution of neutral Hydrogen with the upcoming
SKA-Low could also be used for additional constraints on IGM
properties from images such as those shown in Fig. 3 (see for e.g.
Bianco et al. 2021).

4 DISCUSSION

By analysing various galactic and IGM observables, we note that the interplay of astrophysical and cosmological parameters can lead to similarities as well as differences. For example, in the constrained case, while agreement with observed UVLFs leads to similar SFR for the high-mass galaxies across all models, there are significant differences in the SFR and stellar fraction for low-mass galaxies. On the other hand, while in the unconstrained case, the SFR and stellar fraction show a very similar dependence on stellar mass in all models because the parameters that regulate such quantities are the same, the UVLFs exhibit significant differences.

For the IGM, we note that the 21-cm signal power spectrum at the redshifts of interest for LOFAR is very similar across all four models in the constrained case (bottom panel, Fig. 7). This is due to the choice of astrophysical parameters, which reduces the impact of the cosmological ones, and consequently the differences between the reionization histories in comparison to the unconstrained case. In the unconstrained case, we note significant differences in reionization histories simply because the choice of cosmological parameters changes the overall matter clustering.

These similarities in terms of 21-cm observables have a profound impact on inference modeling, as models with very different cosmological and astrophysical parameters may still produce 21-cm signal power spectra in agreement with current upper limits (for example from Mertens et al. 2020; Trott et al. 2020; HERA Collaboration 2023; Acharya et al. 2024c), or possibly also with an eventual measurement. This is because, while the end of reionization may differ among models, at *z* > 7 the extent and speed of reionization are still largely similar. This means that fixing a priori any of these parameters could exclude viable models. While introducing cosmological and astrophysical parameters as free parameters during inference modelling has been attempted in several earlier studies (see e.g. Santos & Cooray 2006; Liu & Parsons 2016; Kern et al. 2017; Hassan, Andrianomena & Doughty 2020; Schneider, Schaeffer & Giri 2023), these have used only approximate reionization simulation codes due to the significant computational costs involved in building such a high-dimensional parameter space. Indeed, in a more physical approach, it is necessary to run a large number of dark-matter only simulations with different cosmologies, each post-processed with a large set of astrophysical parameters. The problem worsens if one would like to boost the resolution of the simulations to take into account the role of mini halos (see Haiman, Abel & Madau 2001; Iliev, Shapiro & Raga 2005), dwarf galaxies (Wu & Kravtsov 2024), or Lyman Limit Systems (Georgiev et al. 2024; Giri et al. 2024) which likely play a significant role in the EoR. While faster semi-numerical codes and emulators can be employed, they are not ideal tools to model galaxy-scale physics and radiative transfer effects.

Because of this, boosting the resolution of low-resolution simula-
tions using analytic techniques based on smaller boxes with higher
resolution has been proposed (see for example Nasirudin, Iliev &
Ahn 2020; Barsode & Choudhury 2024). Newer methods based on
Machine Learning, such as Generative Adversarial Networks (GANs,
as done in Zhang et al. 2024), bypass the requirement of a higher
resolution simulation. However, it needs to be explored whether

such methods are robust enough to factor in small differences in
individual parameters. In subsequent work with POLAR, we intend to
include resolution-boosting techniques to resolve haloes and galaxies
down by at least two orders of magnitude in mass. The development
and implementation of such techniques will additionally allow us
to create more diverse training sets of power spectra for building
Machine Learning kernels for Gaussian Process Regression based
signal extraction, as shown by Mertens, Bobin & Carucci (2024) and
Acharya et al. (2024a, c).

5 SUMMARY

In this work, we have investigated the impact of different cosmological models on the 21-cm signal and some galactic properties using POLAR, which combines N-body dark matter simulations run with GADGET-4 with the semi-analytic model of galaxy formation L-GALAXIES, and the 1D radiative transfer code GRIZZLY. We have applied the framework to four different cosmological models: a ‘fiducial’ model, which adopts values of the cosmological parameters from Planck Collaboration VI (2020); a ‘h high’ model, with h = 0.733 based on results from studies of Cepheid variables in the host galaxies of 42 Type Ia supernovae (Riess et al. 2022); a ‘σ<sub>8</sub> low’ case, with σ<sub>8</sub> = 0.702 from an anisotropic galaxy clustering measurement analysis done by Tröster et al. (2020); and a ‘σ<sub>8</sub> high’ case, with σ<sub>8</sub> = 0.880 according to recent eROSITA results (Ghirardini et al. 2024). We additionally used the Fixed & Paired approach (Angulo & Pontzen 2016) to suppress cosmic variance by boosting the effective volumes of the simulations. For all the quantities analysed, we took an average of the fixed initial conditions case (which we refer to as the A-series) and its corresponding pair (referred to as the B-series).

We then choose astrophysical parameters in L-GALAXIES such
that the fiducial model matches UV luminosity functions from HST
legacy fields and JWST programmes at z = 10 and 9. In particular, we
focus on parameters affecting the star formation efficiency (α<sub>SF</sub>), star
formation efficiency during galaxy mergers (α<sub>SF,burst</sub>), AGN accretion
rate (k<sub>AGN</sub>), reheating of cold gas by star formation (ε<sub>reheat</sub>, V<sub>reheat</sub>),
and the energy released by each supernova (E<sub>SN</sub>).

To investigate the effect of adopting different values for *h* and σ<sub>8</sub> on other galactic and IGM observables, we use the same astrophysical parameters also for the other three cosmological models. We refer to this as the ‘unconstrained’ case, as by not changing astrophysical parameters for the different cosmological models the resulting UVLFs will not necessarily be consistent with the observed ones. We also built a ‘constrained’ case, where we instead choose the astrophysical parameters such that the UVLFs obtained in each cosmological model are consistent with those observed at *z* = 10 and 9. For this, we increased the star formation efficiency and reduced the energy released per supernova in the σ<sub>8</sub> low model to boost star formation and, in turn, its UVLF. On the other hand, we increased the energy released per supernova for the *h* high and σ<sub>8</sub> high models to suppress star formation and mitigate the impact of increased matter clustering as compared to the fiducial model. For the radiative transfer calculations, we only took into account the stellar sources while modeling the spectral energy distributions and chose a global escape fraction of 12.5 per cent across all redshifts for all four models.

Our results can be summarized as follows:

(i) **Reionization history:** In the unconstrained case, the values of the cosmological parameters strongly influence the matter clustering, which in turn leads to significant differences in how reionization
progresses in each model. While the $\sigma_8$ high model reionizes by $z \approx 6.4$, the $\sigma_8$ low model is only 45 per cent reionized even by $z = 5$. In the constrained case, instead, the impact of astrophysical parameters is significant, and the $\sigma_8$ low model is the first to reionize, at $z = 5.75$, while the $\sigma_8$ high model is only 70 per cent reionized by $z = 5$.

(ii) **Star formation rate (SFR):** In the unconstrained case, no significant differences are observed as the parameters controlling star formation have been kept the same. However, in the constrained case the SFR in low-mass galaxies is impacted by the required energy released per supernova, and thus the SFR is boosted for the $\sigma_8$ low model and suppressed for the $\sigma_8$ high model. In both cases, all models broadly agree with observations at $z = 10$ and 9.

(iii) **Global star formation efficiency (SFE):** As the parameters controlling star formation are the same in the unconstrained case, the global SFE estimated via the stellar fraction at all masses is similar in all models. The stellar fraction is also in agreement with observations and abundance matching estimates. However, in the constrained case, due to a higher value of the parameter regulating the star formation efficiency, the stellar fraction in the $\sigma_8$ low model is also higher. The other three models show minor differences, governed by the differences in the energy released per supernova.

(iv) **21-cm signal power spectrum:** In the unconstrained case, the differences between the four models are not only due to the difference in overdensity but also to those in the neutral fraction and its redshift evolution. However, in the constrained case, the neutral fraction at $z > 8$ is similar in all models, as the impact of matter clustering on the neutral fraction is cancelled out by the impact of setting different values for the astrophysical parameters controlling star formation. Thus, the 21-cm signal power spectrum at higher redshifts is mostly dictated by the overdensity term, and because of this it shows smaller differences among the four models.

Overall, we conclude that jointly varying cosmological and astrophysical parameters can lead to differences in some observables (e.g. low-mass SFR and stellar fraction), while others stay similar (e.g. UVLFs and 21-cm signal power spectrum with the latter remaining well within the expected observational errors despite minor differences). In particular, we note that despite significantly different galactic processes and reionization histories, the 21-cm power spectra are very similar in power across k-bins and in agreement with current observational upper limits (Mertens et al. 2020; Trott et al. 2020; HERA Collaboration 2023; Acharya et al. 2024c). Due to this, when doing inference modelling, it is essential to consider all cosmological and astrophysical parameters as free parameters, with other observational constraints serving as priors. While a limited exploration of astrophysical parameters can be done with excursion set algorithms, a semi-analytic model provides a significantly more rigorous and physical approach to modelling galactic properties. For varying cosmological parameters though, it is necessary to run a large number of N-body simulations to populate the prior parameter space. However, running so many high-resolution simulations would be prohibitively expensive in terms of computational resources and time. To address this issue, in future work we intend to incorporate techniques for boosting the resolution of less costly, low-resolution simulations using either analytic or machine learning techniques. Such an implementation would also additionally allow us to build broader training sets of power spectra for the Machine Learning kernels used with Gaussian Process Regression by

the LOFAR EoR Key Science Project team (as done by Acharya et al. 2024a).

Lastly, to achieve a stronger constraining power on a broadened parameter space, it is essential to jointly employ a variety of galactic observables, e.g. from JWST, Euclid, and upcoming SPHEREx, as well as tomographic imaging of neutral Hydrogen in the IGM with the future SKA-Low.

## ACKNOWLEDGEMENTS

AA thanks the EoR research group at MPA, Volker Springel, and Rüdiger Pakmor for helpful discussions. AA acknowledges Nordita's PhD Fellows Visitor programme for funding a part of this work, and thanks Nordita for their hospitality. Nordita is supported in part by NordForsk. RG acknowledges support from SERB, DST Ramanujan Fellowship no. RJF/2022/000141. GM is supported by the Swedish Research Council project grant 2020-04691_VR. SZ acknowledges the Alexander von Humboldt Foundation for the Humboldt Research Award and the Max Planck Institute for Astrophysics for its hospitality. The post-doctoral contract of IH was funded by Sorbonne Université in the framework of the Initiative Physique des Infinis (IDEX SUPER). LVEK acknowledges the financial support from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (Grant agreement no. 884760, 'CoDEX').

This work used the color map 'arctic' for slice images from the publicly available PYTHON package CMASHER (van der Velden 2020) and the PYTHON package CORECON (Garaldi 2023) for collating neutral fraction observations. Additionally, this work made extensive use of commonly used packages like NUMPY (Harris et al. 2020), MATPLOTLIB (Hunter 2007) and SCIPY (Virtanen et al. 2020).

## DATA AVAILABILITY

The simulation data, and post-analysis scripts used in this work can be shared upon reasonable request to the corresponding author.

## REFERENCES

* Abdurashidova Z. et al., 2022, ApJ, 924, 51

* Acharya A. et al., 2024a, MNRAS, 527, 7835

* Acharya A., Garaldi E., Ciardi B., Ma Q.-b., 2024b, MNRAS, 529, 3793

* Acharya A., Mertens F., Ciardi B., Ghara R., Koopmans L. V. E., Zaroubi S., 2024c, MNRAS, 534, L30

* Adams N. J. et al., 2024, ApJ, 965, 169

* Amon A. et al., 2023, MNRAS, 518, 477

* Anderson L., Pontzen A., Font-Ribera A., Villaescusa-Navarro F., Rogers K. K., Genel S., 2019, ApJ, 871, 144

* Angulo R. E., Pontzen A., 2016, MNRAS, 462, L1

* Arrabal Haro P. et al., 2023a, Nature, 622, 707

* Arrabal Haro P. et al., 2023b, ApJ, 951, L22

* Atek H. et al., 2023, MNRAS, 524, 5486

* Balu S., Greig B., Qiu Y., Power C., Qin Y., Mutch S., Wyithe J. S. B., 2023, MNRAS, 520, 3368

* Barkana R., Loeb A., 2001, Phys. Rep., 349, 125

* Barrera M. et al., 2023, MNRAS, 525, 6312

* Barsode A., Choudhury T. R., 2024, J. Cosmol. Astropart. Phys., 2024, 036

* Basu A., Garaldi E., Ciardi B., 2024, MNRAS, 532, 841

* Becker G. D., Bolton J. S., Madau P., Pettini M., Ryan-Weber E. V., Venemans B. P., 2015, MNRAS, 447, 3402

* Berti M., Spinelli M., Viel M., 2024, MNRAS, 529, 4803

* Bhagwat A., Costa T., Ciardi B., Pakmor R., Garaldi E., 2024, MNRAS, 531, 3406

* Bianco M., Giri S. K., Iliev I. T., Mellema G., 2021, MNRAS, 505, 3982
Bird S., Ni Y., Di Matteo T., Croft R., Feng Y., Chen N., 2022, MNRAS, 512, 3703

Bolton J. S., Haehnelt M. G., Warren S. J., Hewett P. C., Mortlock D. J., Venemans B. P., McMahon R. G., Simpson C., 2011, MNRAS, 416, L70

Bosman S. E. I., Davies F. B., 2024, A&A, 690, A391

Bosman S. E. I. et al., 2022, MNRAS, 514, 55

Bouwens R. J. et al., 2015, ApJ, 803, 34

Bouwens R. J. et al., 2021, AJ, 162, 47

Bouwens R., Illingworth G., Oesch P., Stefanon M., Naidu R., van Leeuwen I., Magee D., 2023a, MNRAS, 523, 1009

Bouwens R. J. et al., 2023b, MNRAS, 523, 1036

Bromm V., Larson R. B., 2004, ARA&A, 42, 79

Carniani S. et al., 2024, Nature, 633, 318

Casavecchia B., Maio U., Péroux C., Ciardi B., 2024, A&A, 689, A106

Chartier N., Wandelt B., Akrami Y., Villaescusa-Navarro F., 2021, MNRAS, 503, 1897

Chornock R., Berger E., Fox D. B., Lunnan R., Drout M. R., Fong W.-F., Laskar T., Roth K. C., 2013, ApJ, 774, 26

Choudhury T. R., Puchwein E., Haehnelt M. G., Bolton J. S., 2015, MNRAS, 452, 261

Choudhury M. et al., 2024, J. Cosmol. Astropart. Physics, 2025, 32

Ciardi B., Madau P., 2003, ApJ, 596, 1

Dark Energy Survey and Kilo-Degree Survey Collaboration, 2023, Open J. Astrophys., 6, 36

Davies F. B. et al., 2018, ApJ, 864, 142

de Belsunce R., Gratton S., Coulton W., Efstathiou G., 2021, MNRAS, 507, 1072

Dijkstra M., Mesinger A., Wyithe J. S. B., 2011, MNRAS, 414, 2139

Dixon K. L., Iliev I. T., Mellema G., Ahn K., Shapiro P. R., 2016, MNRAS, 456, 3011

Eide M. B., Graziani L., Ciardi B., Feng Y., Kakiichi K., Di Matteo T., 2018, MNRAS, 476, 1174

Eide M. B., Ciardi B., Graziani L., Busch P., Feng Y., Di Matteo T., 2020, MNRAS, 498, 6083

Esmerian C. J., Gnedin N. Y., 2021, ApJ, 910, 117

Esmerian C. J., Gnedin N. Y., 2022, ApJ, 940, 74

Esmerian C. J., Gnedin N. Y., 2024, ApJ, 968, 113

Fan X. et al., 2006, AJ, 132, 117

Feldmann R. et al., 2024, MNRAS, 536, 988

Field G. B., 1959, ApJ, 129, 551

Finkelstein S. L. et al., 2015, ApJ, 810, 71

Finlator K., Keating L., Oppenheimer B. D., Davé R., Zackrisson E., 2018, MNRAS, 480, 2628

Fujimoto S. et al., 2023, ApJ, 949, L25

Furlanetto S. R., Oh S. P., Pierpaoli E., 2006, Phys. Rev. D, 74, 103502

Garaldi E., 2023, J. Open Source Softw., 8, 5407

Garaldi E., Kannan R., Smith A., Springel V., Pakmor R., Vogelsberger M., Hernquist L., 2022, MNRAS, 512, 4909

Garaldi E. et al., 2024, MNRAS, 530, 3765

Georgiev I., Mellema G., Giri S. K., 2024, MNRAS, 536, 3689

Ghara R., Choudhury T. R., Datta K. K., 2015, MNRAS, 447, 1806

Ghara R., Mellema G., Giri S. K., Choudhury T. R., Datta K. K., Majumdar S., 2018, MNRAS, 476, 1741

Ghara R. et al., 2020, MNRAS, 493, 4728

Ghara R. et al., 2024, A&A, 687, A252

Ghirardini V. et al., 2024, A&A, 689, A298

Giarè W., Di Valentino E., Melchiorri A., 2024, Phys. Rev. D, 109, 103519

Giri S. K., Mellema G., 2021, MNRAS, 505, 1863

Giri S. K., Mellema G., Aldheimer T., Dixon K. L., Iliev I. T., 2019a, MNRAS, 489, 1590

Giri S. K., D'Aloisio A., Mellema G., Komatsu E., Ghara R., Majumdar S., 2019b, J. Cosmology Astropart. Phys., 2019, 058

Giri S. K., Schneider A., Maion F., Angulo R. E., 2023, A&A, 669, A6

Giri S. K., Bianco M., Schaeffer T., Iliev I. T., Mellema G., Schneider A., 2024, MNRAS, 533, 2364

Gnedin N. Y., 2014, ApJ, 793, 29

Gnedin N. Y., Kaurov A. A., 2014, ApJ, 793, 30

Gorce A., Douspis M., Aghanim N., Langer M., 2018, A&A, 616, A113

Greig B., Mesinger A., Haiman Z., Simcoe R. A., 2017, MNRAS, 466, 4239

Greig B., Mesinger A., Bañados E., 2019, MNRAS, 484, 5094

Greig B., Trott C. M., Barry N., Mutch S. J., Pindor B., Webster R. L., Wyithe J. S. B., 2021a, MNRAS, 500, 5322

Greig B. et al., 2021b, MNRAS, 501, 1

Greig B., Prelogović D., Qin Y., Ting Y.-S., Mesinger A., 2024, MNRAS, 533, 2530

HERA Collaboration et al., 2023, ApJ, 945, 124

Haiman Z., Abel T., Madau P., 2001, ApJ, 551, 599

Han J., Cole S., Frenk C. S., Benitez-Llambay A., Helly J., 2018, MNRAS, 474, 604

Harikane Y. et al., 2022, ApJS, 259, 20

Harikane Y. et al., 2023, ApJS, 265, 5

Harris C. R. et al., 2020, Nature, 585, 357

Hassan S., Andrianomena S., Doughty C., 2020, MNRAS, 494, 5761

Heger A., Woosley S. E., 2002, ApJ, 567, 532

Heintz K. E. et al., 2023a, Science, 384, 890

Heintz K. E. et al., 2023b, Nat. Astron., 7, 1517

Helton J. M. et al., 2024, ApJ, 962, 124

Henriques B. M. B., White S. D. M., Thomas P. A., Angulo R., Guo Q., Lemson G., Springel V., Overzier R., 2015, MNRAS, 451, 2663

Henriques B. M. B., Yates R. M., Fu J., Guo Q., Kauffmann G., Rissawat C., Thomas P. A., White S. D. M., 2020, MNRAS, 491, 5795

Hernández-Aguayo C. et al., 2023, MNRAS, 524, 2556

Hirling P., Bianco M., Giri S. K., Iliev I. T., Mellema G., Kneib J. P., 2024, Astron. Comput., 48, 100861

Hoag A. et al., 2019, ApJ, 878, 12

Hogan C. J., Rees M. J., 1979, MNRAS, 188, 791

Hothi I., Allys E., Semelin B., Boulanger F., 2024, A&A, 686, A212

Hunter J. D., 2007, Comput. Sci. Eng., 9, 90

Hutter A., Dayal P., Yepes G., Gottlöber S., Legrand L., Ucci G., 2021, MNRAS, 503, 3698

Hutter A., Cueto E. R., Dayal P., Gottlöber S., Trebitsch M., Yepes G., 2024, A&A, 694, 15

Iliev I. T., Shapiro P. R., Raga A. C., 2005, MNRAS, 361, 405

Iliev I. T., Mellema G., Ahn K., Shapiro P. R., Mao Y., Pen U.-L., 2014, MNRAS, 439, 725

Jaskot A. E. et al., 2024, ApJ, 973, 111

Jensen H., Laursen P., Mellema G., Iliev I. T., Sommer-Larsen J., Shapiro P. R., 2013, MNRAS, 428, 1366

Jin S. et al., 2023, A&A, 670, L11

Jones T. W. et al., 1998, PASP, 110, 125

Jones G. C. et al., 2024, MNRAS, 536, 2355

Jung I. et al., 2024, ApJ, 967, 73

Kannan R., Garaldi E., Smith A., Pakmor R., Springel V., Vogelsberger M., Hernquist L., 2022, MNRAS, 511, 4005

Kaur H. D., Gillet N., Mesinger A., 2020, MNRAS, 495, 2354

Kern N. S., Liu A., Parsons A. R., Mesinger A., Greig B., 2017, ApJ, 848, 23

Klypin A., Prada F., Byun J., 2020, MNRAS, 496, 3862

Koopmans L., et al., 2015, Proc. Sci., The Cosmic Dawn and Epoch of Reionisation with SKA, SISSA, Trieste, PoS#001

Kostyuk I., Nelson D., Ciardi B., Glatzle M., Pillepich A., 2023, MNRAS, 521, 3077

Kostyuk I., Ciardi B., Ferrara A., 2025, A&A, 695, A32

Kozyreva A., Caputo A., Baklanov P., Mironov A., Janka H.-T., 2025, A&A, 694, A319

La Plante P., Mirocha J., Gorce A., Lidz A., Parsons A., 2023, ApJ, 944, 59

Leethochawalit N. et al., 2023, ApJ, 942, L26

Leung G. C. K. et al., 2023, ApJ, 954, L46

Lewis J. S. W. et al., 2022, MNRAS, 516, 3389

Lidz A., Zahn O., Furlanetto S. R., McQuinn M., Hernquist L., Zaldarriaga M., 2009, ApJ, 690, 252

Liu A., Parsons A. R., 2016, MNRAS, 457, 1864

Long A. S. et al., 2023, ApJ, 970, 14

Longley E. P. et al., 2023, MNRAS, 520, 5016

Looser T. J. et al., 2023, A&A, 697, 16

Ma X. et al., 2018, MNRAS, 478, 1694
Ma X., Quataert E., Wetzel A., Hopkins P. F., Faucher-Giguère C.-A., Kereš D., 2020, MNRAS, 498, 2001

Ma Q.-B., Ciardi B., Eide M. B., Busch P., Mao Y., Zhi Q.-J., 2021, ApJ, 912, 143

Ma Q.-B., Ghara R., Ciardi B., Iliev I. T., Koopmans L. V. E., Mellema G., Mondal R., Zaroubi S., 2023, MNRAS, 522, 3284 (M23)

Ma Q.-B., Chen X.-R., Li M., Guo Q., Ciardi B., Acharya A., Wang X., 2025, ApJ, 986, 14

Madau P., Meiksin A., Rees M. J., 1997, ApJ, 475, 429

Maio U., Viel M., 2023, A&A, 672, A71

Maio U., Péroux C., Ciardi B., 2022, A&A, 657, A47

Maion F., Angulo R. E., Zennaro M., 2022, J. Cosmol. Astropart. Phys., 2022, 036

Mao Y., Tegmark M., McQuinn M., Zaldarriaga M., Zahn O., 2008, Phys. Rev. D, 78, 023529

Mascia S. et al., 2023, A&A, 672, A155

Mason C. A. et al., 2018, ApJ, 857, L11

McGreer I. D., Mesinger A., Fan X., 2011, MNRAS, 415, 3237

McGreer I. D., Mesinger A., D'Odorico V., 2015, MNRAS, 447, 499

McLeod D. J. et al., 2024, MNRAS, 527, 5004

McQuinn M., Zahn O., Zaldarriaga M., Hernquist L., Furlanetto S. R., 2006, ApJ, 653, 815

Mellema G., Iliev I. T., Alvarez M. A., Shapiro P. R., 2006, New A, 11, 374

Mertens F. G. et al., 2020, MNRAS, 493, 1662

Mertens F. G., Bobin J., Carucci I. P., 2024, MNRAS, 527, 3517

Mertens F. G. et al., 2025, A&A, 698, A186

Mesinger A., Furlanetto S., 2007, ApJ, 669, 663

Mesinger A., Furlanetto S., Cen R., 2011, MNRAS, 411, 955

Mesinger A., Aykutalp A., Vanzella E., Pentericci L., Ferrara A., Dijkstra M., 2015, MNRAS, 446, 566

Metha B., Berger S., 2024, PASP, 137, 23

Molaro M., Davé R., Hassan S., Santos M. G., Finlator K., 2019, MNRAS, 489, 5594

Mondal R., Bharadwaj S., Majumdar S., 2017, MNRAS, 464, 2992

Mondal R. et al., 2020, MNRAS, 498, 4178

Mortlock D. J. et al., 2011, Nature, 474, 616

Murray S., 2014, Astrophysics Source Code Library, record ascl:1412.006

Mutch S. J., Geil P. M., Poole G. B., Angel P. W., Duffy A. R., Mesinger A., Wyithe J. S. B., 2016, MNRAS, 462, 250

Nasirudin A., Iliev I. T., Ahn K., 2020, MNRAS, 494, 3294

Ocvirk P. et al., 2016, MNRAS, 463, 1462

Ocvirk P. et al., 2020, MNRAS, 496, 4087

Ono Y. et al., 2012, ApJ, 744, 83

Ota K. et al., 2008, ApJ, 677, 12

Ouchi M. et al., 2010, ApJ, 723, 869

Papovich C. et al., 2023, ApJ, 949, L18

Peacock J. A., 1999, Cosmological Physics, Cambridge Univ. Press, Cambridge

Pentericci L. et al., 2014, ApJ, 793, 113

Planck Collaboration VI, 2020, A&A, 641, A6

Poole G. B., Angel P. W., Mutch S. J., Power C., Duffy A. R., Geil P. M., Mesinger A., Wyithe S. B., 2016, MNRAS, 459, 3025

Qin Y., Mesinger A., Bosman S. E. I., Viel M., 2021, MNRAS, 506, 2390

Riess A. G. et al., 2022, ApJ, 934, L7

Robertson B. E. et al., 2013, ApJ, 768, 71

Robertson B. E. et al., 2023, Nat. Astron., 7, 611

Rosdahl J. et al., 2018, MNRAS, 479, 994

Rubin A. et al., 2016, ApJ, 820, 33

Santos M. G., Cooray A., 2006, Phys. Rev. D, 74, 083517

Santos M., Ferramacho L., Silva M., Amblard A., Cooray A., 2010a, Astrophysics Source Code Library, record ascl:1010.025

Santos M. G., Ferramacho L., Silva M. B., Amblard A., Cooray A., 2010b, MNRAS, 406, 2421

Saxena A., Cole A., Gazagnes S., Meerburg P. D., Weniger C., Witte S. J., 2023, MNRAS, 525, 6097

Schaeffer T., Giri S. K., Schneider A., 2023, MNRAS, 526, 2942

Schauer D., 2002, A&A, 382, 28

Schenker M. A., Ellis R. S., Konidaris N. P., Stark D. P., 2014, ApJ, 795, 20

Schneider A., Schaeffer T., Giri S. K., 2023, Phys. Rev. D, 108, 043030

Schroeder J., Mesinger A., Haiman Z., 2013, MNRAS, 428, 3058

Shaver P. A., Windhorst R. A., Madau P., de Bruyn A. G., 1999, A&A, 345, 380

Sipple J., Lidz A., 2024, ApJ, 961, 50

Smith A., Kannan R., Garaldi E., Vogelsberger M., Pakmor R., Springel V., Hernquist L., 2022, MNRAS, 512, 3243

Sobacchi E., Mesinger A., 2015, MNRAS, 453, 1843

Springel V., White S. D. M., Tormen G., Kauffmann G., 2001, MNRAS, 328, 726

Springel V., Pakmor R., Zier O., Reinecke M., 2021, MNRAS, 506, 2871

Stanway E. R., Eldridge J. J., 2018, MNRAS, 479, 75

Stefanon M., Bouwens R. J., Labbé I., Illingworth G. D., Gonzalez V., Oesch P. A., 2021, ApJ, 922, 29

Sun G., Furlanetto S. R., 2016, MNRAS, 460, 417

Tacchella S., Bose S., Conroy C., Eisenstein D. J., Johnson B. D., 2018, ApJ, 868, 92

Thomas R. M. et al., 2009, MNRAS, 393, 32

Tinker J., Kravtsov A. V., Klypin A., Abazajian K., Warren M., Yepes G., Gottlöber S., Holz D. E., 2008, ApJ, 688, 709

Totani T., Kawai N., Kosugi G., Aoki K., Yamada T., Iye M., Ohta K., Hattori T., 2006, PASJ, 58, 485

Tozzi P., Madau P., Meiksin A., Rees M. J., 2000, ApJ, 528, 597

Trac H., Chen N., Holst I., Alvarez M. A., Cen R., 2022, ApJ, 927, 186

Treu T. et al., 2023, ApJ, 942, L28

Tröster T. et al., 2020, A&A, 633, L10

Trott C. M. et al., 2020, MNRAS, 493, 4711

van der Velden E., 2020, J. Open Source Softw., 5, 2004

Vani A., Ayromlou M., Kauffmann G., Springel V., 2024, MNRAS, 536, 777

Villaescusa-Navarro F. et al., 2018, ApJ, 867, 137

Virtanen P. et al., 2020, Nat. Methods, 17, 261

Vrbanec D. et al., 2016, MNRAS, 457, 666

Watson W. A., Iliev I. T., D'Aloisio A., Knebe A., Shapiro P. R., Yepes G., 2013, MNRAS, 433, 1230

Wu Z., Kravtsov A., 2024, Open J. Astrophys., 7, 56

Zaroubi S., 2013, The Epoch of Reionization. Springer, Berlin, p. 45, https://doi.org/10.1007/978-3-642-32362-1_2

Zhang X., Lachance P., Ni Y., Li Y., Croft R. A. C., Matteo T. D., Bird S., Feng Y., 2024, MNRAS, 528, 281

## APPENDIX A: UV LUMINOSITY FUNCTIONS ACROSS REDSHIFTS

As we tune the astrophysical parameters for the constrained case only to match the UVLFs at $z = 10$ and $9$, it is interesting to explore how the four cosmological models compare to the UVLFs observed at other redshifts. In Fig. A1, we show results at $z = 12, 10, 9, 8, 7, 5$. We note that while the models mildly overestimate observational data at $z = 5$, this could be because astrophysical parameters may need to be evolved with redshift, while here they are kept constant. Interestingly, the $\sigma_8$ low model tends to reproduce the bright end much better at lower redshifts.
<figure><img src="image_11.png" /><figcaption>**Figure A1.** A and B-series averaged UVLFs for the constrained case at $z = 12, 10, 9, 8, 7$, and 5 for the fiducial (cyan), *h* high (red), $\sigma_8$ low (purple), and $\sigma_8$ high (green). We also show the dust attenuated (solid) and unattenuated (dashed) UVLFs and *JWST* and *HST* observations (Bouwens et al. 2015, 2021, 2023a, b; Finkelstein et al. 2015; Harikane et al. 2022, 2023; Leung et al. 2023; Adams et al. 2024; McLeod et al. 2024, grey circles).</figcaption></figure>

## APPENDIX B: FIDUCIAL MODEL WITH

$$F_{\text{esc}} = 25 \%$$

In Fig. B1, we show the redshift evolution of the average of the A and B-series volume-averaged neutral hydrogen fraction $\langle x_{HI} \rangle$ for an escape fraction of $f_{\text{esc}} = 12.5$ per cent (cyan) and 25 per cent (brown). We note that as expected, doubling the escape fraction leads to a faster rate of reionization, getting 50 per cent reionized by $z \approx 7$ for $f_{\text{esc}} = 25$ per cent as compared to $z \approx 6.4$ for $f_{\text{esc}} =$

12.5 per cent. However, both cases are in agreement with observations, except for redshifts $z \le 6$, where the lower $f_{\text{esc}}$ case has slightly higher values. Nevertheless, it is still within the margin of error of observations.

In Fig. B2, we also show the redshift evolution of the normalized 21-cm signal power spectrum ($\Delta_{21\text{cm}}^2$) at $k = 0.15 \text{ hMpc}^{-1}$. We note that just like the other parameters as shown in Fig. 8, the choice of $f_{\text{esc}}$ impacts where $\Delta_{21\text{cm}}^2(k = 0.15 \text{ hMpc}^{-1})$ peaks as well as its overall trend across redshift.
<figure><img src="image_12.png" /></figure>

**Figure B1.** Redshift evolution of the average of the A and B-series volume-averaged neutral hydrogen fraction $\langle x_{HI} \rangle$ for the fiducial model, for an escape fraction of $f_{esc} = 12.5$ per cent (cyan) and 25 per cent (brown). Dotted lines refer to a fit to the curves, which is used for a better estimate of the redshift of reionization when the redshift resolution is too coarse. The vertical grey dashed lines indicate the redshifts observationally relevant for LOFAR ($z = 10.11, 9.16, \text{and } 8.3$), and the black solid line is drawn at $\langle x_{HI} \rangle = 0.5$ to guide the eye. Grey circles are a collection of observational constraints (Fan et al. 2006; Totani et al. 2006; Ota et al. 2008; Ouchi et al. 2010; Bolton et al. 2011; Dijkstra et al. 2011; McGreer et al. 2011, 2015; Mortlock et al. 2011; Ono et al. 2012; Chornock et al. 2013; Jensen et al. 2013; Robertson et al. 2013; Schroeder et al. 2013; Pentericci et al. 2014; Schenker et al. 2014; Choudhury et al. 2015; Mesinger et al. 2015; Sobacchi & Mesinger 2015; Greig et al. 2017, 2019; Davies et al. 2018; Mason et al. 2018; Hoag et al. 2019; Jones et al. 2024, collected in the CORECON module, Garaldi 2023).

<figure><img src="image_13.png" /></figure>

**Figure B2.** Redshift evolution of the A and B-series averaged normalized 21-cm signal power spectrum ($\Delta_{21\text{cm}}^2$) at $k = 0.15 \text{ hMpc}^{-1}$ for the fiducial model, for an escape fraction of $f_{esc} = 12.5$ per cent (cyan) and 25 per cent (brown). The vertical dashed lines indicate the redshifts relevant for LOFAR, i.e. $z = 10.11, 9.16, \text{and } 8.3$.

## APPENDIX C: MAPS OF $\delta T_b$

### C1 B-series maps in the constrained case

In Fig. C1, we show the maps of $\delta T_b$ of the B-series middle slices (the simulations with 'paired' initial conditions in the F&P pair of simulations) for the four cosmological models in the constrained case at $z = 12, 10, 8$, and 6 (from top to bottom). Here, the darkest regions represent the ionized regions with $\delta T_b = 0$. Note that $\delta T_b$ cannot assume negative values, due to the assumption of $T_S \gg T_{\text{CMB}}$. We note that the ionized regions of the A-series middle slices shown in Fig. 3 correspond to the most neutral regions of Fig. C1. This is because of the nature of the Fixed & Paired approach, as regions of matter clustering in the A-series should correspond to voids in the B-series by construction.

In Fig. C2, we show the absolute difference between the maps of $\delta T_b$ of the A and B-series middle slices for the four cosmological models in the constrained case at $z = 12, 10, 8$, and 6 (from top to bottom).

This makes the difference between the fixed simulation and its pair evident, as well as showing that with the progress of reionization, small-scale stochasticity induced by ionized bubbles is a tracer for deviations from the perfect antisymmetry that the system started with. That is, until $z = 8$, the middle slices are effectively random noise maps, as they are primarily governed by the underlying matter distribution. However, by $z = 6$, the impact of the growth of ionized bubbles becomes significant and shows up as features in the difference maps.

### C2 A and B-series maps in the unconstrained case

In Figs C3 and C4, we show the maps of $\delta T_b$ of the A and B-series middle slices for the four cosmological models in the unconstrained case at $z = 12, 10, 8$, and 6 (from top to bottom). Thus, here, we directly see the influence of the choice of cosmological parameters on the 21-cm signal.
<figure><img src="image_14.png" /><figcaption>**Figure C1.** Maps of $\delta T_b$ of the B-series middle slices of single cell thickness (i.e. $\approx 600 \ h^{-1} \ \text{ckpc}$) for the four cosmological models in the constrained case at $z = 12, 10, 8$, and 6 (from top to bottom) corresponding to the A-series slices shown in Fig. 3. Here the dark areas represent the ionized regions with $\delta T_b = 0$. Note that $\delta T_b$ cannot have negative values due to the assumption of $T_S \gg T_{\text{CMB}}$.</figcaption></figure>
<figure><img src="image_15.png" /><figcaption>Figure C2. The absolute difference between the maps of $\delta T_b$ of the A and B-series middle slices shown in Figs 3 and C1.</figcaption></figure>
<figure><img src="image_16.png" /><figcaption>**Figure C3.** Maps of $\delta T_b$ of the A-series middle slices of single cell thickness (i.e. $\approx 600~h^{-1}~\text{ckpc}$) for the four cosmological models in the unconstrained case at $z = 12, 10, 8$, and 6 (from top to bottom). Note that here, due to the higher matter clustering driving up star formation, the $\sigma_8$ high model is already fully ionized by $z = 6$, and the $h$ high model is approaching complete reionization as well. On the other hand, the $\sigma_8$ low case still has large neutral regions.</figcaption></figure>
<figure><img src="image_17.png" /><figcaption>**Figure C4.** Maps of $\delta T_b$ of the B-series middle slices of single cell thickness (i.e. $\approx 600~h^{-1}~\text{ckpc}$) for the four cosmological models in the unconstrained case at $z = 12, 10, 8$, and 6 (from top to bottom).</figcaption></figure>

This paper has been typeset from a TeX/LATEX file prepared by the author.