© The Authors 2023

1. Introduction

Recently, it has been possible to study the creation and evolution of the Universe through various different cosmological experiments, shifting cosmology from mere intellectual speculation to a position that attracts the prestigious nickname of “precision cosmology”. Despite the successful endeavour of creating a cosmological model that matches most of the current observational challenges and constrains the cosmological parameter values with ∼1% precision, some aspects of the so-called “concordance model” require much more thorough investigation. Among them, the nature of the dark sector (dark matter and dark energy), which comprise about 95% of the present composition of the Universe, is among the most interesting open problems in modern astrophysics and cosmology.

The Baryon Acoustic Oscillations from Integrated Neutral Gas Observations (BINGO) telescope, a unique instrument designed to be one of the first radiotelescopes to measure baryon acoustic oscillations (BAOs) in the radio band, may unveil more details about the late evolution of the Universe (Abdalla et al. 2022a). BINGO will detect the integrated signal of the cosmological neutral hydrogen (HI signal) hyperfine transition at 1420 MHz (21 cm) frequency in the redshift interval 0.127 < z < 0.449, corresponding to a redshifted frequency interval 980 < ν < 1260 MHz. This survey will be conducted using a novel technique known as intensity mapping (IM; Peterson et al. 2006), which allows a flux measurement of the signal produced by all HI atoms over large areas of the sky. Combined with the radial dimension offered by the observational frequency band, IM can produce very large surveys covering a significant fraction of the cosmological volume. Cosmological HI datacubes can be used to probe the Universe at lower redshifts, complementing the cosmological information obtained by cosmic microwave background (CMB) temperature and polarization, weak lensing, and supernova (SN) Ia data.

Nevertheless, the very faint radio signal from the 21 cm transition is easily covered by the much stronger foreground emission, that is, mainly the Galactic diffuse component, which is constituted by synchrotron and free-free emissions, but also an extragalactic contribution from unresolved radio sources and the CMB. Mitigating the foreground contamination of the detected signal is essential to recovering the cosmological 21 cm information from the incoming data.

Various methods to separate the astrophysical foreground components from the targeted cosmological signal are proposed in the literature, and most of them have already been tested in the challenging task of CMB data analysis. A good component-separation method greatly contributes to the accurate reconstruction of the 21 cm signal. On the one hand, it works toward decreasing the foreground contamination in the data to a minimum and, on the other hand, prevents the loss of portions of the 21 cm data. An efficient component-separation method (or a combination of complementary methods) reduces the uncertainties arising during the separation process and their propagation into the 21 cm signal power spectrum, avoiding the introduction of bias in the estimation of cosmological parameters.

Some of these methods assume prior knowledge of some foreground properties. Wiener Filtering (Bunn et al. 1994; Tegmark & Efstathiou 1996), for instance, assumes prior information on the frequency dependence of the foreground emission and of the power spectra of the components of sky emission (Delabrouille & Cardoso 2009). Other approaches, such as the Gibbs sampling approaches, assume a parametric model for the emission laws of the sky components (Jewell et al. 2004; Wandelt et al. 2004; Eriksen et al. 2004, 2008; Larson et al. 2007), and fit for the parameters and the amplitude of the various components in each sky pixel. These methods are considered to be “nonblind”, because they rely on prior information about the foregrounds that may not be very well understood.

Alternatively, the so-called “blind” approaches separate components of diverse physical origin in multi-frequency observations, relying only on their statistical independence. These methods include the independent component analysis (ICA; Baccigalupi et al. 2004), which maximizes some measures of the non-Gaussianity of independent sources. Examples are FastICA, first used for foreground removal in CMB datasets (Maino et al. 2002, 2003) and later applied to HI mapping (Wolz et al. 2017); and the spectral matching independent component analysis (SMICA; Delabrouille et al. 2003; Patanchon et al. 2005; Betoule et al. 2009), which uses decorrelation to identify independent components.

The internal linear combination (ILC), and its variants (Tegmark & Efstathiou 1996; Tegmark et al. 2003; Bennett et al. 2003; Eriksen et al. 2004; Saha et al. 2006; Delabrouille et al. 2009; Basak & Delabrouille 2012, 2013; Remazeilles et al. 2011a), are blind component-separation techniques that have been extensively applied to CMB data analysis –in particular on WMAP and Planck survey data– to obtain foreground-cleaned CMB maps. This technique can be adapted to extract maps of other components, such as the 21 cm signal, using an extension called the generalized needlet ILC (Remazeilles et al. 2011b), which is used in this work to separate the 21 cm signal from the foregrounds without using prior information about the latter.

This work tests the performance of the GNILC method in recovering the HI signal in the presence of the different synchrotron models and with different binning in the frequency channels. We present reconstructed HI simulated maps and power spectra in the presence of a combination of various foregrounds and white noise.

The paper is organized as follows. In Sect. 2 we give a brief overview of the instrument. In Sect. 3, we present the cosmological signal, astrophysical foregrounds, and instrumental noise models used in our analysis, with a complementary discussion of the masking process and the description of our simulation plan. Section 4 contains a brief description of GNILC, the foreground removal method used in this work. Section 5 describes the debiasing method used to correct the reconstructed power spectra. In Sects. 6 and 7, we then present the results and the conclusions of this work, respectively.

2. Instrument overview

The BINGO telescope is under construction in Paraiba, Brazil, and will be located at coordinates 7° 2′ 27.6″ S; 38° 16′ 4.8″ W. The site-selection process is described in Peel et al. (2019). BINGO is designed to observe in the frequency interval 980 ≤ ν ≤ 1260 MHz, which corresponds to the redshift interval 0.127 < z < 0.449. BINGO will operate as a transit instrument, covering an instantaneous declination strip of ∼15° in width centered at δ = −15°. Its field of view will cover a declination strip of 15.4° in width measured from center to center, with an angular resolution of FWHM ≈ 40 arcmin at the central frequency of the band (1120 MHz). This accounts for a daily sky coverage of 5320 square degrees. In this work, we consider FWHM = 40′ for the entire BINGO frequency range.

The BINGO optics follows an off-axis crossed-Dragone configuration (Dragone 1978), with a primary 40 m diameter paraboloid and a secondary 34 m diameter hyperboloid. The secondary dish illuminates a focal surface with 28 corrugated horns, each one feeding a polarizer coupled to two magic tees. The optics and the horn design and fabrication are described, respectively, by Abdalla et al. (2022b) and Wuensche et al. (2020). BINGO will operate with a full correlation receiver per horn, connecting two radiometer chains to each magic tee. Radiometers are expected to operate at a nominal system temperature of T_sys ∼ 70 K.

