IDENTIFYING IONIZED REGIONS IN NOISY REDSHIFTED 21 cm DATA SETS

First, let us describe the underlying 21 cm signal and our reionization simulations. The 21 cm signal will be measured through its contrast with the CMB. The brightness temperature contrast between the CMB and the 21 cm line from a neutral hydrogen cloud with neutral fraction x_HI and fractional baryon overdensity δ_ρ is (Zaldarriaga et al. 2004)

$$\begin{equation} \delta T_{\rm{b}} = 28 x_{{\rm H\scriptsize{I}}}(1 + \delta _{\rho })\left(\frac{T_{\rm{S}} - T_{\gamma }}{T_{\rm{S}}} \right) \left(\frac{1 + z}{10} \right)^{1/2}\ {\rm mK}. \end{equation} \tag{ 1 }$$

Here, T_γ denotes the CMB temperature and T_S is the spin temperature of the 21 cm line. Here and throughout we neglect effects from peculiar velocities, which should be a good approximation at the redshifts and neutral fractions of interest (e.g., Mesinger & Furlanetto 2007; Mao et al. 2012). Furthermore, throughout we assume that the spin temperature is globally much larger than the CMB temperature, i.e., we assume that T_S ≫ T_γ. In this case, the 21 cm signal appears in emission and the brightness temperature contrast is independent of T_S. This is expected to be a good approximation for the volume-averaged ionized fractions of interest for our present study, although it will break down at earlier times (e.g., Ciardi & Madau 2003). With these approximations,

$$\begin{equation} \delta T_{\rm{b}} = T_{0} x_{{\rm H\scriptsize{I}}}(1 + \delta _{\rho }), \end{equation} \tag{ 2 }$$

where T₀ = 28[(1 + z)/10]^1/2 mK. Throughout this paper, we refer to the brightness temperature contrast in units of T₀.

In order to simulate reionization we use the "semi-numeric" scheme described in Zahn et al. (2007; see also, e.g., Mesinger et al. 2010 for related work and extensions to this technique). This scheme is essentially a Monte Carlo implementation of the analytic model of Furlanetto et al. (2004b), which is in turn based on the excursion set formalism. The Zahn et al. (2007) algorithm allows us to rapidly generate realizations of the ionization field over large simulation volumes at various stages of the reionization process. The results of these calculations agree well with more detailed simulations of reionization on large scales (Zahn et al. 2007, 2010).

We start by generating a realization of the linear density field in a simulation box with a comoving side length of 1 h⁻¹ Gpc and 512³ grid cells. The ionization field, x_i, is generated following the algorithm of Zahn et al. (2007), assuming a minimum host halo mass of M_min = 10⁸ M_☉, comparable to the atomic cooling mass at these redshifts (Barkana & Loeb 2001). Each halo above M_min is assumed to host an ionizing source, and the ionizing efficiency of each galaxy is taken to be independent of halo mass. In our fiducial model, we adjust the ionizing efficiency so that the volume-averaged ionization fraction is 〈x_i〉 = 0.79 at z_fid = 6.9. We focus most of our analysis on this redshift and on this particular model for the volume-averaged ionized fraction. However, we consider additional redshifts in Section 5, as well as variations around our fiducial ionization history in an effort to bracket current uncertainties in the ionization history (see, e.g., Kuhlen & Faucher-Giguère 2012; Zahn et al. 2012).

From the linear density field and the ionization field, we generate the 21 cm brightness temperature contrast following Equation (2). Using the linear density field here—rather than the evolved nonlinear density field—should be a good approximation for the large scales of interest for our study; we focus on length scales of R ≳ 20 h⁻¹ Mpc and high redshift (z ≳ 6) in subsequent sections.

We mostly consider two concrete examples of upcoming redshifted 21 cm surveys. The first is based on the current, 128-tile version of the MWA (Tingay et al. 2013) and the second is based on an expanded, 500-tile version of the MWA (as described in Lonsdale et al. 2009, and considered in previous work such as McQuinn et al. 2006; Lidz et al. 2008). These two examples are intended to indicate the general prospects for imaging and bubble identification with first- and second-generation 21 cm surveys, respectively. Similar considerations would apply for other experiments, but we choose these as a concrete set of examples. We mainly focus on the 500-tile configuration in this paper because of its greater sensitivity. In Section 5.4, we shift to considering 128-tile configurations, and in Section 6, we consider a LOFAR-style interferometer for comparison. Hereafter, we refer to the 500-tile configuration as the MWA-500 and the 128-tile version as the MWA-128.

Throughout this paper, we work in comoving coordinates described by Cartesian labels (x–y–z), with Fourier counterparts (k_x–k_y–k_z). The Fourier modes can be connected directly with the u–v–ν coordinate system generally used to describe interferometric measurements. Here, u and v describe the physical separation of a pair of antennae in units of the observed wavelength, while ν describes the corresponding observed frequency. The instrument makes measurements for every frequency, ν, in its bandwidth, and for every antenna tile separation, (u, v), sampled by the array. In order to shift to a Fourier space description, the interferometric measurements must first be Fourier-transformed along the frequency direction. With our Fourier convention, the relation between the two sets of coordinates is given by

$$\begin{equation} k_{x} = \frac{2\pi u}{D} \hspace{28.45274pt} k_{y} = \frac{2\pi v}{D} \hspace{28.45274pt} k_z = \frac{2\pi }{\Delta \chi }, \end{equation} \tag{ 3 }$$

where D is the comoving distance to the survey center and Δχ is the comoving distance corresponding to a small difference in observed frequency of Δν (e.g., Liu et al. 2009). For small Δν/ν, we can express Δχ as

$$\begin{equation} \Delta \chi \approx \frac{c (1 + z_{\rm fid})}{H(z_{\rm{fid}})} \frac{|\Delta \nu |}{\nu }, \end{equation} \tag{ 4 }$$

where H(z_fid) is the Hubble parameter at the fiducial redshift and |Δν|/ν is the absolute value of the fractional difference between two nearby observed frequencies.

In order to test the prospects for imaging and bubble identification with the MWA, we must corrupt the underlying 21 cm signal described in Section 2.2 with thermal noise. We do this by generating a Gaussian random noise field in the k-space coordinate system described above, using an appropriate power spectrum. We assume that the co-variance matrix of the thermal noise power is diagonal in k-space. We add the resulting noise field to the underlying 21 cm signal (Equation (2)). The power spectrum of the thermal noise is given by (McQuinn et al. 2006; Furlanetto & Lidz 2007)

$$\begin{equation} P_{\rm{N}}(k,\mu) = \frac{T_{\rm{sys}}^2}{Bt_{\rm{int}}}\frac{D^2\Delta D}{n(k_\perp)}\left(\frac{\lambda ^2}{A_{\rm{e}}} \right)^{2}. \end{equation} \tag{ 5 }$$

