Power spectrum multipoles and clustering wedges during the epoch of reionization

ABSTRACT

We study the viability of using power spectrum clustering wedges as summary statistics of 21 cm surveys during the epoch of reionization (EoR). For observations in a wide redshift range |$z\sim 7{\!-\!}9$| corresponding to a line-of-sight scale of |$\sim 500\,$| Mpc, the power spectrum is subject to anisotropic effects due to the evolution along the light of sight. Information on the physics of reionization can be extracted from the anisotropy using the power spectrum multipoles. Signals of the power spectrum monopole are highly correlated at scales smaller than the typical ionization bubble, which can be disentangled by including higher-order multipoles. By simulating observations of the low-frequency part of the Square Kilometre Array (SKA) Observatory, we find that the sampling of the cylindrical wavenumber |$\boldsymbol {k}$|-space is highly non-uniform due to the baseline distribution, i.e. the distribution of antenna pairs sampling different transverse |$\boldsymbol {k}_\perp$| scales. Measurements in clustering wedges partition the cylindrical |$\boldsymbol {k}$|-space into different radial |$k_\parallel$| scales, and can be used for isolating parts of |$\boldsymbol {k}$|-space with relatively uniform sampling, allowing for more precise parameter inference. Using Fisher Matrix forecasts, we find that the reionization model can be inferred with per cent level precision with |$\sim 120$| h of integration time using SKA-Low. Compared to model inference using only the power spectrum monopole above the foreground wedge, model inference using multipole power spectra in clustering wedges yields a factor of |$\sim 3$| improvement on the marginalized 1D parameter constraints.

1 INTRODUCTION

The formation of luminous objects in the early Universe is key to the evolutionary history of the cosmic large-scale structure. As the new generation of telescopes such as JWST (Gardner et al. 2023) reaches unprecedented sensitivity, for the first time a large number of very early galaxies at |$8 \lesssim z \lesssim 17$| has been observed (e.g. Castellano et al. 2022; Adams et al. 2023; Austin et al. 2023; Finkelstein et al. 2023; Harikane et al. 2023; Santini et al. 2023). The detection of a large population of galaxies hints towards tensions between the observations and the structure formation model in standard |$\Lambda$| cold dark matter (⁠|$\Lambda$|CDM) cosmology (e.g. Yung et al. 2024). However, resolving individual galaxies at very high redshifts requires observations with extreme depths, limiting the survey area of these observations and consequently the amount of information available for inferring cosmology. To compensate this fallback, complementary probes of the early Universe are needed to cover large volumes for cosmological studies of the early Universe.

A promising candidate for probing large cosmological volumes in the early Universe is the cosmic 21 cm line. Neutral hydrogen (H i) has an emission line at the rest wavelength of around 21 cm (Hellwig et al. 1970). H i is the most abundant element in the Universe after recombination (e.g. Wagoner, Fowler & Hoyle 1967), making the 21 cm line an ideal tracer of the very early Universe (e.g. Cooray 2006). After the formation of the first luminous galaxies, H i within the intergalactic medium (IGM) is ionized by the ultraviolet photons emitted by these objects during a period known as the epoch of reionization (EoR; see Furlanetto, Oh & Briggs 2006 for a review). The redshifted 21 cm line can be used to map the distribution of neutral hydrogen in the IGM, and therefore directly probes the EoR (Gnedin & Ostriker 1997).

The 21 cm line signal lies in the radio waveband and can be measured by radio telescopes. However, it is intrinsically faint and requires very sensitive telescopes with large effective collecting area to detect it. The requirements for large survey area and high-angular resolution naturally call for using radio interferometers for 21 cm observations (Madau, Meiksin & Rees 1997). Current experiments targeting the EoR include the Murchison Widefield Array (MWA; Trott et al. 2020), the Hydrogen Epoch of Reionization Array (HERA Collaboration 2023), and the Low-Frequency Array (Mertens et al. 2020). In the future, the low-frequency array of the SKA Observatory (SKAO; Raghunathan et al. 2023), the SKA-Low array, will have the sensitivity to produce high-fidelity tomographic image data for constraining the EoR (Mellema et al. 2013). Measurements of the 21 cm line signal are prone to observational effects. Foreground emissions, i.e. continuum emissions from Galactic and extragalactic sources, are several orders of magnitude brighter than the H i emission (see e.g. Chapman et al. 2012). In order to mitigate the contamination from foregrounds, statistical techniques utilizing the discrepancy in the spectral smoothness of foregrounds and the H i signal are extensively studied. For interferometric observations, the technique of delay transform (Morales & Hewitt 2004), which is Fourier transforming the measured visibility data along the frequency axis into delay space, is widely adopted (e.g. Yoshiura et al. 2021; HERA Collaboration 2023). The continuum foregrounds concentrate in low delay, leaving the high delay modes dominated by the H i signal. Measuring the H i signal by avoiding the low delay modes is referred to as foreground avoidance (Morales et al. 2012). On the other hand, blind signal separation techniques can be used to remove the smooth foregrounds, usually by quantifying the frequency–frequency covariance of the observed signal (e.g. Mertens et al. 2020; Wolz et al. 2022; Cunnington et al. 2023).

Both foreground avoidance and foreground removal techniques rely on the spectral smoothness of the foregrounds. Therefore, systematic effects that break the spectral smoothness of foregrounds can severely contaminate 21 cm observations. These effects include calibration errors (e.g. Barry et al. 2016; Byrne et al. 2019), radio frequency interference (RFI; e.g. Wilensky et al. 2019), antenna cross-coupling (Kern et al. 2020), beam side-lobes (Thyagarajan et al. 2015), and more. As a result, high-fidelity imaging of the IGM is prone to systematics, and current measurements of the EoR focus on obtaining the upper limits of the 21 cm monopole power spectrum (see fig. 30 of HERA Collaboration 2023 and references therein).

As the detection of the 21 cm signal closely intertwines with the mitigation of observational systematics, it is important to choose robust summary statistics of the 21 cm signal to infer cosmology. The spherically averaged power spectrum monopole is the most studied summary statistic for EoR, from theory (e.g. Georgiev et al. 2022), detection methodology (e.g. Liu & Shaw 2020), and model inference (e.g. Greig & Mesinger 2015). However, the power spectrum monopole does not capture the full evolutionary information of the IGM. The growth of ionization bubbles in the IGM mimics a phase transition, where the ionized regions grow and become more interconnected with each other (Furlanetto & Oh 2016; Chen et al. 2019). The distribution of ionized regions is therefore highly non-Gaussian. Parameter inferences from the 1D power spectrum may exhibit biases and degeneracies between the astrophysical parameters (Park et al. 2019). Furthermore, the 21 cm signal is observed along a 3D lightcone, and multiple observational effects introduce anisotropy into the 21 cm signal such as redshift space distortions (Mao et al. 2012) and evolution along the line of sight (Datta et al. 2012). The anisotropy may bias parameter estimation (Greig & Mesinger 2018), which cannot be fully quantified using the spherically averaged power spectrum.

In this work, we aim to provide alternative summary statistics to the spherically averaged monopole power spectrum for EoR. One common approach in the literature is to use higher-order correlation functions for parameter inference, such as the bispectrum (e.g. Shimabukuro et al. 2017; Watkinson, Greig & Mesinger 2022) and the trispectrum (e.g. Cooray, Li & Melchiorri 2008). The modelling of the covariance of higher-order correlation functions requires sophisticated simulations and good understanding of the underlying model (see e.g. Philcox & Ereza 2024 for discussions in the context of galaxy surveys). On the other hand, it has been suggested that tomographic images can be used for alternative summary statistics such as bubble size distributions (e.g. Giri et al. 2018a). Topological analysis on the image data can also be used for fully extracting the underlying information, for example from Minkowski Functionals (e.g. Yoshiura et al. 2017), shape finders (e.g. Bag et al. 2019), and Betti numbers (Giri & Mellema 2021). However, these summary statistics are more susceptible to observational effects. Even simple effects such as Gaussian thermal noise complicate the morphology of the tomographic image data (e.g. Giri, Mellema & Ghara 2018b; Chen et al. 2019). In general, topological quantities cannot be linearly formulated with regard to the data vector, making it difficult to reconstruct the signal and accurately model the covariance of the measurement.

