Bayesian semi-blind component separation for
foreground removal in interferometric 21-cm
observations

Le Zhang, Emory F. Bunn, Ata Karakci, Andrei Korotkov, P.M. Sutter, Peter

T. Timbie, Gregory S. Tucker, Benjamin D. Wandelt

► To cite this version:

Le Zhang, Emory F. Bunn, Ata Karakci, Andrei Korotkov, P.M. Sutter, et al.. Bayesian semi-blind component separation for foreground removal in interferometric 21-cm observations. The Astrophysical Journal Supplement, 2016, 222, pp.3. 10.3847/0067-0049/222/1/3. hal-01554227

HAL Id: hal-01554227

https://hal.science/hal-01554227v1

Submitted on 13 Nov 2025

**HAL** is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
BAYESIAN SEMI-BLIND COMPONENT SEPARATION FOR FOREGROUND REMOVAL
IN INTERFEROMETRIC 21 cm OBSERVATIONS

LE ZHANG<fnref n="1" />, EMORY F. BUNN<fnref n="2" />, ATA KARAKCI<fnref n="3" />, ANDREI KOROTKOV<fnref n="3" />, P. M. SUTTER<fnref n="4" />,<fnref n="5" />,<fnref n="6" />, PETER T. TIMBIE<fnref n="1" />,
GREGORY S. TUCKER<fnref n="3" />, AND BENJAMIN D. WANDELT<fnref n="7" />,<fnref n="8" />,<fnref n="9" />,<fnref n="10" />

Received 2015 May 19; accepted 2015 November 20; published 2016 January 8

ABSTRACT

In this paper, we present a new Bayesian semi-blind approach for foreground removal in observations of the 21 cm signal measured by interferometers. The technique, which we call H1 Expectation-Maximization Independent Component Analysis (HIEMICA), is an extension of the Independent Component Analysis technique developed for two-dimensional (2D) cosmic microwave background maps to three-dimensional (3D) 21 cm cosmological signals measured by interferometers. This technique provides a fully Bayesian inference of power spectra and maps and separates the foregrounds from the signal based on the diversity of their power spectra. Relying only on the statistical independence of the components, this approach can jointly estimate the 3D power spectrum of the 21 cm signal, as well as the 2D angular power spectrum and the frequency dependence of each foreground component, without any prior assumptions about the foregrounds. This approach has been tested extensively by applying it to mock data from interferometric 21 cm intensity mapping observations under idealized assumptions of instrumental effects. We also discuss the impact when the noise properties are not known completely. As a first step toward solving the 21 cm power spectrum analysis problem, we compare the semi-blind HIEMICA technique to the commonly used Principal Component Analysis. Under the same idealized circumstances, the proposed technique provides significantly improved recovery of the power spectrum. This technique can be applied in a straightforward manner to all 21 cm interferometric observations, including epoch of reionization measurements, and can be extended to single-dish observations as well.

Key words: cosmology: observations – diffuse radiation – instrumentation: interferometers – methods: data analysis – methods: statistical

1. INTRODUCTION

One of the main goals of modern cosmology is to measure
the large-scale structure (LSS) of the universe, which encodes
crucial information about cosmological processes. Measuring
these large-scale fluctuations is the primary tool for under-
standing cosmological and astrophysical phenomena and
determining the related parameters. Observations (Smoot
et al. 1992; Hinshaw et al. 2013; Ade et al. 2014) of the
cosmic microwave background (CMB) precisely probe fluctua-
tions induced by primordial matter density perturbations from
the last scattering surface at a redshift of z ~ 1100 to constrain
models of inflation (Peiris et al. 2003; Kinney et al. 2006).
Spectroscopic galaxy redshift surveys such as the 2dF Galaxy
Redshift Survey (Colless et al. 2003), the WiggleZ Dark
Energy Survey (Drinkwater et al. 2010), the Sloan Digital Sky
Survey (York et al. 2000), LAMOST (Wang et al. 2009), and
the Baryon Oscillation Spectroscopic Survey (Anderson
et al. 2013) have attempted to map out the three-dimensional
(3D) structure of galaxies in the universe, providing an even
more precise measurement of cosmological parameters.

Neutral hydrogen (H1) tomography may provide a powerful
alternative to galaxy surveys. The redshifted 21 cm hyperfine
transition line allows for direct measurement of the redshift

(distance) of a source. Using H1 tomography through red-
shifted 21 cm emission could be a relatively inexpensive tool to
probe LSS over enormous volumes of our universe, from
redshift 0 to ultimately as high as ~200, potentially collecting
orders of magnitude more information than either CMB or
galaxy surveys can provide. Cosmological parameters could be
determined with far greater precision (Mao et al. 2008), and the
properties of dark matter, dark energy, neutrinos (Loeb &
Wyithe 2008; Morales & Wyithe 2010; Pritchard & Loeb
2012), and non-Gaussianity (Maldacena 2003; Loeb & Zal-
darriaga 2004; Xu et al. 2015) from inflation could be tightly
constrained. Twenty-one centimeter emission may be the only
way to explore the dark ages and the epoch of reionization
(Madau et al. 1997; Chen & Miralda-Escude 2004; Furlanetto
et al. 2006; McQuinn et al. 2006).

In particular, it is possible to trace the dynamics of the
expansion of the universe and study the evolution of the dark
energy by mapping the large-scale structure of the universe in
three dimensions (Seo & Eisenstein 2003). Several authors
(Peterson et al. 2006; Chang et al. 2008; Loeb & Wyithe 2008;
Wyithe et al. 2008; Morales & Wyithe 2010) have proposed a
new technique, the so-called “21 cm intensity mapping,” for
the observation of baryon acoustic oscillations (BAOs) in the

<footnote n="1">¹Department of Physics, University of Wisconsin, Madison, WI 53706, USA; lzhang263@wisc.edu</footnote>

<footnote n="2">²Physics Department, University of Richmond, Richmond, VA 23173, USA</footnote>

<footnote n="3">³Department of Physics, Brown University, 182 Hope Street, Providence, RI 02912, USA</footnote>

<footnote n="4">⁴Center for Cosmology and Astro-Particle Physics, Ohio State University, Columbus, OH 43210, USA</footnote>

<footnote n="5">⁵INFN—National Institute for Nuclear Physics, via Valerio 2, I-34127 Trieste, Italy</footnote>

<footnote n="6">⁶INAF—Osservatorio Astronomico di Trieste, via Tiepolo 11, I-34143 Trieste, Italy</footnote>

<footnote n="7">⁷Department of Physics, University of Illinois at Urbana-Champaign, 1110 W Green Street, Urbana, IL 61801, USA</footnote>

<footnote n="8">⁸UPMC Univ Paris 06, UMR 7095, Institut d'Astrophysique de Paris, 98 bis, boulevard Arago, F-75014 Paris, France</footnote>

<footnote n="9">⁹CNRS, UMR 7095, Institut d'Astrophysique de Paris, 98 bis, boulevard Arago, F-75014 Paris, France</footnote>

<footnote n="10">¹⁰Department of Astronomy, University of Illinois at Urbana-Champaign, Urbana, IL 61801, USA</footnote>
large-scale structure. By mapping the 21 cm intensity field with an angular resolution of ~10 arcmin across the redshift range from 0 to 3, it may be possible to determine precisely the dark energy equation of state (Bull et al. 2015). This technique has been tested with existing telescopes such as the Green Bank Telescope (Chang et al. 2010; Masui et al. 2013; Switzer et al. 2013) and the Parkes Radio Telescope (Pen et al. 2009). However, to obtain the sensitivity and angular resolution necessary to resolve BAO peaks, dedicated radio telescopes with large numbers of receivers and apertures of the order of 100 m are required. Dedicated single-aperture telescopes (e.g., FAST, Nan et al. 2011) are prohibitively expensive. Instead, interferometers provide an economical way to observe large volumes of the sky and efficiently measure the BAO features. Additionally, interferometry is well-suited to EoR and BAO measurements: visibility naturally measures the Fourier modes of the intensity distribution on the sky, which is directly related to the power spectrum of statistical fluctuations in the H1 distribution. Radio interferometric arrays are already being designed and constructed. Currently, there are several dedicated H1 interferometers that plan to detect the redshifted 21 cm EoR signal: LOFAR<fnref n="11" />, MWA<fnref n="12" />, 21 cm A<fnref n="13" />, and PAPER.<fnref n="14" /> GMRT<fnref n="15" /> is also currently attempting to detect the 21 cm signal from reionization. Furthermore, several 21 cm BAO experiments to conduct low-redshift surveys have been proposed (including BAOBAB, Pober et al. 2013a; BAORadio, Ansari et al. 2012b) or constructed (CHIME<fnref n="16" /> and TianLai, Chen 2012). Next-generation experiments such as SKA<fnref n="17" /> will enable revolutionary progress in 21 cm tomography.

However, H1 tomography is extremely challenging because astrophysical foreground contamination is expected to be at least four orders of magnitude brighter than the 21 cm signal. The success of the statistical detection of the signal relies on robust foreground removal, which has been extensively explored over the past decade. Different methods have different advantages and biases in providing signal estimates. Fortunately, foregrounds are expected to have smooth continuum spectra that vary slowly with frequency, while the 21 cm signal is expected to oscillate dramatically over a short frequency range and to become uncorrelated beyond frequency separations of about 1 MHz (Santos et al. 2005). In principle, one can remove slowly varying components to extract the H1 signal from the contaminants, although this process will always leave some residual contamination. The early ideas use a multi-frequency angular cross-correlation power spectrum to separate foreground from the 21 cm signal (Di Matteo et al. 2002, 2004; Oh & Mack 2003; Cooray & Furlanetto 2004; Zaldarriaga et al. 2004; Santos et al. 2005). An alternative approach is to subtract foreground along each line of sight in the sky or in the uv plane by fitting the frequency dependence with a particular smooth function (e.g., a polynomial) and then using the residuals to determine the 21 cm power spectrum (Santos et al. 2005; Bowman et al. 2006, 2009; McQuinn et al. 2006; Wang et al. 2006; Gleser et al. 2008; Jelić et al. 2008; Harker et al. 2009a, 2010; Liu et al. 2009a; Petrovic & Oh 2011).

Higher-order statistics, like skewness, are also employed to clean foregrounds (Harker et al. 2009b). In addition, Liu et al. (2009b), Trott et al. (2012), and Morales et al. (2006) examined the impact of point-source subtraction residuals on 21 cm fluctuations, and Morales & Hewitt (2004) developed a removal approach based on some symmetries in the 3D Fourier representation of radio interferometric data. More recently, simulations and observations have shown that even if sky foregrounds are spectrally smooth, the chromatic synthesized beam of an interferometer interacting with smooth spectrum foregrounds can still imprint unsmooth spectral features on measured foregrounds, contaminating a “wedge”-like region (the so-called “mode-mixing”) in cylindrical $(k_⊥, k_‖)$ Fourier space (Datta et al. 2010; Morales et al. 2012; Parsons et al. 2012; Trott et al. 2012; Vedantham et al. 2012; Dillon et al. 2013; Pober et al. 2013a, 2013b, 2014; Liu et al. 2014a, 2014b).

In order to subtract foregrounds with minimal prior assumptions on the angular distribution or frequency dependence of foreground components, some *blind* (non-parametric) methods have been extensively studied in the literature. By using only data, those methods can automatically find the correlated features and separate them out into distinct components. A commonly used data analysis technique known as Principal Component Analysis (PCA) has been applied to foreground removal and has been successfully tested on simulated data and real observations (de Oliveira-Costa et al. 2008; Chang et al. 2010; Masui et al. 2013; Switzer et al. 2013). Other non-parametric methods, including the FASTICA algorithm (Hyvärinen 1999; Hyvärinen & Oja 2000; Chapman et al. 2012; Wolz et al. 2014) based on independent component analysis (ICA) and singular value decomposition together with an analysis based on the Karhunen–Loève transform (Shaw et al. 2014), provide encouraging results. The FASTICA algorithm is based on the maximization of non-Gaussianity (neg-entropy), a measure of the deviation of a mixture of signals from a Gaussian distribution. However, the disadvantage of this approach is that the noise leakage into the reconstructed signal is, in principle, difficult to estimate since the signal and noise essentially have the same statistical distribution, which is completely uncorrelated across frequencies. Other methods, such as the Correlated Component Analysis (CCA; Bonaldi & Brown 2015) method, can be referred to as “semi-blind” since they require some prior knowledge. By using second-order statistics, CCA also provides a promising way to separate the 21 cm reionization signal from foreground contamination for SKA measurements. Recently, Vansyngel et al. (2014) provided a unified approach for semi-blind foreground cleaning from multi-frequency CMB data.

