# Revised LOFAR upper limits on the 21-cm signal power spectrum at $z \approx 9.1$ using machine learning and gaussian process regression

Anshuman Acharya<sup>1,★</sup> Florent Mertens<sup>2</sup> Benedetta Ciardi,<sup>1</sup> Raghunath Ghara<sup>3,4</sup> Léon V. E. Koopmans<sup>5</sup> and Saleem Zaroubi<sup>5,6</sup>

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

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

<sup>3</sup>Haverford College, 370 Lancaster Ave, Haverford, PA 19041, USA

<sup>4</sup>Center for Particle Cosmology, Department of Physics and Astronomy, University of Pennsylvania, Philadelphia, PA 19104, USA

<sup>5</sup>Kapteyn Astronomical Institute, University of Groningen, PO Box 800, NL-9700AV Groningen, The Netherlands

<sup>6</sup>Astrophysics Research Centre of the Open University of Israel, Ra'anana 4353701, Israel

Accepted 2024 August 16. Received 2024 August 7; in original form 2024 June 3

## ABSTRACT

The use of Gaussian Process Regression (GPR) for foregrounds mitigation in data collected by the LOw-Frequency ARray (LOFAR) to measure the high-redshift 21-cm signal power spectrum has been shown to have issues of signal loss when the 21-cm signal covariance is misestimated. To address this problem, we have recently introduced covariance kernels obtained by using a Machine Learning based Variational Auto-Encoder (VAE) algorithm in combination with simulations of the 21-cm signal. In this work, we apply this framework to 141 h (≈10 nights) of LOFAR data at $z \approx 9.1$, and report revised upper limits of the 21-cm signal power spectrum. Overall, we agree with past results reporting a 2-$\sigma$ upper limit of $\Delta_{21}^2 < (80)^2 \text{ mK}^2$ at $k = 0.075 \text{ h Mpc}^{-1}$. Further, the VAE-based kernel has a smaller correlation with the systematic excess noise, and the overall GPR-based approach is shown to be a good model for the data. Assuming an accurate bias correction for the excess noise, we report a 2-$\sigma$ upper limit of $\Delta_{21}^2 < (25)^2 \text{ mK}^2$ at $k = 0.075 \text{ h Mpc}^{-1}$. However, we still caution to take the more conservative approach to jointly report the upper limits of the excess noise and the 21-cm signal components.

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

## 1 INTRODUCTION

The eras between redshifts $z \sim 30$ to $z \sim 5$, when the first astro-physical objects formed and started to emit photons that ionized the neutral hydrogen in the intergalactic medium (IGM), are known as the Cosmic Dawn and the Epoch of Reionization (EoR). The brightness fluctuation of the 21-cm line of neutral hydrogen as observed in emission or absorption against the Cosmic Microwave Background, is one of the most promising probes to study this redshift range (e.g. Field 1959; Madau, Meiksin & Rees 1997; Shaver et al. 1999; Tozzi et al. 2000; Ciardi & Madau 2003; Zaroubi 2013). Even a statistical upper limit estimate of the strength of such fluctuations would allow us to constrain models of structure formation in the early Universe (Ghara et al. 2020; Mondal et al. 2020; Greig et al. 2021a, b; Abdurashidova et al. 2022).

While multiple interferometric low-frequency radio telescopes have been designed to search for this signal (e.g. PAPER,<sup>1</sup> MWA,<sup>2</sup>

HERA,<sup>3</sup> NenuFAR<sup>4</sup> and the upcoming SKA<sup>5</sup>), in this work we focus on LOw-Frequency ARray (LOFAR).<sup>6</sup> The most recent and stringent 2-$\sigma$ upper limits on the power of the 21-cm signal from 141 h (≈10 nights) of observation with LOFAR at $z=9.1$, were provided by Mertens et al. (2020, hereafter M20) as $\Delta_{21}^2(k = 0.075 \text{ h Mpc}^{-1}) < (73)^2 \text{ mK}^2$.