In this paper, we argue that the power spectrum multipoles, measured in clustering wedges of the cylindrical |$(\boldsymbol {k}_\perp ,k_\parallel)$|-space, can be used as a robust summary statistic for the EoR. The power spectrum multipole is the 1D average of the cylindrical power spectrum, weighted by the Legendre polynomials of the cosine of the angle to the line-of-sight |$\mu = k_\parallel /|\boldsymbol {k}|$|⁠. It is the simplest extension to the 1D power spectrum monopole, and uses only the two-point statistics of the 21 cm signal. Therefore, techniques such as foreground avoidance and quadratic estimators can be applied to the multipoles in the same way as the monopole, making the power spectrum multipole immediately applicable to the data product of current and near-future surveys. The power spectrum in 2D |$(k_\perp ,k_\parallel)$|-space is found to be the most robust summary statistics for the 21 cm signal (Prelogović & Mesinger 2024), which directly relates to the multipoles. As power spectrum multipoles measure the anisotropy of the 21 cm signal (Majumdar et al. 2016), we expect that they may resolve the problem of parameter degeneracy, and improve the constraining power of 21 cm observations.

Furthermore, the |$\boldsymbol {k}$|-space for the weighted average of multipoles can be divided into multiple regions known as the clustering wedges (e.g. Grieb et al. 2017), based on the values of |$\mu$|⁠. By partitioning the |$\boldsymbol {k}$|-space into clustering wedges, the multipole power spectra can be used to disentangle the uneven sampling of |$\boldsymbol {k}$|-space, originated from the baseline distribution of radio interferometers. We explore the viability of using future SKA-Low observations to measure the 21 cm multipoles, and quantify the constraining power of the power spectrum clustering wedges on the EoR (see Cunnington et al. 2020; Soares et al. 2021; Berti, Spinelli & Viel 2023; Pourtsidou 2023 for post-reionization studies using the power spectrum multipoles).

The rest of this paper is organized as follows. In Section 2, we outline the procedure for simulating mock observations of SKA-Low. The definitions of the power spectrum multipoles and clustering wedges are given in Section 3. We discuss the advantages of using the power spectrum multipoles and clustering wedges comparing to the spherically averaged monopole in Section 4. In Section 5, we present Fisher Matrix formalism for constraints on the EoR model parameters from multipole measurements. The results of the forecast are presented in Section 6. In Section 7, we further clarify and discuss the caveats of the simulations used. We summarize and conclude our findings in Section 8. Throughout this paper, we assume the flat |$\Lambda$|CDM cosmology reported in Planck Collaboration VI (2020).

2 SIMULATION

2.1 The 21 cm signal

In this section, we describe the simulation of the 21 cm signal used in this paper. The 21 cm signal during the EoR is represented by the brightness temperature fluctuation against the background cosmic microwave background (CMB) temperature

where |$T_{S}$| is the gas spin temperature, |$T_\gamma$| is the CMB temperature, z is the redshift of the signal, and |$\tau _{\mu _0}$| is the optical depth at the rest wavelength of the 21 cm signal |$\mu _0$|⁠. In the optically thin approximation |$\tau _{\mu _0}\ll 1$|⁠, the brightness temperature can be written as (e.g. Furlanetto et al. 2006)

where |$x_{\rm H{\small I}}$| is the neutral fraction of the IGM, |$\delta _m$| is the matter overdensity, |$\Omega _m$| is the fraction of matter density in total energy density at present day, |$\Omega _b$| is the fraction of baryon density in total energy density at present day, H is the Hubble expansion rate, |$h = H_0 /(100\, {\rm km\,s^{-1}\,Mpc^{-1}})$| where |$H_0$| is the Hubble parameter at |$z=0$|⁠, and |${\rm d}v_\parallel /{\rm d}r$| is the gradient of the peculiar velocity field along the line of sight.

We use the publicly available software 21cmfast¹ (Mesinger, Furlanetto & Cen 2011; Murray et al. 2020) to simulate the 21 cm signal. 21cmfast uses the excursion set formalism (Furlanetto, Zaldarriaga & Hernquist 2004) to enable fast, seminumerical computations of the reionization process. In particular, it computes the neutral fraction of the IGM by iteratively checking for the maximum radius R around a cell |$\boldsymbol {x}$| where the criterion (e.g. Greig & Mesinger 2015)

can be met. Here, |$\zeta$| is the ionization efficiency parameter to be discussed later, and |$f_{\rm coll}$| is the collapse function describing the fraction of collapsed matter within a spherical radius R residing in haloes with mass larger than |$\bar{M}_{\rm min}$| (e.g. Press & Schechter 1974; Sheth & Tormen 1999). If the above criterion is met, the cell is fully ionized. If not, the ionized fraction of the cell is set to be |$\zeta \, f_{\rm coll}(\boldsymbol {x},z,R_{\rm cell},\bar{M}_{\rm min})$|⁠.

In our simulation settings, the collapse function is tuned by the minimum virial temperature of the star-forming haloes, |$T_{\rm vir}$| (see e.g. Barkana & Loeb 2001). Higher |$T_{\rm vir}$| results in fewer haloes to produce ionizing photons, requiring a larger ionization efficiency |$\zeta$| to ionize the IGM.

The ionization efficiency parameter |$\zeta$| is an accumulation of multiple physical properties of the ionizing sources

where |$f_{\rm esc}$| is the fraction of ionizing photons that escape into the IGM, |$f_*$| is the fraction of galactic gas in stars, |$N_\gamma$| is the number of ionizing photons produced per baryon in stars, and |$\bar{n}_{\rm rec}(\boldsymbol {x},z,R)$| is the average number of times an H i atom recombines. In our simulation settings, we allow inhomogeneous recombination and calculate the recombination number following the procedures in Sobacchi & Mesinger (2014), while keeping the rest of the ionization efficiency |$\xi$| fixed.

As we focus on power spectrum multipoles, the anisotropic effects need to be included in the simulation accurately. The evolution along the observational lightcone produces sizeable anisotropy (e.g. Datta et al. 2012, 2014). Furthermore, redshift space distortions due to the peculiar velocity field also have a major impact on the 21 cm signal (Mao et al. 2012; Majumdar, Bharadwaj & Choudhury 2013). Therefore, we follow Jensen et al. (2013) and simulate the full brightness temperature fluctuation without the optically thin approximation for a 3D lightcone. The simulation voxel is split into subvoxels to apply redshift space distortion effects according to the simulated velocity field, and the perturbed subvoxels are regridded to the simulation box.

Throughout this paper, the simulations have coeval boxes with |$(1600\, {\rm Mpc})^3$| size and |$(2\, {\rm Mpc})^3$| resolution. The |$1600\, {\rm Mpc}$| size corresponds to |$\sim 10\, {\rm deg}$| for |$z\sim 8$|⁠, which is at |$\sim 3\sigma$| of the power beam of SKA-Low discussed later in Section 2.2. The |$2\, {\rm Mpc}$| resolution ensures that the multipoles can be accurately modelled down to |$k\sim 1\, {\rm Mpc^{-1}}$|⁠.

An important aspect of our work is to illustrate the constraining power of the power spectrum multipoles for EoR parameters. Therefore, we choose two different sets of EoR parameters |$({\rm log_{10}[{\it T}_{vir}/K]}, \xi)$|⁠, which produce an almost identical global ionization history as well as a similar spherically averaged 1D power spectrum for comparison. Specifically, we choose |$({\rm log_{10}[{\it T}_{vir}/{\it K}]}, \xi) = (4.7,65),\, (5.1,150)$| to produce two sets of reionization lightcones. This is similar to the ‘faint galaxies’ and ‘bright galaxies’ models presented in Greig & Mesinger (2018). For simplicity, we use the same names for the two models we are considering. Lower values of |$({\rm log_{10}[T_{vir}/K]}, \xi)$| mean that the reionization process is driven by fainter, more numerous sources, and vice versa. An illustration of the differences in the reionization process between the two models is shown in Fig. 1. As seen in Fig. 1, the brightness temperature field indeed has distinct morphologies between the faint and bright models, where the reionization process in the faint model is more homogeneous. The evolution of the global ionization fraction and the 1D monopole power spectrum along the observation lightcone are presented in Fig. 2 for reference.²

We note that, in our simulation, reionization is almost complete at |$z\sim 7.0$|⁠, which is faster than what has been suggested by current measurements (e.g. Bosman et al. 2022; Greig et al. 2024b). As discussed later in Section 2.2, we focus on the redshift range of |$z\sim 7.1{\text {--}}8.8$| and thus a faster reionization produces larger anisotropic effects along the lightcone, which is ideal for presenting the case of comparing the constraining power of the monopole and the multipoles.

2.2 The interferometric observation