Referring to the discussion of Sect. 2 of Wuensche et al. (2022), after 1 year of observations at 60% duty cycle, with 30 frequency (redshift) bins of equal bandwidth and HEALpix resolution N_side = 128 (comparable to the nominal angular resolution of the telescope), BINGO should achieve an estimated sensitivity of 102 μK.

3. Simulated data

This work builds upon and optimizes aspects of previous analyses of the potential of BINGO. We use a simulated data set similar to the one described in Fornazier et al. (2022), and extend that work using different synchrotron models, assessing the capability of clearly recovering the HI signal as a function of the number of channels used in the reconstruction. This choice allows us to assess the level of systematic differences introduced by our lack of knowledge on the foregrounds. It also introduces different samplings of the data and includes ancillary data from the C-Band All-Sky Survey (CBASS) experiment (Jones et al. 2018). We particularly find that this will greatly help our decision regarding the final binning of the data, itself impacting the hardware configuration required by BINGO.

All maps used in this work are generated with the HEALpix package (Górski et al. 2005) with N_side = 256. More details on the components included in the simulated data are found in Sect. 3.1 and on simulation sets in Sect. 3.3.

3.1. Cosmological, astrophysical, and noise components

The HI maps are generated using the Full-Sky Log-normal Astro-Fields Simulation Kit (FLASK) package (Xavier et al. 2016), using the C_ℓs created by the Unified Cosmological Library for C_ℓs (UCLCL) code (McLeod et al. 2017; Loureiro et al. 2019), and with a log-normal distribution with a very small deviation from a pure Gaussian. For more details on generating HI maps using FLASK, see Liccardo et al. (2022).

The foreground simulations are generated using the Planck Sky Model (PSM) software (Delabrouille et al. 2013). The Galactic foregrounds are synchrotron radiation, free-free emission, the anomalous microwave emission (AME), and thermal dust (TD). The extragalactic contaminants include unresolved or faint radio point sources (FRPSs) and the thermal and kinetic Sunyaev Zel’dovich (SZ) effects. Bright radio point sources are not considered, because they can be masked out during the analysis. For more details on foreground models, see Abdalla et al. (2022a), and for details of how the respective maps are generated using PSM, see Fornazier et al. (2022).

To produce the synchrotron component, we used the reprocessed Haslam 408 MHz all-sky map presented in Remazeilles et al. (2015) as a template. This map is then extrapolated to the BINGO frequencies (980–1260 MHz) through a power law given by

(1)

where is the template map defined in the reference frequency ν₀ and in the pixel p; β_s(p) is the spatially variable spectral index. In our simulations, we considered synchrotron maps produced with two models of nonuniform β_s over the sky, namely the Miville-Deschênes et al. (2008) model (hereafter, synchrotron MD), which uses WMAP data at 23 GHz, and the Giardino et al. (2002) model (hereafter, synchrotron GD), which is the result of the combination of the full-sky map of synchrotron emission at 408 from Haslam et al. (1982), the northern-hemisphere map at 1420 MHz from Reich & Reich (1986), and the southern-hemisphere map at 2326 MHz from Jonas et al. (1998).

The spectral index map of the synchrotron MD model has a mean value of −3.00 and a standard deviation of 0.06. These values for the synchrotron GD model are −2.9 and 0.1. Figure 1 shows the synchrotron emission maps generated with the two pixel-dependent spectral-index models described above, and Fig. 2 presents their respective power spectra. Polarized synchrotron emission is not included in this work, despite our understanding that polarization leakage is a critical issue in the data analysis. We are investigating this matter and intend to address the subject in a forthcoming paper.

Fig. 1. Synchrotron emission maps for the MD model (top) and the GD model (bottom), defined in the frequency channel centered at 1115 MHz and with a bandwidth of δν ∼ 9.33 MHz (30 channels configuration). The maps are presented in celestial coordinates and were convolved with a FWHM = 40 arcmin beam. The observed region corresponds to the BINGO covered area with the apodized Galactic mask defined in Sect. 3.2.

Fig. 2. Power spectra referring to the MD (green) and GD (red) synchrotron maps shown in Fig. 1, as well as the power spectrum calculated from the difference between them (blue), defined in the frequency channel centered at 1115 MHz and with a bandwidth of δν ∼ 9.33 MHz (30 channels configuration). The maps were convolved with a FWHM = 40 arcmin beam. The observed region corresponds to the area covered by BINGO with the apodized Galactic mask defined in Sect. 3.2.

Although the thermal emission from dust grains is subdominant in the BINGO frequency band, we chose to include this component in the data in order to obtain a simulated sky that more closely reflects reality. We use the Galactic dust emission maps and a corresponding model obtained from the Planck 2015 data release using the GNILC method (Planck Collaboration Int. XLVIII 2016).

Due to the 28-horn arrangement in the BINGO focal plane (double rectangular scheme, see Wuensche et al. 2022), the observation time is not uniformly distributed over the sky area covered by the instrument. In this case, each horn covers the pixels of a fixed-latitude ring. Therefore, the total observation time of a pixel, and consequently the root mean square (RMS) value of thermal noise per pixel, depend on the latitude. Strong inhomogeneities in the innermost part of the covered region are avoided by repositioning the horns every year of the mission, within a total of t_obs = 5 years, which achieves a better distribution of observation time over the pixel in the covered area (Liccardo et al. 2022). Furthermore, we assume that the noise level map is the same for all frequency channels. The RMS properties of the noise simulations used in this work follow the prescriptions published in Fornazier et al. (2022).

In terms of instrumental systematic uncertainties, we currently only consider the thermal (white) noise contribution over the sky, with an expected system temperature of 70 K. More realistic instrumental uncertainties should also be addressed in a forthcoming paper, such as real RFI emission (as measured on site) and a beam profile derived from the horn beam measurements (Wuensche et al. 2020), the optical design analysis (Abdalla et al. 2022b), and 1/f data measured directly from the receiver.

3.2. The masking process