Here, μ is the cosine of the angle between wavevector k = |k| and the line of sight, so that $k_{\perp } = \sqrt{\vphantom{A^A}\smash{\hbox{$1-\mu ^2$}}} k$ is the transverse component of the wavevector. We assume a system temperature of T_sky = 280[(1 + z)/7.5]^2.3 K (Wyithe & Morales 2007) and a total observing time of t_int = 1000 hr, which is an optimistic estimate for the observing time in one year. At our fiducial redshift of z_fid = 6.9, the comoving distance to the center of the survey is D = 6.42 × 10³ h⁻¹ Mpc. In this equation, λ denotes the observed wavelength of the redshifted 21 cm line, λ = 0.211(1 + z) m, and A_e is the effective area of each antenna tile. We determine A_e by linearly extrapolating or interpolating from the values given in Table 2 of Bowman et al. (2007); the effective area at z_fid = 6.9 is A_e = 11.25 m². We assume that the full survey bandwidth of 32 MHz is broken into individual blocks of bandwidth B = 6 MHz to protect against redshift evolution across the analysis bandwidth (McQuinn et al. 2006). The comoving survey depth corresponding to a B = 6 MHz chunk is ΔD = 69 h⁻¹ Mpc. The n(k_⊥) term describes the configuration of the antenna tiles. More specifically, it is the number density of baselines observing modes with transverse wavenumber k_⊥ (McQuinn et al. 2006). Following Bowman et al. (2006) and McQuinn et al. (2006), we assume that the antenna tiles are initially packed as closely as possible in a dense compact core, and that the number density of antenna tiles subsequently falls off as r⁻² out to a maximum baseline of 1.5 km. The radius of the dense antenna core is set by the requirement that the antenna density falls off as r⁻² outside of the core, and that it integrates to the total number of antennae. For the MWA-500, this gives r_c = 20 m, while for the MWA-128, the core radius is r_c ≈ 8 m. Equation (5) gives the noise power spectrum in units of mK², and so we divide by $T_0^2$ to combine with the simulated 21 cm signal expressed in units of T₀.

Note that the volume of the MWA survey differs somewhat from that of our reionization simulation. In particular, the transverse dimension of the simulation is smaller than that of the MWA by a factor of ∼3, while the simulation is deeper in the line-of-sight direction by about the same factor, as compared with the full MWA bandwidth. However, we remove the long-wavelength modes along the line-of-sight direction to mimic foreground cleaning (Section 2.4), and so we do not, in practice, use the longer line-of-sight scales in our simulation box. As we will see, the ionized regions in the simulation are substantially smaller than the transverse length of the box. Transverse slices should therefore be representative of what the actual MWA will observe from a fraction of its larger field of view. We have checked that the coarser transverse k-space sampling in the simulation compared to in the actual MWA survey does not impact our results.

Next, we need to consider contamination from foreground emission at the frequencies of interest. The relevant foregrounds include diffuse Galactic synchrotron radiation, extragalactic point sources, and Galactic bremsstrahlung radiation. While these foregrounds are many orders of magnitude brighter than the cosmological 21 cm signal, they are expected to individually follow smooth power laws in frequency. Over a sufficiently small frequency range, the summed contributions can also be approximated as following a smooth power law, while the 21 cm signal will vary rapidly. This allows the foregrounds to be removed from the data by, for example, fitting a low-order function along each line of sight and subtracting it. While this procedure is effective at removing foreground contamination, it also removes long-wavelength modes along the line of sight from the signal itself, and hence prevents measuring these modes. Several related methods for foreground removal have been discussed in the literature (e.g., Wang et al. 2006; Harker et al. 2009; Petrovic & Oh 2011; Chapman et al. 2012). In this work, we approximately mimic the degrading effects from foreground removal by subtracting the running mean from the noisy signal along each line of sight, rather than including realizations of the foregrounds in our simulation and excising them with one of the above algorithms. We generally remove the running mean over a bandwidth of 16 MHz, which corresponds to a comoving distance of L_fg = 185 h⁻¹ Mpc  at redshift z_fid = 6.9; we consider the impact of other choices of L_fg in Section 5.3. We defer more detailed models of foreground contamination, and foreground removal algorithms, to future work.

The Wiener filter is the optimal filter for extracting an input signal of known power spectrum when it is corrupted by additive noise, also with known power spectrum. As described in Press et al. (2002), this filter is optimal in that it minimizes the expectation value of the integrated squared error between the estimated signal field and the true signal field. The estimate of the true signal is a convolution of the Wiener filter and the corrupted signal in real space, and so is a product of the two quantities in Fourier space,

$$\begin{equation} \tilde{S}({\bf k}) = C({\bf k})W({\bf k}), \end{equation} \tag{ 6 }$$

where C(k), W(k), and $\tilde{S}({\bf k})$ are the Fourier transforms of the corrupted signal, Wiener filter, and estimated signal, respectively. Requiring that the filter be optimal in the least-square sense results in W(k) taking the form

$$\begin{equation} W(k,\mu) = \frac{P_{\rm{S}}(k)}{P_{\rm{S}}(k) + P_{\rm{N}}(k,\mu)}, \end{equation} \tag{ 7 }$$

where P_S(k) and P_N(k, μ) are the power spectra of the signal and noise, respectively. We note that, while the signal power spectrum is roughly isotropic,¹ the noise power spectrum depends on μ and consequently so does the filter. The filter keeps a unity weighting for k-modes, where P_S(k) ≫ P_N(k) and significantly downweights k-modes, where P_S(k) ≪ P_N(k). This can allow for partial recovery of the original signal, provided that the signal power dominates for some k-modes.

The Wiener filter requires an estimate of the signal power spectrum, P_S(k), and of the total (signal plus noise) power spectrum, P_S(k) + P_N(k, μ), as inputs. These may not be precisely known. However, since the filter is the outcome of a minimization problem—i.e., it minimizes the expected difference between the estimated and true fields—the accuracy of the filter should be insensitive to small changes about its optimal value. In other words, the accuracy of the filter is not expected to change greatly by using estimates of the signal and noise power spectra rather than the true spectra.

Furthermore, we do expect to have an estimate of the underlying signal power spectrum; measuring this statistic is a major goal of redshifted 21 cm surveys. Specifically, the underlying signal power can be estimated by cross-correlating redshifted 21 cm measurements made over two different time intervals (after foregrounds have been removed). The statistical properties of the signal should be identical across the two different time periods, but the thermal noise contributions will be independent. The cross-correlation between two time chunks then provides an unbiased estimate of the signal power (e.g., Liu et al. 2009). Estimates of the noise power spectrum can then be made by subtracting the estimated signal power from the power measured over the entire integration time, which contains both the signal and noise contributions. The Wiener filter does not actually require the noise power spectrum to be known on its own. However, in Section 4 we consider the optimal matched filter, which does have this requirement. Throughout this study, we assume perfect knowledge of the underlying power spectra.