In this section, we present the simulation of visibility data using the configurations of the SKA-Low array. The SKA-Low array will be located at the Murchison Radio-astronomy Observatory³ (MRO), sharing the location with its precursor, the MWA (Tingay et al. 2013). SKA-Low is capable of measuring the 21 cm signal from 50 to 350 MHz, covering |$z\sim 3{\!-\!}25$|⁠. The radio environment at the MRO provides a relatively clean frequency range of 140–200 MHz for EoR measurements (Offringa et al. 2015; Sokolowski, Wayth & Lewis 2015). We choose the subband of 145–175 MHz, corresponding to the redshift range of |$z\sim 7.1{\!-\!}8.8$| towards the later stage of the EoR when more than half of the IGM is ionized (Davies et al. 2018). The channel bandwidth is to set to 200 kHz, corresponding to the line-of-sight resolution of |$k_\parallel \lesssim 1\, {\rm Mpc}^{-1}$|⁠. While the frequency resolution of SKA-Low will be finer than 200 kHz, the non-linear scales |$k_\parallel \gt 1\, {\rm Mpc}^{-1}$| are hard to model and beyond the scope of this work. The 21 cm lightcone is smoothed along the line of sight to match the frequency resolution.

We use the oskar package⁴ (Mort et al. 2010) to simulate the interferometric observations, with the telescope configuration following the latest official SKAO data challenge.⁵ The pointing centre is set to be near the EoR0 field (Lynch et al. 2021) at (RA = 0 deg, Dec. = −30 deg). The u–v coverage of the baseline distribution is simulated with a tracking observation of 12 h with a time resolution of 10 s. The baseline distribution and the power beam are shown in Fig. 3 for reference.

The standard deviation of the Gaussian thermal noise of visibilities follows the radiometer equation (e.g. Condon & Ransom 2016):

where |$T_{\rm sys}$| is the system temperature, |$A_{\rm e}$| is the effective collecting area of the antenna, |$\delta t = 10\,$| s is the time resolution and |$\delta f = 200$| kHz is the channel bandwidth. In this paper, we aim to provide forecasts for power spectrum multipole measurements of a deep observation with coherent averaging of multiple nights. For a total integration time of |$t_{\rm tot}$|⁠, the thermal noise amplitude for each baseline is scaled so that

where |$t_{\rm n} = 12\,$| h is the observation time of one night. In this paper, we follow the anticipated performance of SKA-Low in Braun et al. (2019) to produce the thermal noise. The natural sensitivity of the instrument is set to be |$A_{\rm e}/T_{\rm sys} = [1.052,1.119,1.159, 1.171,1.190]\, {\rm m^2\,K^{-1}}$| at [140,150,160,170,180] MHz, and the amplitudes at each frequency channels are linearly interpolated. The total observation time is chosen to be |$t_{\rm tot} = 120\,$| h. The choice of integration time of 120 h is relatively conservative, as the EoR0 field will potentially be observed for |$\gt 1000$| h. We note that, a large amount of data loss is expected for future surveys due to RFI flagging. It is also common in data analysis to cross-correlate different time blocks to reduce systematic effects (e.g. Abdurashidova et al. 2022; Paul et al. 2023), which effective reduces the integration time by half. Overall, we expect |$t_{\rm tot} = 120\,$| h to be a reasonable lower limit for future SKA-Low observations of one deep field.

As a preliminary study into the detectability and constraining power of EoR multipoles, simulating observational systematics is beyond the scope of our work. In particular, we assume that foreground contamination can be avoided by applying a horizon criterion so that

where the position of the wedge |$c_k$| is determined by a characteristic angular scale |$\theta _0$| so that (Liu, Parsons & Trott 2014)

where |$D_c(z)$| is the comoving distance. We choose |$\theta _0 = \sqrt{\Omega _{\rm beam}}$| following Chen et al. (2023b) which gives |$c_k = 0.28$|⁠, where |$\Omega _{\rm beam}$| is the beam area calculated using the power beam shown in Fig. 3. As the 21 cm signal during the EoR is anisotropic, the power spectrum averaged above the foreground wedge is biased compared to the spherical average (Jensen et al. 2016). It also leads to bias in the multipoles, inducing effective mode-mixing in the 1D multipole power spectra (Raut, Choudhury & Ghara 2018). In this paper, we refer to the average above the foreground wedge as ‘spherical average’ for simplicity. The partition into clustering wedges discussed in Section 4.3 is performed above the foreground wedge instead of using the entire |$\boldsymbol {k}$|-space. Equation (7) can also be written as

where |$\mu \equiv k_\parallel /|\boldsymbol {k}|$| is used for the multipole expansion discussed later in Section 3.

Several observational effects induce foreground leakage into the observation window above the wedge. For example, chromatic data excision due to narrow-band RFI will cause foreground spillover (Wilensky, Hazelton & Morales 2022). Calibration errors also introduce foreground scatter into high delay (e.g. Barry et al. 2016; Byrne et al. 2019). Although not discussed in this paper, we note that the systematic effects can be partially mitigated in the quadratic estimator formalism, for example by including the foreground cleaning operation in the weighting matrix |$\mathbf {R}$| (e.g. Kern & Liu 2021; Chen, Wolz & Battye 2023a; Chen et al. 2023b).

3 POWER SPECTRUM MULTIPOLES IN CLUSTERING WEDGES

In this section, we give an introduction to the power spectrum multipoles in clustering wedges as summary statistics of the EoR. In particular, we discuss the relation between the multipoles and the anisotropy of the 3D power spectrum, and the observational effects in the measured multipole power spectra. These effects lay the foundations for the incentives for using multipoles in wedges, which we discuss later in Section 4.

3.1 Definitions

In the flat-sky and plane-parallel limit, the power spectrum multipole of order |$\ell$| in a clustering wedge can be written as (e.g. Scoccimarro 2015)

where |$k =|\boldsymbol {k}|$| is the mode of the 3D wavenumber vector, |$\mu = k_\parallel /k$|⁠, |$\mu _0$| and |$\mu _1$| are the lower and upper limits of the clustering wedge, |$\mathcal {P}_\ell$| is the Legendre polynomial, |$P_{21}(\boldsymbol {k}) \equiv P_{21}(k,\mu)$| is the anisotropic 21 cm power spectrum and |$w(\boldsymbol {k})$| is the weighting of the 3D power spectrum. The weights are renormalized so that for a specific 1D k-bin, |$\int {\rm d}\mu \, w(k,\mu) = 1$|⁠. We choose the spacing of the k-bins to be logarithmically distributed from |$0.05\, {\rm Mpc^{-1}}$| to |$1\, {\rm Mpc^{-1}}$| with 20 bins. We consider |$\ell = 0,2,4$|⁠, i.e. the monopole, quadrupole, and hexadecapole, in this paper.

For illustration, we also use the ‘dimensionless’ power spectrum:⁶

From equation (11), we can introduce the terminology of the summary statistics. A ‘wedge’ is defined as a region in the |$\boldsymbol {k}$|-space enclosed by the integration limit |$\mu _0 \lt \mu \lt \mu _1$|⁠. For example, the region dominated by the foregrounds, as defined in equation (8), is referred to as the foreground wedge and corresponds to |$\mu _0 = 0$|⁠, |$\mu _1 = \mu _{\rm fg}$|⁠. Consequently, a ‘clustering wedge’ is a region in the |$\boldsymbol {k}$|-space outside the foreground wedge with an arbitrary choice of |$\mu _0$| and |$\mu _1$|⁠, as long as |$\mu _0 \gt \mu _{\rm fg}$|⁠. Throughout this paper, we use ‘spherical average’ to refer to the case where there is no further splitting of |$\boldsymbol {k}$|-space above the foreground wedge. On the contrary, we use ‘monopole/multipoles in clustering wedges’ to refer to the case where the region above the foreground wedge is split into multiple wedges.

3.2 Multipoles and anisotropy

Using the formalism described by equation (11), we first discuss the relation between the power spectrum anisotropy and the multipoles, without the partitioning of clustering wedges and observational effects. In this case, the integration range of |$\mu$| in equation (11) is simply |$\mu _1 = 1$| and |$\mu _0 = 0$|⁠. The power spectrum is uniformly weighted so that |$w(\boldsymbol {k}) = 1$|⁠. We note that the foreground wedge is not included in this subsection and will be discussed later in Section 3.3.

From equation (11), it can be seen that multipoles of different order |$\ell$| are essentially 1D averages of the 3D power spectrum with different weighting schemes, and the weighting schemes are described by the Legendre polynomial |$\mathcal {P}_{\ell }(\mu)$|⁠. As a result, the weighted averages of the power spectrum can be used as a probe for the anisotropy, as illustrated in Fig. 4. For both faint and bright models, it can be seen that the signal power spectrum is highly anisotropic, demonstrated by the contours of constant power spectrum amplitude in the left panels of Fig. 4. The amplitude of the 21 cm power spectrum is visibly higher around |$k_\parallel \sim 1\, {\rm Mpc^{-1}}, k_\perp \sim 0$|⁠. This is due to the fact that the power spectrum is measured along a large lightcone shown in Fig. 1. Along the line-of-sight direction, there is a large evolution of the global ionization fraction of the IGM, which corresponds to an enhanced fluctuation of at large transverse, small radial scales. Meanwhile, the power spectrum amplitude is also higher around |$k_\parallel \sim 0, k_\perp \lesssim 0.5\, {\rm Mpc^{-1}}$|⁠. At these scales, the power spectrum probes the distribution of ionization sources averaged along the line of sight. Overall, for the same |$|\boldsymbol {k}|$|⁠, the power spectrum has higher amplitude at smaller |$\mu$|⁠.