In order to remove the region of the Galactic plane where the foregrounds are more intense, and consequently facilitate the component-separation process, we applied a Galactic mask to the simulated maps. Figure 3 shows the combined emission at 984.7 MHz of all foregrounds considered in our simulations, as well as the region of the sky covered by BINGO. The observed region is a declination strip of ≈15° centered on δ = −15°. We note that the region where the Galactic foregrounds are more intense crosses the area covered by the instrument. The brightness temperature in this intersection region reaches 56 K.

Fig. 3. Description of the masking process. Top left: map with the sum of all foregrounds considered in this work, including the synchrotron MD model (see Sect. 3.1), in the lowest frequency bin of a total of 30 channels centered on 984.7 MHz and limited to 10 K. Top right: result of applying the Galactic mask to the all-foregrounds map. Bottom left: BINGO apodized (5 deg) Galactic mask, preserving a sky fraction of 12.2%. Bottom right: result of applying the BINGO mask to the all-foregrounds map centered at 984.7 MHz. All maps are in celestial coordinates and the dashed lines delimit the BINGO covered area.

We used a Galactic mask that covers 20% of the sky, cutting off the most intense emission region of the Galactic plane. The result of applying this mask to the foregrounds map can be seen in Fig. 3. The maximum temperature within the area covered by BINGO after masking is about 6.5 K.

The BINGO observed region is defined by the feed horn arrangement in the focal plane and by the observation strategy (see Wuensche et al. 2022). The effective masked region is then the intersection between the Galactic mask discussed above and the observed region. To avoid boundary artifacts in calculating the power spectrum of the maps, we use the NaMaster¹ (Alonso et al. 2019) package to produce an apodization of type C2 and width 5° in the mask. The result is shown in the bottom right panel of Fig. 3, with a visible area of 12.2% of the sky.

3.3. Assembling the simulated foreground maps

The set of sky simulations used in this work was created to investigate the efficiency of our component-separation pipeline against simulated foregrounds with different synchrotron models for different numbers of frequency channels and with the addition of foregrounds and noise data from another experiment at a frequency outside the BINGO band.

In order to test the robustness of the method against different foreground models, we created two sets of simulated data. Each set contains a different synchrotron model (MD or GD) in addition to all other components described in the Sect. 3.1. These sets were created with the BINGO project baseline configuration of 30 frequency bins and were called MD30 and GD30.

The sky at low radio frequencies (< 10 GHz) is dominated by astrophysical foregrounds, mainly Galactic synchrotron emission. Independent foreground observations in frequencies outside the BINGO band may improve the characterization of the synchrotron contribution and other components, facilitating its removal. The C-Band All-Sky Survey (CBASS) is an all-sky survey at a frequency of 5 GHz and 1 GHz bandwidth, with a sensitivity ≲0.1mK RMS (per beam) and a resolution of 45 arcmin. CBASS was designed to provide complementary data to the CMB surveys (Jones et al. 2018). Therefore, we also perform a component-separation test, adding the simulated CBASS sky and noise to the MD30 base dataset. We refer to the result of this data combination as MD30+CBASS. The CBASS noise level is ∼437 μK for a HEALpixN_side = 256. The simulated CBASS all-sky foreground and noise emission maps are shown in Fig. 4.

Fig. 4. CBASS 5 GHz all-sky maps. Left: CBASS all-sky foregrounds map, which is a result of the sum of the components described in Table 1. Right: CBASS white-noise map. The temperature scale of the foregrounds map is saturated at ±105 μK. The temperature scale of the noise map is saturated at ±2 × 103 μK and the noise level is ∼437 μK.

3.4. Simulation plan

We generated sets of simulated data with 20, 30, 40, 60, and 80 frequency bins to evaluate the efficiency with which the 21 cm signal is reconstructed as a function of the number of channels used. To perform this test, we used the MD model for the synchrotron component, keeping the other components as described in Sect. 3.1. The created sets were called MD20, MD30, MD40, MD60, and MD80.

The BINGO input maps used in the tests described above are the result of the sum of HI and foregrounds convolved with a 40 arcmin beam plus the estimated BINGO noise. In the case of combining BINGO and CBASS data with the BINGO frequency bins, we added a channel with CBASS foregrounds plus noise. In this configuration, the maps of the cosmological and astrophysical components, both from BINGO and CBASS, are convolved with a beam of 45 arcmin. The CBASS foregrounds map contains the same components considered in the BINGO data (see Sect. 3.1). A summary of our simulation plan, defined by the foregrounds configuration, is provided in Table 1.

4. GNILC method

The GNILC method was developed to be used in CMB experiments (Remazeilles et al. 2011a) and was later adapted to the HI IM (Olivari et al. 2016). The way GNILC functions is by using not only spectral information but also spatial information (angular power spectrum) to, through a constrained principal component analysis (PCA; Murtagh & Heck 1987), discriminate between foregrounds and our targeted signal, considered here as HI plus thermal noise (see Sect. 6.1). GNILC is a blind method, that is, it does not assume prior knowledge about the properties of foregrounds, but only about the power spectrum of the desired signal. Below, we present a summary of the method, as applied to this work, and refer the reader to the original papers cited above for more details.

The sky signal at the frequency ν and pixel p can be modeled as

(2)

where s_ν(p) is the true targeted signal and f_ν(p) is the contribution of the Galactic and extragalactic foregrounds to the data. To adjust component separation to the local conditions of contamination, both in pixel and in harmonic spaces, the method uses a set of bandpass filters defined in the multipole domain ℓ (spherical wavelets), , called needlets. A set of needlets respects the relation , where j is the index associated with the window that filters the signal referring to a certain range of angular scales (or multipole range). Thus, each d_ν(p) is decomposed into maps of , where j varies from zero to the adopted number of needlets minus one, and is the result of the inverse spherical harmonic transform of , according to

(3)

where j defines the angular scales isolated by the respective needlet function. For each pixel p and each j-needlet window, a covariance matrix R^(j)(p) is calculated, whose elements for a frequency pair ν, ν′ are given by

(4)

where is the number of pixels that constitute the domain 𝔇_p centered on the pixel p and whose size is chosen so as to avoid creating artificial anti-correlations between the target signal and the foregrounds. These anti-correlations created when working with small regions of the sky can result in loss of the power of the reconstructed signal and this can be measured by the multiplicative ILC bias b, according to Delabrouille et al. (2009). The number of pixels is therefore adjusted for each j-needlet window according to the choice of the ILC bias b value, which is defined as one of the GNILC input parameters. The relation between the number of pixels in the 𝔇_p domain and the ILC bias is given by