Before applying the Wiener filter to our corrupted simulations, it is useful to estimate the expected signal-to-noise ratio of the filtered maps analytically, using simulated signal power spectra and the noise power spectrum of Equation (5). The expected signal-to-noise ratio of the Wiener-filtered field is $\mathcal{S}_{\rm wf} = \tilde{\sigma }_{\rm S}/\tilde{\sigma }_{N}$, where $\tilde{\sigma }^{2}_{{\rm S(N)}}$ is the filtered signal (noise) variance. The signal and noise variance can in turn be calculated as integrals over their respective power spectra,

$$\begin{equation} \tilde{\sigma }^2_{\rm{S(N)}} = \int \frac{\rm{\it d}^3 k}{(2\pi)^3}|W(k,\mu)|^{2}P_{\rm{S(N)}}(k,\mu). \end{equation} \tag{ 8 }$$

Here, we use P_S(N) to denote the power spectrum of the signal (noise). One can also consider the impact of foreground cleaning here by downweighting modes where the foreground power is large compared to the signal power. In order to consider the dependence of the signal-to-noise ratio on the stage of reionization, we consider simulation outputs in which the volume-averaged ionization fraction is 〈x_i〉 = 0.51, 0.68, 0.79, and 0.89. We consider each of these models at our fiducial redshift of z_fid = 6.9.² Currently, we do not consider still earlier stages of reionization since the prospects for imaging with the MWA-500 are especially poor for lower ionized fractions.

The resulting Wiener filters for the different values of 〈x_i〉 are shown in Figure 1, after integrating over angle μ. In this figure, foreground cleaning has been accounted for by subtracting a running mean along the line of sight, as described in Section 2.4. It is helpful to note, from Equation (7), that the filter is equal to 1/2 for modes where the signal and noise power are equal. The figure suggests that a small range of k-modes with k ≲ 0.1 h Mpc⁻¹ will have signal-to-noise ratio larger than unity for all four ionized fractions considered, although imaging is less promising for the smaller ionized fractions. If the ionized regions are larger than in our fiducial model—as expected if, for example, rarer yet more efficient and more clustered sources dominate reionization (e.g., McQuinn et al. 2007; Lidz et al. 2008)—then the prospects for imaging may improve somewhat. Performing the integrals in Equation (8), while incorporating foreground cleaning, we find that the total signal-to-noise ratio expected for the MWA-500 is $\mathcal{S}_{\rm{wf}} = 0.52$, 0.79, 1, and 1.2 for 〈x_i〉 = 0.51, 0.68, 0.79, 0.89, respectively.

Zoom In Zoom Out Reset image size

With the analytic signal-to-noise ratio estimates as a guide, we apply the Wiener filter to our mock noisy redshifted 21 cm data. The results of these calculations, for a particular slice through the simulation volume, are shown in Figure 2. The side length (1 h⁻¹ Gpc) of each slice is a factor of ∼3 smaller than the transverse dimension of the MWA. One can asses how well the original signal is "recovered" by comparing the top-left panel of the figure which shows the input signal with the bottom-left panel which shows the filtered noisy signal, after mimicking foreground removal. The two panels do not bear a striking resemblance since the average signal-to-noise ratio is only of order unity. Nonetheless, it is encouraging that many of the minima in the filtered noisy signal do indeed correspond to ionized regions in the input signal. Furthermore, we can compare the filtered noisy signal in the bottom-left panel with the top-right panel, which shows filtered pure noise. While these two panels do not look drastically different, they are easily distinguishable from each other given the increased contrast in the filtered noisy signal. Note that the noisy signal—in the absence of any filtering (not shown)—looks like pure noise since the signal to noise per resolution element is extremely low, and so filtering helps significantly. In addition, we see that the filtered noisy signal obtains signal-to-noise values exceeding 6σ, while the statistical significance of the filtered noise does not exceed ∼5σ. Quantitatively, ∼3% (∼0.03%) of the volume in the filtered noisy signal is occupied by pixels with statistical significance greater (in magnitude) than 3σ (5σ). This is expected given that the filtered data cube has an average signal-to-noise ratio of σ_S/σ_N ≈ 1, as anticipated in the analytic calculation of Section 3.1.

Zoom In Zoom Out Reset image size

Comparing the filtered noisy signal and the filtered pure noise, one can see that ionized regions in the underlying signal are diminished if they happen to be coincident with upward fluctuations in the noise, as expected. For example, the ionized region in the bottom-right corner of the unfiltered signal lies very close to a ∼3σ upward fluctuation in the filtered noise and, as a result, appears with weak statistical significance in the filtered noisy signal. Conversely, some of the most statistically significant regions in the filtered noisy signal occur when ionized regions overlap downward noise fluctuations. We can further compare the filtered noisy signal with the filtered noise-less signal, as shown in the bottom-right panel of Figure 2. The filtered noiseless signal is normalized by the standard deviation of the filtered noise so that it can be compared with the signal-to-noise slices in the other panels. This comparison reveals that high significance regions (≳ 5σ) in the filtered noisy signal only line up well with the corresponding regions in the filtered noiseless signal if they are coincident with downward fluctuations in the noise. On its own, the filtered noiseless signal only attains statistical significances of ≲ 4σ_N. Finally, Figure 3 illustrates the impact of foreground cleaning, performed here over a bandwidth of 16 MHz (Section 2.4). Foreground cleaning removes the long-wavelength modes along the line of sight—which is along the vertical axis in the figure—and thereby compresses structures along the line of sight. However, the cleaning process only impacts the long-wavelength line-of-sight modes which still leaves room to image other modes robustly.

Zoom In Zoom Out Reset image size

Note that the slice thickness (8 h⁻¹ Mpc) in Figures 2 and 3 is somewhat arbitrary. However, the Wiener filter smooths out structure on significantly larger scales than this (Figure 1), and so we expect similar results for other values of the slice thickness, provided the slice is thin compared to the cutoff scale of the filter. In practice, of course, one can make many independent maps similar to Figure 2 from the MWA-500 or similar surveys. Collectively, our results mostly confirm previous wisdom; the prospects for imaging with the MWA-500 are limited. Nonetheless, it appears that a signal-to-noise ratio of order unity is achievable, suggesting that the MWA-500 can make low-resolution maps of the reionization process.

The optimal matched filter is suited for the case of a corrupted signal containing a known feature that one would like to extract. The filter acts in Fourier space by cross-correlating the corrupted signal with a template describing the known feature, while downweighting k-modes in the corrupted signal by the noise power. The resulting form of the filter in Fourier space, M(k, μ), is

$$\begin{equation} M(k,\mu) = \frac{T({\bf k})}{P_{\rm{N}}(k,\mu)}, \end{equation} \tag{ 9 }$$

