© The Authors 2025

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 I) in the intergalactic 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 depth (e.g. Planck Collaboration VI 2020), high-redshift (z) Lyα 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α 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 (${\overline{\delta T}}_{\mathrm{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¹; 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²; Parsons et al. 2014; Kolopanis et al. 2019), the Murchison Widefield Array (MWA³; Tingay et al. 2013), the Hydrogen Epoch of Reionization Array (HERA⁴; DeBoer et al. 2017), and the New Extension in Nançay Upgrading LOFAR (NenuFAR⁵; Munshi et al. 2024). Around the start of the next decade, the Square Kilometre Array (SKA)⁶ 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 polarisation 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 ${\overline{\delta T}}_{\mathrm{b}}$ with a minimum at z ≈ 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 ${\overline{\delta T}}_{\mathrm{b}}$ profile around z ∼ 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 Δ², which contains much more information than the global signal. At the current time, such observations have only produced upper limits. A 2σ upper limit of (248)² mK² for k = 0.5 h Mpc⁻¹ at z = 8.6 from the GMRT was the first of its kind (Paciga et al. 2013). Later, PAPER⁷ (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, Δ²(k = 0.075 h Mpc⁻¹, z = 9.1, 141 hours)≈(73)² mK², Δ²(k = 0.34 h Mpc⁻¹, z = 7.9, 1128 hours)≈(21.4)² mK², and Δ²(k = 0.14 h Mpc⁻¹, z = 6.5, 110 hours)≈(43)² mK² 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 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σ upper limit of Δ²(k = 0.076 h Mpc⁻¹, z = 9.1, 140 hours)≈(54.3)² mK². 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 Ω_m = 0.27, Ω_Λ = 0.73, Ω_B = 0.044, and h = 0.7 (Wilkinson Microwave Anisotropy Probem 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 Δ₂₁²(k, z) and associated 1σ error Δ_21, err²(k, z). The best result is at z = 9.1, with a 2σ upper limit of (54.3 mK)² at k = 0.075 h Mpc⁻¹, 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⁵) 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α 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⁸ project PRACE4LOFAR. The dark matter-only simulation was run using the cubep³m (Harnois-Déraps et al. 2013) code with particle mass resolution of 4.05 × 10⁷ M_⊙. The simulation cube has a size of 500 h⁻¹ Mpc, corresponding to a field-of-view of 4.3° ×4.3° at z≈ 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 Mpc⁻¹, requiring a size of at least ≈250 h⁻¹ Mpc.

The density and velocity fields were smoothed onto a 300³ 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 lists consist of halos with more than 25 particles, which leads to a minimum dark matter halo mass ≈10⁹ M_⊙.

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α, 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_⋆ in a dark matter halo with mass M_halo is $M_{\star }=f_{\star }\frac{{\mathrm{\Omega }}_{\mathrm{B}}}{{\mathrm{\Omega }}_{\mathrm{m}}}{M}_{\mathrm{halo}}$, where f_⋆ is the star formation efficiency. In all our simulations we used f_⋆ = 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 (ζ): This parameter regulates the emission rate of ionizing photons escaping into the IGM from a halo⁹. 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}{\mathrm{s}}^{-1}{{\mathrm{M}}_{\odot }}^{-1}$, where we vary log₁₀(ζ) 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}{\mathrm{s}}^{-1}{{\mathrm{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)∝E^−α, where α = 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₁₀(f_X) in the range [–3, 3].

(3) Minimum mass of UV emitting halos (M_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_min contribute to the ionizing photon budget. We vary log₁₀(M_min/M_⊙) 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 ΛCDM cosmology is larger at the low mass end, star formation in halos with mass ≲ 10⁹ M_⊙ 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_min, X): This determines the number density of X-ray emitting halos, as we assume that only halos with a mass larger than M_min, X host X-ray sources. Also for log₁₀(M_min, X/M_⊙) we explore the range [9, 12].

(5) Radio background efficiency (A_r): We assumed the existence of a radio background¹⁰ 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)

$$\begin{array}{c}\hfill T_{\gamma ,\mathrm{eff}}=T_{\gamma }[1+A_r{\left(\frac{{\nu }_{\mathrm{obs}}}{78\mathrm{MHz}}\right)}^{-2.6}].\end{array}$$(1)

Here, T_γ(z) = 2.725 (1 + z) K is the usual CMB temperature and ν_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_γ, eff = T_γ), 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 (δT_b) can be written as (see e.g, Madau et al. 1997; Furlanetto et al. 2006)

$$\begin{array}{cc}\hfill \delta T_{\mathrm{b}}(\mathbf{x},z)& =27x_{\mathrm{HI}}(\mathbf{x},z)[1+{\delta }_{\mathrm{B}}(\mathbf{x},z)](1-\frac{T_{\gamma ,\mathrm{eff}}(z)}{T_{\mathrm{S}}(\mathbf{x},z)})\hfill \\ \hfill & \times \left(\frac{{\mathrm{\Omega }}_{\mathrm{B}}{h}^2}{0.023}\right){\left(\frac{0.15}{{\mathrm{\Omega }}_{\mathrm{m}}{h}^2}\frac{1+z}{10}\right)}^{1/2}\left(\frac{H}{\mathrm{d}{v}_{\mathrm{r}}/\mathrm{d}r+H}\right)\mathrm{mK}.\hfill \end{array}$$(2)

Here, [1 + δ_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α 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 δ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 δT_b, denoted as Δ²(k)¹¹.

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_γ, 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¹² for each redshift:

• ${\overline{x}}_{\text{HII}}$: Volume-averaged ionized fraction.

• ${\overline{T}}_{\text{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 ${\overline{T}}_{\mathrm{K}}$. In general, a higher x_HII threshold results in a larger value of ${\overline{T}}_{\mathrm{K}}$.

• ${\overline{\delta T}}_{\text{b}}$ (mK): Mass-averaged differential brightness temperature.

• f_heat: Volume fraction of heated regions (T_K > T_γ, eff).

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

• $\mathrm{\Delta }{R}_{\text{FWHM}}^{\text{heat}}$ (h⁻¹ 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¹³ (Lewis & Bridle 2002). The likelihood of a set of parameters θ to be excluded at z for an upper limit of the power spectrum Δ₂₁²(k_i, z)±Δ_21, err²(k_i, z), can be defined as (see Ghara et al. 2020, for details)

$$\begin{array}{c}\hfill {\mathcal{L}}_{\mathrm{ex},\mathrm{single}-z}(\mathit{\theta },z)=1-{\displaystyle \prod _i}\frac{1}{2}[1+\mathrm{erf}\left(\frac{{\mathrm{\Delta }}_{21}^2(k_i,z)-{\mathrm{\Delta }}_{\mathrm{m}}^2(k_i,\mathit{\theta },z)}{\sqrt{2}\sigma (k_i,z)}\right)].\end{array}$$(3)

Here, Δ_m²(k, θ, z) is the model power spectrum for a set of parameters θ, while σ is the error that includes contributions from the observational error (Δ_21, err²) as well as from the power spectrum modelling error (Δ_m, err²): $\sigma =\sqrt{{\mathrm{\Delta }}_{21,\mathrm{err}}^4+{\mathrm{\Delta }}_{\mathrm{m},\mathrm{err}}^4}$.

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

$$\begin{array}{c}\hfill {\mathcal{L}}_{\mathrm{ex},\mathrm{joint}-z}(\mathit{\theta })=1-{\displaystyle \prod _{i,j}}\frac{1}{2}[1+\mathrm{erf}\left(\frac{{\mathrm{\Delta }}_{21}^2(k_i,z_j)-{\mathrm{\Delta }}_{\mathrm{m}}^2(k_i,\mathit{\theta },z_j)}{\sqrt{2}\sigma (k_i,z_j)}\right)].\end{array}$$(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⁵ sets of power spectra and IGM parameters at each of the three redshifts by sampling the 5D source parameter space¹⁴. 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, Δ_m, err²(k_i, z) = 0.3 × Δ_m²(k_i, θ, z). This includes the inherent modelling errors due to the approximate nature of the GRIZZLY algorithm¹⁵, as well as the errors introduced by the interpolation scheme. To run the MCMC on the source parameters, we used ten walkers and 10⁶ 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 ℒ_{ex, single − z}(θ, z) and ℒ_{ex, joint − z}(θ).

In this analysis we also use a prior on ${\overline{x}}_{\mathrm{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 (${\overline{x}}_{\mathrm{HII},\max }$) until redshift 6, and fully ionized after that¹⁶. The value ${\overline{x}}_{\mathrm{HII},\max }$ follows from assuming the 1στ measurement, i.e. τ = 0.061 for τ = 0.054 ± 0.007 (Planck Collaboration VI 2020). It represents the largest possible of value of ${\overline{x}}_{\mathrm{HII}}$ at redshift z_i consistent with the τ measurements. We obtain values of ${\overline{x}}_{\mathrm{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 ${\overline{x}}_{\mathrm{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 Δ₂₁²(k, z)−Δ_21, err²(k, z) from the recent LOFAR 1σ 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 × (1 + z)² 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.

Fig. 1. Set of 1000 power spectra randomly chosen from the 105 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 (Ar = 0). The down arrow points refer to the recent LOFAR 1σ upper limit (Δ212(k, z)±Δ21, err2(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 TS = TK. The coloured power spectra are those with a value larger than Δ212(k, z)−Δ21, err2(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.

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 Δ₂₁²(k, z)−Δ_21, err²(k, z) in at least one k-bin. The vertical axes represent the number of disfavoured models out of the total of 10⁵. 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_min and M_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 ζ and f_X, suggesting that the rare ionized or heated regions are large in size.

Fig. 2. Histogram of the disfavoured values of the source parameters for the scenario with Ar = 0. These correspond to the models with power spectra values larger than Δ212(k, z)−Δ21, err2(k, z) in at least one k-bin, and they have a high probability of being ruled out by the recent LOFAR 1σ upper limits (Δ212(k, z)±Δ21, err2(k, z)) from Mertens et al. (2025).

In Fig. 3, we present the distribution of IGM parameters for the same set of models considered in Fig. 2. A significant number of disfavoured models at all redshifts have ${\overline{x}}_{\mathrm{HII}}$ and f_heat smaller than ∼0.6. The characteristic sizes in these models are several tens of h⁻¹ 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 ${\overline{x}}_{\mathrm{HII}}\lesssim 0.1$, f_heat ≲ 0.1 and ${\overline{T}}_{\mathrm{K}}\lesssim 4$ K are strongly influenced by the priors chosen on the source parameters (see Appendix A for details).

Fig. 3. Histogram of the disfavoured values of the IGM parameters for the scenario with Ar = 0. These correspond to the same set of models considered in Fig. 2.

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 ℒ_{ex, single − z}(θ, 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 ℒ_{ex, joint − z}(θ) 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 ζ, M_min, M_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 parameters depend on the chosen prior ranges. We find that the 68% disfavoured credible intervals¹⁷ from the Single-z analysis at z = 9.1 are $M_{\min }\gtrsim 3\times {10}^{10}{\mathrm{M}}_{\odot }$ and $M_{\min ,\mathrm{X}}\gtrsim 6\times {10}^{10}{\mathrm{M}}_{\odot }$, while the limits for ζ 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.

Fig. 4. Posterior distribution of the disfavoured values of the source parameters for the Ar = 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σ and 2σ 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.

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_min and M_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₁₀(M_min) or log₁₀(M_min, X)≳10.5 M_⊙. This suggests that models associated with a low number density of massive UV and X-ray emitting sources are disfavoured. The limits for ζ, on the other hand, differ among the various redshifts. 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.

Fig. 5. Marginalised probability distribution of the disfavoured values of the source parameters for the Ar = 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.

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 ζ) 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 Δ² is dominated by x_HII fluctuations and reaches the largest amplitude for ${\overline{x}}_{\mathrm{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 ${\overline{x}}_{\mathrm{HII}}\lesssim 0.13(0.55)$, ${\overline{T}}_{\mathrm{K}}\lesssim 7.3(21)$ K, f_heat ≲ 0.18 (0.57), $-276(-317)\lesssim {\overline{\delta T}}_{\mathrm{b}}\lesssim -123(-61)$ mK, $R_{\mathrm{peak}}^{\mathrm{heat}}\lesssim 28(40)h^{-1}\mathrm{Mpc}$ and $\mathrm{\Delta }{R}_{\mathrm{FWHM}}^{\mathrm{heat}}\lesssim 46(81)h^{-1}\mathrm{Mpc}$ at 68% (95%) credible intervals for the Single-z analysis case. The PDF of ${\overline{x}}_{\mathrm{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).

Fig. 6. Posterior distribution of the disfavoured values of the IGM parameters for the Ar = 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σ and 2σ 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.

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 δT_b maps, producing large Δ² values at scales 0.076 h Mpc⁻¹ ≲ k ≲ 0.133 h Mpc⁻¹. Unlike the Single-z case, in the Joint-z analysis the IGM parameter’s disfavoured credible intervals show a consistent evolution with redshift. For example, the 95% credible intervals limits on ${\overline{x}}_{\mathrm{HII}}$ are ≲0.85, ≲0.46 and ≲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 (VaryingA_r)

Next, we present the results of the MCMC analysis for the 5D source parameter space, i.e. ζ, 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 ζ ≲ 2.4, $M_{\min }\gtrsim 1.7\times {10}^{10}{\mathrm{M}}_{\odot }$, $M_{\min ,\mathrm{X}}\gtrsim 1.6\times {10}^{10}{\mathrm{M}}_{\odot }$, f_X ≲ 5.2 and A_r ≳ 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 δT_b = 0, while 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 ${\overline{x}}_{\mathrm{HI}}\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.

Fig. 7. Posterior distribution of the source parameters disfavoured values for the Varying Ar 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.

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 ζ 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_γ, eff ≳ 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_γ, eff ≳ 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 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 ${\overline{x}}_{\mathrm{HII}}\lesssim 0.02(0.46)$, ${\overline{T}}_{\mathrm{K}}\lesssim 4.4(44)$ K, f_heat ≲ 0.05 (0.46), $1.3\times {10}^4(2.6\times {10}^4)\gtrsim |{\overline{\delta T}}_{\mathrm{b}}|\gtrsim 158(91)$ mK, $R_{\mathrm{peak}}^{\mathrm{heat}}\lesssim 3(14)h^{-1}\mathrm{Mpc}$ and $\mathrm{\Delta }{R}_{\mathrm{FWHM}}^{\mathrm{heat}}\lesssim 5(25)h^{-1}\mathrm{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).

Fig. 8. Posterior distribution of the IGM parameters of the disfavoured models at z = 9.1 for the Varying Ar 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.

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 ≳ 0.3 h Mpc⁻¹. 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⁵ simulations for A_r = 0. We note that at z ∼ 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 Δ²(k = 0.34 h Mpc⁻¹, z = 7.9, 1128 hours)≈(21.4)² mK², 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), Δ²(k = 0.14 h Mpc⁻¹,  z = 7.8, 51 hours)≈(154)² mK², is still above the extreme EoR models power spectra by a factor of ≈2. The best 2σ result from PAPER is Δ²(k = 0.37 h Mpc⁻¹,  z = 8.4, 135 nights)≈(205)² mK², which is larger than the predicted power spectra by a factor of ≳4.

Fig. 9. Set of 1000 power spectra randomly chosen from the 105 simulated spectra for Ar = 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 TS = TK at those redshifts. The down arrow points refer to the different 2σ upper limits at those redshifts.

At z ∼ 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σ upper limit from GMRT (Paciga et al. 2013) and MWA (Trott et al. 2020) at this redshift are Δ²(k = 0.5 h Mpc⁻¹,  z = 8.6, 50 hours)≈(248)² mK² and Δ²(k = 0.14 h Mpc⁻¹,  z = 8.8, 51 hours)≈(250)² mK², respectively, i.e. more than ≈6 times larger than our extreme models. The best 2σ result from PAPER is Δ²(k = 0.36 h Mpc⁻¹,  z = 8.7, 135 nights)≈(281)² mK², which is larger than the GRIZZLY power spectra by a factor of at least ≈8. We note that the neutral and unheated IGM scenario at z ∼ 9 (dashed curve) hasn’t been ruled out by any of these upper limits.

Similarly to z ∼ 8, also at z ∼ 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σ upper limit result from HERA at this redshift is Δ²(k = 0.36 h Mpc⁻¹,  z = 10.4, 1128 hours)≈(59)² mK², which is marginally able to rule out the neutral and unheated EoR scenario. Compared to z ∼ 8, the number of models ruled out by HERA at z ∼ 10 is smaller. The only other interferometer which has collected data at this redshift is PAPER, with a best 2σ result of Δ²(k = 0.34 h Mpc⁻¹,  z = 9.9, 135 nights)≈(1923)² mK², 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 ∼ 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 ∼ 9 for the allUL case are $M_{\min }\gtrsim 2.4\times {10}^{10}{\mathrm{M}}_{\odot }$ and $M_{\min ,\mathrm{X}}\gtrsim 3.1\times {10}^{10}{\mathrm{M}}_{\odot }$, while the limits for ζ 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 $\mathit{allUL}$ analysis case, the 68% disfavoured credible intervals at z ∼ 9 are ζ ≲ 1.8 and $M_{\min }\gtrsim 1.4\times {10}^{10}{\mathrm{M}}_{\odot }$, while the limits for M_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 ∼ 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.

Fig. 10. Marginalised probability distribution of the disfavoured models’ source parameters at z ∼ 8, 9, and 10 for allUL case. The results are for the Ar = 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.

In Fig. 11 we present the 1D marginalised probability distribution of the disfavoured models’ IGM parameters at z ∼ 9 for the allUL scenario with A_r = 0, which in the Single-z analysis case have ${\overline{x}}_{\mathrm{HII}}\lesssim 0.08(0.54)$, ${\overline{T}}_{\mathrm{K}}\lesssim 7(38)$ K, $-308(-344)\lesssim {\overline{\delta T}}_{\mathrm{b}}\lesssim -134(27)$ mK, $R_{\mathrm{peak}}^{\mathrm{heat}}\lesssim 21(44)h^{-1}\mathrm{Mpc}$ and $\mathrm{\Delta }{R}_{\mathrm{FWHM}}^{\mathrm{heat}}\lesssim 35(104)h^{-1}\mathrm{Mpc}$ at 68% (95%) credible intervals, while f_heat spans the entire range. With the Joint-z analysis we instead obtain ${\overline{x}}_{\mathrm{HII}}\lesssim 0.02(0.42)$, ${\overline{T}}_{\mathrm{K}}\lesssim 6(270)$ K, $-307(-307)\lesssim {\overline{\delta T}}_{\mathrm{b}}\lesssim -176(13)$ mK, f_heat ≲ 0.12(X), $R_{\mathrm{peak}}^{\mathrm{heat}}\lesssim 7(63)h^{-1}\mathrm{Mpc}$ and $\mathrm{\Delta }{R}_{\mathrm{FWHM}}^{\mathrm{heat}}\lesssim 20(151)h^{-1}\mathrm{Mpc}$. While the constraints on ${\overline{x}}_{\mathrm{HII}},f_{\mathrm{heat}},{\overline{T}}_{\mathrm{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.

Fig. 11. Marginalised probability distribution of the IGM parameters of the disfavoured models at z ∼ 9 for allUL and Ar = 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.

Finally, Fig. 12 presents the 1D marginalised probability distribution of the disfavoured IGM parameters at z ∼ 9 for the allUL and Varying A_r scenario. The 68% (95%) credible intervals of the disfavoured models for the Single-z analysis are ${\overline{x}}_{\mathrm{HII}}\lesssim 0.02(0.46)$, f_heat ≲ 0.06 (0.5) ${\overline{T}}_{\mathrm{K}}\lesssim 5(55)$ K, $131(62)\lesssim |{\overline{\delta T}}_{\mathrm{b}}|\lesssim 3812(28184)$ mK, $R_{\mathrm{peak}}^{\mathrm{heat}}\lesssim 5(17)h^{-1}\mathrm{Mpc}$ and $\mathrm{\Delta }{R}_{\mathrm{FWHM}}^{\mathrm{heat}}\lesssim 7(32)h^{-1}\mathrm{Mpc}$. For the Joint-z analysis, these intervals instead are ${\overline{x}}_{\mathrm{HII}}\lesssim 0.02(0.42)$, ${\overline{T}}_{\mathrm{K}}\lesssim 6(371)$ K, $51(13)\lesssim |{\overline{\delta T}}_{\mathrm{b}}|\lesssim 2238(19952)$ mK, f_heat ≲ 0.12(X), $R_{\mathrm{peak}}^{\mathrm{heat}}\lesssim 5(35)h^{-1}\mathrm{Mpc}$ and $\mathrm{\Delta }{R}_{\mathrm{FWHM}}^{\mathrm{heat}}\lesssim 9(76)h^{-1}\mathrm{Mpc}$. Table B2 presents the 68% and 95% 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_r scenario also play an important role in producing similar posterior distributions.

Fig. 12. Same as Fig. 11 but for the Varying Ar scenario.

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_r = 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_r 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 databse 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 (ζ), minimum mass of the UV emitting halos (M_min), minimum mass of X-ray emitting halos (M_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 (${\overline{x}}_{\mathrm{HII}}$), volume-averaged gas temperature of the partially ionized IGM (${\overline{T}}_{\mathrm{K}}$), mass-averaged brightness temperature (${\overline{\delta T}}_{\mathrm{b}}$), volume fraction of the heated region (f_heat), size of the heated regions at which the PDF of the sizes peak ($R_{\mathrm{peak}}^{\mathrm{heat}}$), and FWHM of the PDFs ($\mathrm{\Delta }{R}_{\mathrm{FWHM}}^{\mathrm{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 Δ₂₁²(k, z)−Δ_21, err²(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_min and M_min, X values at z = 9.1 are $\gtrsim 3\times {10}^{10}{\mathrm{M}}_{\odot }$ and $\gtrsim 6\times {10}^{10}{\mathrm{M}}_{\odot }$, respectively, while ζ 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 ${\overline{x}}_{\mathrm{HII}}\lesssim 0.13(0.55)$, ${\overline{T}}_{\mathrm{K}}\lesssim 7.3(21)$ K, f_heat ≲ 0.18 (0.57), $-276(-317)\lesssim {\overline{\delta T}}_{\mathrm{b}}\lesssim -123(-61)$ mK, $R_{\mathrm{peak}}^{\mathrm{heat}}\lesssim 28(40)h^{-1}\mathrm{Mpc}$, and $\mathrm{\Delta }{R}_{\mathrm{FWHM}}^{\mathrm{heat}}\lesssim 46(81)h^{-1}\mathrm{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 ${\overline{x}}_{\mathrm{HII}}$ are ≲0.85, ≲0.46, and ≲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 ζ ≲ 2.4, $M_{\min }\gtrsim 1.7\times {10}^{10}{\mathrm{M}}_{\odot }$, $M_{\min ,\mathrm{X}}\gtrsim 1.6\times {10}^{10}{\mathrm{M}}_{\odot }$, f_X ≲ 5.2, and A_r ≳ 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 ≳ 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 ≳57% of the CMB at 1.42 GHz.

• For such a radio background model, the disfavoured limits on the IGM parameters are ${\overline{x}}_{\mathrm{HII}}\lesssim 0.02(0.46)$, ${\overline{T}}_{\mathrm{K}}\lesssim 4.4(44)$ K, f_heat ≲ 0.05 (0.46), $1.3\times {10}^4(2.6\times {10}^4)\gtrsim |{\overline{\delta T}}_{\mathrm{b}}|\gtrsim 158(91)$ mK, $R_{\mathrm{peak}}^{\mathrm{heat}}\lesssim 3(14)h^{-1}\mathrm{Mpc}$, and $\mathrm{\Delta }{R}_{\mathrm{FWHM}}^{\mathrm{heat}}\lesssim 5(25)h^{-1}\mathrm{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 ∼ 8 and 10, while the LOFAR ones prevail at z ∼ 9. The HERA upper limits rule out a neutral and unheated IGM at z ∼ 8 at a 2σ level. We find that the constraints on ${\overline{x}}_{\mathrm{HII}},f_{\mathrm{heat}}$, and ${\overline{T}}_{\mathrm{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_heat ≲ 0.34, ${\overline{T}}_{\mathrm{K}}\lesssim 160$ K, $R_{\mathrm{peak}}^{\mathrm{heat}}\lesssim 70h^{-1}\mathrm{Mpc}$ and $\mathrm{\Delta }{R}_{\mathrm{FWHM}}^{\mathrm{heat}}\lesssim 110h^{-1}\mathrm{Mpc}$, while we found two regimes of disfavoured credible intervals for ${\overline{x}}_{\mathrm{HII}}$, i.e. ${\overline{x}}_{\mathrm{HII}}\lesssim 0.08$ and $0.45\lesssim {\overline{x}}_{\mathrm{HII}}\lesssim 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⁹ M_⊙ while varying ζ, 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 ${\overline{x}}_{\mathrm{HII}}$, f_heat, and ${\overline{T}}_{\mathrm{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σ upper limit of Δ²(k = 0.075 h Mpc⁻¹,  z = 9.1)≈(25)² mK² 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 ${\overline{x}}_{\mathrm{HII}}\lesssim 0.75$, $27\lesssim {\overline{T}}_{\mathrm{K}}\lesssim 7\times {10}^3$ K, f_heat ≳ 0.62, $4\lesssim {\overline{\delta T}}_{\mathrm{b}}\lesssim 44$ mK, $3\lesssim R_{\mathrm{peak}}^{\mathrm{heat}}\lesssim 160h^{-1}\mathrm{Mpc}$, and $4\lesssim \mathrm{\Delta }{R}_{\mathrm{FWHM}}^{\mathrm{heat}}\lesssim 316h^{-1}\mathrm{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 Δ² 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 Δ² 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 ≳ 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.

Acknowledgments

RG acknowledges support from SERB, DST Ramanujan Fellowship no. RJF/2022/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$\stackrel{`}{\mathrm{e}}$s Grand Centre de Calcul (TGCC) operated by CEA near Paris, France and Marenostrum 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

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

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 ℒex, single − z(θ, z = 9.1) = 0.5 and ℒex, joint − z(θ, z = 9.1) = 0.5.

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 ${\overline{x}}_{\mathrm{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 ${\overline{x}}_{\mathrm{HII}}\lesssim 0.024(0.56)$, ${\overline{T}}_{\mathrm{K}}\lesssim 7.2(618)$ K, no constraints on f_heat, $-315\lesssim {\overline{\delta T}}_{\mathrm{b}}\lesssim -129(18.1)$ mK, $R_{\mathrm{peak}}^{\mathrm{heat}}\lesssim 10.8(83.1)h^{-1}\mathrm{Mpc}$ and $\mathrm{\Delta }{R}_{\mathrm{FWHM}}^{\mathrm{heat}}\lesssim 87(207)h^{-1}\mathrm{Mpc}$ at 68% (95%) credible intervals for the Single-zno-LOFAR analysis case. For the Joint-zno-LOFAR analysis case, the disfavoured models’ 68% (95%) credible intervals are ${\overline{x}}_{\mathrm{HII}}\lesssim 0.019(0.43)$, ${\overline{T}}_{\mathrm{K}}\lesssim 6(562)$ K, no constraints on f_heat, $-315\lesssim {\overline{\delta T}}_{\mathrm{b}}\lesssim -149(18.1)$ mK, $R_{\mathrm{peak}}^{\mathrm{heat}}\lesssim 7(36)h^{-1}\mathrm{Mpc}$ and $\mathrm{\Delta }{R}_{\mathrm{FWHM}}^{\mathrm{heat}}\lesssim 23(93)h^{-1}\mathrm{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 ${\overline{\delta T}}_{\mathrm{b}}$, $R_{\mathrm{peak}}^{\mathrm{heat}}$ and $\mathrm{\Delta }{R}_{\mathrm{FWHM}}^{\mathrm{heat}}$, they are less prominent for ${\overline{x}}_{\mathrm{HII}}$, f_heat and ${\overline{T}}_{\mathrm{K}}$. This suggests that the LOFAR upper limits provide very weak constraints on the latter three IGM parameters. The fact that the posterior for ${\overline{\delta T}}_{\mathrm{b}}$ is affected by the LOFAR upper limits but the other averages (${\overline{x}}_{\mathrm{HII}}$, f_heat and ${\overline{T}}_{\mathrm{K}}$) are not, suggests that the correlations between δ, x_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 ∼ 8 and 10

While the results from the MCMC analysis which includes all available upper limits (allUL case) at z ∼ 9 have been discussed in Sect. 3.3, here we present the results at z ∼ 8 and 10. The 68% disfavoured credible intervals of the source parameters for the A_r = 0 and Single-z scenario at z ∼ 8 are ζ ≲ 1.6 and $M_{\min }\gtrsim 1.4\times {10}^{10}{\mathrm{M}}_{\odot }$, while the limits for M_min, X and f_X span the entire parameter ranges (see Fig. 10). The same intervals for z ∼ 10 are ζ ≲ 5,$M_{\min }\gtrsim 1.7\times {10}^{10}{\mathrm{M}}_{\odot }$, $M_{\min ,\mathrm{X}}\gtrsim 1.6\times {10}^{10}{\mathrm{M}}_{\odot }$, and f_X ≲ 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 ∼ 8 and 10 for the $\mathit{allUL}$ case are significantly different from the LOFAR-only case, as HERA upper limits dominate the MCMC analysis at these two redshifts.

Fig. B.1. Same as Fig. B.2 but for z ∼ 8.

Fig. B.2. Marginalised probability distribution of the disfavoured models’ IGM parameters at z ∼ 10 for the allUL and Ar = 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.

Fig. B.3. Same as Fig. B.1 but for the Varying Ar scenario.

Fig. B.4. Same as Fig. B.2 but for Varying Ar scenario.

The 1D marginalised probability distribution of the disfavoured IGM parameters for the A_r = 0 and allUL case at z ∼ 8 are shown in Fig. B.1 while Fig. B.2 shows the same but for z ∼ 10. The disfavoured models at z ∼ 8 in the allUL, Single-z and A_r = 0 scenario have ${\overline{x}}_{\mathrm{HII}}\lesssim 0.05(0.78)$, ${\overline{T}}_{\mathrm{K}}\lesssim 10(2180)$ K, $-352(-375)\lesssim {\overline{\delta T}}_{\mathrm{b}}\lesssim 3.2(42)$ mK, $R_{\mathrm{peak}}^{\mathrm{heat}}\lesssim 17(169)h^{-1}\mathrm{Mpc}$ and $\mathrm{\Delta }{R}_{\mathrm{FWHM}}^{\mathrm{heat}}\lesssim 112(316)h^{-1}\mathrm{Mpc}$, while f_heat spans the entire range at 68% (95%) credible intervals. The disfavoured models at z ∼ 10 in the allULSingle-z case have ${\overline{x}}_{\mathrm{HII}}\lesssim 0.01(0.35)$, ${\overline{T}}_{\mathrm{K}}\lesssim 5(12)$ K, $-285(-285)\lesssim {\overline{\delta T}}_{\mathrm{b}}\lesssim -241(-117)$ mK, f_heat ≲ 0.03(0.36) $R_{\mathrm{peak}}^{\mathrm{heat}}\lesssim 7(17)h^{-1}\mathrm{Mpc}$ and $\mathrm{\Delta }{R}_{\mathrm{FWHM}}^{\mathrm{heat}}\lesssim 12(36)h^{-1}\mathrm{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 ∼ 8 and 10, as the MCMC analysis is mainly driven by the strong HERA upper limits. The 68% and 95% credible intervals for the Joint-zallUL 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 ∼ 8.

Next we consider the Varying A_r scenario, starting with the Single-z and $\mathit{allUL}$ case. In this case, the disfavoured models at z ∼ 8 show 68% credible intervals on the source parameters as ζ ≲ 1.7,$M_{\min }\gtrsim 1.4\times {10}^{10}{\mathrm{M}}_{\odot }$, while the credible intervals of other parameters span the entire range. The same intervals for z ∼ 10 are ζ ≲ 3.3,$M_{\min }\gtrsim 1.6\times {10}^{10}{\mathrm{M}}_{\odot }$ and f_X ≲ 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 ∼ 8 and 10 is shown in Fig. B.3 and B.4, respectively. The disfavoured models at z ∼ 8 have ${\overline{x}}_{\mathrm{HII}}\lesssim 0.05(0.78)$, ${\overline{T}}_{\mathrm{K}}\lesssim 10(2180)$ K, $28(2)\lesssim |{\overline{\delta T}}_{\mathrm{b}}|\lesssim 1905(14791)$ mK, $R_{\mathrm{peak}}^{\mathrm{heat}}\lesssim 8(138)h^{-1}\mathrm{Mpc}$ and $\mathrm{\Delta }{R}_{\mathrm{FWHM}}^{\mathrm{heat}}\lesssim 16(213)h^{-1}\mathrm{Mpc}$, while f_heat spans the entire range at 68% (95%) credible intervals for the Single-z analysis case. The disfavoured models at z ∼ 10 in the allUL and Varying A_r case have ${\overline{x}}_{\mathrm{HII}}\lesssim 0.01(0.32)$, ${\overline{T}}_{\mathrm{K}}\lesssim 5(27)$ K, $132(91)\lesssim |{\overline{\delta T}}_{\mathrm{b}}|\lesssim 9120(24888)$ mK, f_heat ≲ 0.03(0.31) $R_{\mathrm{peak}}^{\mathrm{heat}}\lesssim 4(11)h^{-1}\mathrm{Mpc}$ and $\mathrm{\Delta }{R}_{\mathrm{FWHM}}^{\mathrm{heat}}\lesssim 5(20)h^{-1}\mathrm{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 B2. 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 ∼ 8.

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

Fig. C.1. A subset of 1000 power spectra simulated with GRIZZLY at z = 9.1. These power spectra correspond to those with Ar = 0. The red and magenta points refer to the recent LOFAR 1σ 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 TS = TK.

Fig. C.2. Posterior distribution of the source parameters at z = 9.1 for the Varying Ar 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σ and 2σ 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.

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σ upper limit of ${\mathrm{\Delta }}^2(k=0.075h{\mathrm{Mpc}}^{-1},z=9.1,141\mathrm{hours})\approx {(25)}^2{\mathrm{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σ upper limits from Mertens et al. (2025, red) and Acharya et al. (2024c, magenta, Δ_{21 − an}² with error Δ_{21 − 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_HI = 1 and T_S = T_K = 2 K, remains above the 1σ 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 Δ_{21 − an}² − Δ_{21 − 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 adopting for the likelihood of a model surviving at z the expression ℒ_{su, single − z}(θ, z) = 1 − ℒ_{ex, single − z}(θ, 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 ζ ≲ 14, $M_{\min }\lesssim 5\times {10}^{11}{\mathrm{M}}_{\odot }$, $M_{\min ,\mathrm{X}}\lesssim 1.5\times {10}^{10}{\mathrm{M}}_{\odot }$, f_X ≳ 100 and A_r ≲ 2 for the Varying A_r scenario. The limits are similar when we use one upper limit at k = 0.075 h Mpc⁻¹.

Fig. C.3. Posterior distribution of the IGM parameters at z = 9.1 for the Varying Ar 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σ and 2σ 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.

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 ${\overline{x}}_{\mathrm{HII}}\lesssim 0.75$, $27\lesssim {\overline{T}}_{\mathrm{K}}\lesssim 7\times {10}^3$ K, f_heat ≳ 0.62, $4\lesssim {\overline{\delta T}}_{\mathrm{b}}\lesssim 44$ mK, $3\lesssim R_{\mathrm{peak}}^{\mathrm{heat}}\lesssim 160h^{-1}\mathrm{Mpc}$ and $4\lesssim \mathrm{\Delta }{R}_{\mathrm{FWHM}}^{\mathrm{heat}}\lesssim 316h^{-1}\mathrm{Mpc}$. The limits on the IGM parameters for the case where we used the upper limit at k = 0.075 h Mpc⁻¹ only are similar to the case where we used three k−bins upper 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).

All Tables

All Figures

Fig. 1. Set of 1000 power spectra randomly chosen from the 105 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 (Ar = 0). The down arrow points refer to the recent LOFAR 1σ upper limit (Δ212(k, z)±Δ21, err2(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 TS = TK. The coloured power spectra are those with a value larger than Δ212(k, z)−Δ21, err2(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.

In the text

	Fig. 2. Histogram of the disfavoured values of the source parameters for the scenario with Ar = 0. These correspond to the models with power spectra values larger than Δ212(k, z)−Δ21, err2(k, z) in at least one k-bin, and they have a high probability of being ruled out by the recent LOFAR 1σ upper limits (Δ212(k, z)±Δ21, err2(k, z)) from Mertens et al. (2025).
In the text

	Fig. 3. Histogram of the disfavoured values of the IGM parameters for the scenario with Ar = 0. These correspond to the same set of models considered in Fig. 2.
In the text

Fig. 4. Posterior distribution of the disfavoured values of the source parameters for the Ar = 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σ and 2σ 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.

In the text

	Fig. 5. Marginalised probability distribution of the disfavoured values of the source parameters for the Ar = 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.
In the text

Fig. 6. Posterior distribution of the disfavoured values of the IGM parameters for the Ar = 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σ and 2σ 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.

In the text

Fig. 7. Posterior distribution of the source parameters disfavoured values for the Varying Ar 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.

In the text

Fig. 8. Posterior distribution of the IGM parameters of the disfavoured models at z = 9.1 for the Varying Ar 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.

In the text

	Fig. 9. Set of 1000 power spectra randomly chosen from the 105 simulated spectra for Ar = 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 TS = TK at those redshifts. The down arrow points refer to the different 2σ upper limits at those redshifts.
In the text

Fig. 10. Marginalised probability distribution of the disfavoured models’ source parameters at z ∼ 8, 9, and 10 for allUL case. The results are for the Ar = 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.

In the text

	Fig. 11. Marginalised probability distribution of the IGM parameters of the disfavoured models at z ∼ 9 for allUL and Ar = 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.
In the text

	Fig. 12. Same as Fig. 11 but for the Varying Ar scenario.
In the text

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 ℒex, single − z(θ, z = 9.1) = 0.5 and ℒex, joint − z(θ, z = 9.1) = 0.5.

In the text

	Fig. B.1. Same as Fig. B.2 but for z ∼ 8.
In the text

Fig. B.2. Marginalised probability distribution of the disfavoured models’ IGM parameters at z ∼ 10 for the allUL and Ar = 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.

In the text

	Fig. B.3. Same as Fig. B.1 but for the Varying Ar scenario.
In the text

	Fig. B.4. Same as Fig. B.2 but for Varying Ar scenario.
In the text

Fig. C.1. A subset of 1000 power spectra simulated with GRIZZLY at z = 9.1. These power spectra correspond to those with Ar = 0. The red and magenta points refer to the recent LOFAR 1σ 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 TS = TK.

In the text

Fig. C.2. Posterior distribution of the source parameters at z = 9.1 for the Varying Ar 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σ and 2σ 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.

In the text

Fig. C.3. Posterior distribution of the IGM parameters at z = 9.1 for the Varying Ar 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σ and 2σ 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.

In the text