(5)

where is the total number of pixels in the j-needlet map, n_ch is the number of frequency channels, and is the number of spherical modes in the j-needlet window.

GNILC uses an estimate of the target signal covariance matrix, , produced with map realizations from prior knowledge of the HI + thermal noise power spectrum to transform and diagonalize the data covariance matrix, R^(j)(p), disentangling foregrounds and signal subspaces on each needlet angular scale. This is the PCA step of the method. Omitting needlet scale (j) and pixel p, the diagonalization of the matrix resulting from the transformation , can be approximated by

(6)

The eigenvalues of that are ≃1 do not contain any relevant power of the foregrounds; that is, the sky emission is dominated by the target signal. The corresponding subset of eigenvectors U_s spans the target signal subspace. The subset of n_fg eigenvectors U_f for which eigenvalues of depart significantly from unity (1 + μ_i ≫ 1) spans the foregrounds subspace. For each j-needlet scale and each pixel p, the statistical Akaike information criterion (AIC) Akaike (1974) is used to select the subset of eigenvectors U_f of the covariance matrix for which eigenvalues are ≫1.

After estimating the dimension of the foregrounds subspace, a -dimensional ILC is performed for each needlet scale j. First, the matrix , whose dimension is n_ch × (n_ch−n_fg), is determined for each pixel p, and then the mixing matrix of the target signal, , is calculated through

(7)

A multidimensional ILC filter (ILC weights matrix) is applied to the data vector d^(j)(p) to make an estimate of the HI signal , according (omitting p and (j)) to

(8)

Finally, each reconstructed j-needlet map, , is transformed to spherical harmonic space and their harmonic coefficients are again band-pass-filtered by the respective j-needlet window, , and the filtered harmonic coefficients are transformed back to maps in pixel space, recreating a single map per j-needlet window, according to

(9)

These maps are then added to give, for each frequency channel ν, the full reconstructed map .

GNILC has two input parameters that control the locations used in the calculation of covariance matrices, both in the pixel and in the harmonic domains: the set of needlets and the ILC bias, and this work tested GNILC with different bias and needlets combinations. The best match for our simulations was the set of cosine-shaped needlets with peaks at ℓ = [0,128,384,767], according to Fig. 5, and the ILC bias b = 0.005, and these were the input parameters used in this work. The GNILC method was tested in this work using simulated data with different configurations and the results are presented in Sect. 6.

Fig. 5. Set of cosine-shaped needlets with peaks located at ℓ = [0, 128, 384, 767].

5. Debiasing procedure

Section 4 describes the GNILC method, which refers to the signal of interest as the component to be recovered. Due to the characteristics of the 21 cm signal and the noise, we chose to recover both as a single component. For more details on this choice, see Sect. 6.1. After obtaining the reconstructed 21 cm plus noise maps, we use a debiasing procedure to reconstruct the power spectrum of the 21 cm signal as a single component. This method consists of estimating the residual noise content and the loss of the HI signal, after passing the data through the ILC filter (see Sect. 4), and correcting their effects on the power spectra of the GNILC output maps. This procedure was also used by Fornazier et al. (2022) in the reconstruction of 21 cm power spectra from simulated data in a configuration of 30 frequency bins. Here, we explore the same method considering different configurations, as presented in Table 1 and Sect. 6.2. A brief description of the technique is presented below.

The debiasing procedure to reconstruct the 21cm power spectra is divided into two steps: (1) estimate the projected noise power spectra, , and debias the GNILC map power spectra from this additive bias; and (2) estimate the multiplicative bias, , and correct the noise-debiased GNILC power spectra from it. These two steps can be summarized by the equation

(10)

where is our final estimate of the 21cm power spectrum at the frequency channel ν.

The additive noise bias is estimated by generating N_realis white noise map realizations for each frequency channel ν and projecting them through the ILC weights matrix (see Eq. (8)) computed in Sect. 4 for the data. The power spectra of the resulting projected noise realizations are then averaged over all realizations.

The multiplicative bias is estimated by generating N_realis realizations of 21 cm signal maps at all frequency channels ν and computing the projected 21 cm signal by applying the ILC weights matrix again to the pure 21 cm map realizations. For each frequency channel ν, we then compute the ratios between the power spectra of the projected 21 cm realizations and the power spectra of the input 21 cm realizations, which we average over all realizations in order to estimate .

The accuracy of the estimation of additive and multiplicative biases depends directly on the number of realizations used. However, the choice of N_realis is not free, but limited by the available computational capacity. The greater the number of channels or realizations used, the longer the debiasing processing time. Considering our computing resources and the settings adopted for the simulated data, we chose to test the debiasing procedure with two different numbers of realizations, as can be seen in Sect. 6.2.

6. Results

Our component-separation pipeline can be divided in two steps: a foreground removal stage, where we use the GNILC method described in the Sect. 4 to recover the 21 cm signal together with the noise from the BINGO simulated data maps; and a debiasing stage, where we use the procedure described in Sect. 5 to obtain the 21 cm power spectra from the GNILC output maps. The BINGO simulated data maps used in this work are the result of the sum of astrophysical foregrounds, HI, and thermal noise, as described in Sect. 3. In the simulations that include an independent foreground observation, we add the CBASS 5 GHz map (foregrounds plus noise) to the set of BINGO simulated data maps. The results obtained in these two steps for the simulation plan in Table 1 are described in the following subsections.

6.1. Reconstruction of 21 cm plus noise maps

When trying to recover the HI signal as a single component using GNILC, it was observed that the reconstructed power spectrum was contaminated by a residual content of thermal noise, mainly at small scales (Olivari et al. 2016), which contributes to the 21 cm component estimation error. As it is possible to obtain the thermal noise characteristics with good precision in an experiment, we chose to reconstruct the HI and noise signals as a single component in a first step, as described in this section. In Sect. 6.2, we use the debiasing procedure to estimate and remove the noise content present in our power spectra in order to be able to reconstruct the pure HI signal.