The morphology of ionization bubbles and their evolution lead to the anisotropy in the power spectrum. As a result, the anisotropic features in the power spectrum differ for the faint and bright models, as seen in Fig. 4. For the bright model, the amplitude of the power spectrum is significantly lower at |$k_\parallel \sim 0, k_\perp \lesssim 0.5\, {\rm Mpc^{-1}}$| comparing to the faint model. As shown in Fig. 1, the bright model gives a more rapid reionization process, with larger but fewer bubbles. Therefore, the fluctuation of the brightness temperature is smaller at small transverse scales for the bright model.

The differences in the anisotropy are then captured in the power spectrum multipoles. In Fig. 5, we show the averaged 1D multipole power spectrum of the 21 cm signal. The monopole power spectrum gives constant weights |$\mathcal {P}_{\ell =0}(\mu)=1$| to the 3D power spectrum before performing the 1D averaging. Therefore, the monopole power spectrum only gives the overall amplitude averaged along a certain scale |$|\boldsymbol {k}|$|⁠. On the other hand, the quadrupole gives positive weights to large |$\mu$| and negative weights to small |$\mu$|⁠. As the 21 cm power spectrum is lower at higher |$\mu$|⁠, the overall quadrupole is negative except for large scales |$k\lesssim 0.05\, {\rm Mpc^{-1}}$|⁠. The hexadecapole weight |$\mathcal {P}_{\ell =4}(\mu)$| is positive at |$\mu \sim 0$|⁠. It first decreases as |$\mu$| increases and the increases, leading to negative values at |$\mu \sim 0.5$| and positive values at |$\mu \sim 1$|⁠. The averaging results in a positive amplitude in the 1D average. Since the power spectrum is less anisotropic for the bright model, the amplitudes of the multipoles are lower comparing to the faint model. This suggests that multipoles can be used to extract information on the reionization model, as we discuss more in Section 4.1.

3.3 Observational effects in multipoles

The 21 cm multipoles reflect the underlying anisotropy of the 3D power spectrum. However, the multipoles are also affected by observational effects. The observational effects influence the summary statistics in different ways, which we outline in this subsection:

• There is mode-mixing along the radial and the transverse direction. Along the line-of-sight, different |$k_\parallel$|-modes are mixed due to bandpass tapering. |$\boldsymbol {k}_\perp$|-modes are also mixed due to the attenuation of the primary beam. In this work, we assume mode-mixing is negligible by using relatively large u-v grids and a renormalization matrix in the quadratic estimator for power spectrum estimation. See Appendix  A for more discussions.

• The foreground wedge forbids measurements in |$\mu \lt \mu _{\rm fg}$|⁠. In the ideal case, the entire |$\boldsymbol {k}$|-space from |$\mu =0$| to |$\mu =1$| is sampled. For spherically averaged multipoles, the measured spectra integrated from |$\mu _0 = \mu _{\rm fg}$| will be biased comparing to the ideal spectra integrated from |$\mu _0 = 0$|⁠.

• The |$\boldsymbol {k}$|-space is not uniformly sampled, with an additional weighting introduced by the baseline distribution. In the ideal case, the weighting of the |$\boldsymbol {k}$|-modes |$w(\boldsymbol {k})$| is uniform so that |$w(\boldsymbol {k}) = 1$|⁠. In real observations, however, the |$\boldsymbol {k}$|-modes will be weighted by the inverse of the noise covariance, so that |$w(\boldsymbol {k}) = \frac{{\rm d}N}{{\rm d}\boldsymbol {k}_\perp }$| where |$\frac{{\rm d}N}{{\rm d}\boldsymbol {k}_\perp }$| is the probability distribution of the baselines.

These observational effects can be disentangled from the underlying 21 cm signal by utilizing the multipoles in clustering wedges, which we discuss in detail in Section 4. For this section, we first establish the intuition of the importance of these effects, by showing the biasing of the power spectrum from the observational effects in Fig. 6. For simplicity, only the faint model is shown here. When foreground avoidance is included, the monopole power spectrum decreases in amplitude. This is due to the fact that the lower |$\mu$| region is excluded where the power spectrum amplitude is higher. For the quadrupole, the foreground wedge corresponds to the region where the Legendre polynomial |$\mathcal {P}_{\ell = 2}(\mu)$| is negative. Therefore, excluding the foreground wedge results in a positive-valued quadrupole power spectrum. For the hexadecapole, when foreground avoidance is not included, the averaging gives positive amplitude since it weights |$\mu \sim 0$| positively. Including the foreground avoidance gives negative amplitude, since |$\mathcal {P}_{\ell = 4}(\mu)$| is negative for |$\mu \sim 0.5$| where the 3D power spectrum has higher amplitude, while positive for |$\mu \sim 1.0$| where the 3D power spectrum is lower.

The uneven |$\boldsymbol {k}_\perp$| sampling further biases the power spectrum multipoles. The exact effects on the multipoles depend on the baseline distribution and |$|\boldsymbol {k}|$|-bins. For SKA-Low, we find that in general, higher values of |$\mu$| are sampled more frequently, which we discuss further in Section 4.3. The monopole power spectrum is therefore lowered as its amplitude is smaller at large |$\mu$|⁠. On the other hand, the higher-order Legendre polynomials are positive at |$\mu \sim 1$|⁠. Therefore, the quadrupole and the hexadecapole become larger when baseline distribution is considered.

In conclusion, the observational effects have significant impact on the measured power spectrum multipoles of the 21 cm signal. In turn, measuring the |$(\boldsymbol {k}_\perp ,k_\parallel)$| dependence of the power spectrum amplitude helps disentangle observational effects with the underlying signal. This provides strong incentives to use power spectrum multipoles in clustering wedges, which we demonstrate in the next section.

4 ADVANTAGES OF MULTIPOLES AND CLUSTERING WEDGES

In this section, we demonstrate the advantages of using multipoles in clustering wedges as summary statistics for EoR surveys. Using power spectrum multipoles allows us to extract information from anisotropy and to probe the evolution of the EoR along the lightcone. Furthermore, the advantage of partitioning the measurements into clustering wedges is illustrated. For simplicity, when not specifically mentioned, the results shown correspond to the simulation using the faint model defined in Section 2.1. The results shown are averages over 10 independent realizations. From here on, all power spectra are calculated with the inverse noise covariance weighted, baseline-sampled 3D bandpowers using the quadratic estimator formalism described in Appendix  A.

In order to understand the information content of the multipoles, we also need to compute the covariance of the power spectra. While the noise covariance can be analytically calculated in terms of its two-point correlation function (see Appendix  A), the signal itself is non-Gaussian and requires the computation of high-order correlation functions (Shaw, Bharadwaj & Mondal 2020) to calculate its covariance. Instead of the quadratic estimator, we use the jackknife method and calculate the signal covariance directly from the simulation lightcone. The lightcone is resampled 25 times along the transverse direction, each with a subarea of |$(320\, {\rm Mpc})^2$| taken out. The conversion from the delay power spectrum to the 21 cm temperature power spectrum is susceptible to the lightcone effects, the plane-parallel approximation, and the treatment of the primary beam. For the results shown in this paper, we use the measured power spectra from visibility data, and ignore the systematic biases of the estimation. We note that, simulation-based inference (e.g. Zhao et al. 2022; Saxena et al. 2023; Greig et al. 2024a) can be used to circumvent these problems, as it does not require an explicit likelihood and takes into account systematic effects through forward-modelling. The detailed treatment of an inference framework to be used on real data is beyond the scope of our work.

4.1 Information on anisotropy

The fluctuation of the 21 cm signal, as described by equation (2), is determined by the multiplication of the ionization field, the spin temperature, the matter density field and the velocity field. While quantities such as the velocity field (e.g. Chapman & Santos 2019) and the spin temperature (e.g. Schaeffer, Giri & Schneider 2024) have a sizeable impact on the 21 cm power spectrum, the features of the fluctuation are mostly influenced by the matter density on large cosmological scales and the ionization field on small scales (Georgiev et al. 2022). Particularly, after the beginning of the percolation of the ionization bubbles, the ionization field becomes highly non-Gaussian (e.g. Shimabukuro et al. 2015) and non-linear (e.g. Hassan et al. 2016). As a result, the evolution of ionization bubbles along the line-of-sight induces anisotropy into the 21 cm power spectrum. The anisotropy can be captured in the cylindrical power spectrum. However, observational effects such as beam attenuation and frequency-dependent sampling create mode-mixing of different |$\boldsymbol {k}$|-modes, making it difficult to accurately measure the signal and its covariance in fine cylindrical |$\boldsymbol {k}$|-bins. Alternatively, power spectrum multipoles can be used for measuring the anisotropy.

