Recovering the H ii region size statistics from 21-cm tomography

Abstract

We introduce a novel technique, called ‘granulometry’, to characterize and recover the mean size and the size distribution of H ii regions from 21-cm tomography. The technique is easy to implement, but places the previously not very well-defined concept of morphology on a firm mathematical foundation. The size distribution of the cold spots in 21-cm tomography can be used as a direct tracer of the underlying probability distribution of H ii region sizes. We explore the capability of the method using large-scale reionization simulations and mock observational data cubes while considering capabilities of Square Kilometre Array 1 (SKA1) low and a future extension to SKA2. We show that the technique allows the recovery of the H ii region size distribution with a moderate signal-to-noise ratio from wide-field imaging (SNR ≲ 3), for which the statistical uncertainty is sample variance dominated. We address the observational requirements on the angular resolution, the field of view, and the thermal noise limit for a successful measurement. To achieve a full scientific return from 21-cm tomography and to exploit a synergy with 21-cm power spectra, we suggest an observing strategy using wide-field imaging (several tens of square degrees) by an interferometric mosaicking/multibeam observation with additional intermediate baselines ( ∼ 2-4 km) in an SKA phase 2.

1 INTRODUCTION

The 21-cm signal arising from the spin flip transition in the neutral hydrogen (H i) atoms promises to unearth an unprecedented amount of information about the cosmic dawn (CD) and the epoch of reionization (EoR). It will help us to resolve many outstanding questions regarding this particular phase of the history of our Universe such as: how did reionization progress with time? What was the morphology of H i distribution in the intergalactic medium (IGM) during the different phases of reionization? Which were the major sources that drive reionization? (for reviews see, e.g. Furlanetto, Oh & Briggs 2006; Morales & Wyithe 2010; Pritchard & Loeb 2012).

The first generation of radio interferometers such as the Giant Metrewave Radio Telescope (GMRT)¹ (Paciga et al. 2013), Low-Frequency Array (LOFAR)² (Yatawatta et al. 2013), Murchison Widefield Array (MWA)³ (Bowman et al. 2013), and Precision Array for Probing the Epoch of Reionization (PAPER)⁴ (Ali et al. 2015) target the detection of this signal through statistical estimators such as the power spectrum and variance. While these statistical estimators ensure an optimum signal-to-noise ratio (SNR) to increase the detection possibility, it will not be possible to characterize the signal completely through these measurements alone, as the signal from the EoR is expected to be highly non-Gaussian in nature (e.g. Ciardi & Madau 2003; Furlanetto, Zaldarriaga & Hernquist 2004b; Bharadwaj & Pandey 2005; Cooray 2005; Iliev et al. 2006; Mellema et al. 2006b; Harker et al. 2009; Mondal et al. 2015; Dixon et al. 2016). As a result, many completely different underlying H i distributions could give rise to the same power spectrum and variance for the signal, which would lead to a large degeneracy and ambiguity in the estimated astrophysical parameters through these statistics. This is why it is essential to quantify and characterize the signal beyond these one- and two-point statistics.

One way to achieve this is through a hierarchy of the higher order correlation functions (such as three- and four-point correlation functions) and their corresponding Fourier analogues (Bharadwaj & Pandey 2005; Pillepich, Porciani & Matarrese 2007; Cooray, Li & Melchiorri 2008; Muñoz, Ali-Haïmoud & Kamionkowski 2015; Mondal et al. 2015; Mondal, Bharadwaj & Majumdar 2017), which can be directly estimated using the observed visibilities in a radio interferometer.

Another alternative to understanding and characterizing the signal more fully is through the images of the 21-cm signal (thus the H i distribution in the IGM) at different stages of reionization. This is more popularly known as ‘21-cm tomography’ (see e.g. Madau, Meiksin & Rees 1997), and it will essentially provide us with a time lapse movie of the state of the H i in a certain portion of the sky during the EoR. As radio interferometry is only sensitive to the fluctuations in the signal, not its absolute value, regions of ionized hydrogen will appear as negative ‘cold spots’ in these images. Thus, studying the properties of cold spots in the 21-cm tomographic data will provide a direct insight into ‘the morphology of reionization’, which is characterized by the shapes, sizes, and spatial distribution of the H ii regions. These studies will help us to directly quantify the ionized volume-filling fraction, the mean size, and the bubble size distribution of the H ii regions.

For this reason, the 21-cm tomographic imaging during the CD and the EoR has been identified as one of the major science goals of the future radio interferometer Square Kilometre Array⁵ (SKA) (Koopmans et al. 2015; Mellema et al. 2015; Wyithe, Geil & Kim 2015). Even currently operating interferometers such as LOFAR and MWA are expected to produce relatively coarser resolution images at a moderate SNR (Zaroubi et al. 2012; Malloy & Lidz 2013). However, the tomographic data cubes still need to be analysed in a statistical sense in order to derive astrophysical parameters.

The data analysis strategy of 21-cm tomographic images is a relatively unexplored field. Previous works have focused on topological analysis of 21-cm tomography through the Minkowski functionals (Gleser et al. 2006; Lee et al. 2008; Friedrich et al. 2011a; Hong et al. 2014; Yoshiura et al. 2016) or by extracting the size, shape, and location of the H ii region using the matched filter technique on the raw observed visibilities and images (Datta, Bharadwaj & Choudhury 2007; Datta et al. 2008; Datta, Bharadwaj & Choudhury 2009; Majumdar et al. 2011; Datta et al. 2012a, 2016; Majumdar, Bharadwaj & Choudhury 2012; Malloy & Lidz 2013). The latter have explored both targeted (where a very bright high-redshift ionizing source has already been identified through observations at other wavelengths) and blind (where the location of the ionizing source is not known) searches for H ii regions on the observed visiblities or on the image planes. For example, Malloy & Lidz (2013) have demonstrated that by using this matched filtering method, it is also possible to extract the size distribution of the 21-cm cold spots from an observational data set. Some additional work has been done on studying size distributions in simulations of reionizations (Friedrich et al. 2011b; Lin et al. 2016), but the application of these techniques to 21-cm tomography has not yet been explored.

One of the main challenges in recovering the statistical measures of the morphology (e.g. shape or size distribution) of the H ii regions from 21-cm tomography is due to the fact that there is no unique definition of the ‘size or shape of an H ii region’. Several methods and definitions of H ii region size and shape have been used to extract the size distribution from simulated 21-cm data (Iliev et al. 2006; McQuinn et al. 2007; Mesinger & Furlanetto 2007; Zahn et al. 2007; Friedrich et al. 2011a; Majumdar et al. 2014; Lin et al. 2016). So far, there does not seem to be a consensus on how these methods can actually trace the underlying size distribution of the |${\rm H\,{\small {II}}}$| regions. Apart from the ambiguity of the definition of an |${\rm H\,{\small {II}}}$| region, the other major challenge in obtaining the size distribution is to identify the ‘true’ ionized cells from all the cold spots present in the 21-cm tomographic data. As radio interferometers can only capture fluctuations in the signal, the tomographic data will have mean subtracted signal represented by positive and negative pixels in the image plane at a certain redshift. A pixel containing neutral gas but having density below the mean density could appear as a negative pixel. In addition, the mean of the 21-cm signal (which corresponds to the ‘zero’ of the signal in that data volume) will also vary with redshift as reionization progresses from its early to late stages. These two intrinsic features of the observed signal will produce additional confusion in identifying true ionized pixels from 21-cm images.

The main goal of this paper is to introduce a novel granulometric analysis method in the context of the EoR 21-cm tomography for the first time. We focus on measuring the |${\rm H\,{\small {II}}}$| region size distribution from 3D data cube and 2D images, hence, directly characterizing the aspect of the morphology of |${\rm H\,{\small {II}}}$| regions. This granulometric technique is based on a well-defined method in mathematical morphology and stochastic geometry (Serra 1983; Dougherty & Lotufo 2003; Chiu et al. 2013), which provides a well-formulated definition and the characterization of ‘size’ of the objects based on the idea of sieving (Matheron 1975). We further investigate the possibility and requirement to recover the |${\rm H\,{\small {II}}}$| region size distribution in 21-cm tomography from future radio interferometric observations.