First, we present the results in the pixel domain for the MD30 case, the BINGO project baseline configuration. Figure 6 shows the input (expected), the reconstructed and the residual 21 cm plus noise maps, all observed near to the central frequency of the BINGO band (1120 MHz). Figure 7 shows the total foreground emission and the respective residual map after the component separation with GNILC. We note the effect of the apodized Galactic mask on the edges of the apparent region of the maps. This result shows the effectiveness of GNILC in removing contaminants: the content of foregrounds in the data is reduced from about 10⁵ to tens of μK. Figure 8 shows a comparison between the real and reconstructed power spectra of HI plus noise, as well as the residual foregrounds after GNILC.

Fig. 6. Simulated 21 cm plus thermal noise input map (top), the GNILC reconstructed map and the respective residuals (middle), and the difference between the two (bottom), corresponding to a channel centered at 1115 MHz in the MD30 configuration. The maps are in celestial coordinates and are covered by the apodized Galactic mask defined in Sect. 3.2. The HI component was convolved with a FWHM = 40 arcmin beam.

Fig. 7. Map containing the sum of all foregrounds described in Sect. 3.1 (top) and their respective residuals after GNILC (bottom) at 1115 MHz in the MD30 configuration (see Table 1). The maps are in celestial coordinates and are covered with the apodized Galactic mask defined in Sect. 3.2. We note the region of extreme declination where the signal is attenuated by the apodized mask.

Fig. 8. Input 21 cm plus noise (black), reconstructed 21 cm plus noise (dashed green), and foregrounds residuals (dashed red) power spectra for frequency channels centered at 0.985 GHz (left), 1.115 GHz (center), and 1.255 GHz (right) for the MD30 configuration. The multipole range considered is 30 ≤ ℓ ≤ 270, with a multipole bin size of Δℓ = 12.

The power spectra are plotted in the range of multipoles 30 ≤ ℓ ≤ 270, which is equivalent to angular scales ∼6° to ∼0.7°. We do not consider very large angular scales, that is, ℓ < 30, in accordance with the limited area of the sky covered by BINGO, or small angular scales, ℓ > 270°, in line with the angular resolution of the instrument, of namely 40 arcmin. Furthermore, as BINGO covers a fraction of the sky, the power spectra of its maps suffer from a loss of angular resolution, given by Δℓ∼180° /γ_max, where γ_max is the maximum extent of the observed area (Ansari & Magneville 2010). In our case, considering the repositioning of the horns during the mission, as described in Wuensche et al. (2022) and Liccardo et al. (2022), we have γ_max = 17.5°. Therefore, in order to respect this limitation and better adapt to the range of multipoles adopted (30 ≤ ℓ ≤ 270), we choose to use a resolution of Δℓ = 12.

We begin by checking the efficiency of GNILC in reconstructing the 21 cm signal considering different models for the dominant component of the simulated data, the Galactic synchrotron emission. To this end, we performed the component separation with the MD30 and GD30 data sets (see Table 1). To quantify the reconstruction efficiency of the maps, we used the Pearson coefficient, defined as

(11)

where and s_ν(p) are the reconstructed and expected HI plus noise maps, respectively, in a given pixel p and frequency channel ν. The and μ_sν terms are the mean value over the pixels of the reconstructed and expected HI plus noise maps, respectively. The Pearson coefficients between the input and reconstructed HI plus noise maps are presented in Fig. 9. The result is seen to vary very little between the two configurations. Furthermore, the average Pearson coefficient over all channels is the same for the two cases, , which shows that GNILC appears to be robust against different synchrotron models in an analysis in the pixel domain.

Fig. 9. Left column: Pearson coefficient ρ calculated for each pair of expected and reconstructed 21 cm plus thermal noise maps. The value in parentheses is the average correlation calculated over all frequency channels. Right column: reconstruction error of the 21 cm signal power spectrum η calculated as an average over the multipole range 30 ≤ ℓ ≤ 270, on each frequency channel. The value in parentheses is the average error calculated over all frequency channels. Top row: comparison between the results of ρ and ν for the MD30 and GD30 configurations (different synchrotron models). Center row: comparison between the results of ρ and ν for the MD30 and MD30+CBASS configurations (inclusion of an independent foreground observation). Bottom row: comparison between the results of ρ and ν for the MD20, MD30, MD40, MD60, and MD80 configurations (different numbers of frequency bins).

To verify the effect of adding an extra channel to the MD30 configuration with the simulated data from the CBASS experiment (see Sect. 3.3), we compared two cases: MD30 and MD30+CBASS. The Pearson coefficients between the GNILC output maps and the expected HI plus noise maps are plotted in Fig. 9. As expected, we see an improvement in the reconstruction of HI plus noise maps when we add the CBASS channel compared to the MD30 case. The average Pearson coefficient calculated across all channels is without CBASS, while it is with the addition of the CBASS channel. Furthermore, the Pearson coefficient shows a more significant improvement in the extreme frequencies of the band.

The accuracy with which GNILC reconstructs the HI plus noise map depends mainly on the number of frequency channels adopted. Theoretically, as we increase the number of frequency bins, GNILC becomes increasingly able to remove astrophysical foregrounds from the data. This is because GNILC is able to adapt the number of degrees of freedom n_fg (see Sect. 4) dedicated to describing the foregrounds to the number of channels. This is done both in the pixel domain and in the harmonic domain, optimizing the reconstruction of the HI plus noise signal (see Olivari et al. 2016). Therefore, we tested the efficiency of the GNILC in reconstructing 21 cm plus noise maps considering 20, 30, 40, 60, and 80 channels. To this end, we used the datasets in the MD20, MD30, MD40, MD60, and MD80 configurations (see Table 1).

The Pearson coefficients for configurations with different numbers of channels are shown in Fig. 9. It is possible to observe that, in general, the Pearson coefficient values increase with an increasing number of channels. However, the points seem to tend towards a limit. The difference in reconstruction efficiency improvement seems to decrease as the number of bins increases. Between the MD60 and MD80 configurations, this difference is already very small, indicating that there is a near-optimal choice of the number of channels, which in our case appears to be around N_ch = 80. The average Pearson coefficient, calculated over all channels, is also directly related to the number of frequency bins. The calculated values range from in the MD20 configuration to in the MD80 case. As expected, an increase in the number of bins gives the method more degrees of freedom to describe and remove the foregrounds, improving the recovery of the HI plus noise signal. However, there seems to be a threshold number of channels above which there is negligible improvement in the reconstruction of the signal of interest.