In this study, we describe HIEMICA, the H1 Expectation-Maximization Independent Component Analysis algorithm, a Bayesian framework to infer power spectra and maps of the H1 signal and foregrounds from “dirty” data cubes. This technique extends the Spectral Matching ICA (SMICA; Delabrouille et al. 2003; Cardoso et al. 2008; Planck Collaboration et al. 2014) approach from two-dimensional (2D) CMB maps to the three-dimensional (3D) 21 cm cosmological signal and can be considered a semi-blind source separation method. The statistical isotropy of the cosmological Gaussian random fields in 3D Fourier space provides a characteristic and unique signature to distinguish them from the statistically non-

<footnote n="11">¹¹ http://www.lofar.org</footnote>

<footnote n="12">¹² http://www.mwatelescope.org</footnote>

<footnote n="13">¹³ http://21 cma.bao.ac.cn</footnote>

<footnote n="14">¹⁴ http://astro.berkeley.edu/~dbacker/eor/</footnote>

<footnote n="15">¹⁵ http://gmrt.ncra.tifr.res.in</footnote>

<footnote n="16">¹⁶ http://chime.phas.ubc.ca</footnote>

<footnote n="17">¹⁷ http://www.skatelescope.org</footnote>
isotropic foregrounds; in the frequency (redshift, or radial) direction, astrophysical foreground sources contain only spectral information, which has no connection with their spatial distributions in the transverse direction.

The purpose of this work is to take a first step toward building a new framework for 21 cm data analysis. Our paper provides a complete description of the mathematical formalism and numerical techniques, which in principle can be applied to any realistic array design. Here, we test the approach for an idealized array; we propose dedicating a future paper to systematic tests of arrays with realistic instrumental effects.

This paper is organized as follows. Section 2 describes the HIEMICA semi-blind component separation algorithm. Section 3 summarizes the details of our performed simulations. Section 4 presents the application of HIEMICA to simulated sky maps and shows the main results. Finally, Section 5 provides our discussion and concluding remarks.

## 2. THE HIEMICA METHOD

### 2.1. Radio Interferometric Measurements

For a multichannel interferometric array, the fundamental observable is the visibility, $V(u, v)$, which is the angular Fourier mode of the intensity fluctuations on the sky measured by a baseline $u$ in a frequency channel $v$. In the flat-sky approximation, the complex visibility can be written simply as the 2D Fourier transform of the beam-modulated sky intensity:

$$V(u, v) = \int d^2\theta A(\theta, v) T(\theta, v) e^{-2\pi i u \cdot \theta}, \quad (1)$$

where $\theta$ denotes the angular sky position, $T(\theta, v)$ is the sky temperature, $A(\theta, v)$ is the primary beam pattern, which determines the observed area of the sky, and the baseline vector $u$ is on the uv-plane perpendicular to the direction of the incoming signals and measures the separation between the two antenna pairs in units of wavelength, i.e., $u = v/c\Delta r$. As the physical separation $\Delta r$ of each pair of antennas is assumed to be fixed, each frequency channel thus gives an independent set of visibility samplings in the uv-plane with frequency-dependent baseline vectors.

In practice, especially for 21 cm tomography, one discretizes the 3D sky into voxels so as to efficiently estimate the intensity distribution $T(\theta, v)$ and associated power spectrum using FFT-based algorithms. If we assume $n_{\perp}$ pixels per 2D image (real or Fourier space) with even spacings in the $\theta_x$ and $\theta_y$ directions (i.e., two transverse directions), and an instrument with $n_{\parallel}$ frequency channels of an equal spectral resolution $\Delta\nu$ to probe the sky signal in a radial direction, then we have a 3D uniform Cartesian grid of size $n_{\perp}n_{\parallel}$. Throughout this paper, the symbols $\perp$ and $\parallel$ denote the two transverse directions and the radial direction, respectively.

Following Sutter et al. (2012), the radio interferometric measurement equation of Equation (1) at the frequency $\nu_i$, for $i = 1, 2, \dots, n_{\parallel}$, can be recast into a discrete operator formalism by

$$y_{\perp}[\nu_i] = I[\nu_i]F_{\perp}A[\nu_i]F_{\perp}^{-1}\tilde{\psi}_{\perp}[\nu_i] + n_{\perp}[\nu_i], \quad (2)$$

where $y_{\perp}[\nu_i] = (y[u_1, \nu_i], \dots, y[u_{n_{\perp}}, \nu_i])^T$ and $\tilde{\psi}_{\perp}[\nu_i] = (\tilde{\psi}[u_1, \nu_i], \dots, \tilde{\psi}[u_{n_{\perp}}, \nu_i])^T$ denote, respectively, the discretized visibilities in the uv-plane and the 2D Fourier “image” of the true angular sky ($\psi_{\perp}[\nu_i]$) of 21 cm and foregrounds at the $i$th

channel. The vector $n_{\perp}[\nu_i] = (n[u_1, \nu_i], \dots, n[u_{n_{\perp}}, \nu_i])^T$ denotes the instrumental noise at the $i$th channel for all of the uv points. Hereafter, bold, upper and lower case letters denote a matrix and a vector, respectively, and the symbol $\sim$ is used for Fourier components. $F_{\perp}$ represents the 2D Fourier transform operator performed in the transverse direction, converting from the real-space domain into the uv-space domain (of course, $F_{\perp}^{-1}$ is its inverse), and $A[\nu_i]$ represents the primary beam pattern of the interferometer measured at the frequency $\nu_i$ and implicitly is a $n_{\perp} \times n_{\perp}$ diagonal matrix. The operator $I[\nu_i]$ specifies the visibility sampling function determined by the distribution of baselines and can also be represented by an $n_{\perp} \times n_{\perp}$ diagonal matrix that has 1s and 0s on the diagonal to indicate whether a given pixel in uv-space has been observed or not.

It is convenient to introduce a vector form for numerical operations. By stacking each vector $\tilde{\psi}_{\perp}[\nu_i]$ as $\tilde{\psi}_{\perp} = [\tilde{\psi}_{\perp}^T[\nu_1], \dots, \tilde{\psi}_{\perp}^T[\nu_{n_{\perp}}]]^T$, Equation (2) can account for all of the frequency channels when rewritten as

$$y = \Phi \tilde{\psi}_{\perp} + n, \quad \Phi = \begin{pmatrix} \Phi_{\perp} & 0 \\ 0 & \Phi_{n_{\parallel}} \end{pmatrix}, \quad (3) \\ \text{with } \Phi_i = I[\nu_i]F_{\perp}A[\nu_i]F_{\perp}^{-1},$$

where $\tilde{\psi}_{\perp}$ is a vector with length $n_{\perp}n_{\parallel}$ and $\Phi$ is a block diagonal matrix, while $y = [y_{\perp}^T[\nu_1], \dots, y_{\perp}^T[\nu_{n_{\perp}}]]^T$ and $n = [n_{\perp}^T[\nu_1], \dots, n_{\perp}^T[\nu_{n_{\perp}}]]^T$ each collect vectors for all frequencies into an $n_{\perp}n_{\parallel}$-element vector.

Note that Equation (3) is presented on a uniform grid in order to enable the FFT operation. Realistic observations from radio interferometers do not sample visibility data on a uniform grid in uv-space. Therefore, analyzing the data on a uniform grid requires an additional process, “gridding,” to map raw visibilities from specified uv points to grids with even spacing. However, the gridding process can distort the signal estimate, especially for small 21 cm signals. Alternatively, one can use the Non-uniform Discrete Fourier Transform (NuDFT; Fessler & Sutton 2003) to directly calculate visibility data at continuous coordinates so as to avoid “gridding” effects. We leave an investigation of the gridding effects and a possible solution by NuDFT to future work.

### 2.2. Linear Mixture Model

The 21 cm signal along the line of sight for each individual pixel is not a smooth function of frequency, while most astrophysical contaminants should have smooth power-law spectral structures. If we assume that the total foreground emission contributing to the observed data in a given frequency channel can be expressed as a linear superposition of a few components and the spectral behavior of each component does not vary across the sky, the true sky temperature in pixel $\theta_k$, for $k = 1, 2, \dots, n_{\perp}$, in the $i$th channel can then be modeled as

$$\psi[\nu_i, \theta_k] = s[\nu_i, \theta_k] + \sum_{j=1}^{n_c} M_{ij}f_j[\theta_k], \quad (4)$$

where $\nu_i$ is the frequency, $\theta_k$ is the pixel position, $n_c$ is the assumed total number of foreground components, $s[\nu_i, \theta_k]$ is the 21 cm signal, and $f_j[\theta_k]$ is the $j$th foreground component. $M_{ij}$ reflects the contribution of the $j$th foreground component to
the ith frequency channel and is assumed to be constant and independent of θ<sub>k</sub>, and thus the matrix **M** is commonly called the “mixing matrix” with dimensions of n<sub>‖</sub> × n<sub>c</sub>. Each column of **M** encodes the frequency dependence of each foreground component.

As the Fourier modes of the isotropic 21 cm signal s(ν, θ) are expected to be mutually independent, expressing its temperature in the Fourier coordinates (η, **u**) is particularly useful. By taking the 3D Fourier transform along the frequency axis and two angular directions, one has

$$
\tilde{s}(\eta, \mathbf{u}) = \int s(\nu, \theta) e^{-2\pi i (\nu\eta + \mathbf{u}\cdot\theta)} d\theta d\nu, \quad (5)
$$

where η is the Fourier dual of the frequency variable ν. Note
that here we perform the “true” Fourier transform of the sky
temperature along the line-of-sight direction that is perpendi-
cular to the uv-plane, rather than a “delay transform.” The
delay transform is the Fourier transform in the frequency of the
spectrum measured by a single baseline, which varies with
frequency (i.e., **u** is a function of ν in Equation (5)) and which
leads to the “mode-mixing” phenomenon (see Liu et al. 2014a
for details).

The foreground components can also be expressed in terms
of their 2D Fourier modes. For the jth component, one has

$$
\tilde{f}_j(\mathbf{u}) = \int f_j(\theta) e^{-2\pi i \mathbf{u}\cdot\theta} d\theta \quad (j = 1, \dots, n_c). \quad (6)
$$

According to the above equations, by performing the angular Fourier transform on Equation (4) and expressing the 21 cm signal in terms of its 3D Fourier modes, Equation (4) can be rewritten in a discretized version as

$$
\tilde{\psi}[\nu_i, u_k] = \sum_{m=1}^{n_{\parallel}} F_{\parallel}^{-1}[i, m] \tilde{s}[\eta_m, u_k] + \sum_{j=1}^{n_c} M_{ij} \tilde{f}_j[u_k], \quad (7)
$$

where *i* denotes the frequency channel number, *k* denotes the
pixel position in the uv-plane, and we have introduced the
operator (matrix) *F*<sub>‖</sub><sup>-1</sup>, with its elements given by Fourier
coefficients *F*<sub>‖</sub><sup>-1</sup>[*i*, *m*] (*i* = 1, ..., *n*<sub>‖</sub>; *m* = 1, ..., *n*<sub>‖</sub>], to repre-
sent the one-dimensional inverse Fourier transform converting
the *η* domain to the frequency domain.

Using matrix multiplications, we can rewrite the linear
system above as the following matrix equation:

$$
\tilde{\Psi} = F_{\parallel}^{-1} \tilde{S} + M \tilde{F}, \quad (8)
$$

