# Constraints on the state of the IGM at $z \sim 8-10$ using redshifted 21-cm observations with LOFAR

R Ghara, S Zaroubi, B Ciardi, G Mellema, S.K Giri, Florent Mertens, M

Mevius, L.V.E Koopmans, I.T Iliev, A Acharya, et al.

► To cite this version:

R Ghara, S Zaroubi, B Ciardi, G Mellema, S.K Giri, et al.. Constraints on the state of the IGM at $z \sim 8-10$ using redshifted 21-cm observations with LOFAR. Astronomy & Astrophysics - A&A, 2025, 699, pp.A109. 10.1051/0004-6361/202554163. hal-05067491

HAL Id: hal-05067491

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

Submitted on 21 Jul 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished 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_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Constraints on the state of the intergalactic medium at z~8-10
using redshifted 21 cm observations with LOFAR

R. Ghara<fnref n="1" />,*, S. Zaroubi<fnref n="2" />,<fnref n="3" />, B. Ciardi<fnref n="4" />, G. Mellema<fnref n="5" />, S. K. Giri<fnref n="6" />,<fnref n="7" />, F. G. Mertens<fnref n="2" />,<fnref n="8" />, M. Mevius<fnref n="9" />,†,
L. V. E. Koopmans<fnref n="2" />,‡, I. T. Iliev<fnref n="10" />,§, A. Acharya<fnref n="4" />,¶, S. A. Brackenhoff<fnref n="2" />,¶, E. Ceccotti<fnref n="2" />,<fnref n="11" />,||, K. Chege<fnref n="2" />,
I. Georgiev<fnref n="3" />,<fnref n="5" />,||, S. Ghosh<fnref n="2" />, I. Hothi<fnref n="8" />,<fnref n="12" />,||, C. Höfer<fnref n="2" />,||, Q. Ma<fnref n="13" />, S. Munshi<fnref n="2" />,||, A. R. Offringa<fnref n="9" />,||, A. K. Shaw<fnref n="14" />,
V. N. Pandey<fnref n="9" />, S. Yatawatta<fnref n="9" />,||, and M. Choudhury<fnref n="15" />,||

Received 17 February 2025 / Accepted 1 May 2025

## ABSTRACT

The power spectra of the redshifted 21 cm signal from the Epoch of Reionization (EoR) contain information about the ionization and thermal states of the intergalactic medium (IGM) and depend on the properties of the sources that existed during that period. Recently, the LOFAR-EoR Key Science Project team has analysed ten nights of LOFAR high-band data and estimated upper limits on the 21 cm power spectrum at redshifts 8.3, 9.1, and 10.1. Here, we used these upper limit results to constrain the properties of the IGM at those redshifts. We focus on the properties of the ionized and heated regions where the temperature is larger than that of the cosmic microwave background (CMB). We modelled the power spectrum of the 21 cm signal with the code GRIZZLY and used a Bayesian inference framework to explore the source parameters for uniform priors on their ranges. The framework also provides information about the IGM properties in the form of derived parameters. We do not include constraints from other observables except for some very conservative limits on the maximum ionization fraction at those redshifts, which we estimated from the CMB Thomson scattering optical depth. In a model that includes a radio background in excess of the CMB, the 95% (68%) credible intervals of dis-favoured models at redshift 9.1 for the chosen priors correspond to IGM states with an averaged ionization and heated fraction below 0.46 (≤0.05), an average gas temperature below 44 K (4 K), and a characteristic size of the heated region of ≤14 *h*⁻<fnref n="1" /> Mpc (≤3 *h*⁻<fnref n="1" /> Mpc). The 68% credible interval suggests an excess radio background that is more than 100% of the CMB at 1.42 GHz, while the 95% cred-ible interval of the radio background efficiency parameter spans the entire prior range. The behaviour of the credible intervals is similar at all redshifts. The models disfavoured by the LOFAR upper limits are extreme, as they are mainly driven by rare and large ionized or heated regions. We find that the inclusion of upper limits from other radio interferometric observations in the Bayesian analysis significantly increases the number of disfavoured EoR models, thus enhancing the disfavoured credible intervals of the IGM parameters, especially those related to the average gas temperature and size distribution of the heated regions. While our constraints are not yet very strong, more stringent upcoming results from 21 cm observations together with the detection of many high-z galaxies, for example with the James Webb Space Telescope, will strengthen understanding of this crucial phase of the Universe.

**Key words.** radiative transfer – galaxies: formation – galaxies: high-redshift – intergalactic medium – cosmology: theory – dark ages, reionization, first stars

## 1. Introduction

The formation of the first stars and galaxies transformed the Universe from a cold and neutral phase to a hot and highly

ionized one. The era when ultraviolet (UV) radiation from the earliest sources ionized the neutral hydrogen (H) in the inter-galactic medium (IGM) is known as the Epoch of Reionization (EoR). While concrete evidence for the existence of this phase transition exists, its details remain unknown. Several observational probes, such as the CMB Thomson scattering optical

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

<footnote n="2">² Kapteyn Astronomical Institute, University of Groningen, PO Box 800, 9700AV Groningen, The Netherlands</footnote>

<footnote n="3">³ ARCO (Astrophysics Research Center), Department of Natural Sciences, The Open University of Israel, 1 University Road, PO Box 808, Ra'anana 4353701, Israel</footnote>

<footnote n="4">⁴ Max-Planck Institute for Astrophysics, Karl-Schwarzschild-Straße 1, 85748 Garching, Germany</footnote>

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

<footnote n="6">⁶ Van Swinderen Institute for Particle Physics and Gravity, University of Groningen, Nijenborgh 4, 9747 AG Groningen, The Netherlands</footnote>

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

<footnote n="8">⁸ LUX, Observatoire de Paris, PSL Research University, CNRS, Sorbonne Université, F-75014 Paris, France</footnote>

<footnote n="9">⁹ Astron, PO Box 2, 7990 AA Dwingeloo, The Netherlands</footnote>

<footnote n="10">¹⁰ Astronomy Centre, Department of Physics and Astronomy, Pevensey II Building, University of Sussex, Brighton BN1 9QH, UK</footnote>

<footnote n="11">¹¹ INAF - Istituto di Radioastronomia, Via P. Gobetti 101, 40129 Bologna, Italy</footnote>

<footnote n="12">¹² Laboratoire de Physique de l'ENS, ENS, Université PSL, CNRS, Sorbonne Université, Universitée Paris Cité, 75005 Paris, France</footnote>

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

<footnote n="14">¹⁴ Department of Computer Science, University of Nevada, Las Vegas, Nevada 89154, USA</footnote>

<footnote n="15">¹⁵ Center for Fundamental Physics of the Universe, Department of Physics, Brown University, Providence 02914, RI, USA</footnote>

<footnote>* Corresponding author: ghara.raghunath@gmail.com</footnote>
depth (e.g. Planck Collaboration VI 2020), high-redshift ($z$) $Ly\alpha$ emitters (e.g. Hu et al. 2010; Morales et al. 2021; Tang et al. 2023; Witten et al. 2024), the Gunn-Peterson trough, and the $Ly\alpha$ damping wings in high-$z$ quasar spectra (e.g. Fan et al. 2006; Becker et al. 2015; Bañados et al. 2018), have provided limited but useful insight about the evolution of the average neutral fraction.

Currently, the redshifted 21 cm signal from neutral hydrogen is the most promising observational probe of the spatial distribution and temporal evolution of H I in the IGM (see e.g. Madau et al. 1997; Furlanetto et al. 2006; Zaroubi 2013). This signal carries a large amount of information on the physical state of the IGM as well as on the properties of the first UV and X-ray emitting sources. Indeed, as it is a line emission, the 21 cm signal permits one to study physical processes at different stages of the EoR.

To measure the redshifted 21 cm signal, several radio observations have been set up or proposed, and they can be divided into two categories. The first type uses a single radio antenna and attempts to measure the time evolution of the sky-averaged brightness temperature ($\delta T_b$) of the 21 cm signal. Into this category fall instruments such as the Experiment to Detect the Global EoR Signature (EDGES; Bowman & Rogers 2010), EDGES2 (Monsalve et al. 2017; Bowman et al. 2018), the Shaped Antenna measurement of the background RAdio Spectrum (SARAS; Patra et al. 2015), SARAS2 (Singh et al. 2017), the Broadband Instrument for Global HydrOgen Reionization Signal (BigHorns; Sokolowski et al. 2015), the Sonda Cosmologica de las Islas para la Deteccion de Hidrogeno Neutro (SciHi; Voytek et al. 2014), the Large Aperture Experiment to detect the Dark Ages (LEDA; Price et al. 2018), the Radio Experiment for the Analysis of Cosmic Hydrogen (REACH; de Lera Acedo et al. 2022), Probing Radio Intensity at High-Z from Marion (Prizm; Philip et al. 2018), and the Mapper of the IGM Spin Temperature (MIST; Monsalve et al. 2024).

The second type uses large radio interferometers to measure the spatial fluctuations of the redshifted H I signal in a statistical sense, mainly through its power spectrum. To this category belong ongoing or completed observations such as the Low Frequency Array (LOFAR<fnref n="1" />; van Haarlem et al. 2013; Patil et al. 2017), the Giant Metrewave Radio Telescope (GMRT; Paciga et al. 2013), the Precision Array for Probing the Epoch of Reionization (PAPER<fnref n="2" />; Parsons et al. 2014; Kolopanis et al. 2019), the Murchison Widefield Array (MWA<fnref n="3" />; Tingay et al. 2013), the Hydrogen Epoch of Reionization Array (HERA<fnref n="4" />; DeBoer et al. 2017), and the New Extension in Nançay Upgrading LOFAR (NenuFAR<fnref n="5" />; Munshi et al. 2024). Around the start of the next decade, the Square Kilometre Array (SKA)<fnref n="6" /> will become available, and its transformative sensitivity should also enable tomographic imaging of the signal (Mellema et al. 2013, 2015; Koopmans et al. 2015; Ghara et al. 2017).

These observations face a range of challenges. Among them is the mitigation of the galactic and extra-galactic foregrounds, which are four to five orders of magnitude brighter than the expected signal (see e.g. Ghosh et al. 2012). However, the smooth frequency dependence of the foreground signal allows for its subtraction (Harker et al. 2009; Bonaldi & Brown

2015; Chapman et al. 2016; Mertens et al. 2018; Hothi et al. 2021), avoidance (Datta et al. 2010; Liu et al. 2014), or suppression (Datta et al. 2007; Ghara et al. 2016). At the same time, an accurate calibration during the data analysis process is required to minimise the artefacts from bright foreground sources (Barry et al. 2016; Patil et al. 2017). Recent improvements in foreground mitigation algorithms (e.g. Liu et al. 2014; Mertens et al. 2018, 2024; Acharya et al. 2024a; Gao et al. 2024), calibration methods (see e.g. Yatawatta & Banerjee 2011; Kern et al. 2019, 2020; Mevius et al. 2022; Gan et al. 2022, 2023), mitigation of radio frequency interference (see e.g. Offringa et al. 2019; Wilensky et al. 2019; Munshi et al. 2025), and ionospheric effects (e.g. Mevius et al. 2016; Edler et al. 2021), as well as improved sky-model building (see e.g. Patil et al. 2016; Ewall-Wice et al. 2017) and inclusion of polarization effects (e.g. Jelić et al. 2010; Spinelli et al. 2018) in the data analysis step, are showing promising results in reducing the effects of systematics and in helping recover the targeted 21 cm signal.

So far only one detection of $\delta T_b$ with a minimum at $z \approx 17$ (Bowman et al. 2018) has been claimed based on an EDGES2 low-band observation. This detection, though, remains disputed (e.g. in Hills et al. 2018; Draine & Miralda-Escudé 2018; Singh & Subrahmanyan 2019; Bradley et al. 2019), and it requires confirmation from other similar observations. A recent observation in the 55-85 MHz band with SARAS3 is inconsistent with the results from Bowman et al. (2018), ruling out the deep $\delta T_b$ profile around $z \sim 17$ with 95% confidence (Singh et al. 2022). Other global 21 cm signal experiments, such as LEDA (Bernardi et al. 2016), EDGES high-band (Monsalve et al. 2017, 2019), and SARAS2 (Singh et al. 2017), have only provided upper limits on the signal amplitude. Results from SARAS and EDGES, however, have started ruling out EoR scenarios and putting constraints on the early source properties, dark matter models, and the strength of the radio background (e.g. Barkana 2018; Fialkov et al. 2018; Muñoz & Loeb 2018; Nebrin et al. 2019; Chatterjee et al. 2019, 2020; Ghara & Mellema 2020; Ghara et al. 2022; Bera et al. 2023; Acharya et al. 2023; Cyr et al. 2024).

This paper focuses on the results of radio interferometric observations of the 21 cm signal power spectrum $\Delta^2$, which contains much more information than the global signal. At the current time, such observations have only produced upper limits. A $2\sigma$ upper limit of $(248)^2 \text{Mpc}^{-1}$ for $k = 0.5 \text{ h Mpc}^{-1}$ at $z = 8.6$ from the GMRT was the first of its kind (Paciga et al. 2013). Later, PAPER<fnref n="7" /> (Parsons et al. 2014; Ali et al. 2015), LOFAR (Patil et al. 2017; Mertens et al. 2020), MWA (Barry et al. 2019; Trott et al. 2020), and HERA (Abdurashidova et al. 2023) produced additional upper limits. To date, $\Delta^2(k = 0.075 \text{ h Mpc}^{-1}, z = 9.1, 141 \text{ hours}) \approx (73)^2 \text{Mpc}^2$, $\Delta^2(k = 0.34 \text{ h Mpc}^{-1}, z = 7.9, 1128 \text{ hours}) \approx (21.4)^2 \text{Mpc}^2$, and $\Delta^2(k = 0.14 \text{ h Mpc}^{-1}, z = 6.5, 110 \text{ hours}) \approx (43)^2 \text{Mpc}^2$ are the best upper limits obtained from LOFAR (Mertens et al. 2020), HERA (Abdurashidova et al. 2023), and MWA (Trott et al. 2020) EoR observations, respectively. Recent observations with the LOFAR-Low Band Antenna array (Gehlot et al. 2019), the Owens Valley Radio Observatory Long Wavelength Array (OVRO-LWA; Eastwood et al. 2019), and NenuFAR (Munshi et al. 2024) have produced upper limits at

<footnote n="1">¹ http://www.lofar.org/</footnote>

<footnote n="2">² http://eor.berkeley.edu/</footnote>

<footnote n="3">³ https://www.mwatelescope.org/</footnote>

<footnote n="4">⁴ https://reionization.org/</footnote>

<footnote n="5">⁵ https://nenufar OBS-nancay.fr/en/homepage-en/</footnote>

<footnote n="6">⁶ http://www.skatelescope.org/</footnote>

<footnote n="7">⁷ The initial strong upper limits from the PAPER collaboration (Beardsley et al. 2016) were later revised to a weaker upper bound after taking signal loss into account in their analysis (Cheng et al. 2018; Kolopanis et al. 2019).</footnote>
higher redshifts. These data have been used to rule out not only some unusual reionization scenarios requiring, for example, atypical cooling mechanisms or the presence of a radio background in addition to the CMB, but also less exotic ones (see e.g. Ghara et al. 2020; Greig et al. 2021; Mondal et al. 2020; Abdurashidova et al. 2022).

Mertens et al. (2025) have recently published new results from the LOFAR team in the redshift range of 8.3–10.1, with a best $2\sigma$ upper limit of $\Delta^2(k = 0.076 h \text{Mpc}^{-1}, z = 9.1, 140 \text{hours}) \approx (54.3)^2 \text{mK}^2$. This is approximately a factor of two improvement at the same k-scale and redshift in comparison to the previous one from Mertens et al. (2020) and thus is expected to rule out even more EoR scenarios than those already excluded in Ghara et al. (2020), Greig et al. (2021), and Mondal et al. (2020). Additionally, the upper limits at the other two redshifts might help in ruling out even more EoR scenarios.

It should be noted that the 21 cm signal measurements do not directly probe astrophysical sources but rather the ionization and thermal state of the IGM. However, presenting the 21 cm signal observables in terms of IGM-related parameters is challenging (see e.g. Mirocha et al. 2022; Ghara et al. 2024; Choudhury et al. 2024). As the 21 cm signal simulation codes use the source parameters as input, previous studies mainly focused on constraining the astrophysical source properties using either Fisher matrices or Bayesian inference techniques (e.g. Ewall, Wice et al. 2016; Park et al. 2019; Cohen et al. 2020). It should be noted that the inference on the astrophysical source properties is limited by the assumptions of the source model used in the 21 cm signal simulation codes. We therefore put less emphasis on the inference results for the source parameters and focus on the IGM properties instead.

In this paper, we use the new LOFAR upper limits to investigate which EoR scenarios are disfavoured within the redshift range of 8.3–10.1. We achieved this by employing the Bayesian inference framework previously used in Ghara et al. (2021) to interpret multi-redshift upper limits from MWA observations (Trott et al. 2020). The main difference between this version of the framework and the one used in Ghara et al. (2020) is that we consider the presence of an excess radio background component in addition to the CMB, similar to what was done in Mondal et al. (2020).

The paper is structured as follows. In Sect. 2, we describe the new LOFAR upper limits used in this study and introduce the basic methodology of our Bayesian interpretation framework. We present our results in Sect. 3, and finally we conclude our study in Sect. 4. We adopted the cosmological parameters $\Omega_m = 0.27$, $\Omega_\Lambda = 0.73$, $\Omega_B = 0.044$, and $h = 0.7$ (Wilkinson Microwave Anisotropy Probe WMAP; Hinshaw et al. 2013), which are the same as used in the N-body simulations employed here.

## 2. Methodology

