A TECHNIQUE FOR FOREGROUND SUBTRACTION IN REDSHIFTED 21 cm OBSERVATIONS

A promising emerging field in cosmology is the proposed use of 21 cm emission from neutral hydrogen to trace the evolution of structure in the universe from redshift 0 to 50. At redshifts up to about 3, by mapping the three-dimensional intensity field with 10 arcmin resolution, it may be possible to precisely measure the expansion history throughout the transition from deceleration to acceleration (Chang et al. 2008; Wyithe et al. 2008; Loeb & Wyithe 2008; Morales & Wyithe 2010). For redshifts near 10, several instruments are under construction to search for 21 cm emission from the Epoch of Reionization (EoR; Morales & Hewitt 2004; Zaroubi & Silk 2005; Backer et al. 2007). Some authors even suggest that observations at redshifts well above 10 may provide very sharp tests of the world model and the composition of the cosmic fluid (Loeb & Zaldarriaga 2004).

Most of the observational effort so far has focused on searching for the EoR signal. The first goal of current EoR programs is to pin down the redshift of reionization. Later measurements will address questions such as the following. What were the first sources of ionizing radiation? What is the power spectrum of the H i structure, and how did it evolve? The first generation of 21 cm EoR observatories may begin to attack this last question by observing the emission of the neutral hydrogen itself. The interpretation of the data in terms of pinning down the UV sources is likely to be much more difficult.

Even at the end of the EoR, there is still enough H i (a few percent of all the hydrogen) in the form of self-shielded clumps to be visible at a brightness temperature contrast of ΔT ∼ 1 mK. We have some information about these structures since they are seen via UV absorption as damped Lyα regions (Wolfe et al. 2005). Observation of 21 cm fluctuations from such systems provides a tool for understanding structure formation in the post-reionization era. These H i regions presumably trace the dark matter and hence can be used to study the evolution of the matter power spectrum. Measurement of the cosmic power spectrum as a function of redshift over the cosmic volume 0.5 < z < 6 could provide cosmic-variance limited determination of cosmological parameters such as the equation of state of the dark energy and the neutrino mass (Loeb & Wyithe 2008).

Multiple lines of observational evidence now indicate that dark energy accounts for ∼70% of the energy density in the universe, but so far we have few clues as to the physics underlying this phenomenon. There are a host of dark energy models and these can be distinguished by their equations of state. A promising approach is to infer the equation of state from the distance–redshift relation in the redshift range of 0 to about 2. In particular, baryon-acoustic oscillations (BAO) can provide a standard ruler for distance determinations. The BAO "wiggles" in the matter power spectrum can be measured efficiently using 21 cm intensity mapping (Chang et al. 2008; Ansari et al. 2008, 2011; Seo et al. 2010).

To achieve these scientific goals a wide variety of innovative 21 cm telescope designs have been proposed, some prototypes built, and some telescopes completed (see Morales & Wyithe 2010, for a review of these instruments and their scientific potential). However, there are many development tasks needed to advance the state of the art of 21 cm cosmology, including calibration strategies, low-noise amplifier development, correlator design, examination of array layouts, radio-quiet site development, and foreground removal.

One of the main difficulties of all 21 cm telescopes is to remove foregrounds which are several orders of magnitude more intense than the H i signal (Morales et al. 2006). The basic idea is that most foregrounds have a smooth power-law frequency dependence in contrast with the H i signal. The exception is Galactic radio recombination lines which occur at known frequencies and can be excised. The main foregrounds are the Galactic synchrotron radiation and extragalactic point sources. Other sources like free–free electron emission (i.e., bremsstrahlung) are much less intense and have power-law spectra similar to the first two components.

Other foreground removal techniques have been proposed as well. Morales et al. (2006) review these techniques and show how to correct the errors they cause in the recovered power spectrum. Liu et al. (2009) show that performing foreground removal in Fourier space eliminates errors that arise from frequency-dependent beam patterns (mode mixing). Liu et al. (2011) have developed a technique for estimating the power spectrum in the presence of foregrounds without introducing noise or bias to the power spectrum.