As the 21-cm signal is seen through foregrounds that are several orders of magnitude stronger than the signal itself, one of the challenges faced by all 21-cm signal experiments is to develop an accurate modelling and removal of the foregrounds. The LOFAR EoR Key Science Project (KSP) team employs Gaussian Process Regression (GPR) to model different components of the observed data, as described by Mertens, Ghosh & Koopmans (2018), Gehlot et al. (2019), and Hothi et al. (2021). In this method, each contributing component of the data is represented by a Gaussian Process, which completely depends on the mean and the covariance between the

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

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

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

<footnote n="3"><sup>3</sup>Hydrogen Epoch of Reionization Array, https://reionization.org/</footnote>

<footnote n="4"><sup>4</sup>New Extension in Nançay Upgrading LOFAR, https://nenufarobs-nancay.fr/en</footnote>

<footnote n="5"><sup>5</sup>Square Kilometre Array, https://www.skao.int/en</footnote>

<footnote n="6"><sup>6</sup>Low-Frequency Array, http://www.lofar.org</footnote>
range of frequencies considered. This covariance can thus be represented as a matrix, which is also referred to as a 'covariance kernel'. However, given the choice of normalization and bias correction in the power-spectra estimation used by M20, Kern & Liu (2021) noted the dependence of the results on the choice of the covariance kernel used for the 21-cm signal so that any misestimation of the kernel could lead to spurious signal gain/loss. This in turn would affect the astrophysical interpretation of the estimated 21-cm signal power spectrum.

To improve the set of 21-cm signal covariance kernels, Mertens, Bobin & Carucci (2024) introduced the concept of a learned covariance kernel, using a Variational Auto-Encoder (VAE) algorithm applied to 21-cm signal simulations. This was employed by Acharya et al. (2024, hereafter A24) after training on GRIZZLY simulations of reionization (Ghara, Choudhury & Datta 2015; Ghara et al. 2018, 2020), and tested against various mock data sets. In this work, we apply the VAE-based covariance kernels trained at $z = 9.16$ by A24 to the same $\approx 10$ nights of data used by M20, and compare its performance with that of the Exponential Matern-class function based kernel employed by M20.

In Section 2, we briefly highlight the covariance kernels used for the different components of the data set; in Section 3, we present two cases of recovery, and compare the obtained upper limits to those presented in M20. In Section 4, we discuss the quality of our results by comparing the residuals on applying our model to the data, and also the limitations of the method. Finally, in Section 5, we give our conclusions.

## 2 METHODOLOGY

GPR (Rasmussen & Williams 2006; Aigrain & Foreman-Mackey 2023) has been used to model radio data in frequency space as noisy observations of the form $y = f(x) + \epsilon$, with $\epsilon$ Gaussian noise vector having variance $\sigma_{\text{noise}}^2$ (see e.g. A24, M20 and Munshi et al. 2024). The vector $f(x)$ can be split into the 21-cm signal ($f_{21}$) and foregrounds. To the latter contribute an intrinsic sky component ($f_{\text{sky}}$) and mode-mixing contaminants ($f_{\text{mix}}$). Further, M20 identified an additional term for bias correction, defining it as an 'excess noise' component ($f_{\text{ex}}$), which corrects for additional systematic noise after the subtraction of foregrounds. While the exact source of this excess is yet to be identified, various possibilities have been discussed in section 6.2 of M20, and the LOFAR EoR KSP team is currently working on testing each one of them. The overall datacube $y$ can be expressed as:

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

For a set of points $x$ (independent parameters in frequency), for each value $x$, $f(x)$ given by a Gaussian Process $GP(m, \kappa)$ is fully defined by its mean $m$ and covariance matrix $\kappa$. The vector $f(x)$ is then fully defined by its mean vector $m$ and covariance matrix vector $K$. The joint distribution for all random variables that share the desired covariance properties is a normal distribution vector ($N$) which can be represented as:

