© ESO 2022

1. Introduction

Over the last century dedicated experiments mapping different large-scale observables in the Universe have increased our understanding of the cosmic history, and established modern observational cosmology as a precise and quantitative science. However, despite great successes in building a cosmological concordance model with tightly constrained parameters, a number of questions regarding the constituents of the Universe are yet to be fully answered. The nature of dark energy, which leads to the observed accelerated expansion of the Universe, is one of the great mysteries in modern cosmology. The Baryon Acoustic Oscillations from Integrated Neutral Gas Observations (BINGO) telescope, which is designed to measure one of the most powerful observables used to characterize dark energy, the Baryon Acoustic Oscillations (BAO), may enlighten this late evolution of the Universe (Abdalla et al. 2022a; Wuensche et al. 2022; Costa et al. 2022).

BINGO will map the integrated sky emission of the neutral hydrogen (H I signal) 21 cm line transition within a redshift interval of 0.127 < z < 0.449. The intensity mapping (IM) technique (Peterson et al. 2006; Battye et al. 2013) thus allows the measurement of the entire H I flux density of a wide patch of the sky at different redshift bins, producing H I maps that can be used as input data to estimate cosmological parameters (Abdalla & Rawlings 2005).

However, a crucial intermediate stage comes before the production of H I maps and the estimation of cosmological parameters. Mitigating the foreground emission in radio sky observations is critical for the reliable recovery of the 21 cm signal, which is much fainter than the diffuse emission from the Galactic interstellar medium (ISM). Several methods have been proposed in the literature to separate the astrophysical foregrounds from the cosmological 21 cm signal, with the aim of accurately reconstructing the power spectrum of the 21 cm signal without biasing the estimation of the cosmological parameters.

Most component separation methods in the literature have been devised to primarily deal with foreground contamination in cosmic microwave background (CMB) data, for which one can rely on the known blackbody frequency spectrum of the CMB to disentangle the signal from the foregrounds through multi-frequency observations. Some of these methods rely on a parametric model for the foregrounds, such as Commander (Eriksen et al. 2004, 2008), a joint CMB and foreground Bayesian fitting method based on Gibbs sampling. Other methods, the so-called blind (or semi-blind) methods, do not rely on any assumption about the frequency dependence of the foregrounds, but mostly exploit statistical correlations between frequency channels to mitigate the foreground contamination. Such blind methods include internal linear combinations (ILCs) such as the needlet ILC (NILC; Delabrouille et al. 2009; Basak & Delabrouille 2012, 2013), a constrained variance minimization implemented on spherical wavelets; the Spectral Matching Independent Component Analysis (SMICA; Delabrouille et al. 2003; Cardoso et al. 2008); and the Correlated Component Analysis (CCA; Bedini et al. 2005; Bonaldi et al. 2006), which use statistical decorrelation to disentangle independent components; the Spectral Estimation Via Expectation Maximization (SEVEM; Fernandez-Cobos et al. 2012), which builds internal foreground templates from different maps between pairs of channels; the Generalized Morphological Component Analysis (GMCA; Chapman et al. 2013), which exploits sparsity to separate CMB and foregrounds; the Independent Component Analysis (ICA), which maximizes some measure of non-Gaussianity (NG) to disentangle independent sources, such as FastICA (Maino et al. 2002); and Bayesian formulations of the ICA (Vansyngel et al. 2016).

In contrast to the CMB, the frequency dependence of the cosmological 21 cm signal is nontrivial, and the emission law somewhat random or decorrelated across frequencies, hence making it more challenging to model 21 cm emission. For this reason, component separation algorithms dedicated to 21 cm data analysis in the literature mostly reduce to foreground subtraction techniques, with a known risk of partial loss of the 21 cm signal during the subtraction. Typical component separation algorithms that have been applied to 21 cm data include principal component analysis (PCA; Alonso et al. 2015; Zuo et al. 2019), independent component analysis (ICA; Chapman et al. 2012; Wolz et al. 2014), and generalized morphological component analysis (GMCA; Chapman et al. 2013; Carucci et al. 2020). For the present analysis we use the generalized needlet ILC (GNILC; Remazeilles et al. 2011b; Olivari et al. 2016), an extension of the blind NILC method that compensates for the lack of information on the frequency dependence of the cosmological signal by some prior information on its spatial statistics (power spectrum), and also exploits here the decorrelation between cosmological signals originating from different redshift bins.

Here we present the reconstructed maps and power spectra of the cosmological 21 cm signal in the presence of various foregrounds and white noise. However, differently from the primordial fluctuations generated by inflationary models, the 21 cm brightness temperature fluctuations depend on the densities, temperatures, and velocity gradients at late cosmic times, hence giving rise to NGs. Higher order statistics are thus necessary to completely characterize the intrinsic non-Gaussian 21 cm field. However, residual foreground contamination can also leave a non-Gaussian imprint in the reconstructed 21 cm field after the component separation step. We can hence use higher order statistics to probe the non-Gaussian features of our signal at different scales and discern them from foreground residuals. The bispectrum, which is the Fourier transform of the three-point correlation function, is well suited to capture NGs, either intrinsic or left over by foreground residuals, in our recovered 21 cm field.

The map that describes the temperature of our Galaxy as a function of sky position does not behave at all like a Gaussian random field; therefore, it should have a very large nonzero bispectrum. We expect this bispectrum to be large compared to the intrinsic NGs produced by the log-normal large-scale structure (LSS) field stated above. Estimating the bispectrum should therefore be a very good test to detect any residual foregrounds in the recovered 21 cm maps.

One of our goals here is to identify spurious non-Gaussian features in the recovered 21 cm maps that would be larger than the intrinsic NG of the 21 cm signal, and which is due to any residual foreground contamination after component separation or greater than the noise present in the simulated data set. If we had such a scenario, it would be clear that residuals from foreground separation have been injected into our output maps, and we would then have a clear tool to identify if this is the case.

Alternatively, it is possible that we might find non-Gaussian features in the bispectrum that are partially due to the 21 cm signal from the above-mentioned nonlinear evolution of matter, but also due to any poorly subtracted foregrounds. In such scenarios, measurements of the bispectrum would have to model the contributions to determine if there are any residual non-Gaussian modes in the final maps that have been left over from the foreground subtraction method. We expect the bispectrum of the 21 cm to be small compared to the non-Gaussian bispectrum signal from the Galaxy, hence this scenario is unlikely.

It is important to note here that for the conditions we use to simulate the 21 cm signal in this paper the cosmological non-Gaussian signal is small, and is compatible with zero within our error bars. In these conditions, the bispectrum output is used in this paper as a double-check for the cleaning procedure and not as a tool to get cosmological information.

This paper is organized as follows. Section 2 describes the simulation codes used on this work, both for foregrounds and the 21 cm signal with non-Gaussian information. Section 3 describes our component separation procedure and debiasing approach, including residual foreground contamination. In Sect. 4 we present the bispectrum module as well as its abilities for pinpointing subtle foreground contamination. In Sect. 5 we present the conclusion of this work.

This paper is the fifth in a series of papers presenting the BINGO project. Companion Paper I is the project paper (Abdalla et al. 2022a), Paper II describes the instrument (Wuensche et al. 2022), Paper III shows the optical design (Abdalla et al. 2022b), Paper IV describes the mission simulations (Liccardo et al. 2022), Paper VI discusses the 21 cm catalog simulations (Zhang et al. 2022), and Paper VII presents the cosmological forecasts for the BINGO telescope (Costa et al. 2022).

2. Simulations

In this paper, we implement the GNILC method on simulated data sets for the BINGO telescope specifications, as given by Wuensche et al. (2022). Simulations are crucial to the analysis toolkit in order to test our entire pipeline: from observed maps to the estimation of the cosmological parameters. We use a HEALPix (Górski et al. 2005) pixelization scheme at N_side = 256, a Gaussian beam with a full width at half maximum (FWHM) of 40′, and celestial coordinates for our simulated maps.

Our simulations of the sky as observed by the BINGO telescope include the non-Gaussian 21 cm signal, several foregrounds, and white noise. The 21 cm power spectrum is generated using the Unified Cosmological Library for C_ℓs (UCLCL) code (McLeod et al. 2017; Loureiro et al. 2019) that enters as an input to the Full-sky Lognormal Astro-fields Simulation Kit (FLASK) code (Xavier et al. 2016), which in turn generates the non-Gaussian simulated 21 cm maps.

