Detecting ionized bubbles around luminous sources during the reionization era using HI 21-cm signal

Here, we briefly describe simulated ionized regions surrounding a very bright quasar at $z=7.1$, and around a collection of galaxies at $z=8.3$ during the EoR. First, we simulated the dark-matter distribution using the $N$-body code²²2https://github.com/rajeshmondal18/N-body developed by [51]. We simulated this within a cubic box of the comoving size of $\sim 527$ Mpc with $2048^{3}$ grids and the same number of particles. The minimum halo mass resolved in this simulation is ${\sim 6.26\times 10^{9}M_{\odot}}$, corresponding to a minimum of $10$ dark matter particles, each with a mass of ${\sim 6.26\times 10^{8}M_{\odot}}$. We note that our simulations do not include the contributions of low-mass halos in the range $\bm{10^{8}M_{\odot}}$ to $\bm{\sim 6.26\times 10^{9}M_{\odot}}$. In reality, these halos are expected to play a significant role in the reionization of the IGM, potentially giving rise to more diffuse, complex, and smaller ionized structures during the EoR. Then we use a Friends-of-Friends ³³3https://github.com/rajeshmondal18/FoF-Halo-finder (FoF) [52] halo finder to identify the collapsed halos that host the galaxies/bright quasars. These galaxies/quasars serve as the sources of ionizing photons during the EoR. Subsequently, the resolution of simulated cubes is degraded to $512^{3}$ grids. This resolution resembles the corresponding angular resolution of the uGMRT and SKA1-Low and has the added benefit of reduced computational requirements for simulation. We then employ a semi-numerical prescription [53, 54, 55] on these outputs to simulate ionization maps and the HI 21-cm differential brightness temperature maps.

The semi-numerical approach compares the average number of ionizing photons $\langle n_{\gamma}\rangle_{R}$ and neutral hydrogen atoms $\langle n_{\rm H}\rangle_{R}$ within a spherical region of radius $R$ around a grid cell. If the ionization condition $\langle n_{\gamma}\rangle_{R}\geq\langle n_{\rm H}\rangle_{R}$ is met for any radius $R$, then the grid cell is reckoned to be ionized, within that region. Otherwise, an ionized fraction value of $x_{\rm HII}=\langle n_{\gamma}\rangle_{\rm grid}/\langle n_{\rm H}\rangle_{\rm grid}$ is assigned to that grid cell. This procedure is repeated for varying radius $R$ up to a maximum value $R_{\rm max}$, over the whole simulation box. It is assumed that the number of ionizing photons from halos, $N_{\gamma}$ contributed to the ionization of the IGM, follows the relation $N_{\gamma}\propto M_{\rm h}$, where $M_{\rm h}$ is the mass of the corresponding halo.

Here, we consider two scenarios. First, we aim to simulate an ionized region around a very bright quasar at a redshift of $7.1$ similar to the one reported in [8]. In the second case, we simulate an ionized region around a collection of galaxies at redshift $8.3$ as reported in [17]. In the first scenario, we assume that such a very bright quasar is hosted by the most massive halo (of mass $M_{\rm h}\sim 6\times 10^{12},M_{\odot})$ found in our simulation [20]. We further assume that the ionizing photon emission rate by the QSO is $\sim 1.3\times 10^{57}\,{\rm s}^{-1}$ [35, 20]. We note that along with the central quasar, many galaxies reside in and outside the QSO ionized regions. These galaxies are also contributing ionizing photons. The mass-averaged neutral Hydrogen fraction resulting from galaxies outside the QSO HII region is $\langle x_{\rm HI}\rangle\sim 0.88$ at a redshift of $z=7.1$. In figure 1 (top left panel), the resulting HI 21cm map is shown with the ionized bubble around the quasar in the center of the map. The ionized region around the central bright QSO is considerably bigger (of radius around $23.5$ comoving Mpc) than all other ionized regions produced by the galaxies.

In the second case, we consider a collection of galaxies near the most massive halo (of mass $M_{\rm h}\sim 3.5\times 10^{12}M_{\odot}$) in our simulation at $z=8.3$. Simulation specifications such as the cubic box size, total number of particles used, and number of grids along a direction are the same as above. The ionizing luminosity of these galaxies, including the most massive halo in our simulation, is adjusted to create an HII bubble with a radius of $\sim 27.8$ Mpc, and the resulting mass-averaged neutral fraction is around $\langle x_{\rm HI}\rangle\sim 0.94$. The ionized bubble from this scenario is shown in the top right panel of figure 1. We see that the ionized bubble in this scenario is slightly aspherical in shape, mainly because it is created by a collection of sources located at different locations inside the bubble.