This paper is organized as follows. Section 2 describes the 21-cm and reionization models. Section 3 introduces granulometry in the context of our data analysis methodology. Section 4 describes our mock data cube under a simple model for radio interferometry. Section 5 presents the theoretical aspects of the granulometric analysis of 21-cm tomography. Sections 6 and 7 present the observational prospects and requirements for measuring the |${\rm H\,{\small {II}}}$| region size distribution using 3D image cubes and 2D image slices from 21-cm tomography. The implications for baseline design and observing strategy are discussed in Section 8. Finally, our conclusions are summarized in Section 9.

2 EOR 21-cM SIGNAL

2.1 Models for 21-cm maps from EoR

To demonstrate the use of granulometric analysis in 21-cm tomography, we employ two models to characterize the distribution of ionized regions: (1) a toy model in which we impose a lognormal bubble size distribution and (2) the results from detailed radiative transfer (RT) simulations.

2.1.1 Lognormal bubble model

Finally, we use Monte Carlo realizations of the lognormal bubble model in a volume V_box = (1 h ^{− 1} cGpc)³ with 256³ cells and 3.9 h ^{− 1} cMpc resolution to create a differential brightness temperature map.

2.1.2 Radiative transfer simulation

To demonstrate our new data analysis methodology for 21-cm tomography, we use snapshots from a full RT simulation performed within the PRACE4LOFAR project.⁶ Reionization is modelled by post-processing an N-body simulation with an RT calculation. The N-body code used is cubep³m (Harnois-Déraps et al. 2013) and the RT code is c2ray (Mellema et al. 2006a). We assume a flat Λcold dark matter cosmology with the Wilkinson Microwave Anisotropy Probe 5 parameters |$h=0.7,\ \Omega _{\rm m}=0.27,\ \Omega _\Lambda =0.73,\ \Omega _{\rm b}=0.044,\ \sigma _8=0.8,\ {\rm and }\,n_{\rm s}=0.96$| (Komatsu et al. 2009), which are also consistent with Planck data (Planck Collaboration XIII 2015).

The simulation we use throughout this paper is similar to the LB1 model in Dixon et al. (2016), where details regarding simulation techniques and methods can be found. However, the volume considered here is a larger 500 h ^{− 1} cMpc on each side for which the N-body simulation was run with 6192³ particles with mass resolution 4.05 × 10⁷ M_⊙ (Dixon et al. in preparation). The minimum dark matter halo mass used in the RT simulation is 1 × 10⁹ M_⊙ (25 particles). The RT simulation was performed on smoothed and gridded density fields consisting of 300³ cells. The ionizing photons from the sources (galaxies) are assumed to linearly scale with the host halo mass such that the number of ionizing photons released into the IGM is g_γ = 1.7 per baryon per 10⁷ yr. This RT simulation of reionization runs from z = 21 to 6.

We use the resulting snapshots of the |${\rm H\,{\small {II}}}$| fraction and the gas density (on a 300³ grid) for the analysis presented in this paper. We select the z = 6.8 snapshot as our fiducial reference redshift unless otherwise stated. At this redshift, the volume-averaged ionized fraction corresponds to 〈x_H II〉_V = 0.40 (see Section 5.3 for a discussion on the dependence of our results on this choice).

3 DATA ANALYSIS METHOD

This section describes basic concepts involved in our data analysis methodology of 21-cm tomography using the granulometric technique. The same analysis method is used both for noiseless and noisy data.

3.1 Granulometry

Granulometry is a technique in mathematical morphology and image analysis that measures a size distribution of objects (Serra 1983; Dougherty & Lotufo 2003). The central idea is based on the concept of sieving (Matheron 1975). This provides a mathematically well-defined measure of the size distribution of a collection of objects, in our case |${\rm H\,{\small {II}}}$| regions and 21-cm cold spots, the latter defined as regions where the differential brightness temperature is less than a certain threshold value. The tool we used for our granulometric analysis is publicly available online.⁷

3.1.1 Binary images and data cubes

3.1.2 Size distribution measurement

The granulometric analysis technique is applied on the binary field. The basic idea is to probe the field with a so-called structuring element of a certain shape. This (loosely speaking) represents the shape of holes in a sieve. For our analysis, we choose the structuring element to be a sphere of radius R (or a disc of radius R when analysing a 2D image), defined as S_R = Π(r − r₀|R), where r₀ is the coordinate of the centre of the structuring element. The mathematical formulation of the concept of sieving corresponds to a morphological opening operation⁸ denoted by a symbol ∘ (e.g. Dougherty & Lotufo 2003). Hereafter, we refer to the morphological opening operation as sieving. A brief introduction is presented in Appendix A.

An advantage of the granulometric measure of the size distribution is that it attempts to recover the true underlying probability distribution function of sizes (Serra 1983). In Section 5.1, we verify this property using an explicit example. In the terminology of Lin et al. (2016), the granulometry method is an unbiased estimator of the true size distribution.⁹ We compare the granulometric method with other size estimators in Section 8.3.

3.1.3 Moments of size distributions and the normalization

3.2 Error estimation

So far, we have only introduced the theory for granulometric analysis. When measuring the |${\rm H\,{\small {II}}}$| region size distributions from observation of 21-cm cold spots, we must also assess the associated error as any measurement comes with uncertainty. The statistical uncertainties of the results come both from the sample variance of the signal and the thermal noise of the instrument. We estimate the error in the measured size distribution of 21-cm cold spots using an ensemble of mock data cubes with noise.

The total error covariance matrix is then given by |$C_{ij}=C_{ij}^N+C_{ij}^S$|⁠. The error is estimated according to the covariance matrix. Note that the sample variance error, which scales with the survey volume as |$\propto V_{\rm survey}^{-1/2}$|⁠, typically dominates. The thermal noise error has a more complicated dependency on the survey volume as it couples with the sample variance error (see Section 7.1), but it is usually a small contribution to the total error budget. We present 1σ error bars in this paper unless otherwise stated.

4 MOCK INTERFEROMETRIC OBSERVATIONS

In this section, we describe how we construct our mock data cubes, which represent the signal as observed with an idealized SKA-like radio interferometer.

4.1 Instrument

We consider EoR experiments with an SKA-like radio interferometer. Our reference array configuration is based on the SKA1-low baseline design documented by Dewdney (2016).¹⁰ SKA1-low has a frequency coverage between 50 and 350 MHz, corresponding to the redshift range 3 ≲ z ≲ 27. The assumed design consists of a total of 512 stations. The final configuration is still being developed but here we assume that approximately half of the stations are distributed within a radius of ∼ 1 km; therefore, we assume that the number of core stations is N_st = 256. We assume a complete uv-coverage within a baseline length of 2 km (we refer to this as the maximum baseline length b_max). The station diameter is taken to be D_st = 35 m. We assume an idealized instrument with an effective area which below the critical frequency ν_crit is the same as the geometric area of the station A_eff = π(D_st/2)², and above it falls off as (ν_crit/ν)² where ν_crit = 110 MHz. The system temperature is given by the sum of the receiver noise, T_rcvr = 40 K, and the sky temperature |$T_{\rm sky}=60(\nu /300{\,\rm MHz})^{-2.55}{\,\rm K}$|⁠. We refer to this idealized instrument as ‘SKA1-low’. The instrumental parameters are summarized in Table 1.

In addition, we consider a future SKA2-like instrument extending the core array of SKA1-low to a 2 km radius, and achieving the complete uv-coverage within the maximum baseline of b_max = 4 km. This improves the angular resolution by a factor of 2, but keeps the field of view (FoV) of the single pointing the same as the SKA1-low. We refer to this instrument as ‘SKA2’.¹¹

4.2 21-cm signal

The 21-cm signal is calculated following equation (2) using the results of the RT simulations described in Section 2.1.2. Our fiducial analysis uses the redshift snapshot at z = 6.8 (corresponding to a frequency of 182 MHz), where the simulation box has an extension of ∼ 4.6 deg on a side at this redshift. The spatial resolution is Δr = 1.67 h ^{− 1} cMpc. Each pixel therefore corresponds to an angular size of 0.92 arcmin. We construct a simulated 21-cm data cube using a coeval snapshot from the RT simulation.¹² Each 2D slice is one pixel width and it corresponds to a bandwidth B = ν₂₁H(z)Δr/[c(1 + z)²] ≈ 0.15 MHz, where H(z) is the Hubble parameter and c is the speed of light.