FLASK can generate fast full-sky simulations of cosmological LSS observables, such as multiple matter density tracers (galaxies, quasars, dark matter halos), CMB temperature anisotropies and weak lensing convergence, and shear fields. UCLCL is a library for computing two-point angular correlation function of various cosmological fields that are related to LSS surveys. It uses the formalism of angular power two-point correlations, and then derives the exact analytical equations for the angular power spectrum of cosmological observables. We describe below the simulated maps of 21 cm emission obtained with FLASK, of foreground components from the PSM code (Delabrouille et al. 2013) and of instrumental noise.

2.1. Cosmological signal

We simulate our non-Gaussian 21 cm map as a log-normal field for the 30 redshift bins observed by the BINGO telescope (see Table 1). The FLASK code produces non-Gaussian fields by applying a transform to an originally Gaussian field in such a way that the transformed field obeys the two-point function originally supplied to the code. This transformation also produces a modification to the one-point function (to produce a log-normal field) chosen by the user and an associated bispectrum that cannot be chosen explicitly by the user. This transformation produces a non-Gaussian signal that is meant to skew the Gaussian field into a log-normal field, especially enhancing the overdensities into higher peaks, which is what is expected by the gravitational evolution of density perturbations. A summary of FLASK characteristics and specific choices for the BINGO series of papers can be found in Paper IV.

2.2. Foregrounds

We use the PSM code (Delabrouille et al. 2013) to generate our simulated foreground maps (see Table 1 for the frequency band distribution; we have simulated bands that start slightly before and extend slightly above the nominal BINGO filter, see Wuensche et al. 2022 for the full instrument description). We consider three Galactic foreground emissions: synchrotron, free-free, and AME (assumed to be spinning dust), which we generate consistently for all the BINGO channels. The thermal dust emission is subdominant in the BINGO frequency range and is neglected in our calculations. Extragalactic foregrounds due to CMB temperature anisotropies and a background of unresolved radio point sources are also a significant source of contamination in the BINGO frequency range, and are included. Figure 1 shows the relative amplitude of all main foreground components as a function of frequency, also illustrating the average frequency scaling of those components of emission in the frequency range of interest for BINGO observations.

Fig. 1. Frequency scaling of the main foreground components simulated with the PSM. The various curves display the standard deviation of the maps of the various components as a function of frequency, at an angular resolution of 40′, and at galactic latitudes greater than 10°. The general shapes of the curves illustrate the average frequency scaling, while the relative amplitudes in the BINGO frequency range show the relative importance of the various components for the detection of 21 cm fluctuations by BINGO. The data points are from the 408 MHz map of Remazeilles et al. (2015), and the 28.4 and 44.1 GHz maps of the LFI instrument on board the Planck satellite.

Figures 2 and 3 show the maps of the components contributing to our simulated sky and their power spectra, respectively. As is evident from Fig. 3, the large foreground contribution to the sky is a challenge for component separation and for the recovery of the 21 cm signal. We need to suppress the foreground power down to the level of the expected 21 cm spectrum, which implies a foreground rejection level better than one part in 10⁴ for the synchrotron, whose power is about eight to nine orders of magnitude above that of the 21 cm signal in the sky region observed by BINGO, after masking the brightest regions around the Galactic plane.

Fig. 2. Simulated maps of the foreground components and the 21 cm cosmological signal in mK units, as observed in the BINGO frequency band 16 (see Table 1): AME (top left), CMB (top right), free-free (middle left), faint radio point sources (FRPS) (middle right), synchrotron (bottom left), and the 21 cm log-normal cosmological signal (bottom right). All maps are shown in celestial coordinates, have a HEALPix resolution of Nside = 256, and are convolved with a 40′ beam.

Fig. 3. Power spectra of the foregrounds and 21 cm signal in the masked and beam-convolved sky (in the ℓ range of 0–300) in the BINGO frequency band 15 (∼1125 MHz; see Table 1). Before computing the power spectrum, all maps were convolved with a 40 arcmin beam, and masked according to the BINGO sky coverage. We also apply a galactic mask described in Sect. 2.3.

In the next sections we briefly describe the five components mentioned earlier. For a more detailed description of the physics behind astrophysical foregrounds, see Abdalla et al. (2022a) and references therein.

2.2.1. Galactic synchrotron

Galactic synchrotron emission originates from interactions between the Galactic magnetic field and relativistic cosmic ray electrons. It is the dominant contaminant for the 21 cm signal in the frequency range covered by BINGO. The synchrotron frequency dependence is modeled as a power law, in units of antenna (Rayleigh-Jeans) temperature, T_sync(ν)∝ν^−β_s. We use as a template the 408 MHz all-sky map produced by Remazeilles et al. (2015), which is extrapolated to BINGO frequencies through a power law considering a non-uniform spectral index, β_s, over the sky (Miville-Deschenes et al. 2008).

In this simulation, we do not include any effect due to polarized synchrotron emission. We made this choice because BINGO will measure the intensity with horns that are specifically sensitive to the circular polarization that has a high attenuation to a linear polarization response. This attenuation is on the order of 40 dB, as simulated in Wuensche et al. (2022) and Abdalla et al. (2022b).

In other words, the BINGO collaboration chose to remove this contamination instrumentally, creating the horns described in Wuensche et al. (2022) and Abdalla et al. (2022b, mainly Fig. 18), whose response is polarized, and using a cross-dragonian optical system that suppresses this type of polarization.

2.2.2. Galactic free-free

Free-free emission is produced by electron-ion interaction in the Galactic ISM. It is, together with synchrotron, the main source of contamination to the 21 cm signal. The free-free emission can be traced by the Hα emission line, since both depend on the emission measure . We simulated free-free emission following Dickinson et al. (2003), using their composite template of Hα and a single spectral emission law, which is uniform over the sky, given by

(1)

where T_e is the electron temperature, T₄ = T_e/10⁴, and a(T, ν) is the Gaunt correction factor given by

This free-free emission law is a slowly varying function of frequency, depending slightly on the electron temperature, which is assumed here to be constant (T_e = 7000 K).

2.2.3. Galactic anomalous microwave emission

Dipole emission from spinning dust grains is currently the most commonly accepted explanation for the so-called dust-correlated anomalous microwave emission (AME). Radio emission is produced when the electric dipoles in small dust grains spin up. This component is mainly observed in the frequency range of 10–60 GHz, and is subdominant in the BINGO frequency range. There is, however, significant uncertainty about the exact emission spectrum.

Here we use the PSM prescription for the spinning dust, as described in Delabrouille et al. (2013), which is modeled using a spinning dust template map extrapolated in frequency using a single emission law. We use a 353GHz thermal dust template obtained using GNILC from the Planck satellite observations (Planck Collaboration Int. XLVIII 2016), which is scaled to the spinning dust emission at 22.8 GHz on the basis of the scaling found in Planck Collaboration XXV (2016). To extrapolate to lower frequencies, the spinning dust emission law is parameterized following the model of Draine & Lazarian (1998), using a mix of 96.2% warm neutral medium and 3.8% reflection nebulae. Our choice for the AME modeling differs somewhat from that in companion Paper IV, which used the Planck GNILC dust optical depth map τ₃₅₃ (Planck Collaboration Int. XLVIII 2016) instead of the GNILC dust intensity map at 353 GHz, and scaled it down from 22.8 GHz to the BINGO frequencies of ∼1 GHz using the publicly available spdust2 code instead of the Draine and Lazarian model as implemented in the PSM. The differences between the two are representative or current uncertainties. They do not matter much as for both models spinning dust emission is subdominant.

2.2.4. Cosmic microwave background

Most of the electromagnetic radiation in the Universe is in the form of a near-isotropic background of thermal radiation originating from the time when free electrons and nuclei in the primordial plasma first combined to form neutral atoms, the cosmic microwave background (CMB). The spectrum of emission is very close to that of a blackbody at an average temperature of T_CMB ≃ 2.726 K, with small temperature fluctuations across the sky, at a level on the order of 100 μK. In the BINGO frequency range, the corresponding brightness fluctuations are comparable in amplitude to those of the 21 cm signal of interest (see Fig. 2).

In our sky simulations, we generate random CMB temperature fluctuations with a harmonic power spectrum based on the best fit spectrum from Planck Collaboration VI (2020). We note, however, that CMB maps from the Planck satellite could also be used to subtract most of this component from the BINGO observations.

2.2.5. Radio point sources