In the above two cases, the targeted HII bubbles are isolated and nearly spherical. This geometry arises because the ionizing radiation from other sources is negligible, and the IGM outside the targeted bubbles remains highly neutral. While these scenarios may be suitable at the onset of reionization, they are likely unrealistic for the later stages of the process. At later stages of reionization, ionized bubbles around very bright sources are surrounded by neighboring ionized regions, leading to an aspherical shape for the targeted bubble. The third scenario accounts for this complexity. The bottom panel of Fig. 1 displays a 2D map of the differential brightness temperature (${\delta T_{b}}$) at redshift $z=7.1$, where the central targeted HII bubble appears highly aspherical. The quasar at the center is the same as the first case. Surrounding it there are numerous other HII regions, unlike the first case. This results in a much lower mean neutral fraction, which is $0.52$ in this scenario.

We want to mention the considerable uncertainty regarding the ionizing photon emissivity rate from QSO, and therefore, the photon emissivity rate and the resulting QSO HII region may vary. The consequences of this uncertainty will be discussed later in the paper.

The astrophysical foregrounds that contaminate the EoR HI 21-cm signal observations comprise various components [56], including diffuse galactic synchrotron emission, radio synchrotron emission from point-like sources, free-free radio emission from diffuse ionized gas, etc.

The diffuse galactic synchrotron emission (DSGE) is one of the most dominant foreground components. We model the angular power spectrum corresponding to the diffuse galactic synchrotron emission using a power-law formulation [57] as,

$${C_{l}}^{M}(\nu)=A_{150}\times{\left(\frac{1000}{l}\right)}^{\beta}{\left(\frac{\nu_{0}}{\nu}\right)^{2\alpha}}.$$

(3.1)

Here, $\nu_{0}=150$ MHz and $\nu=\nu_{0}+\Delta\nu$. The power-law index, $\beta$, is set to $2.34$, and the amplitude $A_{150}=513\,\text{mK}^{2}$ is based on previous work [58]. The parameter $\alpha$ describes how the power spectrum changes with frequency, which is set to $2.8$ [29, 30]. This simplified model allows us to simulate DGSE maps and to study their impact on the detection of individual ionized bubbles in the HI 21-cm maps. We use the following equation to generate Fourier components of the DGSE on each gridded baseline $U$, which is given by

$$\Delta\tilde{T}(U)=\sqrt{\frac{\Omega{C_{l}}^{M}}{2}}[x(U)+iy(U)].$$

(3.2)

Here, $\Omega$ denotes the total solid angle corresponding to the spatial size of simulated maps. $x(U)$ and $y(U)$ are independent Gaussian random variables with zero mean and unit variance. Upon simulating 2D DGSE maps in Fourier space, we apply the inverse Fourier transform on them to obtain the simulated sky map $\delta T_{b}(\bm{\theta})$ at a specific frequency. Subsequently, we use the frequency scaling to obtain the DGSE map at other desired frequencies. Figure 2 shows a single realization of the 2D simulated DGSE map centered at frequency $\nu=175$ MHz corresponding to redshift $z=7.1$. As mentioned earlier, the angular resolution of these maps is $23^{\prime\prime}$, and the strength of temperature fluctuations changes with angular resolution. We note that the mean temperature of the DGSE map is set to zero, as the radio interferometric experiments, such as the uGMRT/SKA1-Low, are not sensitive to the mean due to the lack of the zero spacing baseline. Subsequently, we obtain mock sky maps at multiple frequency channels by combining the HI 21-cm and the foreground maps simulated above.

In radio interferometric observations, the measured quantity is the visibility, denoted as $V(\bm{U},\nu)$, which can be written as the Fourier transform of the product of the specific intensity pattern of the sky, $\delta I_{\nu}(\bm{\theta})$ and $A(\bm{\theta},\nu)$. It is given as,

$$V(\bm{U},\nu)=\int d^{2}\theta\,A(\bm{\theta},\nu)\,\delta I_{\nu}(\bm{\theta})\,e^{2\pi i\bm{\theta}.\bm{U}}.$$

(3.4)

Here, the baseline $\bm{U}=\bm{d}/\lambda$ denotes the antenna separation $\bm{d}$ projected on the plane perpendicular to the line of sight in units of the observed wavelength ($\lambda$) of the signal. $\bm{\theta}$ is a two-dimensional vector on the plane of the sky with the origin at the center of the field of view (FoV), and $A(\bm{\theta},\nu)$ is the beam pattern of the individual antenna. Assuming uniformly illuminated circular aperture of diameter $D$, we model the antenna primary beam pattern as [57],

$$A(\bm{\theta},\nu)=\left[\left(\frac{2\lambda}{\pi\theta D}\right)J_{1}\left(\frac{\pi\theta D}{\lambda}\right)\right]^{2},$$

(3.5)

where $J_{1}$ is the first order Bessel function. The antenna diameter $D=45$m and $35$m for the uGMRT and SKA1-Low respectively. The specific intensity pattern of the sky, $\delta I_{\nu}(\bm{\theta})\equiv\delta T_{\rm b}({\bm{n}},z)$ and can be written as,

