© C. P. Novaes et al. 2022

1. Introduction

Recent cosmological observations, such as distance measurements with type Ia supernovae, observations of cosmic microwave background (CMB) temperature and polarization anisotropies, and the detection of baryon acoustic oscillations (BAO) in large-scale structure spectroscopic surveys, have been key in establishing the current standard cosmological model, ΛCDM. Although a remarkable fit to existing observations, ΛCDM requires 95% of the energy content of the Universe to be in the form of dark matter (∼25%), and a puzzling component that causes the current accelerated expansion of the Universe, dark energy, with approximately 70%. Confidence in the model requires the identification of these two components of the Universe (Abdalla & Marins 2020). Otherwise, the ΛCDM model will be subject to being challenged by alternative interpretations, many of which involve modifications of our theory of gravitation on large scales.

The BAO feature, appearing as a fixed scale in both the CMB and large-scale structure (LSS) data (a geometrical probe), originates in the early Universe, when the temperature was high enough to keep photons and baryons coupled. The BAO are the imprints left by acoustic waves traveling at relativistic speed, generated by the gravitational infall of baryons (and dark matter) into the potential wells of dark matter balanced by the radiation pressure pushing out these baryons from overdense regions. The comoving distance they traveled until recombination when the baryons were released from the drag of the photons defines the BAO scale (Dodelson 2003). This characteristic scale is then imprinted in the LSS distribution and evolves with structure formation as a standard ruler that can be measured as a function of the redshift (Eisenstein et al. 2007a; Seo & Eisenstein 2007). This makes the BAO one of the most powerful and well-established tools for investigating the history of the acceleration of the expansion of the Universe (Weinberg et al. 2013).

The first statistically significant measurements of the BAO feature imprinted on the galaxy distribution were made by Eisenstein et al. (2005), who analyzed the Sloan Digital Sky Survey (SDSS) data, and by Cole et al. (2005), who used the 2dF Galaxy Survey (2dFGS). After this, several other detections of the BAO scale were performed using different biased tracers of the dark matter in several redshift ranges. These include the analyses of SDSS data at various stages of data gathering (Anderson et al. 2012, 2014; Alam et al. 2017, 2021; Ata et al. 2018), in addition to the WiggleZ Dark Energy Survey (DES) (Hinton et al. 2016) and the 6dFGS (Carter et al. 2018) analyses. More recently, data releases of the DES have been analyzed for the detection of the BAO signature through different clustering measurement methods: in configuration, in harmonic and Fourier spaces, and for a two- or three-dimensional distribution of sources (Camacho et al. 2019; Abbott et al. 2019, 2022). These results have used the common template-based fitting method to study the BAO feature from clustering estimates. Alternative approaches for measuring the BAO scale in the angular (two-dimensional) two- and three-point statistics use an empirical parametric fit technique, for example, Sánchez et al. (2011), Carnero et al. (2012), de Simoni et al. (2013), Carvalho et al. (2016, 2020), De Carvalho et al. (2018, 2020), and de Carvalho et al. (2021).

A new and very promising way of measuring BAO is through the detection of structures as traced by the redshifted diffuse 21 cm hyperfine transition line of neutral hydrogen (HI) through the so-called intensity mapping (IM) technique (Chang et al. 2008). Compared to radio surveys using the 21 cm emission to detect individual galaxies (for an example of its use in cosmological analyses, see Avila et al. 2018), limited by the low luminosity of the 21 cm line emission, the IM technique can cover a larger volume of the Universe in a much shorter time, using instruments with a relatively small, and consequently cheaper, collecting area, as discussed and evaluated in Battye et al. (2013) and Abdalla et al. (2022a) (see also Chang et al. 2008; Loeb & Wyithe 2008). The idea is to measure the overall HI brightness temperature field, similarly to what is done for CMB temperature fluctuations, but mapping the Universe as a function of redshift. This is possible given the high abundance of hydrogen, so that a hydrogen map is expected to be a powerful tracer of the underlying total matter content of the Universe. In this sense, through redshift surveys at radio frequencies, the 21 cm IM constitutes a new window for observing the Universe, and for studying the dark sector.

To exploit this new observational window, several radio instruments that are already observing or are still under construction will survey a large volume of the Universe through the 21 cm IM technique and measure the BAO feature. They include the Square Kilometer Array¹ (SKA; SKA Cosmology SWG 2020); MeerKAT² (Santos et al. 2017), a precursor of the SKA; the Five-Hundred-Meter Aperture Spherical Radio Telescope, currently the largest single-dish telescope in the world (FAST; Nan et al. 2011); the Canadian Hydrogen Intensity Mapping Experiment³ (CHIME; Bandura et al. 2014); Tianlai⁴ (Chen 2012); HIRAX (Crichton et al. 2022); and the Baryon Acoustic Oscillations from Integrated Neutral Gas Observations⁵ (BINGO; Battye et al. 2013; Bigot-Sazy et al. 2015; Abdalla et al. 2022a), which is described in the next section.

In this work we evaluate what to expect from future BINGO observations in terms of measuring the transversal (angular) BAO signal. To perform this task, we employ and compare two clustering estimators, the two-point angular correlation function (ACF) in configuration space, and the angular power spectrum (APS) in harmonic space. To our knowledge, this is the first forecasting study of the BAO detection from 21 cm IM simulations using the ACF. The anisotropic two-point correlation function estimator ξ(r_⊥, r_∥) was explored by Avila et al. (2022) for the case of an SKA-like survey employing simulations and by Kennedy & Bull (2021) for a theoretical modeling approach focusing on the MeerKAT survey. Both of them focused on the impact of the instrumental beam smoothing and foreground removal in recovering the BAO scale. Villaescusa-Navarro et al. (2017) also explored the SKA case, assessing the BAO detection through the radial power spectrum, under the presence of instrumental effects and foreground contamination, and the impact of the angular resolution imposed by a single-dish instrument. They reported that the telescope beam will compromise the detection of the isotropic BAO feature at redshifts z ≳ 1, while the radial BAO seems to be robust against the foreground removal. These published results confirm the potential of future 21 cm IM observations for detecting the BAO signature.

In the present work, we use a template-based method to model the APS and ACF estimated from two types of mock realizations, mimicking future BINGO observations, and extract the BAO information from each of them. We use the covariance matrices calculated from the mocks to construct the likelihood corresponding to each template, and using a maximum likelihood estimator, we estimate the parameters of the model for each simulation. We apply this template-fitting procedure to three sets of consecutive bins, so that we can evaluate the BAO detection at three different redshift intervals.

This paper is organized as follows. Section 2 briefly describes the main aspects of the BINGO project. Section 3 introduces the HI clustering theory and the template model constructed for each estimator. Section 4 presents the details of the cosmological 21 cm simulations, the characteristics of the instrumental noise, and the foreground contamination, as well as the foreground cleaning process employed here. Section 5 summarizes our method for the clustering measurements, the BAO fitting process, and the covariance matrix construction. The results from all our analyses are discussed in Sect. 6, and the conclusions are summarized in Sect. 7.

2. The BINGO telescope

The BINGO radio telescope is in construction at a site (for the site selection process, see Peel et al. 2019) located in Paraíba State, northeastern Brazil (latitude: 7°2′27.6″ S; longitude: 38°16′4.8″ W; altitude: 350 to 460 m). Its main scientific goal is to measure the BAO signal imprinted on the 21 cm distribution in the redshift range 0.127 ≤ z ≤ 0.449 (Abdalla et al. 2022a). This goal will be achieved through an IM survey covering 5324 square degrees, using the telescope in sky-transit mode, with a 14.75° wide declination strip centered at a declination δ = −15°.

The telescope has a crossed-Dragone design, with a 40 m diameter primary paraboloid and a 34 m diameter secondary hyperboloid. During BINGO Phase 1, the telescope will operate with a focal plane containing 28 horns (Wuensche et al. 2022; Abdalla et al. 2022b). Each horn is sensitive to circular polarization and is coupled to a correlation receiver with four amplifier chains. They are expected to perform at an expected system temperature T_sys = 70 K. The optical design will produce a suitable angular resolution for the BAO signal, with a full width at half maximum θ_FWHM = 40′ of the beam in the central frequency ν = 1120 MHz. In a second phase, we intend for BINGO to operate with 28 additional horns, totaling 56 horns.

The redshift range covered by BINGO corresponds to the frequency interval 980 − 1260 MHz, which is quite complementary to instruments such as CHIME (400 − 800 MHz; Bandura et al. 2014). A more detailed description of the BINGO project, its scientific goals, and instrument status is available in previous publications of the collaboration (Abdalla et al. 2022a; Wuensche et al. 2022).

3. Modeling the BAO signal

3.1. H I clustering

The intensity mapping of the 21 cm signal measures the brightness temperature of the average HI intensity in a given volume of the Universe, which can be written as function of the redshift, z, as (see, e.g., Battye et al. 2013; Hall et al. 2013)

$$\begin{array}{c}\hfill {\overline{T}}_{\mathrm{HI}}(z)=188h{\mathrm{\Omega }}_{\mathrm{HI}}(z)\frac{{(1+z)}^2}{E(z)}\mathrm{mK},\end{array}$$(1)

where H₀ = 100 h km s⁻¹ Mpc⁻¹ is the Hubble constant, E(z) = H(z)/H₀, and Ω_HI is the HI density parameter. Fluctuations of the HI brightness temperature, δT_HI, are a biased tracer of the dark matter fluctuations, δ. Their Fourier transform can be written in terms of the growth function D(z) as

$$\begin{array}{c}\hfill \delta T_{\mathrm{HI}}(\mathit{k},z)=D(z){\overline{T}}_{\mathrm{HI}}(z)b_{\mathrm{HI}}(z)\delta (\mathit{k},0),\end{array}$$(2)