4.3 Angular and frequency resolution

The angular resolution of the radio interferometric observation is characterized by the maximum baseline of the array as θ_A = λ/b_max radians (see Table 1). For SKA2, the angular resolution is increased by a factor of 2. We implement the angular response of the interferometer by convolving with a Gaussian point spread function (PSF), R(θ), with full width at half-maximum (FWHM) corresponding to θ_A = 2.58(ν/200 MHz) ^{− 1} arcmin.

The frequency resolution is determined by the design of the instrument, and for SKA1-low it is expected to be better than 1 kHz. However, in practice, when analysing the signal a lower frequency resolution is used to increase the SNR. Here, we assume that the data are smoothed in the frequency direction with a Gaussian kernel of exactly the same physical size as the angular PSF. For the chosen redshift, this implies an FWHM of 453 kHz.

4.4 Noise

As explained in Section 3.2, we produce 100 Monte Carlo realizations of the noise cubes. We generate white (Gaussian) noise fields, which have the same spatial (angular) scale and frequency range as the 21-cm data cube (Section 4.2). The rms noise level 〈(ΔT_N)²〉 at the scale of the resolution element is then normalized according to equation (15). Because of the assumption of white noise, the noise power spectrum scales as |$\Delta _{\rm N}^2(k)\propto k^3$|⁠.

To quantify the image quality, we define the SNR for a data cube as the ratio of the rms fluctuations between an image cube and a noise cube on the scale of resolution element |${\rm SNR}(\theta _{\rm A})=\sqrt{\langle (\Delta T_{21})^2\rangle /\langle (\Delta T_{N})^2\rangle }$|⁠.

4.5 Foregrounds

We assume that the various foreground signals are perfectly removed from our data cube. Chapman et al. (2015) discussed the effect of different foreground removal techniques on the reconstructed 21-cm images, showing that good quality reconstructed 21-cm data cubes are in principle obtainable. Studying the impact of foreground residuals on the 21-cm tomographic analysis is beyond the scope of this paper.

5 GRANULOMETRIC ANALYSIS

In this section, we first present the results of granulometric analysis of one constructed and one simulated distribution of |${\rm H\,{\small {II}}}$| regions, as well as of the noiseless 21-cm signals associated with the latter. The goal is to understand the physical properties probed by the granulometric analysis and how well the 21-cm cold-spot size distribution traces the underlying size distribution of the |${\rm H\,{\small {II}}}$| regions.

5.1 A proof of concept: lognormal bubble model

As a proof of concept, we apply the granulometric measurement of size distribution of |${\rm H\,{\small {II}}}$| regions to a Monte Carlo realization of the 2D lognormal bubble model from Section 2.1.1. Fig. 1 shows an example of sieving the morphology of |${\rm H\,{\small {II}}}$| regions, represented by the white areas. The original distribution of |${\rm H\,{\small {II}}}$| regions, i.e. the binary image |$X_{\rm H\,{\small{II}}}(\boldsymbol {r})$|⁠, is shown in the left-hand panel. When it is sieved (X ∘ S_R) with a disc of radius R = 27 h ^{− 1} cMpc (middle panel), the structures smaller than the radius of the disc are removed. A larger disc of radius R = 39 h ^{− 1} cMpc (right-hand panel) sieves most of the structures and only some of the largest |${\rm H\,{\small {II}}}$| regions remain.

The histogram in Fig. 2 shows the differential size distribution measured with the sieving procedure. The differential size distribution clearly picks up the underlying probability distribution function of |${\rm H\,{\small {II}}}$| region sizes used by the lognormal model (dashed curve). The agreement is not perfect because of overlapping |${\rm H\,{\small {II}}}$| regions, which causes a departure from spherical shapes. Nevertheless, the lognormal model proves that the granulometric measure of the size distribution traces the underlying probability distribution of |${\rm H\,{\small {II}}}$| regions. From this, we conclude that the granulometric analysis is a very useful statistical tool to quantify the size distribution of |${\rm H\,{\small {II}}}$| regions with a mathematically well-motivated framework.

5.2 Size distribution of |${\rm H\,{\small {II}}}$| regions: RT simulation

There are physical motivations why the |${\rm H\,{\small {II}}}$| region size distribution follows the form of a modified Schechter function. First, suppose that the ultraviolet (1500 Å) luminosity function of galaxies responsible for driving reionization follows the Schechter function |$\phi (L_{1500})\propto (L_{1500}/L^\ast _{1500})^{\alpha _{\rm L}}\exp (-L_{1500}/L^\ast _{1500})$| with a characteristic luminosity |$L^{\ast }_{1500}$| and a faint-end slope α_L. The ionizing photon luminosity |$\dot{N}_{\rm ion}$| (in units of s ^{− 1}) at <912 Å from a galaxy is given by the product |$\dot{N}_{\rm ion}=f_{\rm esc}\xi _{\rm ion}L_{1500}$| (e.g. Robertson et al. 2013), where f_esc is the escape fraction and ξ_ion is the ratio between the 1500 Å luminosity and the intrinsic ionizing photon production rate of a galaxy. Assuming that all galaxies have their own |${\rm H\,{\small {II}}}$| regions during the early pre-overlapping phase, the bubble number density dn_b(R)/dR per unit radius of |${\rm H\,{\small {II}}}$| regions is given by |$\frac{{\rm d}n_{\rm b}(R)}{{\rm d}R}{\rm d}R=\phi (L_{1500}){\rm d}L_{1500}$|⁠. By estimating the radius of an |${\rm H\,{\small {II}}}$| region by counting photons (cosmological Strömgren sphere), |$R=[\frac{3\dot{N}_{\rm ion}t_{{\rm G}}}{4\pi \bar{n}_{\rm H}(0)}]^{1/3}\propto L_{1500}^{1/3}$|⁠, where t_G is the time interval during which a galaxy is ionizing the IGM and |$\bar{n}_{\rm H}(0)$| is the comoving mean number density of hydrogen atoms. Therefore, the bubble number density scales as |${\rm d}n_{\rm b}(R)/{\rm d}R\propto R^{3\alpha _{\rm L}+2}\exp [-(R/R_\ast )^3]$|⁠, where the characteristic size R_* is given by the ionizing properties of early galaxies, |$R_\ast =[\frac{3f_{\rm esc}\xi _{\rm ion}L^\ast _{1500}t_{{\rm G}}}{4\pi \bar{n}_{\rm H}(0)}]^{1/3}$|⁠. Hence, because dF(<R)/dR ∝ dn_b(R)/dR, we expect the |${\rm H\,{\small {II}}}$| region size distribution to scale as |$R\,{\rm d}F(<R)/{\rm d}R\propto R^{3(\alpha _{\rm L}+1)}\exp [-(R/R_\ast )^3]$|⁠, as long as there is not too much overlap between bubbles. This is indeed a form of the modified Schechter function (16). A second physical motivation comes from the result of the excursion set formalism for reionization, in which the mass function for the ionized bubbles has a form of the modified Schechter function (Furlanetto, Zaldarriaga & Hernquist 2004a).

This relation between the |${\rm H\,{\small {II}}}$| region size distribution and the properties of ionizing sources allows us to provide physical interpretations for the shape of the size distribution. For smaller |${\rm H\,{\small {II}}}$| region sizes, R < R_*, the distribution scales as a power law, |$R\frac{{\rm d}F(>R)}{{\rm d}R}\propto R^\alpha$|⁠. Our RT simulation gives a positive slope α = 0.71, while the above simple expectation gives a negative slope α = −3 for the faint-end slope of α_L = −2. The positive value of the small size end of |${\rm H\,{\small {II}}}$| region size distribution can be interpreted as a result of the well-known characteristic of the reionization process, i.e. that the cosmological |${\rm H\,{\small {II}}}$| regions grow by merging many small ones. Consequently, the small size end of |${\rm H\,{\small {II}}}$| regions is redistributed to larger sizes as reionization progresses. For larger |${\rm H\,{\small {II}}}$| region sizes, R > R_*, the interpretation is more complicated. The distribution shows an exponential drop off, scaling as ∝exp [( − R/R_*)^β], which could be (partially) due to the rapid decline in the population of luminous or clustered ionizing sources responsible for producing large |${\rm H\,{\small {II}}}$| regions. For a robust interpretation, a future study of the relation (if there is any) between the |${\rm H\,{\small {II}}}$| region size distribution and the properties and spatial distribution of ionizing sources is required; none the less, one may postulate that the |${\rm H\,{\small {II}}}$| region size distribution will contain the memory of the properties of ionizing sources.