$$f(x) \sim N(m(x), K(x, x)). \quad (2)$$

Here, the matrix $K$ gives the covariance between the function values at any two points, and can be written as $K_{ij} = \kappa(x_i, x_j, \phi) + \delta_{ij}\sigma_{\text{noise}}^2$, where $\kappa(x_i, x_j, \phi)$ can be optimized by the choice of hyperparameters represented by $\phi$, $\delta_{ij}$ is the Kronecker-delta function and $\sigma_{\text{noise}}^2$ is the noise at $x_i$.

If we assume that the various components are uncorrelated, we can utilize the additive property of matrices. Thus, the covariance

kernels of the right-hand side of equation (1) can be represented as a single covariance kernel $K$ given as:

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

$K_{\text{sky}}$, $K_{\text{mix}}$, and $K_{\text{ex}}$ are modelled using the best-fitting Matern-class functions (Stein 1999) as done by M20 and A24:

$$k_{\text{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 (4)$$

with $\eta_{\text{sky}} = +\infty$, $\eta_{\text{mix}} = 3/2$, and $\eta_{\text{ex}} = 5/2$. Additionally, $l$ is the coherence-scale hyperparameter, with its associated variance given by $\sigma^2$, $r$ is the absolute difference between the frequencies of two sub-bands, $\kappa_{\eta}$ is the modified Bessel function of the second kind, and $\Gamma$ is the Gamma-function. $\eta_{\text{sky}}$, $\eta_{\text{mix}}$, and $\eta_{\text{ex}}$ were obtained by M20 by assuming different values of the hyperparameter $\eta$, and finding the one that maximized the marginal likelihood (or the 'evidence'). This was done by calculating the analytical integral over $f$, which is the log-marginal-likelihood (LML, see section 2.3 in Mertens et al. 2018). While M20 used a gradient-descent-based optimization algorithm for maximizing the LML, A24, instead, introduced an MCMC sampling based approach (Foreman-Mackey et al. 2013) to estimate the hyperparameters by sampling their posterior distributions. This additionally provides a measure of the uncertainty on the hyperparameters.

The noise is modelled based on the Stokes-V visibility difference, according to the methodology laid out in M20. Further, for the 21-cm signal we use the VAE-based kernel at $z = 9.16$ built by A24 by training on GRIZZLY simulations. This VAE-based kernel is fully defined using two hyperparameters $x_1$ and $x_2$ with uninformed flat priors in linear space in the range $[-10, 10]$, and an associated variance. The kernel has been trained on a training set of $\approx 1500$ simulations with four independently variable parameters. While more hyperparameters could be employed, they do not provide any significant improvement, and thus we do not use them to avoid overfitting. A24 also used an additional testing set of $\approx 150$ simulations and found a recovery error of $\lesssim 1$ per cent for wave-modes of $k < 0.43h\text{Mpc}^{-1}$. In this work, the performance of the VAE-based kernel is compared against the results obtained in M20 when using a Matern-class function with $\eta_{21} = 1/2$.

## 3 RESULTS

Following A24, we recover the hyperparameters of the various signal components by applying GPR to the data using an MCMC approach (Foreman-Mackey et al. 2013). We adopt the same broad flat priors for the variances, and uniform priors for the coherence-scale parameters (see Table 1 and section 2.5 of A24). The recovered values are listed in Table 1, where we compare them to those obtained by M20. We find that the hyperparameters for the $f_{\text{sky}}$ and $f_{\text{ex}}$ components differ from the M20 estimates by about $2-4\sigma$, while those for $f_{\text{mix}}$ are in good agreement. This is expected with a more accurate 21-cm kernel, as it reduces the strength of the correlation of hyperparameters of the different components. This allows for better characterisation of the different components in the overall data. We discuss the details below after analysing each component's recovered power spectra.

To have a better understanding of the results quantified above, in Fig. 1, we show the recovered upper limits on the power spectra in two cases. In the first one (top panel), we present the excess noise component (dashed-dotted green line) and the 21-cm signal (blue dashes with downward arrows) separately, as done for the mock data sets in A24, and refer to this as the *individual upper limits case*. We
<table><caption>**Table 1.** Comparison of hyperparameters obtained in M20 and in this work. From left to right the columns refer to the hyperparameters used, the results from M20, the results obtained in this work, and the difference between the two. As the covariance kernels used for the 21-cm signal component are not the same, we just provide the values of the hyperparameters obtained in each case. Lastly, $\sigma_{\text{noise}}^2$ is $74 \times 10^3 \text{ mK}^2$ (see M20).</caption><thead><tr><th>Parameter</th><th>M20</th><th>This work</th><th>Difference</th></tr></thead><tbody><tr><td>$l_{\text{sky}}$</td><td>$47.5^{+3.1}_{-2.8}$</td><td>$38.6^{+2.6}_{-2.6}$</td><td>$2.3\sigma$</td></tr><tr><td>$\sigma_{\text{sky}}^2/\sigma_{\text{noise}}^2$</td><td>$611^{+22}_{-19}$</td><td>$530\pm18$</td><td>$3.0\sigma$</td></tr><tr><td>$l_{\text{mix}}$</td><td>$2.97^{+0.09}_{-0.08}$</td><td>$3.05^{+0.08}_{-0.08}$</td><td>$0.7\sigma$</td></tr><tr><td>$\sigma_{\text{mix}}^2/\sigma_{\text{noise}}^2$</td><td>$50.4^{+2.1}_{-1.9}$</td><td>$48.8^{+2.1}_{-2.0}$</td><td>$0.6\sigma$</td></tr><tr><td>$l_{\text{ex}}$</td><td>$0.26^{+0.01}_{-0.01}$</td><td>$0.32^{+0.01}_{-0.01}$</td><td>$4.2\sigma$</td></tr><tr><td>$\sigma_{\text{ex}}^2/\sigma_{\text{noise}}^2$</td><td>$2.18^{+0.09}_{-0.14}$</td><td>$2.82^{+0.10}_{-0.10}$</td><td>$4.2\sigma$</td></tr><tr><td>$f_{21}$</td><td></td><td></td><td></td></tr><tr><td></td><td>$\eta_{21} = 1/2$</td><td>$x_1 = 0.49^{+1.00}_{-1.00}$</td><td>—</td></tr><tr><td></td><td>$l_{21} > 0.73$</td><td>$x_2 = -0.34^{+0.97}_{-0.97}$</td><td>—</td></tr><tr><td></td><td>$\frac{\sigma_{21}^2}{\sigma_{\text{noise}}^2} < 0.77$</td><td>$\frac{\sigma_{21}^2}{\sigma_{\text{noise}}^2} = 0.01^{+0.09}_{-0.01}$</td><td>—</td></tr></tbody></table>

<figure><img src="image_1.png" /><figcaption>**Figure 1.** Upper limits on the power spectrum of the 21-cm signal from 10 nights of LOFAR observational data at $z \approx 9.1$. Top: individual upper limits case. The upper limits on $\Delta_{21}^2$ (blue dashes with downward arrows) are obtained using the VAE-based kernel separately from the excess noise (dashed-dotted green line). The 2-$\sigma$ confidence of the recovery by the VAE-based kernel is shaded in blue, and noted to effectively extend down to zero. We also plot the upper limits obtained by M20 (crosses with downward arrows), the noise (dashed yellow line), and the 2-$\sigma$ error on the noise (dotted brown). Bottom: joint upper limits case. The upper limits on the power spectrum of the 21-cm signal and excess noise are recovered together (solid blue line and shaded region), as done in M20. The other lines and symbols have the same meaning as in the top panel.</figcaption></figure>

additionally show the 2-$\sigma$ confidence interval on the recovery of the VAE-based kernel, the overall noise (yellow dashed line), and the 2-$\sigma$ error on the noise (brown dotted line). We compare our recovery with the results from M20 shown with magenta crosses. We note that while the upper limits are higher than the 2-$\sigma$ error on the noise, the values approach similar values at the higher k bins (see values listed in Table 2). Further, it is still possible that some fraction of the 21-cm signal may be assigned to the excess noise, and vice-versa. Thus, we also show a second case (bottom panel), which takes a more conservative approach of including the excess noise in the upper limits along with the 21-cm signal (solid blue line), as done in M20. We refer to this as the joint upper limits case.

We note that, in the individual case, the upper limits on the 21-cm signal power spectrum are significantly deeper than those of the combined power spectrum of the excess noise and 21-cm signal found by M20.

Because we depict just the extracted 21-cm signal component, we assume that the bias correction for the excess noise is accurate, with minimal loss of the 21-cm signal. This is an idealized scenario, where the excess noise component is accurately described by the covariance kernel used for it. To verify this, we need a better understanding of the source of the excess noise component. Nevertheless, we report results with our current best choice of the kernel as identified by M20, to test the extent of improvement possible with the inclusion of the VAE-based kernel for the 21-cm signal. However, in Fig. 1, the 2-$\sigma$ confidence interval on the lower side extends to < 1 mK$^2$, and thus is effectively zero. Thus, this confirms that this is just an upper limit and not a detection.

When comparing the joint upper limits cases, we find marginally higher values than those found by M20 for $k \lesssim 0.2 \text{ h cMpc}^{-1}$, and marginally lower for $0.2 \text{ h cMpc}^{-1} \lesssim k \lesssim 0.5 \text{ h cMpc}^{-1}$. This is caused by the power of the excess noise component being mildly higher in the lower k-bins, and mildly lower in the higher k-bins, compared to the results of M20. This shows up as a $4.2\sigma$ difference of the hyperparameters $l_{\text{ex}}$ and $\sigma_{\text{ex}}^2/\sigma_{\text{noise}}^2$ from the results of M20. A closer inspection shows that this occurs because a small fraction of the intrinsic sky component is assigned to the excess noise component at lower k-bins, and vice-versa for the higher k-bins, due to the usage of a more accurate 21-cm kernel. While the difference in the recovered power spectrum for the foregrounds is not significant, the higher contributions at small scales for the excess noise, could allow a better characterisation of the same.

The exact values of the recovered median power spectrum of the 21-cm signal and their upper limits in each k-bin are listed in Table 2 for the joint and individual upper limits case discussed in this work, as well by M20. We also show the 2-$\sigma$ error of the noise (which represents the maximum sensitivity achievable) to compare against the obtained upper limits. We can summarize the results as follows:

(i) The excess noise is confirmed to be the dominant signal component after sky-model and residual foregrounds subtraction, with its recovered power spectrum being about an order of magnitude stronger than the noise.

(ii) While the hyperparameters for the excess noise component differ by $4.2\sigma$ from the results of M20, the recovered power spectrum and 2-$\sigma$ upper limits show only a minor deviation.

(iii) In the case of individual upper limits, the upper limits of the 21-cm signal power spectrum lie above the 2-$\sigma$ error on the noise, and thus are not noise-dominated. However, this is not a detection, because the lower limit on the recovery by the VAE-based kernel is consistent with zero.
<table><caption>**Table 2.** Comparing the results of M20 with those of our joint and individual upper limits of the excess noise component and 21-cm signal power spectra. The first column refers to the *k* bins, and the following to the recovered 21-cm signal power ($\Delta_{21}^2$) and its corresponding 2-$\sigma$ upper limit ($\Delta_{21,UL}^2$) for M20, the joint upper limits (bottom panel in Fig. 1) and the individual upper limits (top panel in Fig. 1) applied in this work. The last column shows the 2-$\sigma$ error on the noise (2-$\sigma$ error<sub>noise</sub>).</caption><thead><tr><th rowspan="2"><i>k</i><br>(<i>h</i> cMpc<sup>-1</sup>)</th><th colspan="2">M20 (<i>f</i><sub>ex</sub> + <i>f</i><sub>21</sub>)</th><th colspan="2">Joint upper limits (<i>f</i><sub>ex</sub> + <i>f</i><sub>21</sub>)</th><th colspan="2">Individual upper limits (<i>f</i><sub>21</sub>)</th><th rowspan="2">2-<i>σ</i> error<sub>noise</sub><br>2Δ<sup>2</sup><sub>noise, err</sub><br>(mK<sup>2</sup>)</th></tr><tr><th>Δ<sup>2</sup><sub>21</sub><br>(mK<sup>2</sup>)</th><th>Δ<sup>2</sup><sub>21,UL</sub><br>(mK<sup>2</sup>)</th><th>Δ<sup>2</sup><sub>21</sub><br>(mK<sup>2</sup>)</th><th>Δ<sup>2</sup><sub>21,UL</sub><br>(mK<sup>2</sup>)</th><th>Δ<sup>2</sup><sub>21</sub><br>(mK<sup>2</sup>)</th><th>Δ<sup>2</sup><sub>21,UL</sub><br>(mK<sup>2</sup>)</th></tr></thead><tbody><tr><td>0.075</td><td>(58.96)<sup>2</sup></td><td>(72.86)<sup>2</sup></td><td>(71.33)<sup>2</sup></td><td>(79.68)<sup>2</sup></td><td>(6.58)<sup>2</sup></td><td>(24.55)<sup>2</sup></td><td>(18.63)<sup>2</sup></td></tr><tr><td>0.100</td><td>(95.21)<sup>2</sup></td><td>(106.65)<sup>2</sup></td><td>(107.05)<sup>2</sup></td><td>(115.77)<sup>2</sup></td><td>(8.56)<sup>2</sup></td><td>(29.43)<sup>2</sup></td><td>(19.05)<sup>2</sup></td></tr><tr><td>0.133</td><td>(142.17)<sup>2</sup></td><td>(153.00)<sup>2</sup></td><td>(158.65)<sup>2</sup></td><td>(167.76)<sup>2</sup></td><td>(10.48)<sup>2</sup></td><td>(37.29)<sup>2</sup></td><td>(23.22)<sup>2</sup></td></tr><tr><td>0.179</td><td>(235.80)<sup>2</sup></td><td>(246.92)<sup>2</sup></td><td>(245.42)<sup>2</sup></td><td>(261.40)<sup>2</sup></td><td>(14.40)<sup>2</sup></td><td>(43.92)<sup>2</sup></td><td>(30.47)<sup>2</sup></td></tr><tr><td>0.238</td><td>(358.95)<sup>2</sup></td><td>(370.18)<sup>2</sup></td><td>(360.08)<sup>2</sup></td><td>(372.28)<sup>2</sup></td><td>(16.29)<sup>2</sup></td><td>(56.32)<sup>2</sup></td><td>(38.70)<sup>2</sup></td></tr><tr><td>0.319</td><td>(505.26)<sup>2</sup></td><td>(520.33)<sup>2</sup></td><td>(472.49)<sup>2</sup></td><td>(490.22)<sup>2</sup></td><td>(19.95)<sup>2</sup></td><td>(73.79)<sup>2</sup></td><td>(58.72)<sup>2</sup></td></tr><tr><td>0.432</td><td>(664.23)<sup>2</sup></td><td>(683.20)<sup>2</sup></td><td>(601.37)<sup>2</sup></td><td>(626.81)<sup>2</sup></td><td>(26.61)<sup>2</sup></td><td>(99.19)<sup>2</sup></td><td>(86.01)<sup>2</sup></td></tr></tbody></table>

We conclude that with these data, performing a bias correction for the excess noise provides an individual upper limit on the 21-cm signal power spectrum that is higher than the 2-$\sigma$ error on the noise. However, the upper limits approach the 2-$\sigma$ error on the noise at higher *k* bins. Thus, lowering the noise threshold with longer observations could allow a greater separation between the upper limits and the noise threshold, and eventually also a detection, provided no correlation between the excess noise and the 21-cm signal is present. For this, the LOFAR EoR KSP team has ≈100 nights of data to analyse and is currently in the process of doing so. We discuss the correlations of the different components and the limitations of the overall model in the subsequent section.

## 4 DISCUSSION

The individual upper limits are greater than the 2-$\sigma$ error on the noise, and thus is not noise dominated. However, it is still possible that the VAE-based kernels provide an incomplete model of the data. In this case, we would be biased towards cases where the model is enough for the 21-cm signal, and the remainder of the data is modelled with the excess and foregrounds components as defined above. For example, better characterisation of the excess noise may indeed lead to broader models for it. In our current model, what we can explore however, is the degeneracy between the VAE-based kernel and the excess noise kernel. A decrease in the correlation between them would indicate that at least the models we use for them are not degenerate. For example, the difference in the intrinsic sky and excess noise component hyperparameters with respect to M20 could be due to a decrease in correlation of the hyperparameters for these components by using the VAE-based kernel for the 21-cm signal. To have a clearer picture of this, it is important to investigate whether the hyperparameters of the different components in our overall model are correlated and whether the overall model is a complete picture of the data.

To assess this, in Fig. 2, we show the corner plot of the hyperparameters, where the purple contours are the 68 per cent, 95 per cent, and 99.7 per cent confidence intervals. We see that hyperparameters of each component are largely uncorrelated with those of other components. In particular, we note that the 21-cm signal hyperparameters are completely uncorrelated from those of the excess noise. This provides greater confidence in the 2-$\sigma$ confidence intervals obtained in the individual upper limits case.

Furthermore, to assess the performance of our model, we plot the residual obtained by subtracting the data cube of the model

<figure><img src="image_2.png" /><figcaption>**Figure 2.** Posterior probability distributions of the GPR model hyperparameters used: $x_1$, $x_2$, the variances $\sigma_{21}^2$, $\sigma_{sky}^2$, $\sigma_{mix}^2$, and the coherence-scales $l_{sky}$, $l_{mix}$, and $l_{ex}$. The purple contours show the 68 per cent, 95 per cent, and 99.7 per cent confidence intervals and the diagonal plots refer to the individual posterior distributions of each of the hyperparameters listed. Note that the hyperparameters of each component are largely uncorrelated with those of other components.</figcaption></figure>

from the observational data cube in Fig. 3, and compare its power spectrum (black solid) to that of the noise data cube (yellow dashed). We find excellent agreement between them within the 2-$\sigma$ error on the noise. To numerically quantify this, we calculate the $\chi^2$ of the residual versus the noise power spectrum, and find this to be ≈6.12. With 7 *k* bins, this gives a reduced $\chi^2$ ≈ 0.87, thus confirming the performance of the model. However, it is still possible that some minor leakage can occur between the different signal components. Thus, better modelling of the excess noise remains necessary to improve the bias correction as more data is added. Furthermore, for an eventual detection, reduction of the noise component through longer observations is also needed.
<figure><img src="image_3.png" /><figcaption>**Figure 3.** Residual (black solid line) obtained by subtracting the full GPR model from the data, and its corresponding 2-σ confidence interval (grey shaded region). The noise is shown as a yellow dashed line with 2-σ error bars. We note an excellent agreement between them.</figcaption></figure>

## 5 CONCLUSIONS

To improve the modelling of the 21-cm signal covariance kernel, Mertens et al. (2024) introduced a Machine Learning based method that uses a VAE based algorithm. This was trained on GRIZZLY simulations (Ghara et al. 2015, 2018, 2020) of the 21-cm signal and tested against mock data sets by A24. In this work, we use the VAE-based kernel trained at $z = 9.16$ to obtain an updated upper limit on the 21-cm signal from 141 h of LOFAR data ($\approx 10$ nights) at $z \approx 9.1$. We compare our results to past efforts, which used analytic functions instead of Machine Learning (M20), and find that they are consistent with errors, although the new upper limits are slightly increased at $k \lesssim 0.2 \ h \ c\mathrm{Mpc}^{-1}$, and mildly decreased for $0.2 \ h \ c\mathrm{Mpc}^{-1} \lesssim k \lesssim 0.5 \ h \ c\mathrm{Mpc}^{-1}$. Overall, we report a 2-σ upper limit of $\Delta_{21}^2 < (80)^2 \mathrm{mK}^2$ at $k = 0.075 \ h \ \mathrm{Mpc}^{-1}$ in the conservative case, when jointly recovering the 21-cm signal and excess noise components.

We also investigate the case in which the upper limits of the 21-cm signal are evaluated including a bias correction for the excess noise component. In this case, we report a 2-σ upper limit of $\Delta_{21}^2 < (25)^2 \mathrm{mK}^2$ at $k = 0.075 \ h \ \mathrm{Mpc}^{-1}$. We confirm that the VAE kernel hyperparameters are uncorrelated from the excess component ones, thus providing confidence for the estimated upper limits. This is strengthened by the fact that the residuals obtained by subtracting the full model from the data are consistent with the noise with $0.87\sigma$. Although the bias correction is promising, we still caution against using the bias-corrected upper limits, as better modelling of the excess noise remains necessary to improve confidence on its separation from the 21-cm signal, as the data improves and the signal to noise increases. In future work, the LOFAR EoR KSP team will be improving the characterisation of the excess noise component and its corresponding covariance kernel, and assess whether an excess noise bias correction can reliably be applied as this work suggests it can. As demonstrated here, this approach would substantially reduce the current upper limits on the 21-cm signal, although this should still be considered as an upper limit and not a detection for the given data, as it is already reaching the sensitivity limit due to the noise, given by the 2-σ error on the noise. This improvement could offer enhanced astrophysical constraints compared to Ghara et al. (2020), allowing for the rejection of a significantly greater number of cold

IGM scenarios and achieving stricter constraints on the sources of X-ray heating. This shall be carried out in future work. Additionally, the team will make use of VAE-based covariance kernels as developed in A24 across multiple redshifts with data from the LOFAR telescope.

## ACKNOWLEDGEMENTS

AA thanks the EoR research group at MPA for helpful discussions. RG acknowledges support from the Kaufman Foundation (gift no. GF01364). LVEK acknowledges the financial support from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (grant agreement no. 884760, 'CoDEX'). SZ thanks Alexander von Humboldt Foundation for the Humboldt Research award and MPA hospitality. Lastly, we would like to thank the referee for their helpful comments that contributed to the improvement of this paper.

## DATA AVAILABILITY

The VAE kernels for $z = 9.16$ can be shared on reasonable request. The PS_EOR package used to run GPR and ML-GPR can be found at https://gitlab.com/flomertens/ps_eor. The LOFAR data analysed in this work can be shared on reasonable request.

## REFERENCES

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

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

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

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

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

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

* 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

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

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

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

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

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

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

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

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

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

* Munshi S. et al., 2024, A&A, 681, A62

* Rasmussen C. E., Williams C. K. I., 2006, Gaussian processes for machine learning. Adaptive computation and machine learning. MIT Press, USA

* 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 New York, NY

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

* Zaroubi S., 2013, Astrophysics and Space Science Library, The Epoch of Reionization. Springer-Verlag, Berlin, 45

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