$$\delta I_{\nu}({\bm{n}},z)=\frac{2k_{B}\nu^{2}}{c^{2}}{\rm\delta}T_{\rm b}({\bm{n}},z).$$

(3.6)

Considering a spherical ionized bubble surrounded by a uniform HI background at the center of the field of view, we can derive the analytical formula for the visibility of the HI 21-cm signal at a particular frequency channel $\nu$ and baseline $\bm{U}$ following [33] as,

$$S_{\rm center}(\bm{U},\nu)=-\pi\bar{I_{\nu}}x_{\rm HI}{\theta_{\nu}}^{2}\left[\frac{2J_{1}(2\pi U\theta_{\nu})}{2\pi U\theta_{\nu}}\right]\Theta\left(1-\frac{\left|\nu-\nu_{c}\right|}{\Delta\nu_{b}}\right),$$

(3.7)

where

$$\bar{I_{\nu}}=2.5\times 10^{2}{\rm Jy/sr}\left(\frac{\Omega_{b}h^{2}}{0.02}\right)\left(\frac{0.7}{h}\right)\left(\frac{H_{0}}{H(z)}\right)$$

(3.8)

is the specific intensity from uniform HI background at redshift $z$, $x_{\rm HI}$ is the mass-averaged neutral hydrogen fraction, $\theta_{R}$ is the angular size of the ionized bubble corresponding to comoving radius $R$, $J_{1}$ is first order Bessel function and $\Delta\nu_{b}$ is ionized bubble size in frequency. Note that a spherical ionization bubble will appear as a circular one at a particular frequency channel. We use this analytical equation to design a filter for detecting the ionized bubbles in the mock HI 21-cm maps obtained from simulations.

To simulate mock uGMRT and SKA1-Low observations, we consider baseline distributions corresponding to both instruments. There are currently $30$ active antennas at uGMRT [59] and the proposed number of antennas in SKA1-Low is $512$ [60]. Table 1 summarises instrument specifications. Figure 3 shows the baseline coverage for an $8$-hour uGMRT and SKA1-Low observations at the $175$ MHz frequency corresponding to $z=7.1$ for the HI 21-cm signal. It also shows the baseline density i.e., the total number of baselines per unit area of circular annulus corresponding to baseline $U$ and $U+dU$. For the uGMRT, baselines vary from $32\lambda$ to $13500\lambda$ and for the SKA1-Low, the baseline range extends from $4\lambda$ to $37000\lambda$ over the same $8$ hour, $175$ MHz observations. As expected, the number density of baselines and their extent are much greater for the SKA1-Low, as it has more antennae, and some of them are placed at large distances compared to the uGMRT. The baseline distribution at $153$ MHz (corresponding to $z=8.3$) is very similar to that at $175$ MHz but will be reduced (refer to Table 1).

We generate the mock visibility data of the HI 21-cm signal and the foregrounds (DGSE and extragalactic point sources) at all simulated baselines at $256$ frequency channels, centered around frequencies $\nu=175$ MHz and $153$ MHz, using the Discrete Fourier Transformation (DFT) of the total simulated maps. For simplicity, we keep the baseline distribution the same in all frequency channels.

Figure 4 shows simulated visibilities corresponding to the 2D central HI 21-cm map around the ionized bubble centered at $z=7.1$. The top left and top right panels show plots of the signal as a function of baseline and frequency respectively, for the uGMRT and the bottom panels display results for the SKA1-Low. We also compare the simulated visibility with the analytical estimation assuming a perfectly spherical ionized bubble of comoving radius $23$ Mpc and surrounded by uniform and fully neutral medium (refer to equation 3.7). The IGM around the ionized bubble in simulated maps is not uniform, but traces the underlying DM distribution and is also filled with many small ionized bubbles. This causes HI density to fluctuate and consequently, the visibility signal also fluctuates as a function of baseline and frequency. We see in Figure 4 that the HI fluctuations are stronger at lower baselines which reflects the fact that HI power spectrum $P_{\rm HI}(k)$ (or equivalently the visibility-visibility correlations) increases at lower $k$-modes or baselines [61, 34]. We further see that the fluctuating HI outside the bubble dominates over the visibilities solely due to the HI 21-cm signal around a spherical ionized bubble which is buried in a uniform HI background. Visibilities are similar for the second scenario at $z=8.3$ for both experiments and we do not show them here explicitly.

To make our study more realistic, we introduce system noise contribution from radio interferometers in our simulated mock visibility data. The system noise contribution to the measured visibility in each baseline and frequency channel is an independent Gaussian random variable with zero mean and the root mean square [33]

$$\sqrt{\langle N^{2}\rangle}=\frac{\sqrt{2}k_{B}T_{\rm sys}}{A_{\rm eff}\sqrt{\Delta\nu\Delta t}},$$