where T(k) is the Fourier profile of the known feature. The filter is optimal in the sense that it maximizes the signal-to-noise ratio in the filtered data cube at the location of the feature being extracted. While the Wiener filter requires an estimate of the signal and total (signal plus noise) power spectra, the matched filter requires a good estimate of the template profile, $T({ \mathchoice{{\hbox{k}}}{{\hbox{k}}}{{\hbox{k}}}{{\hbox{k}}} })$, and the noise power spectrum, P_N(k, μ). For our present application, we would like templates describing the ionized regions. An appropriate choice is not obvious; theoretical models predict that the ionization state of the gas during reionization has a complex, and somewhat uncertain, morphology, with ionized regions of a range of sizes and shapes (Iliev et al. 2006; Zahn et al. 2007; McQuinn et al. 2007). However, we find that the simplest conceivable choice of template filters, corresponding to completely ionized spherical bubbles of varying size, is nonetheless effective at identifying ionized regions with a more realistic and complex morphology. In this case, $T({ \mathchoice{{\hbox{k}}}{{\hbox{k}}}{{\hbox{k}}}{{\hbox{k}}} };R)$ is just the Fourier transform of a spherical top hat of radius R and is given by

$$\begin{equation} T({\bf k};R_{\rm{T}}) = \frac{V}{k^3R_{\rm{T}}^3}\left[ -kR_{\rm{T}}\cos kR_{\rm{T}}+ \sin kR_{\rm{T}}\right], \end{equation} \tag{ 10 }$$

with V denoting the volume of the spherical top hat. Note that the precise normalization of the filter is unimportant since we are mainly interested in the signal-to-noise ratio here, in which case the overall normalization divides out.

It is instructive to first consider an idealized test case that can be treated analytically before applying the matched filter to our full mock 21 cm data sets. In particular, we consider the case of an isolated, spherical, and highly ionized region placed at the origin and embedded in realistic noise. We assume that the neutral fraction exterior to the ionized region is uniform, with a mass-weighted neutral fraction of 〈x_H I(1 + δ_ρ)〉. Ignoring foreground contamination for the moment, the 21 cm signal may be written as

$$\begin{equation} \delta T_{\rm{b}}({\bf x}) - \langle \delta T_{\rm{b}} \rangle = \tilde{B}({\bf x};R_{\rm{B}}) + \tilde{N}({\bf x}), \end{equation} \tag{ 11 }$$

where B(x; R_B) denotes our isolated bubble of radius R_B and N(x) denotes the thermal noise contribution to the signal. We have subtracted off the overall mean brightness temperature, 〈δT_b〉, since this will not be measured in an interferometric observation. The tildes indicate that the spatial average has been removed from each of the underlying signal and noise so that $\tilde{B}({\bf x};R_{\rm{B}})$ and $\tilde{N}({\bf x})$ each have zero mean. In this case $\tilde{B}({\bf x};R_{\rm{B}})$ has an inverted spherical top-hat profile,

$$\begin{eqnarray} \tilde{B}({\bf x};R_{\rm B}) =\left\lbrace \begin{array}{@{}ll@{}} -\langle x_{{\rm H\scriptsize{I}}}(1+\delta _{\rho }) \rangle & {|{\bf x}| < R_{\rm B},} \\ 0 &{\rm otherwise.} \end{array}\right. \end{eqnarray} \tag{ 12 }$$

The Fourier transform of the isolated bubble is hence related to the Fourier transform of our template by $\tilde{B}({\bf k};R_{\rm B}) = -\langle x_{{{\rm H}\scriptsize{I}}}(1+\delta _\rho) \rangle T({\bf k};R_{\rm B})$. Note that we express brightness temperatures in units of T₀ (see Equation (2)), and so all quantities here are dimensionless.

It is straightforward to derive the expected signal-to-noise ratio at the center of the isolated ionized region, and thereby gauge the prospects for bubble detection with a matched filter technique. Let us assume that the radius, R_B, of our template filter is well matched to the true radius of the ionized region. This will maximize the expected signal-to-noise ratio. Neglecting foregrounds for the moment, and using the fact that the thermal noise has zero mean, we find that the signal-to-noise ratio at bubble center for the optimal matched filter is

$$\begin{equation} {\mathcal{S}}(R_{\rm B}) = - \langle x_{{\rm H}\scriptsize{I}}(1 + \delta _\rho) \rangle \left[\int \frac{\rm{\it d}^3 k} {(2\pi)^3}\frac{T^2(k;R_{\rm B})} {P_{\rm N}(k,\mu)}\right]^{1/2}. \end{equation} \tag{ 13 }$$

For our sign convention, in which the template and ionized regions have opposite signs, this quantity is negative—ionized bubbles are regions of low 21 cm signal. The contribution of a Fourier mode to the signal-to-noise ratio depends on the relative size of T²(k; R_B) and P_N(k, μ): modes for which the template is much larger than the noise power contribute appreciably to ${\mathcal {S}}(R_{\rm B})$ while modes dominated by the noise power are not useful. The signal-to-noise ratio depends on the neutral fraction: a larger exterior neutral fraction increases the contrast between an ionized bubble and the exterior, and hence boosts the detectability of the ionized region. We would like to calculate the expected signal-to-noise ratio for ionized regions of different sizes and for various volume-averaged ionization fractions. To do this, we need to connect the volume-averaged ionized fraction with the mass-averaged fraction, 〈x_H I(1 + δ_ρ)〉, which enters into Equation (13). Here, we should incorporate that large-scale overdense regions are generally ionized before typical regions during reionization, i.e., the neutral fraction and overdensity fields are anti-correlated. Defining δ_x = (x_H I − 〈x_H I〉)/〈x_H I〉, we approximate 〈δ_xδ_ρ〉 as fixed at 〈δ_xδ_ρ〉 = −0.25 throughout the reionization process (Lidz et al. 2007).

The results of the signal-to-noise calculation are shown in Figure 4 for the MWA-500 and a LOFAR-type experiment. Here we consider only our fiducial redshift, z_fid = 6.9. The (absolute value of) $\mathcal {S}(R_{\rm{B}})$ is evidently a strongly increasing function of bubble size. This occurs because the thermal noise is a strong function of scale and only the rather large-scale modes are measurable. It is encouraging that the expected signal-to-noise ratio exceeds five, $\mathcal {S}(R_{\rm{B}}) \gtrsim 5$, for a range of radii and neutral fractions. This corresponds to a 5σ detection: "false" bubbles at this significance from downward fluctuations in the noise are highly unlikely, with a fraction of only ∼3 × 10⁻⁷ of pixels in the filtered noise having such a large (negative) significance on their own. For simplicity we neglect the impact of foreground cleaning in this figure: this will degrade the expected signal-to-noise ratios somewhat, as we will consider subsequently (see Sections 4.3 and 5.3).

Zoom In Zoom Out Reset image size