In addition to emission from the diffuse interstellar medium in our own Galaxy, we must take into account emission from the background of distant compact radio sources, both Galactic and extragalactic. The most luminous of these objects can be identified as individual sources, while the rest contribute a diffuse background from the integrated emission of faint objects that cannot be detected individually.

Although not resolved by BINGO, the brightest objects are gathered in a catalog that contains sources characterized by their parameters: type of the source, position on the sky, flux density, and polarization fraction as a function of frequency. In the model we adopt the objects in this catalog are described as a population of point sources. On the other hand, the diffuse background of sources that are not detected individually is gathered in the form of sky background inhomogeneities represented using frequency-dependent brightness fluctuation maps. Although individual sources are not kept in the format of a catalog at each frequency, maps are effectively produced by summing the contribution of a large population of sources.

The population of radio sources is based on observations at 0.8, 1.4, and 4.85 GHz. The emission of each source in the frequency range of interest for BINGO is modeled as a power law, with a distribution of spectral indices for steep and flat radio sources. Details of the radio source model can be found in the description of the Planck Sky Model (Delabrouille et al. 2013).

Sources with flux densities above than the Planck 5σ detection limit in the 30, 70, 353, and 857 GHz are considered “bright” and are not included in the simulation as it is assumed that their contribution will be cut out or fitted out from the BINGO data.

2.3. BINGO sky coverage

After simulating full sky maps containing foregrounds and the 21 cm signal, it is necessary to apply a mask to the maps in order to have the BINGO sky coverage assumed for this work. Our mask excludes pixels in the ranges of latitudes that are not observed by BINGO (i.e., only the sky pixels for which −22.5 ≤ δ ≤ −7.5 are kept).

Inside the region observed by BINGO there are regions where the Galactic emission is too high to hope to detect the 21 cm emission. We generate a Galactic mask using a smoothed version of the intensity map of the total Galactic emission in one of the BINGO channels. This mask is built using the following prescription: All observed pixels are sorted by decreasing Galactic emission amplitude, after smoothing with a 10° Gaussian beam. The brightest 10% are set to 0. The next 30% are attenuated with a cosine apodization function f(x) which smoothly increases from 0 when x = 10% to 1 when x = 40%.

Figure 4 shows the apodized mask map and gnomonic views of the foreground and 21 cm cosmological signals after masking.

Fig. 4. Effects of apodization in the maps. Top: BINGO apodized mask format (5°). The zeros in the scale correspond to non-observed spots and the transition between the observed and non-observed regions is done through the apodization (pixels with values > 0 and ≤1). Bottom left: region centered (ℓ,b)=(270° , − 20° ), where the mask was applied to the 21 cm signal map (gnomonic projection). Bottom right: same configuration as bottom left for the total foreground map.

2.4. BINGO noise simulation

As outlined by the instrument paper, Paper II, this work considers the BINGO Phase 1 configuration with 28 horns, where each horn observes the sky for a long period of time at a given elevation, which corresponds to a fixed declination. We assume that it is possible to change the elevation of the horns after a certain amount of observing time, as defined in Papers III and IV. We also assume that each horn feeds two channels measuring (I + V)/2 and (I − V)/2, respectively, with a system temperature T_sys = 70 K in each channel, where I is the intensity and V is the circular polarization. The value of I is measured by summing the two outputs, and we analyze here only the I maps.

The white noise level in the I signal for a band of total width δν is

(2)

where for a correlation receiver and σ can be related to the minimum detectable flux density for the telescope, see Eq. (3) in Wuensche et al. (2022), and the s is in units of one over δν.

Assuming that the observation time for N_horns horns is uniformly spread across the BINGO survey area, and a pixelized map with total number of pixels , with N_side = 512, the total white noise level per pixel for full observing time τ_obs is

(3)

We assume that for our survey f_sky = 15%, T_sys = 70 K, N_horns = 28, and τ_obs = 1 month gives a noise per pixel at N_side = 512 (for white noise estimation), which corresponds to a noise power in each of the 30 BINGO channels of

(4)

across multipoles ℓ.

We note that in the limit where the number of data samples in the BINGO timeline is much larger than the number of pixels, each pixel value is obtained as the average of data samples uniformly spread over the pixel. The signal in each pixel is thus the integral of the signal coming from beams centered at points distributed over the pixel area. The signal is thus effectively convolved with the beam and with the pixel shape, and the power spectrum is multiplied by the square of the beam transfer function and the pixel window function W_l. The noise spectrum is directly proportional to W_l.

Two complications must be taken into account for an observation with a non-ideal instrument. The first is the presence of low frequency (1/f) noise in the timestreams. We do not consider low frequency noise in this paper and leave it for future analysis. The second is the non-uniformity of the coverage because of the layout of the BINGO horns in the focal plane. This second issue has been solved with the vertical displacement of the horns in the focal plane as described in the optics and simulation Papers III and IV.

This work uses the “double rectangular” horn arrangement described in Abdalla et al. (2022b), with 28 horns. Figure 5 shows an example noise realization for this configuration. Overall, we are able to produce reasonably homogeneous noise maps, apart from edge effects above the minimum and maximum declination covered with this horn arrangement, which should be masked in the final map analysis.

Fig. 5. Noise realizations for the “double retangular” horn arrangements, in Healpix gnomonic coordinates centered in (α = 0° ,δ = −17.5°). Left: white noise realization after multiplying the rms by a Gaussian map, using Nside = 256. Right: rms realization for the double-rectangular configuration (with time spread between five different horn offsets) for two years of observation with a 70 K system temperature at HEALpixNside = 512. (The rms maps are produced at a higher HEALpix pixelization and then degraded to the working resolution discussed in the text). Map of the corresponding rms of the projected white noise part. The color scale is saturated at five times the rms of a map with homogeneous coverage and same sky fraction.

3. Component separation analysis

Observations of the radio sky encompass a mixture of cosmological signals emitted in the early Universe (e.g., CMB and cosmological H I signal), astrophysical sources emitting in the late Universe (e.g., Galactic foregrounds and extragalactic point sources), and instrumental signals (e.g., thermal noise). Component separation is a term that refers to any data processing that tries to disentangle these emissions by exploiting correlations in observations made at separate frequencies, external constraints, and/or physical models of the different sources of emission (Delabrouille & Cardoso 2007).

Component separation techniques can be divided into two categories: parametric methods and blind (or nonparametric) methods. Parametric methods assume a spectral model for the foregrounds, while blind methods use only the observed data to recover the cosmological signal, and therefore do not make any assumption about the foregrounds. We note that in this definition of blind methods only the use of prior information about the cosmological signal is allowed.

The component separation problem is particularly relevant to IM because the observed signal is dominated by astrophysical emission, both from our Galaxy and from extragalactic sources. The removal (or the mitigation) of the astrophysical foreground contamination is a fundamental step in IM data analysis. Given that the spectral signature of the cosmological signal is nontrivial, one of the difficulties for IM component separation is the preservation of the 21 cm signal during foreground removal. Either an excess of foreground residuals in the 21 cm maps or an oversubtraction of the foregrounds from the data will result in erroneous cosmological results.

Several component separation methods have been proposed in the literature to disentangle cosmological H I 21 cm emission and astrophysical foregrounds. Many of them rely on the assumption that foregrounds are spectrally smooth. However, the smoothness assumption may be broken by instrumental effects of the telescope, such as standing waves and calibration uncertainties, although the design of the BINGO telescope has the goal of minimizing these effects.

Slight departures from spectral smoothness of the observed foregrounds may result in biases on the detection of the faint H I signal for some methods such as COMMANDER (Eriksen et al. 2004, 2008) or wp-FIT (Harker et al. 2009), given the huge dynamic range in amplitude between foregrounds and H I emission. In the case of the method GNILC (Olivari et al. 2016) that we adopt in our analysis and describe below, no assumption is made on the spectral shape of the foregrounds. To disentangle the cosmological 21 cm signal from the foreground emission, GNILC relies only on the property that the foreground emission is much more strongly correlated across frequencies than the 21 cm signal is.

3.1. GNILC methodology

In this work we use GNILC, a blind component separation method originally devised for CMB data analysis by Remazeilles et al. (2011a) and extended to 21 cm data analysis by Olivari et al. (2016). The main idea of GNILC is to use prior information on the power spectrum of the cosmological signal to either compensate for the lack of knowledge on the frequency dependence of the targeted signal, as in the case for the 21 cm signal (Olivari et al. 2016), or to overcome spectral degeneracies between components, as in the case for example between cosmic infrared background and thermal dust (Planck Collaboration Int. XLVIII 2016).