For the spherically averaged power spectrum monopole, anisotropic information is lost, leading to obscurity in the interpretation of the amplitude of the signal. The lack of distinguishing power in the 1D monopole is illustrated in the upper panel of Fig. 7. For the cases of increasing |$\xi$| and decreasing |$\log T_{\rm vir}$|⁠, the changes in the power spectrum amplitude are similar. This leads to a potential degeneracy between the model parameters, and indicates the limitations in the constraining power from monopole. On the other hand, as shown in the lower panel of Fig. 7, the differences in the cylindrical power spectrum contain information on the model parameters. The constraining power is then reflected in the higher-order multipoles |$\ell = 2$| and |$\ell = 4$|⁠. Increasing |$\xi$| and decreasing |$\log T_{\rm vir}$| lead to changes in the multipole power spectra with different structures. The differences are visible at all scales in the quadrupole, which can also be seen in the hexadecapole at small scales |$k\gtrsim 0.3\, {\rm Mpc^{-1}}$|⁠.

4.2 Decorrelation of small-scale measurements

The patchiness of ionization bubbles during the EoR breaks the tracing of the brightness temperature field to the density field. As a result, at small scales the temperature power spectrum is highly correlated towards the end of reionization (e.g. Mondal, Bharadwaj & Majumdar 2017). This is due to the fact that the signal is dominated by correlation within the same ionized structure at scales smaller than the typical size of the bubbles.

To illustrate this correlation, we show the correlation matrix of the monopole of the pure 21 cm signals in Fig. 8 for both the faint and bright models. For the faint model, the power spectrum signal is highly correlated at relatively small scales |$k\gtrsim 0.2\, {\rm Mpc^{-1}}$|⁠. The monopole power spectrum is even more correlated for the bright model, especially on the large scales. The small-scale correlation indicates that the fluctuations of the ionization field are dominating the signal. To further establish the connection between the structures of the bubbles and the power spectrum correlation, we show the bubble size distribution for both faint and bright models in the right panel of Fig. 8. The bubble size distribution is calculated using the mean-free-path method (Mesinger & Furlanetto 2007), and only cells that are completely ionized are counted into the bubble sizes. Due to the fact that the reionization process is driven by brighter and less numerous sources in the bright model, there are more large bubbles in the bright model than in the faint model, which can be seen in Fig. 8. The increase in bubble sizes corresponds to the stronger correlation at large scales, verifying our conclusion on the connection between bubble sizes and signal correlation.

As the signal becomes completely correlated at small scales |$k\sim 0.3\, {\rm Mpc^{-1}}$|⁠, the constraining power of the monopole on EoR parameters is severely limited. Although SKA-Low will be able to probe a wide range of scales, the extraction of information from measurements of the monopole is not sufficient. For forthcoming surveys, this issue may be compounded by the fact that shorter baselines are more susceptible to observational systematics (e.g. Amiri et al. 2023; Paul et al. 2023), leading to further signal loss. Alternative summary statistics that can probe small scales are thus called for.

In order to recover information from the small scale measurements, we can utilise higher-order multipoles. As we have shown in Section 3, the multipoles probe the evolution of the signal along the line of sight, subsequently disentangling the signal at smaller scales. In Fig. 9, we show the correlation matrix for the multipole data vector of the pure 21 cm signal including |$\ell = 0,2,4$|⁠. For the quadrupole, the signal completely decorrelates at large scales |$k\lesssim 0.1\, {\rm Mpc^{-1}}$| and becomes less correlated compared to the monopole at |$k\gtrsim 0.1\, {\rm Mpc^{-1}}$|⁠. For the hexadecapole, the correlation is even weaker and the signal only becomes correlated at |$k\gtrsim 0.5\, {\rm Mpc^{-1}}$|⁠. Combined with the fact that the multipoles can distinguish changes in the reionization parameters as shown in Fig. 7, this indicates that including the quadrupole and the hexadecapole will significantly increase the information content of 21 cm surveys.

The cross-correlation between |$P^{\ell = 0}_{21}$| and |$P^{\ell = 4}_{21}$| exhibits a transition from anticorrelation to correlation, as shown in Fig. 9. This is due to the baseline distribution leading to different sampling of the cylindrical k-space at different scales, which we discuss next.

4.3 Non-uniform sampling from baseline distribution

The |$\boldsymbol {k}$|-modes measured in a 21 cm survey are determined by the distribution of the interferometric baselines. For redundant arrays like HERA, specific |$|\boldsymbol {k}_\perp |$| modes will be densely sampled and the k-dependence of the measured power spectrum is mainly contributed by the line-of-sight direction (see e.g. DeBoer et al. 2017). SKA-Low will be able to cover various scales on the transverse plane. However, the sampling of the |$\boldsymbol {k}_\perp$| space is not homogeneous, as the sensitivity of the array will concentrate around specific scales of interest. As a result, the 1D average of multipoles will have highly uneven sampling across the cylindrical |$\boldsymbol {k}$|-space.

In Fig. 10, we show the 1D baseline distribution and the corresponding distribution function of |$\mu$| for each 1D k-bin. The baseline distribution peaks around |$k_\perp \sim 0.1\, {\rm Mpc^{-1}}$|⁠. Therefore, the sampling of |$\mu$| is relatively even for large scales |$k \lesssim 0.1\, {\rm Mpc^{-1}}$|⁠. On the other hand, the baseline distribution decreases sharply after the peak. For relatively small scales |$k \gtrsim 0.1\, {\rm Mpc^{-1}}$|⁠, the contribution into the 1D power spectrum comes mainly from the large |$k_\parallel$| region, skewing the distribution of |$\mu$| towards higher values. The large deviation from uniformly sampled |$\mu$| indicates loss of information, as spherically averaged power spectra will not be able to fully probe the |$k_\parallel$|-dependence of the signal and rather measures certain regions of large |$\mu$|⁠.

To resolve the uneven sampling, we divide the |$\mu$|-region above the foreground wedge into three different clustering wedges, shown in Fig. 10. The foreground wedge is at |$c_k = 0.28$| and equivalently |$\mu = 0.267$|⁠. The |$\boldsymbol {k}$|-space above the foreground wedge is then split into three equally spaced regions of |$\mu$|⁠. The partitioning resolves the uneven sampling by isolating less-sampled regions of |$\mu$|⁠, leading to relatively even sampling in each clustering wedge. As mentioned in Section 4.2, one of the consequences of the uneven sampling is the transition from correlation to anticorrelation between the monopole and the hexadecapole as shown in Fig. 9. As a comparison, we show the correlation matrix for the multipoles of the pure 21 cm signal within each wedge in Fig. 11. The transition disappears when the clustering wedges are used. We note that, due to the non-Gaussianity of the 21 cm signal leading to higher-order correlation (e.g. Majumdar et al. 2018), on small scales the k-modes are correlated between clustering wedges. The full correlation matrix including the cross-correlation between different wedges is shown in Fig. A1 for reference.

Comparing the correlation of multipoles in wedges with the spherical average, we can see that in the wedge with the largest |$\mu$|⁠, the multipole power spectra are less correlated. This is consistent with the fact that the |$\boldsymbol {k}$|-modes with higher |$\mu$| probe small line-of-sight scales, providing more information on the evolution along the lightcone. In the spherical average, at lower k the sample is uniform. At higher k, the sampling becomes heavily skewed towards the third wedge, and the difference in the sampling of |$\mu$| at different k leads to the structure in Fig. 9.

5 PARAMETER CONSTRAINTS

5.1 Fisher matrix

In this section, we forecast the constraining power of power spectrum multipoles for SKA-Low. The information content in summary statistics can be quantified via the Fisher matrix formalism (see e.g. Euclid Collaboration et al. 2020)

where |$\rm tr$| denotes the trace of a matrix, |$\boldsymbol {\mu }$| is the ensemble average of the data vector of the summary statistics, |$\mathbf {C}$| is the data covariance, |$^{\rm T}$| denotes the transpose operation and |$\lbrace \theta \rbrace$| is the parameter set. The inverse of the Fisher matrix gives a lower bound of the error covariance of parameter inference (Rao 1992; Cramér 1999)

where |$\ge$| indicates that |$\mathbf {C}_{\theta }-\mathbf {F}^{-1}$| is positive semi-definite.

The |$1\sigma$| measurement error of a parameter |$\alpha$| can then be written as

The correlation coefficient between two parameters |$\alpha$| and |$\beta$| can be written as