where b_HI(z) is the bias factor as a function of redshift z, and δ(k, 0) is the underlying dark matter distribution at z = 0. In this way, the HI power spectrum is given by $P_{\mathrm{HI}}\approx {[D(z){\overline{T}}_{\mathrm{HI}}(z)b_{\mathrm{HI}}(z)]}^{2}P(k)$, where P(k) is the matter power spectrum at z = 0.

As detailed by Battye et al. (2013) and Seehars et al. (2016), since the 21 cm IM surveys the Universe within tomographic bins of redshift, we can project the three-dimensional quantity $\delta T_{\mathrm{HI}}(\chi (z)\widehat{n})={\overline{T}}_{\mathrm{HI}}(z){\delta }_{\mathrm{HI}}(\chi (z)\widehat{n})$ on the sky by integrating it along the line of sight $\widehat{n}$ as

$$\begin{array}{c}\hfill \delta T_{\mathrm{HI}}(\widehat{n})={\displaystyle \int }\mathrm{d}z\varphi (z){\overline{T}}_{\mathrm{HI}}(z){\delta }_{\mathrm{HI}}(\chi (z)\widehat{n}),\end{array}$$(3)

where ${\delta }_{\mathrm{HI}}(\chi (z)\widehat{n})$ is the density fluctuation of neutral hydrogen at the comoving distance to redshift z, χ(z), and ϕ(z) is the projection kernel (a window function of observation). Here, we assume a top-hat kernel, that is, ϕ(z) = 1/(z_max − z_min) for z_min < z < z_max and ϕ(z) = 0 outside the redshift bin. Decomposing the HI temperature fluctuations in spherical harmonics, we have (see also Sobreira et al. 2011; Costa et al. 2022)

$$\begin{array}{c}\hfill \delta T_{\mathrm{HI}}(\widehat{n})={\displaystyle \sum _{\ell =0}^{\infty }\sum _{m=-\ell }^{\ell }{a}_{\ell m}{Y}_{\ell m}(\widehat{n}),}\end{array}$$(4)

with the a_ℓm harmonic coefficients written as

$$\begin{array}{c}\hfill a_{\ell m}=4\pi i^{\ell }{\displaystyle \int \mathrm{d}z\varphi (z){\overline{T}}_{\mathrm{HI}}(z)\int \frac{{\mathrm{d}}^{3}k}{{(2\pi )}^3}{\delta }_{\mathrm{HI}}(\mathit{k},z)j_{\ell }(k\chi (z))Y_{\ell m}^{\ast }(\mathit{k}),}\end{array}$$(5)

where j_ℓ is the spherical Bessel function. Therefore, from the above equations, we can define the APS of the temperature fluctuations, C_ℓ = ⟨ |a_ℓm|²⟩, as a function of the matter power spectrum, such that (for the analogous case of a galaxy distribution, see also Loureiro et al. 2019)

$$\begin{array}{c}\hfill C_{\ell }^{\mathit{ij}}=\frac{2}{\pi }{\displaystyle \int }\mathrm{d}k{W}_{\mathrm{HI},\ell }^i(k)W_{\mathrm{HI},\ell }^j(k)k^{2}P(k),\end{array}$$(6)

where the indices i and j denote two tomographic bins. For i = j and i ≠ j, we obtain auto- and cross-APS, respectively (see also Costa et al. 2022). The redshift dependence is given by the window function

$$\begin{array}{c}\hfill W_{\mathrm{HI},\ell }^i(k)={\displaystyle \int }\mathrm{d}z{b}_{\mathrm{HI}}(z)\varphi (z){\overline{T}}_{\mathrm{HI}}(z)D(z)j_{\ell }(k\chi (z)).\end{array}$$(7)

In a similar way to the APS, we can also define the ACF, ω(θ), defined as the probability of finding galaxies separated by an angle θ. For a tomographic bin around z, the ACF can be defined as a function of the matter three-dimensional spatial correlation function ξ(r) (Sobreira et al. 2011; Crocce et al. 2011; Sánchez et al. 2011; de Simoni et al. 2013)

$$\begin{array}{cc}\hfill \omega (\theta )& \equiv ⟨\delta T_{\mathrm{HI}}(\widehat{n})\delta T_{\mathrm{HI}}(\widehat{n}+\mathit{\theta })⟩\hfill \end{array}$$(8)

$$\begin{array}{cc}& ={\displaystyle \int \mathrm{d}{z}_{1}f(z_1)\int \mathrm{d}{z}_{2}f(z_2)\xi (r(z_1,z_2,\theta ),\overline{z}),}\hfill \end{array}$$(9)

where $f(z)=b_{\mathrm{HI}}(z)\varphi (z){\overline{T}}_{\mathrm{HI}}(z)$ and r(z₁, z₂, θ) is the radial comoving distance between spatial fluctuations δ_HI at redshifts z₁ and z₂ separated by an angle θ. This equation neglects the time evolution of the matter correlation function inside the bin, so that it can be evaluated at a given $\overline{z}$ (the average redshift of the bin) as

$$\begin{array}{c}\hfill \xi (r,\overline{z})=\frac{1}{2\pi }{\displaystyle \int }\mathrm{d}k{k}^2{j}_0(rk)P(k,\overline{z}),\end{array}$$(10)

where j₀ is the spherical Bessel function of zeroth order and $P(k,\overline{z})$ is the matter power spectrum at redshift $\overline{z}$.

Alternatively, following Sobreira et al. (2011) and Crocce et al. (2011) and using Eq. (4), we can still obtain the ACF as

$$\begin{array}{c}\hfill \omega (\theta )=\langle {\displaystyle \sum _{\ell =0}^{\infty }\sum _{\ell =-\ell }^{\ell }{a}_{\ell m}{a}_{\ell \prime m\prime }{Y}_{\ell m}(\widehat{n})Y_{\ell \prime m\prime }(\widehat{n}+\mathit{\theta })\rangle ,}\end{array}$$(11)

and, as a function of the APS, as

$$\begin{array}{cc}\hfill \omega (\theta )& ={\displaystyle \sum _{\ell m}{C}_{\ell }{Y}_{\ell m}(\widehat{n})Y_{\ell m}(\widehat{n}+\mathit{\theta })}\hfill \end{array}$$(12)

$$\begin{array}{cc}& ={\displaystyle \sum _{\ell }{C}_{\ell }\frac{2\ell +1}{4\pi }{P}_{\ell }(cos\theta ),}\hfill \end{array}$$(13)

where P_ℓ are the Legendre polynomials.

The theoretical auto- and cross-C_ℓs employed in this paper, both as input to the 21 cm log-normal simulations (see Sect. 4.1.1) and to construct the BAO templates, were calculated using the Unified Cosmological Library for C_ℓs (UCLCL) code (Loureiro et al. 2019; McLeod et al. 2017). This code implements Eqs. (6) and (7) using the CLASS Boltzmann code (Lesgourgues 2011; Blas et al. 2011) to estimate the primordial power spectra and transfer function. To account for all processes involved in the evolution of the Universe, we used the window function $W^{\text{Tot},i}(k)=W_{\text{HI},\mathcal{l}}^i(k)+W_{\text{RSD},\mathcal{l}}^i(k)$. In addition to Eq. (7), it also includes a term describing the redshift space distortion (RSD) effect, $W_{\text{RSD},\mathcal{l}}^i(k)$. A detailed description of how the two terms are implemented in the UCLCl code can be found in Loureiro et al. (2019).

3.2. BAO template

The extraction of the BAO features from data clustering estimates is commonly performed by fitting a template model, derived from a parameterization of the matter power spectrum (see, e.g., Anderson et al. 2014, and references in the Introduction). Following the same approach, we write this parameterization as a function of the linear power spectrum, P_lin, and the no-wiggle (no BAO feature) power spectrum, P^nw,

$$\begin{array}{c}\hfill P^{\mathrm{temp}}(k)=[P^{\mathrm{lin}}(k)-P^{\mathrm{nw}}(k)]e^{-k^2{\mathrm{\Sigma }}_{\mathrm{nl}}^2}+P^{\mathrm{nw}}(k).\end{array}$$(14)

We employed the nbodykit⁶ code to obtain both P^lin(k) and P^nw(k), using the transfer functions calculated by the CLASS code and the analytic calculation by Eisenstein & Hu (1998), respectively. The exponential term, $e^{-k^2{\mathrm{\Sigma }}_{\text{nl}}^2}$, takes into account the effect of nonlinear structure growth by damping the signal around the BAO scale. The damping scale, ${\mathrm{\Sigma }}_{\text{nl}}^2=({\mathrm{\Sigma }}_{\perp }^2+{\mathrm{\Sigma }}_{\parallel }^2)/2$, is written in terms of the components along (Σ_∥) and across (Σ_⊥) the line of sight, taking into account the RSD effect. Seo & Eisenstein (2007) have shown that these components can be written as Σ_⊥ = 10.4D(z)σ₈ (prediction for real space) and Σ_∥ = (1 + f)Σ_⊥, where f is the growth rate of cosmic structures (see also the discussion by Chan et al. 2018; Ata et al. 2018). Because we already account for the RSD in the window function $W_{\text{HI},\mathcal{l}}^i$ when projecting P^temp(k) into the APS (Eqs. (6) and (7)), we can consider Σ_nl = Σ_⊥. For our fiducial cosmology (see Sect. 4.1.1), we then fixed the damping scale at Σ_nl = 7.6, 7.1, and 6.7 h⁻¹ Mpc, appropriate values for the lower, intermediate, and higher redshift intervals into which we split the BAO fitting analyses, as discussed below. Then, from the projection of P^temp(k) into C^temp(ℓ), we constructed the template used for the APS analysis,

$$\begin{array}{c}\hfill C(\ell )=B{C}^{\mathrm{temp}}(\ell /\alpha )+{\displaystyle \sum _q{A}_q{\ell }^q,}\end{array}$$(15)

where α, B and A_q are free parameters, the last two of them intended to absorb linear and nonlinear bias effects, noise, and uncertainties in the RSD, in addition to any other difference between the data points and the full shape template, such as those introduced by systematic effects.