Hence, in our case, the only necessary ingredient to GNILC is a theoretical H I 21 cm power spectrum across the redshift bins as a prior to the algorithm: . We also include additional noise information in the GNILC prior, so the prior can be described by a power spectrum given by . No assumption is made about the foregrounds.

The observed data by BINGO in any direction of the sky (or pixel in the map) and at any frequency ν are a mixture of signal, foregrounds, and noise

(5)

where is the 21 cm signal at frequency ν that we aim to recover with GNILC, denotes the total foreground emission at that same frequency, and is the instrumental noise in this channel.

The formalism of GNILC has been described in detail in the literature (Remazeilles et al. 2011a; Planck Collaboration Int. XLVIII 2016; Olivari et al. 2016). We refer the reader to Olivari et al. (2016) for a full description of the method in the context of 21 cm intensity mapping. We summarize it below.

We first define a set of functions in harmonic space, (j = 1, …, 4), called needlet windows, which work as bandpass filters to handle different ranges of angular scales in the maps independently for component separation. The four needlet windows are chosen to satisfy the constraint

(6)

over the full multipole range in order to conserve the total power of the sky emission in the data processing. While choosing more than four needlet windows to partition the multipole range would allow more localized filtering in harmonic space, by the uncertainty principle the support of the filter would in contrast be less compact in pixel space and exceed the size of the small area of sky observed by BINGO. To allow sufficient localization of the filter in pixel space and capture the variations of the foreground contamination inside the BINGO stripe, we have to relax localization in harmonic space by limiting the partitioning to four needlet windows.

The spherical harmonic coefficients of each BINGO channel map d_ℓm(ν) are bandpass-filtered through each needlet window as . An inverse spherical harmonic transform of the bandpass-filtered coefficients thus provides four needlet maps at each frequency. These needlet maps contain only temperature fluctuations of typical angular scales selected by the needlet window.

In the second step, for each needlet scale (j) we compute the n × n data covariance matrix of the needlet maps in each pixel for all pairs of frequencies (ν, ν′) as

(7)

where is a domain of pixels centred around pixel , chosen in such a way to avoid artificial correlations between the signal and the foregrounds at large angular scales where the statistics is poor (the so-called ILC bias, see Delabrouille et al. 2009).

Third, using our prior estimate of the 21 cm (+ noise) power spectrum, , we simulate 21 cm (+ noise) maps for each channel ν. These simulated 21 cm (+ noise) maps receive the same needlet bandpass filtering as the data, thus leading to four needlet maps for each frequency. As in Eq. (7), for each needlet scale (j) we compute the prior signal (+ noise) covariance matrix as

(8)

We note that this prior is independent of the particular realization of the 21 cm signal in the observed sky.

In step four, as evident from Eq. (5), the data covariance matrix (Eq. (7)) must receive contributions from signal, foregrounds, and noise covariance matrices as

(9)

where we omitted the implicit pixel, frequency, and needlet-scale indices to reduce the amount of notations, and defined R_S ≡ R_{21 cm} + R_noise as the covariance matrix of 21 cm signal plus noise.

Hence, the whitened data covariance matrix, defined as , where is the prior covariance matrix (Eq. (8)), reduces to

(10)

where I is the identity matrix. The eigenvalue decomposition of matrix Eq. (10) thus yields

(11)

where D_N collects the m largest eigenvalues departing from unity, U_N the corresponding eigenvectors, and U_S the (n − m) eigenvectors whose eigenvalue is close to unity.

The decomposition in Eq. (11) enables GNILC to identify the signal and foreground subspaces in the data, since the m eigenvalues collected in the m × m diagonal matrix D_N indicate significant power from the foregrounds in the data, while the n − m eigenvalues of matrix that are close to unity correspond to power from the 21 cm signal (+noise) in the data. Hence, the m eigenvectors collected in the n × m matrix U_N form the principal components of the foreground subspace, while the n − m eigenvectors collected in the n × (n − m) matrix U_S form the independent components of the targeted 21 cm signal subspace.

In step five, unlike PCA, the effective dimension m of the foreground subspace is not pre-defined in an ad hoc manner by GNILC, but estimated directly from the data using Eq. (11) and minimizing the Akaike information criterion (AIC), which in the GNILC formalism reduces to solving (e.g., Olivari et al. 2016)

(12)

where μ_i are the eigenvalues of matrix (Eq. (10)). The foreground dimension m, which minimizes the AIC criterion, is denoted m_AIC in the next sections.

The foreground dimension m_AIC is estimated by the AIC locally across the sky and across the angular scales thanks to needlet decomposition. Hence, unlike PCA, GNILC also allows the effective dimension of the foreground subspace to vary with the needlet scale (j) and the position in the sky depending on the local signal-to-foregrounds ratio in the data (Eq. (10)). For the same reason, the n × (n − m_AIC) matrix , which spans the targeted 21 cm signal subspace, varies across the sky and the scales.

In the sixth step, for each needlet scale (j), each pixel , and each frequency ν, an estimate of the 21 cm signal is obtained by applying a multi-dimensional ILC filter to the data

(13)

where the expression for the GNILC weights is given by

(14)

and the signal mixing matrix is given by

(15)

The mixing matrix (Eq. (15)) is the only information needed to implement the GNILC filter (Eq. (14)).

The exact amplitude of the 21 cm prior is not critical for GNILC because it could be multiplied by a constant factor while leaving Eq. (14) unchanged. This is only true as long as this constant factor is not large enough to modify the dimension of the matrix U_S.

Finally, we synthesize the needlet-map estimates of the 21 cm signal to form the complete 21 cm map that includes all scales. First, we compute the spherical harmonic coefficients of the 21 cm map estimates. These coefficients are again bandpass filtered by the needlet windows (this step guarantees the normalization of the maps), and the filtered coefficients are transformed back to real space by inverse spherical harmonic transform. This operation gives one reconstructed 21 cm map per needlet scale and per frequency channel. The reconstructed 21 cm map per needlet scales are finally co-added to give the complete GNILC 21 cm map, , for each frequency channel.

To perform our analysis, we generated a cube with 30 BINGO channels and/or redshift bins for the simulation of the 21 cm signal, with the 21 cm signal having a different seed for each frequency channel. Simulated foregrounds and noise at BINGO frequency channels were then added to the corresponding redshift slices of the 21 cm cube. This cube is labeled as L₀ simulation. We also generated 100 cubes with different realizations of 21 cm signal and noise for our debiasing procedure, which is described in Sect. 3.3.

We apply GNILC to the BINGO channel maps of the L₀ simulation to extract maps of the 21 cm emission at each frequency with reduced foreground contamination, and reconstruct the 21 cm power spectrum for each redshift bin . Error bars on the reconstructed power spectra are computed analytically for each channel ν (e.g., Tristram et al. 2005)

(16)

where ℓ is the central multipole of the bin, Δℓ is the bin size, f_sky is sky fraction of the BINGO survey, and is the binned power spectrum of the GNILC 21 cm map in that channel. The uncertainty calculated by Eq. (16) thus includes cosmic variance from the 21 cm signal plus contributions from residual foregrounds and noise.

3.2. Additive and multiplicative errors

The set of reconstructed 21 cm maps by GNILC across BINGO channels contains residual contamination by foregrounds and noise

(17)

where W s^FG is the residual foreground contribution and W n is the residual noise contribution to the reconstructed 21 cm maps after component separation. These additive errors are in principle minimized by the GNILC filter Eq. (14). In addition, the GNILC filter Eq. (14) is built to ensure that Ws^{21 cm} ≃ s^{21 cm}, so that GNILC fully recovers the 21 cm signal with minimum foregrounds and noise:

(18)

However, a small part of the 21 cm signal can be removed along with the foregrounds by the filtering because the signal and foreground subspaces are not fully orthogonal in the eigenvector decomposition outlined in Eq. (11).

Hence, in practice the reconstructed 21 cm signal may suffer from a small multiplicative error or bias b (i.e., Ws^{21 cm} ≃ b s^{21 cm}) and

(19)

with b < 1. The risk of partial loss of the 21 cm signal is common to all 21 cm foreground removal techniques.

Therefore, the reconstructed power spectrum of the GNILC 21 cm map , at a given frequency ν, may have a small multiplicative error on the 21 cm signal through a suppression factor S_ℓ < 1

(20)

along with additive errors due to projected noise and residual foreground contamination .