Simulations by Bowman et al. (2009) for the Murchison Widefield Array remove the contribution from bright identified sources and then in each R.A.–decl. pixel fit and subtract a polynomial in ν. This two-step strategy allows for the observation of the H i signal at the time of reionization (z ≈ 8). Similar foreground removal simulations have been applied to candidate dark energy observatories; Ansari et al. (2011) show how foregrounds might be separated from the H i signal based on their differing spectral signatures as measured by an interferometric array of dishes or cylinder telescopes.

Below we discuss a different idea for removing foregrounds which is based on using the power spectrum. The physical basis of the idea is that at different frequencies the fluctuations of redshifted H i are not correlated, as we sample completely different regions of H i, while the foreground fluctuations are correlated as they arise from the same galactic inhomogeneities, e.g., turbulent fluctuations as discussed in Cho & Lazarian (2002, 2010). However, unlike Cho & Lazarian (2010) the suggested analysis does not require assumptions about the power spectrum of fluctuations.

In Section 2 we discuss the technique that we propose, in Section 3 the applicability of our technique, and in Section 4 we summarize our results.

2.1. Assumptions

2.2. Basic Idea

We can obtain an accurate shape of the foreground angular power spectrum in the following way. Suppose that we have multi-channel observation data in uv coordinates (left panel of Figure 1). We stack all the uv plane data into a single file (right panel of Figure 1). That is, we calculate ∑^N_{i = 1}T_νi(u, v). (If the observations are done in real space, we first obtain ∑^N_{i = 1}T_νi(θ_x, θ_y) and then we perform a Fourier transform.) In the stacked map, since the foreground signals in different channels are correlated, they are enhanced by a factor of ∼N. However, the noise and the 21 cm signals are enhanced by a factor of $\sqrt{N}$ because they contribute randomly. Therefore, the contrast between the foreground and the noise/21 cm angular spectra will be increased by ∼N times. In this way we can obtain an accurate shape of the foreground power spectrum. We define T_avg(u, v) as

$$\begin{equation} T_{{\rm avg}}(u,v)\equiv \frac{1}{N} \sum _{i=1}^{N} T_{\nu _i}(u,v). \end{equation} \tag{ 1 }$$

What, then, does T_avg(u, v) mean? If T_ν(u, v) varies slowly as a function of ν, we can show that the obtained foreground power spectrum is a good representation for the power spectrum at the central channel (ν = ν_{(N + 1)/2}). Suppose that N is an odd number. Then we can write

$$\begin{equation} T_{\nu _i}^{{\rm for}}(u,v)= T_{\nu _{(N+1)/2}}^{{\rm for}}(u,v)+ a \Delta \nu, \end{equation} \tag{ 2 }$$

where a is a constant, Δν = ν_i − ν_{(N + 1)/2} and the superscript "for" denotes the foreground. We assume $|a \Delta \nu | \ll |T_{\nu _i}^{{\rm for}} |$ and we ignore the second-order quantities. From Equations (1) and (2), we have $\frac{1}{N}\sum _{i=1}^{N} T_{\nu _i}(u,v) \approx T_{\nu _{(N+1)/2}}^{{\rm for}}(u,v),$ which results in

$$\begin{eqnarray} &&C_l^{{\rm avg}} \approx C_l^{\nu _{(N+1)/2}, {\rm for}}, \end{eqnarray} \tag{ 3 }$$

where C^avg_l and C^{ν_i, for}_l are the angular power spectra of T_avg(u, v) and $T_{\nu _i}^{{\rm for}}(u,v)$, respectively.

Zoom In Zoom Out Reset image size

The foreground power spectrum (≈C^avg_l), as well as the foreground map (≈T_avg(u, v)), can be used for other observations. For example, when a future low-noise single channel observation becomes available at ν = ν_{(N + 1)/2}, we can use C^avg_l or T_avg(u, v) to subtract the foreground. It is also possible to use them for future high angular resolution low-noise observations. Cho & Lazarian (2010) argued that when accurate foreground observations are available for a certain range of l, we can extrapolate the foreground spectrum to a range of l where no foreground information is available. Therefore, if a very high angular resolution single-frequency observation is performed in the future, we can extrapolate the foreground spectrum (i.e., C^avg_l) to higher values of l and extract the 21 cm information on the small scales.

2.3. Proposed Techniques