In the limit where the lower bound is achieved, we can use the inverse of the Fisher matrix as the covariance matrix of the parameters which we use for the forecast.

To calculate the partial derivative of the power spectrum multipoles and the data covariance, we run the EoR simulations with small changes to the input parameters, shifting the values of the parameters around the fiducial values. We fix the initial condition when varying the EoR parameters in each realization, and average the calculated derivatives over 10 different realizations, after which we find that convergence has been reached as shown in Appendix  B. We test the choice of step size between 0.1 per cent and 5 per cent and choose it to be 1 per cent to calculate the derivatives. While we find the correlation between the parameters vary slightly, the forecasts for parameter constraints discussed in Section 5 are consistent with different step sizes.

In large-scale structure surveys, it is often assumed that either the mean of the data vector is zero or the derivatives of the data covariance is negligible (e.g. Euclid Collaboration et al. 2020). We find that in our case, both terms on the r.h.s of equation (13) have non-negligible contribution, consistent with the findings in Prelogović & Mesinger (2024).

The total data covariance matrix is the sum of signal covariance and noise covariance. As discussed in Section 4.3, the signal covariance is calculated with the baseline distribution propagated into the power spectra, and the full covariance, including cross-correlation between different wedges as shown in Fig. A1, is used to calculate the Fisher matrix. As mentioned in Section 3, we use the quadratic estimator formalism to calculate the noise covariance, which is then added to the signal covariance. The amplitude of the thermal noise fluctuation is calculated assuming 120 h of integration time following equation (5) as discussed in Section 2.2.

5.2 Detectability of multipoles

Using the total data covariance, we present the forecasts on the measurement errors of power spectrum multipoles for SKA-Low in Fig. 12. To examine which scales can be probed by SKA-Low, we calculate a qualitative upper limit shown as the black dashed lines, above which the measurements provide negligible information. The limits are calculated using the maximum and the minimum of the multipoles in different wedges shown in Fig. 12 and correspond to |$\Delta ^{0,2,4}_{\rm limit} = 10, 50, 50\, {\rm mK^2}$|⁠. For reference, the forecasts for measurement errors for the spherically averaged power spectra are also shown in Fig. 12.

For 120 h of integration time, SKA-Low will be able to measure the multipoles up to |$k\sim 0.2 \, {\rm Mpc^{-1}}$|⁠. Accurately modelling the data covariance is therefore essential, as the signal at small scales is non-Gaussian. Since the baseline distribution peaks around |$k_\perp \sim 0.1 \, {\rm Mpc^{-1}}$|⁠, sensitivities at smaller scales decrease sharply, and there is limited amount of information that can be extracted. For a complete, futuristic EoR survey, measurements down to |$k\sim 1.0 \, {\rm Mpc^{-1}}$| can be reached for |$\sim 10\,000$| h of integration time. We refrain from discussing such a scenario. As discussed later in Section 6, we find that per cent level precision can be reached for |$\sim 120\,$| h. For accurate forecasts of |$\sim 10\,000$| h, the requirements on the accurate modelling of the signal and its covariance will be very stringent, and are beyond the scope of this work.

6 RESULTS

Using the calculated data covariance and the Fisher matrix formalism, we present the forecasts for constraints on EoR parameters for SKA-Low. In order to compare using multipoles as summary statistics with other approaches, we consider four scenarios: Using only the spherically averaged monopole (‘mono + avg’), the spherically averaged multipoles (‘multi + avg’), the monopole in clustering wedges (‘mono + wedge’) and the multipoles in clustering wedges (‘multi + wedge’). Both the faint and the bright models are considered, in order to demonstrate the constraining power for different types of reionization morphology.

In Table 1, we list the forecasts for the |$1\sigma$| measurement errors of the EoR parameters. Apart from the |$1\sigma$| confidence interval, we also show the correlation coefficient of the two parameters as a reference for parameter degeneracy, as well as the Fisher information volume. For both models, we can see that the constraining power of SKA-Low increases significantly by including higher-order multipoles and using clustering wedges. The measurement errors are lower by a factor of |$\sim 2$| when either multipoles or clustering wedges are included. The improvement suggests that a significant amount of information lies in the anisotropy of the 21 cm power spectrum. Combining the power spectrum multipoles with clustering wedges, we find that the measurement errors decrease by a factor of |$\sim 3$| comparing to using spherically averaged monopole. For a survey of |$120\,$| h, the SKA-Low will be able to constrain the ionization efficiency |$\xi$| up to |$\sim 3~{{\ \rm per\ cent}}$| precision and the minimum virial temperature |${\rm log}[T_{\rm vir}]$| up to |$\lesssim 1~{{\ \rm per\ cent}}$| precision.

For the faint model, we find that the parameters are highly degenerate, as illustrated in Fig. 7. Including the multipoles yields a slight improvement in resolving the degeneracy. The lack of improvement can be explained by the lack of information on small scales. The measurements on scales |$k \gtrsim 0.1 \, {\rm Mpc^{-1}}$| are expected to probe small scales along line of sight, which provides information on the evolution of EoR. For the noise level assumed in this work, the small scales cannot be measured with high precision and thus the degeneracy persists. For the bright model, we find no significant correlation between the parameters. We visualize the |$1\sigma$| confidence region of our forecasts in Fig. 13.

The constraints on the EoR parameters are much more stringent for the bright model. This is due to the fact that the bright model adopts extreme fiducial values for the model parameters. As reionization is driven by massive and bright sources in the bright model, the variation in the model parameters produces bigger changes in the power spectra comparing to the faint model.

In conclusion, we find that for SKA-Low, a 21 cm survey with |$120\,$|h of integration time will be able to measure the 21 cm power spectrum with high sensitivity. The visibilities can be summarized into power spectrum multipoles in clustering wedges, and can be used to accurately constrain the EoR parameters.

7 DISCUSSIONS

We have shown in Section 6 that the multipoles in clustering wedges can serve as a robust summary statistics for EoR inference. It is nevertheless important to point out that the forecasts may be subject to the specific choices of simulation configurations in this paper, which we discuss in this section.

As we have discussed in this work, the improvement of parameter constraints from the power spectrum multipoles comparing to the monopole originates from the anisotropy in the power spectrum. We choose the redshift range |$z\sim 7.1{\!-\!}8.8$| towards the end of the reionization, to maximize the effects of bubble morphology and consequently the anisotropy. For the same reason, we also do not split the frequency range into multiple subbands, so that the evolution along the line of sight is maximized. While it is beyond the scope of this work, we note that for the lower frequency range |$50{\!-\!}145\,$| MHz, the power spectrum traces more closely the underlying density and the spin temperature field. As a result, we do not expect the same level of improvements in the constraints on reionization parameters from the multipoles in clustering wedges.

We also comment on the fact that the multipoles in clustering wedges can only be used when wide-field effects in EoR observations are sufficiently mitigated. Throughout this work, we have assumed the foreground wedge to be |$\mu \lt 0.27$|⁠, corresponding to the size of the instrument beam. This leaves a large part of the |$\boldsymbol {k}$|-space for 21 cm measurements, and allows us to perform multipole expansions as well as splitting the k-space into clustering wedges. If the physical horizon is used instead, i.e. the foreground leakage is significant for the whole sky beyond the beam field of view, the 21 cm power spectrum can only be measured at |$\mu \gt 0.96$|⁠. In that case, the multipole expansions and clustering wedges will not be meaningful since there is almost no variation for |$\mu$|⁠. We find that for the spherically averaged monopole, the measurement errors are enlarged by a factor of |$\sim 2.5$| when the foreground wedge is at the physical horizon, as shown in Appendix  C.

We note that the degeneracy of the parameters is discussed in the simplified setting of the EoR simulation, where only two parameters |$\xi$| and |$T_{\rm vir}$| are considered. For simulations with an extended set of model parameters, the degeneracies between parameters, and how multipoles in clustering wedges may be used to resolve them, can be more complicated. As we have shown in Appendix  B, the covariance calculation converges when averaged across 10 realizations. For a more realistic reionization model, the number of realizations needed is likely to increase. The noise covariance is calculated analytically as described in Appendix  A, leading to a total covariance matrix that is numerically stable when inverted. For a more optimistic assumption of the total observation time |$t_{\rm tot} \gg 120\,$| h, the signal covariance will be dominant and demand a large number of realizations for accurate modelling. Accurate description of the degeneracies between a large number of parameters requires Bayesian inference of the posterior instead of the Fisher matrix analysis, which we leave for future work.