(3.9)

where $T_{\rm sys}$ represents the total system temperature, $A_{\rm eff}$ is the effective collecting area of the antenna, and $\Delta\nu$ and $\Delta t$ denote the channel width and correlator integration time, respectively.

For the uGMRT, $T_{\rm sys}=300$ K and $A_{\rm eff}/(2k_{\rm B})=0.33$ K/Jy [59] at frequency $\nu\sim 175$ MHz. The estimated root mean square noise in visibility is approximately $0.64\,{\rm Jy}$ for frequency channel width ($\Delta\nu$) of $62.5$ kHz and integration time ($\Delta t$) of $16$ seconds in each baseline. $T_{\rm sys}=300$ K and $A_{\rm eff}/(2k_{\rm B})=0.348$ K/Jy for the SKA1-Low at the same frequency considering antenna efficiency factor unity. $T_{\rm sys}$ and $A_{\rm eff}/(2k_{\rm B})$ change to $400$ K and $0.348$ K/Jy for the SKA1-Low at $153$ MHz frequency. For the uGMRT, these numbers are $400$ K and $0.33$ K/Jy [59] respectively at $153$ MHz frequency. We note that the rms noise will be reduced by a factor of $\sqrt{t_{\rm obs}/\Delta t}$ where $t_{\rm obs}$ is the total observation time.

Our mock simulation considers $8$ hrs of observations with $16$s and $320$s integration time for the uGMRT and SKA1-Low respectively. This produces a total of $783000$ and $11773440$ baselines for the uGMRT and SKA1-Low respectively. Lower integration time will produce a higher number of baselines for the SKA1-Low and hence will increase the requirement of computational resources which is beyond our reach. The system noise contribution for $8$ hrs of observations will be too high to detect the HI 21-cm signal around the individual ionized bubbles. Therefore, we present our results for $2048$ and $3072$ hrs of uGMRT observations at redshifts $7.1$ and $8.3$, respectively, assuming a total of $256$ and $384$ nights of observations, with each night lasting $8$ hours. Further, for the third scenario, we present the results at redshift $7.1$ for $5120$ hours, assuming $640$ nights of observations. The increased observing time for the third scenario is due to the lower mean neutral hydrogen fraction outside the targeted bubble. Assuming all nights are the same, we reduce the system noise corresponding to an $8$ hrs observation by a factor of $16$, $19.6$, and $25.29$ for the three scenarios, respectively. The resulting mock signal becomes equivalent to $2048$, $3072$, and $5120$ hours of uGMRT observations. Similarly, for the SKA1-Low we present our results for $96$ hrs of observation for the first two cases at redshift $7.1$ and $8.3$ (neutral hydrogen fraction ${\sim 0.9}$) and for $200$ hrs of observation for the third case at redshift $7.1$ (neutral hydrogen fraction ${\sim 0.5}$). Eventually, we obtain mock visibility data with all the signal contributions (HI 21-cm signal, galactic foreground, extragalactic point sources, and system noise contribution) typical of the uGMRT and SKA1-Low observations.

Figure 5 shows the mock visibility data comprising of HI signal, diffused galactic foreground contribution, and system noise as a function of frequency for two baselines $U=34\lambda$ (top left panel) and $U=230\lambda$ (top right panel) for the uGMRT at $175$ MHz. The corresponding mock visibility for the SKA1-Low is shown in the bottom panels (bottom left panel: $U=11\lambda$, and bottom right panel: $U=254\lambda$). They also show diffused galactic foreground contributions (red dashed lines). We see that the diffused galactic foreground dominates over the HI 21-cm signal and system noise at a small baseline (left panel) and, at a large baseline (right panel) the system noise dominates over the other two components. The HI 21-cm signal remains buried under the stronger foreground and system noise contributions.

Figure 6 shows the total mock visibility comprising of HI signal, diffused galactic synchrotron emission, radiation from extra-galactic radio point sources, and system noise as a function of frequency for two baselines ${U=34\lambda}$ (left panel) and ${U=230\lambda}$ (right panel) for the uGMRT at $175$ MHz. We note that Figure 6 shows contribution both from the extragalactic point sources, galactic diffuse synchrotron radiation along with other components unlike Fig. 5 which shows galactic diffuse synchrotron emission only as the foreground component. The radio point sources are modeled using a simplistic power-law form as discussed in subsection 3.2.1. The green lines indicate visibilities due to point sources. We see that at a small baseline (left panel) the diffuse galactic synchrotron emission dominates over the point sources; at a large baseline (right panel), the point sources dominate over the diffuse galactic synchrotron emission. The total mock visibilities are similar for the second scenario at ${z=8.3}$ for both experiments, except for the fact that the system noise and foreground contributions are higher. Therefore, we do not show them here explicitly. The subsequent section discusses the subtraction of the foreground contribution.