Case 1: same power spectrum shape for all channels. In this case, we can easily obtain C^{ν_i, 21 cm}_l. Note that, at this moment, we know only the shape (the l dependence) of C^{ν_i, for}_l. The amplitude can be different in each frequency channel and we do not know it yet. We can determine the amplitude by fitting or by the following method. Suppose that the value of

$$\begin{equation} \sum _{l} \big| \log \big(C_l^{\nu _i}\big) - \log \big(\alpha C_l^{{\rm avg}}\big)\big|^2 \end{equation} \tag{ 4 }$$

has a minimum value at α = α_i. Then we have

$$\begin{equation} C_l^{\nu _i, 21\, {\rm cm}} \approx C_l^{\nu _i} - \alpha _i C_l^{{\rm avg}}. \end{equation} \tag{ 5 }$$

Case 2: linear dependence of T(u, v) on ν. In this case, as we discussed in the previous subsection, we have $C_l^{{\rm avg}} =C_l^{\nu _{(N+1)/2, {\rm for}}}$. The 21 cm spectrum at ν = ν_i is

$$\begin{equation} C_l^{\nu _i, 21\, {\rm cm}}=C_l^{\nu _i}-C_l^{\nu _{i, {\rm for}}}. \end{equation} \tag{ 6 }$$

From Equation (2), we have

$$\begin{equation} C_l^{\nu _i, {\rm for}}=C_l^{\nu _{(N+1)/2}, {\rm for}} + A \Delta \nu, \end{equation} \tag{ 7 }$$

where A = Re[2∑_l/2πT^for_{ν(N + 1)/2}(u, v)a*]. Here "*" denotes the complex conjugate and Re[...] stands for the real part. Therefore, the average spectrum of the 21 cm signal over N channels is

$$\begin{equation} \frac{1}{N}\sum _{i=1}^{N} C_l^{\nu _i,21\, {\rm cm}} =\frac{1}{N} \left(\sum _{i=1}^{N} C_l^{\nu _i} \right) - C_l^{{\rm avg}}. \end{equation} \tag{ 8 }$$

Case 3: quadratic dependence of T(u, v) on ν. In this more general case, we can obtain the average 21 cm spectrum. Suppose that N = even for simplicity. Then we can write

$$\begin{equation} T_{\nu _i}^{{\rm for}}(u,v)=f_0+ a \Delta \nu + b \Delta \nu ^2, \end{equation} \tag{ 9 }$$

where f₀ is a visibility at ν₀ = (ν_N/2 + ν_{1 + N/2})/2 (not observed), a and b are constants, Δν = ν_i − ν₀ = [i − (N + 1)/2]δν, and δν = ν_{i + 1} − ν_i. We note that

$$\begin{eqnarray} | T_{{\rm avg}}(u,v)|^2 &=& \left| f_0 + \frac{b}{2N} \Delta \nu _0^2 \sum _{j=1}^{N/2} (2j-1)^2 \right|^2, \nonumber \\ &=& | f_0 |^2 + \frac{1}{N} \Delta \nu _0^2 \sum _{j=1}^{N/2} (2j-1)^2 {\rm Re} [ f_0 b^* ].\quad\quad \end{eqnarray} \tag{ 10 }$$

On the other hand, we have

$$\begin{eqnarray} && \sum _{i=1}^{N/2} |T_{\nu _i}(u,v)+T_{\nu _{N-i+1}}(u,v)|^2 \mbox{ } = \sum _{j=1}^{N/2} \left| 2f_0 + \frac{b}{2} (2j-1)^2 \delta \nu ^2 \right|^2 \nonumber\\\ms\ms && \quad\quad + \sum _{i=1}^{N/2} \big| T_{\nu _i}^{21\,{\rm cm}} + T_{\nu _{N-i+1}}^{21\,{\rm cm}} \big|^2 = 2N | f_0 |^2 \nonumber \\\ms\ms && \quad\quad + 2 \delta \nu ^2\sum _{j=1}^{N/2} (2j+1)^2 {\rm Re} [ f_0 b^*] + \sum _{i=1}^{N}\big|T_{\nu _i}^{21\,{\rm cm}}\big|^2, \end{eqnarray} \tag{ 11 }$$

where we dropped the cross-correlation terms between the 21 cm and the foreground signals (see Section 3). If we subtract Equation (10) times 2N from Equation (11), we have ∑^N_{i = 1}|T_νi^21 cm|². Therefore, the average spectrum of the 21 cm signal is