For the ACF analyses, the template model was constructed by substituting the same C^temp(ℓ) in Eq. (12) to calculate ω^temp(θ), so that we guarantee the consistency of the two templates, obtaining

$$\begin{array}{c}\hfill \omega (\theta )=B{\omega }^{\mathrm{temp}}(\alpha \theta )+{\displaystyle \sum _q\frac{A_q}{{\theta }^q},}\end{array}$$(16)

where α, B and A_q are free parameters, as in Eq. (15). The last term in both templates, appearing as functions of ℓ and θ, can be written so that the number q of A_q parameters optimizes our fitting pipeline. Section 6.4 shows results from testing different degrees of freedom, that is, different choices for minimum and maximum values to run the q index, for the C_ℓ and ω(θ) templates.

In both cases, the α parameter, the most important parameter of the analysis, is the so-called shift parameter, associated with the change in the BAO peak position with respect to a fiducial cosmology,

$$\begin{array}{c}\hfill \alpha =\frac{D_{\mathrm{A}}(z)/r_{\mathrm{d}}}{{(D_{\mathrm{A}}(z)/r_{\mathrm{d}})}_{\mathrm{fid}}},\end{array}$$(17)

where D_A is the angular diameter distance and r_d is the sound horizon scale at the drag epoch. Then, α characterizes any observed deviation with respect to the model, so that α > 1 (α < 1) indicates a shift of the acoustic peak to smaller (larger) scales. Because the fiducial cosmology adopted to model the BAO templates and to generate the synthetic data are the same, we expect to find α ≈ 1 when fitting them to the APS and ACF clustering computed from the simulations. This allowed us to test the method and to predict what to expect from the BAO analysis of the future BINGO data.

4. Synthetic data preparation

In this section we describe the different sets of simulations we employed in the analyses. Sky maps were produced in the HEALPix pixelization scheme (Gorski et al. 2005), with N_side = 256. The semi-realistic mock data sets (referred to as BINGO-like simulations) include instrumental noise, beams effects, the BINGO sky coverage scan, and simulated foreground signals that will contaminate 21 cm observations in the BINGO frequency range. A foreground cleaning pipeline, outlined in Sect. 4.2, was applied to a small set of BINGO-like simulations, constructed from a subset of the FLASK mocks, to mitigate the effect of the foreground signal on the 21 cm signal and estimate the residual contamination in the recovered maps.

4.1. BINGO-like simulations

4.1.1. Cosmological signal

The 21 cm signal simulations we employed here were produced in two ways. In most of the tests, we used a large set of (fast) log-normal simulations of the 21 cm signal generated by the FLASK⁷ code (Xavier et al. 2016); hereafter FLASK mocks. We also tested the BAO fitting pipeline over a smaller data set based on the density contrast from N-body simulations; hereafter N-body mocks. The two types of simulations are described below.

(a) Fast log-normal distributions. The first set of 21 cm IM simulations consisted of full-sky log-normal realizations generated with FLASK. This is a publicly available code able to produce two- or three-dimensional tomographic realizations of an arbitrary number of random astrophysical fields and reproducing the desired cross-correlations between them. FLASK took as input the auto- and cross-APS, C^ij(ℓ), previously calculated for each of the i and j redshift slices, and created two-dimensional HEALPix maps of correlated log-normal realizations of the projected 21 cm signal in each redshift bin (z-bin). We used the fiducial C^ij(ℓ) computed using the UCLCl code. Cosmological parameters matched those from WMAP 5-year results (Table 2, Five-Year Mean values, of Dunkley et al. 2009), Ω_m = 0.26, Ω_b = 0.044, Ω_Λ = 0.74, and H₀ = 72 km s⁻¹ Mpc⁻¹, for consistency with the N-body mock simulations described below. Although more recent parameter constraints have been obtained with the Planck satellite (Planck Collaboration VI 2020), our pipeline and results do not strongly depend on the exact cosmological parameter values. Here we are not interested in constraining these parameters, but to show the consistency of our method and the detection of BAO.

We generated a total of 1500 log-normal realizations, each of them corresponding to 30 HI maps, one for each tomographic z-bin. Despite our ad hoc choice of 30 channels (δν = 9.33 MHz width each), the BINGO hardware and data analysis pipeline allows different choices for the number of channels. Consequently, their width can be chosen according to what is most suitable for cosmological analyses, for example, as a function of the efficacy of a foreground cleaning procedure (Mericia et al. 2022). Tests of how the BAO measurement is impacted when a different number of channels is used are left for future work.

(b) N-body simulations. A second set of mocks was generated using the friend-of-friends (FoF) halo catalogs from the HR4 N-body simulation (Kim et al. 2015). The HR4 has a box size of 3150 Mpc h⁻¹ and uses WMAP 5-year cosmological parameters.

First, we selected 100 random observer locations in the simulation box to mimic 100 different realizations in just one simulation. We then constructed the light cone catalog from snapshots at z = 0.1, 0.15, 0.2, 0.3, 0.4, and 0.5. According to the designated frequency bin, the halos from the closest two snapshots were selected. These halos were further filtered randomly, and the selection probability was proportional to the radial distance of this bin and the corresponding comoving distance at the redshift of the snapshot. Therefore, the generated light cone halo catalog is a random mixture of the snapshots. We validated the method by comparing the angular power spectrum and halo mass function to the light cone halo catalog provided by HR4. The two results are consistent and meet our requirements.

We constructed the full-sky 21 cm brightness temperature map from the light cone halo catalog following the method described in Zhang et al. (2022). By further considering redshift distortion effects, we constructed 100 realizations of full-sky mock maps, each containing 30 redshift bins. Because BINGO exploits only a strip of the sky, we obtained independent realizations by randomly rotating the full-sky maps before including realistic characteristics of the BINGO observations. Using four independent rotations (the spherical rotations as performed by the healpy package), we generated another 400 mock maps in the BINGO footprint using the 100 original full-sky maps. After performing these simulations, we had 1500 FLASK and 500 N-body mocks on which to test our BAO detection pipelines.

4.1.2. Foreground signals, instrumental effects, and sky coverage

Following previous BINGO papers by Fornazier et al. (2022) and Liccardo et al. (2022), our simulations included the contribution of seven foreground components, including Galactic synchrotron, free-free, thermal dust, and anomalous microwave (AME) emissions, extragalactic thermal and kinetic Sunyaev–Zel’dovich (SZ) effects, and unresolved radio point sources, all generated using a recent version of the Planck Sky Model software (PSM; Delabrouille et al. 2013). The specific configuration of the PSM code that was used to generate our foreground emissions is as described bellow.

The synchrotron emission was based on the 408 MHz all-sky map produced by Remazeilles et al. (2015), extrapolated to BINGO frequencies with a spatially variable spectral index map following the model derived by Miville-Deschênes et al. (2008). Free-free emission was simulated using a template given by the Hα emission map from Dickinson et al. (2003), with a uniform frequency scaling over the sky and slowly varying with frequency. For the AME, which is usually described as radio emission produced by the rapid rotation of electric dipoles associated with small dust grains, we used a high-resolution thermal dust template from Planck observations, scaled to low frequency according to the ratio of the AME and thermal dust as found by Ade et al. (2016), and extrapolated it to BINGO frequencies using a single emission law. A template for the thermal dust emission was obtained by applying the generalized needlet internal linear combination (GNILC) code to Planck 2015 data. This is the component separation code that was also employed here to clean the 21 cm maps (see Sect. 4.2). Dust spectral index and temperature maps, used to extrapolate thermal dust emission across frequencies using a single modified blackbody in each pixel, were obtained from fits of GNILC dust maps at different frequencies. Although AME and thermal dust emission are subdominant at BINGO frequencies, both were taken into account here for completeness.

The extragalactic emission from thermal and kinetic SZ effects was modeled according to prior knowledge of existing galaxy clusters for a number density as a function of mass and redshift predicted by the cosmological model of interest given by Tinker et al. (2008). The unresolved radio point source component was simulated using catalogs of observed sources from 850 MHz and 4.85 GHz. A population of sources below the detection threshold was simulated on the basis of theoretical number counts. A more complete discussion of the foreground simulation for the BINGO frequency can be found in Fornazier et al. (2022) and Liccardo et al. (2022).

We also took the expected thermal (white) noise level per pixel into account by employing a system temperature of 70 K. This level was estimated considering 28 horns, each of them observing at a constant elevation, achieving a survey area of 5324 deg² (a sky fraction of f_sky ∼ 0.13). We assumed five years of observation and the horn arrangement designed for the Phase 1 of observations, as discussed in previous BINGO publications (Wuensche et al. 2022; Abdalla et al. 2022b; Liccardo et al. 2022). For a good sampling of the map in declination, the position of the horns will be shifted to displace their pointing on the sky by a fraction of a beam width in elevation. This will homogenize the observation time over the innermost BINGO area. A detailed description of the noise simulations employed here can be found in Fornazier et al. (2022).

Although this paper does not account for the low-frequency (1/f) noise, we are aware that it is an important complication in real data and can be very detrimental to HI IM data if not efficiently removed (Bigot-Sazy et al. 2015; Li et al. 2021). The 1/f noise is correlated across the frequency band, and when the observed signal is projected onto a map, it appears as stripes, that is, as large-scale spatial fluctuations (Harper et al. 2018). In this sense, because this instrumental effect has the potential to impact our analyses at the lowest redshifts, where the BAO scale is comparable to the 15° stripe of the BINGO coverage, it is being addressed in various working groups across the collaboration and our findings will be presented in upcoming works of the BINGO series. Moreover, we also recall that another instrumental effect relevant for intensity mapping experiments is the leakage of polarized foregrounds into the intensity data (Alonso et al. 2015; Cunnington et al. 2021). Unlike the foregrounds described above, this polarization leakage is not expected to have a smooth dependence on frequency, which makes the cleaning process more complex (Carucci et al. 2020). In this sense, although BINGO is designed to minimize the impact of this effect (Wuensche et al. 2022; Abdalla et al. 2022b), future work will investigate its impact on the foreground cleaning efficiency, as well as on the BAO recovery.