Here, we would like to note that the visibility data used in this study does not contain calibration errors, ionospheric effects, and radio-frequency interference (RFI). A comprehensive investigation of these factors will be addressed in future work.

The foreground contributions are anticipated to be several orders of magnitude higher than the HI 21-cm signal around ionized bubbles during the EoR. Therefore, one significant challenge is reliably subtracting the foreground signal from the observed signal. All foreground mitigation techniques use the fact that the foregrounds are smooth along the frequency [62, 63]. As discussed earlier, the foreground signal in our simulated mock data has also been assumed to be smooth along frequency. To mitigate foreground components, we fit the total mock visibility signal with a $3^{rd}$ order polynomial of observing frequency $\nu$ as follows,

$$f(\nu)=a\nu^{3}+b\nu^{2}+c\nu+d,$$

(4.1)

where parameters $a$, $b$, $c$, and $d$ are constant. Then, we find the best-fit polynomial along each line of sight corresponding to a baseline and subsequently subtract them from the total simulated mock visibility signal.

Figure 7 presents the residual visibility at 175 MHz after the foreground contribution is subtracted from the total mock visibility for the uGMRT (top panels) and the SKA1-Low (bottom panels). Here, we have only considered the diffuse galactic synchrotron emission as the foreground contamination. However, we get a very similar residual signal when we perform similar analysis on the total mock visibilities considering the radio point source contribution, which has been shown in Figure 9. The residual signal contains the HI 21-cm signal, the system noise, and the un-subtracted foreground, if any. To assess the reliability of foreground subtraction, we compare the combined input HI 21-cm signal and system noise with the extracted signal for two different baselines ${U=34\lambda}$ (top left panel) and ${230\lambda}$ (top right panel), as depicted in Figure 8. In the bottom panels of Figure 8, we show a similar comparison for SKA1-Low for two different baselines ${U=11\lambda}$ (bottom left panel) and ${U=254\lambda}$ (bottom right panel). The plots show that the residual visibility closely matches the input visibility data, with only minor deviations. These slight deviations are possibly due to some residual foreground contribution and subtraction of the HI 21-cm signal and the system noise during the foreground subtraction. Consequently, we can conclude that we have mostly subtracted the foreground data from the mock visibility data, making it ready for further analysis, such as matched filtering. We do not present similar results for the second scenario at $z=8.3$ as it is similar to the first scenario.

The residual mock visibility (see Figure 7) contains mainly the HI 21-cm signal around the ionized bubble and the system noise from radio interferometers. We see that the dominance of system noise over the HI 21-cm signal is so pronounced that the latter’s presence is barely discernible. Here, we employ a matched filter technique to maximize the signal-to-noise ratio, thereby enhancing the prospects of detecting the HI 21-cm signal around an ionized bubble. The functional form of the visibility signal corresponding to neutral hydrogen around a spherical ionized bubble is quite well known (refer to eq. 3.7). Additionally, the system noise is anticipated to follow a Gaussian random distribution. These conditions enable us to implement the matched filter-based technique effectively. Below, we briefly summarize the method proposed in [33] and consider the following estimator for detecting ionized bubbles from the HI 21-cm signal, which is given as

$${\hat{E}}=\left[{\sum_{a,b}{S_{f}}^{*}(\bm{U}_{a},\nu_{b})\hat{V}(\bm{U}_{a},\nu_{b})}\right]/\left[{\sum_{a,b}1}\right].$$

(4.2)

Here, $S_{f}(\bm{U},\nu)$ is the filter which is a function of both baseline $\bm{U}$ and frequency $\nu$. $\hat{V}(\bm{U}_{a},\nu_{b})$ is the total observed visibility which consists of the desired HI 21-cm signal, residual foregrounds, system noise, etc. Here, $\bm{U}_{a}$ and $\nu_{b}$ refer to the different baselines and frequency channels in our observations. Here, we note that foreground components are subtracted from the total observed visibilities and we employ the matched filtering technique only on the residual visibility which contains only the HI 21-cm signal and system noise. Figure 3 shows a typical uGMRT and SKA1-Low baseline distribution for $8$ hrs of observations. We note that the estimator is calculated by taking the sum over all baselines and frequency channels. We get one estimator for a given filter. The contribution of the system noise to the estimator is expected to be zero if it is averaged over a large number of independent realizations. However, for a single realization of the system noise, the contribution of it is unlikely to be exactly zero. It would change for different realizations of the system noise. Therefore, we estimate the variance of the estimator due to the system noise using the following equation [33],

$$\left<(\Delta{\hat{E}})^{2}\right>_{NS}={\langle N^{2}\rangle}\left[\sum_{a,b}\left|S_{f}(\bm{U}_{a},\nu_{b})\right|^{2}\right]/\left[\sum_{a,b}\right]^{2}.$$

(4.3)