Here we briefly describe the new LOFAR upper limits on the 21 cm signal power spectrum. We also give an overview of the framework employed to constrain the IGM properties with them.

### 2.1. Observations

The LOFAR-EoR Key Science Project team has analysed data for the North Celestial Pole (NCP) field from 14 nights of observations taken during LOFAR Cycles 0–3. Previously, Mertens et al. (2020) analysed 12 nights (resulting in 141 hours) of NCP data in the frequency range 134–146 MHz, and presented

upper limits on the 21 cm signal power spectrum at $z = 9.1$. In the latest results, the team has added two more nights at 134–146 MHz ($z = [8.73, 9.6]$), and also analysed two additional frequency bands: 122–134 MHz ($z = [9.6, 11.64]$) and 146–158 MHz ($z = [8, 8.73]$), while only the best 10 nights (corresponding to 140 hours) of data for each redshift bin are finally combined. This has resulted in new upper limits on the 21 cm signal power spectrum at $z = 8.3, 9.1$ and 10.1.

The details of the observations and the analysis techniques can be found in Mertens et al. (2025). Here we summarise the final results in Table 1. It shows the mean power spectrum $\Delta_{21}^2(k, z)$ and associated $1\sigma$ error $\Delta_{21,\text{em}}^2(k, z)$. The best result is at $z = 9.1$, with a $2\sigma$ upper limit of $(54.3 \text{ mK})^2$ at $k = 0.075 h \text{Mpc}^{-1}$, which is a factor of ≈2 lower than the one reported in Mertens et al. (2020). This is caused by improvements in the data analysis pipeline (see Mertens et al. 2025 for more details).

### 2.2. Interpretation framework

The framework employed to constrain the IGM properties using the upper limit results from Mertens et al. (2025) relies on a large (~$10^5$) database of reionization scenarios, as well as the associated IGM properties and 21 cm signal power spectra. These theoretical power spectra are then used in a Bayesian analysis in combination with the measured ones to find the exclusion likelihood for a reionization scenario. The same approach was used in Ghara et al. (2021) to interpret upper limits from observations with MWA. We refer the reader to Ghara et al. (2020, 2021) for a detailed description of the framework.

The large database of reionization scenarios and associated 21 cm signals is produced with GRIZZLY (Ghara et al. 2015a, 2018), an independent implementation of the BEARS algorithm (Thomas & Zaroubi 2008; Thomas et al. 2009; Thomas & Zaroubi 2011; Krause et al. 2018). It requires cosmological density and velocity fields, together with the dark matter halo catalogues as input. These are typically obtained from a cosmological N-body simulation. Given a set of astrophysical source parameters, it then employs a one-dimensional radiative transfer scheme to generate ionization fraction, gas temperature and $Ly\alpha$ flux fields. The N-body simulation as well as the source model used in GRIZZLY are described in the following sub-sections.

#### 2.2.1. N-body simulation

The N-body outputs used here are the same employed in Ghara et al. (2020) and were obtained from the PRACE$^8$ project PRACE4LOFAR. The dark matter-only simulation was run using the CUBEP$^3$M (Harnois-Déraps et al. 2013) code with particle mass resolution of $4.05 \times 10^7 \text{ M}_\odot$. The simulation cube has a size of $500 h^{-1} \text{ Mpc}$, corresponding to a field-of-view of $4.3^\circ \times 4.3^\circ$ at $z \approx 9$. This length is enough to cover the largest scale probed by the LOFAR observations, which corresponds to a k-scale of $0.05 h \text{Mpc}^{-1}$, requiring a size of at least ≈250 $h^{-1} \text{ Mpc}$.

The density and velocity fields were smoothed onto a $300^3$ grid (see e.g. Giri et al. 2019a,b). The lists of dark matter halos were obtained on the fly during the N-body simulation using a spherical overdensity halo finder (Watson et al. 2013). The halo

<footnote>$^8$ Partnership for Advanced Computing in Europe: http://www.prace-ri.eu/</footnote>
<table><caption>**Table 1.** Recent LOFAR upper limit on $\Delta^2$ at different k-bins for $z = 8.3, 9.1$, and $10.1$.</caption><thead><tr><th colspan="3">z = 8.3</th><th colspan="3">z = 9.1</th><th colspan="3">z = 10.1</th></tr><tr><th>k (h Mpc<sup>-1</sup>)</th><th>Δ<sup>2</sup><sub>21</sub> (mK<sup>2</sup>)</th><th>Δ<sup>2</sup><sub>21, err</sub> (mK<sup>2</sup>)</th><th>k (h Mpc<sup>-1</sup>)</th><th>Δ<sup>2</sup><sub>21</sub> (mK<sup>2</sup>)</th><th>Δ<sup>2</sup><sub>21, err</sub> (mK<sup>2</sup>)</th><th>k (h Mpc<sup>-1</sup>)</th><th>Δ<sup>2</sup><sub>21</sub> (mK<sup>2</sup>)</th><th>Δ<sup>2</sup><sub>21, err</sub> (mK<sup>2</sup>)</th></tr></thead><tbody><tr><td>0.083</td><td>(47.4)<sup>2</sup></td><td>(31.7)<sup>2</sup></td><td>0.076</td><td>(39.7)<sup>2</sup></td><td>(24.1)<sup>2</sup></td><td>0.076</td><td>(50.7)<sup>2</sup></td><td>(28.5)<sup>2</sup></td></tr><tr><td>0.106</td><td>(77.7)<sup>2</sup></td><td>(31.3)<sup>2</sup></td><td>0.101</td><td>(62.1)<sup>2</sup></td><td>(26.3)<sup>2</sup></td><td>0.101</td><td>(75.5)<sup>2</sup></td><td>(36.3)<sup>2</sup></td></tr><tr><td>0.138</td><td>(110.1)<sup>2</sup></td><td>(35.8)<sup>2</sup></td><td>0.133</td><td>(96.8)<sup>2</sup></td><td>(36.8)<sup>2</sup></td><td>0.133</td><td>(112.4)<sup>2</sup></td><td>(46.8)<sup>2</sup></td></tr><tr><td>0.184</td><td>(158.4)<sup>2</sup></td><td>(46.6)<sup>2</sup></td><td>0.181</td><td>(150.1)<sup>2</sup></td><td>(53.8)<sup>2</sup></td><td>0.181</td><td>(172.4)<sup>2</sup></td><td>(70.0)<sup>2</sup></td></tr><tr><td>0.242</td><td>(234.8)<sup>2</sup></td><td>(62.4)<sup>2</sup></td><td>0.240</td><td>(217.3)<sup>2</sup></td><td>(73.2)<sup>2</sup></td><td>0.240</td><td>(243.4)<sup>2</sup></td><td>(93.1)<sup>2</sup></td></tr><tr><td>0.323</td><td>(303.3)<sup>2</sup></td><td>(67.4)<sup>2</sup></td><td>0.319</td><td>(283.6)<sup>2</sup></td><td>(80.9)<sup>2</sup></td><td>0.323</td><td>(302.8)<sup>2</sup></td><td>(98.3)<sup>2</sup></td></tr><tr><td>0.436</td><td>(338.1)<sup>2</sup></td><td>(80.8)<sup>2</sup></td><td>0.432</td><td>(320.3)<sup>2</sup></td><td>(84.3)<sup>2</sup></td><td>0.434</td><td>(351.2)<sup>2</sup></td><td>(94.4)<sup>2</sup></td></tr></tbody></table>

Notes. The different columns show k-bins, mean value $\Delta^2_{21}(k, z)$ and corresponding $1\sigma$ error $\Delta^2_{21, err}(k, z)$. These are achieved from 140 hours of the NCP field data (see Mertens et al. 2025 for details).

lists consist of halos with more than 25 particles, which leads to a minimum dark matter halo mass $\approx 10^9$ M$_\odot$.

### 2.2.2. GRIZZLY source model

The next step is to employ GRIZZLY (Ghara et al. 2015a, 2018) to generate 3D cubes of the 21 cm signal at the probed redshifts. For this, we used different combinations of astrophysical parameters, which encode the source properties in terms of the emissivity of UV, X-ray, Ly$\alpha$, and radio photons, as each of these affects the spatial fluctuations of the signal (see e.g. Ross et al. 2019; Eide et al. 2020; Reis et al. 2020). Our source model assumes that the stellar mass $M_\star$ in a dark matter halo with mass $M_{\text{halo}}$ is $M_\star = f_\star \frac{\Omega_B}{\Omega_m} M_{\text{halo}}$, where $f_\star$ is the star formation efficiency. In all our simulations we used $f_\star = 0.02$, which is consistent with studies such as Behroozi & Silk (2015) and Sun & Furlanetto (2016).

The source parameters we varied to create our set of simulations are as follows:

(1) Ionization efficiency ($\zeta$): This parameter regulates the emission rate of ionizing photons escaping into the IGM from a halo$^9$. We assume it to be constant per unit stellar mass, or in other words, this emission rate scales linearly with halo mass. The rate is then given by $\dot{N}_i = \zeta \times 2.85 \times 10^{45} \text{ s}^{-1} \text{ M}_\odot^{-1}$, where we vary $\log_{10}(\zeta)$ in the range $[-3, 3]$.

(2) X-ray heating efficiency ($f_X$): This is the parameter that regulates the emission rate of X-ray photons from a halo according to $\dot{N}_X = f_X \times 10^{42} \text{ s}^{-1} \text{ M}_\odot^{-1}$, where we consider the X-ray band to span from 100 eV to 10 keV. Just as the ionizing ultraviolet radiation, the X-ray emission scales linearly with halo mass. We adopt a power-law spectral energy distribution (SED), i.e. $I_X(E) \propto E^{-\alpha}$, where $\alpha = 1.2$ across the X-ray band. The $\dot{N}_X$ value for $f_X = 1$ is in agreement with the SED of high-mass X-ray binaries in local star-forming galaxies for the 0.5–8 keV band (Mineo et al. 2012; Islam et al. 2019). We vary $\log_{10}(f_X)$ in the range $[-3, 3]$.

(3) Minimum mass of UV emitting halos ($M_{\text{min}}$): This parameter determines the number density of UV emitting sources, as we assume that only dark matter halos with a mass larger than $M_{\text{min}}$ contribute to the ionizing photon budget. We vary $\log_{10}(M_{\text{min}}/M_\odot)$ in the range [9, 12]. The lower value is determined by the lowest mass halos resolved in our N-body simulations. While the number density of dark matter halos in a

$\Lambda$CDM cosmology is larger at the low mass end, star formation in halos with mass $\lesssim 10^9$ M$_\odot$ is expected to be significantly suppressed due to thermal feedback, SNe feedback, or inefficient gas accretion (see e.g. Shapiro et al. 1994; Barkana & Loeb 2001; Springel & Hernquist 2003; Sobacchi & Mesinger 2013). These lower mass halos are therefore expected to have a less significant impact on the 21 cm signal power spectrum.

(4) Minimum mass of X-ray emitting halos ($M_{\text{min},X}$): This determines the number density of X-ray emitting halos, as we assume that only halos with a mass larger than $M_{\text{min},X}$ host X-ray sources. Also for $\log_{10}(M_{\text{min},X}/M_\odot)$ we explore the range [9, 12].

(5) Radio background efficiency ($A_r$): We assumed the existence of a radio background$^{10}$ in addition to the CMB. This parameter quantifies the effective brightness temperature of the total radio background as (Fialkov & Barkana 2019; Mondal et al. 2020)

$$T_{\gamma, \text{eff}} = T_{\gamma} \left[ 1 + A_{r} \left( \frac{\nu_{\text{obs}}}{78 \text{ MHz}} \right)^{-2.6} \right]. \quad (1)$$

Here, $T_\gamma(z) = 2.725(1+z)$ K is the usual CMB temperature and $\nu_{\text{obs}}$ is the observation frequency. For simplicity, we adopt a uniform radio background, but the presence of inhomogeneities could enhance the amplitude of the 21 cm signal power spectrum above our estimates (Reis et al. 2020). We vary $A_r$ in the range [0, 416]. The lower bounds correspond to no excess radio background (i.e. $T_{\gamma, \text{eff}} = T_\gamma$), while the upper bounds represent the upper limits of 0.603 K as measured by LWA1 (Dowell & Taylor 2018).

### 2.2.3. Simulating the 21 cm signal power spectrum

The 21 cm signal brightness temperature ($\delta T_b$) can be written as (see e.g. Madau et al. 1997; Furlanetto et al. 2006)

$$\delta T_b(\mathbf{x}, z) = 27 \, x_{\text{HI}}(\mathbf{x}, z) \left[ 1 + \delta_B(\mathbf{x}, z) \left( 1 - \frac{T_{\gamma, \text{eff}}(z)}{T_S(\mathbf{x}, z)} \right) \times \left( \frac{\Omega_B h^2}{0.023} \right) \left( \frac{0.15}{\Omega_m h^2} \frac{1+z}{10} \right)^{1/2} \left( \frac{H}{d v_r / d r + H} \right) \right] \text{mK.} \quad (2)$$

<footnote>$^{10}$ The existence of a radio background of astrophysical or cosmological origin in addition to the CMB is motivated by the results of ARCADE2 (Fixsen et al. 2011) and LWA1 (Dowell & Taylor 2018) observations, which show evidence of excess towards the Rayleigh-Jeans part of the spectrum. In particular, the LWA1 measurement at frequency 40–80 MHz suggests a power-law dependence with spectral index –2.58 ± 0.05 and a temperature of $603^{+102}_{-92}$ mK at 1.42 GHz. We note that the ARCADE2 and LWA1 measurements are debated, as the measurements might suffer from erroneous Galactic modelling (Subrahmanyan & Cowsik 2013).</footnote>

<footnote>$^9$ The amount of ionizing photons reaching the IGM depends on various quantities such as the star formation rate, the stellar initial mass function, the intrinsic ionizing photons emission rate, and the escape fraction. The parameter $\zeta$ combines the effects of all these.</footnote>
Here, $[1 + \delta_B(x, z)]$ is the overdensity at position x and redshift z, extracted from the gridded density fields obtained from the N-body simulation described in Sect. 2.2.1. The quantities $x_{HI}$, $T_S$, H and $dv_r/dr$ are the neutral hydrogen fraction, the spin temperature of hydrogen in the IGM, the Hubble parameter and the velocity gradient along the line of sight, respectively. Given a set of astrophysical source parameters, GRIZZLY generates $x_{HI}$ and gas temperature ($T_K$) maps. Here, we assume $T_S = T_K$ which is expected at the redshifts of interest due to the presence of a strong Ly$\alpha$ background (e.g. Barkana & Loeb 2005; Pritchard & Loeb 2012; Schneider et al. 2021). We include the line-of-sight velocity field effects (so-called redshift space distortions, RSD) to the $\delta T_b$ cubes using the cell movement method (or Mesh-to-Mesh Real-to-Redshift-Space-Mapping scheme; Mao et al. 2012; Ghara et al. 2015b; Ross et al. 2021). Finally, we estimate the dimensionless power spectrum of $\delta T_b$, denoted as $\Delta^2(k)^{11}$.

### 2.2.4. Derived IGM parameters

Similar to Ghara et al. (2020, 2021), our main focus here is on the characterisation of the IGM, as its properties are directly linked to the observed 21 cm signal. However, our framework is not completely independent of source models, as GRIZZLY uses a set of source parameters as input to simulate ionization and temperature fields, which are then used to determine a set of quantities related to the IGM.

When considering the IGM parameters, we defined 'heated regions' as those with $T_K > T_{\gamma,eff}$. To derive the size distribution of these regions, we used the mean free path (MFP) method from Mesinger & Furlanetto (2007), see Ghara et al. (2020) for more details. We then considered the following derived IGM parameters or quantities<fnref n="12" /> for each redshift:

- $\bar{x}_{HII}$: Volume-averaged ionized fraction.

- $\bar{T}_K$ (K): Volume-averaged gas temperature of neutral regions, where we consider the gas to be neutral if $x_{HII} < 0.5$. We note that choosing a different value would change the estimate of $\bar{T}_K$. In general, a higher $x_{HII}$ threshold results in a larger value of $\bar{T}_K$.

- $\delta\bar{T}_b$ (mK): Mass-averaged differential brightness temperature.

- $J_{\text{heat}}$: Volume fraction of heated regions ($T_K > T_{\gamma,eff}$).

- $R_{\text{peak}}^{\text{heat}}$ ($h^{-1}$ Mpc): Characteristic size of the heated regions, corresponding to the size at which the distribution of heated regions size peaks.

- $\Delta R_{\text{FWHM}}^{\text{heat}}$ ($h^{-1}$ Mpc): Full width at half maximum (FWHM) of the PDF of the size distribution of heated regions.

### 2.2.5. Bayesian inference

To explore the source parameter space as given in Table 2, we use the publicly available Markov chain Monte Carlo (MCMC) sampler code COBAYA<fnref n="13" /> (Lewis & Bridle 2002). The likelihood of a set of parameters $\theta$ to be excluded at z for an upper limit

of the power spectrum $\Delta_{21}^2(k, z) \pm \Delta_{21,\text{err}}^2(k, z)$, can be defined as (see Ghara et al. 2020, for details)

$$ \mathcal{L}_{\text{ex,single-z}}(\theta, z) = 1 - \prod_i \frac{1}{2} \left[ 1 + \text{erf} \left( \frac{\Delta_{21}^2(k_i, z) - \Delta_m^2(k_i, \theta, z)}{\sqrt{2}\sigma(k_i, z)} \right) \right]. \quad (3) $$

Here, $\Delta_m^2(k, \theta, z)$ is the model power spectrum for a set of parameters $\theta$, while $\sigma$ is the error that includes contributions from the observational error ($\Delta_{21,\text{err}}^2$) as well as from the power spectrum modelling error ($\Delta_{m,\text{err}}^2$): $\sigma = \sqrt{\Delta_{21,\text{err}}^4 + \Delta_{m,\text{err}}^4}$.