Because all simulations were generated as full-sky maps, we selected the BINGO region using an appropriate mask that accounted not only for the expected sky coverage, but also cut out a region with strong Galactic foreground signal. To avoid the impact of sharp edges in the masked region, the mask was apodized with a cosine square transition of 5° using the NaMaster code (Alonso et al. 2019, see also Sect. 5.1.1). A detailed description of the procedure to create this mask is provided in Mericia et al. (2022). Figure 1 shows an illustrative example of a 21 cm IM simulation in the lowest redshift bin 0.127 < z < 0.138 in the BINGO region.

Fig. 1. Illustrative example of a FLASK simulation at the lowest redshift bin (0.127 < z < 0.138), convolved with a Gaussian beam of θFWHM = 40 arcmin, after selecting the sky region (apodized mask) considered in all the analyses. The color indicates the brightness temperature in μK. The Mollweide projection corresponds to celestial coordinates.

We also accounted for the angular smoothing effect introduced by the BINGO beam, with θ_FWHM = 40 arcmin, which we approximated by a Gaussian of the same width for all the frequency bins. Each simulated sky map was convolved with such a Gaussian beam before thermal noise was added. In addition, we investigated the impact of θ_FWHM scaled with the frequency ν according to θ_FWHM(ν) = θ_FWHM(ν₀) × ν₀/ν, where ν₀ = 1120 MHz and θ_FWHM(ν₀) = 40 arcmin (Bigot-Sazy et al. 2015). Although in both cases we approximated the BINGO beam by a Gaussian, there will be contribution from structures such as the side lobes, discussed in the companion papers Wuensche et al. (2022) and Abdalla et al. (2022b), as well as the possible coupling of the horns, which is still to be investigated.

4.2. Foreground cleaning process

One of the main challenges of 21 cm IM observations is the efficient subtraction of the foreground contamination, whose amplitude can be up to 10⁴ larger than the HI signal. Even when an efficient cleaning technique is employed, a residual contamination is expected to remain in the observations, and therefore it is important to evaluate their impact on the cosmological analyses. In this work, we used the GNILC method (Remazeilles et al. 2011), a nonparametric component separation technique, to perform the foreground cleaning process. This method has been tested in previous work and showed an excellent performance when applied to 21 cm IM observations (Olivari et al. 2016) and, more specifically, to the BINGO-simulated data (see Liccardo et al. 2022; Fornazier et al. 2022; Mericia et al. 2022).

The temperature signal, x(p), represented by a vector of dimension n_ch (where n_ch is the number of BINGO channels) for each sky pixel p, can be modeled as

$$\begin{array}{c}\hfill \mathit{x}(p)=\mathit{s}(p)+\mathit{n}(p)+\mathit{f}(p),\end{array}$$(18)

where s(p) is the 21 cm signal, n(p) is the instrumental noise, and f(p) is the total foreground signal. All contributions to x(p) have the same dimension n_ch. To distinguish foreground emission from 21 cm signals, GNILC exploits the fact that foreground emissions are highly correlated between frequencies channels and thus effectively rely on a few independent spectral degrees of freedom, while the 21 cm signal, probing shells at different redshifts, is only weakly correlated between frequency bands. Using a prior (theoretical template) of the 21 cm and thermal noise power spectra, GNILC evaluates the effective dimension of the foreground subspace locally in pixel and harmonic space, as identified by a principal component analysis of empirical channel data covariance matrices, computed on a redundant basis (a tight frame) of a type of wavelets on the sphere (called “needlets”). With this approach, GNILC does not use any prior information on the foreground components, except for the safe assumption that their emission between frequencies is strongly correlated. Then, GNILC projects the data x(p) out of the foreground subspace.

The 21 cm plus noise signal was then reconstructed for each wavelet scale by applying an ILC filter W to the data, $\widehat{\mathit{s}}=\mathit{Wx}$, in the subspace orthogonal to that of the foregrounds. The ILC filter was constructed in such a way as to preserve the 21 cm signal while filtering out the foreground contamination. The complete reconstructed 21 cm plus noise map at each frequency was finally obtained by adding the maps corresponding to each wavelet scale. A detailed discussion of the use of the GNILC method for 21 cm foreground cleaning can be found in Olivari et al. (2016) and Fornazier et al. (2022).

Finally, we emphasize that the reconstructed 21 cm maps contain a residual contribution from the foreground components along with the noise,

$$\begin{array}{c}\hfill \widehat{\mathit{s}}\simeq \mathit{s}+\mathit{Wn}+\mathit{Wf},\end{array}$$(19)

where Wn and Wf are the residual noise and foreground contribution to the reconstructed 21 cm signal, respectively. In the case of simulations for which we know the exact input foreground contamination, as is the case here, we can then estimate its residual contribution over the recovered cosmological signal. Here, we used this facility to include a realistic foreground residual contribution to the simulations to evaluate its impact on the BAO fitting process, avoiding the need to apply GNILC to each one of them (see Sect. 6.3).

5. Method

This section summarizes the method we employed to extract the BAO information from BINGO-simulated HI maps. We describe (1) how the APS and ACF were measured from the simulated maps, (2) how the covariance matrices were calculated, and (3) the fitting method we used to estimate the α parameter from each of the simulations.

5.1. Clustering measurements

5.1.1. Angular power spectrum

The partial-sky coverage is a common situation for CMB and galaxy surveys, whether due to the observation strategy or to cuts of regions with strong foreground contamination. The APS computation in these cases is not as straightforward as it would be in full-sky analyses, where C_ℓs are computed by expanding the radiation field in spherical harmonics and averaging the a_ℓm coefficients,

$$\begin{array}{c}\hfill C_{\ell }=\frac{1}{2\ell +1}{\displaystyle \sum _{m=-\ell }^{m=+\ell }}{|a_{\ell m}|}^2.\end{array}$$(20)

The effect of a mask when calculating the a_ℓm coefficients over partial sky is to mix (or couple) different multipoles so that the estimated APS, ${\widehat{C}}_{\ell }$, can be written in terms of the true one, C_ℓ, through

$$\begin{array}{c}\hfill {\widehat{C}}_{\ell }={\displaystyle \sum _{\ell \prime }}{\mathcal{M}}_{\ell \ell \prime }{C}_{\ell \prime },\end{array}$$(21)

where ℳ_ℓℓ′ is the multipole mixing matrix, which depends on the mask geometry (Hivon et al. 2002).

Here we estimated the APS using the so-called pseudo-C_ℓ formalism as implemented in the NaMaster code⁸ (Alonso et al. 2019). The mixing matrix was estimated analytically, according to the shape of the mask, and was used to obtain an unbiased estimate of C_ℓ (for examples of application of this formalism and/or this code to different cases, see, e.g., Brown et al. 2005; Moura-Santos et al. 2016; Loureiro et al. 2019; Nicola et al. 2020; Anderson et al. 2022).

However, the loss of information due to the mask usually makes the mixing matrix ill-conditioned and, for this reason, singular, with no direct inversion. To deal with this, NaMaster allows the usage of bandbowers, calculating a binned and invertible version of ℳ_ℓℓ′. Here, we considered a linear binning, calculating $C_{\mathcal{l}}^i$ at bins of width Δℓ, centered at ℓ. This procedure also helps making the corresponding covariance matrix more diagonal. Different Δℓ widths were tested so that we were able to find the more appropriate width that was to be used in the more realistic analyses, composing the fiducial configuration, defined when running the several robustness tests.

5.1.2. Angular correlation function

We measured the ACF from each of the simulations using the following equation (de Simoni et al. 2013):

$$\begin{array}{c}\hfill {\omega }^i(\theta )=\frac{{\sum }_{pp\prime }\delta T_{\mathrm{p}}\delta T_{p\prime }w{t}_{\mathrm{p}}w{t}_{p\prime }}{{\sum }_{pp\prime }w{t}_{\mathrm{p}}w{t}_{p\prime }},\end{array}$$(22)

with δT_p = T_b − ⟨T_b⟩, where T_b is the HI brightness temperature at each pixel p, and θ is the angular separation between pixels p and p′. The parameter wt_p is associated with the observed sky region (mask information; Fig. 1), receiving value 1 when pixel i contains valid data and 0 when its area has not been observed or if it is located in the Galactic region and is removed from the analyses. The ω(θ), as given by Eq. (22), were calculated at equally spaced θ values, with a bin width Δθ, using the public code TreeCorr⁹ (Jarvis et al. 2004). In the same way as for the APS analyses, the best Δθ width for our fiducial configuration was defined after a few different values were tested.

5.2. Covariance matrix

We used N = 1500 FLASK and N = 500 N-body mock realizations to estimate the covariance matrices as

$$\begin{array}{c}\hfill {\mathtt{C}}_{\mathit{km}}^{\mathit{ij}}=\frac{1}{N-1}{\displaystyle \sum _{n=1}^N}(X_k^{i,n}-{\overline{X}}_k^i)(X_m^{j,n}-{\overline{X}}_m^j),\end{array}$$(23)

where k and m indices ran over the different ℓ or θ bins, and i and j indicate the z-bin; $X_k^i$ is the clustering measurement, $C_{\mathcal{l}}^i$ or ωⁱ(θ), and ${\overline{X}}_k^i$ is the respective average over the simulations.

It is well known that the statistical noise in the covariance matrix estimated from simulations introduces a bias on its inverse, ${[{\mathtt{C}}_{\mathit{km}}^{\mathit{ij}}]}^{-1}$, the precision matrix. The bias is mainly dependent on the number N of simulations used to construct ${\mathtt{C}}_{\mathit{km}}^{\mathit{ij}}$ and the number of entries N_p of the data vectors X. However, as proposed by Hartlap et al. (2007) and as used, for example, in the analyses of Camacho et al. (2019), Ata et al. (2018), and Anderson et al. (2014), an unbiased version of the precision matrix can be obtained by rescaling it as