Figure 9 shows that, in general, the Pearson coefficient is lower at the extreme frequencies of the band due to the edge effect, as observed by Wolz et al. (2014). This is because the adjustment for the extreme channels is restricted to only one side, while for the inner channels this is done by both sides (Harker et al. 2009). A possible solution to this is to extend the frequency coverage beyond the nominal values, as proposed by Alonso et al. (2015). Furthermore, at the lower frequencies of the band, the removal of foregrounds is also affected by the greater intensity of synchrotron emission. Finally, an increase in the Pearson coefficient is observed in relation to the extreme frequencies for the MD30+CBASS case and for the configurations with a greater number of frequency channels, such as the MD60 and MD80.

Figure 10 shows the Pearson coefficient as a function of the number of channels, calculated in channels centered at lower, central, and higher frequencies, within the BINGO band (980–1260 MHz). This figure more clearly reveals the stabilization trend in the efficiency of the reconstruction of the 21 cm plus noise maps with the increase in the number of channels. This effect is more evident in the central frequency channels, where the Pearson coefficients are already higher than in the other bins. This result reinforces the conclusion mentioned above regarding our ability to obtain a near-optimal reconstruction in the MD80 configuration.

Fig. 10. Left: Pearson coefficients calculated for each pair of expected and reconstructed HI plus thermal noise maps. Right: 21 cm power spectrum reconstruction error, calculated as an average over the multipole range 30 ≤ ℓ ≤ 270. Both parameters are calculated for simulated data with a different number of channels (MD20, MD30, MD40, MD60 and MD80 configurations), in channels centered at three different frequency channels: lowest frequencies (channels centered on 987, 984.7, 983.5, 982.3, and 981.7 MHz), central frequencies (channels centered on 1113, 1115.3, 1116.5, 1117.7, and 1118.2 MHz), and highest frequencies (channels centered on 1253, 1255.3, 1256.5, 1257.7, and 1258.2 MHz).

So far we have presented the GNILC performance in the reconstruction of the 21 cm plus noise maps. In the following section, we show the results of using the debiasing procedure described in the Sect. 5 to recover the noiseless 21 cm power spectra.

6.2. Noiseless 21 cm power spectrum reconstruction

In this section, we present the results of noiseless 21 cm signal recovery in the multipole domain using the debiasing procedure presented in Sect. 5. The power spectra of the GNILC output maps () obtained in the previous section are used here as one of the inputs for the bias-correction method (see Eq. (10)). To measure the reconstruction error of the 21 cm power spectrum, we define the index η(ν) given by the absolute mean difference between the recovered and real , normalized by the real power spectrum, according to

(12)

This index is calculated for each frequency channel ν and for a multipole range 30 ≤ ℓ ≤ 270, as defined in Sect. 6.1.

To perform the debiasing for all cases proposed in Table 1 on a reasonable timescale, we chose to work with a low number of realizations. Therefore, we initially adopted a base number N_realis = 10 and after correcting the bias for all the configurations, we performed a test increasing this amount for the configuration with the best results.

Initially, we performed the debiasing on the power spectra of the GNILC reconstructed maps for the cases with different synchrotron models. Figure 9 shows the error η of power spectrum reconstruction of the 21 cm signal in each frequency channel of the MD30 and GD30 configurations. The mean reconstruction error for the case MD30 is 7.2% and for the case GD30 is 6.4%. This difference is due to the ability of GNILC to adapt to the spatial variation of foregrounds in order to describe and remove them (see Sect. 4). Thus, maps generated with different models (MD and GD) may require different numbers of degrees of freedom to represent the contaminants in a given region of the sky and to reconstruct the target signal.

Figure 11 shows the power spectra of the foreground emission in the MD30 and GD30 configurations and their respective residuals after component separation with GNILC. At scales ℓ < 100, GNILC can be seen to remove more foregrounds from the GD30 configuration than from the MD30. Furthermore, these residuals are preserved after the debiasing step, contributing to the error in estimating the power spectrum of the 21 cm signal.

Fig. 11. Power spectra of foreground emission in MD30 (red) and GD30 (blue) configurations, as well as their respective residuals present in the signal reconstructed with the GNILC (dashed red and dashed blue). The power spectra are plotted for frequency channels centered at 0.985 GHz (left), 1.115 GHz (center), and 1.255 GHz (right). The multipole range considered is 30 ≤ ℓ ≤ 270, with a multipole bin size of Δℓ = 12.

Next, we present the results of the debiasing procedure in its basic project configuration (MD30 case), considering the addition of a channel with simulated CBASS map to the BINGO data. Figure 9 shows the error η of power spectrum reconstruction of the 21 cm signal in each frequency channel. The reconstruction error for the MD30 and MD30+CBASS cases is 7.2% and 5.7%, respectively. As expected, the inclusion of an extra channel with CBASS data increases the number of degrees of freedom available to describe the foregrounds, improving the reconstruction of the 21 cm HI signal.

Figure 9 shows the Pearson coefficient ρ calculated between the expected and reconstructed 21 cm plus thermal noise at each frequency channel, and the reconstruction error η of the 21 cm signal in the harmonic domain in each frequency channel, considering simulated data with different numbers of channels: MD20, MD30, MD40, MD60, and MD80. An improvement in the cosmological signal reconstruction can be observed with an increase in the number of channels because the number of dimensions available to describe the components increases, as discussed in Sects. 4 and 6.1. The smallest average error over all channels is , obtained with 80 frequency channels.

Figure 10 (right) shows the power spectrum reconstruction error calculated in three different channels, corresponding to the lowest, central, and highest frequencies within the BINGO band (980–1260 MHz). We note a reduction in the estimation error with an increasing number of channels, in addition to a stabilization trend of this reduction. This effect is more evident in the lowest frequency channels, where the estimation errors are already higher for smaller numbers of channels. This result corroborates the findings presented in Sect. 6.1 regarding obtaining a near-optimal reconstruction with 80 frequency bins.

To evaluate the effect of varying the number of realizations on the 21 cm spectrum reconstruction error, we reprocessed the debiasing procedure with a greater number of HI and noise realizations. To this end, we chose the MD80 configuration, the case with the best results in the previous analysis, and performed the bias correction with N_realis = 50 realizations. We then compared with the previous results presented in Fig. 12. The 21 cm signal reconstruction error for the case MD80 with N_realis = 10 is 4.7% and with N_realis = 50 is 3.0%. As expected, the reconstruction of the 21 cm signal is better when we use more HI and thermal noise realizations to estimate the additive and multiplicative bias present in the GNILC reconstructed 21 cm plus noise power spectrum. This is because with more realizations the method improves the estimation of the reconstructed spectra bias. Our objective here is not to optimize the number of realizations, which would require much more computational processing time than we had available, but to perform a sensitivity analysis of the debiasing method with this parameter.

