# 21-cm Signal from the Epoch of Reionization: A Machine Learning upgrade to Foreground Removal with Gaussian Process Regression

Anshuman Acharya, Florent Mertens, Benedetta Ciardi, Raghunath Ghara, Léon Koopmans, Sambit Giri, Ian Hothi, Qing-Bo Ma, Garrelt Mellema, Satyapan Munshi

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

To cite this version:

Anshuman Acharya, Florent Mertens, Benedetta Ciardi, Raghunath Ghara, Léon Koopmans, et al. 21-cm Signal from the Epoch of Reionization: A Machine Learning upgrade to Foreground Removal with Gaussian Process Regression. Monthly Notices of the Royal Astronomical Society, 2024, 527 (3), pp.7835-7846. 10.1093/mnras/stad3701. hal-04344205

HAL Id: hal-04344205

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

Submitted on 4 Nov 2024

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# 21-cm signal from the Epoch of Reionization: a machine learning upgrade to foreground removal with Gaussian process regression

Anshuman Acharya<sup>1,★</sup> Florent Mertens<sup>1,2</sup> Benedetta Ciardi,<sup>1</sup> Raghunath Ghara<sup>1,3</sup>

Léon V. E. Koopmans<sup>1,4</sup> Sambit K. Giri<sup>1,5</sup> Ian Hothi,<sup>2,6</sup> Qing-Bo Ma,<sup>7,8</sup> Garrelt Mellema<sup>1,9</sup> and Satyapan Munshi<sup>4</sup>

Accepted 2023 November 23. Received 2023 October 20; in original form 2023 May 17

## ABSTRACT

In recent years, a Gaussian process regression (GPR)-based framework has been developed for foreground mitigation from data collected by the LOw-Frequency ARray (LOFAR), to measure the 21-cm signal power spectrum from the Epoch of Reionization (EoR) and cosmic dawn. However, it has been noted that through this method there can be a significant amount of signal loss if the EoR signal covariance is misestimated. To obtain better covariance models, we propose to use a kernel trained on the GRIZZLY simulations using a Variational Auto-Encoder (VAE)-based algorithm. In this work, we explore the abilities of this machine learning-based kernel (VAE kernel) used with GPR, by testing it on mock signals from a variety of simulations, exploring noise levels corresponding to ≈10 nights (≈141 h) and ≈100 nights (≈1410 h) of observations with LOFAR. Our work suggests the possibility of successful extraction of the 21-cm signal within 2σ uncertainty in most cases using the VAE kernel, with better recovery of both shape and power than with previously used covariance models. We also explore the role of the excess noise component identified in past applications of GPR and additionally analyse the possibility of redshift dependence on the performance of the VAE kernel. The latter allows us to prepare for future LOFAR observations at a range of redshifts, as well as compare with results from other telescopes.

**Key words:** methods: data analysis – techniques: interferometric – dark ages, reionization, first stars – cosmology: observations.

## 1. INTRODUCTION

The Epoch of Reionization (EoR) follows the period of cosmic dawn, when the first stars, galaxies, black holes, and other astrophysical objects formed. These objects began to radiate photons that ionized the neutral gas across the Universe. Understanding this period, which spreads over redshifts z ≈ 5–15, is crucial to learn more about the nature, timing, and mechanism of the formation and evolution of the aforementioned astrophysical objects, as well as their impact on the interstellar medium and intergalactic medium (IGM) surrounding them (Ciardi & Ferrara 2005; Morales & Wyithe 2010; Pritchard & Loeb 2012; Furlanetto 2016; Liu & Shaw 2020). From observations of the Gunn–Peterson troughs of high-z quasars (Becker et al. 2001; Fan et al. 2006) and the optical depth for Thomson scattering of the cosmic microwave

background (CMB) radiation (Planck Collaboration XLVII 2016), we can deduce that most reionization took place in the range 6 ≲ z ≲ 10, with recent observations suggesting an end of reionization at z < 6 (see e.g. Becker et al. 2015 and Bosman et al. 2022).

There exist multiple indirect probes to study this period. For example, the evolution of the observed Lyman-α emitter luminosity function at z > 6 (Clément et al. 2012; Schenker et al. 2013), and Lyman-α absorption profiles to distant quasars (Mortlock 2016; Greig et al. 2017; Davies et al. 2018). However, the most sensitive probe to study the EoR is through the fluctuations of the redshifted 21-cm line of neutral Hydrogen against the CMB (Hogan & Rees 1979; Madau, Meiksin & Rees 1997; Shaver et al. 1999; Tozzi et al. 2000; Ciardi & Madau 2003; Zaroubi 2013). A statistical detection of the strength of these 21-cm brightness temperature fluctuations can allow us to constrain our models of the early Universe and the formation of the first stars and galaxies. For this, a number of interferometric low-frequency radio telescopes

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

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

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

<footnote n="4"><sup>4</sup>Kapteyn Astronomical Institute, University of Groningen, PO Box 800, Groningen NL-9700AV, the Netherlands</footnote>

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

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

<footnote n="7"><sup>7</sup>School of Physics and Electronic Science, Guizhou Normal University, Guiyang 550001, P. R. China</footnote>

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

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

<footnote>*E-mail: anshuman@mpa-garching.mpg.de</footnote>
have been designed, such as LOw-Frequency ARray (LOFAR),<fnref n="1" /> Hydrogen Epoch of Reionization Array (HERA),<fnref n="2" /> Murchison Wide-field Array,<fnref n="3" /> and Precision Array to Probe EoR,<fnref n="4" /> which over the years have been providing increasingly tighter upper limits. For example, The HERA Collaboration et al. (2022) reported $\Delta^2(k = 0.34 h \text{ Mpc}^{-1}) \le 457 \text{ Mpc}^2$ at $z = 7.9$ and $\Delta^2(k = 0.36 h \text{ Mpc}^{-1}) \le 3496 \text{ Mpc}^2$ at $z = 10.4$ from 94 nights of observation, and Mertens et al. (2020, hereafter M20) reported $\Delta^2(k = 0.075 h \text{ Mpc}^{-1}) < 5329 \text{ Mpc}^2$ from 141 h ($\approx 10$ nights) of observation with LOFAR at $z = 9.1$.

One of the major challenges faced in the detection of the 21-cm signal is the fact that it is buried under foregrounds (synchrotron and free-free emissions from the Milky Way and other galaxies) that are several orders of magnitude stronger. To address this issue, The HERA Collaboration et al. (2022) uses the 'foreground avoidance' technique (Kerrigan et al. 2018; Morales et al. 2019) by focusing on regions in Fourier space which are mostly foreground free, while the LOFAR EoR Key Science Project (KSP) team uses foreground modelling and removal, which allows the maximization of scales explored, as well as boosts the sensitivity up to an order of magnitude (Pober et al. 2014).

The most stringent constraints obtained with LOFAR data were presented in M20, where Gaussian process regression (GPR, as described in Mertens, Ghosh & Koopmans 2018; Gehlot et al. 2019; Hothi et al. 2021) was used for hyperparameter optimization with different Matern-class functions (equation 6, see below) chosen as covariance kernels for modelling different components of the observed data and then recovering the fitted data cube. However, Kern & Liu (2021) pointed out some issues with this approach. Primarily, they found that given the choice of normalization and bias correction in the power-spectra estimation used in M20, misestimation of the covariance kernel for the EoR signal could lead to significant signal loss. This can have severe ramifications on the astrophysical interpretations of the estimated 21-cm signal power spectrum. Further, they show that alternative choices for normalization and weighting schemes could reduce the dependence on the choice of covariance priors, thus reducing its impact on the estimation of the 21-cm signal. However, here we focus on improving the covariance prior, while in future works we plan to explore other normalization and bias correction schemes to further upgrade the overall analysis pipeline.

To improve the covariance kernel, we refer to Mertens, Bobin & Carucci (2023). They propose a machine learning (ML)-based approach to GPR, where the covariance kernel for the 21-cm signal is obtained by implementing a ML-based algorithm that learns from simulations. The results obtained can then be compared against runs of the same simulation code, to constrain the physical parameters used in it. As the covariance kernel is trained over a range of physical parameters, this would significantly reduce chances of misestimation, and thus it can be reliably used to derive astrophysical parameters necessary for the same simulation code to generate similar power spectra.

In this paper, we use GRIZZLY (Ghara, Choudhury & Datta 2015; Ghara et al. 2018, 2020) for generating the training, test, and validation data sets. This code has been employed previously (see Ghara et al. 2020) to constrain astrophysical parameters based on

the results obtained in M20. As this combines N-body simulations with one-dimensional (1D) radiative transfer, it is more physically precise than seminumerical algorithms, while not being as computationally expensive as codes that use 3D radiative transfer (while still performing reasonably well, as shown in Ghara et al. 2018). Here, we test the performance of this ML-based kernel versus covariance kernels used with GPR in previous work. In Section 2, we discuss a range of simulations used to generate mock 21-cm data sets, as well as introduce the ML-trained 21-cm kernel. We also provide a short introduction to GPR. We report the results and comparisons between kernel performances in Section 3. Finally, we discuss the role of the excess noise component found in M20 and the overall performance of the ML-based kernel in Section 4. In a companion paper, we will apply the new pipeline to 10 nights of LOFAR data at $z = 9.1$, as was done in M20.

## 2. METHODOLOGY

In this section, we introduce the pipeline used to implement GPR to recover the 21-cm signal from mock data sets, comparing the performance of an ML-based kernel with kernels used in M20.

### 2.1 Simulations of the 21-cm signal

The 21-cm differential brightness temperature relative to the CMB for any patch of the IGM is given by (see Furlanetto, Oh & Pierpaoli 2006)

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