Here, $\langle N^{2}\rangle$ is the variance of the system noise for single visibility and can be calculated using eq. 3.9. The symbol $\sum_{a,b}$ denotes summation over all baselines and frequency channels. We note that the variance to the estimator depends on the chosen filter, system noise, and baseline distribution but not on the desired signal. The signal-to-noise ratio (SNR) of a particular detection (or non-detection) can be calculated as,

$${\rm SNR}=\frac{\langle\hat{E}\rangle}{\sqrt{\left<(\Delta{\hat{E}})^{2}\right>_{\rm NS}}}.$$

(4.4)

In this work, we consider ${\rm SNR}\gtrsim 5$ as the threshold for detection. We would like to mention that the SNR reaches its maximum when the selected filter exactly (or closely) matches with the HI 21-cm signal which is contaminated by the system noise. Therefore, the chosen filter has the same functional form as the HI 21-cm signal around an ionized bubble (refer to eq. 3.7). It also has one free parameter i.e., the radius of the ionized bubble $R_{\rm b}$. We assume that the location of the bubble is known a priori from other observations. To find a filter that matches the targeted signal, one needs to explore the entire possible ranges of all parameters (here the size of the targeted bubble). The particular set of parameters yielding the maximum SNR will reveal the size of the ionized bubble in the observed HI 21-cm signal.

Here, we present our findings regarding the detectability of ionized bubbles around bright QSO/galaxies in HI 21-cm maps, employing the matched filter technique. We assume that the targeted ionized bubble is located at the antenna phase center. Thus, we simulate the targeted ionized bubble at the center of the simulation cube. Subsequently, we apply the matched filtering technique to the residual visibility obtained after subtracting foreground contaminants. As explained above the residual visibility contains the HI 21-cm signal, system noise, and unsubtracted foregrounds. First, we calculate the estimator ${\hat{E}}$ using equation 4.2. Because the actual size of the targeted ionized bubble is not known in real observations we calculate the estimator for various filter sizes ranging from $5$ to $35$ Mpc. Subsequently, we estimate the variance of the estimator $\left<(\Delta{\hat{E}})^{2}\right>_{\rm NS}$ due to the system noise using equation 3.9 and the signal to noise ratio (SNR, refer to eq. 4.4) for the same range of filter sizes. Figure 10 shows the resulting SNR as a function of filter size $R_{\rm f}$ for $2048$ hrs of uGMRT observation (left panel) and $96$ hrs of SKA1-Low observations (right panel) at $z=7.1$.

The ten different curves in Figure 10 represent the SNR for ten independent realizations of the system noise. We see that for both the uGMRT and the SKA1-Low, the SNR peaks at filter size $\sim 23$ Mpc, with minor deviations in different noise realizations. As per the principle of the matched filter technique, the SNR peaks when the filter matches the actual targeted signal (i.e., the HI 21-cm signal around the ionized bubble) in the simulated HI maps. Therefore, we can conclude that the extracted radius of the ionized bubble in the simulated HI 21-cm map is around $\sim 23$ Mpc. We want to note that the actual radius of the ionized bubble measured directly from the simulated map is $23.5$ Mpc, which is very close to the extracted ones using the matched filter technique. The location of the peak in the SNR curve shifts around the mean value for different realizations, which we shall discuss in detail in the next section.

Further, we see that peak SNR values in all the ten realizations are $\gtrsim 4$ for $2048$ hrs of observations with the uGMRT and $\sim 10$ for $96$ hrs of observations with the SKA1-Low. We assume a total bandwidth of $16$ MHz and $6.25$ MHz while subtracting the foreground for the uGMRT and SKA1-Low, respectively. The choice of smaller bandwidth for the SKA1-Low is solely due to the limited computational resources available. Use of the higher bandwidth for the SKA1-Low would have subtracted the HI 21-cm signal less during the foreground subtraction process, thus increasing the SNR. Subsection 5.4 discusses the impact of foreground subtraction in more detail. The higher SNR obtained in our results indicates that it is possible to significantly detect individual ionized bubbles in the HI 21-cm maps using the uGMRT with $2048$ hrs observations and using the SKA1-Low with only around $\sim 100$ hrs of observations.

To assess the significance of the above result, we estimate the SNR for $100$ independent realizations of the system noise for the uGMRT and $20$ independent realizations for SKA1-Low, while keeping the HI 21-cm signal and foreground components the same. Figure 10 shows the SNR for some of these realizations for the uGMRT (left panel) and SKA1-Low (right panel) at $z=7.1$. The left panel of Figure 11 shows the corresponding mean SNR and $1-\sigma$ spread for all $100$ realizations of the uGMRT, while the right panel shows the same for the SKA1-Low with $20$ realizations. We see that all the realizations show a peak in the SNR. In principle, the peak should appear at the actual ionized bubble size of $23.5$ Mpc. However, the peak location varies in the range of $18-30$ Mpc and $22-24$ Mpc for the uGMRT and SKA1-Low, respectively, for different realizations. It can be seen in Figure 12, which shows a histogram of peak locations for $100$ realizations for the uGMRT.