$$\begin{equation} \frac{1}{N}\sum _{i=1}^{N} C_l^{\nu _i,21\, {\rm cm}}=\frac{1}{N}\left(\sum _{i=1}^{N/2} D_l^{\nu _i} \right) - 2C_l^{{\rm avg}}, \end{equation} \tag{ 12 }$$

where Dνil is the angular spectrum of $T_{\nu _i}(u,v)+T_{\nu _{N-i+1}}(u,v)$.4

Our technique critically depends on the assumption of slow variation of the $T_{\nu _i}^{{\rm for}}(u,v)$ over frequency channels. This assumption may need further investigation. However, we believe that this is a good assumption due to the following reasons. First of all, we know that the major foreground, which is the synchrotron foreground, is caused by MHD turbulence (see the statistical description of fluctuations in Lazarian & Pogosyan 2012). The variations with frequency from an elementary synchrotron emitting volume are given by ∝n₀B²_⊥ν^{(1 − α)/2} where B_⊥ is the magnetic field strength perpendicular to the line of sight and the number density of cosmic rays has a power-law distribution n(E)dE = n₀E^−αdE ∼ n₀E⁻³dE. Therefore, we expect the spectral components to change as ∝(ν/ν₀)⁻¹. This dependence can be compensated over the frequency interval of Δν = (N − 1)δν and the residual emission, e.g., arising from other subdominant components of the foreground, may be treated with the suggested technique. The dependence of the other foreground components, e.g., free–free, dust, and spinning dust (see Lazarian & Finkbeiner 2003 for a review) on the frequency is weak and therefore the linear/quadratic expansion in frequencies in Equation (2)/Equation (9) should be also valid. If we observe, for example, 610 MHz, the 21 cm signal is believed to be decorrelated if Δν > 0.5 MHz (see, for example, Ghosh et al. 2011). If we use 12 channels, we have Δν/ν₀ = 5.5/610 ≈ 0.01. Therefore, any third or higher order effects should be very small. For synchrotron, since (ν/ν₀)⁻¹ = ∑^∞_{n = 0}(1 − ν/ν₀)ⁿ, the ratio of the third to the second term of the expansion is ∼O(ν/ν₀ − 1) ∼ 0.01. For free–free the frequency dependence is weaker and truncation error will be smaller.

In this paper we assume that the cross-correlation between the foreground and the 21 cm signals is negligible. This needs further verification. In general, we can write C_l = C^for_l + C^21 cm_l + C_l^cc, where

$$\begin{equation} C_l^{{\rm cc}}={\rm Re} \left[ \frac{2}{N_{{\rm sum}}} \sum _{ l/2\pi } T^{{\rm for}}(u,v)T^{21\,{\rm cm}}(u,v)^* \right] \end{equation} \tag{ 13 }$$

is the cross-correlation term.

When the 21 cm signal is γ times smaller than the foreground one (i.e., T^21 cm ∼ T^for/γ), the cross-correlation term C^cc_l can be larger than C^21 cm_l if the number of modes used for the summation in Equation (13) is less than ∼γ².

However, if we use Equation (12) or (8), we may circumvent this difficulty. Let us consider the more general case of Equation (12). The cross-correlation term that will appear additionally in the right-hand side of Equation (12) is

$$\begin{eqnarray} &&C_l^{{\rm cc}} = {\rm Re} \left[ \frac{2}{N_{{\rm sum}}} \sum _{ l/2\pi } (\mathcal {A}-\mathcal {B}) \right], \end{eqnarray} \tag{ 14 }$$

$$\begin{eqnarray} &&\mathcal {A}=\frac{1}{N}\sum _{i=1}^{N/2} \left[ (2 f_0+2b\Delta \nu ^2)\big(T_{\nu _i}^{21\,{\rm cm}} +T_{\nu _{N-i+1}}^{21\,{\rm cm}}\big)^* \right], \end{eqnarray} \tag{ 15 }$$

$$\begin{eqnarray} &&\mathcal {B}= \frac{2}{N^2} \left(\sum _{i=1}^{N}(f_0+b \Delta \nu ^2) \right) \left(\sum _{i=1}^{N} T_{\nu _i}^{21\,{\rm cm}*} \right), \end{eqnarray} \tag{ 16 }$$