$$\begin{array}{c}\hfill {[{\mathtt{C}}_{\mathit{km}}^{\mathit{ij}}]}^{-1}\to \frac{N-N_{\mathrm{p}}-2}{N-1}{[{\mathtt{C}}_{\mathit{km}}^{\mathit{ij}}]}^{-1}.\end{array}$$(24)

The rescaled version of the precision matrix was employed in all our analyses.

5.3. Parameter inference

The fit of the template models to the estimated clustering measurement goes through two (recursive) steps: one of them is to fit the nuisance parameters, B and A_q, appearing linearly in both template models, and the other step is to use the maximum likelihood estimator (MLE) to evaluate α, as described in Eqs. (15) and (16). Following Chan et al. (2018) and Camacho et al. (2019), we performed a least-square fit to find the best-fit A_q values, and, fixing them, we repeated the procedure to obtain the best-fit B value, in both cases by minimizing the χ² defined as

$$\begin{array}{c}\hfill {\chi }^2(\lambda )={\displaystyle \sum _{k,m}[X_k^i-x_k^i(\lambda )]{[{\mathtt{C}}_{\mathit{km}}^{\mathit{ij}}]}^{-1}[X_m^j-x_m^j(\lambda )],}\end{array}$$(25)

where λ = {α, B, A_q} and $x_k^i$ is the template model for the ith z-bin and the kth ℓ/θ bin. We then maximized the likelihood function,

$$\begin{array}{c}\hfill L(\alpha )\propto exp(-{\chi }^2/2),\end{array}$$(26)

now dependent on only one parameter, since χ(λ) is reduced to χ(α), allowing us to estimate the best-fit α value (see also Anderson et al. 2014). In practice, this value was obtained through a grid search, so that the least-square fit of the nuisance parameters, A_q and B, was repeated for each α value we evaluated.

Taking advantage of a tomographic approach, we were able to combine redshift bins so that the significance of the BAO detection can increase. We fit one α parameter through a joint analysis applying the MLE over a set of N_z consecutive z-bins. In particular, we fit the nuisance parameters, B and A_q, to each z-bin individually, and, fixing their values, we find one best-fit α for each set of N_zz-bins.

The 1σ error for the α estimates, σ_α, is defined as the deviation from the maximum likelihood point, at ${\chi }_{\text{min}}^2$, by Δχ² = 1. It is worth noting that this Δχ² = 1 rule is valid in the case of a Gaussian likelihood, L(p), for only one parameter, in our case, p = α (see Chan et al. 2018, and references therein, for details). In principle, this can be employed in our analyses, as commonly seen in literature (e.g., Ata et al. 2018; Carter et al. 2018; Abbott et al. 2019, 2022), but we have to be careful with deviations from Gaussian distribution, which can lead to σ_α values that underestimate the error bars. In the next section, we discuss and evaluate the validity of this rule for each clustering statistic used here.

6. BAO fitting results

In this section we present our findings from analyzing FLASK and N-body mock realizations. We start with an ideal case, analyzing pure 21 cm simulations, only accounting for the BINGO sky coverage and fixed 40 arcmin beam (Sect. 6.1). The contribution of thermal noise and foreground residuals are then considered one at a time, so that we can investigate the impact of each one separately, as well as the impact of a redshift-dependent beam size (Sects. 6.2 and 6.3). We also present the results from an extensive (but not exhaustive) list of robustness tests evaluating several aspects of the method employed here to extract the BAO signal from BINGO-like simulations, using FLASK mocks as the cosmological signal (Sect. 6.4). These tests allow finding the most appropriate configuration (hereafter called fiducial configuration), optimizing the analysis. In particular, the fiducial configuration is defined by

• (a)

  Angular and multipole binning: Δθ = 0.5° and Δℓ = 10,

• (b)

  A_q parameters for C_ℓ and ω(θ) templates: q = −1, 0, 1, 2 and q = 0, 1, 2, respectively.

The 500 N-body mocks were analyzed using only this fiducial configuration (Sect. 6.5).

The choice for the range of scales, θ, and ℓ intervals is also important. Because the BINGO coverage is limited to a declination strip of ∼15°, the C_ℓ determination at the smallest multipoles (ℓ ≲ 15) is compromised. Then, applying the BAO fitting to the (high signal-to-noise ratio) mean C_ℓ from the 1500 FLASK mocks, we chose ℓ_min ≈ 32, for which α estimate is the least biased with the smallest σ_α error. In this way, the choice of the multipole range can also account for any artifact introduced by the mask geometry, which would be more evident in the mean clustering measurement. Details about the fitting procedure over the mean APS and ACF are presented in the next subsections.

For all the z-bins we considered the same minimum multipole, but the maximum multipole was chosen according to the multipole ranges at which the BAO wiggles are concentrated, which are more spread out over higher multipoles the higher the redshift (Fig. 2). In this sense, after a few tests over the mean C_ℓ, we chose 141 ≲ ℓ_max ≲ 401 for lower to higher redshifts, so that we used only the range containing the BAO wiggles, avoiding to increase the number of data points with multipoles that do not bring information about the BAO signal (see discussion in Sect. 5.2).

Fig. 2. Theoretical APS calculated for three of the 30 tomographic bins linearly spaced in frequency, whose average redshifts are ${\overline{z}}_5=0.166$, ${\overline{z}}_{15}=0.263$, and ${\overline{z}}_{25}=0.377$. The BAO wiggles appear at higher multipoles as the redshift increases. The dashed black line shows the white noise amplitude. The expected contribution from foreground residuals at the highest z-bin (${\overline{z}}_{30}=0.442$; the one with the highest amplitude of foreground residual) is represented by the gray line. The dotted lines show the effect introduced by the instrumental beam with θFWHM = 40′.

A similar procedure was used to choose the θ range for the ACF analyses. Applying the BAO fitting pipeline over the mean ω(θ) measurements from the FLASK mocks, we chose for each z-bin the range of scales that was large enough to encompass slightly more than the full width of the BAO peak. Our choices for the θ ranges extend from 10.5° ⪅ θ ⪅ 21.0° for the smallest z-bin to 2.5° ⪅ θ ⪅ 8.0° for the highest z-bin.

The θ and ℓ ranges used for each z-bin in all our analyses are presented in Table 1. The large number of z-bins demanded a large number of tests to find the more appropriate angular and mutipole ranges for each estimator. We therefore only show our final choice in each case. We use the double tilde to refer to minimum and maximum θ and ℓ values because although these numbers are exact for the fiducial values of Δθ and Δℓ, they are slightly different when testing other binning schemes. For such cases, we chose the range of scales that was as close as possible to the fiducial configuration. This is shown in Table 1. All tests for choosing the best angular and multipole ranges were applied to the BINGO-like simulations because our goal was to optimize our analyses for future BINGO observations.

6.1. Fitting clustering measurements from 21 cm only simulations

Employing the fiducial configuration, we measured the BAO signal by jointly fitting three sets of N_z = 10 consecutive z-bins, the lower, intermediate, and higher z-bins (redshift ranges of Δz ≈ 0.09, 0.11, and 0.12 widths), estimating only one α parameter for each set. The advantage of a tomographic approach, combining different z-bins, is to improve the error bars and statistical significance of the measurements. The nuisance parameters, one parameter B and four (three) parameters A_q, were fit to the C_ℓ (ω(θ)) estimated for each of the ten z-bins individually, so that for the fiducial configuration, we had a total of 10 × 5(4) = 50(40) nuisance parameters. The α fitting results obtained from analyzing the 21 cm only simulations (accounting for BINGO sky coverage and beam) are shown in the first part of Tables 2 and 3 for the ACF and APS estimators, respectively.

We started to test our BAO fitting pipeline over the high signal-to-noise ratio mean clustering measurements, that is, the mean C_ℓ and ω(θ) from the 1500 FLASK mocks. Our aim in fitting a measurement with such a high signal-to-noise ratio is not only to help decide about angular and multipole ranges and their binning widths, Δℓ and Δθ, but also to be able to identify any possible systematic bias on α estimates. Hereafter, we use α^m and ${\sigma }_{\alpha }^m$ to denote the results of fitting the mean C_ℓ or ω(θ), and α and σ_α the results for individual mocks. The fitting over the mean clustering measurement required dividing the original covariance matrix by the total number of simulations, but this scaling would not provide a realistic estimate of σ_α for one mock realization. Because likelihoods can be well approximated by Gaussians for high signal-to-noise ratio measurements of the BAO feature, we therefore followed Abbott et al. (2019) and divided the original covariance matrix by 10 to obtain the likelihood and fit the α^m parameter. The corresponding 1σ uncertainty was calculated as described in the previous section and was then normalized to obtain ${\sigma }_{\alpha }^m\to \sqrt{10}{\sigma }_{\alpha }^m$. Using a different normalization factor, for example, 20 instead of 10, leads to similar error amplitude. The ${\alpha }^m\pm {\sigma }_{\alpha }^m$ results are presented in the last column of Tables 2 and 3. Using the fiducial configuration to fit the mean clustering measurements from the 21 cm only simulations, we find a prediction of 24%, 10%, and 5% errors for the lower, intermediate, and higher z-bins, respectively. The same procedure was employed to all tests fitting the mean clustering measurements, namely, varying the fiducial configuration or/and including contaminant signals, whose results are shown in the last column of Tables 2–5 and are discussed in the next subsections. In general, α^m estimates are consistent with unity, although they present a small bias, which in all the cases is well below the statistical uncertainty.