3.3. Debiasing the power spectrum from noise bias and 21 cm signal loss

Using the BINGO specifications, we generated 100 realizations of white noise maps (1 ≤ i ≤ 100) for each channel ν. We found that using 100 simulations was enough for a suitable bias subtraction to be obtained for the purposes of this paper; however, we cannot guarantee that a covariance arising from such simulations would not have errors that are large enough for the purposes of parameter fitting. It is beyond the scope of this paper to assess the covariance errors from any lack of simulations.

We compute the projected noise realizations by applying the GNILC filter Eq. (14) of the fixed sky realization L₀ to the white noise map realizations:

(21)

We compute the power spectra of the projected noise realizations (Eq. (21)) (i.e., ) for each realization (i) and each frequency channel ν. We then average over all the N_SIM realizations to get an estimate of the projected noise power spectrum in Eq. (20):

(22)

Finally, we subtract the estimated projected noise power spectrum (Eq. (22)) to the GNILC power spectrum (Eq. (20)) as

(23)

The resulting power spectrum is thus corrected for the noise bias, so that

(24)

where in the second line of Eq. (24) we neglect the residual foreground power that is already strongly mitigated by GNILC.

We then need to correct the reconstructed 21 cm power spectrum for the multiplicative bias as

(25)

where is an estimate of the 21 cm suppression factor S_ℓ, which we compute as follows.

Using the prior on the 21 cm power spectrum, we generate 100 pure 21 cm map realizations (1 ≤ i ≤ 100) for each channel ν, which we then pass through the GNILC filter of the fixed sky realization L₀ to obtain the projected 21 cm map realizations

(26)

For each realization (i) we compute the ratio of the power spectra of the projected to the input 21 cm realizations

(27)

in other words, the suppression factors for all realizations (i) and frequency channels ν. By averaging over all realizations we then get an overall estimate of the 21 cm suppression factor for a given channel

(28)

which we use to renormalize the noise-debiased GNILC power spectrum as following the prescription in Eq. (25).

The top panels of Figs. 6–8 show, for three different BINGO channels, the reconstructed 21 cm power spectrum after foreground cleaning by GNILC and debiasing. The green bins correspond to the prior signal (+noise) power spectrum used in the analysis. We can see that for each channel the recovered 21 cm power spectrum (yellow bins) matches relatively closely the power spectrum of the input 21 cm signal (blue bins) of the L₀ simulation across multipoles. The difference between the recovered 21 cm power spectrum and the input 21 cm power spectrum is shown in the bottom right panels of each figure, highlighting an unbiased recovery of the 21 cm signal at multipoles 20 ≲ ℓ ≲ 800. In the bottom left panels of each figure we show our estimate of the suppression factor after foreground cleaning, showing a 2–6% loss (depending on multipole and channel) of the 21 cm signal before correction.

Fig. 6. Different aspects concerning the debias procedure for BINGO channel 10 (1070–1081 MHz). Top: reconstructed 21 cm power spectrum (yellow dots) for the BINGO channel 10 (1070−1081 MHz) in logarithmic scale (left) and linear scale (right), after foreground cleaning with GNILC and debiasing. Bottom right: difference between the reconstructed and input 21 cm signal power spectra. Bottom left: estimate of the suppression factor on the 21 cm signal across multipoles.

Fig. 7. Different aspects concerning the debias procedure for the BINGO channel 15 (1125–1136 MHz). Top: reconstructed 21 cm power spectrum (yellow dots) for the BINGO channel 15 (1125–1136 MHz) in logarithmic scale (left) and linear scale (right), after foreground cleaning with GNILC and debiasing. Bottom left: estimate of the suppression factor on the 21 cm signal across multipoles. Bottom right: difference between the reconstructed and input 21 cm signal power spectra.

Fig. 8. Different aspects concerning the debias procedure for the BINGO channel 20 (1180−1190 MHz). Top: reconstructed 21 cm power spectrum (yellow dots) for the BINGO channel 20 (1180–1190 MHz) in logarithmic scale (left) and linear scale (right), after foreground cleaning with GNILC and debiasing. Bottom left: estimate of the suppression factor on the 21 cm signal across multipoles. Bottom right: difference between the reconstructed and input 21 cm signal power spectra.

Our results show that our debiasing procedure applied to GNILC is quite successful in removing the additive noise bias and correcting for the small multiplicative bias on the 21 cm signal, in addition to performing efficient foreground cleaning. We have obtained similar results for the other BINGO channels. It is beyond the scope of this paper to investigate the cosmological implications of these projections and biases; however, this will be incorporated into the final BINGO pipeline so that we take into account the effects that such residual biases have on the estimation of cosmological parameters.

In the next sections we investigate the residual foreground contamination in the reconstructed 21 cm maps by looking at the power spectrum of the projected foreground components (Sect. 4.2) and the bispectrum of the reconstructed 21 cm maps (Sect. 4).

4. Bispectrum analysis as a test of residual foreground contamination

The core science aim of BINGO is to use the angular power spectrum of the fluctuation of the redshifted H I 21 cm radiation in order to measure the BAO and redshift-space distortions (Abdalla et al. 2022a). However, using higher order statistics to characterize the data can bring important information that cannot be probed with the power spectrum alone.

If the signal is Gaussian, the power spectrum contains all the information. However, NGs (i.e., deviations in our maps from Gaussian statistics) will be present in our data at all redshifts, but most notably in the low redshift 21 cm signal. They might come from the early-time evolution of the Universe (what we call primordial NG), will be imprinted onto the 21 cm maps by gravity itself (Bernardeau et al. 2002; de Putter 2018), and/or will get imprinted on the 21 cm maps when these maps are cleaned via residuals from the Galactic foreground distribution.

From the CMB anisotropy measurements (Planck Collaboration XIII 2016) we know that our Universe evolved from adiabatic initial fluctuations that are very close to Gaussian. Inflation is the most accepted paradigm for the early evolution of the Universe that generated such initial fluctuations. There are a plethora of models of inflation, and also alternatives to inflation models, that explain the evolution of the Universe at early times and that are in accordance with current cosmological observations. However, these models might predict different observables like primordial NGs (PNGs; for reviews, see Bartolo et al. 2004; Liguori et al. 2010; Chen 2010). In this way NGs represent a distinctive signature of these models, and measuring them will allow us to exclude and differentiate models and might teach us about the physics that happened at earlier times of our Universe. The presence of PNGs in the initial fluctuations that seed the density perturbations will also be imprinted in the 21 cm anisotropies measurement at late times.

However, PNGs are not the only source of NGs in the late time 21 cm H I signal. As we already discussed, measuring the cosmological 21 cm signal is a daunting task since the 21 cm H I data are dominated by foregrounds coming from cosmological and astrophysical sources. These foregrounds together with radio frequency interference might be limiting factors of 21 cm experiments, if not cleaned and mitigated correctly. Foreground cleaning techniques and RFI mitigation are some of the main techniques applied to 21 cm maps in order to recover the cosmological 21 cm H I signal. However, these techniques might introduce NG features in the residual maps, even if NGs are not present initially. In addition, any residual foreground contamination left over in the reconstructed maps of the 21 cm radiation would also create a large non-Gaussian imprint on the signal. These systematic effects are inherent in the 21 cm residual maps that will be analyzed since foreground removing is always necessary.

Given that the 21 cm H I data at late times has a NG component, we need to use higher order statistics to characterize the data and estimate these signatures. The use of the bispectrum in the 21 cm data analysis is not so unusual (see, e.g., Cunnington et al. 2021; Jolicoeur et al. 2021; Durrer et al. 2020). In Cunnington et al. (2021) there is a good description of this technique over the years. In this work we follow mainly the description found in Liguori et al. (2010), Komatsu (2001) and Smith & Zaldarriaga (2011) in order to develop the tool for identifying the NG. For the contour plot and bispectrum analysis we follow Regan & Shellard (2010).

As we mention above, the bispectrum output is used as a double-check for the cleaning procedure; it does not provide cosmological information on f_nl, g_nl (the third- and fourth-order amplitudes, respectively, of non-Gaussianity; see Komatsu 2001 for more details), or any other parameter.

In this sense, no cosmological analysis of information related to these parameters took part here. In addition, we are not arguing that the bispectrum of the 21 cm radiation is actually zero in the configurations measured, but we are arguing that in the lognormal simulations that we are producing it is close enough to zero for the assumptions made in this paper to be valid.

4.1. Angular bispectrum