Unlike our previous inference studies, here we also consider the combination of the three upper limits into a joint likelihood:

$$ \mathcal{L}_{\text{ex,joint-z}}(\theta) = 1 - \prod_{i,j} \frac{1}{2} \left[ 1 + \text{erf} \left( \frac{\Delta_{21}^2(k_i, z_j) - \Delta_m^2(k_i, \theta, z_j)}{\sqrt{2}\sigma(k_i, z_j)} \right) \right]. \quad (4) $$

We note that multiple upper limits can be combined into a single likelihood only when the redshift evolution of the source parameters is well understood. As here we assume the source parameters to be redshift-independent, some caution should be applied in the interpretation of the joint analysis. In the future, we plan to extend it by employing the code POLAR (Ma et al. 2023; Acharya et al. 2024b), which includes a redshift evolution in the modelling of the source parameters.

Similar to Ghara et al. (2021), our framework does not directly use GRIZZLY in the MCMC analysis. Instead, we first employ GRIZZLY to generate $10^5$ sets of power spectra and IGM parameters at each of the three redshifts by sampling the 5D source parameter space<fnref n="14" />. Then, in the MCMC analysis we employ a linear interpolation scheme (Weiser & Zarantonello 1988; Virtanen et al. 2020) over the 5D source parameter space, which uses the previously generated power spectra and IGM parameters. We also consider a conservative 30% modelling error, $\Delta_{m,\text{err}}^2(k_i, z) = 0.3 \times \Delta_m^2(k_i, \theta, z)$. This includes the inherent modelling errors due to the approximate nature of the GRIZZLY algorithm<fnref n="15" />, as well as the errors introduced by the interpolation scheme. To run the MCMC on the source parameters, we used ten walkers and $10^6$ steps, and ensure that all the MCMC chains converge well before the final step. We adopt a uniform prior in log-scale on the explored parameter ranges (see Table 2). COBAYA produces the list of corresponding derived IGM parameters using the same linear interpolation scheme.

As the upper limits towards the largest k values exceed the values for all of the power spectra in our simulation database (see next section), they do not yield any constraints. Therefore, only the three smallest k-bins (see Table 1) are used in the MCMC analysis to estimate $\mathcal{L}_{\text{ex,single-z}}(\theta, z)$ and $\mathcal{L}_{\text{ex,joint-z}}(\theta)$.

<footnote n="11">¹¹ We note that in our previous studies (Ghara et al. 2020, 2021) we ignored the redshift space distortions, which results in underestimating the 21 cm signal power spectrum, especially when it is dominated by density fluctuations.</footnote>

<footnote n="12">¹² Unlike in Ghara et al. (2020), here we do not consider separate probability distribution functions (PDFs) for the H II regions. In the absence of X-ray heating, the H II regions are equivalent to the heated regions, as the gas temperature within them is ~$10^4$ K, i.e. much larger than $T_{\gamma,eff}$.</footnote>

<footnote n="13">¹³ https://cobaya.readthedocs.io/en/latest/index.html</footnote>

<footnote n="14">¹⁴ Each source parameter, $s_i$, is sampled by linearly binning the range $[\log_{10}(s_{i,min}), \log_{10}(s_{i,max})]$.</footnote>

<footnote n="15">¹⁵ A comparison study of GRIZZLY with the 3D radiative transfer scheme C²RAY for different reionization scenarios assuming $T_S \gg T_\gamma$ showed an excellent match of the power spectra, having a less than 10% difference (Ghara et al. 2018). However, the level of agreement across the entire parameter space is not known. Also, we do not have robust accuracy estimates of GRIZZLY for scenarios with $T_S$ fluctuations. Thus we used a conservative modelling error which also covers other uncertainties such as an interpolation error in the MCMC analysis. We note that the chosen modelling error is still significantly smaller than the $1\sigma$ error on the LOFAR upper limits at the large scales.</footnote>
<table><caption>**Table 2.** Overview of the source parameters used in GRIZZLY with their explored ranges.</caption><thead><tr><th>Source Parameters</th><th>Description</th><th>Explored range</th></tr></thead><tbody><tr><td>ζ</td><td>Ionization efficiency</td><td>[10<sup>-3</sup>, 10<sup>3</sup>]</td></tr><tr><td>M<sub>min</sub></td><td>Minimum mass of UV emitting halos</td><td>[10<sup>9</sup> M<sub>☉</sub>, 10<sup>12</sup> M<sub>☉</sub>]</td></tr><tr><td>M<sub>min,X</sub></td><td>Minimum mass of X-ray emitting halos</td><td>[10<sup>9</sup> M<sub>☉</sub>, 10<sup>12</sup> M<sub>☉</sub>]</td></tr><tr><td>f<sub>x</sub></td><td>X-ray heating efficiency</td><td>[10<sup>-3</sup>, 10<sup>3</sup>]</td></tr><tr><td>A<sub>r</sub></td><td>Excess radio background efficiency</td><td>[0, 416]</td></tr></tbody></table>

Notes. We note that the parameter space exploration was done on the log-spaced grids for each parameter.

In this analysis we also use a prior on $\bar{x}_{HII}$, derived from the measured Thomson scattering optical depth towards the CMB. For each of the three redshifts $z_i$ we assume a stepwise reionization history, in which the Universe is neutral until $z_i$, partially ionized at a level ($\bar{x}_{HII,max}$) until redshift 6, and fully ionized after that<fnref n="16" />. The value $\bar{x}_{HII,max}$ follows from assuming the $1\sigma$ $\tau$ measurement, i.e. $\tau = 0.061$ for $\tau = 0.054 \pm 0.007$ (Planck Collaboration VI 2020). It represents the largest possible of value of $\bar{x}_{HII}$ at redshift $z_i$ consistent with the $\tau$ measurements. We obtain values of $\bar{x}_{HII,max} = 1, 0.79$ and $0.57$ at $z = 8.3, 9.1$ and $10.1$, respectively. The same approach was used in Ghara et al. (2020, 2021). We note that in a realistic reionization scenario, the evolution of $\bar{x}_{HII}$ is expected to be gradual.

It should be noted that the Bayesian analysis explores the source parameters rather than the IGM parameters, as our framework uses the former as input. The source parameter values obtained from the MCMC chains are then used to generate lists of corresponding IGM parameter values. Finally, these lists are employed to obtain the posterior distribution of the IGM parameters. Because of this procedure, the constraint on the IGM parameters might be significantly affected by the priors on the source parameters. We discuss this in Appendix A.

## 3. Results

Before describing the quantitative analysis of our results, we present in Fig. 1, as an illustration, a subset of the 21 cm signal power spectra with $A_r = 0$ at $z = 8.3, 9.1$ and $10.1$. We find that at each redshift some of the power spectra (indicated with coloured lines) have values larger than $\Delta_{21}^2(k, z) - \Delta_{21,\text{err}}^2(k, z)$ from the recent LOFAR $1\sigma$ upper limits from Mertens et al. (2025, indicated by the down arrow points) in at least one k-bin. For this reason, they have a high probability of being ruled out. Overall, we find that 10%, 13%, and 5% of our simulated power spectra fall into this category at $z = 8.3, 9.1$ and $10.1$, respectively. Here it should be kept in mind that these fractions depend on the priors chosen for the source parameters. On the other hand, the power spectra in grey have a low exclusion probability as they remain below that limit. The thick dashed lines correspond to the power spectra of a completely neutral and unheated IGM, obtained assuming constant values of $x_{HI} = 1$ and $T_S = T_K = 0.021 \times (1+z)^2$ K. We note that the observed upper limits are still unable to rule out this scenario. We estimate that for $z = 9.1$ at least 2–3 times more integration time would be required to achieve that.

<figure><img src="image_2.png" /><figcaption>**Fig. 1.** Set of 1000 power spectra randomly chosen from the $10^5$ simulated power spectra. Panels from top to bottom refer to $z = 8.3, 9.1$, and $10.1$, respectively. These correspond to the scenario with no additional radio background other than the CMB ($A_r = 0$). The down arrow points refer to the recent LOFAR $1\sigma$ upper limit ($\Delta_{21}^2(k, z) \pm \Delta_{21,\text{err}}^2(k, z)$) from Mertens et al. (2025). As a reference, the dashed line corresponds to the power spectrum for a completely neutral IGM with no X-ray heating, and with $T_S = T_K$. The coloured power spectra are those with a value larger than $\Delta_{21}^2(k, z) - \Delta_{21,\text{err}}^2(k, z)$ in at least one k-bin, i.e. they have a high probability of being ruled out, while those in grey have a low exclusion probability.</figcaption></figure>

As a further illustration, Fig. 2 shows the distribution of source parameters for models with a high probability of being ruled out, i.e. which have a value larger than $\Delta_{21}^2(k, z) - \Delta_{21,\text{err}}^2(k, z)$ in at least one k-bin. The vertical axes represent the number of disfavoured models out of the total of $10^5$. We only consider the $A_r = 0$ scenario and each source parameter range has been logarithmically binned into eight bins. We note that a significant fraction of these disfavoured models has high values of $M_{\text{min}}$ and $M_{\text{min,X}}$, indicating that they represent scenarios in which the ionized or heated regions are rare. A large fraction of models also has high values of $\zeta$ and $f_X$, suggesting that the rare ionized or heated regions are large in size.

In Fig. 3, we present the distribution of IGM parameters for the same set of models considered in Fig. 2. A significant number

<footnote n="16">¹⁶ The exact redshift of the end of reionization ($z_{\text{end}}$) is debated. While the Gunn-Peterson trough in high-$z$ quasar spectra suggests $z_{\text{end}} \approx 6$ (e.g. Fan et al. 2006; Mortlock et al. 2011; Venemans et al. 2015; Bañados et al. 2018), Ly$\alpha$ forest observations at lower redshifts suggest a late reionization (see e.g. Becker et al. 2018; Eilers et al. 2018).</footnote>
<figure><img src="image_3.png" /><figcaption>**Fig. 2.** Histogram of the disfavoured values of the source parameters for the scenario with $A_r = 0$. These correspond to the models with power spectra values larger than $\Delta_{21}^2(k, z) - \Delta_{21,\text{err}}^2(k, z)$ in at least one k-bin, and they have a high probability of being ruled out by the recent LOFAR $1\sigma$ upper limits ($\Delta_{21}^2(k, z) \pm \Delta_{21,\text{err}}^2(k, z)$) from Mertens et al. (2025).</figcaption></figure>

of disfavoured models at all redshifts have $\bar{x}_{\text{HII}}$ and $f_{\text{heat}}$ smaller than ~0.6. The characteristic sizes in these models are several tens of $h^{-1}$ Mpc, indicating that the ionized or heated regions are both rare and large. At the same time, the volume occupied by them remains less than 50% of the total volume. We note that a large number of disfavoured models with $\bar{x}_{\text{HII}} \le 0.1$, $f_{\text{heat}} \le 0.1$ and $\bar{T}_{\text{K}} \le 4$ K are strongly influenced by the priors chosen on the source parameters (see Appendix A for details).

In the following sections, we present the results of our MCMC analysis. Following the same approach adopted in Ghara et al. (2021), we first perform the MCMC using the likelihood $\mathcal{L}_{\text{ex,single-z}}(\theta, z)$ as given in Eq. (3), which accounts for upper limits at a single redshift only. We denote this as ‘Single-z’ analysis. We also run an MCMC analysis where we use the joint likelihood $\mathcal{L}_{\text{ex,joint-z}}(\theta)$ as shown in Eq. (4), which combines all the upper-limit results. We denote this as ‘Joint-z’ analysis. We note again that as our source parameters are assumed to be redshift independent, the results of the Joint-z analysis should be interpreted with some caution. We first consider the case of no additional radio background and follow it by the case in which we allow for an additional radio background. Finally we describe how the inclusion of upper limits on the power spectrum from other 21 cm experiments impacts the results.

### 3.1. Analysis without excess radio background ($A_r = 0$)

For $A_r = 0$ our parameter space is 4D, defined by $\zeta, M_{\text{min}}, M_{\text{min},X}$ and $f_X$. This scenario corresponds to one in which the only radio background at 1420 MHz is the CMB.

Figure 4 shows the posterior distribution of the GRIZZLY source parameter space at $z = 9.1$. Here we highlight that the likelihood used in the MCMC represents the exclusion probability of a model. We also note that the 68% and 95% credible interval contours on the two-dimensional contour panels have remained open, indicating that the constraints on the source

<figure><img src="image_4.png" /><figcaption>**Fig. 3.** Histogram of the disfavoured values of the IGM parameters for the scenario with $A_r = 0$. These correspond to the same set of models considered in Fig. 2.</figcaption></figure>

parameters depend on the chosen prior ranges. We find that the 68% disfavoured credible intervals<fnref n="17" /> from the Single-z analysis at $z = 9.1$ are $M_{\text{min}} \ge 3 \times 10^{10} \text{ M}_\odot$ and $M_{\text{min},X} \ge 6 \times 10^{10} \text{ M}_\odot$, while the limits for $\zeta$ and $f_X$ span the entire parameter ranges. The 95% disfavoured credible intervals of the source parameters span the entire parameter ranges for both the Single-z and Joint-z analysis. The 68% credible intervals for the Joint-z analysis case are similar to those of the Single-z case. One reason behind the similarity could be that the Joint-z analysis results are likely to be biased by the strongest upper limit results, which are at $z = 9.1$.

The marginalised probability distributions of the disfavoured source parameter values at all redshifts are shown in Fig. 5. The related credible intervals are listed in Table 3 for both the Single-z and Joint-z cases. The credible intervals for $M_{\text{min}}$ and $M_{\text{min},X}$ for the three Single-z cases are qualitatively consistent with each other, as they disfavour models at the highest end of the prior range, i.e. $\log_{10}(M_{\text{min}})$ or $\log_{10}(M_{\text{min},X}) \ge 10.5 \text{ M}_\odot$. This suggests that models associated with a low number density of massive UV and X-ray emitting sources are disfavoured. The limits for $\zeta$, on the other hand, differ among the various redshifts.

<footnote n="17">¹⁷ A credible interval is a parameter range with a certain probability (e.g. 68%) of containing the true value of the parameter. Our constraints on the source and IGM parameters do not represent models ruled out with a high significance, but rather models with a high probability to be excluded. These interval limits are also highly dependent on the prior chosen on the source parameter ranges.</footnote>
<figure><img src="image_5.png" /><figcaption>**Fig. 4.** Posterior distribution of the disfavoured values of the source parameters for the $A_r = 0$ scenario. These are disfavoured by the LOFAR upper limits from Mertens et al. (2025) at $z = 9.1$. Red refers to the case when each redshift is considered separately (Single-z) at $z = 9.1$, while blue represents the joint MCMC analysis results (Joint-z). The contour levels in the two-dimensional contour plots refer to the $1\sigma$ and $2\sigma$ credible intervals of the models disfavoured by the LOFAR upper limit at this redshift. The diagonal panels represent the corresponding marginalised probability distributions of each parameter.</figcaption></figure>

For $z = 10.1$ models towards the higher end of our prior interval are disfavoured, but for $z = 8.3$ it is rather models at the lower end of our prior interval which are found to be less likely. The same trend is seen for the $f_X$ parameter. This can be understood from the difference in the number density of halos between these two redshifts. Scenarios more likely to be disfavoured require a high contrast, so the ionized or heated regions should not be too small and few or too large and many. With a smaller halo density at $z = 10.1$, relatively large efficiencies are needed to produce a large contrast population of ionized or heated regions, but at the larger halo density for $z = 8.3$ such high efficiencies would ionize or heat too much of the IGM, yielding a low contrast, and instead lower efficiencies are needed to create a large contrast distribution of ionized or heated regions.

Indeed, as also seen in our previous studies (Ghara et al. 2020, 2021, 2024) only a few IGM states can produce very large amplitudes for the large-scale power spectrum. One of these is a patchy reionization model with cold neutral regions. Such a model requires reionization by highly luminous UV emitting sources (large value of $\zeta$) with a low number density (large value of $M_{\min}$), combined with a small value of $f_X$ to keep the neutral regions cool. In such a scenario, the large-scale $\Delta^2$ is dominated by $x_{HII}$ fluctuations and reaches the largest amplitude for $\bar{x}_{HI} \sim 0.5$ (see e.g. Ghara et al. 2024). A second possibility is a patchy heating scenario, which requires rare and bright X-ray sources, implying large values for both $f_X$ and $M_{\min}$.

Figure 6 shows the posterior distributions of the disfavoured IGM parameter values at $z = 9.1$. The corresponding credible intervals are listed in Table 4. The disfavoured models have $\bar{x}_{HII} \lesssim 0.13$ (0.55), $T_K \lesssim 7.3$ (21) K, $f_{\text{heat}} \lesssim 0.18$ (0.57), $-276$ (-317) $\lesssim \delta T_b \lesssim -123$ (-61) mK, $R_{\text{peak}}^{\text{heat}} \lesssim 28$ (40) $h^{-1}$ Mpc

<figure><img src="image_6.png" /><figcaption>**Fig. 5.** Marginalised probability distribution of the disfavoured values of the source parameters for the $A_r = 0$ scenario. The corresponding models are disfavoured by the LOFAR upper limits from Mertens et al. (2025) at $z = 9.1$. The black, red, and green curves refer to the Single-z case at $z = 8.3$, 9.1 and 10.1, respectively, while the blue curves are for the Joint-z analysis case.</figcaption></figure>