In addition, Tables 2–5 present the same summary statistics that were obtained applying the BAO fitting pipeline to each mock for each of the three redshift ranges. In particular, there we show the average shift parameter ⟨α⟩, the average error in α estimates ⟨σ_α⟩, and two measures of the dispersion of the α distribution, the symmetric error around ⟨α⟩ encompassing 68% of the mocks (less sensitive to the tails) σ₆₈, and the common standard deviation of the distribution σ_std, as well as the average ⟨χ²⟩ (and degrees of freedom, d.o.f.). All these quantities were calculated from the fraction N_s of the total 1500 mocks, which correspond to the mocks fulfilling our selection criteria, namely, those whose α ± σ_α interval are inside the range [0.6, 1.4] and that have ${\chi }^2>{\chi }_{\text{nw}}^2$, where ${\chi }_{\text{nw}}^2$ is obtained from fitting a non-BAO template model to the clustering estimates (see Sect. 6.6). This fraction N_s is also presented in the tables, as is the fraction N_d of the mock whose α estimate is inside the range [0.8, 1.2]. We consider these as the mocks with a BAO detection. Figures 3 and 4 show for ACF and APS the distribution of the α parameter, the σ_α error, and the χ², estimated from the N_s fraction of the mocks, respectively; the blue curves show the 21 cm only results.

Fig. 3. Histogram distribution of the α parameter, σα error, and χ2 obtained by applying the BAO fitting pipeline to each of the FLASK simulations for the ACF estimator using the fiducial configuration. The rows from top to bottom show the results for each redshift range, the lower (1−10 z-bins), the intermediate (11−20 z-bins), and the higher redshifts (21−30 z-bins). The average redshift for each range is presented in the first panel of each row. In all the plots, the solid blue, dashed red, and dotted green lines show the distribution of the parameters resulting from analyzing the 21 cm only simulations (21 cm), the noisy 21 cm simulations (21 cm + WN), and the BINGO-like simulations (adding also the foreground residual; 21 cm + WN + FG), respectively. All these results correspond to the Ns fraction of the mocks. See Table 2.

The detection criterion of having the full interval α ± σ_α inside the prior range [0.8, 1.2] (e.g., Ata et al. 2018; Chan et al. 2018) has commonly been employed in the literature. However, we can observe from these analyses that the distribution of the α parameter in general seems to be more Gaussian than what we find in our analyses (left panels in Figs. 3 and 4). A possible reason for this is the lower redshift we considered here, where the contribution from nonlinear effects is stronger. However, results from Villaescusa-Navarro et al. (2017), who investigated the BAO detection from 21 cm signal for the SKA case for the redshift range 0.35 < z < 3.05, also showed α distributions with clear deviations from a perfect Gaussian (but note the smaller number of simulations employed there, namely, 100). As pointed out by Chan et al. (2018) and Abbott et al. (2022), a natural consequence from having a (approximate) Gaussian distribution is a reasonable concordance among the three different error measurements, that is, ⟨σ_α⟩∼σ₆₈ ∼ σ_std. This is not our case, as can be seen from Tables 2 and 3 (as well as from Tables 4 and 5), indicating that ⟨σ_α⟩ is not meaningful or representative of the error in the α measurements for individual mock realizations (see also Figs. 3 and 4). Our results show that compared to ⟨σ_α⟩, the errors given by σ₆₈ are overestimated by ∼18% to 33%, from the smallest to the highest z-bins when the ACF is used, but is underestimated by ∼50% for APS. Moreover, comparing the expected ${\sigma }_{\alpha }^m$ obtained fitting the mean C_ℓ and ω(θ) to the average ⟨σ_α⟩, we find a reasonable agreement for intermediate and higher z-bins, but not for the lower z-bins, in particular, for the ACF estimator. In addition, although σ_std agrees better with σ₆₈, we still find non-negligible differences among them, which confirms our α distributions as non-Gaussian (a larger σ_std indicates non-Gaussian tails). These reasons motivated our choice of using the α values instead of α ± σ_α ranges, which belong to the interval [0.8, 1.2] as a criterion for a BAO detection (defining the N_d fraction), as well as our choice of using the 68% spread of the α distributions, σ₆₈, as the representative error in our measurements.

Fig. 4. Same as Fig. 4, but for the APS estimator. All distributions are again constructed from the corresponding Ns fraction of the mocks. See Table 3.

Comparing the results obtained for the three redshift ranges, we confirm that the two clustering estimators perform better at higher redshifts, as expected. We recall that the lower the redshift, the larger the BAO angular scale, θ_BAO, which for the BINGO redshift range, 0.127 < z < 0.449, corresponds to 18° ≳θ_BAO ≳ 7° for the fiducial cosmology. This means that at the smallest redshifts, the θ_BAO is larger than the 15° stripe of the BINGO coverage, leading to worse statistics at these z-bins. Figure 5 compares APS and ACF results obtained from 21 cm only simulations (first row), showing the improvement in their performance for intermediate and higher redshifts, as well as a high correlation among their results. Comparing the performance from each of the clustering estimators, we found that they show different sensitivities to specific redshift ranges. Our results show that C_ℓ measurements lead to smaller σ₆₈ errors for intermediate and higher redshift ranges, namely, 20% and 25% smaller σ₆₈ when compared to ω(θ). Furthermore, both N_s and N_d fractions of the mocks are larger for the APS at lower and intermediate redshifts, and comparable at higher redshifts. For the N_d fraction, the APS has ∼1.8 higher probability of detecting BAO at lower redshifts at 21 cm only mocks compared to the ACF.

Fig. 5. Comparison of the α values estimated using the ACF and APS estimators (black dots). The columns from left to right show estimates from lower to higher redshift ranges, with average ⟨z⟩ = 0.172, ⟨z⟩ = 0.270, and ⟨z⟩ = 0.386. The first row shows results for the 21 cm only analyses, and the middle and bottom rows show results including white noise only and white noise + foreground residual, respectively. The red dots in the bottom row show α estimates from the reconstructed (foreground cleaned) maps.

Finally, comparing the two estimators, we find in summary an apparently better performance of the APS over the ACF. This conclusion is mainly based on the following aspects of the ACF results: the slightly larger σ₆₈ uncertainties at intermediate and higher redshifts; the quite larger bias on the average ⟨α⟩, especially at lower z-bins, corresponding to 0.95σ₆₈; and the smaller N_s and N_d fractions at small and intermediate z-bins.

We also evaluated the z-bin width and the way in which the 30 z-bins were combined to estimate the shift parameter. First, the BAO fitting pipeline was applied to three ranges of redshift composed by N_z = 12, 10, 8 consecutive z-bins with δν = 9.33 MHz, so that α is estimated over three redshift ranges of same total width each (Δz ≈ 0.10). Then, we tested a tomographic binning of the BINGO frequency band into 30 z-bins linearly spaced in redshift (δz ≈ 0.01), instead of frequency, but still estimating α for three redshift ranges of N_z = 10 z-bins each (same total redshift width). We find that in general, both cases lead to results comparable to those from Tables 2 and 3, although each redshift range and clustering estimator has its own most advantageous configuration. In this sense, for simplicity, and because previous BINGO papers employ a redshift binning linear in frequency (Liccardo et al. 2022; Fornazier et al. 2022; Costa et al. 2022), we continued to use z-bins with δν = 9.33 MHz width and the BAO fitting process considering three sets of N_z = 10 z-bin for both estimators.

6.2. Impact of including instrumental noise

Taking all specifications summarized in Sect. 4.1.2 into account, we estimated the white noise level at each pixel in the BINGO region, σ_pix, considering the same resolution as was used to generate the 21 cm simulations, N_side = 256. From this σ_pix map, we generated as many realizations of the corresponding white noise map as necessary by multiplying it by random values defined by a Gaussian distribution of zero mean and unitary standard deviation. We added a different noise realization to each of the 1500 FLASK mocks and repeated the same process of calculating C_ℓ and ω(θ) from each of them.

White noise adds a constant term N_ℓ at all multipoles to the APS term, which for BINGO dominates the 21 cm signal as a function of redshift bin for ℓ ≥ 200 − 350 (Fig. 2). Then, we debiased the C_ℓ from each noisy 21 cm simulation by subtracting the expected N_ℓ amplitude, computed from the theoretical noise level. In contrast, white noise affects the ACF over all angular scales, so that this debiasing cannot be considered. However, because the ω(θ) was calculated by averaging over δT_pδT_p′ from all pairs of pixels p and p′ separated by θ (Eq. (16)), a natural consequence is a reduction of the noise.

Our results show that the presence of thermal noise increases the σ₆₈ uncertainty by 19% and 90% for intermediate and higher z-bins when the ACF is used, while at lower z-bins, this error surprisingly diminishes by 14% (similarly to the respective bias on ⟨α⟩). The noise impact is more expressive for APS analyses, with an increase of 9%, 51%, and 122% in σ₆₈ from lower to higher z-bins. However, although APS results are more affected by the noise, it still has a smaller σ₆₈ uncertainty for higher z-bins, while for intermediate z-bins, the two estimators have comparable errors. The second row of Fig. 5 shows a comparison of the best-fit α from ω(θ) and C_ℓ estimators. It is clearly much more spread due to the presence of noise. In addition, for both estimators, the N_s and N_d fractions decrease due to the presence of thermal noise, with the greatest impact over N_s, decreasing by 10% for the ACF estimator at lower z-bins. Again surprising, N_d increases by 6% for ACF at small z-bins after including noise. In addition, the number of simulations with a BAO detection decreases more strongly the higher the redshift. In summary, the presence of thermal noise has the main impact on higher z-bins, for which the BAO feature appears at smaller angular scales (spread to larger multipoles; Fig. 2), where noise starts to dominate.

6.3. Foreground residual contamination