where $x_{H1}$ is the fraction of neutral hydrogen, $\delta_B$ is the fractional overdensity of baryons, $T_S$ is the spin temperature, $T_{CMB}$ is the temperature of the CMB photons at that redshift, $\Omega_m$ is the total matter density, $\Omega_B$ is the baryon density, $z$ is the redshift, and $h$ is the Hubble constant in units of $100 \text{ km s}^{-1} \text{ Mpc}^{-1}$. In this equation, the parameters affecting the large-scale fluctuations of $\delta T_b$ are $x_{H1}$, $T_S$ and $\delta_B$.

We consider a variety of simulations to generate mock 21-cm differential brightness temperature maps as discussed below. In Section 2.1.1, we employ maps generated using GRIZZLY, where we focus on variations tied to fluctuations in $x_{H1}$ and $T_S$, while assuming that the fluctuations due to $\delta_B$ are small and can thus be ignored. In Section 2.1.2, we do not make this assumption and employ maps generated using the reionization simulation code CRASH. In Appendix A, we also consider the additional case of using 21CMFAST (Mesinger & Furlanetto 2007; Greig & Mesinger 2015) to generate the 21-cm differential brightness temperature maps.

#### 2.1.1 GRIZZLY simulations

GRIZZLY (Ghara, Choudhury & Datta 2015; Ghara et al. 2018, 2020) employs a 1D radiative transfer scheme in combination with cosmological density fields and halo catalogues obtained from an N-body simulation to produce brightness temperature maps of the 21-cm signal at different redshifts for a given source model. The algorithm has been shown to reproduce results similar to those obtained with 3D radiative transfer schemes with the same N-body simulation, while being at least $10^5$ times faster (Ghara et al. 2018). Because of this, we can run a large number of GRIZZLY simulations without

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

<footnote n="2">²https://reionization.org/</footnote>

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

<footnote n="4">⁴http://eor.berkeley.edu</footnote>
the process being too computationally expensive. Furthermore, it has a wide range of physical parameters that can be varied, thus allowing us to explore the role of different physical processes in generating different models of the 21-cm signal. The density fields, velocity fields, and the halo lists used in this work are obtained from the same N-body simulation (500 h⁻¹ cMpc box length, 6912³ particles, with a mass resolution of 4.05 × 10⁷ M☉) which was used in Ghara et al. (2020). In this study, we consider two major GRIZZLY models presented in Sections 3.1 and 3.2 of Ghara et al. (2020). Similar to their implementation, we use four physical parameters to model the sources: the ionization efficiency (ζ), the minimum mass of the ultraviolet (UV) emitting haloes (Mmin), the minimum mass of the X-ray emitting haloes (MminX), and the X-ray heating efficiency (fX). The emission rate of ionizing photons and X-rays per unit stellar mass from a halo are ζ × 2.85 × 10⁴⁵ and fX × 10⁴² s⁻¹ M☉⁻¹, respectively. Further, the X-ray spectral index α is fixed at 1.2, as done in Ghara et al. (2020). Lastly, they note that all other IGM properties can be derived from these parameters, and thus using just these to define the simulation is sufficient. The properties of the two models adopted are listed below:

(i) **xH1 fluctuation dominated model:** Here, we assume a uniform Ly α background strong enough to allow the spin temperature T_S to be fully coupled to the gas temperature T_K. Further, we adopt the following parameters: ζ = 7.0, Mmin = MminX = 10⁹ M☉, and fX = 100, which makes the gas temperature (and in turn, T_S) significantly high compared with T_CMB due to strong X-ray heating. This assumption of T_S ≫ T_CMB ensures that δT_b becomes insensitive to the (1 - T_CMB/T_S) term from equation (1). Thus, all variability of δT_b is tied to the fluctuation of the neutral hydrogen fraction xH1.

(ii) **T_S fluctuation dominated model:** In this case, while we continue to have the assumption of a strong, uniform Ly α background to ensure coupling of T_S and T_K, we change our parameters to relax the condition of T_S ≫ T_CMB. This is done by reducing the X-ray heating and ionization efficiency. Thus, we adopt the following parameters: ζ = 3.0, Mmin = 10⁹ M☉, MminX = 10¹⁰ M☉, and fX = 1. This allows for greater variability tied to T_S, with regions of partial reionization and heating forming in the IGM.

### 2.1.2 CRASH simulations

As a reference, we also use the simulations of reionization described in Eide et al. (2018, 2020) and Ma et al. (2021). These are obtained by post-processing the large-scale, high-resolution hydrodynamical simulation MASSIVE BLACK-II (Khandai et al. 2015; box length 100 h⁻¹ cMpc, 2 × 1792³ gas and dark matter particles, corresponding to a resolution of 2.2 × 10⁶ and 1.1 × 10⁷ h⁻¹ M☉, respectively) with the multifrequency 3D radiative transfer code CRASH (Ciardi et al. 2001; Maselli, Ciardi & Kanekar 2009; Graziani, Maselli & Ciardi 2013; Graziani, Ciardi & Glatzle 2018; Glatzle, Ciardi & Graziani 2019). Here, we make use of the 'Stars' simulation run (which includes only stellar type sources) to generate the mock 21-cm signal data at z = 9.18. We refer the reader to the original papers for more detailed information on the simulations.

## 2.2 Gaussian process regression

GPR (Rasmussen & Williams 2006; Aigrain & Foreman-Mackey 2023) can be used to model a noisy observation y = f(x) + ε, with ε Gaussian noise having variance σ²noise. This is achieved by modelling the Gaussian process as a joint probability distribution for y = {yᵢ}ᵢ=1,...,N, as f = f(x), which is fully defined by its mean

vector (m) and covariance matrix (K, also called covariance 'kernel') as

$$f \sim \mathcal{N}(m(\mathbf{x}), K(\mathbf{x}, \mathbf{x})), \quad (2)$$

for a set of points x (independent parameters). Here, the covariance matrix K gives the covariance between the function values at any two points and can be written as Kᵢⱼ = κ(xᵢ, xⱼ, φ) + δᵢⱼσᵢ², where κ(xᵢ, xⱼ, φ) can be optimized by the choice of hyperparameters represented by φ, and δᵢⱼ is the Kronecker-delta function.

When applying it to radio data to extract the 21-cm signal, we split this function into a foreground component, fᵢₗ, and the 21-cm signal, fᵢ₁, giving

$$y = f_{\mathrm{fg}}(x) + f_{\mathrm{21}}(x) + \epsilon \quad (3)$$

where y is the observed data and x is the frequency. Next, following M20, we further split the foreground component into the intrinsic sky emission component (fsky), which comes from the confusion-limited extragalactic sources and from the Milky Way, and the mode-mixing contaminants component fmix, which has contributions from the instrument chromaticity and calibration errors. Beyond the foreground, we also model the noise (represented by ε in equation 3) using estimates of the noise variance for ≈10 nights of observation from M20. In addition to this, M20 found a significant spectrally-correlated residual, and thus we inject this 'excess noise' component (fex) into our model as well. This gives an updated version of equation 3:

$$y = f_{\mathrm{sky}}(x) + f_{\mathrm{mix}}(x) + f_{\mathrm{ex}}(x) + f_{\mathrm{21}}(x) + \epsilon. \quad (4)$$

For the sake of simplicity and utilizing the additive property of matrices, we can reduce equation (4) to y = f(x) + ε, and represent f(x) with its corresponding covariance kernel K (from equation 2) as

$$K = K_{\mathrm{sky}} + K_{\mathrm{mix}} + K_{\mathrm{noise}} + K_{\mathrm{ex}} + K_{\mathrm{21}}. \quad (5)$$

M20 modelled each of these kernels using the best possible fit Matern-class functions (equation 6, Stein 1999):

$$K_{\mathrm{Matern}}(r) = \sigma^2 \frac{2^{1-\eta}}{\Gamma(\eta)} \left(\frac{\sqrt{2\eta r}}{l}\right)^{\eta} \kappa_{\eta} \left(\frac{\sqrt{2\eta r}}{l}\right). \quad (6)$$

Note that in the Matern-class function, r is the absolute difference between the frequencies of two sub-bands, κη is the modified Bessel function of the second kind, and Γ is the Gamma-function.

M20 obtained the best possible fit Matern-class function by taking different values of the hyperparameter η, maximizing the marginal likelihood (also known as the evidence), and obtaining estimates for the coherence-scale hyperparameter l and the variance σ². Then, for each kernel, they chose the η that led to the highest evidence by calculating the analytical integral over f which is the log-marginal-likelihood (LML, see Section 2.3 in Mertens, Ghosh & Koopmans 2018) and choosing the kernel that maximizes its value. For calculating the hyperparameters (listed in the second column of Table 1), M20 used a gradient-descent-based optimization algorithm for maximizing the LML.

## 2.3 ML-trained 21-cm kernel