In order to estimate the number of bubbles that can be detected from these curves, we need to consider how many bubbles there are of different sizes, i.e., we need to fold in an estimate of the bubble size distribution. In particular, while the contrast of an ionized region increases with the neutral fraction, large ionized bubbles become increasingly scarce for larger values of the neutral fraction. For instance, we can consider the model bubble size distributions in Figure 4 of Zahn et al. (2007). This figure indicates that bubbles of radius larger than 30 h⁻¹ Mpc are exceedingly rare for neutral fractions larger than 〈x_H I〉>0.5, with only the tail end of the distribution extending past 25 h⁻¹ Mpc. However, bubbles of this size are relatively common later in reionization. Since Figure 4 indicates that only bubbles with R ≳ 30 h⁻¹ Mpc exceed $\mathcal {S}(R_{\rm{B}}) \gtrsim 5$, this suggests that bubble detection is feasible for the MWA-500 after the universe is more than ∼50% ionized, but that it will be difficult to use this method at earlier stages of the reionization process. Also, note again that the calculation here neglects the effects of foreground cleaning. However, we find that incorporating foreground cleaning only has a small effect on bubbles of this size (≲ 30 h⁻¹ Mpc, Section 5.3). Bubble detection will also be challenging once the universe is less than ∼10%–20% neutral, owing mostly to the reduced contrast between the bubbles and typical regions. If the ionized bubbles at a given stage of the EoR are larger than in the model of Zahn et al. (2007), then the prospects for bubble detection will be enhanced. We refer the reader to McQuinn et al. (2007) for a quantitative exploration of the bubble size distribution across plausible models for the ionizing sources.

Finally, it is interesting to consider a LOFAR-style interferometer, as discussed further in Section 6. This is shown as the red dotted curve in Figure 4. The expected $\mathcal {S}(R)$ exceeds that of the MWA-500 for small bubble radii, before flattening off at larger radii. This occurs because the LOFAR-style interferometer has more collecting area per baseline, but a larger minimum baseline. This makes it more sensitive to the smaller ionized regions, but less sensitive to larger ones.

While the signal-to-noise curves in this toy case provide a useful guide, we should keep in mind their limitations. First, it considers only the case of a single isolated ionized region. Next, we consider here only the signal to noise at the bubble center, while an ionized region will typically have a strong (negative) signal to noise over much of its volume. Specifically, for bubbles with 10 Mpc h⁻¹ ≲ R_B ≲ 60 Mpc h⁻¹, we find that the signal-to-noise value at the edge of bubble is roughly half that at its center. This can help significantly with detection. Finally, we consider only the average signal-to-noise ratio here. In practice, the signal-to-noise ratio in a filtered map may fluctuate significantly around this average, as we will see.

With the estimates of the previous section as a rough guide, we now apply the matched filter to our noisy mock redshifted 21 cm data. In order to illustrate the results of passing our mock data through a matched filter, we start by examining simulated signal-to-noise fields for a single template radius of R_T = 35 h⁻¹ Mpc. This template radius corresponds to the typical size of the ionized bubbles we believe we can detect (see Figure 4). A representative slice through the simulation is shown in Figure 5. The results look promising, with signal-to-noise ratios comparable to the values anticipated in the idealized calculation of Figure 4. Although the Wiener filter provides the best overall map, or data cube, one can still detect individual features at greater significance by applying a matched filter. Comparing with Figure 2, it is clear that the Wiener filter is passing more small-scale structure than the matched filter shown here. This results in the signal-to-noise ratio being larger (in absolute value) for the matched filter than for the Wiener filter. In particular, we find values of the signal-to-noise ratio that are as low as ∼ − 10 in the matched filter data cube, a significant improvement over the global minimum of ∼ − 6 for the Wiener filter. Moreover, we can compare the filtered noisy signal in the bottom-left panel with the filtered pure noise field in the top-right panel. They differ by more than in the case of the Wiener filter. Indeed, the very low signal-to-noise ratio regions (shown in dark blue/purple in the bottom-left panel) line up fairly well with ionized regions in the top-left panel. This is especially apparent when comparing the filtered noisy signal to the filtered noiseless signal, as shown in the bottom-right panel. For the slice shown, almost all of the significant features in the filtered noiseless signal are preserved in the noisy case. Figure 6 shows the impact of foreground cleaning: as in the case of the Wiener filter (Figure 3), this compresses structures along the line of sight and reduces the overall signal-to-noise ratio in the data cube. The signal-to-noise ratio is still significant enough, however, to robustly identify ionized regions.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

As with the Wiener filter, and in what follows subsequently, we show slices of 8 h⁻¹ Mpc thickness. This choice is arbitrary, but we expect similar results provided the slice thickness is small compared to the radius of the template filter. It is important to keep in mind, however, that the full data cube will consist of many separate slices of this thickness. Also note that the transverse dimension of the MWA-500 is larger than that of our simulation box by a factor of ∼3, and so these slices represent only ∼1/9 of the MWA field of view.

These results are promising, but they are for a single filtering scale, and so we can do significantly better by considering a range of template radii, and looking for extrema in the resulting signal-to-noise fields. In particular, we proceed to apply a sequence of filters with template radii up to R_T ⩽ 75 h⁻¹ Mpc —see Section 4.5 for a justification of this maximum—across the simulation volume. We assign the minimum (most negative) signal-to-noise value obtained over the range of template radii to each simulation pixel and use this to construct a new field. The position of local minima in this field is chosen to be the centers of candidate bubbles, and each such bubble is assigned a radius according to the scale of the template filter that minimizes its signal to noise. We focus on minimum values since ionized regions are expected to appear as regions of low 21 cm signal. All candidate bubbles whose central signal to noise is lower than −5 are considered to be detected ionized regions.

We find it important to apply one additional criterion to robustly identify ionized regions. The criterion is that a low signal-to-noise region on scale R_B must additionally be low in signal to noise at all smaller smoothing scales, R_T < R_B. This guards against the possibility that a detected bubble will be centered on neutral material that is nevertheless surrounded by ionized hydrogen. A region like this will have a high (least negative) signal to noise when filtered on small scales and then dive down (gaining statistical significance) when filtered on scales containing the surrounding ionized material. We discard such spurious bubbles by requiring that the field is low on all smaller smoothing scales. The only downside to this procedure is that it occasionally discards true ionized regions whose center happens to coincide with a significant upward noise fluctuation. Overall, however, it improves the quality of detected bubbles (Section 4.4). This cut also requires a threshold choice; we reject candidate regions if their signal-to-noise ratio crosses above a threshold $\mathcal {S}_{\max }$ at any smoothing scale less than R_T. After trying several thresholds, we found the most effective choice to be $\mathcal {S}_{\max } = -1\sigma$. In principle, one might use the full curve of signal-to-noise ratio versus template radius for each candidate bubble to help verify the detection and determine the properties of the bubble. In practice, we found that individual signal-to-noise curves are noisy and difficult to incorporate into our analysis and so we do not consider this possibility further in what follows.

We apply this algorithm to the mock redshifted 21 cm data and identify 220 ionized regions across the simulation volume (which is different than the MWA survey volume, as we will discuss subsequently). A representative example of a detected bubble is shown in Figure 7. The circle in the figure identifies the detected bubble size and the location of its center in both the filtered noisy signal (top-left and top-right panels), as well as in the input signal (bottom-left and bottom-right panels). The algorithm has convincingly identified an ionized region. The detected bubble overlaps a small fraction (10%) of neutral material in the input signal. Although this particular ionized bubble is well identified, most of the ionized regions in the signal will escape detection. This is because the significance levels of the detected bubbles are not that high, and an ionized region generally needs to be coincident with a downward fluctuation in the noise to pass our significance threshold. For example, consider the larger ionized region below and to the left of the detected region in the bottom-left panel of Figure 7. This region, while larger and therefore more detectable on average than the identified bubble, happens to coincide with a large upward noise fluctuation and hence fails to cross the significance threshold. While we cannot identify all of the large ionized regions in the noisy mock data, we can robustly identify some regions; this may still be quite valuable.