5.3 Relation between |${\rm H\,{\small {II}}}$| regions and cold spots in 21-cm tomography

Observationally, we can only analyse 21-cm tomographic data. In tomographic images, the |${\rm H\,{\small {II}}}$| regions appear as cold spots, i.e. regions devoid of 21-cm emission. As only temperature fluctuations are observable in radio interferometry, those regions appear to have negatives values in the 21-cm images. Thus, intuitively the size distribution of cold spots could be directly translated into that of |${\rm H\,{\small {II}}}$| regions. However, a complication arises from the contribution of the density fluctuations to the 21-cm signal as low-density neutral regions can also give rise to negative values in the images (see equation 2). We would therefore like to answer the following question: how well do the volume-filling factor and the size distribution of 21-cm cold spots trace the volume-filling factor and size distribution of |${\rm H\,{\small {II}}}$| regions?

5.3.1 Volume-filling factors of |${\rm H\,{\small {II}}}$| regions and 21-cm cold spots

The 21-cm cold spots can encode the statistics of voids in matter density, as well as of ionized regions. To interpret the cold spots as a signature of |${\rm H\,{\small {II}}}$| regions, we should understand the nature of the void size distribution.

Fig. 4 shows the volume-filling factors of 21-cm cold spots (black solid line), |${\rm H\,{\small {II}}}$| regions (red dashed), and voids (blue dotted). The filling factors of 21-cm cold spots and |${\rm H\,{\small {II}}}$| regions are measured from the simulation box. Similarly, the filling factor of voids is calculated from the negative excursion sets δ_b(r) < 0 (underdense regions) of the density fluctuation field.

The figure shows that in our simulation the volume-filling factor of 21-cm cold spots traces that of |${\rm H\,{\small {II}}}$| regions well for 6 < z < 7. However, when the filling factor of |${\rm H\,{\small {II}}}$| regions drops below 0.4, the filling factor of 21-cm cold spots starts to deviate from that of |${\rm H\,{\small {II}}}$| regions. At z > 9, it approaches the filling factor of voids, Q_void.¹⁵ We refer to the difficulty in interpreting 21-cm cold spots due to the contributions from both |${\rm H\,{\small {II}}}$| regions and voids as the void–|${\rm H\,{\small {II}}}$| regions confusion.

The transition around a filling factor of 0.4 can be understood as follows. The level of the void contamination is controlled by the probability that the negative excursion of density fluctuation passes below the threshold given by the mean 21-cm signal 〈δT₂₁〉 = T₀〈x_H I(1 + δ_b)〉 ≈ T₀(1 − Q_H ii). Because the 21-cm signal from the neutral parts of the IGM is δT₂₁ = T₀(1 + δ_b), in order for the density fluctuations to appear as cold spots in the mean subtracted (observed) 21-cm signal, i.e. ΔT₂₁ = δT₂₁ − 〈δT₂₁〉 < 0, the density fluctuations must be lower than the threshold value, δ_b < −Q_H ii. The rms density fluctuations in our simulation are σ_b ≈ 0.37[(1 + z)/8]⁻¹. The volume-filling factor of |${\rm H\,{\small {II}}}$| regions Q_H ii ≈ 0.4 corresponds to this rms density fluctuation level. Below this value of Q_H ii, we expect a large contamination from the density fluctuations to 21-cm cold spots since many 1σ_b fluctuations can pass below the threshold level ≈−Q_H ii. This is the reason why the assumption that the 21-cm cold-spots filling factor is a good tracer of |${\rm H\,{\small {II}}}$| regions breaks down around Q_H ii ≈ 0.4(≈σ_b).

However, this is not a fundamental limitation of the granulometric analysis of tomographic data; the apparent confusion is associated with the somewhat arbitrary choice of the threshold value |$\Delta T^{\rm th}_{21}=0$| when creating a binary field. In a future improved granulometric analysis, the use of multiple threshold values should allow us to mitigate the apparent confusion to a large extent. In this introductory paper on granulometry however, we adopt a single threshold value |$\Delta T^{\rm th}_{21}=0$|⁠. In the following, we discuss the factors controlling the void–|${\rm H\,{\small {II}}}$| region confusion to help its mitigation in future work.

5.3.2 Size distributions of |${\rm H\,{\small {II}}}$| regions and 21-cm cold spots

Fig. 5 (top panel) shows the size distributions of 21-cm cold spots dQ₂₁(<R)/dR, |${\rm H\,{\small {II}}}$| regions dQ_H ii(<R)/dR, and voids dQ_void(<R)/dR measured by the granulometric analysis of the RT simulation at z = 6.8. At this redshift, the void–|${\rm H\,{\small {II}}}$| region confusion is small. The size distribution of voids is measured by applying the same granulometric analysis to the negative excursion sets of the density fluctuation field, δ_b(r) < 0. The bottom panel shows the difference between the size distributions of 21-cm cold spots and |${\rm H\,{\small {II}}}$| regions. Note that the size distributions shown in Fig. 5 (and hereafter unless otherwise stated) are normalized to the volume-filling factor, i.e. dQ(<R)/dR ≡ QdF(<R)/d R.¹⁶

As suggested by the results on the filling factors in the previous section, the size distribution of 21-cm cold spots traces that of |${\rm H\,{\small {II}}}$| regions very well for Q_H ii > 0.4. The slight deviation seen at R ≲ 5 h ^{− 1} cMpc occurs because of the contamination by voids. This leads to a slight underestimation of the mean size of |${\rm H\,{\small {II}}}$| regions inferred from the mean size of 21-cm cold spots (by ∼10 per cent). The voids remain a minor contaminant to the overall 21-cm cold-spot size statistics of 21-cm tomography affecting the results only by approximately 10 per cent.

Fig. 6 is the equivalent of Fig. 5 at z = 7.6 (〈x_H II〉_V = 0.16). As expected, we see that the contamination from voids dominates over the signature of |${\rm H\,{\small {II}}}$| regions in the 21-cm cold-spot size distribution. As noted in Section 5.3.1, this is because the underlying filling factor of |${\rm H\,{\small {II}}}$| regions is small during the first half of reionization.

The large difference between the shapes of the size distributions of |${\rm H\,{\small {II}}}$| regions and voids provides a way to avoid the void–|${\rm H\,{\small {II}}}$| region confusion. The void size distribution is always confined within R ≲ 5 h ^{− 1} cMpc. In fact, the size distribution of the excursion sets of a Gaussian random field can be well understood and predicted given a priori knowledge of the matter power spectrum (Bardeen et al. 1986; Bond & Efstathiou 1987; Sheth & van de Weygaert 2004). Therefore, we can test the robustness of the identification of cold spots as |${\rm H\,{\small {II}}}$| regions against the null hypothesis of void size statistics.

Overall, we conclude that the size distribution of 21-cm cold spot traces that of |${\rm H\,{\small {II}}}$| regions very well during the second half of reionization. For the first half of reionization, the interpretation of 21-cm cold spots becomes increasingly difficult because of large contamination from voids for our canonical choice of threshold |$\Delta T_{12}^{\rm th}=0$|⁠. We note that this may be mitigated by choosing lower values for the threshold.

6 RECOVERY OF H ii REGION SIZE DISTRIBUTION: 3D DATA SETS

So far we have only considered the case of a pure simulated 21-cm signal. We now turn our attention to the prospects for recovering the |${\rm H\,{\small {II}}}$| region size distribution from 21-cm tomography using a SKA-like radio interferometer.

The recovery of the size distributions of 21-cm cold spots and |${\rm H\,{\small {II}}}$| regions from real-world radio interferometric observations is subjected to noise, instrumental response, foregrounds, and observing strategy. We therefore ask the question: what are the requirements to recover the |${\rm H\,{\small {II}}}$| region size distribution from 21-cm tomographic data observed with the SKA?