The limitations of GPR as pointed out by Kern & Liu (2021) mainly boil down to the choice of covariance kernel for the 21-cm signal. While the choice of hyperparameters allows a variety of functions to be accessed, the same function might not work equally well across the k-space. Having a function obtained by employing ML trained on power spectra of simulations where the sources of ionization
<table><caption>**Table 1.** We use the results of M20 as our input (2nd column) for the hyperparameters (1st column) to generate a mock data set composed of the $f_{sky}$, $f_{mix}$, and $f_{ex}$ components, in addition to the $x_{H1}$ ($T_S$) fluctuation dominated models from GRIZZLY to generate the 21-cm signal component at $z = 9.16$. For the recovery using GPR, we keep the value of $\eta$ fixed for the $f_{sky}$, $f_{mix}$, and $f_{ex}$ components, and estimate the median and the 68 per cent confidence intervals for the coherence scales and variances from the MCMC-based LML maximization (4th to 7th columns) using the priors as in M20 (3rd column), i.e. we either use linear scale uniform priors $U$, or logarithmic scale indicated as $-$. As the latter is over several orders of magnitude, it has not been listed and can be assumed to be an uninformed prior from $-\infty$ to $+\infty$. For $\approx 10$ nights of observation, we use $\sigma_{\text{noise}}^2 = 74 \times 10^3 \text{ mK}^2$ (as obtained in M20). For $\approx 100$ nights of observation, we scale it down by a factor of 2 for the $f_{sky}$ component and by a factor of 20 for the rest of the components (see Section 2.4). The values for $\sigma_{\text{noise}}^2$ in each case are listed in the first row for each component. Note that $\sigma_{\text{noise}}^2$ has been scaled to the Field of View of the GRIZZLY simulations, which is $3.03^{\circ} \times 3.03^{\circ}$.</caption><thead><tr><th rowspan="2">Hyperparameter</th><th rowspan="2">Input value from M20</th><th rowspan="2">Prior</th><th colspan="2">x<sub>H1</sub> fluctuation dominated</th><th colspan="2">T<sub>S</sub> fluctuation dominated</th></tr><tr><th>10 nights</th><th>100 nights</th><th>10 nights</th><th>100 nights</th></tr></thead><tbody><tr><td>σ<sup>2</sup><sub>noise</sub> (× 10<sup>3</sup> mK<sup>2</sup>)</td><td></td><td></td><td>74</td><td>37</td><td>74</td><td>37</td></tr><tr><td>η<sub>sky</sub></td><td>+∞</td><td></td><td>+∞</td><td>+∞</td><td>+∞</td><td>+∞</td></tr><tr><td>l<sub>sky</sub></td><td>47.5</td><td>U(10, 100)</td><td>43.44 ± 4.02</td><td>47.76 ± 1.64</td><td>52.00 ± 6.91</td><td>45.84 ± 1.67</td></tr><tr><td>σ<sup>2</sup><sub>sky</sub>/σ<sup>2</sup><sub>noise</sub></td><td>611</td><td>—</td><td>665.17 ± 29.52</td><td>579.91 ± 24.11</td><td>598.99 ± 27.13</td><td>633.43 ± 28.23</td></tr><tr><td>σ<sup>2</sup><sub>noise</sub> (× 10<sup>3</sup> mK<sup>2</sup>)</td><td></td><td></td><td>74</td><td>3.7</td><td>74</td><td>3.7</td></tr><tr><td>η<sub>mix</sub></td><td>3/2</td><td></td><td>3/2</td><td>3/2</td><td>3/2</td><td>3/2</td></tr><tr><td>l<sub>mix</sub></td><td>2.97</td><td>U(1, 10)</td><td>2.99 ± 0.07</td><td>2.97 ± 0.08</td><td>3.08 ± 0.08</td><td>2.90 ± 0.08</td></tr><tr><td>σ<sup>2</sup><sub>mix</sub>/σ<sup>2</sup><sub>noise</sub></td><td>50.4</td><td>—</td><td>54.80<sup>+2.53</sup><sub>-2.42</sub></td><td>46.05<sup>+2.14</sup><sub>-2.04</sub></td><td>55.82<sup>+2.64</sup><sub>-2.52</sub></td><td>52.12<sup>+2.58</sup><sub>-2.46</sub></td></tr><tr><td>σ<sup>2</sup><sub>noise</sub> (× 10<sup>3</sup> mK<sup>2</sup>)</td><td></td><td></td><td>74</td><td>3.7</td><td>74</td><td>3.7</td></tr><tr><td>η<sub>ex</sub></td><td>5/2</td><td></td><td>5/2</td><td>5/2</td><td>5/2</td><td>5/2</td></tr><tr><td>l<sub>ex</sub></td><td>0.26</td><td>U(0.1, 0.8)</td><td>0.26 ± 0.003</td><td>0.26 ± 0.005</td><td>0.27 ± 0.01</td><td>0.26 ± 0.01</td></tr><tr><td>σ<sup>2</sup><sub>ex</sub>/σ<sup>2</sup><sub>noise</sub></td><td>2.18</td><td>—</td><td>2.15<sup>+0.05</sup><sub>-0.05</sub></td><td>1.91<sup>+0.11</sup><sub>-0.10</sub></td><td>2.25<sup>+0.07</sup><sub>-0.07</sub></td><td>2.32<sup>+0.18</sup><sub>-0.17</sub></td></tr><tr><td>σ<sup>2</sup><sub>noise</sub> (× 10<sup>3</sup> mK<sup>2</sup>)</td><td></td><td></td><td>74</td><td>3.7</td><td>74</td><td>3.7</td></tr><tr><td>x<sub>1</sub></td><td>—</td><td>—</td><td>-0.10 ± 0.99</td><td>-0.23 ± 1.02</td><td>-0.41 ± 0.92</td><td>-0.07 ± 1.00</td></tr><tr><td>x<sub>2</sub></td><td>—</td><td>—</td><td>-0.17 ± 0.98</td><td>0.42 ± 0.77</td><td>0.92 ± 0.59</td><td>-0.06 ± 0.31</td></tr><tr><td>σ<sup>2</sup><sub>21</sub>/σ<sup>2</sup><sub>noise</sub></td><td>—</td><td>—</td><td>0.04<sup>+0.21</sup><sub>-0.03</sub></td><td>0.17<sup>+0.37</sup><sub>-0.12</sub></td><td>0.25<sup>+0.33</sup><sub>-0.14</sub></td><td>1.70<sup>+1.07</sup><sub>-0.65</sub></td></tr></tbody></table>

are modelled using parameters that sample a wide range of values, allows greater flexibility, and reduces chances of misestimation. Further, it allows for a direct comparison with physical quantities, as we can reliably derive the source parameters necessary to generate simulations that produce the power spectra estimated by the ML-based kernel.

To achieve this, Mertens, Bobin & Carucci (2023) use a Variational Auto-Encoder (VAE, Kingma & Welling 2013, 2019) algorithm. Simply put, an Auto-Encoder (AE, Goodfellow, Bengio & Courville 2016) is an unsupervised neural network which compresses data by reducing the number of independent parameters used to describe it into what is referred to as a 'latent space' of hyperparameters. Thus, it is primarily used for data compression by filtering out independent parameters that are deemed to be unnecessary because they only slightly affect data recovery. This is a two steps process, where the first step of reducing the number of independent parameters into the latent space is called *encoding*, while the step of recovering the data given the latent space parameters is called *decoding*. Instead of taking an input of just a set of parameters $a_1, \dots, a_n$, a VAE (Pinheiro Cinelli et al. 2021) uses probability distributions of each parameter, thus allowing to interpolate in the latent space, and to generate a large range of new samples of reconstructed data (in our case 21-cm signal models), which are not limited by the data that the encoder was trained on. However, this also means that the VAE has an inherently larger error than an AE.

Thus, reconstructing the training data using the decoder with the latent space generated by the VAE as input would not be an exact match to the original training data. However, the VAE is designed to optimize a trade-off between reconstruction accuracy and the fidelity of the latent space representation, by minimizing the Kullback-Leibler (KL) divergence loss (Kullback & Leibler 1951), which is a measure of the divergence between the distribution of reconstructed data and the training data. So while the reconstructed

data and training data might not be an exact match, if their overall distribution is similar (i.e. divergence is minimized), the training is considered successful. Thus, in the training of a VAE algorithm, the following sources of error exist:

(i) **Encoder**: the error due to sampling from a distribution for each independent parameter to build a latent space. Sampling from a distribution is expected to be noisier than choosing point values.

(ii) **Decoder**: the error due to deriving the independent parameters, given some value of the latent space parameters. As this does not involve sampling from distributions, the contribution to the total error is expected to be smaller.

Mertens et al. (2023) show that GPR can be used to estimate the values of the latent space parameters, after which the decoder of the VAE kernel can be used to estimate the independent parameters, and in turn to obtain the recovered 21-cm signal. In this case we would then use only the decoder part of the VAE algorithm. However, as minimizing KL divergence loss requires training the encoder and decoder together, we proceed as follows.

We start by using two hyperparameters, $x_1$ and $x_2$, and an associated variance, as the latent space parameters. We train the VAE algorithm on a data set of $\approx 1500$ simulations (the training set), and validate it against an independent data set of $\approx 150$ simulations (the validation set). We refer to this as the *VAE kernel*. The training and validation data sets are generated by running GRIZZLY simulations with a wide range of values for the parameters introduced in Section 2.1.1. We sample them randomly at the same redshift as the targeted 21-cm signal ($z = 9.16, 8.30, \text{and } 10.11$) and in the following ranges: $\zeta = [-3, 3]$, $M_{\text{min}} = [9, 12]$, $M_{\text{minX}} = [9, 12]$, and $f_X = [-3, 2]$ in the log space. We choose these ranges to be significantly broader than necessary (i.e. the performance of the VAE kernel does not show any appreciable difference even if they are reduced by
multiple orders of magnitude), and also note that the performance of the kernel is not impacted if we fix one of the parameters to a specific value and vary the remaining three. This ensures that the sampling range does not induce any major bias. Next, to train the VAE we use 2000 iterations, as we find that the KL divergence loss and the reconstruction loss for both the training and validation data sets stabilize after ≈500 iterations. The reconstruction loss is defined as the total error made when using the encoder to obtain specific values of the latent space parameters, and then employing the decoder to retrieve the data from those parameters. We then evaluate the ratio between the output and input data as a function of the wave number $k$, and finally calculate the median ratio, which is ≈1 across all $k$-bins. However, as discussed earlier, the 68 per cent confidence interval is significant, being ≈10 per cent for $k = [0.06, 0.12]$ $h$ Mpc$^{-1}$, ≈35 per cent for $k = [0.12, 0.43]$ $h$ Mpc$^{-1}$, and ≈27 per cent for $k = [0.43, 1.11]$ $h$ Mpc$^{-1}$ at $z = 8.30, 9.16$, and 10.11.