We see that, although in most cases, the peak is in the range of $22-26$ Mpc for the uGMRT, some extreme realizations show peaks at higher/lower filter radii. The peak value of the SNR for the uGMRT also changes in the range of $\sim 3-7$. Due to much lower system noise and higher sensitivity, the SKA1-Low should be able to detect the ionized bubble with reduced observing time and higher SNR. The estimated bubble size is also highly accurate for the SKA1-Low, which we do not explicitly show here. We note that the above analysis considers only the DGSE as the foreground contaminant. However, we have repeated the analysis including contributions from both the DGSE and radio point sources as foregrounds. We find very similar results to those presented here. Therefore, we do not explicitly show the outcomes of the combined foreground contamination here.

Unlike the ionized bubble at $z=7.1$, mainly created by a bright QSO, the bubble at $z=8.3$ is produced by a collection of galaxies as indicated in [17]. This ionized bubble is relatively more aspherical in shape than at $z=7.1$. Figure 13 presents the SNR as a function of filter size $R_{\rm f}$ for 3072 hrs of uGMRT (left panel) and $96$ hrs of SKA1-Low (right panel) observations at $z=8.3$.

Different curves are for different independent noise realizations. The system temperature is higher at $z=8.3$ compared to $z=7.1$. Therefore, the SNR decreases at a higher redshift for the same observation time. We see that for both the uGMRT and the SKA1-Low, the SNR peaks at a filter size $\sim 28$ Mpc, with minor deviations in different noise realizations. It is very close to the actual radius of the ionized bubble measured directly from the simulated map, which is $27.8$ Mpc.

Further, we see that peak SNR values in all ten realizations are $\gtrsim 4$ for the uGMRT and $\sim 8$ for the SKA1-Low. Like the previous case, we assume a total bandwidth of $16$ MHz and $6.25$ MHz while subtracting the foreground for the uGMRT and SKA1-Low, respectively. We see that it is possible to reliably detect individual ionized bubbles in the HI 21-cm maps using the uGMRT with $3072$ hrs of observation time and the SKA1-Low with only $\sim 100$ hrs of observations at $z=8.3$.

The left panel of Figure 14 shows the mean SNR along with the $1-\sigma$ spread for $100$ noise realizations of the uGMRT, while the right panel shows the same for the SKA1-Low with $20$ realizations at $z=8.3$. Further, the peak location varies between $24$ and $36$ Mpc for the uGMRT for different realizations (Figure 13), which can be seen in Figure 15, which shows a histogram of peak locations for $100$ realizations for the uGMRT. The peak value of the SNR for the uGMRT also changes in the range of $\sim 4-6$ and $6.5-8.5$ for the SKA1-Low.

Here, we present our findings for the third case, where the targeted ionized bubble is surrounded by other ionized regions which result in a lower neutral hydrogen fraction of ${x_{\rm HI}=0.52}$. Further, the targeted bubble at the center also became aspherical. A lower neutral hydrogen fraction means that the background HI 21-cm signal is weaker, which results in a lower SNR value. We use the same matched filtering technique as described earlier but increased the observation time to $5120$ hours for the uGMRT and $200$ hours for the SKA1-Low, considering the lower SNR and increased difficulty in detecting ionized regions at a lower neutral fraction. Figure 16 shows the SNR as a function of filter size ${R_{\rm f}}$ for $5120$ hours of uGMRT observation (left panel) and $200$ hours of SKA1-Low observations (right panel) at $z=7.1$. As in previous cases, each curve represents a different noise realization. We see that for the uGMRT, the SNR peaks at a filter size ${\sim 25}$ Mpc, and for the SKA1-Low, it peaks at ${\sim 23}$ Mpc, with minor deviations in different noise realizations.

Further, we see that peak SNR values in all ten realizations are ${\gtrsim 3}$ for the uGMRT and ${\sim 11}$ for the SKA1-Low. Like the previous cases, we assume a total bandwidth of $16$ MHz and $6.25$ MHz while subtracting the foreground for the uGMRT and SKA1-Low, respectively. We see that it is possible to reliably detect individual ionized bubbles in the HI 21-cm maps using the uGMRT with $5120$ hours of observation time and the SKA1-Low with only $\sim 200$ hours of observations at $z=7.1$ with significantly much lower neutral hydrogen fraction value.

The left panel of Figure 17 shows the mean SNR along with the ${1-\sigma}$ spread for $100$ noise realizations of the uGMRT, while the right panel shows the same for the SKA1-Low with $20$ realizations at $z=7.1$ for ${x_{\rm HI}=0.52}$. Further, the peak location varies between $16$ and $32$ Mpc for the uGMRT for different realizations (Figure 16), which can be seen in Figure 18, that shows a histogram of peak locations for $100$ realizations for the uGMRT. The peak value of the SNR for the uGMRT also changes in the range of $\sim 3-5$ and $10.5-12.5$ for the SKA1-Low.