$$\begin{eqnarray} &&| \mathcal {A-B} | \lesssim 2|b| \Delta \nu _{{\rm max}}^2 \big|T_{\nu _i}^{21\,{\rm cm}}\big| /\sqrt{N}, \end{eqnarray} \tag{ 17 }$$

where the variables are defined similarly (i.e., Δν = ν_i − ν₀, N = even, ...), "$\mathcal {A}$" and "$\mathcal {B}$" are from the first and the second terms on the right-hand side of Equation (12), respectively, Δν_max = ν_N − ν₀, and we dropped "(u, v)" for simplicity. Note that the terms containing f₀ cancel out.

Therefore, if

$$\begin{equation} 2b\Delta \nu _{{\rm max}}^2\big/\sqrt{NN_{{\rm sum}}} < \big|T_{\nu _i}^{21\,{\rm cm}}\big|, \end{equation} \tag{ 18 }$$

we can neglect the cross-correlation term. If the cross-correlation term is indeed larger than the 21 cm spectrum C^21 cm_l for any reason, our technique can reveal the cross-correlation spectrum. Even in this case, it may not be difficult to infer the shape of C^21 cm_l from the shape of C^cc_l.

Our technique is different from standard polynomial-fitting approaches (e.g., McQuinn et al. 2006; Liu et al. 2009). First, although we assume a polynomial dependence of the foregrounds on frequency, we do not need the actual fitting stage. Therefore, our technique can be used as a synergistic technique to previous standard suggestions. Second, since C_l∝∑|T²|, we expect that the high-order frequency dependence of T may cancel out during the summation. Therefore, in the power spectrum we may have a smoother signal which is much easier to constrain. Third, our technique can give useful insight into removing noise. Suppose that noise signals are larger than or comparable to the 21 cm signals and that the 21 cm signals at ν and ν' are correlated if |ν − ν'| ≲ Δν = ν_{i + 1} − ν_i. In this case our technique will give an average spectrum of the noise and the 21 cm signals, ∼C^noise_l + C_l^21 cm. Now, we add M − 1 new channels between each ν_i and ν_{i + 1} (i = 1, ..., N), so that there are a total of NM channels. Here we assume that the total bandwidth stays roughly the same and the bandwidth of each channel is now ∼M times narrower than before. Then we apply our technique again. Note that the 21 cm signals are correlated in adjacent channels, while the noise signals are not. The newly obtained average spectrum is now proportional to ${\sim} C_l^{{\rm noise}}+C_l^{21\,{\rm cm}}/\sqrt{M}$. Therefore, by comparing the two results, we may extract the spectrum of the 21 cm signals. For example, if the noise and the 21 cm signals have different dependence on l (or k), we can apply the technique in Cho & Lazarian (2010).

We have described a new technique to remove the foreground from redshifted 21 cm observations. We have assumed that the foreground signals change slowly as the frequency changes: we have assumed up to a quadratic dependence of $T_{\nu _i}^{{\rm for}}(u,v)$ on frequency. We obtain the foreground spectrum by simply adding all the observed $T_{\nu _i}(u,v)$. If the observations are done in real space, we add all the $T_{\nu _i}(\theta _x,\theta _y)$ in real space first and then we perform the Fourier transform. When we add the $T_{\nu _i}(u,v)$ of all channels, the foreground spectrum goes up by a factor of ∼N² because they are highly correlated. However, the 21 cm spectrum goes up by a factor of ∼N because the signals in different channels contribute randomly. This way, we can obtain an accurate shape of the foreground power spectrum. Then, we obtain the 21 cm power spectrum by subtracting the foreground power spectrum obtained this way. We have described how to obtain the average 21 cm spectrum (Equations (8) and (12)). We have derived the condition for neglecting the cross-correlation term (Equation (18)).

J.C.'s work was supported by the Basic Science Research Program through the National Research Foundation of Korea (NRF) funded by the Ministry of Education, Science, and Technology (2011-0012081). A.L. acknowledges the support of the NSF Center for Magnetic Self-Organization, the NASA grant NNX11AD32G, and the NSF grant AST 0808118. P.T. is partly supported by NSF grant AAG 0908900.