We note a similar reconstruction error for training sets of sizes ranging from 1000 to 5000 simulations (with the validation set scaling as ~10 per cent of the training set in size), while the errors become worse when using less than 1000 simulations. We also tried to used three hyperparameters, but saw no significant improvement in the recovery error. We thus choose to use two hyperparameters to avoid overfitting.

A high reconstruction error was to be expected, as it includes also the error due to the encoding process. While this is not required for our purpose, we do need to evaluate the exact contribution from the decoder. For this, we create a testing set of ≈150 simulations, and ensure that it also includes 'extreme' models from GRIZZLY (which we define as cases where the power spectrum differs by at least an order of magnitude from the mean of the power spectra from the training and validation data sets), along with some injected stochastic noise. We then apply GPR for hyperparameter optimization (using MCMC, as discussed in Section 2.5) to estimate $x_1, x_2$, and the associated variance. This is then used with the decoder to obtain the recovered signal. The median ratio of the recovered and input data across all $k$-bins is again ≈1. But now we obtain a 68 per cent confidence interval of ≈0.5 per cent for $k = [0.06, 0.12]$ $h$ Mpc$^{-1}$, ≈1 per cent for $k = [0.12, 0.43]$ $h$ Mpc$^{-1}$, and ≈3 per cent for $k = [0.43, 1.11]$ $h$ Mpc$^{-1}$ at $z = 8.28, 9.16$, and 10.11, respectively. By comparing with the total error estimated above, we note that even where decoder error is high, i.e. for the highest $k$-bins, it is still a minor contributor. Therefore, as the overall error is ≲5 per cent, we accept the VAE kernel as a reliable ML-based kernel for the 21-cm signal, and proceed to use it with GPR for signal recovery.

## 2.4 Generating mock data sets

To build a mock data set in the gridded $U-V$ cube domain of radio observations, we derive a full data set $y$ by adding each term on the right-hand side of equation (4). For this, we adopt the values of $\eta, l$, and $\sigma^2$ for $K_{sky}, K_{mix}$, and $K_{ex}$ given in the second column of Table 1 to generate their corresponding $f_{sky}, f_{mix}$, and $f_{ex}$. These values are obtained from the results of M20, where they used $\eta_{sky} = +\infty$ for the intrinsic sky, $\eta_{mix} = 3/2$ for the mode-mixing contaminants, and $\eta_{ex} = 5/2$ for the excess noise to fix Matern-class functions for each of these components, and then employed GPR to obtain the coherence-scale hyperparameter $l$ and its associated variance by adopting $\eta_{21} = 1/2$ for the 21-cm signal. Note that while $f_{sky}$ and $f_{mix}$ are not independent quantities, their recovery with GPR treats them as such. In this paper, we also generate them as independent components to build our mock data sets (unlike in real data, where the mode-mixing component depends on the true sky signal). Thus, the overall quality of the

recovered $f_{sky}$ and $f_{mix}$ is better from the mock data sets, than from real data. However, the effect on the accuracy of recovery is not expected to be severe, as the impact of not factoring in the inter-dependency is insignificant, as compared with the overall power of these components. That is, even without assuming their inter-dependence, GPR can recover them with reasonably high accuracy.

Next, we also assume the value $\sigma^2_{noise} \sim 74 \times 10^3$ mK$^2$ (from M20) to simulate the noise component $\epsilon$ for ≈10 nights of observation with LOFAR scaled to the Field of View corresponding to the GRIZZLY simulations, which is 3.03° × 3.03°. The variance for the intrinsic sky-emission, mode-mixing, and excess noise components scales by the noise variance as indicated in Table 1.

We also consider the case of ≈100 nights of observation, to provide estimates for future observations of similar duration with LOFAR. In this case, the noise variance is expected to be reduced by a factor of 10, assuming all effects to scale uniformly from ≈10 to ≈100 nights of data. However, it should be noted that the data used in M20 was plagued with issues such as Radio Frequency Interference (RFI) flagging and bad ionospheric conditions, and it would be of little scientific value to assume the same conditions to persist for a 10 times longer observational duration. For this reason, we assume ideal conditions (e.g. picking observation nights with good ionospheric conditions) and assume that the thermal noise is reduced by a factor of 20 instead, obtaining $\sigma^2_{noise} \sim 3.7 \times 10^3$ mK$^2$. As we expect the excess noise to behave in a similar fashion, we also reduce its variance $\sigma^2_{ex}$ by a factor of 20. M20 noted that the mode-mixing contaminants were not decreasing when integrating over more nights of data. However, the parts of this component due to effects of ionosphere and calibration errors are uncorrelated from night to night, and thus are expected to decrease with longer integration, leading to a reduction of the mode-mixing term. Additionally, this should also decrease because of the improved UV-coverage. Here, we thus consider an ideal scenario to understand how the performance of the VAE kernel compares with previously used kernels when also the mode-mixing term's variance $\sigma^2_{mix}$ is reduced by a factor of 20. However, we scale down the intrinsic sky component only by a factor of 2. This is justified because (i) the confusion noise limit is set by the resolution of the LOFAR telescope, so that unresolved point sources cannot be subtracted even when integrating over more nights; (ii) modelling an increasing number of point sources is not a trivial task, and doing it accurately does not seem feasible in the short-term. Thus, we limit ourselves to an assumption of modest improvements. Because of this, we assume that the variance of $f_{sky}$ is scaled by a thermal noise variance of $\sigma^2_{noise} \sim 37 \times 10^3$ mK$^2$, corresponding to a factor of 2 reduction from the noise variance in the ≈10 nights case.

Following the procedure outlined above, we generate two sets of foreground, noise, and excess noise components, i.e. for ≈10 and ≈100 nights of observation. We then add the 21-cm signal (see Section 2.1) to obtain a complete mock data set.

The input power spectra for the overall mock data set using GRIZZLY is shown in Figs 1 and 2 for the $x_{H1}$ and $T_{S}$ fluctuation dominated case, respectively. The ratio between the 21-cm signal and the $1\sigma$ uncertainty of the noise power spectrum is a measure of the signal-to-noise ratio (SNR). Thus, for $k$-bins where the 21-cm signal's power is greater than the $1\sigma$ error of the noise, there is a chance of detectability with $1\sigma$ confidence. From Fig. 1, we can thus conclude that for ≈10 nights of observation a detection of the $x_{H1}$ fluctuation dominated 21-cm signal is unlikely, as the SNR is ≪1 for all $k$-bins. However, the chances of detectability improve for ≈100 nights of observation, as SNR ≈ 1. In Fig. 2, we see that detecting the $T_{S}$ fluctuation dominated signal should be possible also for ≈10
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Power spectrum of the mock data set (solid line) generated using the $x_{H1}$ fluctuation dominated model from GRIZZLY. It consists of the foreground component (intrinsic sky + mode-mixing contaminants; topmost dashed), the excess noise (dashed-dotted), the noise (middle dashed), and the 21-cm signal at $z = 9.16$ (bottom-most thick dashed). We also show the $1\sigma$ upper limit (dotted) achievable if the data set is thermal noise dominated, i.e. any signal below this line has SNR < 1. The top (bottom) panel refers to a case with the noise corresponding to 10 (100) nights of observation.</figcaption></figure>

nights of observations, as SNR is >1 for all k-bins, while for ≈100 nights it is ≳ 10, assuring a detection as long as the covariance kernel used to model the 21-cm signal is correctly estimated.

While we first focus on the 21-cm signal at $z = 9.16$ in order to make a direct comparisons with M20, we also consider $z = 8.30$ and 10.11 to prepare for future LOFAR results at these redshifts.

## 2.5 Recovery with MCMC

As discussed in Section 2.2, M20 used a gradient-descent method to maximize the LML. However, in this paper we instead use MCMC sampling (Foreman-Mackey et al. 2013) to estimate the hyperparameters by sampling their posterior distributions. The benefit of MCMC sampling is that it allows us to also have a measurement of the uncertainty on the hyperparameters, which can then be propagated. We build the posterior distributions adopting flat uniform priors with broad ranges as used in M20. For the coherence scales, we provide a smaller range for the uniform prior (U) to improve convergence time, as done by M20. Thus, as listed in the third column of Table 1, we use a range of (10, 100) for $l_{sky}$ and of (1, 10) for $l_{mix}$. To ensure that the converged value for $l_{ex}$ remains in the $1\sigma$ confidence interval, we used a range of (0.1, 0.8) rather than of (0.2, 0.8) as in M20. While we note that a Gamma prior for the variances of the different component leads to faster convergence, we still adopt flat priors in the logarithmic scale over several orders of magnitude (thus effectively an uninformed prior) to minimize the chances of bias.

In addition to this, we note that the coherence scale and the variance for the 21-cm signal are dependent on the baseline length. So, while to recover the intrinsic sky, mode mixing and excess noise components we continue to employ the same Matern-class functions used to

<figure><img src="image_4.png" /><figcaption>**Figure 2.** As Fig. 1 for the $T_S$ fluctuation dominated model.</figcaption></figure>

generate them (i.e. with the same $\eta$ values of the input, see Table 1), when using Matern-class functions to recover the 21-cm signal, we modify the hyperparameter $l$ and the variance as

$$l = \frac{l_0}{1 + 0.001l_{\alpha}l_0(u - u_{\min})} \quad \text{and} \quad \sigma^2 = \sigma_0^2\sigma_{\text{norm}}^2 \left(\frac{u}{u_{\min}}\right)^{\sigma_{\alpha}^2}, \quad (7)$$