6.1 Effect of the field of view: sample variance error

The first requirement addresses the FoV or sky coverage. Fig. 7 shows the effect of a finite FoV on the observed size distributions of 21-cm cold spots measured from the granulometric analysis of our mock image cube of 4.6² deg² FoV with frequency width 44 MHz (red dashed lines), and 2.3² deg² FoV with frequency width 22 MHz (black solid), using the SKA1-low with the limiting case of noise-free data. The frequency width is chosen to match with the size of the FoV so that the comoving lengths of the line-of-sight and perpendicular directions are equal. The error bars correspond to 1σ uncertainty due to the sample variance (Section 3.2). The red and black curves are the best-fitting modified Schechter functions, and the vertical dotted line indicates the FWHM of the angular resolution. Note however that the cubes we have analysed have a large-frequency width which implies that the real data will be subject to the light-cone effect (Barkana & Loeb 2006; Datta et al. 2012b). The results in this section should therefore be interpreted as four (one) pointed observations of an FoV of 4.6² deg² (2.3² deg²) with a frequency width of approximately 10 MHz, for which the light-cone effect can most likely be neglected (Datta et al. 2014).

The figure shows that the smaller FoV (2.3² deg²) does not introduce any significant systematic bias compared to our fiducial FoV of 4.6² deg². The slight under(over)estimation of the large(small) sizes of 21-cm cold spots remains within 10 per cent–20 per cent at most, which is within the 1σ error due to sample variance. The sample variance error increases for larger 21-cm cold-spot sizes because the number of large |${\rm H\,{\small {II}}}$| regions within an FoV is smaller. Note that below the angular resolution scale the error appears to be negligibly small. This is because the angular smoothing artificially induces any cold spots at that scale to have always the same size, and therefore the FoV-to-FoV variation is suppressed.

No systematic bias will be introduced as long as the FoV is large enough to include the largest scale of |${\rm H\,{\small {II}}}$| regions (the theoretical prediction for the characteristic size is of order of tens of comoving Mpc, e.g. ∼ 0.5 deg for ∼ 50 h ^{− 1} cMpc. See, for example Wyithe & Loeb 2004; Geil et al. 2015; Dixon et al. 2016). Therefore, the single-pointing FoV of SKA1-low ( ∼ 3 deg) is large enough to measure the 21-cm cold-spot size distribution with a reasonable sample variance error when the 3D image cube is used.

6.2 Effect of angular resolution: smoothing bias

The second requirement comes from the finite angular resolution of the instrument, which strongly affects the recovery and interpretation of the 21-cm cold-spot and |${\rm H\,{\small {II}}}$| region size distributions.

Fig. 8 shows the effect of finite angular resolution on the 21-cm cold-spot size distribution at redshift z = 6.5 and 6.8 for a noiseless image cube with a 4.6² deg² FoV and a frequency depth of 44 MHz, using an SKA1-low (blue squares) and SKA2 (green triangles) configuration. The frequency resolution is set to match the angular resolution in comoving coordinates. The size distribution of the pure simulated 21-cm tomography (red circles) represents the case of an ideal, fully resolved observation. The vertical lines indicate the angular resolutions of SKA1-low (2.8 arcmin, blue dotted) and SKA2 (1.4 arcmin, green dashed).

The angular resolution introduces a systematic bias in the observed size distribution of 21-cm cold spots. Both at z = 6.5 and 6.8, the shape of the size distributions is systematically shifted towards larger sizes. This is somewhat counterintuitive as an obvious expectation is the suppression of the small end of a size distribution. This systematic bias occurs because the angular smoothing (or Gaussian PSF) actually mixes and merges many 21-cm cold spots (including both |${\rm H\,{\small {II}}}$| regions and voids). As a result, many small cold spots become a large one. We call this important effect the smoothing bias. Since the granulometric analysis is performed on the binary field satisfying ΔT₂₁ < 0, smoothing does not always remove small cold spots. Instead, because of the mixing of cold spots, intrinsically small cold spots re-populate the larger part of a size distribution. Quantitatively, the impact of the smoothing bias is shown in the bottom panels of Fig. 8, where the difference between the simulation and the mock observations is plotted. The smoothing bias is larger for a larger angular smoothing scale. Importantly, the figure shows that the smoothing bias impacts all scales of the size distribution, even when the angular resolution scale is below the scale of interest.

The impact of a finite angular resolution is less severe when the characteristic size or mean size of the underlying (true) 21-cm cold-spot size distribution is much larger than the angular resolution, θ_A ≪ R_*/D_A(z) [or ≪〈R〉/D_A(z)], where D_A(z) is the angular diameter distance to redshift z. For example, at 〈x_H ii〉_V = 0.40, the angular resolution causes a systematic bias in size distribution by ∼50(20) per cent fractional error for SKA1-low (SKA2), while at 〈x_H ii〉_V = 0.63, the bias is only about ∼10 per cent. We note that a chromatic PSF will introduce an extra complication (Vedantham, Udaya Shankar & Subrahmanyan 2012) and possibly a systematic bias. However, addressing this is beyond the scope of this paper.

This analysis suggests that the effect of angular resolution must be taken into account when interpreting the 21-cm cold-spot size distribution. If an unbiased measurement (to the accuracy of order 10 per cent) of the size distribution is required, we should correct the angular resolution bias by a: (i) ‘baseline-design’ approach, where we ensure that the maximum baseline is large enough that the angular resolution bias in the size distribution and the characteristic size of cold spots remains sufficiently small, or a (ii) ‘forward-modelling’ approach, where we create a large suite of models and mock data cubes and directly fit the simulated size distribution in the observation using a Markov Chain Monte Carlo approach. However, as we will show in Section 8, a sufficiently high angular resolution must be achieved for both methods to work.

6.3 Effect of thermal noise: splitting bias

The third requirement is on the thermal noise. Fig. 9, showing a 2D image from the data cube, visually illustrates the effect of noise on the cold spots in 21-cm tomography. The thermal noise (middle panel) acts as an additive noise, and contaminates the underlying 21-cm cold spots (left-hand panel) by distorting the shape and size of the |${\rm H\,{\small {II}}}$| regions (right-hand panel).

The effect of thermal noise on the observed 21-cm cold-spot size distribution is shown in Fig. 10, which plots the measured size distributions from a mock 3D image cube with SKA1-low (21-cm+noise, black solid), noiseless 21-cm signal (21 cm, red dashed), and noise cube (noise, blue dotted). The 3σ error due to the thermal noise (inner error bars) and the total uncertainty including the 3σ sample variance error (outer error bars) are shown. The left-hand (right) panel corresponds to the two mock observations with SKA1-low having an rms noise level of 4.35 mK (2.0 mK) per resolution element for an interferometer with a 2 km maximum baseline.

First, as shown in the left-hand panel in Fig. 10, an SKA1-low with σ_N = 4.35 mK permits us to measure the 21-cm cold-spot size distribution using a low SNR 3D image cube (black solid) within ≲50 per cent of the noiseless case (red dashed).

The thermal noise, however, introduces another important systematic bias. The size distribution measured from the low SNR image cube is shifted systematically towards lower sizes beyond the statistical error. This bias occurs because the noise splits up the cold spots into smaller ones when the positive excursion sets of the noise fluctuation occur inside 21-cm cold spots. We thus call this systematic bias the splitting bias. The bottom panels in Fig. 10 indeed show that originally large 21-cm cold spots are split and re-populate the small end of the size distribution.

To understand the mechanism of the splitting bias, Fig. 11 shows the 1D profiles of the brightness temperature contrast for the same three data sets as above. The cold spots B and D are split into two separate ones in the mock observation. This occurs when the positive excursion sets of noise exceed the contrast of the underlying 21-cm cold spots, ΔT_coldspot ≈ 5-10 mK. Because the noise is a Gaussian random field, the probability that the positive excursion of noise exceeds the 21-cm cold-spot contrast is |$P(>\Delta T_{\rm coldspot})=\frac{1}{2}[1-{\rm erf}(\frac{\Delta T_{\rm coldspot}}{\sqrt{2}\sigma _{\rm N}})]\approx 1-12\hbox{\,per\,cent}$| for ΔT_coldspot ≈ 5-10 mK, where erf is the error function. This is a noticeable effect for σ_N = 4.35 mK (as shown in the left-hand panel of Fig. 10) where the ∼2σ peaks of noise contaminate the signal inducing the splitting bias. On the other hand, this means that only high sigma peaks act as contaminants to the granulometric measurement of the size distribution. This is the reason why a low SNR image cube still permits a reasonable recovery of the underlying size distribution of 21-cm cold spots.