The last step was to take residual foreground into account, as expected for BINGO. For this we added the most important foreground signals contributing to the BINGO frequency range along with the thermal noise (Sect. 4.1.2) to a 21 cm realization. Using the GNILC code, we performed component separation to reconstruct the noisy 21 cm signal as discussed in Sect. 4.2. Moreover, from Eq. (19), using the ILC filter and the known foreground contamination, we estimated the amplitude of the foreground residual in each z-bin, Wf. This process was repeated for ten different realizations of the 21 cm signal, each of them contaminated by a different realization of the thermal noise, but all with the same foreground signal contribution. These ten different estimates of the foreground residual were used to include the expected contaminant signal to each of the 1500 mock realizations (the BINGO-like simulations). The noise realizations are still different from one mock to another, while the foreground residual was repeated every ten mocks. This avoided the necessity of applying the (computationally expensive) component separation process to all the simulations.

Given the 1500 semi-realistic BINGO-like simulations, we calculated the C_ℓ and ω(θ) from each of them and repeated the BAO fitting process using the fiducial configuration. We again debiased C_ℓ by subtracting the expected N_ℓ term. The results are shown in Tables 2 and 3, as well as in Figs. 3 and 4, for ACF and APS estimators, respectively. The last row of Fig. 5 also shows how the foreground residual can affect the concordance among the two estimators. In general, this contamination has a less significant impact, with negligible changes in the ⟨σ_α⟩, σ₆₈, and σ_std uncertainties compared to a case including noise (similar conclusions were drawn by Villaescusa-Navarro et al. 2017). Still, both N_s and N_d fractions decrease, especially at higher z-bins using ACF, and another ∼10% of the 1500 mocks no longer satisfy the selection and detection criteria. We note also that for all the redshift ranges, the best-fit α^m from the mean ACF shows a slight shift to values lower than those found when only noise was included. Because this shift persists at some redshift ranges when the fiducial configuration was changed (Table 2) and also when the N-body mocks were used (Table 4), this suggests that the foreground residual might cause this systematic bias. For the APS, on the other hand, the results do not lead to the same conclusion. In any case, even though this shift is observed from ω(θ) results after the foreground residual contribution is included, it is quite small with respect to the error estimates, that is, small enough to be of no concern.

In order to evaluate the validity of including the expected foreground residual to the mocks, instead of applying the GNILC code to each of them, we added the BAO fitting results from the ten reconstructed maps to the last row of Fig. 5. The red dots appear to agree well with the black dots. The middle and last panels do not show all ten red dots because not all the reconstructed maps provided best-fit α values in the range [0.6, 1.4] for (at least) one of the estimators.

Additionally, we tested the impact of a frequency-dependent beam size, with the BINGO θ_FWHM taking values from ∼35′ to ∼45′ from low to high z-bins. In this case, the component separation was applied after the resolution of each frequency map (21 cm, along with thermal noise and foreground signals) was converted into the common θ_FWHM ∼ 45′, the largest beam size. The foreground residual for each z-bin was estimated as before. It was added to a realization after each frequency map, containing the 21 cm signal (FLASK mocks), convolved with the corresponding beam size, and the thermal noise, was converted into the common θ_FWHM. The BAO fitting results from this set of 1500 mock realizations are presented in the last part of Tables 2 and 3. These results show that for the two clustering estimators, a frequency-dependent beam size has a negligible impact on the BAO analyses for lower and intermediate redshift ranges. At higher z-bins, the bias on ⟨α⟩ for ω(θ) and C_ℓ is slightly larger, but still quite smaller than the error estimates, and in addition, the N_s and N_d fractions for C_ℓ and ω(θ) slightly decrease and increase, respectively. We find that the foreground residuals have a larger amplitude in the higher z-bins than when the beam size is fixed (Fig. 2), which mainly contributes at larger multipoles (a deeper investigation of the impact of the frequency-dependent beam size in the foreground cleaning is necessary, but beyond the scope of this paper and will be considered in future work). Still, even with a larger contribution, the foreground residual does not seem to have a significant impact on BAO results. This allows us to keep the previous conclusions.

6.4. Robustness tests

Here we test the aspects of the fiducial configuration, that is, we consider variations of (a) the angular (Δθ) and multipole (Δℓ) binning and (b) the q index in the templates, Eqs. (15) and (16). These robustness tests were performed over the BINGO-like simulations because they are the closest to what BINGO will observe. We varied only one aspect of the fiducial configuration at a time. This aspect is listed in the first column of Tables 2 and 3, where the results of all robustness tests are summarized.

From testing the Δθ width, Table 2 shows that the different binning leads to small changes at distinct statistics. The most discrepant result concerns the larger bias on α^m obtained from the mean ω(θ) at all the redshift ranges for Δθ = 0.25°. Moreover, although Δθ = 0.40° shows slightly better results for some statistics, for example, larger N_s and N_d fractions at some redshift ranges, we find a smaller bias on ⟨α⟩ at higher z-bins for Δθ = 0.50°. In this sense, and in order to have a more diagonal covariance matrix, we chose to use Δθ = 0.50° as the fiducial configuration. Similar conclusions can be drawn from the tests on the ω(θ) template, showing that some choices of q can help to improve different statistics. However, by comparing all template tests, we note the large reduced ⟨χ²⟩ for all redshift ranges, using q = −2, 0, 1, and the large bias on ⟨α⟩ and α^m for the ACF at small z-bins for q = 0, 1. In general, cases q = 0, 1, 2 and q = −1, 0, 1, 2 present the better results, especially in terms of the fraction N_d of mocks with BAO detection. The first case produces larger N_s and N_d for lower z-bins, while the second case provides larger fractions at higher z-bins. Given the smaller bias amplitude in ⟨α⟩ at higher z-bins, we chose q = 0, 1, 2 as part of the fiducial configuration.

Results from testing Δℓ width, presented in Table 3, show a good concordance between the different uncertainties from all three cases considered, but a clearly poorer performance for Δℓ = 20. While this larger bin width gives larger N_s and N_d fractions, it also introduces large bias on ⟨α⟩ and α^m for most of the redshift ranges. When Δℓ = 10 and 15 are compared, the second case leads to slightly larger N_s and N_d fractions, but also to a larger bias on ⟨α⟩ and α^m, as well as a larger ⟨χ²⟩, for the lower z-bins. It is therefore more appropriate to use Δℓ = 10 as part of the fiducial configuration. Finally, testing the different templates to fit C_ℓ measurements, we again searched for the one with the smallest bias and uncertainty, combined with the largest N_s and N_d. In this sense, q = 0, 1, 2 and q = −1, 0, 1, 2 have the best results, but the latter gives a larger fraction of the mocks with BAO detection, N_d, and also larger N_s. We therefore chose it as part of the fiducial configuration.

From all the robustness tests, considering each redshift range individually, we might conclude that each of them would require a slightly different configuration. However, for simplicity, we chose to use the same θ/ℓ binning and ω(θ)/C_ℓ templates for all three redshift ranges, that is, the (a) and (b) aspects of the fiducial configuration, as discussed.

6.5. Tests of N-body mocks

An additional test of the BAO fitting pipeline was its application on the set of 500 N-body mocks, generated through a completely different method, as described in Sect. 4.1.1. We employed the fiducial configuration to analyze the 21 cm only simulations and the BINGO-like simulations, constructed using the N-body mocks. We used the same ten foreground residual maps as before and a fixed beam size θ_FWHM = 40 arcmin. The summary statistics from these analyses are presented in Tables 4 and 5 for ACF and APS, respectively.

From the 21 cm only analyses, we find small differences with respect to what is obtained from the FLASK mocks, more expressive for the APS estimator, especially the smaller amplitude of ⟨χ²⟩ and fractions N_s and N_d for all the redshift ranges. However, these differences can be partially explained by the number of available N-body mocks: one-third of the number of FLASK mocks. The smaller number of simulations implies, for example, a lower signal-to-noise ratio from the mean clustering measurements, providing a poorer fitting over them (note the slightly large bias on α^m from small and intermediate z-bins for both estimators). Additionally, for the APS (ACF), the number of d.o.f., increases (decreases) with the redshift, making the covariance matrix and its inverse the more biased the higher (lower) the redshift, explaining the discrepancy in the ⟨χ²⟩ and the uncertainty estimates. In this sense, we can argue that the slight differences between other summary statistics from N-body and FLASK mocks seems to be of statistical origin, although we also need to take into account that the two types of mocks were generated following very different methods.

The impact of introducing noise and foreground residuals to the N-body mocks is very similar to that obtained using the FLASK mocks. The greater impact appears in intermediate and high z-bins for the same reasons as previously discussed (see Sects. 6.2 and 6.3).

6.6. Wiggle versus no-wiggle template fitting

In order to evaluate the goodness of fit of the ω(θ) and C_ℓ templates to the corresponding measurements from each mock, we compared the minimum χ², evaluated at the best-fit α_bf value, to that obtained by using a no-BAO template, ${\chi }_{\text{nw}}^2$. To do this, we constructed a template model for each estimator (the same Eqs. (15) and (16)) using C^temp(ℓ) obtained by imposing P^temp(k) = P^nw, instead of using the parameterization from Eq. (14). The square root of the difference among the two quantities, $S=\sqrt{\mathrm{\Delta }{\chi }^2}=\sqrt{{\chi }^2-{\chi }_{\mathrm{nw}}^2({\alpha }_{\mathrm{bf}})}$, given in term of number of σ, is commonly used to evaluate the statistical significance of the BAO signal detection (see, e.g., Villaescusa-Navarro et al. 2017; Camacho et al. 2019; Abbott et al. 2022). This means that the BAO model (baryon signature) is preferred by Δχ² compared to a model without the BAO.