where $l_0$ and $\sigma_0^2$ are the coherence-scale parameter and variance used in M20 for baseline length $u$ (where $u_{min}$ is the minimum baseline length), but now we introduce the additional parameters $l_\alpha$ and $\sigma_\alpha^2$ to fully define the coherence-scale hyperparameter and the variance. Further, $\sigma_{\text{norm}}^2$ is chosen such that the mean of the variance over all baselines is $\sigma_0^2$. The two additional parameters (i.e. $l_\alpha$ and $\sigma_\alpha^2$) thus allow us to encode the dependence on the baseline length into the covariance kernel for the 21-cm signal.

To recover the 21-cm signal component, we use Matern-class functions with two specific values of $\eta$: the Exponential Matern-class function with $\eta = 1/2$ (which, as shown in M20, maximizes the evidence), and the Matern32 function with $\eta = 3/2$. We then compare their performance for recovery with the VAE kernel using GPR. For the hyperparameters $x_1$ and $x_2$, we again take an uninformed flat prior in the linear space.

## 3. RESULTS

Here, we discuss a qualitative comparison between the results which are shown in Figs 3 and 4 using the three aforementioned kernels (Exponential, Matern32, and VAE). We then analyse the estimated coherence-scale hyperparameter and variance values for each of the components of the mock data sets defined in Section 2.1.1 in Table 1. Further, we qualitatively compare the results obtained by using the full simulations of reionization (Section 2.1.2) rather than the mock data sets generated with GRIZZLY. Lastly, we explore the role of redshift on the performance of the VAE kernel, by testing cases at $z = 8.30$ and 10.11, and by comparing them against the results obtained for $z = 9.16$.
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Power spectrum recovered using the Exponential (dash-dotted line) and Matern32 (dotted) Matern-class functions-based covariance kernels, and the VAE-based kernel (solid), together with the $x_{H1}$ fluctuation dominated model 21-cm signal (darker dashed) and noise (lighter dashed) at $z = 9.16$. The $2\sigma$ uncertainty on the recovered signal for each kernel is shown as a shaded area in the corresponding colours. The top (bottom) panel refers to a case with the noise corresponding to 10 (100) nights of observation. We also plot the estimated upper limits of power at $k = 0.075 \ h\ Mpc^{-1}$ from M20 (cross). Note that this value can be significantly higher than the signal due to more complex noise in real data.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 4.** As Fig. 3 for the $T_S$ fluctuation dominated model.</figcaption></figure>

## 3.1 $x_{H1}$ fluctuation dominated model

In Fig. 3, we compare the results from the three kernels (VAE in blue solid, Exponential in orange dash-dotted, and Matern32 in green dotted) in recovering the power spectrum of the 21-cm signal at $z = 9.16$ (in grey dashed). We also show the $2\sigma$ uncertainty on the recovery from different kernels to compare their performance. We note that if the lower bound of the $2\sigma$ uncertainty of the recovered signal is below the uncertainty on the thermal noise (as shown in Fig. 1) the recovery qualifies as an upper limit, otherwise it is referred to as a detection. Based on this, we note that in this case, the recovery from all kernels is going to provide upper limits, as the thermal noise uncertainty is higher than the 21-cm signal for ≈10 nights, and comparable with it for ≈100 nights.

We note that for ≈10 nights of observations, while the VAE kernel has uncertainty bands wider than the Matern-class function-based kernels, the input 21-cm signal is contained within its constrained region. Thus, while the recovered signal for the three kernels is comparable, the VAE kernel is robust enough to compensate for the overestimation and to contain the signal within the $2\sigma$ limits of the error, although it is still an upper limit rather than a constrained detection. As discussed in Section 2.4, the reason for this non-detection when using Matern-class functions-based kernels and broad uncertainty bands for the VAE kernel is due to the low SNR, which is <0.1 across all k-bins. This, however, improves to an average SNR of ≈1 for ≈100 nights of observations, for which, as expected, we obtain tighter constraints and an improved prediction of the actual value. We see, though, that the recovered signal from Matern-

class function-based kernels is still overpredicted. In particular, for $k > 0.5 \ h\ Mpc^{-1}$ the Exponential and Matern32 kernels are unable to contain the signal even in their $2\sigma$ uncertainty bands. On the other hand, the VAE kernel contains the signal in its $2\sigma$ uncertainty bands across all k-bins, despite the recovered signal being about an order of magnitude higher than the input signal. The VAE kernel also does a much better job in recovering the overall shape of the power spectrum. By comparing to the estimated power at $k = 0.075 \ h\ Mpc^{-1}$ from M20 (cross), it is clear that in this case the VAE kernel is also capable of significantly improving the 21-cm upper limits estimate. However, we highlight that the VAE kernel applied to real data is still likely to provide upper limits higher than what has been shown here, because of the more complex noise, and thus the improvement provided by the kernel might be lesser.

In the 4th and 5th columns of Table 1, we show the MCMC estimates for the coherence-scale hyperparameter and the variance obtained by applying GPR to the input power spectrum of the mock data set (in indigo, Fig. 1). The covariance kernels for the intrinsic sky, mode mixing, and excess noise components are the same as those used to generate them (i.e. the value of $\eta$ is fixed), while we adopt the VAE kernel for the 21-cm signal. Note that the variances for all components are scaled down by the corresponding value of $\sigma_{\text{noise}}^2$. This is equal to $74 \times 10^3 \ \text{mK}^2$ for all components for ≈10 nights of observations (see M20), while for ≈100 nights this corresponds to scaling down by a factor of 2 (i.e. $\sigma_{\text{noise}}^2 = 37 \times 10^3 \ \text{mK}^2$) for the intrinsic sky component, and by a factor of 20 (i.e. $\sigma_{\text{noise}}^2 = 3.7 \times 10^3 \ \text{mK}^2$) for other components (as discussed in Section 2.4). From the MCMC estimates, we note that the measurement of the coherence scale for the $f_{\text{sky}}$ and $f_{\text{mix}}$ improves from ≈10 to ≈100 nights of observation. However, the variance estimates do not show an improvement, and even slightly worsen for the excess noise component. The estimates of $x_1$ and $x_2$ and their associated variance $\sigma_{21}^2$ agree within error limits for both cases.
## 3.2 $T_S$ fluctuation dominated model

Here, we perform a comparison between covariance kernels for the spin temperature fluctuation model, similarly to what done in the previous section. The results are shown in Fig. 4. As seen in Fig. 2 and discussed in Section 2.4, the SNR is larger than 1 also for ≈10 nights of observations, suggesting better chances of detectability and smaller uncertainty ranges. Indeed, all three covariance kernels contain the signal within the $2\sigma$ uncertainty bands around their recovered signal, and the uncertainty range for the VAE kernel is ≈2 orders of magnitude smaller than the equivalent case for the $x_{H1}$ fluctuation dominated model. As discussed in the previous section, we compare the lower bounds of the recovered signal with the thermal noise uncertainty from Fig. 2 to classify the recovery as a detection or an upper limit.

For ≈100 nights of observations, the SNR is ≫10 across all k-bins, so that the recovered signal is expected to reproduce the input signal with a significantly narrower $2\sigma$ uncertainty range, provided that the covariance kernel chosen is a reliable estimate of the covariance of the input 21-cm signal. Indeed, from Fig. 4 we note that the VAE kernel fully recovers the signal with less than one order of magnitude uncertainty. However, the Exponential kernel contains the 21-cm signal only in the lowest k-bins and shows significant bias in the estimated power at higher k-bins. Its broad $2\sigma$ uncertainty shows that the recovery just provides upper limits even in the low k-bins. On the other hand, the Matern32 kernel performs better, and provides a successful detection, albeit with broader uncertainty ranges than the VAE kernel recovery. This suggests that the Exponential kernel is definitely not a good match for the covariance of the input 21-cm signal, and, as highlighted by Kern & Liu (2021), would lead to significant errors in the estimated physics, if used. The Matern32 kernel is better; however, the VAE kernel improves upon it even further. This problem with the Exponential kernel appears in the ≈100 and not in the ≈10 nights of observation due to the similarity of power and shape of the excess-noise and 21-cm signal components. Thus, a covariance kernel which is not a reliable estimate of the covariance of the input 21-cm signal would not be able to distinguish between the two, and may either ignore both equally, or identify one over the other purely by chance. It can also be argued that the only reason for any ‘detection’ at low k-bins using the Exponential kernel could possibly just be the detection of the excess noise component, wrongly interpreted as the 21-cm signal one.

The recovered values for the coherence scales and variances for the intrinsic sky, mode mixing, and excess noise components, as well as those for the hyperparameters $x_1$ and $x_2$ and associated variance for the 21-cm signal are listed in the 6th and 7th columns of Table 1 along with their 68 per cent confidence intervals. As expected, we find an improvement in recovery of the input values for ≈100 nights of observation in comparison to ≈10 nights, particularly for the coherence-scale hyperparameter.

While we note better estimates for $l_{sky}$ and $\sigma_{sky}^2/\sigma_{noise}^2$ for ≈100 nights of observations in both 21-cm signal models, in the $x_{H1}$ fluctuation dominated model the input $\sigma_{sky}^2/\sigma_{noise}^2$ is not included within the estimate error of the recovered values. A similar behaviour is observed in the recovery of the variance for $f_{mix}$ and $f_{ex}$ in the $x_{H1}$ fluctuation dominated model for ≈100 nights of observations. Lastly, we also note that the $x_1$ and $x_2$ hyperparameters and associated variance for the 21-cm signal in both models agree within the error estimates. While the estimated variance for ≈100 nights of observations is higher, it also has a broad 68 per cent confidence interval.

<figure><img src="image_7.png" /><figcaption>Figure 5. As Fig. 3 for the CRASH simulation of reionization.</figcaption></figure>

## 3.3 CRASH simulations