We wish to compute higher angular correlation functions for the first-order brightness temperature fluctuations. We focus here on the three-point correlation function or on its Fourier transform, called the bispectrum, given that higher order correlators are usually subdominant. Since our aim is to determine the angular bispectrum, similarly to what was done for the angular power, we decompose the brightness temperature fluctuation in spherical harmonics

(29)

where the hats denote unit vectors. Given this, the angular three-point correlation function with the flat sky approximation is given by (Liguori et al. 2010; Komatsu 2001; Smith & Zaldarriaga 2011)

(30)

where is the angular bispectrum and B_ℓ1ℓ2ℓ3 is the averaged angular bispectrum given by

(31)

The matrix denotes the Wigner-3j symbol, representing the azimuthal angle dependence of the bispectrum, which is invariant under permutations. It describes three angular momenta that form a triangle L₁ + L₂ + L₃ = 0, where m₁ + m₂ + m₃ = 0, which implies that the matrix is only nonzero if the triangle conditions are satisfied: |ℓ_i − ℓ_j|≤ℓ_k ≤ ℓ_i + ℓ_j. The angular correlation function is invariant under parity, which implies that ℓ₁ + ℓ₂ + ℓ₃ is even. In this way, B_ℓ1ℓ2ℓ3 is only nonvanishing if the above triangle and parity conditions are met.

We note here that we neglect the contributions from the shot noise in the theoretical and simulated estimation of the bispectrum and the power spectrum. We assume that the contribution from the shot noise is small assuming that it arises from the number density of galaxies which emit in H I (Olivari et al. 2018). It is beyond the scope of this analysis to check if this assumption is accurate and this term can be neglected or if there are conditions under which this term can be large; however, there are simulations that show that this can in fact be larger than the value suggested by Olivari et al. (2018) and in Sect. 4 of Zhang et al. (2020).

The angular three-point correlation is invariant under rotations, thus the bispectrum can be written as

(32)

where b_ℓ1ℓ2ℓ3 is the reduced bispectrum, which is a real and symmetric function of ℓ₁, ℓ₂, and ℓ₃, and is the Gaunt integral defined as

(33)

The Gaunt integral obeys the conditions mentioned above. As all the dependencies on the Wigner-3j symbol appears only in the Gaunt integral, it is easier to study the physical properties of the bispectrum using b_l1l2l3. We can then rewrite the angular averaged bispectrum B_ℓ1ℓ2ℓ3 in terms of the reduced bispectrum b_ℓ1ℓ2ℓ3:

(34)

The bispectrum can form different triangles depending on the relation between ℓ₁,  ℓ₂,  ℓ₃. The type of triangle describes the shape of the bispectrum, with the following known shapes (Jeong & Komatsu 2009): equilateral ℓ₁ = ℓ₂ = ℓ₃, squeezed ℓ₁ ≂ ℓ₂ ≫ ℓ₃, isosceles ℓ₂ = ℓ₃, elongated ℓ₁ = ℓ₂ + ℓ₃, and folded ℓ₁ = 2ℓ₂ = 2ℓ₃. We use certain combinations of the above-defined bispectrum in later sections to show how its measurements can impose strict tests to the nature and level of the foreground residuals in the simulations and data extraction we presented at the beginning of this paper.

For the purposes of this work, different than the forms cited above, it is useful to define a specific subset of bispectrum values where ℓ₁ + ℓ₂ + ℓ₃ = N, where the sum of the ℓs equals a given factor. We call this shape “equisize”. This configuration corresponds to triangles of roughly similar size in ℓ space. We define the quantity

(35)

which is the sum of the entries of the nonzero bispectrum measurements, where the configurations obey ℓ₁ + ℓ₂ + ℓ₃ = N from the Kronecker delta. We note that the specific case where N is divisible by three, one of the above-mentioned configurations will correspond to the equilateral configuration where ℓ₁ = ℓ₂ = ℓ₃.

4.2. GNILC foreground residuals

The AIC criterion Eq. (12) is used by GNILC to determine the effective dimension of the foreground subspace (i.e., the number of principal components of matrix Eq. (10)) favored by the data among the class of models 1 ≤ m ≤ n, depending on the local signal-to-foreground ratio (Eq. (10)) across the sky and the scales. Without the AIC penalty the maximum likelihood solution would tend to overestimate the dimension of the foreground subspace (i.e., m ≲ n, where n is the number of available channels), and thus would tend to remove the 21 cm signal from the data along with the foreground contamination. The nominal setting for GNILC is to use the AIC selected value m_AIC as the default value.

In this subsection and in the next we purposefully run GNILC with sub-optimal settings for m_AIC in order to have sub-optimal foreground subtraction to quantify our ability to assess the quality of the foreground subtraction method via the bispectrum, as defined in the previous subsection. We aim to establish that we have tools available to measure and assess the systematic effects in such residual maps as obtained with the GNILC method, as described in this paper, although this approach can also be used for other component separation methods.

Because of the finite number of available frequency channels, but with a high dimensionality of both foreground and 21 cm components, component separation methods face a trade-off between two directions: removing as much foreground contamination as possible from the data without losing much 21 cm signal, or conserving as much 21 cm signal as possible at the cost of accepting more residual foreground contamination. The AIC enables GNILC to find a sweet spot in this trade-off. However, the n-dimensional space of the data is not a direct sum of the foreground and 21 cm subspaces, hence the frontier between the two subspaces is not strict.

In order to test the robustness of our results, in this and the following subsections we test GNILC with three different options for the dimension parameter of the foreground subspace: m_AIC − 1, m_AIC as selected by the AIC, and m_AIC + 1 (i.e., the AIC selected dimension increased by one). It is reasonable to assume that by increasing the number of dimensions by one or two, we would have increasingly more aggressive foreground mitigation strategies, which in turn require a larger correction factor for the 21 cm signal loss, but also a smaller projection of foregrounds in the respective 21 cm reciprocal space.

When we artificially impose one less dimension for the foreground subspace in GNILC (i.e., m_AIC − 1) the power spectrum of the recovered 21 cm map in channel 15 shows higher residual foreground contamination at all multipoles. Similarly, when we impose an extra dimension, the residuals decrease.

In contrast, as shown in Fig. 9, the more aggressive foreground strategy increases the loss of 21 cm signal across multipoles, while the default AIC selection m_AIC guarantees minimal loss of 21 cm signal. We find, however, that the signal loss on the 21 cm field decreases slowly, ranging from 2% to 10% (depending on multipoles) in the ranges of strategies that we have analyzed.

Fig. 9. Suppression factor of the 21 cm signal for different values of the foreground dimension: (mAIC − 1) (blue), mAIC (yellow), and (mAIC + 1) (green).

In addition, in Fig. 10 we present the correlation between the phases of the recovered and input 21 cm maps for the three cases studied here. Increasing the default dimension of the foreground subspace clearly reduces the scatter due to foregrounds in the correlation. However, the correlation coefficient for pure 21 cm (i.e., no noise) also degrades because of the increasing loss of signal according to Fig. 9.

Fig. 10. Phase comparison between the reconstructed 21 cm signal and the prior (21 cm signal added to white noise) for different values of the foreground dimension: mAIC − 1 (top), mAIC (middle), and mAIC + 1 (bottom). Shown here are results for channel 15.

Figure 11 shows the angular power spectrum of the residual foreground components projected in the recovered 21 cm signal for the three different options of the foreground dimension parameter in GNILC. The left column shows the results for noise-free simulations, while the right column shows the results in the presence of white noise and the use of 21 cm plus the noise prior by GNILC.

Fig. 11. Angular power spectra of the projected foreground components for the default dimension (mAIC − 1) as selected by AIC (top), (mAIC) (middle), and (mAIC + 1) (bottom) in the absence of noise (left column) or in the presence of white noise (right column) for channel 15. Here the residuals are calculated with mask and convolved with a 40 arcmin beam.

In the case of a less conservative foreground subtraction than that suggested by the AIC method (i.e., a dimension of m_AIC − 1), we see in the top left panel of Fig. 11 that there is significant residual contamination from free-free and synchrotron emission at ℓ ≲ 40. However, in the case of negligible noise, this configuration does manage to recover the majority of the power spectrum in the intermediate ℓ range 40 < ℓ < 200. The overall residual foreground contamination is negligible compared to the 21 cm signal over the larger range of multipoles 40 ≲ ℓ ≲ 500, which is of interest for BAO measurements.