<table><caption>**Table 3.** Constraints on the source parameters derived using the LOFAR upper limits from Mertens et al. (2025) for the $A_r = 0$ scenario.</caption><thead><tr><th>Source<br>Parameters</th><th>z = 8.3<br>(Single-z)</th><th>z = 9.1<br>(Single-z)</th><th>z = 10.1<br>(Single-z)</th><th>Joint-z</th></tr></thead><tbody><tr><td>log<sub>10</sub>(ζ)</td><td>[-3, 0.94]</td><td>X</td><td>[-0.53, 3]</td><td>X</td></tr><tr><td>log<sub>10</sub>(M<sub>min</sub>/M<sub>☉</sub>)</td><td>[10.45, 12]</td><td>[10.49, 12]</td><td>[10.23, 11.7]</td><td>[10.45, 12]</td></tr><tr><td>log<sub>10</sub>(M<sub>min,x</sub>/M<sub>☉</sub>)</td><td>[10.78, 12]</td><td>[10.76, 12]</td><td>[10.5, 11.7]</td><td>[10.7, 12]</td></tr><tr><td>log<sub>10</sub>(f<sub>X</sub>)</td><td>[-3, 1.54]</td><td>X</td><td>[0.6, 3]</td><td>X</td></tr></tbody></table>

**Notes.** The first column shows the list of source parameters, the following three constraints obtained from the Single-z analysis at three redshifts, and the last one those obtained with the Joint-z analysis. All these are disfavoured limits at 68% credible intervals level. The ‘X’ sign means that the credible intervals span the entire prior range of a parameter. The 95% disfavoured credible intervals span the entire prior ranges for all cases and thus are not shown here.

and $\Delta R_{\text{FWHM}}^{\text{heat}} \lesssim 46$ (81) $h^{-1}$ Mpc at 68% (95%) credible intervals for the Single-z analysis case. The PDF of $\bar{x}_{HII}$ indicates that a significant fraction of the disfavoured models are nearly fully neutral ones, where the large-scale power spectrum is mainly determined by fluctuations in $T_s$. As previously mentioned, the derived constraints on the IGM parameters might be significantly affected by the chosen priors on the source parameters (see Appendix A for details).

Table 4 also includes the results on the IGM parameters obtained from the MCMC analysis performed at the other two redshifts. The Single-z analysis gives fairly similar results at all redshifts, for both the 68% and 95% credible intervals levels. Indeed, the disfavoured IGM states are those which cause a high contrast in the $\delta T_b$ maps, producing large $\Delta^2$ values at scales $0.076 h \text{ Mpc}^{-1} \lesssim k \lesssim 0.133 h \text{ Mpc}^{-1}$. Unlike the Single-z case, in the Joint-z analysis the IGM parameter’s disfavoured credible
<figure><img src="image_7.png" /><figcaption>**Fig. 6.** Posterior distribution of the disfavoured values of the IGM parameters for the $A_r = 0$ scenario. These are disfavoured by the LOFAR upper limits from Mertens et al. (2025) at $z = 9.1$. Red refers to the case when $z = 9.1$ is considered separately (Single-z), while blue represents the joint redshift analysis obtained at $z = 9.1$ (Joint-z). The contour levels in the two-dimensional contour plots refer to the $1\sigma$ and $2\sigma$ credible intervals of the models disfavoured by the LOFAR upper limit at this redshift. The diagonal panels represent the corresponding marginalised probability distributions of each parameter.</figcaption></figure>

intervals show a consistent evolution with redshift. For example, the 95% credible intervals limits on $\bar{x}_{HII}$ are $\le 0.85$, $\le 0.46$ and $\le 0.17$ at $z = 8.3, 9.1$ and $10.1$, respectively. The main reason behind this evolution is that in this case, the disfavoured limits of $M_{min}$ and $M_{min,X}$ are basically the same at all redshifts, while the number density of dark matter halos contributing to ionization or heating increases towards lower redshift for a fixed $M_{min}$ and $M_{min,X}$. This causes a more efficient ionization and heating process towards the lower redshift for disfavoured models from the Joint-z analysis.

## 3.2. Analysis with excess radio background (Varying $A_r$)

Next, we present the results of the MCMC analysis for the 5D source parameter space, i.e. $\zeta$, $M_{min}$, $M_{min,X}$, $f_X$ and $A_r$, the same as above with the addition of the excess radio background. Table 2 lists the priors on these parameters.

Figure 7 shows the posterior distribution of the five parameters for the disfavoured models at $z = 9.1$. For the Single-z approach, the 68% disfavoured credible intervals limits are $\zeta \le 2.4$, $M_{min} \ge 1.7 \times 10^{10} M_\odot$, $M_{min,X} \ge 1.6 \times 10^{10} M_\odot$, $f_X \le 5.2$ and $A_r \ge 4.6$ (see Table 5). The 95% disfavoured credible intervals of the source parameters span the entire parameter ranges for all three Single-z and Joint-z analyses. The physical reasons behind the high amplitude of the large-scale power spectrum are the same as those discussed in the previous section. Additionally, a higher value of $A_r$ increases the signal strength and thus the power spectrum amplitude. We note that the effect of $A_r$ is model dependent. In fact, as for $A_r = 0$, also in this scenario the signal from fully ionized gas is zero because $\delta T_b = 0$, while

<table><caption>**Table 4.** Constraints on the IGM parameters derived using the LOFAR upper limits from Mertens et al. (2025) for the $A_r = 0$ scenario.</caption><thead><tr><th>IGM<br>Parameters<br>(MCMC<br>Model)</th><th>Credible<br>intervals<br>[%]</th><th>$z = 8.3$</th><th>$z = 9.1$</th><th>$z = 10.1$</th></tr></thead><tbody><tr><td>$\bar{x}_{HII}$<br>(Single-z)</td><td>68<br>95</td><td>[0, 0.15]<br>[0, 0.57]</td><td>[0, 0.13]<br>[0, 0.55]</td><td>[0, 0.2]<br>[0, 0.47]</td></tr><tr><td>$\bar{x}_{HII}$<br>(Joint-z)</td><td>68<br>95</td><td>[0, 0.24]<br>[0, 0.85]</td><td>[0, 0.06]<br>[0, 0.46]</td><td>[0, 0.01]<br>[0, 0.17]</td></tr><tr><td>$f_{\text{heat}}$<br>(Single-z)</td><td>68<br>95</td><td>[0, 0.21]<br>[0, 1]</td><td>[0, 0.18]<br>[0, 0.57]</td><td>[0, 0.22]<br>[0, 0.47]</td></tr><tr><td>$f_{\text{heat}}$<br>(Joint-z)</td><td>68<br>95</td><td>[0, 0.49]<br>[0, 0.9]</td><td>[0, 0.15]<br>[0, 0.54]</td><td>[0, 0.05]<br>[0, 0.21]</td></tr><tr><td>$\bar{T}_{\text{K}}$ (K)<br>(Single-z)</td><td>68<br>95</td><td>[1.8, 7.7]<br>[1.8, 50]</td><td>[2.2, 7.3]<br>[2.2, 21]</td><td>[2.6, 12.8]<br>[2.6, 22]</td></tr><tr><td>$\bar{T}_{\text{K}}$ (K)<br>(Joint-z)</td><td>68<br>95</td><td>[1.8, 28.5]<br>[1.8, 977]</td><td>[2.2, 7.2]<br>[2.2, 35.6]</td><td>[2.6, 5]<br>[2.6, 15]</td></tr><tr><td>$\delta T_b$ (mK)<br>(Single-z)</td><td>68<br>95</td><td>[-300, -124]<br>[-345, -23]</td><td>[-276, -123]<br>[-317, -61]</td><td>[-222, -113]<br>[-283, -80]</td></tr><tr><td>$\delta T_b$ (mK)<br>(Joint-z)</td><td>68<br>95</td><td>[-264, -4]<br>[-313, 10]</td><td>[-308, -141]<br>[-337, -36]</td><td>[-308, -220]<br>[-289, -116]</td></tr><tr><td>$R_{\text{heat,peak}}^{\text{Mpc}}$<br>(Single-z)</td><td>68<br>95</td><td>[0, 26.3]<br>[0, 56.4]</td><td>[0, 28.1]<br>[0, 39.5]</td><td>[0, 20.6]<br>[0, 35.5]</td></tr><tr><td>$R_{\text{heat,peak}}^{\text{Mpc}}$<br>(Joint-z)</td><td>68<br>95</td><td>[0, 68]<br>[0, 205]</td><td>[0, 25]<br>[0, 41]</td><td>[0, 14]<br>[0, 23]</td></tr><tr><td>$\Delta R_{\text{FWHM}}^{\text{Mpc}}$<br>(Single-z)</td><td>68<br>95</td><td>[6, 47]<br>[3, 118]</td><td>[9, 46]<br>[4, 81]</td><td>[11, 38]<br>[3.5, 79]</td></tr><tr><td>$\Delta R_{\text{FWHM}}^{\text{Mpc}}$<br>(Joint-z)</td><td>68<br>95</td><td>[0, 135]<br>[0, 331]</td><td>[0, 40]<br>[0, 85]</td><td>[0, 21]<br>[0, 33]</td></tr></tbody></table>

**Notes.** The first column shows the list of IGM parameters derived from both the Single-z and Joint-z analysis, the second refers to the credible intervals of the disfavoured models, while the last three to the constraints obtained from the analysis at three redshifts.

the strength of the signal from neutral or partially ionized gas increases for larger values of $A_r$. This causes a different enhancement in the power spectrum amplitude between a scenario with a large fraction of ionized regions (such as a patchy reionization scenario) and a neutral scenario (e.g. a patchy heating scenario with $\bar{x}_{HII} \sim 1$). This is one of the major reasons behind the differences in the posterior distribution of the source parameters compared to the previous case.

Table 5 shows the 68% credible intervals on the disfavoured values of the source parameters at all redshifts for the Single-z analysis, and at $z = 9.1$ for the Joint-z one. The credible intervals of the Single-z analysis are similar in trend for all redshifts, except for $M_{min,X}$ which remains unconstrained at $z = 10.1$. The disfavoured credible intervals on $\zeta$ and $f_X$ indicate smaller values towards lower redshifts. The 68% credible interval limits on disfavoured $A_r$ values for the Single-z analysis correspond to $T_{\gamma,\text{eff}} \ge 60, 55$ and $95$ K at $z = 8.3, 9.1$ and $10.1$, respectively, which are equivalent to an excess radio background larger than 140, 100, 207% of the CMB at 1.42 GHz. The Joint-z analysis posteriors for the source parameters show a trend similar to those from the Single-z cases, but closer to the $z = 9.1$ Single-z analysis results. The 68% credible intervals on the disfavoured $A_r$ values correspond to $T_{\gamma,\text{eff}} \ge 43.2$ K at $z = 9.1$ for the Joint-z analysis, meaning that an excess radio background higher than 57% of the CMB at 1.42 GHz is disfavoured. On the other hand, the 95% credible intervals on the disfavoured $A_r$ value span the entire prior range of $A_r$ for both the Single-z and Joint-z cases.
<figure><img src="image_8.png" /><figcaption>**Fig. 7.** Posterior distribution of the source parameters disfavoured values for the *Varying A<sub>r</sub>* scenario. The corresponding models are disfavoured by the LOFAR upper limits from Mertens et al. (2025) at *z* = 9.1. Red refers to the case when each redshift is considered separately (*Single-z*) at *z* = 9.1, while blue represents the joint MCMC analysis results (*Joint-z*). The contour levels in the two-dimensional contour plots refer to the 1σ and 2σ credible intervals of the disfavoured models. The diagonal panels represent the corresponding marginalised probability distributions of each parameter.</figcaption></figure>

<table><caption>**Table 5.** Constraints on the source parameters derived using the LOFAR upper limits from Mertens et al. (2025) for the *Varying A<sub>r</sub>* scenario.</caption><thead><tr><th>Source<br>Parameters</th><th>z = 8.3<br>(Single-z)</th><th>z = 9.1<br>(Single-z)</th><th>z = 10.1<br>(Single-z)</th><th>Joint-z</th></tr></thead><tbody><tr><td>log<sub>10</sub>(ζ)</td><td>[-3, 0.08]</td><td>[-3, 0.4]</td><td>[-3, 0.61]</td><td>[-3, 0.36]</td></tr><tr><td>log<sub>10</sub>(M<sub>min</sub>/M<sub>☉</sub>)</td><td>[10.25, 12]</td><td>[10.23, 12]</td><td>[10.2, 11.7]</td><td>[10.23, 12]</td></tr><tr><td>log<sub>10</sub>(M<sub>min,X</sub>/M<sub>☉</sub>)</td><td>[10.27, 12]</td><td>[10.2, 12]</td><td>X</td><td>[10.2, 12]</td></tr><tr><td>log<sub>10</sub>(f<sub>X</sub>)</td><td>[-3, 0.46]</td><td>[-3, 0.72]</td><td>[-3, 0.98]</td><td>[-3, 0.85]</td></tr><tr><td>log<sub>10</sub>(A<sub>r</sub>)</td><td>[0.9, 2.62]</td><td>[0.66, 2.62]</td><td>[0.86, 2.62]</td><td>[0.42, 2.62]</td></tr></tbody></table>

**Notes.** The first column shows the list of source parameters, the following three constraints obtained from the *Single-z* analysis at three redshifts, and the last one those obtained with the *Joint-z* analysis at *z* = 9.1. All these are disfavoured limits at 68% credible intervals level. The ‘X’ sign means that the credible intervals span the entire prior range of a parameter. The 95% disfavoured credible intervals span the entire prior ranges for all four cases.

Figure 8 shows the IGM parameters’ posterior distributions for the disfavoured models obtained from the MCMC analysis at *z* = 9.1. The 68% (95%) credible intervals for the *Single-z* analysis case at *z* = 9.1 are $\bar{x}_{HII} \lesssim 0.02$ (0.46), $\bar{T}_K \lesssim 4.4$ (44) K, $f_{\text{heat}} \lesssim 0.05$ (0.46), $1.3 \times 10^4$ (2.6 × 10⁴) $\gtrsim |\delta T_b| \gtrsim 158$ (91) mK, $R_{\text{peak}}^{\text{heat}} \lesssim 3$ (14) $h^{-1}$ Mpc and $\Delta R_{\text{FWHM}}^{\text{heat}} \lesssim 5$ (25) $h^{-1}$ Mpc. The ionization and the thermal states of the disfavoured models are similar to those in the scenario with $A_r = 0$, and correspond to extreme models with rare and large ionized or heated regions. We repeat the caveat that the constraints on the IGM parameters might be significantly impacted by the chosen priors on the source parameters (see Appendix A for details).

<figure><img src="image_9.png" /><figcaption>**Fig. 8.** Posterior distribution of the IGM parameters of the disfavoured models at *z* = 9.1 for the *Varying A<sub>r</sub>* scenario. Red refers to the case when each redshift is considered separately (*Single-z*), while blue represents the joint MCMC analysis obtained at *z* = 9.1 (*Joint-z*). The contour levels in the two-dimensional contour plots refer to the 1σ and 2σ credible intervals of the models disfavoured by the LOFAR upper limit from Mertens et al. (2025) at this redshift. The diagonal panels represent the corresponding marginalised probability distributions of each parameter.</figcaption></figure>

Table 6 shows the disfavoured limits for all three *Single-z* and *Joint-z* analyses at 68% and 95% credible intervals. In the *Single-z* cases, the disfavoured limits on IGM parameters for all redshifts are similar, indicating that the IGM ionization and thermal states producing high amplitudes of the large-scale power spectra are similar. As in the case with $A_r = 0$, we also see a consistent redshift evolution of the IGM parameters’ disfavoured limits in the *Joint-z* analysis case.

### 3.3. Analysis including results from other interferometers

Besides LOFAR, several other radio interferometers, such as GMRT, PAPER, MWA and HERA, have produced upper-limits on the 21 cm signal power spectrum in the range $z = 8-10$. In this section we consider all available upper limits to derive constraints on the IGM physical properties in this redshift range. Except for LOFAR and MWA, the other upper limits are obtained in a limited small-scale range with $k \gtrsim 0.3 h \, \text{Mpc}^{-1}$. Combining all available upper limits in a MCMC analysis needed to consider a broader $k$ range compared to the earlier cases.