The power spectra resulting from the recovery using the three kernels in the case of the CRASH simulation are shown in Fig. 5. As now the input 21-cm signal has a power which lies in between the two GRIZZLY models, this translates into an intermediate SNR across k-bins. Due to this, we are able to successfully contain the input signal in the $2\sigma$ ranges of the recovered signals in most k-bins for ≈10 nights of observation with all kernels, but given the thermal noise uncertainty power, it is still classified as an upper limit. However, for ≈100 nights, we note that while the VAE kernel does an excellent job of recovering the signal with narrow $2\sigma$ uncertainty bands, they still indicate that the recovery is an upper limit. On the other hand, the Matern-class functions-based kernels underestimate the signal and do not contain the input 21-cm signal within $2\sigma$ uncertainty bands for some k-bins, despite them being significantly larger. When comparing with the thermal noise uncertainty, we note that they only provide upper limits for the signal in most k-bins. This result is similar to that with the $T_S$ fluctuation dominated model with GRIZZLY, and thus is due to the same reasons discussed in Section 3.2.

## 3.4 Redshift dependence

To evaluate the performance of the VAE kernel at different redshifts, we use the simulations at $z = 8.30, 9.16$, and $10.11$ of both GRIZZLY models, and compare the performance of the VAE kernel for noise levels equivalent to ≈100 nights of observations. It should be noted that we use the VAE kernels trained for the respective redshifts to avoid making the assumption of the kernels being redshift agnostic.

To analyse our recovery technique, we check how *accurately* it recovers the input signal, and how *precise* the results it reports are. For this purpose, we use the quantities defined below, with their values listed in Table 2:
<table><caption>**Table 2.** Average bias, $\langle \text{PS}_{\text{rec}}/\text{PS}_{\text{in}} \rangle_k$, average scaled uncertainty, $\langle \text{Err}_{\text{rec}}/\text{PS}_{\text{in}} \rangle_k$, and average z-score, $\langle \text{z-score} \rangle_k$, at various redshifts. These values are obtained when the signal recovery is done employing the VAE kernel with the $x_{H1}$ (left) and the $T_S$ (right) fluctuation dominated model and $\approx 100$ nights of observations.</caption><thead><tr><th rowspan="2">z</th><th colspan="5">x<sub>H1</sub> fluctuation dominated</th><th colspan="5">T<sub>S</sub> fluctuation dominated</th></tr><tr><th>⟨PS<sub>rec</sub>/PS<sub>in</sub>⟩<sub>k</sub></th><th>⟨Err<sub>rec</sub>/PS<sub>in</sub>⟩<sub>k</sub></th><th>z-score<sub>min</sub></th><th>⟨z-score⟩<sub>k</sub></th><th>z-score<sub>max</sub></th><th>⟨PS<sub>rec</sub>/PS<sub>in</sub>⟩<sub>k</sub></th><th>⟨Err<sub>rec</sub>/PS<sub>in</sub>⟩<sub>k</sub></th><th>z-score<sub>min</sub></th><th>⟨z-score⟩<sub>k</sub></th><th>z-score<sub>max</sub></th></tr></thead><tbody><tr><td>8.30</td><td>3.90</td><td>14.88</td><td>+ 0.46</td><td>+ 0.74</td><td>+ 0.85</td><td>0.69</td><td>0.50</td><td>-4.00</td><td>-1.36</td><td>-1.25</td></tr><tr><td>9.16</td><td>8.57</td><td>21.12</td><td>+ 0.73</td><td>+ 0.82</td><td>+ 0.88</td><td>0.90</td><td>0.58</td><td>-0.79</td><td>-0.60</td><td>+ 0.23</td></tr><tr><td>10.11</td><td>29.29</td><td>78.22</td><td>+ 0.87</td><td>+ 1.02</td><td>+ 1.16</td><td>0.34</td><td>0.83</td><td>-3.50</td><td>-1.96</td><td>-0.92</td></tr></tbody></table>

<figure><img src="image_8.png" /><figcaption>**Figure 6.** Top panel: Input 21-cm signal for the $x_{H1}$ (dashed lines) and the $T_S$ (solid) fluctuation dominated model are shown for $z = 8.30$ (blue), 9.16 (orange), and 10.11 (green). Middle panel: Recovered 21-cm signal and its associated uncertainty divided by the input 21-cm signal to give the bias and scaled uncertainty (the average values over k-bins for these quantities are listed in Table 2) for the $x_{H1}$ fluctuation dominated model. For comparison, the line of zero bias (i.e. recovered signal perfectly matching the input signal) is also shown. Bottom panel: Same as the middle panel, but for the $T_S$ fluctuation dominated model.</figcaption></figure>

(i) *Average bias*, $\langle \text{PS}_{\text{rec}}/\text{PS}_{\text{in}} \rangle_k$: this is the average bias given as the ratio between the recovered ($\text{PS}_{\text{rec}}$) and input 21-cm signal ($\text{PS}_{\text{in}}$) power spectra, averaged across all k-bins. An accurate recovery has a value close to 1, with higher (lower) values indicating an over(under)estimation of the input 21-cm signal. We note that the deviation of average bias from 1 increases significantly with redshift for the $x_{H1}$ fluctuation dominated model. However, the trend is not so clear for the $T_S$ fluctuation dominated signal, as the power and slope of the excess noise component match closely those of the input 21-cm signal at $z = 8.30$ and $k \le 0.2$ h Mpc$^{-1}$, which makes differentiating between them more difficult. However, we still get an average bias of $\approx 0.70$, suggesting only a minor underestimation, which is primarily observed at the lowest k-bins (see Fig. 6).

(ii) *Average scaled uncertainty*, $\langle \text{Err}_{\text{rec}}/\text{PS}_{\text{in}} \rangle_k$: the $2\sigma$ uncertainty (given by the shaded region in Figs 3, 4, 5, and referred to as $\text{Err}_{\text{rec}}$ hereafter) on the recovered power spectrum ($\text{PS}_{\text{rec}}$) gives

the precision of recovery. However, except for cases of extremely poor recovery, the absolute value of $\text{Err}_{\text{rec}}$ generally depends on the absolute value of the corresponding $\text{PS}_{\text{rec}}$. Thus, to compare different recoveries for a given $\text{PS}_{\text{in}}$, we need to convert it into a unitless quantity. We first tried to do this by calculating the ratio $\text{Err}_{\text{rec}}/\text{PS}_{\text{rec}}$ as a measure of the precision of the recovery. However, this ratio contains no information on the accuracy of $\text{PS}_{\text{rec}}$. Indeed, one could have a precise recovery, i.e. a low $\text{Err}_{\text{rec}}/\text{PS}_{\text{rec}}$, but an inaccurate $\text{PS}_{\text{rec}}$, i.e. a bias which deviates significantly from 1 (see point above). This ratio, then, is not useful, as it does not quantify the overall quality of recovery. To overcome this issue, we use the ratio $\text{Err}_{\text{rec}}/\text{PS}_{\text{in}}$ instead. As the magnitude of $\text{Err}_{\text{rec}}$ depends on that of $\text{PS}_{\text{rec}}$, it also carries the information of the bias in recovery. Further, as we divide by $\text{PS}_{\text{in}}$, the scaled uncertainty becomes independent of the specific $\text{PS}_{\text{in}}$ being recovered. This allows us to use it to compare between cases with different $\text{PS}_{\text{in}}$, such as at different redshifts (see top panel of Fig. 6 where the power spectrum varies due to the time evolution of ionizing bubbles) and different physical models (as discussed in Sections 2.1). We use this generalized comparison, and call it the *scaled* uncertainty for each k-bin. Averaging this across all k-bins provides a handy quantity to compare the precision of recovery for input 21-cm signals with different physical properties. For example, in Fig. 6 we observe that in the $x_{H1}$ fluctuation dominated model, while small-scale variability due to partial reionization is restricted, it still has variability tied to the large-scale distribution of neutral Hydrogen, which increases for lower redshifts. This boosts the power at large scales, corresponding to the low k-bins. The same can also be seen in the $T_S$ fluctuation dominated model, although its overall power is boosted, as it allows small-scale variability in $\delta T_b$ as well. Using the scaled uncertainty, we can compare the precision of recovery across redshift for both cases. The difference in $\langle \text{Err}_{\text{rec}}/\text{PS}_{\text{in}} \rangle_k$ is quite significant for the $x_{H1}$ fluctuation dominated model, going from $\approx 15$ at $z = 8.30$ to $\approx 80$ at $z = 10.11$.

(iii) *z-score*, $\langle \text{z-score} \rangle_k$: The z-score (Kirch 2008) is a popular quantity to evaluate quality of recovery. In our case, it is defined as $\frac{\text{PS}_{\text{rec}}-\text{PS}_{\text{in}}}{\sigma}$ or $\frac{\text{PS}_{\text{rec}}-\text{PS}_{\text{in}}}{\text{Err}_{\text{rec}}/2}$ at each k-bin, and it measures how much the recovered signal deviates from the input 21-cm signal, in units of standard deviation of the recovered signal. The z-score is thus a more explicit method to combine into a single quantity the information provided by the bias along with that of the uncertainty. The only possible issue is that $\sigma(= \text{Err}_{\text{rec}}/2)$ depends on $\text{PS}_{\text{rec}}$ and thus their ratio would mask the accuracy of recovery as discussed in the point above (see below for an example case). Further, we note that we cannot just report an average of z-scores across all k-bins, as the distribution of z-scores is not necessarily Gaussian. Indeed, we find that while it is approximately Gaussian for the $x_{H1}$ fluctuation dominated model, this is not the case for the $T_S$ fluctuation dominated model. Thus, we report the minimum and maximum z-scores (z-score$_{min}$ and z-score$_{max}$) along with the average
($\langle z\text{-score} \rangle_k$). When $\langle z\text{-score} \rangle_k > (<)0$, its exact value quantifies the extent of over(under)prediction. We note that in the $x_{H1}$ fluctuation dominated model, the average z-score worsens with increasing redshift, consistently with the behaviour of the average bias and average scaled uncertainty. This trend is not detected for the $T_S$ fluctuation dominated signal, due to the same reasons discussed above for the average bias in (i). We also note that at $z = 10.11$, $\langle z\text{-score} \rangle_k \approx -2$ and $+1$ in the $T_S$ and $x_{H1}$ fluctuation dominated model, respectively, naively suggesting a better recovery for the latter. This, though, is not correct, but simply a consequence of the very broad error bars and the inverse proportionality of $\langle z\text{-score} \rangle_k$ with the error. This reasoning exposes the limitation of the z-score. Indeed, by looking at $\langle \text{PS}_{\text{rec}}/\text{PS}_{\text{in}} \rangle_k$ and $\langle \text{Err}_{\text{rec}}/\text{PS}_{\text{in}} \rangle_k$ for recovery of the $x_{H1}$ fluctuation dominated signal (see Table 2), we note that the deviation from zero bias (obtained when $\langle \text{PS}_{\text{rec}}/\text{PS}_{\text{in}} \rangle_k = 1$) and zero uncertainty (when $\langle \text{Err}_{\text{rec}}/\text{PS}_{\text{in}} \rangle_k = 0$) is significantly higher than for the $T_S$ fluctuation dominated model. In fact, these numbers suggest a better quality of recovery in the latter case, which is understandable as the SNR in this model is higher. Thus, while we report the z-score due to its popularity, we recommend using the average bias and scaled uncertainty for evaluating the quality of recovery.