Finally, we comment on the relation between the multipoles in clustering wedges and the cylindrical power spectrum. The multipoles and clustering wedges can be viewed as alternative ways of binning the |$\boldsymbol {k}$|-space, as opposed to binning in terms of |$\boldsymbol {k}_\perp$| and |$k_\parallel$|⁠. In theory, the cylindrical power spectrum should have the complete information content of the two-point statistics of the EoR signal. However, the cylindrical binning of |$\boldsymbol {k}$|-space results in a signal covariance matrix that is hard to model. The fine resolution in |$\boldsymbol {k}$|-space requires a large number of realizations for each parameter set for convergence. For example, Prelogović & Mesinger (2024) uses 130 realizations for each sample of the parameter space. Moreover, many observational effects will introduce mode-mixing to the power spectrum, such as primary beam, bandpass tapering and chromatic sampling rate. All of these effects further complicate the calculation of the signal covariance, which are beyond the scope of this work to quantify. The multipoles, on the other hand, are calculated in 1D k-bins, and therefore are less susceptible to these effects.

8 CONCLUSIONS

In this paper, we explore the advantages of using power spectrum multipoles and clustering wedges as summary statistics over using the spherically averaged monopole power spectrum for EoR measurements. In particular, we comprehensively demonstrate where additional information arises in the power spectrum multipoles and in the partition of the cylindrical |$\boldsymbol {k}$|-space. We then present forecasts for parameter constraints using multipole power spectra and clustering wedges for SKA-Low observations of redshift range |$z\sim 7.1-8.8$|⁠.

We confirm the importance of extracting information from power spectrum anisotropy. Through simulating mock 21 cm observations for future SKA-Low surveys using different model parameters for the EoR, we verify that the 1D monopole power spectrum is not optimal in distinguishing physical models of reionization, as it averages over the entire |$\boldsymbol {k}$|-space without retaining the anisotropic information of the 21 cm signal. We point out that power spectrum multipoles provide a direct summary of the anisotropic 21 cm cylindrical power spectrum. There is visible improvement in the constraining power from the amplitude of the power spectra when multipoles are included. Small-scale measurements of multipoles probe small-scale evolution along the line of sight and contain information for constraining EoR parameters.

EoR signals at small scales are dominated by the morphology of the ionization bubbles. The resulting monopole power spectrum is therefore highly correlated, especially for the latter stage of reionization considered in this work. We quantify the covariance of the signal with a large suite of simulations of the reionization lightcone, and find that small scales |$k\gtrsim 0.2\, {\rm Mpc^{-1}}$| are almost completely correlated for the monopole. We verify that this correlation is indeed caused by the growth of ionization bubbles, as simulations with larger and fewer bubbles show stronger correlation at large scales. To our knowledge, we are the first to point out that this correlation can be disentangled by including the quadrupole and hexadecapole power spectra. The higher-order multipoles probe line-of-sight scales along which reionization evolve rapidly, providing information beyond the compressed 1D fluctuation.

The baseline distribution of a radio interferometer is not uniform. As a result, we find that for SKA-Low, the sampling in the cylindrical |$\boldsymbol {k}$|-space is skewed heavily towards large values of |$\mu$| for scales |$k\gtrsim 0.2\, {\rm Mpc^{-1}}$| when averaged into 1D power spectra. Therefore, measurements at large k effectively only probe a small region of |$\boldsymbol {k}$|-space with high |$k_\parallel$|⁠. This calls for the partitioning of |$\boldsymbol {k}$|-space into clustering wedges. Within each clustering wedge, the sampling of |$\mu$| becomes relatively uniform. The power spectrum multipoles can then be used to measure the |$k_\parallel$|-dependence of the signal, optimally extracting information on the EoR. Accurate calculation of the signal covariance will be crucial for model inference, as there are significant contributions from higher-order correlations in the covariance matrix.

In order to estimate the constraining power of the multipoles and clustering wedges for the EoR, we perform Fisher matrix forecasts to quantify the error covariance of reionization parameters assuming |$120\,$| h of observation with SKA-Low. We find that on large scales |$k\lesssim 0.1\, {\rm Mpc^{-1}}$|⁠, the multipole power spectra can be measured with relatively high precision. As a result, SKA-Low will be able to give stringent constraints on EoR parameters, with the measurement errors reaching |$\lesssim 1~{{\ \rm per\ cent}}$|⁠. Comparing to only using the spherically averaged monopole, the power spectrum multipoles in clustering wedges yield a factor of |$\sim$|3 improvement on parameter constraints.

It is worth noting that the strong correlation between the EoR parameters does not disappear when multipoles and clustering wedges are used. We find that for small variations of the EoR parameters, the differences in the multipoles are mostly in relatively small scales. Due to the lack of sensitivity for scales |$k\gtrsim 0.1\, {\rm Mpc^{-1}}$|⁠, the degeneracy between EoR parameters is only slightly resolved by including the multipoles. This is due to the relatively sparse distribution of long baselines in SKA-Low, and our conservative assumption of |$120\,$|hrs of observation time. For a complete 21 cm survey with |$\sim 10,000\,$|hrs of integration time, SKA-Low will be able to probe into the small scales |$k\sim 1\, {\rm Mpc^{-1}}$| where the multipoles can be used to further break the degeneracy.

Our work provides strong incentives for using power spectrum multipoles in clustering wedges in forthcoming EoR measurements. In addition to the advantages mentioned above, the multipoles are less susceptible to various observational systematic effects. Quantifying the covariance of the multipoles in different wedges from visibility data and using a conventional model inference framework with the multipoles are within the reach of the anticipated data quality in the near future. Robust model inference with the power spectrum multipoles and clustering wedges will provide crucial insights into the physics driving the cosmic reionization.

ACKNOWLEDGEMENTS

ZC and AP are funded by a UKRI Future Leaders Fellowship (grant MR/X005399/1; PI: Alkistis Pourtsidou). The computation underlying this work is performed on the cuillin cluster located at the Institute for Astronomy, University of Edinburgh. Apart from aforementioned packages, this work also uses pytorch (Paszke et al. 2019), numpy (Harris et al. 2020), scipy (Virtanen et al. 2020), astropy (Astropy Collaboration 2018), camb (Lewis, Challinor & Lasenby 2000), casa (CASA Team 2022), and matplotlib (Hunter 2007).

For the purpose of open access, the author has applied a Creative Commons Attribution (CC BY) licence to any Author Accepted Manuscript version arising from this submission.

DATA AVAILABILITY

The scripts for simulating reionization lightcones are publicly available as a GitHub repository⁷ The scripts for visibility simulation and the subsequent Fisher matrix analysis will be shared on reasonable request to the corresponding author.

Footnotes

REFERENCES

APPENDIX A: QUADRATIC ESTIMATOR FOR POWER SPECTRUM MULTIPOLES

In this section, we discuss the quadratic estimator formalism used in this paper for power spectrum multipoles. The quadratic estimator for cosmological power spectra is well studied in the context of CMB (e.g. Tegmark 1997; Tegmark & de Oliveira-Costa 2001), galaxy clustering (e.g. Tegmark et al. 1998), and more recently in interferometric 21 cm observations (e.g. Liu & Tegmark 2011; Kern & Liu 2021). For this work, we follow the formalism in Kern & Liu (2021) and Chen et al. (2023a), and specify our estimator in the context of the delay power spectrum from visibilities. We then generalize it to power spectrum multipoles.

A radio interferometer measures the sky emission with baselines, i.e. pairs of antennas, as visibilities (Smirnov 2011)

where |$I(l,m,f)$| is the flux density of the sky signal at frequency f and angular coordinate on the sky |$(l,m)$|⁠, and |$A(l,m,f)$| is the power response of the instrument. For a pair of antennas, the distance vector between them, denoted as |$\boldsymbol {b}$|⁠, can be written in East-North-Up (ENU) coordinate frame as

where |$\rm T$| denotes the transpose operation and |$\lambda$| is the observing wavelength.

For each instantaneous baseline, the visibilities across the frequency channels can be Fourier transformed so that

For the 21 cm signal, the coordinates |$(u,v)$| reflect the angular scale in Fourier space while the delay time |$\eta$| reflects the Fourier mode of density fluctuation along the line of sight:

where |$\boldsymbol {u} = (u,v)^{\rm T}$|⁠.

For a bandpower in the |$\alpha ^{\rm th}$||$\boldsymbol {k}$|-bin of a multipole moment |$\ell$|⁠, which we denote as |$p^\ell _\alpha$|⁠, the quadratic estimator can be written as

where |$\mathbf {V}$| is the visibility data vector, |$\dagger$| denotes the Hermitian transpose, |$\mathbf {E}_\alpha ^\ell$| is an estimation matrix discussed later, and |$\hat{b}^\ell _\alpha$| is the bias correction term. |$C_{\rm H{\small I}}$| is a factor that converts the flux density to brightness temperature unit and renormalizes the survey volume (Parsons et al. 2014)