In the case of m_AIC + 1 we can see in Fig. 11 that the residual foreground contributions, plotted in celestial coordinates, show little difference in comparison with m_AIC; however, the loss of the 21 cm signal in the reconstruction further increases in all channels, as highlighted in Fig. 9.

Figure 12 shows maps of the residual foreground contributions to the reconstructed 21 cm map from each individual component for the three different choices of the dimension of the foreground subspace, and for either noise-free or noisy simulations. The first two rows show larger residual foreground contamination for the m_AIC − 1 case compared to the more aggressive versions used in GNILC in the middle (m_AIC) and bottom (m_AIC + 1) rows. We can see from these results that there is a trade-off between the aggressiveness of the choice of m where aggressive versions lose signal but have lower systematics, and the AIC technique implemented in GNILC is doing a good job at determining this specific dimension.

Fig. 12. Gnomonic view of residual foreground components for dimension mAIC − 1, mAIC and mAIC + 1, with and without noise in the simulation.

4.3. BINGO pipeline: Bispectrum estimation of foreground residual contamination

We calculate the bispectrum of the reconstructed 21 cm maps as well as the foreground maps and projections using Eqs. (33) and (35). Since the bispectrum measures the non-Gaussian statistics intrinsic to the field, it is expected to return B_ℓ values that are compatible with zero for any Gaussian maps, including the noise input maps we have generated.

The reconstructed 21 cm maps by GNILC contain non-Gaussian information from both the intrinsic 21 cm fluctuations, whose distribution is log-normal given that they were simulated with FLASK, and the residual foreground contamination after component separation. The bispectrum module should therefore be able to detect non-Gaussian features in the input 21 cm maps generated by FLASK and in the recovered 21 cm obtained by GNILC, provided it has a sufficient signal-to-noise ratio.

There are several configurations we can select in order to estimate the non-Gaussian nature of a map via a bispectrum. In this analysis, as outlined in Sect. 4.1, we select configurations where ℓ₁ + ℓ₂ + ℓ₃ = N, which we call equisize configurations for the bispectrum (Regan & Shellard 2010). For the tests in this section we used a maximum value of N = 30 and N = 60 for the configurations.

In Fig. 13 we show the results for the B_ℓ values from the input of our simulations (21 cm + noise + foregrounds) in the top row of the plots, from the 21 cm maps of the FLASK simulations plus the noise in the second row of the plot, and for the simulated noise maps for the same simulation in the third row of the results. The plot shows the results for three different frequency bins in our simulation. We note that there is clearly a significant signal in the input foreground maps, which is natural and to be expected.

Fig. 13. Contour charts for the equisize configuration of the bispectrum in three different channels: 10 (first column), 15 (second column), and 20 (third column). Each line is related to each case analyzed: total foregrounds (first row), 21 cm + white noise (second row), and white noise (third row).

This signal in the bispectrum is around 16 orders of magnitude larger in the bispectrum than the signals obtained in the second and third rows. We expect the signal from the white noise configuration to have an expectation value of zero as they contain no non-Gaussianities, whereas the signal in the second row should have a non-Gaussian signal in it, although reasonably small given that this signal should come only from the log-normal transformation within FLASK, which is meant to produce a slightly skewed one-point function. This means that, in essence, this will enhance some peaks of the density field according to the log-normal distribution and will have a reasonably small amount of non-Gaussianities.

In Fig. 14 we plot the same contours for the bispectrum that are plotted in Fig. 13, but for the GNILC output maps obtained with a less aggressive value of m_AIC − 1, which we know is a configuration that is sub-optimal in terms of foreground subtraction from the previous sections, as well as the GNILC output for m_AIC, which is the preferred value chosen by the AIC algorithm. We can see that the bispectrum results for the value of m_AIC are comparable to the values obtained in the input signal. Although it is impossible from these plots alone to indicate if the GNILC run with m_AIC has significant residuals in the bispectrum, we can certainly state that if there are any, they are at most of the same level of the noise inputs. However, we can clearly assert that the level of the signal in the first row, when the code is run with m_AIC − 1, has a measurable bispectrum above the bispectrum of the noise and the bispectrum of the 21 cm signal.

Fig. 14. Contour charts for the considered configurations of the bispectrum (where ℓ1 + ℓ2 + ℓ3 = 30) in three different channels: 10 (first column), 15 (second column), and 20 (third column). Each line is related to each case analyzed: GNILC output using mAIC − 1 (first row) and GNILC output using mAIC (second row). A very similar pattern arises from the results with mAIC − 1 and the first row for the measurements with the input foregrounds in Fig. 13. These patterns are similar, albeit with a much smaller amplitude, which indicates a severe reduction of foreground levels, but not a complete removal, which is directly identified in our analysis.

This indicates that we do have a detection of foreground residuals, and hence have an independent way to check if the GNILC results are compatible with an efficient foreground subtraction. Although this check was made for GNILC, we can use such methods to test and reliably check if there are any residual foregrounds in any such maps that are produced by other methods than GNILC. We can also test optimal values for the ILC bias, which is also a parameter choice within GNILC that we have fixed to be 0.01 in this work.

In order to show this detection in a clearer way, we have also plotted Eq. (35) as a function of channel number in Fig. 15. In this figure we can see a clear detection of the residuals, which are shown to be well measured by our bispectrum test for channels below 5 and above 15.

Fig. 15. Residual values for Bℓ1 + ℓ2 + ℓ3 = N for N = 30 and N = 60, for the GNILC reconstructed maps as a function of the channel number. We plot different values of mAIC (top and bottom figures), and the residuals are a few orders of magnitude smaller than those presented in Fig. 14, independently of which mAIC is chosen for our bispectrum test for sub-optimal foreground extractions parameters.

Interestingly, the reconstruction for sub-optimal values of m_AIC are reasonable in the regions of channels from 5 to 10, which are exactly the channels where the suppression factor S_ℓ for the choice of m_AIC is close to the suppression factor for m_AIC − 1, whereas the suppression factor for the m_AIC case is closer to the m_AIC + 1 case above channel 15 and below channel 5. This clearly shows that there is an interplay between the dimension and aggressiveness of GNILC as a function of the channel number, which is dictated by the data, and is encapsulated in the measurements of the bispectrum shown in this plot.

Finally, we plot the summed bispectrum again in Fig.16 increasing the scale so that we can see all the foreground projections in our simulations. Given the smaller nature of each component, we plot the absolute value of B_{ℓ1 + ℓ2 + ℓ3 = N} and compare their magnitudes. We can see a channel-dependent structure in CMB (green curve) and AME (blue curve) foreground residuals, even though both are placed significantly below the noise level estimated by the bispectrum analysis (red dashed curve).

Fig. 16. Bℓ1 + ℓ2 + ℓ3 = N for N = 30, for the GNILC reconstructed maps as a function of the channel number. The results are plotted for different projected components from our foregrounds as well as the noise and 21 cm input field, and show that the residuals can be compared to the original values expected within the 21 cm field and noise realization.

On the other hand, the level of the synchrotron residual (black curve) is equivalent to the white noise (red dashed curve) and 21 cm (orange curve) projected components. Our tests indicate that GNILC is able to remove foregrounds including any non-Gaussian residual all the way down to the noise level injected into the data. However, we cannot conclude that these non-Gaussian residuals are removed below that noise level.

5. Conclusions

In this paper we have presented an analysis of the GNILC method to remove foregrounds in a simulated scenario of BINGO Phase 1 observations. This analysis is based upon the simulations of the foregrounds and white noise, as expected from the electronics, and does not include the 1/f noise contribution. This paper discusses in more detail the tools that were used in Paper IV, complementing the time domain analysis described in that work. It is based upon the project, instrument, and optical descriptions presented in Papers I, II and III.

We used the simulated maps described in Sect. 2 combined with white noise realizations and foreground emission to produce a realistic sky, as should be observed by BINGO. The simulated maps were also masked to remove the Galactic plane contribution, so that the subsequent analysis was performed in a sky map in which we expect low contamination from our Galaxy.

Section 3 described the component separation performed by GNILC on the simulated sky, as well as the additional steps to calibrate the residual bias left by the GNILC filter. We conclude that the recovered H I power spectrum is compatible with our input simulations within our noise levels and therefore should meet our scientific requirements.

The bispectrum module described in Sect. 4 was used to check for the presence of a non-Gaussian signal in the output H I maps, which might indicate they still contain a significant level of foreground residuals.