The trends in various quantities discussed above are linked to the physical nature of the 21-cm power spectrum and its redshift evolution. The drop in the SNR with increasing redshift [due to a decrease in signal power as shown in Fig. 6 and explained in (ii)], leads to a worsening of the average bias and scaled uncertainty, especially for the $x_{H1}$ fluctuation dominated model. As already mentioned, when the excess noise and the input 21-cm signal have similar power and slope (as at $z = 8.30$ for the $T_S$ fluctuation dominated case), we observe limitations in the capability of differentiating among the two, but the effects are minor and the trends of average bias, scaled uncertainty, and z-score for the $x_{H1}$ fluctuation dominated model are also observed in the $T_S$ fluctuation dominated model when going from $z = 9.16$ to $z = 10.11$.

Thus, we find that the VAE kernel does not add significant biases, with its recovery and associated uncertainty largely scaling with the physical properties of the 21-cm signal.

# 4. DISCUSSION

## 4.1 Role of the excess noise component

In M20, it was noted that the excess noise component was poorly constrained, and thus the combined excess noise and 21-cm signal components were jointly recovered, as separating them was not statistically justifiable. Thus, it is important to understand how well constrained the excess noise component has to be, in order to separate the 21-cm signal from it.

To explore this, we looked at the $x_{H1}$ fluctuation dominated model of the 21-cm signal at $z = 9.16$. We note that in the $\approx 100$ nights case for the chosen excess noise component, the power spectrum recovered with the VAE kernel is slightly overestimated, with the $2\sigma$ bands on both sides of it spanning $\approx 2$ orders of magnitude. We then generated a range of input excess noise components (by varying either the coherence-scale hyperparameter or the variance), and analysed the effects on the recovery of the 21-cm signal power spectrum using the VAE kernel by looking at the average bias, scaled uncertainty and z-score, as defined in Section 3.4. The results are reported in Table 3, where $f_{\text{var}}$ is the factor by which the coherence-scale hyperparameter and the variance from the results of M20 were scaled. We see that, overall, the average bias and scaled uncertainty are reduced when increasing the coherence-scale hyperparameter or

<table><caption>**Table 3.** Average bias, $\langle \text{PS}_{\text{rec}}/\text{PS}_{\text{in}} \rangle_k$, average scaled uncertainty, $\langle \text{Err}_{\text{rec}}/\text{PS}_{\text{in}} \rangle_k$, and average z-score, $\langle z\text{-score} \rangle_k$, for the $x_{H1}$ fluctuation dominated model with $\approx 100$ nights of observations at $z = 9.16$. These values are obtained when the signal recovery is done employing the VAE kernel and multiplying the coherence-scale hyperparameter, $l_{\text{ex}}$, or the variance, $\sigma_{\text{ex}}^2$, listed in the second column of Table 1 (and as used in M20) by a factor $f_{\text{var}}$.</caption><thead><tr><th>f<sub>var</sub></th><th>&lang;PS<sub>rec</sub>/PS<sub>in</sub>&rang;<sub>k</sub></th><th>&lang;Err<sub>rec</sub>/PS<sub>in</sub>&rang;<sub>k</sub></th><th>z-score<sub>min</sub></th><th>&lang;z-score&rang;<sub>k</sub></th><th>z-score<sub>max</sub></th></tr></thead><tbody><tr><td></td><td></td><td></td><td><i>l</i><sub>ex</sub></td><td></td><td></td></tr><tr><td>0.5</td><td>12.88</td><td>31.37</td><td>+0.88</td><td>+1.04</td><td>+1.19</td></tr><tr><td>0.75</td><td>9.61</td><td>26.75</td><td>+0.77</td><td>+0.89</td><td>+0.97</td></tr><tr><td>1.0</td><td>8.57</td><td>21.12</td><td>+0.73</td><td>+0.82</td><td>+0.88</td></tr><tr><td>1.25</td><td>8.04</td><td>27.40</td><td>+0.75</td><td>+0.82</td><td>+0.89</td></tr><tr><td>1.5</td><td>7.77</td><td>15.65</td><td>+0.78</td><td>+0.84</td><td>+0.90</td></tr><tr><td>1.75</td><td>5.64</td><td>14.76</td><td>+0.81</td><td>+0.87</td><td>+0.93</td></tr><tr><td>2.0</td><td>4.55</td><td>19.38</td><td>+0.81</td><td>+0.85</td><td>+0.94</td></tr><tr><td></td><td></td><td></td><td><i>&sigma;</i><sub>ex</sub><sup>2</sup></td><td></td><td></td></tr><tr><td>1.00</td><td>8.57</td><td>21.12</td><td>+0.73</td><td>+0.82</td><td>+0.88</td></tr><tr><td>0.75</td><td>6.64</td><td>32.81</td><td>+0.75</td><td>+0.80</td><td>+0.88</td></tr><tr><td>0.50</td><td>4.77</td><td>16.89</td><td>+0.66</td><td>+0.71</td><td>+0.78</td></tr><tr><td>0.25</td><td>3.24</td><td>9.13</td><td>+0.58</td><td>+0.62</td><td>+0.67</td></tr></tbody></table>

decreasing the variance, while no substantial difference is observed in the average z-score. This is possibly because the bias and scaled uncertainty decrease at the same rate, and thus their effects roughly cancel out.

## 4.2 Overall performance of the VAE kernel

Usually it is expected that the recovery of the 21-cm signal from an overall data set with foreground components, noise and signal, is not possible if SNR<1. In Section 3.1, we have indeed shown that Matern-class functions-based kernels are unable to contain the input 21-cm signal within their $2\sigma$ uncertainty bands when SNR < 1. However, the VAE kernel is not only able to do so, but also to recover the overall shape of the power spectrum, as seen in the top panel of Fig. 3.

Further, as highlighted by Kern & Liu (2021), misestimation of the covariance kernel can significantly hamper signal detection given the currently used normalization and bias correction schemes. This means that as the Matern-class functions model the 21-cm signal only approximately, their results can be significantly biased for more complex models of the 21-cm signal, as given example by the $T_S$ fluctuation dominated model of GRIZZLY and the 21-cm signal model from the CRASH simulations. Indeed, in these cases the Matern-class functions-based kernels fail to recover the signal also for noise levels equivalent to $\approx 100$ nights of observation, even when $\langle \text{SNR} \rangle_k \approx 10$. The VAE kernel does not suffer from such a limitation and performs well also when used with an input signal from the CRASH simulations. This shows that the VAE covariance kernel is a more robust estimate of the covariance of the 21-cm signal, and can successfully report a detection within $2\sigma$ uncertainty regardless of the exact physical properties of the observed 21-cm differential brightness temperature power spectrum. Lastly, we note that it performs well across all redshifts analysed here. This reconfirms the robustness of the VAE kernel in constraining the 21-cm signal, with an increase/decrease in uncertainty tied to the $\langle \text{SNR} \rangle_k$ of the signal itself.

Overall, the $2\sigma$ uncertainty bands given by the VAE kernel contain the signal in all cases discussed here. We consider the recovery limit of the VAE kernel in terms of SNR averaged over k-bins to be that from $\approx 10$ nights of observation in the $x_{H1}$ fluctuation dominated model, i.e. $\langle \text{SNR} \rangle_k = 5 \times 10^{-2}$. For values lower than this, we do
not expect the VAE kernel-based recovery to contain the input signal within its $2\sigma$ uncertainty bands across $k = [0.05, 1.00] h \text{ Mpc}^{-1}$. We explore one such case in Appendix A, and indeed find that the recovery does not contain the input signal within its $2\sigma$ uncertainty bands, but instead provides upper limits. Yet it still outperforms the Matern-class functions-based kernels in recovering the shape of the power spectrum, spread of uncertainty on recovery, as well as the reported upper limits.

## 4.3 Limitations

In this work, we reduce the possibility of biases in the EoR covariance kernel by incorporating a more physically motivated covariance. We showcase the robustness of the generated VAE covariance kernel by testing it against not just mock 21-cm signals obtained with GRIZZLY, but with signals generated using very different frameworks as shown in Section 3.3 and in Appendix A. This has also been demonstrated in Mertens, Bobin & Carucci (2023).