Note that an opposite effect also occurs: a creation of spurious cold spots due to the negative excursion sets of noise (for example see the cold spot C in Fig. 11). However, because these spurious cold spots are always small (R ≲ 5 h ^{− 1} cMpc), the net effect of noise is a systematic bias shifting the size distribution towards smaller sizes.

At a noise level of σ_N = 2.0 mK (i.e. SNR = 3.1), shown in the right-hand panel of Fig. 10, the splitting bias becomes much smaller. The probability that contaminating noise peaks occurs for 21-cm cold spots with ΔT_coldspot ≈ 5-10 mK is only P(>ΔT_coldspot) < 1 per cent.

The level of the systematic bias entering into the cold-spot size distribution measurement does not scale linearly with the noise level (and hence, with the integration time). There is a jump in the quality of the recovery of the 21-cm cold-spot size distribution as a function of the rms noise/integration time. When the rms noise passes the threshold required by the expected brightness temperature contrast of 21-cm cold spots, the quality of the measurement jumps up. Above this threshold, increasing the integration time gradually improves the quality, but the gain is not as pronounced.

The statistical uncertainty of the thermal noise is negligibly small and is shown as the small inner error bars in Fig. 10. This is because in a large FoV (4.6 × 4.6 deg²) the variation of cold-spot sizes due to the thermal noise tends to average out. As the large end of the size distribution has fewer cold spots to average out the noise error, the statistical error bars increase slightly at larger sizes. On the other hand, the error bars become negligibly small for smaller cold spots as there are many samples.

7 RECOVERY OF H ii REGION SIZE DISTRIBUTION: 2D DATA SETS

In the previous section, we have considered a data cube that contains the full 3D information about the 21-cm signal. In real-world observations, such full 3D information representing a fixed redshift may not be available due to the light-cone effect and a bad quality or noise contamination in some frequency slices which must be discarded from the final analysis. Therefore, in this section, we consider more restricted observations of 2D images. We examine the prospects and requirements for recovering the |${\rm H\,{\small {II}}}$| region size distribution from a 2D tomographic image. Choosing only one frequency slice may be too drastic. However, we use this case to illustrate the effects of reducing the frequency width of the 21-cm data set.

The qualitative effects of the FoV, angular resolution, and noise on 2D images are the same as the case for 3D data cubes, i.e. (1) a finite FoV introduces a statistical error due to sample variance, but no systematic bias; (2) a limited angular resolution introduces a smoothing bias; and (3) thermal noise introduces a splitting bias, but a small statistical error. However, quantitatively the situation differs because the sample variance becomes extremely large for a 2D image and dominates the errors. Therefore, in the following sub-sections, we examine in detail the role of sample variance in 2D images.

7.1 Effect of the FoV for 2D images

The error on the 2D size distribution measurement is sample variance dominated. Fig. 12 shows the effect of a finite FoV on the 2D size distribution of 21-cm cold spots obtained from 2D images. 2D radii, R_⊥, of cold spots are defined as the apparent perpendicular lengths of the cold spots’ radii in the 2D image. We analyse a single-frequency image which corresponds to a channel width of 0.45 MHz. For a better recovery of the size distribution (black solid curve), the granulometric analysis of a larger FoV 21-cm image (red solid histogram, 4.6² deg² FoV) is preferred. Averaging over many 2.3² deg² FoV single-pointing images (black dashed curve) works equally well. The bottom panel shows that the sample variance error of 4.6² deg² FoV data is smaller than for the 2.3² deg² FoV data by a factor of 2, because of the factor of 4 increase in the area probed. Although the statistical error due to sample variance differs only by a factor of 2, the large end (R_⊥ ≳ 10 h ^{− 1} cMpc) of the size distribution is difficult to determine by a single-pointing FoV image. A large variation (due to the difficulty in sampling many large cold spots) could produce a catastrophic error in the measured 2D size distribution (see blue dashed histogram). While such a big error can also occur for large FoV data, the probability is much lower. Therefore, we conclude that increasing the area of the sky probed by increasing the FoV by an interferometric mosaicking or multibeaming (multiple EoR windows) technique will be important when analysing single-frequency 21-cm images (see Section 8 for discussion).

Having a large effective FoV is also important for controlling the thermal noise uncertainty. The statistical uncertainty due to the thermal noise is smaller for a large FoV data as the thermal noise error couples with the FoV. This is because the noise randomly modifies the shape and size of 21-cm cold spots, introducing a statistical uncertainty in the measured size distribution (see Section 6.3). For a fixed rms noise level, a larger FoV allows the random modifications to be averaged out because there are many cold spots. As the number of larger cold spots is smaller in a small FoV data, the noise error increases for larger cold-spot sizes.

Therefore, having a large sky coverage by mosaicking/multibeaming is also advantageous for reducing the statistical error bars due to the thermal noise. Note that, however, a large sky coverage does not reduce the splitting bias (Section 6.3), which is only reduced by having a lower absolute value of the rms noise level after a longer integration time.

7.2 Performance of a successful recovery

Finally, we test the performance of the granulometric measurement of the 21-cm cold-spot size distribution once all the requirements are satisfied. Fig. 13 shows the expected results for SKA1-low about the granulometric measurement of the 21-cm cold spots using a mock 4.6² deg² FoV image at z = 6.8 after interferometric mosaic imaging by patching many single-pointing data. The image is formed with a 0.45 MHz channel width. The black points show the mock measurements and the black curve shows the best-fitting modified Schechter function. The red curve and histogram refers to an ideal case in the absence of noise. The blue histogram shows the spurious noise cold-spot size distribution.

With an rms noise σ_N = 4.35 mK (SNR = 1.4 at the resolution element), as shown in the left-hand panel, the splitting bias due to thermal noise is still noticeable on the measured size distribution (black). The splitting bias is within the statistical uncertainties when both sample variance and noise errors are included. For the SKA1-low data with σ_N = 2.0 mK,¹⁷ the splitting bias becomes negligibly small. Thus, successfully recovering the 21-cm cold-spot size distribution at 2 km baseline resolution level is possible with the SKA1-low when the mosaicking/multibeaming technique is used to increase the effective sky coverage.

Fig. 14 shows the expected improvement for SKA2 if it is extended to longer intermediate-scale baselines. The green points are the mock SKA2 measurements, and the green curve with triangles indicates the best-fitting modified Schechter function. As already discussed in Section 6.2, we should note that the finite angular resolution introduces a fundamental instrumental limitation for SKA1-low (the best-fitting modified Schechter function is shown as the blue curve with squares). The smoothing bias systematically shifts the measured size distribution to sizes larger than in the case of a fully resolved 21-cm signal (red curve with open circles). The bias could be as large as a factor of ∼2 for SKA1-low, which exceeds the total statistical uncertainties. For SKA2, assuming the complete uv-coverage extends to 4 km baselines, the data with σ_N ≈ 2 mK noise at 1.4 arcmin resolution reduce the smoothing bias. At R_⊥ > 6 h ^{− 1} cMpc, i.e. well above the angular resolution smoothing scale (FWHM ≈ 2 h ^{− 1} cMpc), the SKA2 will be able to recover the underlying true cold-spot size distribution (red curve) well within the statistical uncertainties using a 4.6² deg² FoV and SNR ≈ 3 image. At this quality, the sample variance error is the dominant source of uncertainty. Thus, increasing sky coverage is necessary to measure the size distribution better. For the smaller end of the size distribution, increasing angular resolution is necessary to prevent a large smoothing bias.

8 DISCUSSIONS

In this section, we discuss the implications and the importance of measuring the |${\rm H\,{\small {II}}}$| region size distribution from 21-cm tomography to understand the EoR in a wider context. A possible observing strategy and baseline design are also discussed.