where $\tilde{\psi}[\nu_i, u_k]$ is the $(i, k)$th entry of the matrix $\tilde{\Psi}$, $\tilde{s}[\eta_m, u_k]$ is the $(m, k)$th entry of the matrix $\tilde{S}$, and $\tilde{f}_j[u_k]$ is the $(j, k)$th entry of the matrix $\tilde{F}$. Applying a “vectorization” operator $(\text{vec}\,\cdot)$ that converts the matrix into a column vector by stacking the columns into a long column vector, one can express $\tilde{\psi}_{\perp}$ in Equation (3) with $\tilde{\Psi}^T$ by

$$
\tilde{\psi}_{\perp} = \mathrm{vec}[\tilde{\Psi}^T] = \mathrm{vec}[(F_{\parallel}^{-1}\tilde{S})^T] + \mathrm{vec}[(M\tilde{F})^T]. \quad (9)
$$

Substituting Equation (9) into Equation (3), and using the identity $\mathrm{vec}[\mathbf{PQ}] = (\mathbf{Q}^T \otimes \mathbf{I}_m)\mathrm{vec}[\mathbf{P}]$ for any matrices $\mathbf{P}_{m \times n}$ and $\mathbf{Q}_{n \times p}$, where $\otimes$ denotes the Kronecker product and $\mathbf{I}_m$ is an $(m \times m)$ identity matrix, yields

$$
\tilde{\psi}_{\perp} = (\mathbf{B} \otimes \mathbf{I}_{n_{\perp}}) \mathbf{x}, \quad \mathbf{B} = (\mathbf{F}_{\parallel}^{-1} \mathbf{M}), \quad (10)
$$

where **x** = (tilde{s}^T, tilde_f^T)^T with tilde_s = vec[tilde{s}^T], tilde_f = vec[tilde{f}^T], I_n⊥ is the n⊥ × n⊥ identity matrix, and **B** is a partitioned matrix with dimensions of n_‖ × (n_‖ + n_c). Explicitly, tilde_s and tilde_f are vectors of lengths n_⊥n_‖ and n_⊥n_c, respectively. They are constructed by stacking tilde_s[η_m] and tilde_f_j as tilde_s = (tilde{s}^T[η_1], ..., tilde{s}^T[η_m))^T and tilde_f = (tilde{f}_1^T, ..., tilde{f}_{n_c}^T)^T, where the vectors tilde_s[η_m] and tilde_f_j, each with length n_⊥, denote an image of the 3D Fourier modes of the 21 cm signal at η_m and an image of the 2D Fourier modes of the jth foreground component, respectively, i.e., tilde_s[η_m] = (tilde{s}[u_1, η_m], ..., tilde{s}[u_n⊥, η_m]) and tilde_f_j = (tilde{f}_j[u_1], ..., tilde{f}_j[u_n⊥]).

By inserting Equation (10) into Equation (3), the measure-
ment equation finally can be expressed as

$$
\mathbf{y} = \mathbf{H}\mathbf{x} + \mathbf{n}, \quad \mathbf{H} = \mathbf{\Phi}(\mathbf{B} \otimes \mathbf{I}_{n_{\perp}}). \quad (11)
$$

2.3. ICA Assumption

The ICA assumption is that the data can be considered as a
linear mixture of a set of statistically mutually independent
components (ICs). The cosmological 21 cm signal is expected
to be well approximated by an isotropic Gaussian random field
and uncorrelated with foregrounds. If we also assume that the
diffuse foregrounds consist of several statistically ICs, each of
them also being an isotropic Gaussian random field with zero
mean, as defined in Equation (10), the covariance matrix of **x**
then becomes a diagonal matrix, namely,

$$
\mathbf{C} = \langle \mathbf{x} \mathbf{x}^{\dagger} \rangle = \begin{pmatrix} \langle \tilde{s} \tilde{s}^{\dagger} \rangle & 0 \\ 0 & \langle \tilde{f} \tilde{f}^{\dagger} \rangle \end{pmatrix}. \quad (12)
$$

The unknown diagonal matrices 〈 ˜s ˜s $^{†}$〉 and 〈 ˜f ˜f $^{†}$〉 are to be
estimated from the observed data. They determine the 3D
power spectrum of the H I signal and the angular power spectra
of the ICs, defined through

$$
\langle \tilde{s}[u_k, \eta_m] \tilde{s}^*[u_{k'}, \eta_{m'}] \rangle = P_{H1}(u_k, \eta_m) \delta_{kk'} \delta_{mm'} \quad (13)
$$