Figure 9 shows those 2σ upper limits along with 1000 GRIZZLY power spectra randomly chosen from the $10^5$ simulations for $A_r = 0$. We note that at $z \sim 8$ only LOFAR (Mertens et al. 2025) and HERA (Abdurashidova et al. 2023) upper limits can rule out some EoR models at the 2σ level. In comparison to LOFAR, a significantly larger number of models is ruled out by the HERA best 2σ upper limit $\Delta^2(k = 0.34 h \, \text{Mpc}^{-1}, z = 7.9, 1128 \, \text{hours}) \approx$
<table><caption>**Table 6.** Constraints on the IGM parameters derived using the LOFAR upper limits from Mertens et al. (2025) for the Varying $A_r$ scenario.</caption><thead><tr><th>IGM<br>Parameters<br>(MCMC<br>Model)</th><th>Credible<br>intervals</th><th>z = 8.3</th><th>z = 9.1</th><th>z = 10.1</th></tr></thead><tbody><tr><td rowspan="2">x&#772;<sub>HII</sub><br>(Single-z)</td><td>68%</td><td>[0, 0.03]</td><td>[0, 0.02]</td><td>[0, 0.01]</td></tr><tr><td>95%</td><td>[0, 0.47]</td><td>[0, 0.46]</td><td>[0, 0.36]</td></tr><tr><td rowspan="2">x&#772;<sub>HII</sub><br>(Joint-z)</td><td>68%</td><td>[0, 0.06]</td><td>[0, 0.02]</td><td>[0, 0.004]</td></tr><tr><td>95%</td><td>[0, 0.8]</td><td>[0, 0.43]</td><td>[0, 0.16]</td></tr><tr><td rowspan="2">f<sub>heat</sub><br>(Single-z)</td><td>68%</td><td>[0, 0.07]</td><td>[0, 0.05]</td><td>[0, 0.03]</td></tr><tr><td>95%</td><td>[0, 0.49]</td><td>[0, 0.46]</td><td>[0, 0.34]</td></tr><tr><td rowspan="2">f<sub>heat</sub><br>(Joint-z)</td><td>68%</td><td>[0, 0.27]</td><td>[0, 0.07]</td><td>[0, 0.01]</td></tr><tr><td>95%</td><td>[0, 0.88]</td><td>[0, 0.53]</td><td>[0, 0.2]</td></tr><tr><td rowspan="2">T&#772;<sub>K</sub><br>(Single-z)</td><td>68%</td><td>[1.8, 5.0]</td><td>[2.2, 4.4]</td><td>[2.6, 4.8]</td></tr><tr><td>95%</td><td>[1.8, 6.7]</td><td>[2.2, 4.4]</td><td>[2.6, 4.1]</td></tr><tr><td rowspan="2">T&#772;<sub>K</sub><br>(Joint-z)</td><td>68%</td><td>[1.8, 16]</td><td>[2.2, 5.2]</td><td>[2.6, 4.5]</td></tr><tr><td>95%</td><td>[1.8, 977]</td><td>[2.2, 101]</td><td>[2.6, 186]</td></tr><tr><td rowspan="2">|δT<sub>bl</sub>| (mK)<br>(Single-z)</td><td>68%</td><td>[174, 11481]</td><td>[158, 13182]</td><td>[396, 11220]</td></tr><tr><td>95%</td><td>[93, 22908]</td><td>[91, 25703]</td><td>[144, 29174]</td></tr><tr><td rowspan="2">|δT<sub>bl</sub>| (mK)<br>(Joint-z)</td><td>68%</td><td>[25, 3715]</td><td>[135, 13301]</td><td>[177, 16218]</td></tr><tr><td>95%</td><td>[5, 19952]</td><td>[35, 28840]</td><td>[158, 28840]</td></tr><tr><td rowspan="2">R<sup>heat</sup><sub>peak</sub><br>(h<sup>-1</sup> Mpc)<br>(Single-z)</td><td>68%</td><td>[0, 3.0]</td><td>[0, 3.0]</td><td>[0, 3.0]</td></tr><tr><td>95%</td><td>[0, 14.4]</td><td>[0, 14]</td><td>[0, 11.2]</td></tr><tr><td rowspan="2">R<sup>heat</sup><sub>peak</sub><br>(h<sup>-1</sup> Mpc)<br>(Joint-z)</td><td>68%</td><td>[0, 6]</td><td>[0, 3.0]</td><td>[0, 3.0]</td></tr><tr><td>95%</td><td>[0, 105]</td><td>[0, 16.2]</td><td>[0, 8.5]</td></tr><tr><td rowspan="2">ΔR<sup>heat</sup><sub>FWHM</sub><br>(h<sup>-1</sup> Mpc)<br>(Single-z)</td><td>68%</td><td>[2, 7]</td><td>[0, 5]</td><td>[0, 4]</td></tr><tr><td>95%</td><td>[0, 27]</td><td>[0, 25]</td><td>[0, 18]</td></tr><tr><td rowspan="2">ΔR<sup>heat</sup><sub>FWHM</sub><br>(h<sup>-1</sup> Mpc)<br>(Joint-z)</td><td>68%</td><td>[0, 10]</td><td>[0, 6]</td><td>[0, 4]</td></tr><tr><td>95%</td><td>[0, 152]</td><td>[0, 30]</td><td>[0, 13.5]</td></tr></tbody></table>

**Notes.** The first column lists the IGM parameters derived from both the Single-z and Joint-z analysis, the second refers to the credible intervals of the disfavoured models, while the last three contain the constraints obtained from the analysis at three redshifts.

(21.4)$^2$ Mpc$^2$, including the scenario of a completely neutral and unheated IGM (dashed line). This is consistent with the findings of Abdurashidova et al. (2022). The best upper limit from MWA (Trott et al. 2020), $\Delta^2(k = 0.14 h \text{ Mpc}^{-1}, z = 7.8, 51 \text{ hours}) \approx (154)^2 \text{ Mpc}^2$, is still above the extreme EoR models power spectra by a factor of $\approx 2$. The best $2\sigma$ result from PAPER is $\Delta^2(k = 0.37 h \text{ Mpc}^{-1}, z = 8.4, 135 \text{ nights}) \approx (205)^2 \text{ Mpc}^2$, which is larger than the predicted power spectra by a factor of $\ge 4$.

At $z \sim 9$, the LOFAR upper limits from Mertens et al. (2025) are the best available results, and remain the only ones which rule out a set of extreme EoR scenarios. The best $2\sigma$ upper limit from GMRT (Paciga et al. 2013) and MWA (Trott et al. 2020)

<figure><img src="image_10.png" /><figcaption>**Fig. 9.** Set of 1000 power spectra randomly chosen from the $10^5$ simulated spectra for $A_r = 0$ scenario. Panels from top to bottom refer to $z = 8.3, 9.1$, and $10.1$, respectively. The dashed lines correspond to the power spectrum for a completely neutral IGM with no X-ray heating and with $T_S = T_K$ at those redshifts. The down arrow points refer to the different $2\sigma$ upper limits at those redshifts.</figcaption></figure>

at this redshift are $\Delta^2(k = 0.5 h \text{ Mpc}^{-1}, z = 8.6, 50 \text{ hours}) \approx (248)^2 \text{ Mpc}^2$ and $\Delta^2(k = 0.14 h \text{ Mpc}^{-1}, z = 8.8, 51 \text{ hours}) \approx (250)^2 \text{ Mpc}^2$, respectively, i.e. more than $\approx 6$ times larger than our extreme models. The best $2\sigma$ result from PAPER is $\Delta^2(k = 0.36 h \text{ Mpc}^{-1}, z = 8.7, 135 \text{ nights}) \approx (281)^2 \text{ Mpc}^2$, which is larger than the GRIZZLY power spectra by a factor of at least $\approx 8$. We note that the neutral and unheated IGM scenario at $z \sim 9$ (dashed curve) hasn't been ruled out by any of these upper limits.

Similarly to $z \sim 8$, also at $z \sim 10$ only LOFAR and HERA upper limits can rule out some extreme EoR models, as shown in the bottom panel of Fig. 9. The best $2\sigma$ upper limit result from HERA at this redshift is $\Delta^2(k = 0.36 h \text{ Mpc}^{-1}, z = 10.4, 1128 \text{ hours}) \approx (59)^2 \text{ Mpc}^2$, which is marginally able to rule out the neutral and unheated EoR scenario. Compared to $z \sim 8$, the number of models ruled out by HERA at $z \sim 10$ is smaller. The only other interferometer which has collected data at this redshift is PAPER, with a best $2\sigma$ result of $\Delta^2(k =$
<figure><img src="image_11.png" /></figure>

**Fig. 10.** Marginalised probability distribution of the disfavoured models' source parameters at $z \sim 8, 9$, and 10 for allUL case. The results are for the $A_r = 0$ scenario. The black, red, and green curves refer to the Single-z case at $z = 8.3, 9.1$ and 10.1, respectively, while the blue curves are for the Joint-z analysis case. The solid curves represent results derived from the LOFAR upper limits only, while the dotted curves refer to cases where we have used all available upper limits.

0.34 $h$ Mpc$^{-1}$, $z = 9.9$, 135 nights) $\approx (1923)^2$ mK$^2$, which is several orders of magnitude higher than our theoretical predictions.

Next, we performed an MCMC analysis on the disfavoured models for the $A_r = 0$ scenario, including all available upper limit results at $z \sim 8, 9$, and 10. In this case, we included one k-bin corresponding to the best upper limit results from each telescope in addition to the three LOFAR k-bins considered in our previous analysis. We denote this as the allUL scenario. The 1D posterior distribution of the disfavoured source parameters for allUL and $A_r = 0$ is shown (dotted lines) in Fig. 10 along with the LOFAR only case (solid). The 68% disfavoured credible intervals from the Single-z analysis at $z \sim 9$ for the allUL case are $M_{\text{min}} \ge 2.4 \times 10^{10}$ M$_\odot$ and $M_{\text{min},X} \ge 3.1 \times 10^{10}$ M$_\odot$, while the limits for $\zeta$ and $f_X$ span the entire parameter ranges. The 95% disfavoured credible intervals span the entire parameter ranges for both the Single-z and Joint-z analysis. For the Joint-z and allUL analysis case, the 68% disfavoured credible intervals at $z \sim 9$ are $\zeta \le 1.8$ and $M_{\text{min}} \ge 1.4 \times 10^{10}$ M$_\odot$, while the limits for $M_{\text{min},X}$ and $f_X$ span the entire parameter ranges, similarly to the Single-z analysis. The presence of small-scale power spectra upper limits in the MCMC analysis, causes a significant change in the marginalised probability distribution of the source parameters with respect to the LOFAR-only case. The minimum change is observed at $z \sim 9$, as the LOFAR upper limits at the large-scales are the best ones and thus dominate the MCMC analysis. Instead, at the other two redshifts the posteriors are mainly determined by the HERA upper limits.

In Fig. 11 we present the 1D marginalised probability distribution of the disfavoured models' IGM parameters at $z \sim 9$ for the allUL scenario with $A_r = 0$, which in the Single-

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

**Fig. 11.** Marginalised probability distribution of the IGM parameters of the disfavoured models at $z \sim 9$ for allUL and $A_r = 0$ scenario. The red and blue curves refer to the Single-z and Joint-z cases, respectively. The solid curves represent results derived from the LOFAR upper limits only, while the dotted refer to cases where we have used all available upper limits.

$z$ analysis case have $\bar{x}_{\text{HII}} \le 0.08$ (0.54), $\bar{T}_K \le 7$ (38) K, $-308$ ($-344$) $\le \bar{\delta T}_b \le -134$ (27) mK, $R_{\text{peak}}^{\text{heat}} \le 21$ (44) $h^{-1}$ Mpc and $\Delta R_{\text{FWHM}}^{\text{heat}} \le 35$ (104) $h^{-1}$ Mpc at 68% (95%) credible intervals, while $f_{\text{heat}}$ spans the entire range. With the Joint-z analysis we instead obtain $\bar{x}_{\text{HII}} \le 0.02$ (0.42), $\bar{T}_K \le 6$ (270) K, $-307$ ($-307$) $\le \bar{\delta T}_b \le -176$ (13) mK, $f_{\text{heat}} \le 0.12(X)$, $R_{\text{peak}}^{\text{heat}} \le 7$ (63) $h^{-1}$ Mpc and $\Delta R_{\text{FWHM}}^{\text{heat}} \le 20$ (151) $h^{-1}$ Mpc. While the constraints on $\bar{x}_{\text{HII}}$, $f_{\text{heat}}$, $\bar{T}_K$ look similar for the allUL and LOFAR-only cases, significant differences are seen for the other parameters. Indeed, due to the presence of the small-scale upper limits with a large error, more reionization scenarios with smaller emission/ionized regions are disfavoured in the allUL scenario.

Finally, Fig. 12 presents the 1D marginalised probability distribution of the disfavoured IGM parameters at $z \sim 9$ for the allUL and Varying $A_r$ scenario. The 68% (95%) credible intervals of the disfavoured models for the Single-z analysis are $\bar{x}_{\text{HII}} \le 0.02$ (0.46), $f_{\text{heat}} \le 0.06$ (0.5) $\bar{T}_K \le 5$ (55) K, 131 (62) $\le |\bar{\delta T}_b|$ $\le 3812$ (28184) mK, $R_{\text{peak}}^{\text{heat}} \le 5$ (17) $h^{-1}$ Mpc and $\Delta R_{\text{FWHM}}^{\text{heat}} \le 7$ (32) $h^{-1}$ Mpc. For the Joint-z analysis, these intervals instead are $\bar{x}_{\text{HII}} \le 0.02$ (0.42), $\bar{T}_K \le 6$ (371) K, 51 (13) $\le |\bar{\delta T}_b|$ $\le 2238$ (19952) mK, $f_{\text{heat}} \le 0.12(X)$, $R_{\text{peak}}^{\text{heat}} \le 5$ (35) $h^{-1}$ Mpc and $\Delta R_{\text{FWHM}}^{\text{heat}} \le 9$ (76) $h^{-1}$ Mpc. Table B2 presents the 68% and 95%
<figure><img src="image_13.png" /><figcaption>**Fig. 12.** Same as Fig. 11 but for the *Varying A*<sub>*r*</sub> scenario.</figcaption></figure>

credible intervals at all three redshifts for both the *Single-z* and *Joint-z* analysis. Contrary to Fig. 11, here we see that the constraints on all our IGM parameters look similar to the *allUL* and *LOFAR-only* case in the *Single-z* analysis. The difference on the posterior distributions for the *Joint-z* case as shown in Fig. 12 is also not as significant as found in Fig. 11. While the *LOFAR* *z* = 9.1 upper limits, being the best ones at this redshift, try to keep the MCMC results similar for the *allUL* and *LOFAR-only* *Single-z* case, the effect of source parameters’ priors in the *Varying A*<sub>*r*</sub> scenario also play an important role in producing similar posterior distributions.

We have also performed an MCMC analysis for the other two redshifts, the results of which can be found in Appendix B. This includes the constraints on the IGM parameters for the *allUL* case in Tables B.1 and B2. Figures B.1 and B.2 show the comparison between the *allUL* and *LOFAR-only* *A*<sub>*r*</sub> = 0 scenario for redshifts ~ 8 and 10, respectively. Clearly, the 1D posterior distributions are biased by the HERA upper limits, making the *allUL* and *LOFAR-only* PDFs different for most of the IGM parameters. This is also true for Figs. B.3 and B.4, which show the same comparison but for the *Varying A*<sub>*r*</sub> scenario for redshifts ~ 8 and 10, respectively.

**4. Summary and discussion**

Recently, Mertens et al. (2025) analysed 140 hours of *LOFAR* data for redshifts 8.3, 9.1, and 10.1, and produced new upper limits on the dimensionless spherically averaged power spectrum of the 21 cm signal (see Table 1). Here, we have used those upper

limits to infer the physical state of the IGM as well as the properties of the ionizing sources at those redshifts. As the upper limit results are still weak, we have derived constraints on the disfavoured rather than favoured models of reionization.

Our framework uses a database of hundreds of thousands of reionization models simulated with the GRIZZLY code, which employs a one-dimensional radiative transfer approach. The models are obtained by changing the following parameters associated with source properties: ionization efficiency ($\zeta$), minimum mass of the UV emitting halos ($M_{\text{min}}$), minimum mass of X-ray emitting halos ($M_{\text{min},X}$), X-ray heating efficiency ($f_X$), and the amplitude of a radio background in excess of the CMB ($A_r$). In addition to the power spectra of the 21 cm signal, various other quantities were derived from the reionization models to characterise the ionization and thermal state of the IGM, such as the volume-averaged ionization fraction ($\bar{x}_{\text{HII}}$), volume-averaged gas temperature of the partially ionized IGM ($\bar{T}_{\text{K}}$), mass-averaged brightness temperature ($\delta\bar{T}_{\text{b}}$), volume fraction of the heated region ($f_{\text{heat}}$), size of the heated regions at which the PDF of the sizes peak ($R_{\text{peak}}^{\text{heat}}$), and FWHM of the PDFs ($\Delta R_{\text{FWHM}}^{\text{heat}}$). We then used an MCMC framework to derive constraints on the source parameters and IGM-related quantities. We considered two different scenarios, one in which there is no excess radio background and another in which there is one characterised by a parameter $A_r$. We performed one MCMC analysis considering each redshift independently (referred to as the *Single-z* analysis case) and another in which we derived a joint likelihood for the three redshifts (*Joint-z* analysis). The constraints obtained on the disfavoured source and IGM parameters are based on the chosen uniform priors of the source parameters as presented in Table 2. We highlight that our source parameters are assumed to be redshift independent, and thus the results of the *Joint-z* analysis should be interpreted with caution. The main findings of the paper can be summarised as follows.

- More than 10% of the simulation database with $A_r = 0$ are well above $\Delta_{21}^2(k, z) - \Delta_{21,\text{err}}^2(k, z)$ in at least one k-bin at all three redshifts, indicating that these models have a high probability of being ruled out.

- The models disfavoured by both the *Single-z* and the *Joint-z* analysis for the scenario with $A_r = 0$ typically have a small number density of bright UV and X-ray emitting sources at all redshifts. More specifically, in the *Single-z* analysis, the 68% credible intervals on the disfavoured $M_{\text{min}}$ and $M_{\text{min},X}$ values at $z = 9.1$ are $\ge 3 \times 10^{10} M_{\odot}$ and $\ge 6 \times 10^{10} M_{\odot}$, respectively, while $\zeta$ and $f_X$ remain unconstrained (see Table 3 for the constraints at all redshifts). We find that the 95% disfavoured credible intervals of the source parameters span the entire prior ranges. The MCMC analysis results in the *Joint-z* case are significantly biased by the upper limits at $z = 9.1$, which are the strongest among the three redshifts considered here.