Zoom In Zoom Out Reset image size

It is also clear that the underlying ionized regions are manifestly non-spherical, creating an ambiguity as to what the appropriate "radius" of the region is. Focusing on the bottom-right panel in Figure 7, we could imagine the size being reasonably described by a radius ∼50% larger, so as to enclose more of the nearby ionized material. However, our method naturally favors radii causing little overlap with neutral material at these size scales. Therefore, an ionized region like the one shown in Figure 7 is more likely to be detected as several small ionized regions than one large one, although both characterizations seem reasonable.

Figure 8 gives a further example of how the algorithm identifies bubbles, and some of the ambiguities that can result. This figure shows an example of an irregular, yet contiguous, ionized region that is detected as more than one ionized bubble. Here, we show spatial slices through the center of the middle sphere, marked with a solid circle, which happens to intersect neighboring ionized bubbles, whose cross sections are shown as dashed circles. Hence, our algorithm generally represents large, irregularly shaped, yet contiguous, regions as multiple ionized bubbles.

Zoom In Zoom Out Reset image size

It is important to emphasize further the difference between the simulated results shown here and the idealized test case of the isolated bubble shown in the previous section. In particular, we consider here the application of matched filters to the 21 cm signal during a late phase of reionization in which many ionized regions, with a broad size distribution, fill the survey volume: the ionized regions are not isolated bubbles in a sea of partly neutral material. When applying a matched filter of template radius R_T around a point, the values of the field at many neighboring pixels impact the filtered field at the point in question. It is hence possible that a filtered pixel is affected by several distinct neighboring ionized regions. Indeed, this can result in even neutral regions having low signal-to-noise ratios provided they are surrounded by many nearby ionized regions. For instance, in the low noise limit, any region with volume-averaged neutral fraction lower than the cosmic mean would pass our significance threshold. To guard against this type of false detection, we implemented the requirement that a candidate bubble has low signal to noise for all template radii smaller than the detected radius. Another possibility might be to treat small ionized regions as an additional noise term in the filter. However, in practice, our attempts along these lines introduced an additional level of model dependence without significantly increasing the quality of the detected bubbles. Ultimately, it is important to keep in mind that the signal-to-noise values quoted here reflect only the likelihood that a value arises purely from noise, and so they are not strictly indicative of the quality of the detected bubbles.

We hence turn to describe the characteristics of the detected ionized regions, and to quantify the method's level of success in detecting ionized bubbles. To do this, we calculate the fractional overlap of each detected bubble with ionized material in the underlying signal. Additionally, we estimate how many ionized bubbles should be detectable across the entire MWA-500 survey volume.

The matched filter technique finds 220 bubbles across our simulation volume. However, the algorithm for determining bubble positions and sizes allows for bubbles to occupy overlapping areas, as shown in Figure 8. We find that ∼55% of the detected bubbles have some overlap with another bubble, although only ∼15% of the total volume occupied by detected bubbles is occupied by more than one. Regardless, ∼96% of the detected ionized bubbles have an average ionized fraction larger than x_i = 0.79, which is the volume-averaged ionization fraction of the simulation box at the redshift under consideration. Furthermore, ∼42% of the detected bubbles have an ionized fraction greater than x_i = 0.9. The lowest ionization fraction of a detected bubble is x_{i, min} = 0.77, just slightly below the volume-averaged ionization fraction of the simulation. In total, we detect nine bubbles whose ionized fractions are less than the average ionization fraction of the box. Inspection reveals that these regions happen to be coincident with significant (⩽−3σ) downward noise fluctuations.

In Figure 9, we plot the volume-averaged ionized fraction within each of our detected bubbles against the detected bubble radius for a significance threshold of 5σ. For comparison, we show the 1σ spread in the ionized fraction enclosed by randomly distributed spheres of the same size.³ The spread in ionization of the randomly distributed spheres around the box average ionized fraction, 〈x_i〉 = 0.79, decreases with increasing radius; this reflects the drop-off in the power spectrum of the ionization field toward large scales. Most of the detected bubbles are significantly more ionized than random regions, as expected, indicating a significant success level. There are a few poor detections which result mostly from downward noise fluctuations. There is a small overall decrease in the ionized fraction of detected bubbles larger than R_B ≳ 40 h⁻¹ Mpc, suggesting that we may no longer be detecting individual ionized regions here. These regions may potentially be distinguished from isolated bubbles by examining the signal-to-noise ratio as a function of template radius closely, as we discuss in Section 4.5.

Zoom In Zoom Out Reset image size

We can estimate the number of ionized regions detectable in the MWA-500 by scaling from our simulation volume to the MWA survey volume. At z_fid = 6.9, for an ionized fraction of 〈x_i〉 = 0.79, we expect to find 140 bubbles in a B = 6 MHz chunk of the MWA, over its entire field of view of ∼770 deg². About 135 (60) of these detected bubbles are expected to have ionized fractions larger than 79% (90%). This estimate comes from simply scaling our simulation volume (which is deeper than the MWA bandwidth) to a 6 MHz portion of the MWA survey volume. Analyzing the MWA data over a 6 MHz chunk is meant to guard against redshift evolution: the full bandwidth of the survey is B = 32 MHz and so the prospects for bubble detection across the full survey are even better than this estimate suggests. The precise gain will be dependent on how rapidly the bubble size distribution evolves across the full survey bandwidth. One caveat with our estimate, however, is that B = 6 MHz corresponds to only ∼70 h⁻¹ Mpc . This is comparable to the size of our larger bubbles, and so analyzing chunks this small might weaken our ability to detect large bubbles. This effect is not incorporated in our scaling estimate, which simply takes the ratio of the MWA survey volume and our simulation volume. In practice, one can perform the bubble extraction for different analysis bandwidths to help ensure robust detections.

By increasing the significance threshold for detection, one might hope to improve the quality of the detected bubbles. However, we find that this leads to the detected bubbles being larger on average rather than being more ionized. This results from the fact that, at a fixed ionized fraction, a larger bubble will have a greater (in magnitude) signal-to-noise value than a smaller bubble. By increasing the significance threshold, some smaller, highly ionized bubbles are no longer detected, while larger, less-ionized bubbles are. Specifically, we find that by increasing the threshold to 6σ (7σ) the number of detected bubbles drops to ∼100 (∼50) with 99% (100%) being more ionized than the box, but only 34% (33%) being more than 90% ionized. Hence, for any of these thresholds we robustly detect many bubbles, although the detailed success rate and overall number vary somewhat with the precise choice of threshold. Note that lowering the significance threshold has the undesirable effect of increasing the number of false positives. Ultimately, with real data one should explore a range of threshold choices and examine the impact.