$$
\langle \tilde{f}_j(u_k) \tilde{f}_j^*(u_k') \rangle = \delta_{jj'} \delta_{kk'} C_j^j (\ell = 2\pi|\mathbf{u}|), \quad (14)
$$

where $P_{\text{H1}}(u_k, \eta_m)$ denotes the 3D power spectrum of the 21 cm signal as a function of **u** and η, and $C_j^j(\ell)$ denotes the angular power spectrum of the jth foreground component as a function of the multipole *ℓ*. The relation *ℓ* = 2π|**u**|, where **u** is determined by the pixel position *u*<sub>*k*</sub>, has been used in the flat-sky approximation. Note that physical sources such as the Galactic synchrotron and free-free emissions actually do have a non-zero cross-correlation and their spatial distributions clearly appear to be non-isotropic. Nevertheless, as we will demonstrate, the ICA assumption can be regarded as an effective decomposition of sources and does not appear to affect our ability to remove foregrounds.

2.3.1. Noise

We also assume that the instrument noise is an uncorrelated Gaussian distribution and can be obtained from a reasonable noise model through

$$
\langle \tilde{n}(u_k, \nu_i) \tilde{n}^*(u_{k'}, \nu_{i'}) \rangle = P_N(u_k, \nu_i) \delta_{kk'} \delta_{ii'}. \quad (15)
$$

The covariance matrix thus has the form of a known block
diagonal matrix (with **n** defined in Equation (3)), where each
block is also diagonal:

$$N = \langle nn^{\dagger} \rangle = \begin{pmatrix} N_1 & & & 0 \\ 0 & \ddots & & N_{n|} \\ & & \ddots & \\ & & & N_{n|} \end{pmatrix}. \quad (16)$$

Here, $N_i$ denotes the $n_⊥ \times n_⊥$ covariance matrix of the noise map in the uv-plane at the ith frequency channel, which is also diagonal.

## 2.4. The Expectation-Maximization (EM) Algorithm

The goal of our analysis is to identify and separate the components from visibility data that contains mixtures of foregrounds and signal. Using the Bayesian framework proposed by Snoussi et al. (2002) for CMB data, we will show that, without any assumption on priors, by performing a semi-blind ICA for which only the noise covariance matrix $N$ is well known, one can successfully separate the components by jointly estimating the covariance matrix $C$ and the mixing matrix $M$, and can accurately estimate the power spectrum of the $H_1$ signal.

In Bayesian inference, information about unknown parameters $\theta$ that we want to estimate from the data $y$ is expressed in the form of a posterior probability distribution. Using Bayes' theorem, it can be computed as $p(\theta|y) \propto L(y|\theta)p(\theta)$, where $L(y|\theta)$ is the likelihood and $p(\theta)$ is the prior distribution of $\theta$; $\theta = (C, M)$ in our case. If we assume flat priors for the mixing matrix and the power spectra of all of the sources, and assume that they are uncorrelated and independent, then the prior reduces to $P(\theta) = P(C)P(M) \propto 1$. Thus, exploring the observed data posterior $p(y|\theta)$ is equivalent to exploring the likelihood. Following Snoussi et al. (2002), given the data model in Equation (11), $y = Hx + n$, the mixing matrix $M$ and the covariance matrix $C$ can be estimated by maximizing the likelihood function. For independent and Gaussian sources, the likelihood can be expressed as

$$-2 \ln L(y|\theta) \propto \ln|HCH^\dagger + N| + \text{Tr}[(HCH^\dagger + N)^{-1}yy^\dagger]. \quad (17)$$

Unfortunately, in general, such a likelihood evaluation is computationally intractable when applied to large data sets with realistic interferometric observations. However, in the analogous case of CMB imaging observations<fnref n="18" />, this likelihood can be approximated by the SMICA-likelihood (Delabrouille et al. 2003) to measure a “spectral mismatch” in the Fourier domain between the empirical covariance matrices of the data and their ensemble averages, which depend on the estimated parameters. In practice, minimization of the spectral mismatch is equivalent to maximizing the likelihood and can be achieved with the EM algorithm. The EM algorithm has been applied to estimate the CMB power spectrum and reconstruct the CMB map from multi-frequency microwave maps (Delabrouille et al. 2003). The EM algorithm is an elaborate technique to find the maximum-likelihood estimate of the parameters when it is intractable to directly evaluating the likelihood function. The calculation is simplified by assuming the existence of additional, missing parameters. It is an iterative algorithm to

repeatedly solve a tractable complete-data problem instead of solving a difficult incomplete-data problem.

Here, we briefly summarize the EM algorithm when applied to our model. The likelihood $p(y|\theta)$ first can be obtained by marginalizing the joint distribution $p(y, x|\theta)$ over the missing data $x$ as

$$p(y|\theta) = \int p(y, x|\theta)dx. \quad (18)$$

The key idea is that the EM algorithm does not maximize $p(y|\theta)$ directly; instead, it maximizes the so-called EM-functional as follows:

$$\begin{aligned} Q(\theta|\theta^n) &= E[\ln p(y, x|\theta)|y, \theta^n] \\ &= \int \ln p(y, x|\theta)p(x|y, \theta^n)dx, \end{aligned} \quad (19)$$

which is the expected value of the complete-data log-likelihood with respect to the missing data $x$ given the observed data $y$ and the current parameter estimates $\theta^n$. The evaluation of this expectation is called the Expectation (E)-step of the algorithm.

In our model, the prior distribution of the complex Fourier modes $x$ is assumed to be Gaussian with zero mean,

$$p(x|C) \propto \exp\left(-\frac{1}{2}x^{\dagger}C^{-1}x\right), \quad (20)$$

where the diagonal covariance matrix $C$ is defined in Equation (12). Within a Bayesian framework, the joint posterior distribution $p(y, x|\theta)$ can be expressed by $p(y, x|\theta) = p(y|x, \theta)p(x|\theta)$, yielding

$$-2 \ln p(y, x|\theta) = \ln|N| + (y - Hx)^{\dagger}N^{-1}(y - Hx) + \ln|C| + x^{\dagger}C^{-1}x + \text{cst.} \quad (21)$$

Using Bayes' rule, the conditional probability distribution function $p(x|y, \theta^n)$ for the signal given the data is also the Gaussian,

$$p(x|y, \theta^n) \propto \exp\left(-\frac{1}{2}(x - \hat{x})^\dagger\Sigma^{-1}(x - \hat{x})\right), \quad (22)$$

where

$$\hat{x} = [H^{\dagger}N^{-1}H + C^{-1}]^{-1}H^{\dagger}N^{-1}y, \quad (23)$$

$$\Sigma = (H^{\dagger}N^{-1}H + C^{-1})^{-1}, \quad (24)$$

where $\hat{x}$ is the so-called Wiener-Filtered (WF) map and $\Sigma$ is the corresponding covariance $\Sigma = \langle(x - \hat{x})(x - \hat{x})^{\dagger}\rangle$. The solution for $\hat{x}$ is the general map-making problem in cosmology. This WF map can be computed efficiently by implementing a preconditioned conjugate-gradient method that allows one to iteratively reach the solution in a tractable amount of computation time.

Using Equations (21) and (22), we integrate out $x$ to derive the expression of Equation (19), given by

$$\begin{aligned} Q(\theta|\theta^n) = \text{cst.} & -\frac{1}{2}\ln|N| - \frac{1}{2}\ln|C| - \frac{1}{2}\text{Tr}[C^{-1}\hat{R}_{xx}^\dagger] \\ & -\frac{1}{2}\text{Tr}[N^{-1}(\hat{R}_{yy} - H\hat{R}_{yx}^\dagger - \hat{R}_{yx}H^\dagger + H\hat{R}_{xx}H^\dagger)], \end{aligned} \quad (25)$$

<footnote n="18">¹⁸ The convolution with a frequency-dependent beam response in real space becomes a simple product in Fourier space.</footnote>
where

$$ \hat{\mathbf{R}}_{yy} = yy^{\dagger}, \quad \hat{\mathbf{R}}_{yx} = y\hat{x}^{\dagger}, \quad \hat{\mathbf{R}}_{xx} = \Sigma + \hat{x}\hat{x}^{\dagger}. \qquad (26) $$

Note here that $\hat{\mathbf{R}}_{yy}$, $\hat{\mathbf{R}}_{yx}$, and $\hat{\mathbf{R}}_{xx}$ only depend on $\theta^n$ rather than $\theta$.

The second step, called the *Maximization (M)-step*, updates the parameters by maximizing the expectation we computed in the previous *E-step*. These two steps are repeated as necessary. This procedure is guaranteed to increase the likelihood $p(y|\theta)$ monotonically with successive iterations. In order to obtain the parameter $\theta^{n+1}$ at iteration $n+1$, we solve the gradient equation with respect to $M$ and $C$ to maximize the functional $Q$. To do so, let us first introduce some notations that allow us to express the derivative simply.

Let $\hat{x}$ and $\Sigma$ be partitioned as follows:

$$ \hat{x} = (\hat{s}^T, \hat{f}^T)^T, \qquad (27) $$

with

$$ \hat{s} = (\hat{s}_{\perp}^T [\eta_1], \dots, \hat{s}_{\perp}^T [\eta_{n_{\parallel}}])^T, \quad \hat{f} = (\hat{f}_1^T, \dots, \hat{f}_{n_c}^T)^T. \qquad (28) $$

Here, $\hat{x}$ is equally divided into $(n_{\parallel} + n_c)$ subvectors, each with length $n_{\perp}$, $\hat{s}_{\perp}$ [$\eta_i$] denotes the WF map of the 21 cm signal at $\eta_i$, and $\hat{f}_j$ denotes the WF map of the $j$th foreground component. Similarly,

$$ \Sigma = \begin{pmatrix} \Sigma^{\bar{s}\bar{s}} & \Sigma^{\bar{s}\bar{f}} \\ \Sigma^{\bar{f}\bar{s}} & \Sigma^{\bar{f}\bar{f}} \end{pmatrix} \qquad (29) $$

with

$$ \Sigma^{\bar{s}\bar{s}} = \begin{pmatrix} \Sigma_{11}^{\bar{s}\bar{s}} & \cdots & \Sigma_{1n_{\perp}}^{\bar{s}\bar{s}} \\ \vdots & \ddots & \vdots \\ \Sigma_{n_{\perp}1}^{\bar{s}\bar{s}} & \cdots & \Sigma_{n_{\perp}n_{\perp}}^{\bar{s}\bar{s}} \end{pmatrix} \qquad \Sigma^{\bar{s}\bar{f}} = \begin{pmatrix} \Sigma_{11}^{\bar{s}\bar{f}} & \cdots & \Sigma_{1n_c}^{\bar{s}\bar{f}} \\ \vdots & \ddots & \vdots \\ \Sigma_{n_{\perp}1}^{\bar{s}\bar{f}} & \cdots & \Sigma_{n_{\perp}n_c}^{\bar{s}\bar{f}} \end{pmatrix} \\ \Sigma^{\bar{f}\bar{f}} = \begin{pmatrix} \Sigma_{11}^{\bar{f}\bar{f}} & \cdots & \Sigma_{1n_c}^{\bar{f}\bar{f}} \\ \vdots & \ddots & \vdots \\ \Sigma_{n_c1}^{\bar{f}\bar{f}} & \cdots & \Sigma_{n_c n_c}^{\bar{f}\bar{f}} \end{pmatrix}, \qquad (30) $$

where

$$ \Sigma_{ij}^{\bar{s}\bar{s}} = \langle (\bar{s}[\eta_i] - \hat{s}[\eta_i])(\bar{s}[\eta_j] - \hat{s}[\eta_j])^{\dagger} \rangle, \qquad (31) $$

$$ \Sigma_{ij}^{\bar{s}\bar{f}} = \langle (\bar{s}[\eta_i] - \hat{s}[\eta_i])(\bar{f}[\eta_j] - \hat{f}[\eta_j])^{\dagger} \rangle, \qquad (32) $$

$$ \Sigma_{ij}^{\bar{f}\bar{f}} = \langle (\bar{f}[\eta_i] - \hat{f}[\eta_i])(\bar{f}[\eta_j] - \hat{f}[\eta_j])^{\dagger} \rangle, \qquad (33) $$

each with dimensions of $n_{\perp} \times n_{\perp}$.

The gradient equation then reads

$$ = \frac{\partial Q}{\partial \mathbf{C}} = -(\mathbf{C}^{-1}) + \mathbf{C}^{-1}\hat{\mathbf{R}}_{xx}\mathbf{C}^{-1} = 0 \qquad (34) $$

$$ \frac{\partial Q}{\partial M_{ij}} = \mathrm{Tr}[N_i^{-1}(\mathbf{\Phi}_i \hat{f}_j y_i^\dagger + (\mathbf{\Phi}_i \hat{f}_j y_i^\dagger)^\dagger)] \\ - \frac{\partial \mathrm{Tr}[N^{-1} \mathbf{H} \hat{\mathbf{R}}_{xx} \mathbf{H}^\dagger]}{\partial M_{ij}} = 0, \qquad (35) $$

where the derivative of the last term can be well approximated by

$$ \frac{\partial \mathrm{Tr}[N^{-1} \mathbf{H} \hat{\mathbf{R}}_{xx} \mathbf{H}^{\dagger}]}{\partial M_{ij}} \approx \sum_{k=1}^{n_c} M_{ik} \mathrm{Tr}[N_i^{-1} \mathbf{\Phi}_i (\Sigma_{kj}^{\bar{f}\bar{f}} + \hat{f}_k \hat{f}_j^{\dagger} + (\Sigma_{kj}^{\bar{f}\bar{f}} + \hat{f}_k \hat{f}_j^{\dagger})^{\dagger}) \mathbf{\Phi}_i^{\dagger}], \qquad (36) $$

where $\mathbf{\Phi}_i = \mathbf{I}[\nu_i] \mathbf{F}_{\perp} \mathbf{A}[\nu_i] \mathbf{F}_{\perp}^{-1}$, as defined in Equation (2), and we have neglected the terms that include $\Sigma^{\bar{s}\bar{f}}$ and $\Sigma^{\bar{s}\bar{s}}$, although the complete expression can be obtained in a straightforward manner. This is an excellent approximation for 21 cm foregrounds since they are several orders of magnitude brighter than the 21 cm signal.<fnref n="19" /> By making this approximation, we are able to dramatically speed up the computation. According to Equations (34), (35), and (37), one can establish the *update scheme* as follows:

$$ \mathbf{M}_{i,:} \leftarrow (\rho_i^{\bar{f}y}[1], \dots, \rho_i^{\bar{f}y}[n_c]) \times \begin{pmatrix} \rho_i^{\bar{f}\bar{f}}[1, 1] \cdots \rho_i^{\bar{f}\bar{f}}[1, n_c] \\ \vdots & \ddots & \vdots \\ \rho_i^{\bar{f}\bar{f}}[n_c, 1] \cdots \rho_i^{\bar{f}\bar{f}}[n_c, n_c] \end{pmatrix}^{-1}, \qquad (37) \\ C_{ii} \leftarrow (\hat{\mathbf{R}}_{xx})_{ii}, \qquad (38) $$

where $\mathbf{M}_{i,:}$ denotes the $i$th row of the mixing matrix, and

$$ \rho_i^{\bar{f}\bar{f}}[k, j] = \mathrm{Re}\{\mathrm{Tr}[N_i^{-1} \mathbf{\Phi}_i (\Sigma_{kj}^{\bar{f}\bar{f}} + \hat{f}_k \hat{f}_j^{\dagger}) \mathbf{\Phi}_i^{\dagger}]\} \qquad (39) $$

$$ \rho_i^{\bar{f}y}[j] = \mathrm{Re}\{\mathrm{Tr}[N_i^{-1} \mathbf{\Phi}_i \hat{f}_j y_i^{\dagger}]\}. \qquad (40) $$

One can verify that by setting $\mathbf{\Phi}_i$ as an identity matrix, for $i = 1, 2, \dots, n_{\parallel}$, such that no instrumental effects are present, this is in complete agreement with the structure of the expression for the mixing matrix in the literature (Snoussi et al. 2002; Delabrouille et al. 2003).

The problem here is evaluating the matrix $\Sigma$ from its inverse as given by Equation (24) in order to obtain the covariance matrix $C$ from Equation (38). Since the matrix $\Sigma$ has dimensions of $n_{\perp}$ ($n_{\parallel} + n_c$) $\times n_{\perp}$ ($n_{\parallel} + n_c$), which are typically of the order of $10^6 \times 10^6$ in a 21 cm survey, such matrix inversion is completely computationally prohibitive. However, since the power spectrum is only determined by the diagonal components of $\Sigma$, not by the off-diagonals, by using the same trick presented in the Gibbs sampling techniques (Jewell et al. 2004; Wandelt et al. 2004) for CMB data analysis, one can verify that the ensemble-averaged solution for the following linear equation has the desired covariance. By solving

$$ [\mathbf{H}^{\dagger} \mathbf{N}^{-1} \mathbf{H} + \mathbf{C}^{-1}] \mathbf{z}_1 = \mathbf{H}^{\dagger} \mathbf{N}^{-1/2} \mathbf{z}_1 + \mathbf{C}^{-1/2} \mathbf{z}_2, \qquad (41) $$

where the real vectors $\mathbf{z}_1$ and $\mathbf{z}_2$ are of lengths $n_{\perp}n_{\parallel}$ and $n_{\perp}(n_c + n_{\parallel})$, respectively, with elements drawn from a standard normal distribution, one can find that $\Sigma = \langle \mathbf{z}\mathbf{z}^{\dagger} \rangle$, where $\langle \cdot \rangle$ represents the average over all of the solutions for a large number of realizations of $\mathbf{z}_1$ and $\mathbf{z}_2$. This algorithm can rapidly

<footnote n="19">¹⁹ Note that in the case of CMB, the evaluation of the mixing matrix not only depends on CMB foregrounds but also on the CMB signal itself since they are roughly comparable.</footnote>
provide accurate approximate solutions in massively parallel computers. Therefore one can obtain the diagonal components of $\hat{\mathbf{R}}_{xx}$ in Equation (38) efficiently by $(\hat{\mathbf{R}}_{xx})_{ii} = \langle |z_i|^2 \rangle + |\hat{x}_i|^2$, where $i$ runs over all of the elements. In the same way as in $\hat{x}$, the vector $z$ with length $(n_\parallel + n_c)n_\perp$ can also be equally split into $(n_\parallel + n_c)$ subvectors, each with length $n_\perp$, as $z = ((z_i^\parallel)^T, \dots, (z_{n_\parallel}^\parallel)^T, (z_i^\perp)^T, \dots, (z_{n_c}^\perp)^T)^T$, corresponding to contributions from the 21 cm signal, denoted by $z_i^\parallel$ for $i = 1, \dots, n_\perp$, and the ICs, denoted by $z_i^\perp$ for $i = 1, \dots, n_c$. Note that, in general, the noise could be correlated, resulting in non-zero off-diagonal components in the noise covariance matrix, $N$. Solving Equation (41) will require multiplication by a dense $n_\perp n_\parallel \times n_\perp n_\parallel$ inverse noise covariance matrix, with a computational scaling of $O(n_\perp^2 n_\parallel^2)$. In the most general case, inverting $N$ is a significant computational challenge. However, the inverse noise covariance matrix needs only to be computed once and stored, and then can be used many times to compute Equations (39)–(41).

Using the above-mentioned trick, the trace term in Equation (39) can be computed simply by the ensemble-averaged dot products of two vectors as

$$ \mathrm{Tr}[N_i^{-1} \mathbf{\Phi}_i (\mathbf{\Sigma}_{kj}^{\hat{\mathbf{f}}_j \hat{\mathbf{f}}_k} + \hat{\mathbf{f}}_j \hat{\mathbf{f}}_k^\dagger) \mathbf{\Phi}_i^\dagger] \\ = \langle (\mathbf{\Phi}_i z_j^\parallel)^\dagger (N_i^{-1} \mathbf{\Phi}_i z_k^\parallel) \rangle + (\mathbf{\Phi}_i \hat{\mathbf{f}}_j)^\dagger N_i^{-1} (\mathbf{\Phi}_i \hat{\mathbf{f}}_k), \quad (42) $$

where the vectors $z_j^\parallel$ and $z_k^\parallel$ represent the subvectors of $z$ solved by Equation (41). The above equation also applies to the case of correlated noise.

Moreover, in each EM iteration, we fix the magnitude of each column of $M$ to unity and adjust the corresponding power spectra of the ICs accordingly, similar to Delabrouille et al. (2003), to break the degeneracy and keep the scale-invariant product $\sum_j M_{ij} M_{i'j} \langle \tilde{f}_j [u_k] \tilde{f}_k^* [u_k] \rangle = \sum_j M_{ij} M_{i'j} C_j^i [u_k]$ unchanged for an arbitrary pixel $u_k$ at any frequency channels $i$ and $i'$, as seen from Equation (7).

After updating $M$ and $C$, since the quantities $\ell^2 C_i^\ell$ for the ICs and $k_\parallel k_\perp^2 P_{H1}$ for the $H1$ signal are expected to vary more slowly than $C_i^\ell$ and $P_{H1}$ themselves, it is more appropriate to perform bin averages after each iteration and update the associated elements in $C$ while ensuring that the number of observations is smaller than the number of estimated parameters in $C$ and $M$. The simulation shows that using band-averaged spectra can significantly increase the convergence speed and obtain much more stable spectra with high accuracy in iterations. According to Equation (13), the bin averages then yield the following estimates:

$$ C_j^i(\mathbf{u}) \leftarrow \left( \sum_{\mathbf{u} \in D_q} |\mathbf{u}|^2 \right)^{-1} \sum_{\mathbf{u} \in D_q} |\mathbf{u}|^2 C_j^i(\mathbf{u}) \quad (43) $$

$$ P_{H1}(\mathbf{u}, \eta) \leftarrow \left( \sum_{\substack{|\mathbf{u}| \in D_q \\ \eta \in D_z}} \eta |\mathbf{u}|^2 \right)^{-1} \sum_{\substack{|\mathbf{u}| \in D_q \\ \eta \in D_z}} \eta |\mathbf{u}|^2 P_{H1}(\mathbf{u}, \eta), \quad (44) $$

where $D_q$ is the set of $\mathbf{u}$ values contributing to bin $q$ and $D_z$ is the set of $\eta$ values contributing to bin $z$. Alternatively, if the 21 cm signal is expected to be highly isotropic, then it is also

appropriate to perform bin averages by averaging over spherical shells of constant $|k|$ to obtain the spherically averaged $H1$ power spectrum. The bin-averaging schemes mentioned above are especially necessary at the beginning of iterations with a poor initial guess since bin averaging can smooth out any unreasonable values and highly suppress foreground contamination that is expected to significantly increase in strength toward low $k_\parallel$ and would result in a substantial overestimate of the $H1$ power spectrum. One has to keep using the bin-averaging process until the derived amplitude of the $H1$ signal in each bin converges to a reasonable range.

In addition, although the bin widths can be chosen arbitrarily, an appropriate choice can reduce the correlation between the band-power estimates while retaining accurate detection of the structure of the power spectrum. The minimum bin width can be approximated by the typical correlation length in Fourier space. For an interferometer, one may choose a bin width of perpendicular modes, $\Delta \mathbf{u}$, greater than the characteristic width of the Fourier transformed primary beam pattern $A(\mathbf{x})$, and a bin width of parallel modes, $\Delta \eta$, roughly greater than the inverse of the bandwidth.

The EM algorithm iterates until there is no significant change in the likelihood. The stopping criterion of the iterative process is here set in terms of the relative change of the bin-averaged $H1$ power spectrum in the last iterations, i.e., $|P_{H1}^{n+1} - P_{H1}^n|/|P_{H1}^n| < 10^{-3}$ for all of the bins, and a typical number of iterations is a few hundred.

## 2.5. The Initial Guess

Since the EM algorithm is a hill-climbing approach, the searching procedure may fail to reach a global maximum and instead converge on a local maximum if the likelihood functions are not convex. In this study, we adapt a commonly used strategy to solve this problem. We try many different initial values varied in reasonable ranges and choose the solution that has the highest converged likelihood value. After a careful investigation, we chose the following algorithm for initialization. For the initial value of $M$, we perform PCA analysis (de Oliveira-Costa et al. 2008) to obtain the $n_c$ eigenvectors associated with the largest eigenvalues for the frequency-frequency covariance matrix estimated by averaging over all of the pixels of the data, and use those vectors as the $n_c$ columns of the mixing matrix $M$. After the mixing matrix is initialized, the power spectra $C_j^i$ are chosen to be the corresponding diagonal elements in $\hat{\mathbf{x}}(\hat{\mathbf{x}})^\dagger$, which is solved using Equation (23) and setting $C^{-1} = 0$ such that there is no prior information about $C$; Equation (23) reduces to a standard map-making equation. After initializing $M$ and $C_j^i$, we assume as prior information about the $H1$ signal a flat power spectrum with an amplitude comparable to the noise level. The simulations show that the $H1$ signal reconstruction is quite insensitive to the initial guess for its power spectrum as long as it is not unreasonably large.

## 3. SIMULATIONS

We perform simulations to generate dirty sky map data cubes. These data cubes include realistic models for the cosmological 21 cm signal, several diffuse foreground components, and instrumental noise. As this is just the first
test of HIEMICA, we assume for simplicity that the primary beam is unity for all frequency channels and assume complete uv coverage. We leave detailed investigations about more realistic observations to future work.

Because an interferometer only measures temperature fluctuations around the mean and is insensitive to the mean value of the brightness temperature, we set the mean of the 3D H1 signal and the means of the foreground and noise maps at each frequency channel to zero.

The 3D dirty sky map is simulated in a box with 64³ pixels (i.e., 64 pixels per side), covering a 30° × 30° sky patch and spanning over 780–880 MHz at intervals of 1.59 MHz, corresponding to redshifts between 0.82 and 0.61. In comoving coordinates, this box corresponds to about 1341 × 1341 × 606 Mpc³ and the size of each pixel is about 21 × 21 × 9.5 Mpc³. The code was run using the best-fit cosmological parameters from the Planck measurements (Ade et al. 2014).

### 3.1. H1 Signal

In cosmology, the power spectrum is typically represented in ($k_{\perp}$, $k_{\parallel}$) comoving coordinates. If the observed frequency band is small enough (i.e., probing a small range in redshifts) and one uses the flat-sky approximation, then there is a linear mapping between these variables:

$$ \mathbf{u} = \frac{k_{\perp} D_c(z)}{2\pi}; \quad \eta \approx \frac{c(1+z)^2}{2\pi H_0 \nu_{21} E(z)} k_{\parallel}, \qquad (45) $$

where $E(z) \equiv \sqrt{\Omega_m(1+z)^3 + \Omega_{\Lambda}}$, $\nu_{21}$ is the rest frequency of the 21 cm line, $D_c$ is the transverse comoving distance, $z$ is the redshift of the observation, $H_0$ is the Hubble parameter, $c$ is the speed of light, and $\Omega_m$ and $\Omega_{\Lambda}$ are the normalized matter and dark energy density, respectively. The angular wavenumber $k_{\perp}$ and the parallel wavenumber $k_{\parallel}$ are the components of the wavenumber $\mathbf{k}$ perpendicular and parallel to the line-of-sight direction, respectively. Therefore, by inserting Equation (45) into Equation (13), we obtain the relation of the power spectrum defined under different coordinates:

$$ P_{H1}(k_{\perp}, k_{\parallel}) = \frac{c(1+z)^2}{H_0 \nu_{21} E(z)} P_{H1}(\mathbf{u}, \eta). \qquad (46) $$

The 21 cm brightness temperature and the corresponding 3D power spectrum can be written as

$$ P_{H1}(k) = (\bar{T}_{H1}(z))^2 b^2(k, z) D^2(z) P_{cdm}(k) \qquad (47) $$

$$ \bar{T}_{21}(z) \approx 0.084 \, \text{mK} \frac{(1+z)^2}{E(z)} \frac{\Omega_B}{0.044} \frac{f_{H1}(z)}{0.01}, \qquad (48) $$

where $b(k, z)$ is the bias parameter and $P_{cdm}$ is the cold dark matter power spectrum at the present day. $\Omega_B$ is the baryon density fraction. $D(z)$ is the growth factor for dark matter perturbations defined such that $D(0) = 1$. For the purpose of this paper, we assume $b = 1$ over redshift and scale. For simplicity, we also neglect the effects of redshift-space distortions caused by the peculiar velocities of H1 clouds and galaxies in the H1 power spectrum since the ICA-based approach is insensitive to the detailed shape of the power spectrum. The H1 mass fraction is expected to increase with redshift and we assume a linear dependence: $f_{H1}(z) = 0.008(1+z)$. The matter power spectrum $P_{cdm}$ with

<table><caption>**Table 1**

Foreground Model Parameters for the Angular Power Spectrum $C_{\ell}(\nu, \nu')$ Used in Equation (51)</caption><thead><tr><th></th><th>A (mK<sup>2</sup>)</th><th>β</th><th>α</th><th>ξ</th><th>Δα</th></tr></thead><tbody><tr><td>Extragalactic point sources</td><td>57.0</td><td>1.1</td><td>2.07</td><td>1.0</td><td>0.2</td></tr><tr><td>Extragalactic free-free</td><td>0.014</td><td>1.0</td><td>2.10</td><td>35</td><td>0.03</td></tr><tr><td>Galactic synchrotron</td><td>700</td><td>2.4</td><td>2.80</td><td>4.0</td><td>0.15</td></tr><tr><td>Galactic free-free</td><td>0.088</td><td>3.0</td><td>2.15</td><td>35</td><td>0.03</td></tr></tbody></table>

the BAO can be parametrized according to the following simple empirical formula:

$$ \frac{P_{cdm}(k_{\perp}, k_{\parallel})}{P_{\text{ref}}(k_{\perp}, k_{\parallel})} = 1 + \qquad (49) $$

$$ A k \exp(-(k/\tau)^{\alpha}) \sin\left(2\pi\sqrt{\frac{k_{\perp}^2}{k_{\text{BAO}\perp}^2} + \frac{k_{\parallel}^2}{k_{\text{BAO}\parallel}^2}}\right), \qquad (50) $$

where $k = \sqrt{k_{\perp}^2 + k_{\parallel}^2}$, and the parameters $A$, $\alpha$, and $\tau$ are adjusted according to the formula presented in Eisenstein & Hu (1998). $P_{\text{ref}}(k)$ is the smooth “no-wiggles” power spectrum at $z=0$, which can be computed from the fitting formula given by Eisenstein & Hu (1998).

The parameters $k_{\text{BAO}\perp}$ and $k_{\text{BAO}\parallel}$ are the sinusoidal scales in the radial and transverse directions in $k$ space. We choose the following values for these parameters used in this paper: $A = 1.0$, $\tau = 0.1 \, h \, \text{Mpc}^{-1}$, $\alpha = 1.4$, and $k_{\text{BAO}\perp} = k_{\text{BAO}\parallel} = 0.060 \, h \, \text{Mpc}^{-1}$.

### 3.2. Foregrounds

We model the foregrounds as isotropic random Gaussian fields described by angular power spectra $C_{\ell}(\nu, \nu')$ based on Santos et al. (2005). In this paper, we assume that bright resolved point sources have been removed accurately and consider four dominant diffuse components: Galactic synchrotron emission, Galactic and extragalactic free-free emission, and extragalactic radio point sources. The angular power spectrum of each source takes the following generic form:

$$ C_f(\ell, \nu, \nu') = A \left(\frac{1000}{\ell}\right)^{\beta} \left(\frac{\nu_f^2}{\nu\nu'}\right)^{2\alpha} \exp\left(-\frac{\ln^2(\nu/\nu')}{2\xi^2}\right), \qquad (51) $$

where $\nu_f$ is the reference frequency with $\nu_f = 130$ MHz. We list the parameters of the foreground models used in this paper in Table 1. Based on such models, Figure 1 shows the angular power spectrum of each foreground contribution for one realization. The corresponding maps are shown in Figure 2.

Following Davies et al. (2006), the frequency index $\alpha$ of real-world foregrounds varies slightly across the sky. For example, in the case of synchrotron emission, the direction-dependent spectral index reflects variations of the relativistic electron density and Galactic magnetic field. In this paper, for each foreground component, we assume the indices in different directions to be Gaussian distributed with a mean of $\alpha$ and an rms of $\Delta\alpha$ as shown in Table 1, consistent with model parameters chosen in Liu & Tegmark (2012).

### 3.3. Noise

For each frequency channel, we assume Gaussian noise in the uv-plane with zero mean and that the noise maps at two
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Angular power spectra derived from the simulated maps of the 21 cm signal (blackthick) and foregrounds at $z = 0.71$ ($\nu = 830$ MHz) with a bandwidth of $\Delta\nu = 1.59$ MHz. The foreground components include Galactic synchrotron radiation (red-solid), extragalactic point sources (blue-dashed), and Galactic free-free (green-dotted) and extragalactic free-free (brown long-dashed).</figcaption></figure>

different frequencies are uncorrelated. We make noise realizations in the uv-plane for a given noise power spectrum for each frequency. The thermal noise of the measurement in units of brightness temperature is given by (McQuinn et al. 2006; Ansari et al. 2012a)

$$ P_N(\mathbf{u}, \nu) = \langle |\tilde{n}(\mathbf{u}, \nu)|^2 \rangle = \left(\frac{\lambda^2}{A_e}\right)^2 \frac{T_{\text{sys}}^2}{\Delta\nu \Delta t n_b} \\ \approx 1.9 \times 10^{-7} \, \text{mK}^2(1 + z_\nu)^4 \left(\frac{T_{\text{sys}}}{50 \, \text{K}}\right)^2 \\ \times \left(\frac{10 \, \text{m}}{A_e}\right) \left(\frac{10}{n_b}\right) \left(\frac{1 \, \text{MHz}}{\Delta\nu}\right) \left(\frac{30 \, \text{days}}{\Delta t}\right), \quad (52) $$

where $z_\nu$ is the corresponding redshift of the frequency $\nu$, $T_{\text{sys}}$ is the system temperature, $n_b$ is the total number of redundant baselines within that $u_k$ pixel with integration time $\Delta t$ for each baseline, $\Delta\nu$ is the width of the frequency channel, and $A_e$ is the effective area of each individual antenna in the array.

In order to keep the results as general as possible, we do not consider a specific configuration of the array; the main goal in this paper is not to discuss instrumental effects in 21 cm power spectrum measurements. Instead, we consider two simple scenarios in our simulation to illustrate the effectiveness of the presented foreground-cleaning approach. One scenario assumes a currently achievable noise level of $\sqrt{P_N} \approx 1 \times 10^{-3}$ mK for each pixel of the uv-plane, and the other one assumes a noise level of about $2.2 \times 10^{-4}$ mK for future measurements. The former has an averaged signal-to-noise ratio (S/N) over the whole data cube of about 1 while the latter is about 5.

## 4. APPLICATIONS TO SIMULATED DATA

In this section, we show our results obtained by applying the HIEMICA algorithm to the simulated data. There are two main findings, namely, the map reconstruction and the power spectrum recovery. As the fluctuations of the cosmological 21 cm signal are expected to be isotropic, we focus on the recovery of the averaged 3D power spectrum in spherical shells of constant $|k|$. However, redshift-space distortions will slightly break the spherical symmetry, and we will investigate the consequences of this in a future work. The estimated power spectra are obtained by averaging over the results from 10 independent realizations of the simulated data cubes and the associated statistical errors are obtained from their dispersions. Each data cube is for a fixed sky patch that includes the same simulated foregrounds combined with independent realizations of the instrumental noise and the H1 signal.

### 4.1. The Independent Components

Figure 3 shows the four WF ICs at $\nu = 830$ MHz found by HIEMICA. The WF components are calculated by Equation (23), and the corresponding map of the $i$th component at the frequency bin $\nu$ is obtained by multiplying such a component by the $(i, \nu)$th entry of the mixing matrix, $M_{\nu,i}$. As can be seen based on Figures 2–4, we find that only the first component is compatible with the input Galactic synchrotron contribution, while the others have no corresponding physical foregrounds and have unphysical spectral and spatial behaviors. These unphysical behaviors are essentially due to the fact that the realistic astrophysical components have significant cross-correlations while our ICA assumes mutual independence among the source signals. Thus, each separated component is just a mixture of all of the source signals, but the sum of these ICs is essentially the same as the total astrophysical signal, as seen in Figures 5(a) and (b). Meanwhile, the WF maps of the 21 cm signal (bottom panel of Figure 5) show the encouraging result that the 21 cm signal, even in the case of S/N = 1, is nearly perfectly recovered across the map by HIEMICA.

There is a slight difference between the low- and high-S/N cases (Figures 5(d) and (e)) at small scales, but we find no apparent difference at large scales. The recovered amplitudes of the fluctuations for S/N = 1 are smaller than those for S/N = 5 by several tens of percent, which is essentially due to the noise suppression in WF maps, that is, $\propto S/(S+N)$, where S/N tends to be low for high-k modes. However, the large-scale fluctuations of the 21 cm signal in both cases are much larger than the noise amplitudes, resulting in almost identical large-scale fluctuations.

### 4.2. Determination of the Number of ICs

When analyzing data with HIEMICA, we have to choose the number of Independent foreground Components (ICs) to be estimated. Our statistical framework for ICA can provide a rigorous determination of this number by using an empirical statistical approach, the Likelihood Ratio Test (LRT). The LRT provides an objective criterion for selecting between possible models by using the ratio of the likelihoods on the data. Although in practice evaluating likelihood functions is a hard task, we can instead apply HIEMICA several times with successively increasing $N_c$ until the recovered power spectra are essentially unchanged when more components are added into the analysis.
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Simulated diffuse foreground components for a 30° × 30° patch of sky at $\nu = 830$ MHz. We have subtracted the mean value (i.e., DC mode) to reflect an interferometric observation.</figcaption></figure>

The LRT statistic is given by $\text{LRT} = 2\ln \mathcal{L}_s - 2\ln \mathcal{L}_c$, in which $\ln \mathcal{L}_s$ and $\ln \mathcal{L}_c$ denote the logs of the likelihood functions for a relatively simpler model $s$ (the so-called null hypothesis) with fewer parameters and a more complex model $c$ (alternative hypothesis), respectively. Asymptotically, the LRT statistic follows a $\chi^2$ distribution with degrees of freedom equal to the difference in the number of free parameters between the two models. In our case, in the analysis of the data with 64 frequency channels and 30ℓ-bins, adding an additional foreground component would increase the degrees of freedom by 64 + 30 = 94.

The comparison of the fitting model between the models for $n_c = 1$ and 2 (i.e., labeled as 1 → 2 ICs in Table 2) and $n_c = 2$ and 3 (as 2 → 3 ICs) shows that, from the values of LRT, $\ln L$ increases significantly as the number of ICs increased up to 3. The significance values ($p \simeq 0$) for S/N of 1 and 5 indicate that only one or two foreground components would result in an unsuccessful fit and the data strongly prefer three foreground components. Moreover, especially in the case of high S/N, there is no statistically significant difference when using more components, given the larger $p$ values ($p > 0.01$). Here, all of the values of LRT and $p$ are the median over 10 independent simulations. Therefore, we recommend using three components in the reconstructed foreground model for our specific simulation parameters in Table 1. It is worth noting that, as mentioned by Alonso et al. (2015), the optimal number of ICs in fact strongly depends on the spectral smoothness of true foregrounds, characterized by a frequency correlation length $\xi$ as defined in Equation (51). In fact, a longer coherence length implies a smoother frequency spectrum for a physical foreground component. Consequently, a smaller value of $\xi$ may

require more ICs to successfully model such physical components and remove them accordingly.

### 4.3. 3D Power Spectrum Results

Figure 6 shows the rapid convergence of the resulting H1 power spectrum when successively increasing the number of ICs. It can be seen that the recovered power spectrum for an IC number of 2 is slightly different from those for IC numbers of 3 and 4, but the 3 and 4 components result in almost the same power spectrum, which is consistent with the varying trend in the LRT. This can also be confirmed by comparing the reconstructed WF ICs in Figure 3. As the amplitude of the fluctuations of the second component (see Figure 3(b)) is indeed about 20 times larger than the input 21 cm signal (see Figure 5(c)), the removal of the second component is therefore necessary. Furthermore, since the averaged amplitude of the fluctuations of the third component (see Figures 3(c) and 4) is about 10 times smaller than the 21 cm signal, the effects of such a component in signal recovery are expected to be less important but can still lead to notable improvements at large scales as the corresponding 3D Fourier modes are not significantly smaller than those of the H1 signal. Moreover, the amplitude of fluctuations of the fourth component (see Figure 3(d)) is about two times smaller than the third one, leading to rather small effects in signal reconstruction.

In the comparison between the estimated signal power spectra for S/N = 1 and 5 in Figure 6, we note that for $k \lesssim 0.2 \text{ Mpc}^{-1}$ over which the averaged S/N > 1 for both cases, the impact of noise at these scales tends to be negligibly small since the recovered power spectra follow the true spectrum very closely and there are no apparent differences
<figure><img src="image_4.png" /><figcaption>**Figure 3.** Wiener-filtered maps of independent components (ICs) reconstructed by applying HIEMICA to the simulated data cube with S/N = 5 for the same sky patch and frequency as in Figure 2. Note that we do not expect the recovered maps to explicitly correspond to the input maps since these ICs are assumed to be mutually independent while the true physical components may have correlations.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Figure 4.** Top panel: HIEMICA-derived coefficients of each column of the mixing matrix, representing the frequency dependence of the ICs, where the solid red line is for the first component, the dotted blue is for the second, the dashed green is for the third, and the long-dashed brown is for the fourth. When applying HIEMICA to simulated data cubes, we have assumed $n_c = 4$ and S/N = 5. Middle panel: the rms brightness temperatures of ICs calculated from the derived mixing matrix and their angular power spectra. Bottom panel: the rms brightness temperatures of the input physical foreground components including Galactic synchrotron radiation (solid red), extragalactic point sources (dashed blue), Galactic free-free (dotted green), and extragalactic free-free (long-dashed-brown). As mentioned in Figure 3, the temperature fluctuations of the ICs are not necessarily the same as the physical components.</figcaption></figure>

between the two cases. On the other hand, the recovered power spectrum for S/N = 1 begins gradually to overestimate the true one at $k > 0.2 \text{ Mpc}^{-1}$ by up to 12%, but for S/N = 5 the recovered spectrum is almost same as the true one with only 2% deviation. The additional power in the former case is primarily because the large noise level causes a slight underestimate of the true foreground contributions. The “excess power” is barely significant; the simulation truth is still contained in the $1\sigma$ error bars.

In Figures 6(c) and (d), we zoom in on the region with BAO features to illustrate clearly the dependence of the power spectrum estimate on the number of ICs as well as on the noise amplitude. We can see that too few components, such as $n_c = 2$ in both cases of S/N = 1 and 5, result in an overestimate of the power spectrum at scales of $0.1 \lesssim k \lesssim 0.2 \text{ Mpc}^{-1}$, since those two components cannot completely describe the spectral and spatial properties of the foreground. As the number of ICs increases to 3 and 4, the resulting power spectra decrease and both rapidly converge to the true values. As the S/N drops by a factor of 5, the measurement errors become larger since such errors can be approximately estimated by $\sigma(k) \propto \sqrt{P_{H1}(k) + P_{N}(k)}$. More importantly, one can see that if $n_c = 3$ or 4, the recovered $H_1$ spectrum is well within the $2\sigma$ error bars of the estimates.

Note that we have applied a correction to convert the prediction intervals on the underlying power spectrum into error bars that approximate Bayesian power spectrum inference in the limit of low noise and small numbers of modes. The resulting upper and lower error bars are asymmetric. The detailed analytical estimates of such a correction are presented in the appendix. As a consequence, compared with the
<figure><img src="image_6.png" /><figcaption>**Figure 5.** Same as Figure 3, but for the total foregrounds and the 21 cm signal to clearly illustrate that a successful separation was achieved by the HIEMICA-cleaning process. There is an underestimate of the 21 cm signal in the case of S/N = 1, which is caused by noise contamination as a Wiener-filtered map $\propto S/(S + N)$. An almost perfect recovery occurs over all scales for S/N = 5.</figcaption></figure>

<table><caption>**Table 2**
Comparison of Models by the Likelihood Ratio Test (LRT) for the Number of Independent Foreground Components</caption><thead><tr><th rowspan="2">Model</th><th rowspan="2">Δ dof</th><th>LRT</th><th>LRT</th><th rowspan="2">p-value (S/N = 1)</th><th rowspan="2">p-value (S/N = 5)</th></tr><tr><th>(S/N = 1)</th><th>(S/N = 5)</th></tr></thead><tbody><tr><td>1 → 2 ICs</td><td>94</td><td>17558</td><td>20872</td><td>0<sup>a</sup></td><td>0<sup>a</sup></td></tr><tr><td>2 → 3 ICs</td><td>94</td><td>196</td><td>280</td><td>&lt;10<sup>-6a</sup></td><td>&lt;10<sup>-6a</sup></td></tr><tr><td>3 → 4 ICs</td><td>94</td><td>120</td><td>125</td><td>0.036</td><td>0.018</td></tr><tr><td>4 → 5 ICs</td><td>94</td><td>96</td><td>108</td><td>0.42</td><td>0.15</td></tr></tbody></table>

**Note.**

<sup>a</sup> *p* < 0.01, typically used to assess the significance of LRT.

uncorrected errors, the upper error bars increased by about 60%, 40%, and 20%, respectively, in the lowest three *k*-bins and the associated lower error bars decreased by about 40%, 15%, and 10%. The differences are less than a few percent in the higher *k*-bins due to the large number of modes contained within those bins.

We can also see that the recovered cosmological signal decreases approximately monotonically with increasing number of ICs, especially for the low-*k* modes ($k \lesssim 0.1$ Mpc$^{-1}$). This result may be interpreted as follows: (1) since each recovered foreground component is a mixture of all the source signals, including the cosmological signal, a small amount of cosmological signal can leak into the reconstructed ICs, strongly biasing and reducing the amplitude of the derived cosmological signal; (2) the angular power spectrum of the dominant foreground component, Galactic synchrotron emission, has the largest amplitude for small *l* (see Figure 1),

resulting in a strong contamination in the cosmological signal at large scales in the transverse direction while increasing the risk of leakage; (3) the performance of the ICs-based reconstruction also depends on the number of measured Fourier modes in each *k*-bin and, in the simulated data, the lowest and the second lowest *k*-bins contain about 500 and 50 times fewer modes than those in the high-*k* bins, respectively.

As discussed before, using the LRT, the three-component foreground model used in HIEMICA gives rise to the successful recovery of the H1 power spectrum. In order to evaluate the recovery performance, we present the derived power spectrum in Figure 7 compared with that from the existing commonly used PCA technique. We can see that after applying PCA to subtract the three dominant foreground components from the data cube, the residual apparently underestimates the H1 power spectrum at large scales, which seems to be an over-subtraction of the foregrounds. Compared with the input ones, the amplitudes in the first two *k*-bins are smaller by factors of 14.1 and 3.7 if S/N = 1, and by factors of 16.1 and 4.1 if S/N = 5, respectively, whereas the ratios of the HIEMICA-derived power spectra and the input spectra in those two bins are only about 91%-92% if S/N = 1 and 91%-94% if S/N = 5, respectively. In addition, the recovered amplitudes in those two *k*-bins are consistent with the input values within $2\sigma$. The HIEMICA method thus does a much better job at foreground removal and cosmological signal recovery at large scales.

More promisingly, we find that the HIEMICA approach provides an unbiased estimate of the H1 power spectrum for the high-*k* modes with $k \gtrsim 0.1$ Mpc$^{-1}$ where S/N tends to be low, whereas the power spectrum of the PCA-derived residual
<figure><img src="image_7.png" /><figcaption>**Figure 6.** Dependence of the 21 cm signal recovery on the number of independent components, for data with S/N = 1 (left) and 5 (right). Upper panels: we show the spherically averaged three-dimensional power spectra of the simulated 21 cm signal (black), noise (dotted black), reconstructed 21 cm signal for the HIEMICA algorithm run with the assumption of $n_c = 2$ (dotted red), 3 (dashed green), and 4 (long-dashed blue), respectively. Vertical bars indicate the $1\sigma$ errors estimated from 10 realizations and the horizontal bars the bin width of $\Delta k = 0.0129$. The statistical uncertainties are typically smaller than the symbol sizes and are mostly invisible. Lower panels: same as upper panels, but showing the H1 power spectrum divided by a smoothed spectrum without baryon acoustic oscillations to single out those oscillations.</figcaption></figure>

in the noise-dominant regime seems to deviate strongly from the true signal, resulting in an overestimate of H1 measurements.

Since the isotropic H1 signal has the same statistical property as that of the assumed instrumental noise, commonly used ICA-based approaches have trouble breaking the degeneracy, leaving residuals consisting of the reconstructed 21 cm signal, noise, and fitting errors. As proposed in Chapman et al. (2012)

and Bonaldi & Brown (2015), after foreground cleaning by CCA or FASTICA, one has to manually subtract the noise power spectrum from that of the residual data cube and obtain a comparable amplitude with the true signal. However, those previous studies show that such a scheme can cause an overfitting problem and lead to some negative or zero-closed amplitudes of the residual power spectrum in noise-dominated k modes, implying that a leakage of H1 signal into the
<figure><img src="image_8.png" /><figcaption>**Figure 7.** Same as Figure 6, but for a comparison of 21 cm power spectrum recovery based on the PCA technique (dotted magenta) and the HIEMICA approach (long-dashed gren), for the data with S/N = 1 (left) and S/N = 5 (right). The spherically averaged three-dimensional power spectra of the simulated 21 cm signal (solid black), noise (dotted black), HIEMICA-derived (assumed $n_c = 3$) 21 cm signal, and PCA-derived residuals from projecting out the first three dominant eigen components based on the frequency-frequency covariance of data.</figcaption></figure>

subtracted foreground components is likely to be serious in regions with low S/N. Some other schemes, such as weighting maps with inverse noise variance to build the frequency covariance used in PCA (Alonso et al. 2015), could mitigate such noise contamination.

In contrast to the above-mentioned standard PCA- and ICA-based approaches, the proposed HIEMICA method in fact provides an unbiased and efficient estimate of Gaussian components, derived from the maximum-likelihood principle. If the noise covariance matrix *N* is exactly known, then this method can optimally estimate the underlying 21 cm signal and foreground covariance matrices. It avoids any artificial effects that can occur when manually subtracting the noise power spectrum from the residuals after the foreground-cleaning processes. Our simulations confirm this, showing that the ensemble-averaged recovered H I power spectrum converges to the true value even in low S/N regions where the noise is about an order of magnitude higher.

We note that our semi-blind approach requires exact knowledge of the noise in order to infer the 21 cm power spectrum accurately; an incorrect noise covariance matrix can bias the signal estimate. However, a relatively lower accuracy in modeling the noise covariance matrix may be sufficient for the power spectrum estimation if the uncertainty in the noise level is not greater than the signal. In general, the corresponding impacts on the estimated parameters can be evaluated by a Monte-Carlo (MC) approach. But qualitatively, the level of the uncertainty in the estimated 21 cm power spectrum is about of the same order of magnitude as that of the noise. A simple way to understand this is to consider that in the absence of

instrumental effects, if the noise is assumed to be uncorrelated, the value of $P_{H1}$ which maximizes the likelihood function can be approximately estimated by $P_{H1} \approx P_D - P_F - P_N$, where $P_D$, $P_F$, and $P_N$ are the power spectra of the data points, the foreground, and noise, respectively. As the total amplitude of the foreground fluctuations is significantly larger than that of the noise fluctuations, one can expect that the estimate of the foreground power spectrum roughly remains unchanged when slightly varying the input noise amplitude. Thus assuming an input noise with an amplitude either greater or smaller than the true one would lead to either an underestimate or an overestimate of the 21 cm power spectrum, i.e., $\Delta P_{H1} \approx -\Delta P_N$. In order to check the above estimate, we test the effect of using an input noise power spectrum with 10%–50% larger or smaller amplitudes than the true one and find that the resulting changes $\Delta P_{H1}$ are indeed consistent with our expected variations ($-\Delta P_N$) within a factor of 2 over all scales. In addition, incorrect noise estimates would also bias the estimate for the PCA method, since the debiasing MC simulations would give the wrong results if the noise level were misspecified.

In more realistic cases, an experiment may exhibit a low level of cross-correlated Gaussian or non-Gaussian noise, both of which are unknown. In the former case, similar to the above qualitative estimate, we expect that the small level of correlation will slightly bias the signal estimate but not significantly affect it. This is because, as mentioned above, (1) the H I power spectrum is approximately determined by the difference among the diagonal components of the data point, the foreground, and noise model, not by their off-diagonal
ones, and (2) the estimate of the foreground is unchanged or insensitive to the noise model. In the latter case, as the variance of the non-Gaussian noise is non-negative, the H1 power spectrum is thus likely to be slightly overestimated. We will investigate the impact of noise misspecification in a future work.

Also, this blind component separation algorithm makes no assumptions about the spectral information of the foreground and can naturally include all of the instrumental effects so as to be able to coherently estimate all of the source components and the associated errors. More detailed comparisons between HIEMICA and other approaches for data in the presence of realistic instrumental effects will be performed in future work.

## 5. CONCLUSIONS

In this paper, we present a non-parametric source separation algorithm, HIEMICA, which is an extension to the 3D version of the SMICA method proposed by Snoussi et al. (2002) and Delabrouille et al. (2003) for source identification in noisy mixtures of CMB maps. The HIEMICA algorithm is a fully Bayesian framework used to infer the 3D power spectrum and maps of the underlying H1 signal and the spatial power spectra and the frequency dependence (i.e., mixing matrix) of uncorrelated foreground components. We adapt the EM algorithm to efficiently maximize the likelihood function for unknown parameter estimates. As the statistical properties of the cosmological signal are significantly different from those of each astrophysical source, the spectral matching method can, in principle, blindly separate the H1 signal from highly foreground-contaminated maps.

The simulations show that HIEMICA is able to successfully reconstruct the H1 3D power spectrum across all scales and is much more robust than the PCA approach, which removes some of the 21 cm signal at large scales ($k \lesssim 0.1$ Mpc$^{-1}$) and overestimates the H1 power spectrum in the noise-dominated regions where noise leaks into the estimated signal. To evaluate the impact of the number of ICs in signal recovery, we adapt the LRT method to rigorously assess the likelihood fit so as to determine the optimal number of foreground components used in the HIEMICA analysis. We find for the number of ICs higher than 3, a number comparable to that of the simulated astrophysical components, the reconstruction does not gain significant improvements and converges toward a stable result consistent with the input H1 power spectrum.

Although our results are quite promising, the simulated observations in this study are idealized with no instrumental effects. As instrument effects such as the frequency-dependent uv sampling and primary beam can significantly complicate the foreground removal and bias the signal estimate, it is important to test HIEMICA in more detail in the context of specific experiment set-ups in future work.

We acknowledge use of the open-source PETSc library (-Balay et al. 1997, 2014a, 2014b) and FFTW (Frigo & Johnson 2005). L.Z. thanks Flavien Vansyngel for useful comments on early drafts. L.Z. and P.T. acknowledge support from NSF awards IIS 1250720 (“Big Data” program) and AST 1216525. P.M.S. is supported by the INFN IS PD51 “Indark.” We thank an anonymous referee for valuable suggestions which helped us to significantly improve this paper.

## APPENDIX

Here, we present a derivation of an approximate expression for the posterior distribution of the signal power spectrum in terms of the sampling distribution of the estimator so as to calculate Bayesian credible intervals of the inferred parameter. The sampling distribution of the estimator in our case can be derived from an MC distribution of the Maximum Likelihood Estimator (MLE) of the signal power spectrum P.

The power spectrum is simply the variance of the Fourier mode coefficients. For simplicity, if we assume that all of the Fourier modes are independent and the variance of all the Fourier mode coefficients within the bin is the same, the posterior distribution of P given data $d_i$ for $i = 1, \dots, n$ can be expressed in an ideal noise-free case as

$$p(P|d) \propto \mathcal{L}(d|P)p(P) \propto P^{-n/2}e^{-\frac{1}{2}\sum_{i=1}^{n} d_i^2/P} P^{-1}, \quad (53)$$

where the non-informative Jeffreys prior (Jeffreys 1961), $p(P) = P^{-1}$, has been applied on the power spectrum. As we know, the MLE $\hat{P} = \frac{1}{n} \sum_{i=1}^{n} d_i^2$ is unbiased over data sets drawn from

$$\langle \hat{P} \rangle = P \quad (54)$$

and has variance

$$\langle (\Delta \hat{P})^2 \rangle = \frac{2}{n} P^2. \quad (55)$$

If we want to construct an approximate Bayesian posterior from MC samples from the sampling density of the MLE $\hat{P}$ and the effective number of degrees of freedom is unknown, then we can determine it from the mean and variance of $\hat{P}$ using

$$\tilde{n} = 2 \frac{\langle \tilde{P} \rangle}{\langle (\Delta \hat{P})^2 \rangle}. \quad (56)$$

Assuming that noise is subdominant, inserting Equation (56) into Equation (53), we can use $\hat{P}$ and the effective number of degrees of freedom to approximate the posterior for P as

$$p(P|d) \propto e^{-\frac{1}{2}\tilde{n}\hat{P}/P} / P^{\tilde{n}+2/2}, \quad (57)$$

from which the asymmetric credible interval can be derived.

## REFERENCES

Ade, P. A. R., Aghanim, N., Armitage-Caplan, C., et al. for the Planck Collaboration 2014, A&A, 571, A16

Alonso, D., Bull, P., Ferreira, P. G., & Santos, M. G. 2015, MNRAS, 447, 400

Anderson, L., Aubourg, E., Bailey, S., et al. 2013, MNRAS, 427, 3435

Ansari, R., Campagne, J. E., Colom, P., et al. 2012a, A&A, 540, A129

Ansari, R., Campagne, J.-E., Colom, P., et al. 2012b, CRPhy, 13, 46

Balay, S., Abhyankar, S., Adams, M. F., et al. 2014a, PETSc Users Manual, Tech. Rep. ANL-95/11—Revision 3.5 (Argonne National Laboratory)

Balay, S., Abhyankar, S., Adams, M. F., et al. 2014b, PETSc Web page, http://www.mcs.anl.gov/petsc

Balay, S., Gropp, W. D., McInnes, L. C., & Smith, B. F. 1997, in Modern Software Tools in Scientific Computing, ed. E. Arge, A. M. Bruaset, & H. P. Langtangen (Basel: Birkhäuser)

Bonaldi, A., & Brown, M. L. 2015, MNRAS, 447, 1973

Bowman, J. D., Morales, M. F., & Hewitt, J. N. 2006, ApJ, 638, 20

Bowman, J. D., Morales, M. F., & Hewitt, J. N. 2009, ApJ, 695, 183

Bull, P., Ferreira, P. G., Patel, P., & Santos, M. G. 2015, ApJ, 803, 21

Cardoso, J.-F., Martin, M., Delabrouille, J., Betoule, M., & Patanchon, G. 2008, ISTSP, 2, 735

Chang, T., Pen, U., Bandura, K., & Peterson, J. B. 2010, Natur, 466, 463
Chang, T.-C., Pen, U.-L., Peterson, J. B., & McDonald, P. 2008, *PhRvL*, **100**, 091303

Chapman, E., Abdalla, F. B., Harker, G., et al. 2012, *MNRAS*, **423**, 2518

Chen, X. 2012, *IJMPS*, **12**, 256

Chen, X.-L., & Miralda-Escude, J. 2004, *ApJ*, **602**, 1

Colless, M., Peterson, B. A., Jackson, C., et al. 2003, *astro-ph/0306581*

Cooray, A., & Furlanetto, S. R. 2004, *ApJL*, **606**, L5

Datta, A., Bowman, J. D., & Carilli, C. L. 2010, *ApJ*, **724**, 526

Davies, R. D., Dickinson, C., Banday, A. J., et al. 2006, *MNRAS*, **370**, 1125

Delabrouille, J., Cardoso, J.-F., & Patanchon, G. 2003, *MNRAS*, **346**, 1089

de Oliveira-Costa, A., Tegmark, M., Gaensler, B. M., et al. 2008, *MNRAS*, **388**, 247

Dillon, J. S., Liu, A., & Tegmark, M. 2013, *PhRvD*, **87**, 043005

Di Matteo, T., Ciardi, B., & Miniati, F. 2004, *MNRAS*, **355**, 1053

Di Matteo, T., Perna, R., Abel, T., & Rees, M. J. 2002, *ApJ*, **564**, 576

Drinkwater, M. J., Jurek, R. J., Blake, C., et al. 2010, *MNRAS*, **401**, 1429

Eisenstein, D. J., & Hu, W. 1998, *ApJ*, **496**, 605

Fessler, J. A., & Sutton, B. P. 2003, *ITSP*, **51**, 560

Frigo, M., & Johnson, S. G. 2005, *IEEEP*, **93**, 216

Furlanetto, S., Oh, S. P., & Briggs, F. 2006, *PhR*, **433**, 181

Gleser, L., Nusser, A., & Benson, A. J. 2008, *MNRAS*, **391**, 383

Harker, G., Zaroubi, S., Bernardi, G., et al. 2009a, *MNRAS*, **397**, 1138

Harker, G., Zaroubi, S., Bernardi, G., et al. 2010, *MNRAS*, **405**, 2492

Harker, G. J. A., Zaroubi, S., Thomas, R. M., et al. 2009b, *MNRAS*, **393**, 1449

Hinshaw, G., Larson, D., Komatsu, E., et al. 2013, *ApJS*, **208**, 19

Hyvärinen, A. 1999, Fast and Robust Fixed-Point Algorithms for Independent Component Analysis, *IEEE Trans. On Neural Networks*, **10**, 626

Hyvärinen, A., & Oja, E. 2000, *NN*, **13**, 411

Jeffreys, H. 1961, Theory of Probability (3rd ed.; Oxford: Oxford Univ. Press)

Jelić, V., Zaroubi, S., Labropoulos, P., et al. 2008, *MNRAS*, **389**, 1319

Jewell, J., Levin, S., & Anderson, C. H. 2004, *ApJ*, **609**, 1

Kinney, W. H., Kolb, E. W., Melchiorri, A., & Riotto, A. 2006, *PhRvD*, **74**, 023502

Liu, A., Parsons, A. R., & Trott, C. M. 2014a, *PhRvD*, **90**, 023018

Liu, A., Parsons, A. R., & Trott, C. M. 2014b, *PhRvD*, **90**, 023019

Liu, A., & Tegmark, M. 2012, *MNRAS*, **419**, 3491

Liu, A., Tegmark, M., Bowman, J., Hewitt, J., & Zaldarriaga, M. 2009a, *MNRAS*, **398**, 401

Liu, A., Tegmark, M., & Zaldarriaga, M. 2009b, *MNRAS*, **394**, 1575

Loeb, A., & Wyithe, J. S. B. 2008, *PhRvL*, **100**, 161301

Loeb, A., & Zaldarriaga, M. 2004, *PhRvL*, **92**, 211301

Madau, P., Meiksin, A., & Rees, M. J. 1997, *ApJ*, **475**, 429

Maldacena, J. M. 2003, *JHEP*, **0305**, 013

Mao, Y., Tegmark, M., McQuinn, M., Zaldarriaga, M., & Zahn, O. 2008, *PhRvD*, **D78**, 023529

Masui, K. W., Switzer, E. R., Banavar, N., et al. 2013, *ApJL*, **763**, L20

McQuinn, M., Zahn, O., Zaldarriaga, M., Hernquist, L., & Furlanetto, S. R. 2006, *ApJ*, **653**, 815

Morales, M. F., Bowman, J. D., & Hewitt, J. N. 2006, *ApJ*, **648**, 767

Morales, M. F., Hazelton, B., Sullivan, I., & Beardsley, A. 2012, *ApJ*, **752**, 137

Morales, M. F., & Hewitt, J. 2004, *ApJ*, **615**, 7

Morales, M. F., & Wyithe, J. S. B. 2010, *ARA&A*, **48**, 127

Nan, R., Li, D., Jin, C., et al. 2011, *IJMPD*, **20**, 989

Oh, S. P., & Mack, K. J. 2003, *MNRAS*, **346**, 871

Parsons, A. R., Pober, J. C., Aguirre, J. E., et al. 2012, *ApJ*, **756**, 165

Peiris, H. V., Komatsu, E., Verde, L., et al. 2003, *ApJS*, **148**, 213

Pen, U.-L., Staveley-Smith, L., Peterson, J. B., & Chang, T.-C. 2009, *MNRAS*, **394**, L6

Peterson, J. B., Bandura, K., & Penx, U. L. 2006, in Proc. 41st Recontres de Moriond, 2006 Contents and structure of the universe, 283

Petrovic, N., & Oh, S. P. 2011, *MNRAS*, **413**, 2103

Planck Collaboration, Ade, P. A. R., Aghanim, N., et al. 2014, *A&A*, **571**, A12

Pober, J. C., Liu, A., Dillon, J. S., et al. 2014, *ApJ*, **782**, 66

Pober, J. C., Parsons, A. R., Aguirre, J. E., et al. 2013a, *ApJL*, **768**, L36

Pober, J. C., Parsons, A. R., DeBoer, D. R., et al. 2013b, *AJ*, **145**, 65

Pritchard, J. R., & Loeb, A. 2012, *RPPh*, **75**, 086901

Santos, M. G., Cooray, A., & Knox, L. 2005, *ApJ*, **625**, 575

Seo, H.-J., & Eisenstein, D. J. 2003, *ApJ*, **598**, 720

Shaw, J. R., Sigurdson, K., Pen, U.-L., Stebbins, A., & Sitwell, M. 2014, *ApJ*, **781**, 57

Smoot, G. F., Bennett, C. L., Kogut, A., et al. 1992, *ApJL*, **396**, L1

Snoussi, H., Patanchon, G., Macías-Pérez, J. F., Mohammad-Djafari, A., & Delabrouille, J. 2002, in AIP Conf. Series, Vol. 617, Bayesian Inference and Maximum Entropy Methods in Science and Engineering, ed. R. L. Fry (New York: AIP) p. 125

Sutter, P. M., Wandelt, B. D., & Malu, S. S. 2012, *ApJS*, **202**, 9

Switzer, E. R., Masui, K. W., Bandura, K., et al. 2013, *MNRAS*, **434**, L46

Trott, C. M., Wayth, R. B., & Tingay, S. J. 2012, *ApJ*, **757**, 101

Vansyngel, F., Wandelt, B. D., Cardoso, J.-F., & Benabed, K. 2014, *arXiv:1409.0858*

Vedantham, H., Udaya Shankar, N., & Subrahmanyan, R. 2012, *ApJ*, **745**, 176

Wandelt, B. D., Larson, D. L., & Lakshminarayanan, A. 2004, *PhRvD*, **70**, 083511

Wang, X., Chen, X., Zheng, Z., et al. 2009, *MNRAS*, **394**, 1775

Wang, X., Tegmark, M., Santos, M. G., & Knox, L. 2006, *ApJ*, **650**, 529

Wolz, L., Abdalla, F. B., Blake, C., et al. 2014, *MNRAS*, **441**, 3271

Wyithe, S., Loeb, A., & Geil, P. 2008, *MNRAS*, **383**, 1195

Xu, Y., Wang, X., & Chen, X. 2015, *ApJ*, **798**, 40

York, D. G., Adelman, J., Anderson, J. E., Jr., et al. 2000, *AJ*, **120**, 1579

Zaldarriaga, M., Furlanetto, S. R., & Hernquist, L. 2004, *ApJ*, **608**, 622