where |$\lambda$| is the observing frequency, |$k_{\rm B}$| is the Boltzmann constant, X is the comoving distance of the 21 cm signal at the observing frequency, B is the frequency bandwidth of the observation, and |$\Omega _{\rm ps}$| is the power-squared beam area defined as

Y is the comoving length scale per frequency interval

where z is the redshift of the 21 cm signal at the observing frequency, and |$H(z)$| is the Hubble parameter.

The data vector comprises visibilities measured at each baseline at each time step in each frequency channel. For an ideal experiment where the data vector consists only of H i signal and the beam response is unitary, the 21 cm power spectrum for |$\boldsymbol {k}_\alpha$| can be written as (McQuinn et al. 2006; Parsons et al. 2012, 2014)

The power beam of the instrument is multiplied with the signal on the sky, and therefore convolves the signal in Fourier space. The mode-mixing can be written in the data vector as

where |$\mathcal {M}^\ell _{\alpha \beta }$| is the mode-mixing matrix that can be calculated given the power beam, and |$\tilde{V}^\ell _\alpha$| is the delay-transformed data vector for multipole |$\ell$|⁠:

|$\mathcal {F}_\ell$| is the Fourier transformation kernel that transforms the visibility to multipole overdensity. In this paper, we consider |$\ell = (0,2,4)$|⁠. |$\mathcal {F}_0$| is simply the Discrete Fourier Transform kernel. For higher moments, in the flat-sky and plane-parallel limit we have (Scoccimarro 2015)

where |${\rm diag}[]$| denotes the operation that diagonalizes a vector into a matrix, |$\boldsymbol {\hat{k}}_\parallel$| is a vector so that |$\boldsymbol {\hat{k}}_\parallel ^\alpha = \boldsymbol {\mu }^\alpha = k_\parallel ^\alpha /|\boldsymbol {k}_\alpha |$|⁠.

Similar to equation (A10), we can then write down the estimator when the measurement is ideal:

where |$\boldsymbol {w}_\alpha$| is a selection vector with the |$\alpha ^{\rm th}$| element being 1 and all other elements 0.

We can then factorize the estimator of equation (A6) in a way that relates to the Fourier multipole overdensity

where the bias term is neglected for now, |${\rm tr}[]$| denotes the trace of a matrix, and we have defined a mode-mixed Fourier space estimator

Assuming that the 21 cm field does not correlate for different |$\boldsymbol {k}$|-modes, the covariance matrix of the ideal signal can be written as

where |$\beta$| loops over every measured |$\boldsymbol {k}$|-mode and |$p_{\beta }^\ell$| is the power spectrum multipole at |$\boldsymbol {k}_\beta$|⁠.

The window function of the estimator is then

To remove the bias from noise and foregrounds in the estimator, simply note |$\langle \hat{p}^\ell _\alpha \rangle = C_{\rm H{\small I}}\, {\rm tr}[\mathbf {E}^\ell _\alpha \mathbf {C}]$| where |$\mathbf {C}$| is the data covariance and the bias correction term is therefore

where |$\mathbf {C}_{\rm N}$| and |$\mathbf {C}_{\rm fg}$| are the noise and foreground covariance, respectively. Throughout this paper, we adopt the foreground avoidance method (Morales et al. 2012) and only estimate the power spectrum at |$\boldsymbol {k}$|-modes where the foreground power is negligible.

The estimator can be empirically factorized as

where |$\mathbf {S}^\ell$| is a renormalization matrix to be determined later, and |$\mathbf {R}$| is a weighting matrix that can include various data analysis techniques such as frequency tapering, inverse noise covariance weighting, and foreground cleaning.

Substituting equation (A22) into equation (A20), the window function becomes

Therefore, we can form the quantity

so that

where we utilize the fact that |${\rm diag}[\boldsymbol {w}_\alpha ]$| is commutable with any matrix.

The formulation of the window function in equation (A25) is determined by the choice of |$\mathbf {R}$|⁠, which then determines |$\mathbf {H}^\ell$|⁠, and the subsequent choice of |$\mathbf {S}^\ell$|⁠. For example, a common choice of |$\mathbf {S}^\ell$| is the diagonal matrix |$\mathbf {S}^\ell _{\alpha \beta } = \delta ^{\rm K}_{\alpha \beta } (\mathbf {H}^\ell _{\alpha \beta })^{\rm -1}$| (e.g. Hamilton 1997), where |$\delta ^{\rm K}$| is the Kronecker delta. |$\mathbf {S}^\ell = (\mathbf {H}^\ell)^{-1}$| decorrelates the estimated bandpowers. Throughout this paper, we choose |$\mathbf {S}^\ell = (\mathbf {H}^\ell)^{-1/2}$| to balance the measurement errors and correlation of different k-bins.

Assuming that the effects of mode-mixing are negligible in the signal covariance, the covariance matrix of the estimator can be written as

where |$\mathbf {C}_{\rm s}$| is the signal covariance. This is different from the case in e.g. the CMB, due to the fact that thermal noise in visibility is complex. The above equation also assumes that the u–v coverage of the interferometer baselines is complete.

While the estimator formalism presented in this section is generic and takes into account the entire data vector, constructing the entire data covariance matrix is too computationally expensive for the scope of our work. We implement the estimator at each u–v grid to compute the 3D band powers and covariances, and then perform the 1D binning assuming that different u–v grids do not correlate. Since we are only using the quadratic estimator to calculate noise covariance in this work, this simplification does not impact our results. A full detection methodology from visibility data to multipoles is left for future work.

For reference, we show the full signal correlation matrix in Fig. A1, calculated using the jackknife method discussed in Section 3, and the full noise correlation matrix in Fig. A2, calculated using the quadratic estimator discussed above. The amplitude of the total data covariance matrices is shown in Section 5.2.

APPENDIX B: CONVERGENCE TESTS

In this section, we discuss the convergence of our Fisher matrix analysis. The calculation of Fisher matrix depends on the calculation of the derivatives of the summary statistics as well as the covariance with regard to the model parameters. However, the growth of ionization structures during EoR is nonlinear and non-Gaussian, leading to potential numerical instabilities in the jackknife method we use to calculate the covariance of the summary statistics. To resolve the potential instability, we average over 10 independent realizations for each set of input parameter values. We evaluate the sufficiency of the number of realizations in this section.

To test the convergence of the Fisher matrix, we perform a jackknife test by excluding a certain number of realizations at a time. The rest of the realizations are then used to calculate the Fisher matrix. We then calculate the mean and the standard deviation of the measurement errors among all possible combinations. For example, to evaluate convergence of the number of realizations |$N_{\rm realization} = 2$|⁠, we exclude |$10-N_{\rm realization} = 8$| realizations resulting in the total of 45 combinations.

The resulting relation between the number of realizations and the projected measurement errors of the parameters is shown in Fig. B1. For simplicity, only the scenario of using the multipoles in clustering wedges is shown. The error bars of the measurement errors are linearly extrapolated from |$N_{\rm realization} = 1{\!-\!}9$| to |$N_{\rm realization} = 10$|⁠. The measurements errors first decrease then plateau with increased number of realizations. When the number of realizations is small, the measurement errors include numerical instabilities which enlarge the error bars. As shown, the values for |$\sigma _{\xi }$| and |$\sigma _{{\rm log}[T_{\rm vir}]}$| plateau after around five realizations, suggesting that the results have converged for both models.

Furthermore, it is important to test that degeneracy between the parameters is consistent across realizations. To illustrate the convergence of the parameter correlation, we show the covariance of the parameter measurement for |$N_{\rm realization} = 9$| in Fig. B2. As shown, for the faint model, the correlation between the parameters is consistent across all combinations, indicating that the convergence has been reached. The bright model, while also consistent, shows larger variations across different combinations. As we have discussed in Section 6, the parameter values are more extreme for the bright model, resulting in larger instabilities. Nevertheless, the parameters are not strongly correlated in the bright model, and the variations do not affect the conclusions reached in this work.

APPENDIX C: FORECASTS FOR THE FOREGROUND WEDGE AT THE HORIZON

As mentioned in Section 7, we cannot perform multipole expansion or clustering wedge partition if the foreground wedge is at the physical horizon. For reference, we calculate the Fisher matrix for the spherically averaged monopole for the horizon case, and present the results here.

In Fig. C1, we show the forecasts for the spherically averaged monopole for the horizon case, compared against the case of foreground wedge inside the primary beam field of view. For the faint model, the |$1\sigma$| confidence interval increases by a factor of |$\sim 2.5$|⁠. For the bright model, the increase in error bars of the parameters is around a factor of 2. We note that, since most of the sampling from the baselines centre around |$\mu \sim 1$| as shown in Fig. 10, the loss in sensitivity in the monopole is not as severe as the |$\mu \gt 0.96$| wedge may suggest.