8.1 Synergy between 21-cm power spectra and tomography

Mellema et al. (2015) and Koopmans et al. (2015) speculated about the synergy between a 21-cm power spectrum analysis and tomography. Here, we present an explicit example to support their arguments.

Fig. 15 shows two noiseless 21-cm tomographic images with SKA1-low using the RT simulation and a Gaussian random field which has a 21-cm power spectrum identical to the one from the RT simulation. By eye, we can discern clear morphological differences between the two images. However, given that their power spectra are identical, how can we tell whether they are two different realizations of the same reionization process or intrinsically different models?

Fig. 16 shows that the 21-cm cold-spot size distribution measurement from 21-cm tomography using the granulometric analysis has the capability to break this degeneracy in 21-cm power spectrum measurements. This clearly illustrates, although for an extreme example, that tomographic observations add valuable information to the measurement of the 21-cm signal.

Of course, the information from the power spectrum and tomography are complementary to each other. For example, the 21-cm power spectrum contains information about the clustering of |${\rm H\,{\small {II}}}$| regions, which is not captured by size distribution. Furthermore, the higher order multipoles of the power spectrum can potentially quantify the strength and nature of correlation between the underlying matter distribution and the |${\rm H\,{\small {II}}}$| region distribution (Majumdar et al. 2016a,b). Therefore, we expect a clear synergy between radio interferometric observations of 21-cm power spectra and imaging tomography. Hence, we place an emphasis on optimizing and balancing the design of future 21-cm experiments for both power spectrum and tomography to achieve an optimal scientific return.

8.1.1 Role of intermediate baselines

How well tomographic data can distinguish different models will of course depend on the (angular) resolution of the images, which is set by the length of the longest baselines used in the construction of the tomographic data. Fig. 17 shows the effect of angular resolution on the 21-cm power spectrum and the cold-spot size distribution measured from 21-cm tomography. It shows the size distributions and power spectra measured with resolutions ≈ 5.70, 2.85, and1.43 arcmin. These correspond to maximum baseline lengths of 1 km (red), 2 km (blue), and 4 km (green) at z = 6.8. Note that all angular resolutions used are smaller than or comparable to the characteristic size [ ∼ 10 h ^{− 1} cMpc ( ∼ 5.5 arcmin)] of |${\rm H\,{\small {II}}}$| regions.

For a mock observation using a maximum baseline of 1 km, the size distributions of both the RT simulation and the Gaussian random field start to resemble each other, as they converge to the same shape with decreasing angular resolution. In fact, the only difference is due to the volume-filling factor of cold spots. A low angular resolution (large smoothing scale) dictates the shape of the size distribution dF(<R)/dR [and dQ(<R)/dR]. The signature of underlying |${\rm H\,{\small {II}}}$| region size distribution is erased, or not significant. While the particular level of angular resolution required to distinguish models differs depending on the stage of reionization, high angular resolution is clearly necessary for differentiating the shapes of size distributions. The power spectra are obviously indistinguishable by design at all angular resolutions.

Therefore, the merit of intermediate baselines and high angular resolution for this particular data set is clear from Fig. 17. Here, the measurement of the |${\rm H}\,{\small {II}}$| region size distribution can only break the degeneracy in the 21-cm power spectrum measurements when baselines longer than 2 km are used. Of course, this conclusion depends on the actual size distribution of ionized regions. For later stages, where the size distribution peaks at larger sizes, maximum baselines of 1 km would still suffice.

8.2 Observing strategy and baseline design

Overall, we suggest that a gradual roll-out of longer intermediate baselines in the SKA2 phase will be beneficial to quantify the reionization morphology from 21-cm tomography. Since the granulometric size distribution measurement can work with a moderate SNR (∼3) imaging, it is more important to sample a large range of k modes in order to eliminate the systematic biases than to increase the SNR of a limited range of large-scale k modes (which may be preferred for power spectra analysis). In this regard, a tiered radio survey (Koopmans et al. 2015) combining a deep narrow field (a single beam pointing) and a shallow wide field (multibeaming) could also benefit the measurement of |${\rm H\,{\small {II}}}$| region size distributions from 21-cm tomography. The tiered observation would deliver a layer of images with a sufficient SNR and FoV for a given angular scale for the |${\rm H\,{\small {II}}}$| size distribution measurement. Therefore, using the deep-narrow- and shallow-wide-field images, we may be able to separately measure the small and large ends of the size distribution, which could make effective use of the integration time and compromise with other science cases such as the power spectrum measurement. In addition, we expect a sweet spot in the baseline configuration, as the splitting bias from noise and the smoothing bias from finite angular resolution have an opposite effect on the size distribution. Further investigation of optimal observing strategy and baseline design for 21-cm tomography will certainly benefit the SKA EoR science cases. Naturally, the first aim should be a detection of a direct 21-cm image at large scales with compact dense core stations, and a large-scale pilot imaging is certainly valuable for testing the imaging algorithm/calibration/foreground removal for 21-cm tomography. Once the pilot phase is over, a larger gain will be achieved by extending to longer intermediate baselines, which would provide a robust physical interpretation and the unbiased size distribution of |${\rm H\,{\small {II}}}$| regions.

8.3 Comparison with other size measures

The literature contains several other algorithms for characterizing the size distribution of |${\rm H\,{\small {II}}}$| regions during reionization, e.g. (1) volumes of topologically connected regions, or the equivalent radii of these (a.k.a. the friends-of-friends method, Iliev et al. 2006; Shin, Trac & Cen 2008); (2) the largest radius at which the average ionized fraction around a pixel is above a certain threshold (e.g. 90 per cent, spherical average method, Zahn et al. 2007; McQuinn et al. 2007); (3) distance from a random pixel in |${\rm H\,{\small {II}}}$| regions to the edge of the |${\rm H\,{\small {II}}}$| region containing the pixel (mean-free path method, Mesinger & Furlanetto 2007); (4) equivalent radius of Watershed basin (Watershed algorithm, Lin et al. 2016).

In the terminology introduced by Lin et al. (2016), the mean-free path and spherical average methods are both biased and diffusive in that even in the case of an ensemble of equally sized, non-overlapping, spherical |${\rm H\,{\small {II}}}$| regions, they produce a range of sizes (Lin et al. 2016). The other two methods, just as granulometry, would reproduce the correct, unbiased size distribution.

The friends-of-friends method focuses on the volume of connected regions, which means that due to the percolative nature of reionization, already at around 10 per cent ionization most of the ionized volume will be part of one large region. It therefore is more suitable to test the percolation process (Furlanetto & Oh 2016).

The watershed method was recently proposed by Lin et al. (2016). It is commonly used in image processing and it segments the data into discrete bubbles. Lin et al. (2016) showed it to produce good bubble size distributions. However, unlike granulometry, it requires an additional tunable parameter to suppress Poisson noise.

The reason that so many different methods have been proposed is because there is no unique way to capture with a simple size distribution the complex morphology of ionized regions during reionization. Friedrich et al. (2011a) and Lin et al. (2016) compared these different methods and explored the similarities and differences between their results. Malloy & Lidz (2013) showed that the matched filter technique can be used to measure the bubble size distribution although the associated systematic bias is still to be understood. Each of these methods can in principle be used to derive its answer for the size distribution in real data and use this as a metric to compare to simulation results. However, because they each use a different approach we should not compare the results between different methods. Future studies will reveal which method(s) are most useful to constrain reionization parameters from 21-cm tomographic data.

9 CONCLUSIONS

Using 21-cm tomography, one can directly measure a fundamental quantity of the reionization process, i.e. the size distribution of the |${\rm H\,{\small {II}}}$| regions, in a model-independent way. The central questions that we have dealt with in this paper are – what can we learn from 21-cm tomography using a next generation radio interferometer such as the SKA? What are the observational requirements for achieving a good science return from such a 21-cm tomography? We summarize our main conclusions as follows.