We noticed that during the process of foreground subtraction using a polynomial along the frequency axis, a small part of the HI 21-cm signal also gets subtracted along with the foreground. It reduces the SNR to some extent. We further notice that the amount of subtraction of the HI 21-cm signal depends on the frequency bandwidth used during the foreground subtraction process. The fraction of HI 21-cm signal subtracted is more for smaller bandwidth compared to the case when larger frequency bandwidth is used. Consequently, the drop in the SNR is more significant for smaller bandwidths. Figure 19 shows the SNR as a function of the bandwidth for $2048$ hrs of observations with uGMRT at $z=7.1$. We see that the SNR, which is lower than $2$ at $6$ MHz bandwidth, increases to $4.5$ when $16$ MHz bandwidth is assumed. Here, we note that the SNR reaches $6.5$ for the above observation when no foreground is used and the process of foreground subtraction is skipped. Therefore, the SNR is reduced by $30$ percent when the bubble size is $\sim 23.5$ Mpc, and $16$ MHz bandwidth is assumed for foreground subtraction. Using higher bandwidth during the foreground subtraction would enhance the SNR and prospects of ionized bubble detection. However, increased bandwidth will increase the computational cost during the analysis. We find similar results for the SKA1-Low but do not explicitly show them here. The impact of foreground subtraction on the SNR will also depend on the foreground subtraction methods, and therefore, a thorough investigation is required.

We see that the visibility corresponding to the HI 21-cm signal around a spherical ionized bubble in the IGM is proportional to the mass-averaged neutral hydrogen fraction $x_{\rm HI}$ of the IGM outside the bubble (refer to eq. 3.7). One can further show that [36],

$${\rm SNR}=Cx_{\rm HI},$$

(5.1)

where $C$ depends on the bubble size, redshift, and various observational parameters such as the $T_{\rm sys}$, total number of antennae, observing time, antenna collecting area, and baseline distribution. All these quantities are constant at a fixed redshift. One can use this relation to place constraints on $x_{\rm HI}$. Now, we focus on the ionized bubble around the QSO at $z=7.1$ and $2048$ hrs of observation time with uGMRT. The SNR value changes for different noise realizations in Figure 10 (left panel). The average SNR value peaks at $4.6$ (Figure 11) at the filter size 23 Mpc for the uGMRT at $z=7.1$. Further, we see that $1\sigma$ uncertainty in the SNR value at the same filter size is 0.73. Therefore, we estimate the hydrogen neutral fraction using $x_{\rm HI}=({\rm SNR}\pm 1\sigma_{\rm SNR})/C$, where the SNR and $1\sigma_{\rm SNR}$ is measured at the peak location.

To calculate $C$, we simulate an entirely neutral and uniform HI 21-cm signal around a spherical ionized bubble of radius 23 Mpc, i.e., $x_{\rm HI}=1$ outside the ionized bubble. It allows us to estimate $C$ using eq. 5.1. Now, after adding the foreground component only with the HI 21-cm signal, we follow the entire analysis and calculate the SNR as a function of filter size without adding system noise from the radio interferometers. The SNR, in this case, peaks at 23 Mpc and is supposed to be equal to $C$ (following eq. 5.1) as the $x_{\rm HI}=1$ outside the ionized bubble in this case. We use this value of $C$ to estimate the mean $x_{\rm HI}$ and the $1\sigma$ uncertainty associated with it. The estimated neutral fraction is $x_{\rm HI}=0.93\pm 0.15$ at $z=7.1$ for $2048$ hrs of observations with uGMRT. It is very close to the actual $x_{\rm HI}$ in the simulation, which is $0.88$. The estimated $1\sigma$ upper limit of $x_{\rm HI}$ exceeds unity which is, of course, unrealistic. The SKA1-Low places much tighter constraints, which is $x_{\rm HI}=0.92\pm 0.055$ using only $96$ hrs of observations for the same scenario. Similarly we also studied this for the third case, where the neutral hydrogen fraction is lower in the simulation, discussed in subsection 5.3. We find the estimated neutral hydrogen fraction to be ${x_{\rm HI}=0.49\pm 0.09}$ at $z=7.1$ for $5120$ hours of observations with the uGMRT, which is very close to the actual neutral hydrogen fraction ${x_{\rm HI}=0.52}$.

The above result demonstrates the prospects of constraining the mass-averaged neutral fraction using the technique discussed in this work. The SKA1-Low is expected to put much tighter constraints on the $x_{\rm HI}$. A thorough analysis is required considering different bubble sizes and hydrogen neutral fractions at different redshifts for various instruments, and we leave this exercise for future work.