Fig. 12. Reconstruction error η of the 21 cm signal power spectrum calculated as an average over a multipole range, defined here as 30 ≤ ℓ ≤ 270, on each frequency channel of a dataset with MD80, considering two different numbers of realizations (Nrealis = 10 and Nrealis = 50). The value in parentheses is the average error over all channels in each dataset. The mean error for the Nrealis = 10 case is and for the Nrealis = 50 case is .

To make our final estimate of the 21 cm power spectrum, we repeat the entire component separation procedure presented so far N_realis times, considering different N_realis realizations of the GNILC input maps (different realizations of HI and noise plus foregrounds). Each debiasing procedure was performed with N_realis independent realizations of HI and noise. Finally, we took the debiased power spectra and calculated the mean and the standard deviations (error bars). We carried out this entire procedure for N_realis = 10 and N_realis = 50 using the data in the MD80 configuration. Figure 13 shows the 21 cm all-sky real and estimated power spectra for two different N_realis and in three different frequency channels. It can be observed that the spectrum relative to N_realis = 50 is closer to the real power spectrum in all the frequencies presented. The best result is obtained for N_realis = 50 at 1258 MHz, where the average power spectrum reconstruction error is .

Fig. 13. Real (black) and reconstructed (green for Nrealis = 10 and red for Nrealis = 50) 21 cm signal all-sky power spectra, with respective error bars for the MD80 case (see Table 1). The power spectra plotted refer to the channels centered at 982 MHz (top), 1118 MHz (center), and 1258 MHz (bottom).

It is expected that a better estimate of the 21 cm power spectrum can be made with a larger N_realis. However, due to our available computational capacity, we limit the number of realizations to N_realis = 50. The complete set of simulations for the MD80 configuration took two and a half months to be ready. For this, we used 56 cores of the processor Intel Xeon Gold 5120 2.20 GHz and 512 GB of RAM.

7. Conclusions

This paper presents two new results relevant to the BINGO data analysis and hardware configuration before commissioning time. The component-separation procedure described here is based on GNILC, which was introduced by Liccardo et al. (2022) and Fornazier et al. (2022) as part of the BINGO collaboration, and is divided in two steps: first, we recover the HI + thermal noise maps applying GNILC to the BINGO simulated data; and second, we reconstruct the noiseless 21 cm power spectra, passing the first step results through a debiasing procedure. For our analysis, as default, we used the BINGO project baseline configuration, MD30, and ten HI and noise realizations in debiasing, that is, N_realis = 10.

Our primary findings show that GNILC is robust against different foreground-emission models. We tested two synchrotron models, finding that the reconstructed HI plus noise maps in both cases show no significant differences in the pixel domain. In the harmonic domain, after the debiasing procedure, we obtained a mean HI power spectra reconstruction error over all frequency channels in the 30 ≤ ℓ ≤ 270 interval of 7.2% for the MD30 case and 6.4% for the GD30 case. The difference between the foreground residuals of the two models is evident for ℓ < 100. This is due to the GNILC ability to adapt the foregrounds removal not only to different regions of the sky, but also to different angular scales.

We also confirm, as expected, that an extra channel with a simulated sky from the CBASS experiment improves our 21 cm power-spectrum reconstruction. Considering the MD30 configuration, the mean error in the power-spectrum reconstruction decreases from 7.2% without CBASS information to 5.7% when we include the CBASS 5 GHz channel in the simulated data (MD30+CBASS case).

Our second important finding is a near-optimal reconstruction of the HI signal in a 80-frequency-channel configuration for the BINGO experiment. In this case, we obtained a mean error in our power spectrum reconstruction of 3% in the multipole interval of 30 ≤ ℓ ≤ 270. This result was obtained by comparing the component-separation results with simulated data with different numbers of frequency bins. These foregrounds-removal runs indicate a stabilization trend in the reduction of the error on our estimation of the 21 cm signal with an increasing number of channels. More precisely, we find that the recovery quality of the HI signal begins to stabilize at N_ch = 60 to N_ch = 80, suggesting that it would not be worth using more than 100 redshift bins in the component-separation process. The definition of the optimal number of frequency (redshift) bins will have a significant impact on the hardware configuration, helping us to define the number of binning channels needed to preserve information from the raw data in < 1 ms sampling mode.

We also tested the effect of increasing the number of realizations of HI and noise used in the debiasing procedure on the quality of the 21 cm signal estimate. The mean power spectrum reconstruction error – calculated over all channels – obtained with N_realis = 10 is and with N_realis = 50 is . These results indicate, as expected, that a greater number of HI and noise realizations allows a better estimation and correction of additive and multiplicative biases of the GNILC output power spectra. For an optimal choice of N_realis, a more detailed analysis is necessary.

Finally, we repeated the entire component-separation process for 10 and 50 different realizations of the BINGO simulated input data (sky + noise) in the MD80 configuration and obtained our final estimates of the reconstructed 21 cm spectra. We compared these results with the power spectrum of an independent 21 cm full-sky realization, considered here as the real map of HI. Our estimated spectra are in good agreement with the input spectra (independent realization). As expected, with 50 realizations, it is possible to obtain a more accurate reconstruction of the HI signal than with 10 realizations. The mean reconstruction error calculated in the range 30 ≤ ℓ ≤ 270 is equal to 6.3%, 3.4%, and 2.8% for channels centered at 982 MHz, 1118 MHz, and 1258 MHz.

Our results indicate that the debiasing procedure described in this work should work efficiently in the BINGO data-analysis pipeline, suggesting that, despite recovering the signal with good efficiency in the harmonic space covered by BINGO, a more detailed analysis of the debiasing procedure should be carried out in future works, particularly because we did not include the 1/f contribution to the overall noise in this analysis.

Acknowledgments