• Granulometric analysis of the 21-cm tomographic data allows the measurement of the |${\rm H\,{\small {II}}}$| region size distribution.

   We have introduced a novel technique, called ‘granulometry’, to quantify the |${\rm H\,{\small {II}}}$| region size distribution in a mathematically well-formulated framework, but which is in practice rather simple to implement. The technique attempts to trace the underlying probability distribution function of the |${\rm H\,{\small {II}}}$| region sizes. This places the previously not so well-defined concept of ‘the morphology of reionization’ on a firm mathematical foundation. Using mock interferometric observations of RT simulations, we have shown that the granulometric analysis of 21-cm tomographic data allows us to recover the |${\rm H\,{\small {II}}}$| region size distribution with the SKA; the theoretical and observational systematics, requirements, and observing strategies have also been examined in detail.

   The measured size distribution is well described by a modified Schechter function. Even with our simplest application of the granulometric analysis, the 21-cm cold spots work as an excellent tracer of |${\rm H\,{\small {II}}}$| regions during the second half of reionization (i.e. |$Q_{{\rm H\,{\small {II}}}} \ge 0.4$|⁠). The size statistics of cold spots can therefore be directly interpreted as that of the |${\rm H\,{\small {II}}}$| regions. For the first half of reionization, attention must be paid to the possible confusion between |${\rm H\,{\small {II}}}$| regions and the voids of density fluctuations in the 21-cm cold-spot size distribution. Although this is not a fundamental limitation of the granulometric analysis, an improved method to separate them for a physical interpretation of the |${\rm H\,{\small {II}}}$| region size statistics must be used in future work.

• Observational requirements to recover the |${\rm H\,{\small {II}}}$| region size distribution from 21-cm tomography are attainable with the SKA. An ideal observing strategy would be a moderate SNR, wide-field mosaic/multibeamed imaging with additional longer intermediate baselines ( ∼ 2-4 km).

   To measure the cold-spot size distribution using 21-cm tomography of an actual radio interferometric observation, we considered two limiting cases: (1) the measurement from a 3D image cube and (2) the measurement from a 2D image slice. The fundamental requirements for recovering the true cold-spot size distribution is to have a high angular resolution (i.e. intermediate baselines of ∼ 2-4 km) and low sample variance (i.e. large sky coverage), both of which are achievable by the SKA by employing an appropriate observational strategy. Each of the requirements and systematics are discussed below.

   A finite angular resolution, even if it is below the characteristic scale of |${\rm H\,{\small {II}}}$| regions, introduces the most important source of systematic bias (‘smoothing bias’) for the measured 21-cm cold-spot size distribution at all scales. This smoothing bias skews the true size distribution to larger sizes. Although the current baseline distribution of SKA1-low (maximum 2 km baseline) limits the bias below approximately ∼50 per cent, a higher angular resolution is desirable for an accurate unbiased measurement. Controlling the smoothing bias is a key to successfully recover the statistical characterization of sizes and morphology of |${\rm H\,{\small {II}}}$| regions.

   The sample variance due to a finite FoV and frequency sampling is a main source of statistical uncertainties. Analysis of a 3D image cube suffers little from the finite FoV because of the large-frequency samples. The data from a single-pointing observation of SKA1-low allows us to measure the cold-spot size distribution with a small statistical uncertainty. On the other hand, for the analysis of a 2D image slice from 21-cm tomography, the sample variance must be reduced to small enough values for a reasonable measurement. A single-pointing observation suffers from a large sample variance in this case. In fact, because of the propagation of sample variance to the statistical uncertainty of the thermal noise, lowering the sample variance is also advantageous to reduce the statistical error. Mosaicking/multibeaming techniques for 21-cm imaging tomography will be desirable to increase the sky coverage (i.e. to lower sample variance).

   The thermal noise is not a major obstacle. This can be already manageable with SKA1-low. A moderate SNR (≲3) of 21-cm images still permits the recovery the cold-spot size distribution within ∼50 per cent. The required rms noise level is ∼ 4 mK, which should be achievable by the SKA1-low data. The largest observational systematics is the bias, instead of a statistical error, caused by the thermal noise (‘splitting bias’). The splitting bias artificially skews the measured size distribution to smaller sizes. However, its effect becomes negligible even for SNR ≈ 3 imaging (with the rms noise level of ∼ 2 mK). A long, but manageable, integration time reduces the effect of noise to a negligible level.

• Requirements for the science beyond power spectra and the synergy between 21-cm power spectrum and tomography.

Our results put a significant importance on having intermediate baselines in radio interferometers. This important factor should be taken into account for the design of a future extension of SKA1-low to SKA2, and other radio telescopes. The intermediate baselines are needed to fully exploit the synergy between 21-cm power spectrum and tomography. For science beyond power spectra, the availability of intermediate baselines is a pre-requisite for distinguishing the size distributions and morphology of |${\rm H\,{\small {II}}}$| regions across different reionization models (Section 8.1).

The intermediate baselines are necessary for reducing a fundamental instrumental limitation by the angular smoothing bias. In addition, they may be beneficial for overcoming a more immediate challenge of calibration such as the scintillation noise due to ionospheric turbulence. Provided that the future extension of SKA will invest on additional longer intermediate baselines, a promising observing strategy for 21-cm tomography is an interferometric mosaicking/multibeaming imaging.

Finally, the granulometric analysis introduced in this paper is only the tip of the iceberg of the entire spectrum of tools in mathematical morphology and stochastic geometry. They provide powerful means for quantifying the morphology of reionization based on a firm mathematical foundation and theory. Armed with this foundation, the synergy between 21-cm power spectrum and tomography provides many opportunities for directly probing the reionization morphology and extending 21-cm science beyond power spectrum.

Acknowledgments

KK acknowledge support from Richard Ellis and the European Research Council Advanced Grant FP7/669253. KK would also like to thank for the hospitality offered by the Department of Astronomy at Stockholm University during one visit, where this work was initiated. Some of the work was also conducted during the workshop supported by the Munich Institute for Astro- and Particle Physics (MIAPP) of the DFG cluster of excellence ‘Origin and Structure of the Universe’. SM would like to acknowledge financial assistance from the European Research Council under ERC grant number 638743-FIRSTDAWN. We acknowledge PRACE for awarding us computational time under PRACE4LOFAR grants 2012061089 and 2014102339 and access to resource Curie based in France at CEA and to resource SuperMUC at LRZ. This work was supported by the Science and Technology Facilities Council (grant numbers ST/F002858/1 and ST/I000976/1); and the Southeast Physics Network (SEPNet).

REFERENCES

APPENDIX A: BASIC OPERATIONS IN MATHEMATICAL MORPHOLOGY

Mathematical morphology (Matheron 1975; Serra 1983) defines a set of operations to formulate an algebra of shapes. In this appendix, we present an intuitive explanation of these basic morphological operations using diagrams. For a mathematically rigorous treatment, the reader is referred to Matheron (1975), Serra (1983), Dougherty & Lotufo (2003), and Chiu et al. (2013).

Suppose that a binary field (image), X, is probed by a symmetric structuring element, S. There are four elemental operations as follows.

The Minkowski addition (dilation),¹⁸ denoted by a symbol ⊕, is defined as the union of a binary field, X, and a structuring element, S, as the centre of the structuring element is moved inside the binary field. Fig. A1 (top) shows how this operation works. The Minkowski addition enlarges the original binary field with a structuring element.

The Minkowski subtraction (erosion), denoted by a symbol ⊖, is a dual to the Minkowski addition. It is defined as the intersection of a binary field, X, with the centre of a structuring element, S, as the structuring element moves inside the binary field. Fig. A1 (bottom) shows how the Minkowski subtraction works. Only the parts that can fit the structuring element inside the binary field remain. The Minkowski subtraction shrinks the original binary field with a structuring element.

The morphological opening, denoted by a symbol ∘, is then defined by a consecutive operation: Minkowski subtraction followed by Minkowski addition. The morphological opening of a binary field, X, by a structuring element, S, is expressed as X ∘ S ≡ (X ⊖ S) ⊕ S. Fig. A2 (top) shows the initial and final results of the morphological opening operation. A step-by-step algebra of shapes defined in this framework of mathematical morphology is shown in Fig. A2 (bottom).

The morphological closing is a dual operation to the morphological opening. Although here we do not use this operation, we define it for completeness. The morphological closing, denoted by a symbol •, is defined by a consecutive operation: Minkowski addition followed by Minkowski subtraction, X • S ≡ (X ⊕ S) ⊖ S.