It is worth mentioning one further detail of our algorithm. In the previous section, we set the maximum template radius considered at R_{T, max} = 75 h⁻¹ Mpc, without justification. In fact, we have a sensible and automated way for arriving at this choice. We discuss this procedure briefly here.

A good candidate ionized region should in fact obey three criteria. First, it should have a large (negative) signal-to-noise ratio, so that it is unlikely to result from a noise fluctuation. Second, the signal-to-noise ratio should be small for all trial radii smaller than the optimal template radius, as discussed in Section 4.3. Finally, the total signal must itself be small in an absolute sense. In the limit of low noise, anything less neutral than average would qualify as a bubble by the first criterion, and so this third criterion may then become important for robustly identifying bubbles. This low noise limit is relevant for the MWA-500 only on very large smoothing scales, where the noise averages down significantly.

Since this third criterion becomes important only on very large smoothing scales here, we use it only to set the maximum template radius considered. Without this third consideration, our algorithm generally identifies a few excessively large ionized bubbles, but this can be easily understood and avoided as follows. Consider, for the moment, the 21 cm brightness temperature field in the absence of noise and foregrounds. Let us further work in units of T₀ (Equation (2)), and remove the average brightness temperature contrast across the data cube. In this case, the signal inside a highly ionized bubble is expected to be −〈x_H I(1 + δ_ρ)〉. If we now spherically average the field on scales smaller than the bubble, the value at bubble center will not change from this value, −〈x_H I(1 + δ_ρ)〉. Once the smoothing scale becomes larger than the bubble scale, however, surrounding neutral material will increase the value of the filtered field at bubble center. Hence, if the filtered field becomes everywhere larger than −〈x_H I(1 + δ_ρ)〉 on some smoothing scale, it is clear that no larger ionized bubbles exist within the data cube. This suggests that we can set the maximum template radius by requiring that the filtered noisy signal reaches sufficiently small values, at some locations across the data cube, for there to still plausibly be completely ionized regions. Since the presence of noise only increases the variance, this should provide a conservative estimate of the maximum size of the ionized regions. In practice, we need to choose a threshold criterion without assuming prior knowledge of the neutral fraction. Here, we set the maximum template radius to be the smallest smoothing scale at which the filtered noisy field everywhere exceeds −〈x_H I(1 + δ_ρ)〉 ⩾ −0.075. This corresponds to the expected contrast at 〈x_HI〉 = 0.1, assuming 〈δ_xδ_ρ〉 = −0.25, and yields a maximum template radius of R_{T, max} = 70 h⁻¹ Mpc , 73 h⁻¹ Mpc , 75 h⁻¹ Mpc , and 75 h⁻¹ Mpc  for 〈x_i〉 = 0.51, 0.68, 0.79, and 0.89, respectively. The precise threshold used here, −0.075, is somewhat arbitrary but this choice is only being used to set the maximum template radius considered.⁴

In order to consider bubble finding at earlier and later stages of the EoR, we run our matched filter on simulation outputs with volume-averaged ionized fractions of 〈x_i〉 = 0.51, 0.68, and 0.89. As discussed in Section 3.1 these outputs are actually at slightly different redshifts, but we generate the 21 cm field as though they were at z_fid = 6.9. As far as bubble detection is concerned, varying the ionized fraction leads to two, generally competing, effects. First, the bubbles grow as reionization proceeds. This tends to boost detection, since it is only the large-scale modes that are detectable over the thermal noise. Second, however, the contrast between an ionized region and a typical volume of the universe is reduced as reionization proceeds. This makes bubble detection more difficult. Both of these effects are quantified in the idealized isolated bubble case in Figure 4. It is also clear that the ideal ionized fraction for bubble detection will be somewhat survey dependent. As already illustrated in Figure 4 and discussed further in Section 6, a LOFAR-type interferometer will perform better when the ionized regions are smaller.

We find that the matched filter is capable of detecting ionized regions for each of the ionized fractions studied. In Figure 10, we show histograms of the detected bubble size distributions for each ionized fraction. Since we preferentially detect large ionized regions, we do not expect these distributions to be representative of the true underlying bubble size distributions. For example, in Figure 10, the size distribution peaks around ≳ 40 h⁻¹ Mpc  for the case of x_i = 0.79, despite volume-weighted size distribution peaking around ∼30 h⁻¹ Mpc  in Figure 4 of Zahn et al. (2007) at roughly the same ionized fraction. Nonetheless, the histograms illustrate a general shift from smaller to larger detected bubble radius as the ionized fraction increases. By applying the matched filter to several redshift bins, one can potentially observe precisely this evolution with the MWA-500. This would complement studies of the 21 cm power spectrum evolution over the same redshift range (e.g., Lidz et al. 2008). From the histograms, one can see that—of the models shown—the best ionized fraction for bubble detection is 〈x_i〉 = 0.79. This is apparently near the sweet spot for the MWA-500 where the bubbles are large enough in the model for detection, but the contrast with typical regions is still sufficiently large.

Zoom In Zoom Out Reset image size

The average ionized fraction within detected bubbles varies significantly across the different ionized fractions considered. Specifically, the percentage of detected bubbles that are more than 90% ionized is 0%, 15%, 43%, and 91% for 〈x_i〉 = 0.51, 0.68, 0.79, and 0.89, respectively. However, in each case the percentage of detected bubbles with ionized fraction larger than the (global) volume-averaged ionization fraction is fixed at ≳ 95%. At first glance, one aspect of these results may appear to be in tension with the calculations of Section 4.2, where we estimated that bubble detection would be unsuccessful for neutral fractions larger than 〈x_H I〉 ≳ 0.5. However, this estimate considered the detection of isolated bubbles. Inspection reveals that the detected bubbles at 〈x_i〉 = 0.51 each correspond to clusters of smaller ionized regions. Evidently, these appear as a single larger ionized region after convolving with the template filter and downweighting the noisy short-wavelength modes. In practice, it may be possible to distinguish this case from that of an isolated bubble by analyzing the trend of signal-to-noise ratio versus trial template radius. The signal-to-noise ratio is expected to grow more rapidly with radius (before reaching the bubble scale) for a truly isolated bubble.

We now consider how the prospects for bubble detection diminish if reionization occurs earlier and the observations are focused on the corresponding redshifts. In particular, we examine the case that our model with an ionized fraction of 〈x_i〉 = 0.79 is observed at a higher redshift. We focus on this case since this ionized fraction appears close to optimal for bubble detection. Aiming for only a rough estimate here, we consider the prospects for detecting an R_B = 40 h⁻¹ Mpc bubble.

Although several different factors in the noise power spectrum of Equation (5) scale with redshift, the dominant scaling is with the sky temperature. The noise power scales as $P_{\rm{N}} \propto T_{\rm{sky}}^{2}$, and the sky temperature follows T_sky∝ν^−2.6∝(1 + z)^2.6. Therefore, we expect the signal-to-noise ratio of a detected bubble to fall off with increasing observation redshift roughly as