The BINGO project is supported by FAPESP grant 2014/07885-0. E.J.M. acknowledges CAPES for the Ph.D. fellowship. C.A.W. acknowledges a CNPq grant 313597/2014-6. T.V. acknowledges CNPq Grant 308876/2014-8. L.S. is supported by the National Key R&D Program of China (2020YFC2201600) and NSFC grant 12150610459. A.A.C acknowledges financial support from the National Natural Science Foundation of China (grant 12175192). The authors thank the HEALPix creators for the HEALPix package (http://healpix.sourceforge.net; Górski et al. 2005). C.P.N. also acknowledges FAPESP for financial support through grant 2019/06040-0.

References

All Tables

All Figures

Fig. 1. Synchrotron emission maps for the MD model (top) and the GD model (bottom), defined in the frequency channel centered at 1115 MHz and with a bandwidth of δν ∼ 9.33 MHz (30 channels configuration). The maps are presented in celestial coordinates and were convolved with a FWHM = 40 arcmin beam. The observed region corresponds to the BINGO covered area with the apodized Galactic mask defined in Sect. 3.2.

In the text

Fig. 2. Power spectra referring to the MD (green) and GD (red) synchrotron maps shown in Fig. 1, as well as the power spectrum calculated from the difference between them (blue), defined in the frequency channel centered at 1115 MHz and with a bandwidth of δν ∼ 9.33 MHz (30 channels configuration). The maps were convolved with a FWHM = 40 arcmin beam. The observed region corresponds to the area covered by BINGO with the apodized Galactic mask defined in Sect. 3.2.

In the text

Fig. 3. Description of the masking process. Top left: map with the sum of all foregrounds considered in this work, including the synchrotron MD model (see Sect. 3.1), in the lowest frequency bin of a total of 30 channels centered on 984.7 MHz and limited to 10 K. Top right: result of applying the Galactic mask to the all-foregrounds map. Bottom left: BINGO apodized (5 deg) Galactic mask, preserving a sky fraction of 12.2%. Bottom right: result of applying the BINGO mask to the all-foregrounds map centered at 984.7 MHz. All maps are in celestial coordinates and the dashed lines delimit the BINGO covered area.

In the text

	Fig. 4. CBASS 5 GHz all-sky maps. Left: CBASS all-sky foregrounds map, which is a result of the sum of the components described in Table 1. Right: CBASS white-noise map. The temperature scale of the foregrounds map is saturated at ±105 μK. The temperature scale of the noise map is saturated at ±2 × 103 μK and the noise level is ∼437 μK.
In the text

	Fig. 5. Set of cosine-shaped needlets with peaks located at ℓ = [0, 128, 384, 767].
In the text

Fig. 6. Simulated 21 cm plus thermal noise input map (top), the GNILC reconstructed map and the respective residuals (middle), and the difference between the two (bottom), corresponding to a channel centered at 1115 MHz in the MD30 configuration. The maps are in celestial coordinates and are covered by the apodized Galactic mask defined in Sect. 3.2. The HI component was convolved with a FWHM = 40 arcmin beam.

In the text

	Fig. 7. Map containing the sum of all foregrounds described in Sect. 3.1 (top) and their respective residuals after GNILC (bottom) at 1115 MHz in the MD30 configuration (see Table 1). The maps are in celestial coordinates and are covered with the apodized Galactic mask defined in Sect. 3.2. We note the region of extreme declination where the signal is attenuated by the apodized mask.
In the text

	Fig. 8. Input 21 cm plus noise (black), reconstructed 21 cm plus noise (dashed green), and foregrounds residuals (dashed red) power spectra for frequency channels centered at 0.985 GHz (left), 1.115 GHz (center), and 1.255 GHz (right) for the MD30 configuration. The multipole range considered is 30 ≤ ℓ ≤ 270, with a multipole bin size of Δℓ = 12.
In the text

Fig. 9. Left column: Pearson coefficient ρ calculated for each pair of expected and reconstructed 21 cm plus thermal noise maps. The value in parentheses is the average correlation calculated over all frequency channels. Right column: reconstruction error of the 21 cm signal power spectrum η calculated as an average over the multipole range 30 ≤ ℓ ≤ 270, on each frequency channel. The value in parentheses is the average error calculated over all frequency channels. Top row: comparison between the results of ρ and ν for the MD30 and GD30 configurations (different synchrotron models). Center row: comparison between the results of ρ and ν for the MD30 and MD30+CBASS configurations (inclusion of an independent foreground observation). Bottom row: comparison between the results of ρ and ν for the MD20, MD30, MD40, MD60, and MD80 configurations (different numbers of frequency bins).

In the text

Fig. 10. Left: Pearson coefficients calculated for each pair of expected and reconstructed HI plus thermal noise maps. Right: 21 cm power spectrum reconstruction error, calculated as an average over the multipole range 30 ≤ ℓ ≤ 270. Both parameters are calculated for simulated data with a different number of channels (MD20, MD30, MD40, MD60 and MD80 configurations), in channels centered at three different frequency channels: lowest frequencies (channels centered on 987, 984.7, 983.5, 982.3, and 981.7 MHz), central frequencies (channels centered on 1113, 1115.3, 1116.5, 1117.7, and 1118.2 MHz), and highest frequencies (channels centered on 1253, 1255.3, 1256.5, 1257.7, and 1258.2 MHz).

In the text

Fig. 11. Power spectra of foreground emission in MD30 (red) and GD30 (blue) configurations, as well as their respective residuals present in the signal reconstructed with the GNILC (dashed red and dashed blue). The power spectra are plotted for frequency channels centered at 0.985 GHz (left), 1.115 GHz (center), and 1.255 GHz (right). The multipole range considered is 30 ≤ ℓ ≤ 270, with a multipole bin size of Δℓ = 12.

In the text

Fig. 12. Reconstruction error η of the 21 cm signal power spectrum calculated as an average over a multipole range, defined here as 30 ≤ ℓ ≤ 270, on each frequency channel of a dataset with MD80, considering two different numbers of realizations (Nrealis = 10 and Nrealis = 50). The value in parentheses is the average error over all channels in each dataset. The mean error for the Nrealis = 10 case is and for the Nrealis = 50 case is .

In the text

	Fig. 13. Real (black) and reconstructed (green for Nrealis = 10 and red for Nrealis = 50) 21 cm signal all-sky power spectra, with respective error bars for the MD80 case (see Table 1). The power spectra plotted refer to the channels centered at 982 MHz (top), 1118 MHz (center), and 1258 MHz (bottom).
In the text