- For $A_r = 0$, the disfavoured limits on the IGM parameters at $z = 9.1$ are $\bar{x}_{\text{HII}} \le 0.13$ (0.55), $\bar{T}_{\text{K}} \le 7.3$ (21) K, $f_{\text{heat}} \le 0.18$ (0.57), $-276$ ($-317$) $\le \delta\bar{T}_{\text{b}} \le -123$ ($-61$) mK, $R_{\text{peak}}^{\text{heat}} \le 28$ (40) $h^{-1}$ Mpc, and $\Delta R_{\text{FWHM}}^{\text{heat}} \le 46$ (81) $h^{-1}$ Mpc at the 68% (95%) credible interval level (see Table 4 for details), indicating extreme reionization models with rare and large ionized or heated regions. The *Single-z* analysis gives similar credible intervals at all redshifts, showing that only a few IGM physical states can produce very high amplitudes of the large-scale power spectrum. As for the source parameters, the results of the *Joint-z* analysis for the IGM parameters are mainly biased by the $z = 9.1$ upper limits.
They show a consistent redshift evolution. For example, the 95% credible interval limits on $\bar{x}_{HII}$ are $\le 0.85$, $\le 0.46$, and $\le 0.17$ at $z = 8.3$, 9.1, and 10.1, respectively.

- When also considering an excess radio background, the models disfavoured by both the Single-z and the Joint-z analysis have a low number density of bright UV and X-ray emitting sources at all redshifts. More specifically, in the Single-z analysis, the 68% credible intervals on the source parameters at $z = 9.1$ are $\zeta \le 2.4$, $M_{\min} \ge 1.7 \times 10^{10} M_{\odot}$, $M_{\min,X} \ge 1.6 \times 10^{10} M_{\odot}$, $f_X \le 5.2$, and $A_r \ge 4.6$ (see Table 5 for the constraints at all redshifts). We find that the 95% disfavoured credible intervals span the entire prior ranges of the source parameters for both the Single-z and Joint-z analysis. We also find that the 68% credible intervals are $A_r \ge 7.9$ and 7.2 at $z = 8.3$ and 10.1, respectively, suggesting that an excess radio background that is at least 140%, 100%, and 207% of the CMB at 1.42 GHz is disfavoured at $z = 8.3$, 9.1, and 10.1, respectively. As for $A_r = 0$, the Joint-z analysis is biased by the upper limits at $z = 9.1$. In terms of the excess radio background, in this case we find that the 68% limits on the disfavoured models' $A_r$ at $z = 9.1$ correspond to an excess radio background that is $\ge 57\%$ of the CMB at 1.42 GHz.