$$\begin{equation} \mathcal {S}(z) = \mathcal {S}(z_{\rm{fid}}) \left(\frac{1+z_{\rm{fid}}}{1+z}\right)^{2.6}. \end{equation} \tag{ 14 }$$

This indicates the signal-to-noise ratio for a bubble detected with a signal to noise of $\mathcal {S}(z_{\rm{fid}})$ at redshift z_fid = 6.9, if the bubble were instead observed at redshift z. A relatively large bubble with R_B ≈ 40 h⁻¹ Mpc has a typical signal-to-noise ratio at bubble center of $\mathcal {S} \approx 4$ at our fiducial redshift. This value is found by incorporating foreground cleaning into the corresponding curve in Figure 4. According to Equation (14), the signal-to-noise value will be reduced to a significance of $\mathcal {S} \approx 2$ (1) at z = 9.3 (12.5). The bubble will, in fact, be more detectable than implied by this one number—the signal-to-noise ratio at bubble center—since an ionized region should have low signal to noise over much of its volume. From this, we conclude that bubble detection should be feasible with the MWA-500 if our fiducial ionized fraction occurs later than z ≲ 9 or so, but that the prospects are rather limited in the case of significantly earlier reionization. A range of recent work in the literature, however, suggests that reionization is unlikely to complete so early. See, for example, Figure 9 from the recent study of Kuhlen & Faucher-Giguère (2012) which combines Lyα forest data (Fan et al. 2006), measurements of the Thomson optical depth from WMAP (Komatsu et al. 2011), and measurements of the Lyman break galaxy luminosity function (Bouwens et al. 2011). Hence, the prospects for bubble detection appear good for the MWA-500.

Next, we consider the impact of variations around our standard foreground cleaning model. As discussed previously (Section 2.4), our standard assumption is that the impact of foreground cleaning can be approximately mimicked by removing the running mean, over a bandwidth of B = 16 MHz, across each line of sight. The optimal foreground cleaning strategy avoids "overfitting" by removing the smoothest possible function over the largest possible bandwidth, in order to preserve the underlying signal as much as possible. It also avoids "underfitting" by ensuring that foreground residuals do not excessively contaminate the signal. Liu & Tegmark (2012), for example, find that 21 cm foregrounds can be removed to one part in 10⁵ or 10⁶ by subtracting roughly three modes over ∼32 MHz of bandwidth. This should have a fairly similar impact to our fiducial cleaning model, but we would expect a bit more degradation in this case. A detailed investigation would add foreground contamination into our mock data cubes, and explore the impact of various cleaning algorithms directly.

Here, we instead check how our results change for slightly more and less aggressive foreground cleaning. In particular, we remove the running mean over each of B = 32 MHz, 8 MHz, and 6 MHz and rerun our bubble finding algorithm (at z_fid = 6.9, 〈x_i〉 = 0.79). This has little impact on our results for the case of B = 32 MHz and 8 MHz. In particular, the number of identified bubbles varies by less than 10%–15% and the quality of detected bubbles decreases slightly for the more aggressive cleaning model, and improves slightly in the most optimistic case. Specifically, for B = 16 MHz, 96% (43%) of bubbles have ionized fractions exceeding x_i = 0.79 (0.9); for B = 8 MHz the corresponding numbers are 90% (32%); and for B = 32 MHz the same numbers are 96% (62%). For the most aggressive case of B = 6 MHz, we find that the number of detected bubbles drops significantly to ∼50% of the number detected in the fiducial case. The larger degradation in the B = 6 MHz results because the foreground cleaning scale has become comparable to the scale of many of the detected bubbles. Despite this drop, we find the quality of the detected bubbles to be about the same as with the B = 8 MHz model, with 92% (36%) exceeding x_i = 0.79 (0.9). The distribution of detected bubble sizes also depends on the bandwidth used for foreground cleaning, with smaller bandwidths generally corresponding to smaller detected bubble sizes when B ≲ 16 MHz. This behavior is consistent with work by Petrovic & Oh (2011) who found that a bubble whose size is comparable to the scale for foreground cleaning should have its contrast reduced. Even in the most aggressive case considered, however, we still detect many ionized regions robustly. While these estimates are encouraging, a more detailed study is warranted. It may also be advantageous to estimate the power spectrum of the foregrounds, and incorporate this as an additional noise term in Equations (7) and (9) for each of the Wiener filter and the matched filter, respectively.

So far our analysis has focused on the MWA-500, which is meant to represent a second-generation 21 cm survey. In the near term, it is timely to consider the prospects for a 128-tile version of the MWA (the MWA-128) which is ramping up to take data in the very near future. This should be significantly less sensitive, since the number of baselines scales as the number of antenna tiles squared.

In order to generate thermal noise representative of the MWA-128, we start by considering a similar antenna distribution as for the MWA-500. In particular, we assume that all of the antenna tiles are packed as closely as possible within a core of radius 8 m and that the antenna distribution subsequently falls off as r⁻² out to a maximum baseline of 1.5 km. After comparing the thermal noise power spectrum in this configuration with that in Beardsley et al. (2012), we find that our noise power is larger by up to a factor of a few. This could possibly be due to our approximation of a smooth antenna distribution being less valid for the MWA-128, or to the fact that our analytic formula for the noise power spectrum does not incorporate a full treatment of rotation synthesis. In an effort to bracket somewhat the impact of the detailed antenna distribution, we further consider a configuration where all antenna tiles are packed as closely as possible in a dense core of radius ∼25 m. This resembles the "supercore" configuration considered in Lidz et al. (2008) for the power spectrum.

The results of applying the optimal matched filter for a single template radius of 35 h⁻¹ Mpc are shown in Figure 11 for the r⁻² tile distribution. This shows that the sensitivity is much lower than for the MWA-500, as expected. It is much more difficult to distinguish the filtered noisy signal (bottom-left panel) from the filtered pure noise (top-right) panel here than in Figure 5. Most of the significant, dark blue regions in the filtered noisy signal correspond simply to low noise regions. However, applying the detection algorithm we do nonetheless detect seven bubbles across a volume equivalent to a 6 MHz chunk of the MWA-128 survey. The success is generally lower than in the case of the MWA-500: here 75% of detected bubbles exceed the average ionized fraction of the box, while ∼42% exceed x_i = 0.9. In the supercore configuration, we find slightly higher significance levels (up to 6.9σ) but the identified regions generally correspond to several large clustered ionized regions, rather than a single ionized bubble. Altogether, the algorithm identifies 10 ionized regions across the MWA survey in the supercore configuration, but the identified regions have a lower overall quality than in the r⁻² configuration.

Zoom In Zoom Out Reset image size

Our conclusion is that bubble detection is only marginally possible with the MWA-128. While the results are unlikely to be very compelling, it is worth applying the matched filter to the first-generation surveys as an initial test. Even a few weakly identified candidate bubbles would provide compelling targets for follow-up observations. Another possibility is to focus on targeted searches around known bright sources for the MWA-128 (e.g., Wyithe & Loeb 2004; Datta et al. 2012).