We also conclude that the bispectrum module is able to recognize if a non-Gaussian pattern is present in the output maps and that, for the BINGO Phase 1 configuration, we are able to reduce such residuals below the noise levels detected by the bispectrum verification. We also found that the residuals are clearly identified in the bispectrum analysis in cases where suboptimal foreground cleaning strategies are used in place of the nominal GNILC method. This can be a valuable tool for testing and verification of the quality of the foreground subtraction steps during the BINGO data analysis.

Acknowledgments

The BINGO project is supported by FAPESP grant 2014/07885-0. K.S.F.F. thanks São Paulo Research Foundation (FAPESP) for financial support through grant 2017/21570-0. Support from CNPq is gratefully acknowledged (E.A.). F.B.A. acknowledges the UKRI-FAPESP grant 2019/05687-0, and FAPESP and USP for Visiting Professor Fellowships where this work has been developed. L.S. is supported by the National Key R&D Program of China (2020YFC2201600). E.J.M. acknowledges the support by CAPES. M.R. acknowledges funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement No 725456, CMBSPEC). R.G.L. thanks CAPES (process 88881.162206/2017-01) and the Alexander von Humboldt Foundation for the financial support. C.P.N. thanks São Paulo Research Foundation (FAPESP) for financial support through grant 2019/06040-0. A.R.Q., L.B., and M.V.S. acknowledge PRONEX/CNPq/FAPESQ-PB (Grant no. 165/2018). V.L. acknowledges the postdoctoral FAPESP grant 2018/02026-0. T.V. acknowledges CNPq Grant 308876/2014-8. C.A.W. acknowledges a CNPq grant 2014/313.597. J.Z. was supported by IBS under the project code, IBS-R018-D1. A.A.C. acknowledges financial support from the China Postdoctoral Science Foundation, grant number 2020M671611. B. W. and A.A.C. were also supported by the key project of NNSFC under grant 11835009. M.P. acknowledges funding from a FAPESP Young Investigator fellowship, grant 2015/19936-1. Some of the results in this paper have been derived using the HEALPix package (Górski et al. 2005), and also Python language, including the packages healpy, astropy, numpy, and matplotlib.

References

All Tables

All Figures

Fig. 1. Frequency scaling of the main foreground components simulated with the PSM. The various curves display the standard deviation of the maps of the various components as a function of frequency, at an angular resolution of 40′, and at galactic latitudes greater than 10°. The general shapes of the curves illustrate the average frequency scaling, while the relative amplitudes in the BINGO frequency range show the relative importance of the various components for the detection of 21 cm fluctuations by BINGO. The data points are from the 408 MHz map of Remazeilles et al. (2015), and the 28.4 and 44.1 GHz maps of the LFI instrument on board the Planck satellite.

In the text

Fig. 2. Simulated maps of the foreground components and the 21 cm cosmological signal in mK units, as observed in the BINGO frequency band 16 (see Table 1): AME (top left), CMB (top right), free-free (middle left), faint radio point sources (FRPS) (middle right), synchrotron (bottom left), and the 21 cm log-normal cosmological signal (bottom right). All maps are shown in celestial coordinates, have a HEALPix resolution of Nside = 256, and are convolved with a 40′ beam.

In the text

	Fig. 3. Power spectra of the foregrounds and 21 cm signal in the masked and beam-convolved sky (in the ℓ range of 0–300) in the BINGO frequency band 15 (∼1125 MHz; see Table 1). Before computing the power spectrum, all maps were convolved with a 40 arcmin beam, and masked according to the BINGO sky coverage. We also apply a galactic mask described in Sect. 2.3.
In the text

Fig. 4. Effects of apodization in the maps. Top: BINGO apodized mask format (5°). The zeros in the scale correspond to non-observed spots and the transition between the observed and non-observed regions is done through the apodization (pixels with values > 0 and ≤1). Bottom left: region centered (ℓ,b)=(270° , − 20° ), where the mask was applied to the 21 cm signal map (gnomonic projection). Bottom right: same configuration as bottom left for the total foreground map.

In the text

Fig. 5. Noise realizations for the “double retangular” horn arrangements, in Healpix gnomonic coordinates centered in (α = 0° ,δ = −17.5°). Left: white noise realization after multiplying the rms by a Gaussian map, using Nside = 256. Right: rms realization for the double-rectangular configuration (with time spread between five different horn offsets) for two years of observation with a 70 K system temperature at HEALpixNside = 512. (The rms maps are produced at a higher HEALpix pixelization and then degraded to the working resolution discussed in the text). Map of the corresponding rms of the projected white noise part. The color scale is saturated at five times the rms of a map with homogeneous coverage and same sky fraction.

In the text

Fig. 6. Different aspects concerning the debias procedure for BINGO channel 10 (1070–1081 MHz). Top: reconstructed 21 cm power spectrum (yellow dots) for the BINGO channel 10 (1070−1081 MHz) in logarithmic scale (left) and linear scale (right), after foreground cleaning with GNILC and debiasing. Bottom right: difference between the reconstructed and input 21 cm signal power spectra. Bottom left: estimate of the suppression factor on the 21 cm signal across multipoles.

In the text

Fig. 7. Different aspects concerning the debias procedure for the BINGO channel 15 (1125–1136 MHz). Top: reconstructed 21 cm power spectrum (yellow dots) for the BINGO channel 15 (1125–1136 MHz) in logarithmic scale (left) and linear scale (right), after foreground cleaning with GNILC and debiasing. Bottom left: estimate of the suppression factor on the 21 cm signal across multipoles. Bottom right: difference between the reconstructed and input 21 cm signal power spectra.

In the text

Fig. 8. Different aspects concerning the debias procedure for the BINGO channel 20 (1180−1190 MHz). Top: reconstructed 21 cm power spectrum (yellow dots) for the BINGO channel 20 (1180–1190 MHz) in logarithmic scale (left) and linear scale (right), after foreground cleaning with GNILC and debiasing. Bottom left: estimate of the suppression factor on the 21 cm signal across multipoles. Bottom right: difference between the reconstructed and input 21 cm signal power spectra.

In the text

	Fig. 9. Suppression factor of the 21 cm signal for different values of the foreground dimension: (mAIC − 1) (blue), mAIC (yellow), and (mAIC + 1) (green).
In the text

	Fig. 10. Phase comparison between the reconstructed 21 cm signal and the prior (21 cm signal added to white noise) for different values of the foreground dimension: mAIC − 1 (top), mAIC (middle), and mAIC + 1 (bottom). Shown here are results for channel 15.
In the text

	Fig. 11. Angular power spectra of the projected foreground components for the default dimension (mAIC − 1) as selected by AIC (top), (mAIC) (middle), and (mAIC + 1) (bottom) in the absence of noise (left column) or in the presence of white noise (right column) for channel 15. Here the residuals are calculated with mask and convolved with a 40 arcmin beam.
In the text

	Fig. 12. Gnomonic view of residual foreground components for dimension mAIC − 1, mAIC and mAIC + 1, with and without noise in the simulation.
In the text

	Fig. 13. Contour charts for the equisize configuration of the bispectrum in three different channels: 10 (first column), 15 (second column), and 20 (third column). Each line is related to each case analyzed: total foregrounds (first row), 21 cm + white noise (second row), and white noise (third row).
In the text

Fig. 14. Contour charts for the considered configurations of the bispectrum (where ℓ1 + ℓ2 + ℓ3 = 30) in three different channels: 10 (first column), 15 (second column), and 20 (third column). Each line is related to each case analyzed: GNILC output using mAIC − 1 (first row) and GNILC output using mAIC (second row). A very similar pattern arises from the results with mAIC − 1 and the first row for the measurements with the input foregrounds in Fig. 13. These patterns are similar, albeit with a much smaller amplitude, which indicates a severe reduction of foreground levels, but not a complete removal, which is directly identified in our analysis.

In the text

Fig. 15. Residual values for Bℓ1 + ℓ2 + ℓ3 = N for N = 30 and N = 60, for the GNILC reconstructed maps as a function of the channel number. We plot different values of mAIC (top and bottom figures), and the residuals are a few orders of magnitude smaller than those presented in Fig. 14, independently of which mAIC is chosen for our bispectrum test for sub-optimal foreground extractions parameters.

In the text

	Fig. 16. Bℓ1 + ℓ2 + ℓ3 = N for N = 30, for the GNILC reconstructed maps as a function of the channel number. The results are plotted for different projected components from our foregrounds as well as the noise and 21 cm input field, and show that the residuals can be compared to the original values expected within the 21 cm field and noise realization.
In the text