- For such a radio background model, the disfavoured limits on the IGM parameters are $\bar{x}_{HII} \le 0.02$ (0.46), $\bar{T}_K \le 4.4$ (44) K, $f_{\text{heat}} \le 0.05$ (0.46), $1.3 \times 10^4$ ($2.6 \times 10^4$) $\ge |\delta\bar{T}_b|$ $\ge 158$ (91) mK, $R_{\text{peak}}^{\text{heat}} \le 3$ (14) $h^{-1}$ Mpc, and $\Delta R_{\text{FWHM}}^{\text{heat}} \le 5$ (25) $h^{-1}$ Mpc at $z = 9.1$ for the Single-z case at the 68% (95%) credible interval level (see Table 6 for all redshifts' constraints). This indicates that the disfavoured ionization and thermal states are similar to those in the scenario with $A_r = 0$ and correspond to extreme reionization models with rare and large ionized or heated regions. Similarly to the $A_r = 0$ scenario, we also find a consistent redshift evolution of the IGM parameters' disfavoured limits in the Joint-z analysis for the Varying $A_r$ case.

- The inclusion of upper limits from other radio interferometric observations (allUL case) in our study significantly increases the disfavoured parameter space of both the source and IGM parameters. We find that the HERA upper limits dominate the Bayesian analysis at $z \sim 8$ and 10, while the LOFAR ones prevail at $z \sim 9$. The HERA upper limits rule out a neutral and unheated IGM at $z \sim 8$ at a $2\sigma$ level. We find that the constraints on $\bar{x}_{HII}$, $f_{\text{heat}}$, and $\bar{T}_K$ are similar for the allUL and LOFAR-only cases, while those on the other parameters are significantly different.

The constraints derived in this paper improve on those obtained in Ghara et al. (2020) using the earlier LOFAR upper limits at $z = 9.1$ from Mertens et al. (2020). Indeed, in the absence of a radio background in excess of the CMB, the 95% credible intervals on the disfavoured models’ IGM parameters obtained in Ghara et al. (2020) were $f_{\text{heat}} \le 0.34$, $\bar{T}_K \le 160$ K, $R_{\text{peak}}^{\text{heat}} \le 70$ $h^{-1}$ Mpc and $\Delta R_{\text{FWHM}}^{\text{heat}} \le 110$ $h^{-1}$ Mpc, while we found two regimes of disfavoured credible intervals for $\bar{x}_{HII}$, i.e. $\bar{x}_{HII} \le 0.08$ and $0.45 \le \bar{x}_{HII} \le 0.62$. We note that there are some differences between the analysis done here and in Ghara et al. (2020). For example, the latter fixed $M_{\min} = 10^9 M_{\odot}$ while varying $\zeta$, $M_{\min,X}$ and $f_X$ with different prior ranges, and they did not include redshift space distortions in their simulated power spectra, resulting in under-predictions of the large-scale power spectrum for several scenarios. The improvement in the present version of our interpretation framework, as well as in the observed upper limits, has enabled us to rule out more reionization models and produce

better posterior distributions for both the disfavoured source and IGM parameters.

Previously, Mondal et al. (2020) constrained the excess radio background as 0.1–9.6% of the CMB at 1.42 GHz using the LOFAR upper limits at $z = 9.1$ from Mertens et al. (2020). We note that a direct comparison to these constraints is not possible, as they were obtained by marginalizing models that are allowed under those upper limits, while our results are derived from the marginalisation of the excluded models. In Ghara et al. (2021), an excess radio background of at least 25% of the CMB at 1.42 GHz was excluded using MWA upper limits in the range $z = 6.5–8.7$ from Trott et al. (2020) with 68% credible intervals on $A_r$.

The MCMC results presented in this study are dependent on the prior ranges of the source parameters. This is true not only for the disfavoured models but also for the allowed ones. Since the volume in parameter space is set by the priors, their impact can be very large, particularly if the interested scale is unknown and a flat prior in log space is employed. Therefore, the results strongly depend on the lower and upper limits set on the model parameters. In this work, we find that the constraints on $\bar{x}_{HII}$, $f_{\text{heat}}$, and $\bar{T}_K$ are significantly affected by the chosen priors on the source parameters (see Appendix A). To improve on this, observational data on high-$z$ galactic properties, such as the luminosity function from the James Webb Space Telescope (JWST; e.g. Donnan et al. 2022; Castellano et al. 2023), can be used to increase our understanding of the prior ranges of the source parameters and thus to improve the constraints on the source and IGM parameters of the excluded models.

The upper limits on the 21 cm signal power spectrum used in this study contain excess noise and thus are conservative in nature. This is also applicable to the previous LOFAR upper limits from Patil et al. (2017), Mertens et al. (2020), and Ceccotti et al. (2025). The source of this excess noise is still uncertain (see e.g. Gan et al. 2022). Recently, Acharya et al. (2024c) considered applying a bias correction for excess noise as well, and they reported a $2\sigma$ upper limit of $\Delta^2(k = 0.075 h \text{ Mpc}^{-1}, z = 9.1) \approx (25)^2 \text{ mK}^2$ based on 141 hours of LOFAR data used in Mertens et al. (2020). In Appendix C, we perform an inference study on the Acharya et al. (2024c) upper limits and find that the 95% credible intervals of the IGM parameters of the allowed models are $\bar{x}_{HII} \le 0.75$, $27 \le \bar{T}_K \le 7 \times 10^3$ K, $f_{\text{heat}} \ge 0.62$, $4 \le \delta\bar{T}_b \le 44$ mK, $3 \le R_{\text{peak}}^{\text{heat}} \le 160$ $h^{-1}$ Mpc, and $4 \le \Delta R_{\text{FWHM}}^{\text{heat}} \le 316$ $h^{-1}$ Mpc for the Varying $A_r$ and Single-z scenario at $z = 9.1$. As in this scenario, the upper limits are stronger, and the results are on the allowed rather than the excluded models, so they cannot be directly compared to the results presented in the rest of the paper. We also caution the reader that the robustness of the Acharya et al. (2024c) upper limits needs further investigation, as also discussed in the original paper. For this reason, we have focused our attention on the constraints obtained from the weaker upper limits of Mertens et al. (2025).

We also note that in principle, an inference study such as the one discussed here should be included in the data analysis pipeline together with the foreground mitigation rather than as a separate step. This, however, is beyond the scope of this paper but will be considered in future studies.

It should be noted that the current LOFAR upper limits on $\Delta^2$ are rather large, and thus the constraints obtained on the disfavoured IGM and source parameters are not yet very strong. However, this study illustrates the potential of the 21 cm signal observation as a probe of the EoR, as even before an actual detection of the signal becomes available, more stringent upper limits
on $\Delta^2$ would provide stronger constraints on the IGM parameters during the first billion years of the Universe. Additionally, other probes of this era such as the CMB, the global 21 cm signal, and high-$z$ sources (e.g. galaxies, quasars, and fast radio bursts) will strengthen the understanding of the IGM during this crucial period in the Universe's history. In particular, the detection of many faint galaxies at $z \ge 6$ with JWST, complemented by observations with the Atacama Large Millimetre Array and the future European Extremely Large Telescope, will improve our understanding of the properties of early galaxies. This information will be fundamental to improving the prediction of the 21 cm signal through the calibration of the source model with these observations (see e.g. the GRIZZLY extension POLAR; Ma et al. 2023) and will thus provide a much deeper understanding of the EoR sources and IGM properties.

**Acknowledgements.** RG acknowledges support from SERB, DST Ramanujan Fellowship no. RJF2022/000141. We acknowledge that the N-body simulation used in this paper has been achieved using the PRACE Research Infrastructure resources Curie based at the Très Grand Centre de Calcul (TGCC) operated by CEA near Paris, France and Marenström based in the Barcelona Supercomputing Center, Spain (PRACE under PRACE4LOFAR grants 2012061089 and 2014102339 as well as under the Multi-scale Reionization grants 2014102281 and 2015122822). SKG is supported by NWO grant number OCENW.M.22.307. LVEK, SAB, KC, SG, CH and SM acknowledge 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"). FGM acknowledges the support of the PSL Fellowship. EC acknowledge support from the Centre for Data Science and Systems Complexity (DSSC), Faculty of Science and Engineering at the University of Groningen, and the Ministry of Universities and Research (MUR) through the PRIN project 'Optimal inference from radio images of the epoch of reionization'. SZ acknowledge support by the Israel Science Foundation (grant no. 1388/24). GM's research is supported by the Swedish Research Council project grant 2020-04691_VR. MC's research is supported by the Center for Fundamental Physics of the Universe (CFPU) Postdoctoral fellowship, Brown University, Providence, RI.

## References

* Abdurashidova, Z., Aguirre, J. E., Alexander, P., et al. 2022, ApJ, 924, 51

* Abdurashidova, T. H. C. Z., Adams, T., Aguirre, J. E., et al. 2023, ApJ, 945, 124

* Acharya, S. K., Cyr, B., & Chluba, J. 2023, MNRAS, 523, 1908

* Acharya, A., Mertens, F., Ciardi, B., et al. 2024a, MNRAS, 527, 7835

* Acharya, A., Ma, Q. b., Giri, S. K., et al. 2024b, MNRAS, submitted [arXiv:2410.11620]

* Acharya, A., Mertens, F., Ciardi, B., et al. 2024c, MNRAS, 534, L30

* Ali, Z. S., Parsons, A. R., Zheng, H., et al. 2015, ApJ, 809, 61

* Bañados, E., Venemans, B. P., Mazzucchelli, C., et al. 2018, Nature, 553, 473

* Barkana, R. 2018, Nature, 555, 71

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

* Barkana, R., & Loeb, A. 2005, ApJ, 626, 1

* Barry, N., Hazelton, B., Sullivan, I., Morales, M. F., & Pober, J. C. 2016, MNRAS, 461, 3135

* Barry, N., Wilensky, M., Trott, C. M., et al. 2019, ApJ, 884, 1

* Beardsley, A. P., Hazelton, B. J., Sullivan, I. S., et al. 2016, ApJ, 833, 102

* Becker, G. D., Bolton, J. S., Madau, P., et al. 2015, MNRAS, 447, 3402

* Becker, G. D., Davies, F. B., Furlanetto, S. R., et al. 2018, ApJ, 863, 92

* Behroozi, P. S., & Silk, J. 2015, ApJ, 799, 32

* Bera, A., Ghara, R., Chatterjee, A., Datta, K. K., & Samui, S. 2023, JApA, 44, 10

* Bernardi, G., Zwart, J. T. L., Price, D., et al. 2016, MNRAS, 461, 2847

* Bonaldi, A., & Brown, M. L. 2015, MNRAS, 447, 1973

* Bowman, J. D., & Rogers, A. E. E. 2010, Nature, 468, 796

* Bowman, J. D., Rogers, A. E. E., Monsalve, R. A., Mozdzen, T. J., & Mahesh, N. 2018, Nature, 555, 67

* Bradley, R. F., Tauscher, K., Rapetti, D., & Burns, J. O. 2019, ApJ, 874, 153

* Castellano, M., Fontana, A., Treu, T., et al. 2023, ApJ, 948, L14

* Ceccotti, E., Offringa, A. R., Koopmans, L. V. E., et al. 2025, A&A, 696, A56

* Chapman, E., Zaroubi, S., Abdalla, F. B., et al. 2016, MNRAS, 458, 2928

* Chatterjee, A., Dayal, P., Choudhury, T. R., & Hutter, A. 2019, MNRAS, 487, 3560

* Chatterjee, A., Dayal, P., Choudhury, T. R., & Schneider, R. 2020, MNRAS, 496, 1445

* Cheng, C., Parsons, A. R., Kolopanis, M., et al. 2018, ApJ, 868, 26

* Choudhury, M., Ghara, R., Zaroubi, S., et al. 2024, arXiv e-prints [arXiv:2407.03523]

* Cohen, A., Fialkov, A., Barkana, R., & Monsalve, R. A. 2020, MNRAS, 495, 4845

* Cyr, B., Chluba, J., & Acharya, S. K. 2024, Phys. Rev. D, 109, L121301

* Datta, K. K., Bharadwaj, S., & Choudhury, T. R. 2007, MNRAS, 382, 809

* Datta, A., Bowman, J. D., & Carilli, C. L. 2010, ApJ, 724, 526

* DeBoer, D. R., Parsons, A. R., Aguirre, J. E., et al. 2017, PASP, 129, 045001

* de Lera Acedo, E., de Villiers, D. I. L., Razavi-Ghods, N., et al. 2022, Nat. Astron., 6, 984

* Donnan, C. T., McLeod, D. J., Dunlop, J. S., et al. 2022, MNRAS, 518, 6011

* Dowell, J., & Taylor, G. B. 2018, ApJ, 858, L9

* Draine, B. T., & Miralda-Escudé, J. 2018, ApJ, 858, L10

* Eastwood, M. W., Anderson, M. M., Monroe, R. M., et al. 2019, AJ, 158, 84

* Edler, H. W., de Gasperin, F., & Rafferty, D. 2021, A&A, 652, A37

* Eide, M. B., Ciardi, B., Graziani, L., et al. 2020, MNRAS, 498, 6083

* Eilers, A.-C., Davies, F. B., & Hennawi, J. F. 2018, ApJ, 864, 53

* Ewall-Wice, A., Hewitt, J., Mesinger, A., et al. 2016, MNRAS, 458, 2710

* Ewall-Wice, A., Dillon, J. S., Liu, A., & Hewitt, J. 2017, MNRAS, 470, 1849

* Fan, X., Strauss, M. A., Becker, R. H., et al. 2006, AJ, 132, 117

* Fialkov, A., & Barkana, R. 2019, MNRAS, 486, 1763

* Fialkov, A., Barkana, R., & Cohen, A. 2018, Phys. Rev. Lett., 121, 011101

* Fixsen, D. J., Kogut, A., Levin, S., et al. 2011, ApJ, 734, 5

* Furlanetto, S. R., Oh, S. P., & Briggs, F. H. 2006, Phys. Rep., 433, 181

* Gan, H., Koopmans, L. V. E., Mertens, F. G., et al. 2022, A&A, 663, A9

* Gan, H., Mertens, F. G., Koopmans, L. V. E., et al. 2023, A&A, 669, A20

* Gao, L. Y., Koopmans, L. V. E., Mertens, F. G., et al. 2024, arXiv e-prints [arXiv:2412.16853]

* Gehlot, B. K., Mertens, F. G., Koopmans, L. V. E., et al. 2019, MNRAS, 488, 4271

* Ghara, R., & Mellema, G. 2020, MNRAS, 492, 634

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

* Ghara, R., Datta, K. K., & Choudhury, T. R. 2015b, MNRAS, 453, 3143

* Ghara, R., Choudhury, T. R., & Datta, K. K. 2016, MNRAS, 460, 827

* Ghara, R., Choudhury, T. R., Datta, K. K., & Choudhuri, S. 2017, MNRAS, 464, 2234

* Ghara, R., Mellema, G., Giri, S. K., et al. 2018, MNRAS, 476, 1741

* Ghara, R., Giri, S. K., Mellema, G., et al. 2020, MNRAS, 493, 4728

* Ghara, R., Giri, S. K., Ciardi, B., Mellema, G., & Zaroubi, S. 2021, MNRAS, 503, 4551

* Ghara, R., Mellema, G., & Zaroubi, S. 2022, JCAP, 2022, 055

* Ghara, R., Shaw, A. K., Zaroubi, S., et al. 2024, A&A, 687, A252

* Ghosh, A., Prasad, J., Bharadwaj, S., Ali, S. S., & Chengalur, J. N. 2012, MNRAS, 426, 3295

* 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., et al. 2019b, JCAP, 2019, 058

* Greig, B., Mesinger, A., Koopmans, L. V. E., et al. 2021, MNRAS, 501, 1

* Harker, G., Zaroubi, S., Bernardi, G., et al. 2009, MNRAS, 397, 1138

* Harnois-Déraps, J., Pen, U.-L., Iliev, I. T., et al. 2013, MNRAS, 436, 540

* Hills, R., Kulkarni, G., Meerburg, P. D., & Puchwein, E. 2018, Nature, 564, E32

* Hinshaw, G., Larson, D., Komatsu, E., et al. 2013, ApJS, 208, 19

* Hothi, I., Chapman, E., Pritchard, J. R., et al. 2021, MNRAS, 500, 2264

* Hu, E. M., Cowie, L. L., Barger, A. J., et al. 2010, ApJ, 725, 394

* Islam, N., Ghara, R., Paul, B., Choudhury, T. R., & Nath, B. B. 2019, MNRAS, 487, 2785

* Jelić, V., Zaroubi, S., Labropoulos, P., et al. 2010, MNRAS, 409, 1647

* Kern, N. S., Parsons, A. R., Dillon, J. S., et al. 2019, ApJ, 884, 105

* Kern, N. S., Parsons, A. R., Dillon, J. S., et al. 2020, ApJ, 888, 70

* Kolopanis, M., Jacobs, D. C., Cheng, C., et al. 2019, ApJ, 883, 133

* Koopmans, L., Pritchard, J., Mellema, G., et al. 2015, Proceedings of Advancing Astrophysics with the Square Kilometre Array - PoS(AASKA14), AASKA14 (Sissa Medialab)

* Krause, F., Thomas, R. M., Zaroubi, S., & Abdalla, F. B. 2018, New Astron., 64, 9

* Lewis, A., & Bridle, S. 2002, Phys. Rev. D, 66, 103511

* Liu, A., Parsons, A. R., & Trott, C. M. 2014, Phys. Rev. D, 90, 023019

* Ma, Q.-B., Ghara, R., Ciardi, B., et al. 2023, MNRAS, 522, 3284

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

* Mao, Y., Shapiro, P. R., Mellema, G., et al. 2012, MNRAS, 422, 926

* Mellema, G., Koopmans, L. V. E., Abdalla, F. A., et al. 2013, Exp. Astron., 36, 235

* Mellema, G., Koopmans, L., Shukla, H., et al. 2015, Advancing Astrophysics with the Square Kilometre Array (AASKA14), 10

* Mertens, F. G., Ghosh, A., & Koopmans, L. V. E. 2018, MNRAS, 478, 3640

* Mertens, F. G., Mevius, M., Koopmans, L. V. E., et al. 2020, MNRAS, 493, 1662

* Mertens, F. G., Bobin, J., & Carucci, I. P. 2024, MNRAS, 527, 3517
Mertens, F. G., Mevius, M., Koopmans, L. V. E., et al. 2025, A&A, 698, A186

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

Mevius, M., van der Tol, S., Pandey, V. N., et al. 2016, Radio Sci., 51, 927

Mevius, M., Mertens, F., Koopmans, L. V. E., et al. 2022, MNRAS, 509, 3693

Mineo, S., Gilfanov, M., & Sunyaev, R. 2012, MNRAS, 419, 2095

Mirocha, J., Muñoz, J. B., Furlanetto, S. R., Liu, A., & Mesinger, A. 2022, MNRAS, 514, 2010

Mondal, R., Fialkov, A., Fling, C., et al. 2020, MNRAS, 498, 4178

Monsalve, R. A., Rogers, A. E. E., Bowman, J. D., & Mozdzen, T. J. 2017, ApJ, 835, 49

Monsalve, R. A., Fialkov, A., Bowman, J. D., et al. 2019, ApJ, 875, 67

Monsalve, R. A., Altamirano, C., Bidula, V., et al. 2024, MNRAS, 530, 4125

Morales, A. M., Mason, C. A., Bruton, S., et al. 2021, ApJ, 919, 120

Mortlock, D. J., Warren, S. J., Venemans, B. P., et al. 2011, Nature, 474, 616

Muñoz, J. B., & Loeb, A. 2018, Nature, 557, 684

Munshi, S., Mertens, F. G., Koopmans, L. V. E., et al. 2024, A&A, 681, A62

Munshi, S., Mertens, F. G., Koopmans, L. V. E., et al. 2025, A&A, 697, A203

Nebrin, O., Ghara, R., & Mellema, G. 2019, JCAP, 2019, 051

Offringa, A. R., Mertens, F., & Koopmans, L. V. E. 2019, MNRAS, 484, 2866

Paciga, G., Albert, J. G., Bandura, K., et al. 2013, MNRAS, 433, 639

Park, J., Mesinger, A., Greig, B., & Gillet, N. 2019, MNRAS, 484, 933

Parsons, A. R., Liu, A., Aguirre, J. E., et al. 2014, ApJ, 788, 106

Patil, A. H., Yatawatta, S., Zaroubi, S., et al. 2016, MNRAS, 463, 4317

Patil, A. H., Yatawatta, S., Koopmans, L. V. E., et al. 2017, ApJ, 838, 65

Patra, N., Subrahmanyan, R., Sethi, S., Udaya Shankar, N., & Raghunathan, A. 2015, ApJ, 801, 138

Philip, L., Abdurashidova, Z., Chiang, H. C., et al. 2018, J. Astron. Instrum., submitted [arXiv:1806.09531]

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

Price, D. C., Greenhill, L. J., Fialkov, A., et al. 2018, MNRAS, 478, 4193

Pritchard, J. R., & Loeb, A. 2012, Rep. Progr. Phys., 75, 086901

Reis, I., Fialkov, A., & Barkana, R. 2020, MNRAS, 499, 5993

Ross, H. E., Dixon, K. L., Ghara, R., Iliev, I. T., & Mellema, G. 2019, MNRAS, 487, 1101

Ross, H. E., Giri, S. K., Mellema, G., et al. 2021, MNRAS, 506, 3717

Schaeffer, T., Giri, S. K., & Schneider, A. 2024, Phys. Rev. D, 110, 023543

Schneider, A., Giri, S. K., & Mirocha, J. 2021, Phys. Rev. D, 103, 083025

Shapiro, P. R., Giroux, M. L., & Babul, A. 1994, ApJ, 427, 25

Singh, S., & Subrahmanyan, R. 2019, ApJ, 880, 26

Singh, S., Subrahmanyan, R., Udaya Shankar, N., et al. 2017, ApJ, 845, L12

Singh, S., Jishnu, N. T., Subrahmanyan, R., et al. 2022, Nat. Astron., 6, 607

Sobacchi, E., & Mesinger, A. 2013, MNRAS, 432, L51

Sokolowski, M., Tremblay, S. E., Wayth, R. B., et al. 2015, PASA, 32, e004

Spinelli, M., Bernardi, G., & Santos, M. G. 2018, MNRAS, 479, 275

Springel, V., & Hernquist, L. 2003, MNRAS, 339, 312

Subrahmanyan, R., & Cowsik, R. 2013, ApJ, 776, 42

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

Tang, M., Stark, D. P., Chen, Z., et al. 2023, MNRAS, 526, 1657

Thomas, R. M., & Zaroubi, S. 2008, MNRAS, 384, 1080

Thomas, R. M., & Zaroubi, S. 2011, MNRAS, 410, 1377

Thomas, R. M., Zaroubi, S., Ciardi, B., et al. 2009, MNRAS, 393, 32

Tingay, S. J., Goeke, R., Bowman, J. D., et al. 2013, PASA, 30, 7

Trott, C. M., Jordan, C. H., Midgley, S., et al. 2020, MNRAS, 493, 4711

van Haarlem, M. P., Wise, M. W., Gunst, A. W., et al. 2013, A&A, 556, A2

Venemans, B. P., Bañados, E., Decarli, R., et al. 2015, ApJ, 801, L11

Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nat. Methods, 17, 261

Voytek, T. C., Natarajan, A., Jáuregui García, J. M., Peterson, J. B., & López-Cruz, O. 2014, ApJ, 782, L9

Watson, W. A., Iliev, I. T., D'Aloisio, A., et al. 2013, MNRAS, 433, 1230

Weiser, A., & Zarantonello, S. E. 1988, Math. Comput., 50, 189

Wilensky, M. J., Morales, M. F., Hazelton, B. J., et al. 2019, PASP, 131, 114507

Witten, C., Laporte, N., Martin-Alvarez, S., et al. 2024, Nat. Astron., 8, 384

Yatawatta, S., & Banerjee, P. 2011, 2011 XXXth URSI General Assembly and Scientific Symposium (IEEE Piscataway)

Zaroubi, S. 2013, in The First Galaxies, eds. T. Wiklind, B. Mobasher, & V. Bromm, Astrophysics and Space Science Library, 396, 45
<figure><img src="image_14.png" /><figcaption>**Fig. A.1.** Marginalised probability distribution of the IGM parameters of the disfavoured models for a constant exclusion likelihood. The red and blue curves refer to the *Single-z* and *Joint-z* cases at *z* = 9.1, respectively. The solid curves represent results from the main text, while the dashed refers to cases where we have used $\mathcal{L}_{\text{ex,single-z}}(\theta, z = 9.1) = 0.5$ and $\mathcal{L}_{\text{ex,joint-z}}(\theta, z = 9.1) = 0.5$.</figcaption></figure>

## Appendix A: Effect of the source parameter priors on the IGM parameters constraints

In the study described in the main text we use source parameters as input for the MCMC analysis. After this the list of source parameter values is used to produce a corresponding list of IGM parameters, which is then employed to constrain the IGM parameters themselves. This method might introduce non-uniform posterior distributions of the IGM parameters even when the LOFAR upper limits are not considered in the analysis. In this section, we aim to investigate the effect of the source parameters' priors on the constraints on the IGM parameters, by adopting the following approach. We fix to 0.5 the exclusion likelihood shown in Eqs. (3) and (4), independently from the LOFAR upper limits, but still use the priors on the maximum value of $\bar{x}_{\text{HII}}$, and we denote this as the *no-LOFAR* scenario. In this case, the posterior distribution of the source parameters is expected to be independent of the choice of the constant likelihood, while the posterior distribution of the disfavoured IGM

parameter values is entirely determined by the chosen priors on the source parameters.

We consider the four source parameters for the $A_r = 0$ scenario from Sect. 3.1, but with a constant likelihood. The dashed curves in Fig. A.1 refer to the resulting marginalised probability distribution of the IGM parameters of the disfavoured models, for both the *Single-z* and *Joint-z* analysis. The distribution shows the effect of the source priors on the posteriors for the IGM parameters. The disfavoured models have $\bar{x}_{\text{HII}} \lesssim 0.024$ (0.56), $\bar{T}_K \lesssim 7.2$ (618) K, no constraints on $f_{\text{heat}}$, $-315 \lesssim \bar{\delta T}_b \lesssim -129$ (18.1) mK, $R_{\text{peak}}^{\text{heat}} \lesssim 10.8$ (83.1) $h^{-1}$ Mpc and $\Delta R_{\text{FWHM}}^{\text{heat}} \lesssim 87$ (207) $h^{-1}$ Mpc at 68% (95%) credible intervals for the *Single-z* no-LOFAR analysis case. For the *Joint-z* no-LOFAR analysis case, the disfavoured models' 68% (95%) credible intervals are $\bar{x}_{\text{HII}} \lesssim 0.019$ (0.43), $\bar{T}_K \lesssim 6$ (562) K, no constraints on $f_{\text{heat}}$, $-315 \lesssim \bar{\delta T}_b \lesssim -149$ (18.1) mK, $R_{\text{peak}}^{\text{heat}} \lesssim 7$ (36) $h^{-1}$ Mpc and $\Delta R_{\text{FWHM}}^{\text{heat}} \lesssim 23$ (93) $h^{-1}$ Mpc. For comparison, the solid curves show the probability distribution of the IGM parameters of the disfavoured models obtained in Section 3.1. Any differences between the solid and dashed curves are due to the effect of using the LOFAR upper limits. Such differences are clearly visible for $\bar{\delta T}_b$, $R_{\text{peak}}^{\text{heat}}$ and $\Delta R_{\text{FWHM}}^{\text{heat}}$, they are less prominent for $\bar{x}_{\text{HII}}$, $f_{\text{heat}}$ and $\bar{T}_K$. This suggests that the LOFAR upper limits provide very weak constraints on the latter three IGM parameters. The fact that the posterior for $\bar{\delta T}_b$ is affected by the LOFAR upper limits but the other averages ($\bar{x}_{\text{HII}}$, $f_{\text{heat}}$ and $\bar{T}_K$) are not, suggests that the correlations between $\delta$, $x_{\text{HII}}$ and $T_K$ play important roles in determining the 21 cm signal power spectra (see e.g. Schaeffer et al. 2024).

## Appendix B: Inference including all available upper limit results at $z \sim 8$ and 10

While the results from the MCMC analysis which includes all available upper limits (*allUL* case) at $z \sim 9$ have been discussed in Sect. 3.3, here we present the results at $z \sim 8$ and 10. The 68% disfavoured credible intervals of the source parameters for the $A_r = 0$ and *Single-z* scenario at $z \sim 8$ are $\zeta \lesssim 1.6$ and $M_{\text{min}} \gtrsim 1.4 \times 10^{10}$ M$_\odot$, while the limits for $M_{\text{min},X}$ and $f_X$ span the entire parameter ranges (see Fig. 10). The same intervals for $z \sim 10$ are $\zeta \lesssim 5$, $M_{\text{min}} \gtrsim 1.7 \times 10^{10}$ M$_\odot$, $M_{\text{min},X} \gtrsim 1.6 \times 10^{10}$ M$_\odot$, and $f_X \lesssim 4$. In this case, the source parameters' 95% disfavoured credible intervals span the entire parameter ranges for both the redshifts. The 1D marginalised probability distribution of the source parameters at $z \sim 8$ and 10 for the *allUL* case are significantly different from the LOFAR-only case, as HERA upper limits dominate the MCMC analysis at these two redshifts.

The 1D marginalised probability distribution of the disfavoured IGM parameters for the $A_r = 0$ and *allUL* case at $z \sim 8$ are shown in Fig. B.1 while Fig. B.2 shows the same but for $z \sim 10$. The disfavoured models at $z \sim 8$ in the *allUL*, *Single-z* and $A_r = 0$ scenario have $\bar{x}_{\text{HII}} \lesssim 0.05$ (0.78), $\bar{T}_K \lesssim 10$ (2180) K, $-352$ ($-375$) $\lesssim \bar{\delta T}_b \lesssim 3.2$ (42) mK, $R_{\text{peak}}^{\text{heat}} \lesssim 17$ (169) $h^{-1}$ Mpc and $\Delta R_{\text{FWHM}}^{\text{heat}} \lesssim 112$ (316) $h^{-1}$ Mpc, while $f_{\text{heat}}$ spans the entire range at 68% (95%) credible intervals. The disfavoured models at $z \sim 10$ in the *allUL*, *Single-z* case have $\bar{x}_{\text{HII}} \lesssim 0.01$ (0.35), $\bar{T}_K \lesssim 5$ (12) K, $-285$ ($-285$) $\lesssim \bar{\delta T}_b \lesssim -241$ ($-117$) mK, $f_{\text{heat}} \lesssim 0.03$ (0.36) $R_{\text{peak}}^{\text{heat}} \lesssim 7$ (17) $h^{-1}$ Mpc and $\Delta R_{\text{FWHM}}^{\text{heat}} \lesssim 12$ (36) $h^{-1}$ Mpc at 68% (95%) credible intervals. As expected, the disfavoured IGM parameter space is much larger in the *allUL* case compared to the LOFAR-only case at $z \sim 8$ and 10, as the
<table><caption>**Table B.1.** Constraints on the IGM parameters for the $A_r = 0$ and allUL scenario for both the Single-z and Joint-z analysis at three redshifts.</caption><thead><tr><th>IGM<br>Parameters<br>(MCMC<br>Model)</th><th>Credible<br>intervals</th><th>z = 8.3</th><th>z = 9.1</th><th>z = 10.1</th></tr></thead><tbody><tr><td rowspan="2">x&#773;<sub>HII</sub><br>(Single-z)</td><td>68%</td><td>[0, 0.05]</td><td>[0, 0.08]</td><td>[0, 0.01]</td></tr><tr><td>95%</td><td>[0, 0.78]</td><td>[0, 0.54]</td><td>[0, 0.35]</td></tr><tr><td rowspan="2">x&#773;<sub>HII</sub><br>(Joint-z)</td><td>68%</td><td>[0, 0.06]</td><td>[0, 0.02]</td><td>[0, 0.01]</td></tr><tr><td>95%</td><td>[0, 0.79]</td><td>[0, 0.42]</td><td>[0, 0.15]</td></tr><tr><td rowspan="2">f<sub>heat</sub><br>(Single-z)</td><td>68%</td><td>[X]</td><td>[X]</td><td>[0, 0.03]</td></tr><tr><td>95%</td><td>[X]</td><td>[X]</td><td>[0, 0.36]</td></tr><tr><td rowspan="2">f<sub>heat</sub><br>(Joint-z)</td><td>68%</td><td>[X]</td><td>[0, 0.12]</td><td>[0, 0.03]</td></tr><tr><td>95%</td><td>[X]</td><td>[X]</td><td>[X]</td></tr><tr><td rowspan="2">T&#773;<sub>K</sub><br>(Single-z)</td><td>68%</td><td>[1.8, 10]</td><td>[2.2, 7]</td><td>[2.6, 5]</td></tr><tr><td>95%</td><td>[1.8, 2180]</td><td>[2.2, 38]</td><td>[2.6, 12]</td></tr><tr><td rowspan="2">T&#773;<sub>K</sub><br>(Joint-z)</td><td>68%</td><td>[1.8, 13]</td><td>[2.2, 6]</td><td>[2.6, 5]</td></tr><tr><td>95%</td><td>[1.8, 1737]</td><td>[2.2, 270]</td><td>[2.6, 79]</td></tr><tr><td rowspan="2">δT&#773;<sub>b</sub><br>(Single-z)</td><td>68%</td><td>[-352, 3.2]</td><td>[-308, -134]</td><td>[-285, -241]</td></tr><tr><td>95%</td><td>[-375, 42]</td><td>[-344, 27]</td><td>[-285, -117]</td></tr><tr><td rowspan="2">δT&#773;<sub>b</sub><br>(Joint-z)</td><td>68%</td><td>[-349, 7]</td><td>[-307, -176]</td><td>[-285, -239]</td></tr><tr><td>95%</td><td>[-374, 45]</td><td>[-307, 13]</td><td>[-285, -15]</td></tr><tr><td rowspan="2">R<sup>heat</sup><sub>peak</sub><br>(h<sup>-1</sup> Mpc)<br>(Single-z)</td><td>68%</td><td>[0, 17]</td><td>[0, 21]</td><td>[0, 7]</td></tr><tr><td>95%</td><td>[0, 169]</td><td>[0, 44]</td><td>[0, 17]</td></tr><tr><td rowspan="2">R<sup>heat</sup><sub>peak</sub><br>(h<sup>-1</sup> Mpc)<br>(Joint-z)</td><td>68%</td><td>[0, 21]</td><td>[0, 7]</td><td>[0, 5]</td></tr><tr><td>95%</td><td>[0, 169]</td><td>[0, 63]</td><td>[0, 66]</td></tr><tr><td rowspan="2">ΔR<sup>heat</sup><sub>FWHM</sub><br>(h<sup>-1</sup> Mpc)<br>(Single-z)</td><td>68%</td><td>[6, 112]</td><td>[6, 35]</td><td>[2, 12]</td></tr><tr><td>95%</td><td>[3, 316]</td><td>[3, 104]</td><td>[3, 36]</td></tr><tr><td rowspan="2">ΔR<sup>heat</sup><sub>FWHM</sub><br>(h<sup>-1</sup> Mpc)<br>(Joint-z)</td><td>68%</td><td>[0, 21]</td><td>[8, 20]</td><td>[0, 11]</td></tr><tr><td>95%</td><td>[0, 295]</td><td>[3, 151]</td><td>[0, 66]</td></tr></tbody></table>

MCMC analysis is mainly driven by the strong HERA upper limits. The 68% and 95% credible intervals for the Joint-z allUL case are shown in Table B.1. The numbers show that the results for the Joint-z MCMC analysis are mainly driven by the HERA upper limits at $z \sim 8$.

Next we consider the Varying $A_r$ scenario, starting with the Single-z and allUL case. In this case, the disfavoured models at $z \sim 8$ show 68% credible intervals on the source parameters as $\zeta \le 1.7, M_{min} \ge 1.4 \times 10^{10} M_{\odot}$, while the credible intervals of other parameters span the entire range. The same intervals for $z \sim 10$ are $\zeta \le 3.3, M_{min} \ge 1.6 \times 10^{10} M_{\odot}$ and $f_X \le 5.5$, while the other two parameters' credible intervals span the entire ranges. The 1D marginalised probability distribution of the disfavoured IGM parameters for this case at $z \sim 8$ and 10 is shown in Fig. B.3 and B.4, respectively. The disfavoured models at $z \sim 8$ have $\bar{x}_{HII} \le 0.05 (0.78), \bar{T}_K \le 10 (2180) K$,

<figure><img src="image_15.png" /><figcaption>**Fig. B.1.** Same as Fig. B.2 but for $z \sim 8$.</figcaption></figure>

28 (2) $\le |\delta T_b| \le 1905$ (14791) mK, $R_{\text{peak}}^{\text{heat}} \le 8$ (138) $h^{-1}$ Mpc and $\Delta R_{\text{FWHM}}^{\text{heat}} \le 16$ (213) $h^{-1}$ Mpc, while $f_{\text{heat}}$ spans the entire range at 68% (95%) credible intervals for the Single-z analysis case. The disfavoured models at $z \sim 10$ in the allUL and Varying $A_r$ case have $\bar{x}_{HII} \le 0.01 (0.32), \bar{T}_K \le 5$ (27) K, 132 (91) $\le |\delta T_b| \le 9120$ (24888) mK, $f_{\text{heat}} \le 0.03(0.31) R_{\text{peak}}^{\text{heat}} \le 4$ (11) $h^{-1}$ Mpc and $\Delta R_{\text{FWHM}}^{\text{heat}} \le 5$ (20) $h^{-1}$ Mpc at 68% (95%) credible intervals for the Single-z analysis case.

The 68% and 95% credible intervals for the Joint-z, allUL case under the Varying $A_r$ scenario are shown in Table B.2. Similarly to the $A_r = 0$ scenario above, we find that the Joint-z MCMC analysis is also mostly driven by HERA upper limits at $z \sim 8$.

## Appendix C: Inference using Acharya et al. (2024c) upper limit results

The LOFAR upper limits considered in the main text, as well as those obtained by Mertens et al. (2020), contain excess noise and thus are conservative limits in nature. Recently, Acharya et al. (2024c) have re-analysed the 10 nights of data used in Mertens et al. (2020) with a different approach, and reported a much stronger $2\sigma$ upper limit of $\Delta^2(k =$
<figure><img src="image_16.png" /></figure>

**Fig. B.2.** Marginalised probability distribution of the disfavoured models' IGM parameters at $z \sim 10$ for the *allUL* and $A_r = 0$ scenario. The red and blue curves refer to the *Single-z* and *Joint-z* cases, respectively. The solid curves represent results when we consider LOFAR upper limits only, while the dotted curves refer to cases where we have used all available upper limit results from different radio interferometric observations.

0.075 $h$ Mpc$^{-1}$, $z = 9.1$, 141 hours) $\approx (25)^2$ mK$^2$ (see Table C1). This was obtained by assuming an accurate bias correction for the excess noise by subtracting the excess component. As this approach increases the risk of signal loss, more testing is required to verify the robustness of the results. Nevertheless, here we use these upper limits to explore their impact on the derived constraints for the physical state of the IGM.

Figure C.1 shows 1000 randomly chosen models from the GRIZZLY simulation database, together with the LOFAR $1\sigma$ upper limits from Mertens et al. (2025, red) and Acharya et al. (2024c, magenta, $\Delta^2_{21\text{-an}}$ with error $\Delta^2_{21\text{-an,err}}$). We see that the red dashed curve, corresponding to the power spectrum of a completely neutral and unheated IGM at $z = 9.1$ obtained assuming constant values of $x_{\text{HI}} = 1$ and $T_{\text{S}} = T_{\text{K}} = 2$ K, remains above the $1\sigma$ upper limits from Acharya et al. (2024c) at least for one k-bin. Unlike the results presented in the main text, here we find that all models are above $\Delta^2_{21\text{-an}} - \Delta^2_{21\text{-an,err}}$.

When performing the same MCMC analysis discussed in the main text to identify the disfavoured models, we find that the marginalised probabilities are high over the full parameter space, indicating that a large number of models are disfavoured. For this reason, differently from the analysis performed in the main text, here we choose to study the properties of the favoured models by

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

**Fig. B.3.** Same as Fig. B.1 but for the *Varying $A_r$* scenario.

adopting for the likelihood of a model surviving at $z$ the expression $\mathcal{L}_{\text{su,single-z}}(\theta, z) = 1 - \mathcal{L}_{\text{ex,single-z}}(\theta, z)$.

We perform the *Single-z* analysis following two approaches, the first of which is the same one adopted in the main text, i.e. we consider the upper limits at the three smallest k-bins. In the second case, instead, we consider only the smallest k-bin upper limit, which is the strongest, and thus it is the one that mainly affects the constraints on the source and IGM parameters. Using only one k-bin also have the advantage of breaking the correlation originating from the approach adopted in Acharya et al. (2024c) when using upper limits at different k-bins.

Figure C.2 shows the posterior distribution of the source parameters for the *Varying $A_r$* scenario at $z = 9.1$ obtained from the MCMC analysis using the upper limits from Acharya et al. (2024c) on the allowed models. The 68% credible intervals limits of the allowed models are $\zeta \lesssim 14$, $M_{\text{min}} \lesssim 5 \times 10^{11}$ M$_\odot$, $M_{\text{min,X}} \lesssim 1.5 \times 10^{10}$ M$_\odot$, $f_X \gtrsim 100$ and $A_r \lesssim 2$ for the *Varying $A_r$* scenario. The limits are similar when we use one upper limit at $k = 0.075 h$ Mpc$^{-1}$.

Figure C.3 shows the posterior distribution of the IGM parameters of the favoured models obtained from the *Single-z* analysis with the upper limits from Acharya et al. (2024c) at $z = 9.1$. The 95% credible intervals of the allowed models' IGM parameters are $\bar{x}_{\text{HI}} \lesssim 0.75$, $27 \lesssim \bar{T}_{\text{K}} \lesssim 7 \times 10^3$ K, $f_{\text{heat}} \gtrsim 0.62$, $4 \lesssim \delta\bar{T}_{\text{b}} \lesssim 44$ mK, $3 \lesssim R_{\text{peak}}^{\text{heat}} \lesssim 160$ $h^{-1}$ Mpc and $4 \lesssim \Delta R_{\text{FWHM}}^{\text{heat}} \lesssim 316$ $h^{-1}$ Mpc. The limits on the IGM parameters for the case where we used the upper limit at $k = 0.075 h$ Mpc$^{-1}$ only are similar to the case where we used three k-bins upper
<figure><img src="image_18.png" /><figcaption>**Fig. B.4.** Same as Fig. B.2 but for *Varying A*<sub>*r*</sub> scenario.</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>**Fig. C.1.** A subset of 1000 power spectra simulated with GRIZZLY at $z = 9.1$. These power spectra correspond to those with $A_r = 0$. The red and magenta points refer to the recent LOFAR $1\sigma$ upper limits from Mertens et al. (2025) and Acharya et al. (2024c), respectively. As a reference, the red dashed line corresponds to the power spectrum for a completely neutral IGM with no X-ray heating and $T_S = T_K$.</figcaption></figure>

limits. We note that, similar to the scenarios considered in the main text, these constraints are also highly dependent on the chosen priors on the source parameters (see Appendix A).

<figure><img src="image_20.png" /><figcaption>**Fig. C.2.** Posterior distribution of the source parameters at $z = 9.1$ for the *Varying A*<sub>*r*</sub> scenario. These models are favoured by the LOFAR upper limits from Acharya et al. (2024c). Magenta indicates the case in which upper limits in three *k*-bins are considered, while teal represents the case when only one *k*-bin is used. The contour levels in the two-dimensional contour plots refer to the $1\sigma$ and $2\sigma$ credible intervals of the models favoured by the LOFAR upper limit at this redshift. The diagonal panels represent the corresponding marginalised probability distributions of each parameter.</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>**Fig. C.3.** Posterior distribution of the IGM parameters at $z = 9.1$ for the *Varying A*<sub>*r*</sub> scenario. These posterior distributions correspond to the models that are favoured by the LOFAR upper limits from Acharya et al. (2024c). Magenta indicates the case in which upper limits in three *k*-bins are considered, while teal represents the case when only one *k*-bin is used. The contour levels in the two-dimensional contour plots refer to the $1\sigma$ and $2\sigma$ credible intervals of the models favoured by the LOFAR upper limits at this redshift. The diagonal panels represent the corresponding marginalised probability distributions of each parameter.</figcaption></figure>
<table><caption>**Table B2.** Constraints on the IGM parameters at three different redshifts for the *Varying A*<sub>r</sub> and *allUL* scenario for both the *Single-z* and *Joint-z* analysis.</caption><thead><tr><th>IGM<br>Parameters<br>(MCMC<br>Model)</th><th>Credible<br>intervals</th><th>z ~ 8</th><th>z ~ 9</th><th>z ~ 10</th></tr></thead><tbody><tr><td rowspan="2">x&#773;<sub>HII</sub><br>(Single-z)</td><td>68%</td><td>[0, 0.05]</td><td>[0, 0.02]</td><td>[0, 0.01]</td></tr><tr><td>95%</td><td>[0, 0.78]</td><td>[0, 0.46]</td><td>[0, 0.32]</td></tr><tr><td rowspan="2">x&#773;<sub>HII</sub><br>(Joint-z)</td><td>68%</td><td>[0, 0.06]</td><td>[0, 0.02]</td><td>[0, 0.004]</td></tr><tr><td>95%</td><td>[0, 0.8]</td><td>[0, 0.42]</td><td>[0, 0.16]</td></tr><tr><td rowspan="2">f<sub>heat</sub><br>(Single-z)</td><td>68%</td><td>X</td><td>[0, 0.06]</td><td>[0, 0.03]</td></tr><tr><td>95%</td><td>X</td><td>[0, 0.5]</td><td>[0, 0.31]</td></tr><tr><td rowspan="2">f<sub>heat</sub><br>(Joint-z)</td><td>68%</td><td>[X]</td><td>[0, 0.12]</td><td>[0, 0.03]</td></tr><tr><td>95%</td><td>[X]</td><td>[X]</td><td>[X]</td></tr><tr><td rowspan="2">T&#773;<sub>K</sub><br>(Single-z)</td><td>68%</td><td>[1.8, 10.0]</td><td>[2.2, 5]</td><td>[2.6, 5]</td></tr><tr><td>95%</td><td>[1.8, 2180]</td><td>[2.2, 55]</td><td>[2.6, 27]</td></tr><tr><td rowspan="2">T&#773;<sub>K</sub><br>(Joint-z)</td><td>68%</td><td>[1.8, 20]</td><td>[2.2, 6]</td><td>[2.6, 6]</td></tr><tr><td>95%</td><td>[1.8, 2089]</td><td>[2.2, 371]</td><td>[2.6, 110]</td></tr><tr><td rowspan="2">|δT&#773;<sub>b</sub>|<br>(mK)<br>(Single-z)</td><td>68%</td><td>[28, 11905]</td><td>[131, 3812]</td><td>[132, 9120]</td></tr><tr><td>95%</td><td>[2, 14791]</td><td>[62, 28184]</td><td>[91, 24888]</td></tr><tr><td rowspan="2">|δT&#773;<sub>b</sub>|<br>(mK)<br>(Joint-z)</td><td>68%</td><td>[18, 1584]</td><td>[51, 2238]</td><td>[77, 2630]</td></tr><tr><td>95%</td><td>[3, 15848]</td><td>[12, 19952]</td><td>[19, 31622]</td></tr><tr><td rowspan="2">R<sup>heat</sup><sub>peak</sub><br>(h<sup>-1</sup> Mpc)<br>(Single-z)</td><td>68%</td><td>[0, 8.0]</td><td>[0, 5]</td><td>[0, 4.0]</td></tr><tr><td>95%</td><td>[0, 138]</td><td>[0, 17]</td><td>[0, 11]</td></tr><tr><td rowspan="2">R<sup>heat</sup><sub>peak</sub><br>(h<sup>-1</sup> Mpc)<br>(Joint-z)</td><td>68%</td><td>[0, 20]</td><td>[0, 5]</td><td>[0, 3.0]</td></tr><tr><td>95%</td><td>[0, 138]</td><td>[0, 35]</td><td>[0, 15]</td></tr><tr><td rowspan="2">ΔR<sup>heat</sup><sub>FWHM</sub><br>(h<sup>-1</sup> Mpc)<br>(Single-z)</td><td>68%</td><td>[0, 16]</td><td>[0, 9]</td><td>[0, 5]</td></tr><tr><td>95%</td><td>[0, 213]</td><td>[0, 76]</td><td>[0, 20]</td></tr><tr><td rowspan="2">ΔR<sup>heat</sup><sub>FWHM</sub><br>(h<sup>-1</sup> Mpc)<br>(Joint-z)</td><td>68%</td><td>[0, 20]</td><td>[0, 9]</td><td>[0, 6]</td></tr><tr><td>95%</td><td>[0, 218]</td><td>[0, 76]</td><td>[0, 30]</td></tr></tbody></table>

<table><caption>**Table C1.** Upper limit on Δ² as reported in Acharya et al. (2024c) at z = 9.1 for different k-bins.</caption><thead><tr><th>k (h Mpc<sup>-1</sup>)</th><th>Δ<sup>2</sup><sub>21-an</sub> (mK<sup>2</sup>)</th><th>Δ<sup>2</sup><sub>21-an, err</sub> (mK<sup>2</sup>)</th></tr></thead><tbody><tr><td>0.076</td><td>(5.98)<sup>2</sup></td><td>(15.42)<sup>2</sup></td></tr><tr><td>0.101</td><td>(8.16)<sup>2</sup></td><td>(18.7)<sup>2</sup></td></tr><tr><td>0.133</td><td>(10.77)<sup>2</sup></td><td>(20.89)<sup>2</sup></td></tr><tr><td>0.181</td><td>(15.24)<sup>2</sup></td><td>(25.13)<sup>2</sup></td></tr><tr><td>0.240</td><td>(17.6)<sup>2</sup></td><td>(24.47)<sup>2</sup></td></tr><tr><td>0.319</td><td>(20.88)<sup>2</sup></td><td>(28.15)<sup>2</sup></td></tr><tr><td>0.436</td><td>(27.12)<sup>2</sup></td><td>(41.63)<sup>2</sup></td></tr></tbody></table>

**Notes.** The different columns show k-bins, mean value Δ²<sub>21-an</sub>(k, z) and corresponding 1σ error Δ²<sub>21-an, err</sub>(k, z).