To evaluate the effective significance of our BAO detection over BINGO-like simulations, we considered the complete BINGO redshift range, fitting all the 30 z-bins together instead of three subsets of them. Because the number of d.o.f. increases significantly in this case, namely, 256 and 604 for ACF and APS, respectively, an unbiased estimate of the covariance matrix and its inverse requires a larger set of simulations, which directly influences the χ² estimate (see the discussion in Sect. 6.5). To do this, we rotated the original 1500 FLASK mocks six times, similarly to the procedure used to enlarge the number of N-body mocks, as described in Sect. 4.1.1, achieving a total of 10 500 mocks. The clustering statistics measured from the 10 500 BINGO-like simulations were used to estimate the covariance matrix, while the Δχ² quantity was calculated only for the 1500 original mocks. The distribution of $\sqrt{\mathrm{\Delta }{\chi }^2}$ from the corresponding N_s ∼ 90% and 94% of the mocks for ACF and APS, with average significance of ⟨S⟩ = 1.48σ and 1.61σ, respectively, are presented in Fig. 6. From the N_s fraction of mocks, we find ∼14% and ∼24% of them with $S=\sqrt{\mathrm{\Delta }{\chi }^2}>2\sigma$, for ACF and APS, respectively. Confirming the slightly better performance of the APS over the ACF estimator, the first seems to provide a BAO detection with a higher significance than the last. Moreover, although the $\sqrt{\mathrm{\Delta }{\chi }^2}$ distribution is shown only for the BINGO-like simulations, it is impacted by systematic effects similarly to other summary statistics. Compared to 21 cm only case, $\sqrt{\mathrm{\Delta }{\chi }^2}$ is significantly shifted to lower values due to the inclusion of thermal noise, while the contribution of foreground residual has a negligible impact.

Fig. 6. Histogram distribution of the significance, $S=\sqrt{\mathrm{\Delta }{\chi }^2}$, of the BAO signal, in terms of number of σ, estimated for each of the FLASK simulations (i.e., the corresponding Ns fraction of mocks) considering each clustering estimator. These results are obtained by evaluating the BINGO-like simulations combining all the 30 z-bins.

7. Conclusions

The BAO feature measured from galaxy distribution is currently recognized as one of the most important cosmological probes. Its dependence on the evolution and the components of the Universe makes it a powerful tool for investigating the physical nature of dark energy. Since its first detection, the BAO scale has been measured using different matter tracers in several redshift ranges, and now the radio 21 cm IM experiments, such as BINGO, plan to provide a complementary route for studying the BAO.

We assessed the BINGO telescope potential in its phase 1 operation to detect the transversal BAO signal. To do this, we developed a template fitting pipeline with two clustering estimators, ω(θ) and C_ℓ, to extract the BAO information from two different sets of simulations, FLASK and N-body mocks. Our analyses accounted for the BINGO sky coverage and beam effect and evaluated the impact of including thermal noise and foreground residual. The fitting procedure to measure the shift parameter, α, from the mocks was performed over three sets of consecutive z-bins using a maximum likelihood estimator. Most of our analyses were performed employing the fiducial configuration (aspects a and b defined in Sect. 6), optimized to be applied to all redshift ranges. Our main results are summarized below.

• We found that both estimators perform better at higher redshifts, but show different sensitivities to specific redshift ranges. The bias amplitude in ⟨α⟩, the σ₆₈ uncertainty, and the selection and detection fractions, N_s and N_d, suggest that the APS slightly outperforms the ACF estimator, although the ACF presents slightly smaller σ₆₈ at lower z-bins. The APS also leads to a greater average significance, ⟨S⟩, than the ACF. This behavior is observed by analyzing 21 cm only simulations and remains after contaminant signals and redshift-dependent beam size are included (see Tables 2–5).

• Accounting for the impact of including thermal noise and foreground residual one at a time, we found that the first affects the two clustering estimators more. It affects the results more significantly the higher the redshift. It enlarges the σ₆₈ error significantly, in particular for the APS estimator, for which it reaches 122% in higher z-bins. The reason for this is that the BAO wiggles appear at higher multipoles (smaller angular scales) the higher the redshift. At these scales, the thermal noise starts to dominate. The inclusion of the foreground residual, on the other hand, does not seem to have a significant impact in our results.

• Robustness tests of the fiducial configuration, considering different Δθ and Δℓ binning and alternative template models (q indexes), suggested that a more appropriate choice would be to select an optimal configuration for each of the three redshift ranges.

• In the fiducial configuration, the BAO fitting pipeline applied to 500 N-body mocks (Tables 4 and 5) showed a few aspects that differed from those in the analysis of FLASK mocks. This is caused mainly by the smaller number of N-body mocks.

• We also evaluated how preferred a BAO model is when compared to a fit of a non-BAO model, and found an average significance for the BAO detection of ∼1.48σ for ACF and ∼1.48σ for APS, combining all 30 z-bins. A fraction of 14% and 24% of the FLASK mocks for ACF and APS, respectively, provide a > 2σ BAO detection.

We emphasize that some aspects of the analysis presented here might benefit from the inclusion of more realistic choices of instrumental issues, such as the structure of the beam and instrumental effects such as 1/f noise and polarization leakage, or even aspects of the method that can be improved and deserve further study. These include a test of alternative component separation procedures, a different number of z-bins for the BAO fitting (see the results from the BINGO paper, Mericia et al. 2022), and a different fitting method from MLE.

Moreover, it is well known (Eisenstein et al. 2007b) that the nonlinear gravitational evolution of the Universe smears out the acoustic signature by inducing a damping and broadening on the BAO peak in addition to shifting its position so that the lower the redshift, the more difficult an accurate measurement of the BAO features. Because BINGO will survey a large volume of the Universe at relatively low redshifts, we also expect to improve our BAO fitting results, in particular at lower redshifts, by considering a reconstruction procedure (Obuljen et al. 2017).

In summary, we conclude that intermediate and higher redshift intervals are the most promising for measuring the BAO scale, with a probability of detection of more than ∼80% (∼70%) with the APS (ACF) estimator. These numbers are obtained using a fixed fiducial configuration and can be significantly improved by the choice of an optimal configuration for each redshift interval. Although we showed that systematic effects have a non-negligible impact on α estimates, our analyses of a semi-realistic scenario (BINGO-like mocks) confirmed that BINGO should be able to successfully detect the BAO signal at radio frequencies.

Acknowledgments

The BINGO project is supported by São Paulo Research Foundation (FAPESP) grant 2014/07885-0. C.P.N. would like to thank Edilson de Carvalho, Armando Bernui, Henrique Xavier, and Hugo Camacho for very enlightening and useful discussions. C.P.N. also acknowledges FAPESP for financial support through grant 2019/06040-0. J.Z. acknowledges support from the Ministry of Science and Technology of China (grant Nos. 2020SKA0110102). R.G.L. thanks CAPES (process 88881.162206/2017-01) and the Alexander von Humboldt Foundation for the financial support. L.S. is supported by the National Key R&D Program of China (2020YFC2201600). A.A.C. acknowledges financial support from the National Natural Science Foundation of China (grant 12175192). L.B., F.A.B, A.R.Q., and M.V.S. acknowledge PRONEX/CNPq/FAPESQ-PB (Grant no. 165/2018). E.F.: The Kavli IPMU is supported by World Premier International Research Center Initiative (WPI), MEXT, Japan. We thank an anonymous referee for her/his very insightful comments. This research made use of astropy (Astropy Collaboration 2018), healpy (Zonca et al. 2019), numpy (Van Der Walt et al. 2011), scipy (Virtanen et al. 2020) and matplotlib (Hunter 2007).

References

All Tables

All Figures

	Fig. 1. Illustrative example of a FLASK simulation at the lowest redshift bin (0.127 < z < 0.138), convolved with a Gaussian beam of θFWHM = 40 arcmin, after selecting the sky region (apodized mask) considered in all the analyses. The color indicates the brightness temperature in μK. The Mollweide projection corresponds to celestial coordinates.
In the text

Fig. 2. Theoretical APS calculated for three of the 30 tomographic bins linearly spaced in frequency, whose average redshifts are ${\overline{z}}_5=0.166$, ${\overline{z}}_{15}=0.263$, and ${\overline{z}}_{25}=0.377$. The BAO wiggles appear at higher multipoles as the redshift increases. The dashed black line shows the white noise amplitude. The expected contribution from foreground residuals at the highest z-bin (${\overline{z}}_{30}=0.442$; the one with the highest amplitude of foreground residual) is represented by the gray line. The dotted lines show the effect introduced by the instrumental beam with θFWHM = 40′.

In the text

Fig. 3. Histogram distribution of the α parameter, σα error, and χ2 obtained by applying the BAO fitting pipeline to each of the FLASK simulations for the ACF estimator using the fiducial configuration. The rows from top to bottom show the results for each redshift range, the lower (1−10 z-bins), the intermediate (11−20 z-bins), and the higher redshifts (21−30 z-bins). The average redshift for each range is presented in the first panel of each row. In all the plots, the solid blue, dashed red, and dotted green lines show the distribution of the parameters resulting from analyzing the 21 cm only simulations (21 cm), the noisy 21 cm simulations (21 cm + WN), and the BINGO-like simulations (adding also the foreground residual; 21 cm + WN + FG), respectively. All these results correspond to the Ns fraction of the mocks. See Table 2.

In the text

	Fig. 4. Same as Fig. 4, but for the APS estimator. All distributions are again constructed from the corresponding Ns fraction of the mocks. See Table 3.
In the text

Fig. 5. Comparison of the α values estimated using the ACF and APS estimators (black dots). The columns from left to right show estimates from lower to higher redshift ranges, with average ⟨z⟩ = 0.172, ⟨z⟩ = 0.270, and ⟨z⟩ = 0.386. The first row shows results for the 21 cm only analyses, and the middle and bottom rows show results including white noise only and white noise + foreground residual, respectively. The red dots in the bottom row show α estimates from the reconstructed (foreground cleaned) maps.

In the text

	Fig. 6. Histogram distribution of the significance, $S=\sqrt{\mathrm{\Delta }{\chi }^2}$, of the BAO signal, in terms of number of σ, estimated for each of the FLASK simulations (i.e., the corresponding Ns fraction of mocks) considering each clustering estimator. These results are obtained by evaluating the BINGO-like simulations combining all the 30 z-bins.
In the text