However, biases are still possible, especially if the true signal, and thus its covariance, is fundamentally different from what we obtain with our simulation codes. One way to further minimize this bias is to use mock data obtained from different codes to train the covariance kernel, which we plan to do in further upgrades of our pipeline. In the future, we will also investigate methods to reduce the dependence on the prior by using different normalization and bias correction schemes as suggested by Kern & Liu (2021).

## 5. SUMMARY

The LOFAR EoR KSP team strives for a successful detection of the 21-cm signal from the EoR at $z \approx 7-11$. Past theoretical models indicated that at least 1000 h of observation would be necessary to lead to a successful detection (Mertens, Ghosh & Koopmans 2018), while M20 provided upper limits using 141 h ($\approx 10$ nights) of observation. In this respect, an optimal choice of the covariance kernel for the 21-cm signal component is crucial. Indeed, as shown in Kern & Liu (2021), given the currently used normalization and bias correction scheme, a mismatch between the adopted and the actual covariance kernel of the 21-cm signal can induce a significant signal loss, which can in turn lead to incorrect astrophysical interpretations from any 'successful' detection.

To improve the choice of the 21-cm signal covariance kernel, Mertens, Bobin & Carucci (2023) introduce a ML method which employs a VAE-based algorithm. As the training done using VAE is not limited by the form of the specific function (as in the case of Matern-class functions), nor by the kernels of the training data sets (as in the case of a simple AE), it allows to reproduce the covariance kernel of the 21-cm signal with a greater flexibility. This is showcased in terms of the robustness of the VAE-based kernel's performance in comparison with previously used kernels based on Matern-class functions.

We show that the result on using the VAE kernel is able to contain the input 21-cm signal within its $2\sigma$ uncertainty band in all cases explored where $\langle \text{SNR} \rangle_k \gtrsim 5 \times 10^{-2}$. It is also usually better than the results from Matern-class functions-based covariance kernels in recovering the overall shape of the power spectrum of the signal. A key result in this paper is that Matern-class functions-based kernels are unable to recover the 21-cm signal for the $T_S$ fluctuation dominated model even for $\approx 100$ nights of observation, for which $\langle \text{SNR} \rangle_k \approx 10$, while a recovery with the VAE kernel is successful. A similar result is obtained also with a

21-cm signal generated using the CRASH simulations, thus clearly indicating that the Matern-class functions-based kernels do not correctly match the covariance of more complex models of the signal. Thus, this analysis suggests that using the VAE kernel can mitigate to a significant extent the issues highlighted by Kern & Liu (2021) given no change to the normalization and bias correction schemes.

Further, we show that the behaviour of the VAE kernel is consistent across all redshifts of interest, with changes in its performance strongly correlating with the neutral hydrogen distribution, which changes the strength of the resultant power spectrum, and thus the $\langle \text{SNR} \rangle_k$. This suggests that the VAE-based kernel can be used for any choice of redshift without additional correction factors, making the algorithm developed here directly applicable to LOFAR data at $z \approx 8$ to 10, whose results can then be compared with results from telescopes like HERA.

We also explore the effects that the properties of the excess noise component identified in M20 have on the recovery of the 21-cm signal. As expected, we find that having a higher coherence scale or a lower variance for the components leads to better recovery.

In companion papers we will apply the VAE kernel to the $\approx 10$ nights of LOFAR data used in M20, and explore the range of theoretical models which are consistent with the upper limits provided by the VAE kernel, as done in Ghara et al. (2020). Applying the VAE kernel to observations much longer than $\approx 10$ nights requires a significant improvement in the modelling of the intrinsic sky component, which would eventually be limited by the confusion noise due to the angular resolution of LOFAR. Further improvements, such as noise mitigation, can be implemented by choosing data from nights with better ionospheric conditions and lesser contribution from RFI flagging. All these aspects are currently being explored by the LOFAR EoR KSP team.

## ACKNOWLEDGEMENTS

AA thanks the EoR research group at MPA, as well as Hitesh Kishore Das and Abinaya Swaruba Rajamuthukumar for helpful discussions. We would also like to thank the referee for their helpful comments that helped improve this paper significantly.

## DATA AVAILABILITY

The VAE kernels generated for $z = 8.30, 9.16$, and $10.11$ shall be made public on publication. The PS_EOR package used to run GPR and ML-GPR can be found at https://gitlab.com/flomertens/ps_eor. The underlying training, testing, and validation data sets generated using GRIZZLY, and the mock 21-cm signal data generated using GRIZZLY and CRASH can be shared upon reasonable request.

## REFERENCES

* Aigrain S., Foreman-Mackey D., 2023, ARA&A, 61, 329

* Becker R. H. et al., 2001, AJ, 122, 2850

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

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

* Ciardi B., Ferrara A., 2005, Space Sci. Rev., 116, 625

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

* Ciardi B., Ferrara A., Marri S., Raimondo G., 2001, MNRAS, 324, 381

* Clément B. et al., 2012, A&A, 538, A66

* Davies F. B. et al., 2018, ApJ, 864, 142
Eide M. B., Graziani L., Ciardi B., Feng Y., Kakiichi K., Di Matteo T., 2018, MNRAS, 476, 1174

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

Fan X. et al., 2006, AJ, 131, 1203

Foreman-Mackey D., Hogg D. W., Lang D., Goodman J., 2013, PASP, 125, 306

Furlanetto S. R., 2016, The 21-cm Line as a Probe of Reionization. Springer International Publishing, Cham, p. 247

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

Gehlot B. K. et al., 2019, MNRAS, 488, 4271

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

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

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

Glatzle M., Ciardi B., Graziani L., 2019, MNRAS, 482, 321

Goodfellow I., Bengio Y., Courville A., 2016, Deep Learning. MIT Press, USA

Graziani L., Maselli A., Ciardi B., 2013, MNRAS, 431, 722

Graziani L., Ciardi B., Glatzle M., 2018, MNRAS, 479, 4320

Greig B., Mesinger A., 2015, MNRAS, 449, 4246

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

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

Hothi I. et al., 2021, MNRAS, 500, 2264

Kern N. S., Liu A., 2021, MNRAS, 501, 1463

Kerrigan J. R. et al., 2018, ApJ, 864, 131

Khandai N., Di Matteo T., Croft R., Wilkins S., Feng Y., Tucker E., DeGraff C., Liu M.-S., 2015, MNRAS, 450, 1349

Kingma D. P., Welling M., 2013, arXiv:1312.6114

Kingma Diederik P., Welling Max 2019, Foundations and Trends® in Machine Learning, 12, 307–392.

Kirch W. ed. 2008, z-Score. Springer, Netherlands, Dordrecht, p. 1484

Kullback S., Leibler R. A., 1951, Ann. Math. Stat., 22, 79

Liu A., Shaw J. R., 2020, PASP, 132, 062001

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

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

Maselli A., Ciardi B., Kanekar A., 2009, MNRAS, 393, 171

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

Mertens F. G. et al., 2020, MNRAS, 493, 1662 (M20)

Mertens Florent G., Bobin Jérôme Carucci Isabella P 2024, MNRAS, 527, 3517–3531.

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

Morales M. F., Wyithe J. S. B., 2010, ARA&A, 48, 127

Morales M. F., Beardsley A., Pober J., Barry N., Hazelton B., Jacobs D., Sullivan I., 2019, MNRAS, 483, 2207

Mortlock D., 2016, Quasars as Probes of Cosmological Reionization. Springer International Publishing, Cham, p. 187

Pinheiro Cinelli L., Araújo Marins M., Barros da Silva E. A., Lima Netto S., 2021, Variational Autoencoder. Springer International Publishing, Cham, p. 111

Planck Collaboration XLVII 2016, A&A, 596, A108

Pober J. C. et al., 2014, ApJ, 782, 66

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

Rasmussen C. E., Williams C. K. I., 2006, Gaussian Processes for Machine Learning (Adaptive Computation and Machine Learning). MIT Press, USA

Schenker M. A. et al., 2013, ApJ, 768, 196

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

Stein M. L., 1999, Interpolation of Spatial Data: Some Theory for Kriging. Springer Science and Business Media, Berlin

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

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

Zaroubi S., 2013, The Epoch of Reionization. Springer, Berlin, Heidelberg, p. 45

## APPENDIX A: INVESTIGATING A LOW-SNR CASE

We use the seminumerical code 21CMFAST (Mesinger & Furlanetto 2007; Greig & Mesinger 2015) to generate a mock 21-cm signal at $z = 9.1$ within a box of length 400 cMpc. We adopt the values $\zeta = 30$ per cent, $R_{mfp} = 15$ cMpc, and $T_{vir}^{feed} = 5 \times 10^4$ K (see Greig & Mesinger 2015 for more details).

The power spectra recovered using the three kernels are shown in Fig. A1. As in this model the input 21-cm signal is much weaker than in those discussed in the main text, for $\approx 10$ nights of observation ($\langle (\text{SNR})_k \rangle \approx 1.6 \times 10^{-2}$) the input signal is never contained within the $2\sigma$ uncertainty bounds of the recovered signal with the exception of the highest $k$-bins, where the SNR is higher. When comparing with the thermal noise uncertainty, we find that all three kernels manage to provide upper limits, with the VAE kernel also recovering the overall shape. A detection is not possible because of the low SNR, with the uncertainty on the thermal noise being almost two orders of magnitude larger than the input 21-cm signal. For $\approx 100$ nights ($\langle (\text{SNR})_k \rangle \approx 0.3$), the signal is within the $2\sigma$ uncertainty bounds for all kernels, while still providing upper limits and not a detection because the SNR is still low despite the increased integration time. The VAE kernel has the tightest constraints on the uncertainty. While the VAE kernel overestimates the input signal, it does an excellent job in recovering its overall shape.

<figure><img src="image_9.png" /><figcaption>**Figure A1.** As Fig. 3 for 21CMFAST simulations of reionization.</figcaption></figure>

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

© The Author(s) 2023.

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