ALL-SKY INTERFEROMETRY WITH SPHERICAL HARMONIC TRANSIT TELESCOPES

Mapping the universe with the 21 cm line of neutral hydrogen will revolutionize our view of the universe. It holds the promise of unraveling the mysteries of dark energy (Chang et al. 2008; Loeb & Wyithe 2008), unveiling the epoch of reionization (Furlanetto et al. 2006), and perhaps even extending our view of the cosmos out far enough to shine light on the primordial dark ages (Loeb & Zaldarriaga 2004). Rapidly probing large volumes of the universe requires new large wide-field telescopes along with powerful new digital processing hardware such as GMRT,⁴ LOFAR,⁵ MWA,⁶ MITEoR (Zheng et al. 2013), PAPER,⁷ BAOBAB,⁸ BAORadio,⁹ BINGO (Battye et al. 2013), CHIME,¹⁰ EMBRACE/EMMA (Kant et al. 2011), and Tianlai,¹¹ In recent years it has become increasingly clear that new methods of interpreting and analyzing the data from these revolutionary new instruments will be necessary to realize their scientific potential (Myers et al. 2003; Tegmark & Zaldarriaga 2009; Parsons & Backer 2009; Liu et al. 2010; Liu & Tegmark 2011; Parsons et al. 2012; Dillon et al. 2013).

We describe here the spherical harmonic transit telescope, a new paradigm for analyzing wide-field transit telescopes in the spherical harmonic domain, which is naturally suited to mapping the 21 cm universe. Any telescope with fixed pointing observes the sky transit through its field of view. The periodic rotation about the poles over the course of a sidereal day creates a linear correspondence between time, t, and azimuthal angle, ϕ. We obtain a simple mapping between the observed data and a linear combination of the spherical harmonic coefficients a_lm of the sky at fixed angular wavenumber m, mediated by the angular response of each element. In the following, we elaborate and make precise this basic idea in the context of wide-field interferometers, including the radial (frequency) direction.

This formalism diverges sharply from traditional characterizations of radio interferometry, which are better suited to observations with a narrow field of view, often assuming tracking of a particular source of interest, and exploit Fourier transform mapping between the sky and the uv-plane. What we describe here is an all-sky formalism for describing interferometry that naturally incorporates the observable modes on the sky—the spherical harmonics. Previous work has described interferometry on the full sky (Kim 2007; McEwen & Scaife 2008), however, in the domain of transit telescopes there are many computational advantages we can exploit.

The foremost challenge for any 21 cm mapping experiment is separating the cosmological signal from astrophysical contaminants which are 10³–10⁵ times larger (Furlanetto et al. 2006; Morales & Wyithe 2010). Conceptually this is simple—the primary foreground sources (diffuse synchrotron emission from the Galaxy and emission from extragalactic point sources) are smooth as a function of frequency, while the 21 cm signal decorrelates quickly as each frequency corresponds to a different radial slice of the universe. To remove foregrounds one just needs to model and remove the smooth frequency component from their observations. Unfortunately, in practice, the large dynamic range between the amplitude of the foregrounds and the 21 cm signal makes several real-world effects extremely problematic. While the properties of the cosmological 21 cm signal are thought to be well understood, the astrophysical foregrounds are poorly constrained at the small angular and frequency scales that will be probed by forthcoming 21 cm experiments—a successful technique should be robust to uncertainties in foreground modeling. Of course, these experiments will help characterize the properties of real-world foregrounds. More troublesome is the phenomenon of angular-frequency mode-mixing: in any real experiment the shape of the beam on the sky will vary with the observed frequency (Liu et al. 2009). This mode-mixing makes simple frequency-only foreground methods ineffective in practice. We show below that the spherical harmonic transit telescope formalism can naturally address the issues of model uncertainty and mode-mixing, and enable efficient and effective discrimination of the 21 cm signal from obscuring foregrounds.

Any foreground removal method aims to find a subset of the data within which there is significantly more 21 cm signal than astrophysical foregrounds. However, in the presence of mode-mixing, it is not obvious how to select a basis which separates the two components—what we would like is a method which can automatically generate this. Just such a technique exists in the form of the Karhunen-Loève (KL) transform. In this paper, we show how the m-mode formalism, described below, makes the use of the KL transform computationally feasible. The result is a remarkably effective and robust filter for rejecting bright foregrounds and we demonstrate its effectiveness using realistic simulations of the radio sky and a simple fiducial interferometer configuration.

In Section 2, we introduce the all-sky formalism that is the basis for this technique. In Section 3, we discuss the map-making process in the spherical harmonic transit telescope paradigm. In Section 4, we discuss how to best represent statistics of the cosmological 21 cm signal and foregrounds in the measurement basis, and in Section 5 we discuss how the KL transform can be used to detect faint signals in the presence of bright foregrounds. In Section 6, we quantify the information lost due to foreground removal using Fisher Analysis, and estimate errors on power spectra. We present our conclusions in Section 7. In Appendix A, we discuss the signal and foreground models we employ. In Appendix B, we describe how we create realistic simulations of radio emission.

In this section, we introduce the m-mode formalism, a new description of the measurement process for transit interferometers.

In radio interferometry, a visibility V_ij is the instantaneous correlation between two feeds F_i and F_j. We will assume that we can take a linear combination of the signal from a dual polarization antenna, with no cross-polarization or polarization leakage, such that we are sensitive only to the total intensity (Stokes I) part of the sky. The fully polarized extension to this work is also a tractable problem, we address this in a subsequent paper, J. R. Shaw et al. (in preparation; for a treatment for a non-transit telescope, see Kim 2007). At any instant, a visibility is given by

$$\begin{eqnarray} V_{ij} & = \left\langle F_i F_j^* \right\rangle \nonumber \\ & = \frac{1}{\Omega _{ij}}\int d^2\hat{\boldsymbol{n}}\, A_i(\hat{\boldsymbol{n}}) A_j^*(\hat{\boldsymbol{n}}) e^{2 \pi i \hat{\boldsymbol{n}}\cdot \boldsymbol{u}_{ij}} T(\hat{\boldsymbol{n}}), \end{eqnarray} \tag{ 1 }$$

where $\boldsymbol{u}_{ij} = (\boldsymbol{r}_i - \boldsymbol{r}_j) / \lambda$ is the spatial separation between the two feeds divided by the observed wavelength (that is the separation in the uv-plane), $\hat{\boldsymbol{n}}$ is the position on the celestial sphere, and $A_i(\hat{\boldsymbol{n}})$ gives the primary beam of feed i. In the above we have normalized our visibilities such that they are temperature-like, and we have defined them in terms of the brightness temperature T = λ²I/2k_b instead of the total intensity I. The quantity $\Omega _{ij} = \sqrt{\Omega _i \Omega _j}$ is the geometric mean of the individual beam solid angles

$$\begin{equation} \Omega _i = \int |A_i(\hat{\boldsymbol{n}}) |^2 d^2\hat{\boldsymbol{n}}, \end{equation} \tag{ 2 }$$

which also gives the effective antenna area A_effΩ = λ². This ensures that for a sky with uniform brightness temperature T, the auto-correlation of an antenna V_ii = T with our definition.

As the Earth turns, both the primary beams and the baseline separations rotate relative to the celestial sphere. This means that the measured visibilities change periodically with the sidereal day. We take this into account by explicitly including the dependence on the azimuthal angle ϕ and by averaging over each sidereal day.

The measured visibilities are also corrupted by instrumental noise for which we add a noise term n_ij(ϕ). We assume the noise is stationary such that its statistics are independent of ϕ. Rewriting Equation (1) in terms of a transfer function B_ij leaves the measured visibility as

$$\begin{equation} V_{ij}(\phi) = \int \! d^2\hat{\boldsymbol{n}}\, B_{ij}(\hat{\boldsymbol{n}}; \phi) T(\hat{\boldsymbol{n}}) + n_{ij}(\phi), \end{equation} \tag{ 3 }$$

where the transfer function is

$$\begin{equation} B_{ij}(\hat{\boldsymbol{n}}; \phi) = \frac{1}{\Omega _{ij}}A_i(\hat{\boldsymbol{n}}; \phi) A_j^*(\hat{\boldsymbol{n}}; \phi) e^{2 \pi i \hat{\boldsymbol{n}}\cdot \boldsymbol{u}_{ij}(\phi)} \;. \end{equation} \tag{ 4 }$$

Taking advantage of the periodicity in ϕ, we Fourier transform the system

$$\begin{eqnarray} V^{ij}_m &= \int \frac{d\phi }{2\pi } V_{ij}(\phi) e^{-i m \phi } \end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray} & = \sum _{l m^{\prime }}\int \frac{d\phi }{2\pi } B^{ij}_{l m^{\prime }}(\phi) a_{l m^{\prime }} e^{-i m \phi }+ n^{ij}_m, \end{eqnarray} \tag{ 6 }$$

where to proceed to the second line we have inserted the spherical harmonic expansions of both the sky and the beam transfer function

$$\begin{eqnarray} T(\hat{\boldsymbol{n}}) & = \sum _{lm} a_{lm} Y_{lm}(\hat{\boldsymbol{n}}) \;, \end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray} B_{ij}(\hat{\boldsymbol{n}}; \phi) & = \sum _{lm} B^{ij}_{lm}(\phi) Y_{lm}^*(\hat{\boldsymbol{n}}) \;. \end{eqnarray} \tag{ 8 }$$

Note that we have defined $B^{ij}_{lm}$ relative to the conjugate spherical harmonic in order to simplify later notation. As the ϕ dependence simply rotates the functions about the Earth's polar axis, the transfer function at any ϕ is trivially $B^{ij}_{lm}(\phi) = B^{{ij}}_{lm}(\phi \!=\!0) e^{i m \phi }$. Combined with the exponential factor in the integral, this simply generates the Kronecker delta $\delta _{mm^{\prime }}$, and we find

$$\begin{equation} V^{ij}_m = \sum _{l} B^{ij}_{l m} a_{l m}+ n^{ij}_m \;. \end{equation} \tag{ 9 }$$

This gives a simple description of how the observed sky maps into the measured data given a telescope design (which is contained in the beam transfer matrices $B^{ij}_{lm}$). This transformation does not mix m-modes on the sky, and can therefore be performed on an m-by-m basis—for any particular m and frequency ν, the measured visibilities are simply a projection of the l-modes on the sky for the measured m. As the optical system is of a finite size, this limits both the l and m to which the telescope is sensitive, ensuring we only need to consider a finite number of degrees of freedom, both measured (V_m) and on the sky (a_lm).

In fact, while the positive and negative m-modes may be independent measurements, they are still observations of the same sky—by transforming the conjugate $V_{-m}^*$ and using that $a_{lm} = a^*_{l,-m}$ for a real field we see that

$$\begin{equation} V^{ij *}_{-m} = \sum _l (-1)^m B^{ij *}_{l, -m} a_{lm} + n^{ij *}_{-m} \;. \end{equation} \tag{ 10 }$$

In light of this, we will change our notation such that we are considering only the actual degrees of freedom on the sky. Let us separate out the positive and negative m parts by defining

$$\begin{eqnarray} B_{l m}^{ij, +} = B^{ij}_{l m} &\quad n^{ij, +}_{m} & = n^{ij}_{m} \end{eqnarray} \tag{ 11 }$$

$$\begin{eqnarray} B_{l m}^{ij, -} = (-1)^m B^{ij *}_{l, -m} &\quad n^{ij, -}_m & = n^{ij *}_{-m} \end{eqnarray} \tag{ 12 }$$

which is valid for m ⩾ 0. Additionally, to prevent double counting the m = 0 measurement, we need to set $B_{l 0}^{ij -} = n_0^{ij -} = 0$. This gives a modified version of Equation (14)

$$\begin{equation} V_{m}^{ij,\pm } = \sum _l B_{l m}^{ij, \pm } a_{lm} + n_{m}^{ij, \pm } \;. \end{equation} \tag{ 13 }$$

For brevity of notation, we will introduce a label α which indexes both the positive and negative m parts of all of the included feed pairs ij, such that any particular α specifies exactly the values of ij, ± (exactly how α is packed is unimportant). This gives

$$\begin{equation} V^\alpha _m = \sum _{l} B^\alpha _{l m} a_{l m}+ n^\alpha _m. \end{equation} \tag{ 14 }$$

The beam transfer matrices above can be written in an explicit matrix notation

$$\begin{equation} \left({{\bf B}}_m\right)_{(\alpha \nu) (l \nu ^{\prime })} = B^{\alpha, \nu }_{lm} \delta _{\nu \nu ^{\prime }}, \end{equation} \tag{ 15 }$$

where the row index labels all combinations of baseline (α) and frequency (ν), whereas the column index is over all multipoles (l) and frequencies (ν'). Similarly we can define vectors for the visibilities and harmonic coefficients

$$\begin{equation} \left(\boldsymbol{v}_m\right)_{(\alpha \nu)} = V^{\alpha, \nu }_m \, \quad \left(\boldsymbol{a}_m\right)_{(l \nu)} = a_{lm}^\nu. \end{equation} \tag{ 16 }$$

From here onward we will drop the subscript m denoting the spherical harmonic order, all the equations below are valid for any m. This allows us to rewrite Equation (14) as

$$\begin{equation} \boldsymbol{v}= {{\bf B}}\, \boldsymbol{a}+ \boldsymbol{n}. \end{equation} \tag{ 17 }$$

This simple linear description of the measurement process of a transit telescope is extremely powerful. By reducing it down to a linear mapping between a finite number of degrees of freedom, it allows us to apply the standard tools of signal processing. In the subsequent sections we apply it to solve two challenging problems in 21 cm radio astronomy.

In astronomy, being able to transform our measured signal into an accurate map of the sky is essential. Whilst in this paper we explicitly avoid this process for our analysis, preferring to carry it out directly in the data space, maps are still needed for visualization and cross-checking. Map-making with interferometric data is generally a complicated process performed by algorithms such as CLEAN (Högbom 1974) and its derivatives. This is especially true with wide fields of view where mosaicking and w-projection are generally required. However, the m-mode formalism makes the map-making process on the full sky conceptually simple.

First, we assume that the instrumental noise n follows a complex Gaussian distribution with covariance ${{\bf N}}= \langle \boldsymbol{n}\boldsymbol{n}^\dagger \rangle$, and the different frequency channels are independent.¹² For stationary noise, the m-modes are uncorrelated and the likelihood function of the observed sky for a single m and frequency ν is

$$\begin{equation} p(\boldsymbol{v}\vert \boldsymbol{a}) = \frac{1}{\left|\pi {{\bf N}}\right|} \exp {(- \left(\boldsymbol{v}- {{\bf B}}\boldsymbol{a}\right)^\dagger {{\bf N}}^{-1} \left(\boldsymbol{v}- {{\bf B}}\boldsymbol{a}\right))}, \end{equation} \tag{ 18 }$$

where the vector a contains all of the harmonic coefficients for the given m.

To estimate the sky corresponding to a given set of visibilities, we will look for a maximum likelihood solution $d p / d \boldsymbol{a}= 0$. In particular, we want to find the value of a that minimizes

$$\begin{equation} \chi ^2 = |{{\bf N}}^{-\frac{1}{2}} \boldsymbol{v}- ({{\bf N}}^{-\frac{1}{2}} {{\bf B}})\boldsymbol{a}|^2. \end{equation} \tag{ 19 }$$

The matrix N^−1/2 represents any factorization such that (N^−1/2)†N^−1/2 = N⁻¹. Provided N contains no noiseless modes, it is positive-definite and so this factorization must exist. The maximum likelihood solution is given by the Moore–Penrose pseudo-inverse ¹³

$$\begin{equation} \hat{\boldsymbol{a}} = ({{\bf N}}^{-\frac{1}{2}} {{\bf B}})^+ {{\bf N}}^{-\frac{1}{2}} \boldsymbol{v}\;, \end{equation} \tag{ 20 }$$

where the superscript + denotes the pseudo-inverse.

Depending on the number of baselines measured and the maximum l we are sensitive to, the problem may either be over- or under-constrained. In either regime, the Moore–Penrose pseudo-inverse gives a solution; in the former case this reduces to the standard map-making equation $\hat{\boldsymbol{a}} = ({{\bf B}}^\dagger {{\bf N}}^{-1} {{\bf B}})^{-1} {{\bf B}}^\dagger {{\bf N}}^{-1} \boldsymbol{v}$, and in the latter case selects the solution which also minimizes $\left|\hat{\boldsymbol{a}}\right|^2$, effectively setting the unconstrained degrees of freedom to zero.

As both distinct frequencies and m-modes are independent, map-making for a set of full sky observations is a case of collating the estimates for each individual ν and m.

For intensity mapping experiments, our data has three components: the 21 cm signal which we are ultimately trying to extract, the foregrounds, and instrumental noise. Understanding the two-point statistics of the data is of paramount importance to our analysis—not only do the correlations of the signal encode most of the cosmological information that we are interested in (see Appendix A), but to efficiently extract this we require knowledge of the two-point statistics of all three components. Here we write the linear relationship between these two-point statistics of the data, and how they are related to the underlying physical correlations.

The statistics of instrumental noise live in the visibility space, the basis of our measurements. However, the other components are naturally represented on the sky, and must be projected into this space using Equation (14). The lowest non-zero moment of the visibilities is their covariance

$$\begin{eqnarray} C_{(\alpha \nu m); (\alpha ^{\prime } \nu ^{\prime } m^{\prime })} &=& \left\langle V^m_{\alpha \nu } V^{m *}_{\alpha ^{\prime } \nu ^{\prime }} \right\rangle \nonumber \\ &=& \sum _{l l^{\prime }} B^{\alpha \nu }_{l m} \left\langle a_{l m \nu }^{*} a_{l^{\prime } m^{\prime } \nu ^{\prime }} \right\rangle B^{\alpha ^{\prime } \nu ^{\prime } *}_{l^{\prime } m^{\prime }} + \left\langle n_{(\alpha \nu m)} n^*_{(\alpha ^{\prime } \nu ^{\prime } m^{\prime })} \right\rangle \;. \nonumber\\ \end{eqnarray} \tag{ 21 }$$

This is the covariance between all measured degrees of freedom: baselines, frequencies, and m-modes. For the experiments listed in Section 1, we expect ≳ 10³, ∼10², and 10³, respectively. This gives matrices of dimension ≳ 10⁸, too large to be tackled with current technology, both in terms of computation and storage.

Instead, let us make an approximation that will dramatically reduce this complexity. If we think of the sky as a statistically isotropic random field, its two-point statistics become dramatically simpler

$$\begin{equation} \langle a_{l m \nu ^{\prime }} a_{l^{\prime } m^{\prime } \nu ^{\prime }}^{*} \rangle = C_l(\nu, \nu ^{\prime }) \delta _{l l^{\prime }} \delta _{m m^{\prime }}, \end{equation} \tag{ 22 }$$

and importantly, they are automatically uncorrelated in the m index. This means that the full signal covariance Equation (21) is block diagonal and thus allows us to calculate all statistics on an m-by-m basis. For a specific m-mode

$$\begin{equation} C_{(\alpha \nu); (\alpha ^{\prime } \nu ^{\prime })} = \sum _{l} B^{\alpha \nu }_{l} B^{\alpha ^{\prime } \nu ^{\prime } *}_{l^{\prime }} C_l(\nu, \nu ^{\prime }) + N_{(\alpha \nu); (\alpha ^{\prime } \nu ^{\prime })} \;, \end{equation} \tag{ 23 }$$

where we have dropped all the m-indices, and N is the power spectrum on the instrumental noise. As the number of m-modes we are sensitive to is usually ≳ 10³, assuming statistical isotropy saves at least a factor of a million in computation and a thousand in storage. This ability to efficiently perform calculations incorporating the full statistics opens up new avenues for the data analysis of transit instruments. Synchrotron emission from our Galaxy clearly violates this assumption of statistical isotropy though, as we will demonstrate, this does not appear to affect our analysis and in particular our ability to clean foregrounds.

In matrix notation

$$\begin{equation} {{\bf C}}= {{\bf B}}{{\bf C}}_{\rm sky} {{\bf B}}^\dagger + {{\bf N}}\;. \end{equation} \tag{ 24 }$$

where we will split C_sky into independent 21 cm signal and foreground parts C_sky = C₂₁ + C_f.

The statistical models used for each component are chosen to be appropriate for the frequency ranges of interest. In the fiducial example that follows this is 400–600 MHz, corresponding to z ∼ 1–2 for the cosmological signal. These are described in Appendix A.

To clean our data we simply aim to find a subset within which there is significantly more 21 cm signal than astrophysical foregrounds. However, in the presence of mode-mixing there is no immediately apparent representation in which to perform this. This basis can be found using the KL transform (often called the signal-to-noise eigendecomposition), which has a long history in cosmology (e.g., Bond 1995; Vogeley & Szalay 1996; Tegmark et al. 1997) and has been used for the analogous problem of E/B mode separation for polarization of the cosmic microwave background (Lewis et al. 2002; Bunn et al. 2003). This transform simultaneously diagonalizes both the signal and foreground covariance matrices, generating a set of modes with no foreground or signal correlations. This makes comparing the amount of signal and foreground power in each mode trivial.

To reduce the risk of foreground uncertainties biasing our analysis, we will prioritize the removal of foreground contaminated modes at the expense of cosmological signal. In contrast, the instrumental noise is well understood, and we should be able to dig deeper into this contaminant to extract useful cosmological information with little risk. In practice, this means we will start with a filter which aggressively removes foregrounds only; subsequently we will add the instrumental noise back in, which will allow us to compress the data by completely removing noise-dominated modes, while retaining those with a small fraction of signal.

This requires models for the statistics of both the signal and the foregrounds. The signal is modeled as a simple Gaussian random field for the 21 cm emission, whereas the foreground model includes both the synchrotron emission from our Galaxy, and the contribution from a background of extragalactic point sources. The details of both are discussed in Appendix A.

The KL transform seeks to find a linear transformation of the data $\boldsymbol{v}^{\prime } = {{\bf P}}\boldsymbol{v}$ such that the 21 cm signal S = BC₂₁B† and foreground F = BC_fB† covariance matrices are jointly diagonalized. That is

$$\begin{equation} {{\bf S}}\rightarrow {{\bf S}}^{\prime } = {{\bf P}}{{\bf S}}{{\bf P}}^\dagger = {{\bf \Lambda }}\;, \end{equation} \tag{ 25 }$$

and

$$\begin{equation} {{\bf F}}\rightarrow {{\bf F}}^{\prime } = {{\bf P}}{{\bf F}}{{\bf P}}^\dagger = {{\bf I}}\;, \end{equation} \tag{ 26 }$$

where Λ is a diagonal matrix, and I is the identity. In this diagonal basis we can simply compare the amount of power expected in each mode by the ratio of the diagonal elements (this is given by the corresponding entries of Λ), and identify the regions of the space with low foreground contamination (large entries in Λ).

This transformation can be found by solving the generalized eigenvalue problem ${{\bf S}}\boldsymbol{x}= \lambda {{\bf F}}\boldsymbol{x}$ (for our purposes the routine ZHEGVD from LAPACK is sufficient). This gives a set of eigenvectors x, and corresponding eigenvalues λ. Writing the eigenvectors in a matrix P, row-wise, gives the transformation matrix to diagonalize the covariances. The eigenvalues λ corresponding to each eigenvector give the diagonal matrix Λ.

To isolate the 21 cm signal, we want to select modes with eigenvalue (signal-to-foreground power; S/F) greater than some threshold (see Figure 1). To project into this basis we define the matrix P_s which contains only the rows from P corresponding to eigenvalues greater than the threshold s.

Zoom In Zoom Out Reset image size

For most of the analysis we can work directly in the eigenbasis. However, for visualizing our results, we want to be able to transform back to the sky (by way of the measured visibilities). To project back into the higher dimensional space we simply generate the full inverse P⁻¹ and remove columns corresponding to the rejected modes (we denote this matrix $\bar{{{\bf P}}}_{s}$). This is equivalent to projecting into the full eigenbasis, zeroing the foreground contaminated modes, and then using the full-inverse P⁻¹.

For further analysis, we must include all noise terms, both foregrounds and instrumental. Writing the total noise contribution as N_all = F + N, the matrix in the truncated basis is

$$\begin{eqnarray} {{\bf N}}^{\rm all} \rightarrow {{\bf N}}^{\rm all}_s & =& {{\bf P}}_s \left({{\bf F}}+ {{\bf N}}\right){{\bf P}}_s^\dagger\\ \ms & =& {{\bf I}}+ {{\bf P}}_s {{\bf N}}{{\bf P}}_s^\dagger. \end{eqnarray} \tag{ 27 }$$

As the transformed instrumental noise matrix will not remain diagonal, this gives a correlated component between all our modes. However, as it is useful if our modes are uncorrelated, we make a further KL transformation on the foreground removed signal ${{\bf S}}_s = { {\boldsymbol \Lambda }}_s$, and total noise ${{\bf N}}^{\rm all}_s$ covariance matrices. For computational and storage efficiency we apply a further cut-off to include only modes with a signal-to-total-noise ratio greater than the cut-off value t. We denote this as projection matrix $\tilde{{{\bf Q}}}$. For notational convenience we will write the total transformation in terms of a single matrix R = Q_tP_s, having chosen suitable values for the two cut-offs s and t. As above, we will define an inverse $\bar{{{\bf R}}} = \bar{{{\bf P}}}_s \bar{{{\bf Q}}}_t$, which remains orthogonal to the removed space. Quantities in this final basis we denote with tildes, for example a visibility mapped into this basis is $\tilde{\boldsymbol{v}} = {{\bf R}}\boldsymbol{v}$, and a covariance is $\tilde{{{\bf C}}}= {{\bf R}}{{\bf C}}{{\bf R}}^\dagger$.

5.1. Cylinder Example

While this method works with any transit telescope, to illustrate the foreground removal process we will simulate a cylinder telescope, such as CHIME or the Pittsburgh Cylinder Telescope (Bandura 2011). These are transit interferometers composed of multiple parabolic cylinders where each only focuses in the East–West direction. This gives a long and thin primary beam on the sky, extending nearly from horizon to horizon in the North–South direction but which is only around 1°wide East–West.

Feeds are spaced along the axis of each cylinder—when correlated these provide resolution in the N–S direction. Correlations between cylinders enhance the E–W resolution. The telescope operates as a transit telescope such that the entire visible sky is observed once per sidereal day.

For a cylinder uniformly illuminated by a particular feed, near the axis the beam pattern is a sinc function in the E–W direction, and uniform in the N–S direction (Wilson et al. 2009, chapter 6). To extend this off-axis, we modulate by projected area of the telescope giving

$$\begin{equation} A^2(\hat{\boldsymbol{n}}) = \rm {sinc}^2{\left(\pi \, \hat{\boldsymbol{n}}\cdot \hat{\boldsymbol{u}}\frac{W}{\lambda } \right)} \, \Theta \left(\hat{\boldsymbol{n}}\cdot \hat{\boldsymbol{z}}\right)\, \hat{\boldsymbol{n}}\cdot \hat{\boldsymbol{z}} \end{equation} \tag{ 29 }$$

where W is the cylinder width and $\hat{\boldsymbol{z}}$ is a unit vector pointing to the zenith and $\hat{\boldsymbol{u}}$ is a unit vector pointing East in the ground-plane. The step function Θ masks out the regions where the sky is below the horizon, and the final factor $\hat{\boldsymbol{n}}\cdot \hat{\boldsymbol{z}}$ accounts for the projected area of the telescope.

We model a two cylinder telescope observing the sky with 64 frequency channels from 400–600 MHz. Each cylinder is 15 m wide and has 60 feeds regularly spaced by 0.25 m (with the feeds lining up E–W between cylinders). The telescope is located at a latitude of 45°. These specifications correspond to a slightly smaller half-bandwidth version of the CHIME pathfinder telescope being constructed at the Dominion Radio Astrophysical Observatory.

The noise covariance is diagonal for all m, frequencies, and baselines. For small m-modes with m ≪ 1/(2πΔϕ) (where Δϕ is the angular integration time), the noise variance is

$$\begin{equation} N^{ij}_m = \frac{T_{{\rm sys},i}(\nu) T_{{\rm sys},j}(\nu)}{4 \pi N_{\rm day} t_{\rm sid} \Delta \nu } \;, \end{equation} \tag{ 30 }$$

where T_{sys, i} is the system temperature of a single polarization of feed i, N_day is the number of sidereal days observed, t_sid is the length of a sidereal day, and Δν is the width of the frequency channel. As we combine the two polarizations into a single unpolarized signal, this reduces the noise power spectrum by a factor of two. For this example T_sys = 50 K, and we assume two full years of observation (that is, 730 complete sidereal days).

In Figure 1 we show the spectrum of S/F eigenvalues for the telescope. The KL mode distribution of S/F has an extremely rapidly rising spectrum so that the information retained (approximately the number of modes) is quite insensitive to the cut threshold s for values between 10⁻² and 10².

To demonstrate the foreground removal process, we simulate time-streams from separate realizations of the signal and foregrounds using Equation (17), and project them through the filtering process to make maps. The visibilities are filtered using $\boldsymbol{v}_{\rm clean} = \bar{{{\bf R}}} {{\bf R}}\boldsymbol{v}$, and then are turned into a three-dimensional (3D) map using Equation (20). In Figure 2 we show the original simulation, the map made from the unfiltered visibilities, and the map made from the foreground filtered visibilities. The simulated signal and foreground maps are described in Appendix B. Note that the foreground maps are not simply realizations of the model used to generate the foreground filter—unlike the input model they are both non-Gaussian and anisotropic. Figure 2 clearly illustrates how the foreground amplitude is dramatically reduced by the process, whilst the signal retains its overall character. Though the foreground residuals are clearly highest in the Galactic center, even these are significantly lower than the filtered signal.

Zoom In Zoom Out Reset image size

In the previous section, we have demonstrated that the KL transform gives an effective method for removing foregrounds. Though a visual inspection of Figure 2 suggests that the 21 cm signal is largely untouched, we would like to be able to quantify how much useful information remains. In this section, we will make use of the Fisher Information matrix, a measure of the information contributed to our knowledge of a set of parameters for a given observation, by telling us how much the covariance of those parameters is improved (see Dodelson 2003, chapter 11 for an overview). We will use this to forecast power spectrum errors, for the same telescope, with and without foreground removal.

After projection into the reduced eigenbasis, let us assume that the remaining modes follow a complex Gaussian distribution with zero mean. This assumption should be reasonable provided we have successfully removed the modes containing any significant foreground contribution. In this case the contribution of the m^th mode to the ab-element of the Fisher Information matrix for a set of parameters p_a is

$$\begin{equation} F_{ab}^{\scriptscriptstyle \left(m\right)} = \rm {tr}{(\tilde{{{\bf C}}}_a \tilde{{{\bf C}}}^{-1} \tilde{{{\bf C}}}_b \tilde{{{\bf C}}}^{-1})} \;. \end{equation} \tag{ 31 }$$

where $\tilde{{{\bf C}}}_a = \partial \tilde{{{\bf C}}}/ \partial p_a$, and the common factor of 1/2 does not appear as we are dealing with a complex Gaussian distribution. Though in the constructed eigenbasis $\tilde{{{\bf C}}}= \tilde{{{\bf \Lambda }}}+ {{\bf I}}$ is diagonal, $\tilde{{{\bf C}}}_a$ can have off-diagonal elements. Again this process is performed on a per-m basis. As there is no coupling between them, the total Fisher information is simply the sum over all m-modes

$$\begin{equation} F_{ab} = \sum _m F_{ab}^{\scriptscriptstyle \left(m\right)} \;. \end{equation} \tag{ 32 }$$

For a set of parameters p_a that we are trying to determine, the inverse of the Fisher matrix is the lowest order approximation to their covariance.

In this work, we will focus on forecasting the errors on the shape of the matter power spectrum P(k) whilst keeping all other cosmological parameters fixed. Such forecasting has been performed using the uv-plane in Seo et al. (2010) and Ansari et al. (2012).

We parameterize the power spectrum in terms of a linear summation of different basis functions

$$\begin{equation} P(k) = \sum _a p_a P_a(k) \;. \end{equation} \tag{ 33 }$$

In Appendix A, we describe how to project this quantity into the angular power spectrum of 21 cm fluctuations that we use to calculate the visibility correlations. For simplicity, each of our bands is simply equal to the input power spectrum within a fixed k-band, and zero outside, such that the fiducial model is p_a = 1.

For the band-powers p_a that we are trying to estimate, the matrices $\tilde{{{\bf C}}}_a$ are simply the projection of the basis functions $P_a(\boldsymbol{k})$ into the eigenbasis. Starting from the angular power spectra C_{a; l}(ν, ν'), corresponding to each of the basis functions $P_a(\boldsymbol{k})$ (using Equation (A2))

$$\begin{equation} \tilde{{{\bf C}}}_{a} = {{\bf R}}{{\bf B}}{{\bf C}}_{21,a} {{\bf B}}^\dagger {{\bf R}}^\dagger \;. \end{equation} \tag{ 34 }$$

In practice, explicitly calculating the $\tilde{{{\bf C}}}_a$ in this way is computationally expensive, we instead use a Monte Carlo technique. We can form the estimator $\hat{q}_a = \tilde{\boldsymbol{x}}^\dagger \tilde{{{\bf C}}}^{-1} \tilde{{{\bf C}}}_a \tilde{{{\bf C}}}^{-1} \tilde{\boldsymbol{x}}$, which has the property that its covariance $\langle \hat{q}_a \hat{q}_b \rangle - \left\langle \hat{q}_a \right\rangle \left\langle \hat{q}_b \right\rangle = F_{ab}$ (Padmanabhan et al. 2003). This means we can estimate the F_ab by averaging over realizations of $\tilde{\boldsymbol{x}}$. For details, see Dillon et al. (2013).

In Figure 3 we plot the power spectrum errors for two cases: in the presence of foregrounds that have been cleaned using our method and without foregrounds at all. In the case without foregrounds, F = 0 and we only perform the final KL transform to diagonalize the signal and instrumental noise. For the foregrounds we have cleaned modes with S/F < 10 and additionally have removed modes with a small ratio of signal-to-total-power. This corresponds to setting s = 10 and t = 0.01. This is a clear demonstration of the effectiveness of the technique—it reduces our sensitivity on large scales as we would expect (as the removed foreground are smooth on large scales), while only slightly reducing our ability to constrain the small scale power spectrum.

Zoom In Zoom Out Reset image size

In this paper, we have introduced a powerful formalism for describing the measurement process of transit telescopes (either interferometric or otherwise). It is a natural formalism to describe interferometry on the full sky—sidestepping the standard complications that arise when dealing with wide-field interferometric data such as mosaicking and w-projection. A spherical harmonic transit telescope allows for compact and computationally efficient representations of the data and its statistics, which enable new ways of approaching important problems like map-making and foreground removal.

Using the m-mode formalism and approximating the foregrounds as statistically isotropic, allows the powerful KL transformation to be used, automatically finding the basis in which the astrophysical foregrounds and 21 cm signal are maximally separated. The KL approach would be computationally impossible otherwise and is a key advantage of the m-mode formalism. Using this technique we can take the full 3D data set into account and overcome the mode-mixing problem. The filters we construct are highly effective and robust, a fact we have demonstrated by propagating through realistically simulated 21 cm and foreground timestreams. In our fiducial example, shown in Figure 2, peak-to-peak foreground amplitude was reduce by a factor of ∼2 × 10⁷ leaving the peak-to-peak amplitude of the 21 cm signal around 20 times brighter that the foreground residuals.

We have also used this formalism to produce realistic forecasts for the power spectrum constraints from a fiducial 21 cm cylinder interferometer. We have demonstrated that foreground cleaning does not significantly degrade 21 cm power spectrum estimates on baryon acoustic oscillation (BAO) scales and below compared to a hypothetical foreground-free measurement. We anticipate that the spherical harmonic transit telescope formalism will be a powerful tool that can be applied to inform experimental design and test the interplay between real-world systematics and design constraints on 21st century 21 cm science. We will explore this further in J. R. Shaw et al. (2013, in preparation).

We model the 21 cm signal and foregrounds as isotropic fields described by an angular power spectrum C_l(ν, ν'). We base our models on Santos et al. (2005), although we will only include the Galactic synchrotron and extragalactic point source contributions. Both these contributions are assumed to take the form

$$\begin{eqnarray} &&C_l(\nu, \nu ^{\prime }) = A \left(\frac{l}{100}\right)^{-\alpha } \! \left(\frac{\nu \nu ^{\prime }}{\nu _0^2} \right)^{-\beta } \! e^{ - \frac{1}{2 \xi ^2_l} \ln ^2{(\nu / \nu ^{\prime })}} \,. \end{eqnarray} \tag{ A1 }$$

As the models are calibrated for observations of the reionization epoch, we need to transform them into the higher frequencies we are concerned with. We list the parameters for both these models in Table 1.

For the point source model, which is based on the results of Di Matteo et al. (2002), we change the pivot frequency ν₀ from 150 MHz to 408 MHz and also rescale the amplitude in order to raise the maximum flux of unsubtracted sources from 0.1 mJy to 0.1 Jy.

The galactic synchrotron model we use is not only calibrated for low frequencies, but also high Galactic latitudes. As we will measure large fractions of the sky, we take this into account by changing the A and angular power law index β to be consistent with the angular power spectrum of the 408 MHz Haslam map for Galactic latitudes |b| > 5° from La Porta et al. (2008).

We model the 21 cm brightness temperature as being a biased tracer of the underlying matter fluctuations. These fluctuations are naturally characterized by the angular power spectrum (Lewis & Challinor 2007; Datta et al. 2007). However exact calculation of this quantity requires double-integration over highly oscillatory functions, instead we use the flat-sky approximation from Datta et al. (2007)

$$\begin{equation} C_l(z, z^{\prime }) = \frac{1}{\pi \chi \chi ^{\prime }} \int _0^\infty \!\! dk_\parallel \cos {\left(k_\parallel \Delta \chi \right)} P_{T_b}(\boldsymbol{k}; z, z^{\prime }), \end{equation} \tag{ A2 }$$

where χ and χ' are the comoving distances to redshift z and z'. Their difference is denoted by Δχ = χ − χ'. The vector k has the components k_∥ and $l / \bar{\chi }$ in the directions parallel and perpendicular to the line of sight ($\bar{\chi }$ is the mean of χ and χ'). This approximation is accurate to the 1% level for l > 10 (Datta et al. 2007).

We model the 21 cm brightness temperature power spectrum $P_{T_b}$ as

$$\begin{equation} P_{T_b}(\boldsymbol{k}; z, z^{\prime }) = \bar{T}_b(z) \bar{T}_b(z^{\prime }) (b + f \mu ^2)^2 P_m(k; z, z^{\prime }), \end{equation} \tag{ A3 }$$

where P_m(k; z, z') = P(k)D₊(z)D₊(z') is the real-space matter power spectrum, D₊(z) is the growth factor normalized such that D₊(0) = 1, b is the bias, and the growth rate f = dln D₊/dln a, the logarithmic derivative of the growth factor D₊. We assume that the bias b = 1 at all redshifts. The mean brightness temperature is assumed to take the form

$$\begin{eqnarray} \bar{T}_b(z) &=& 0.3 \left(\frac{\Omega _{{\rm H}\,\scriptsize {\rm I}}}{10^{-3}}\right)\nonumber \\ &&\times \left(\frac{\Omega _m + (1+z)^{-3}\Omega _\Lambda }{0.29} \right)^{-1/2} \left(\frac{1+z}{2.5}\right)^{1/2} \,{{\rm mK}} \quad\quad \end{eqnarray} \tag{ A4 }$$

given in Chang et al. (2008). We assume that the neutral hydrogen fraction takes a value $\Omega _{{\rm H}\,\scriptsize {\rm I}} ={5\times 10^{-4}}$ (Masui et al. 2013).

B.1. Galactic Synchrotron

In order to test our methods we require simulated maps of the Galactic emission from our own Galaxy in the range 400–1400 MHz with 1 MHz resolution. Though there are maps at both 800 MHz and 1420 MHz, the only public all-sky radio survey in this range is the 408 MHz Haslam map (Haslam et al. 1982). However, the Global Sky Model (de Oliveira-Costa et al. 2008) is based on a compilation of maps from 10 MHz to 94 GHz. We use the Global Sky Model to generate maps at both 400 MHz and 1420 MHz, and use these to estimate an effective spectral index at each location on the sky

$$\begin{equation} \alpha (\hat{\boldsymbol{n}}) = \frac{\log {T_{1420}(\hat{\boldsymbol{n}})} - \log {T_{400}(\hat{\boldsymbol{n}})}}{\log {1420} - \log {400}} \;. \end{equation} \tag{ B1 }$$

By combining this with the Haslam map¹⁵ we can extrapolate to simulate a map of the sky at any desired set of frequencies

$$\begin{equation} T_{\rm base}(\hat{\boldsymbol{n}}, \nu) = T_{408}(\hat{\boldsymbol{n}}) \left(\frac{\nu }{{408}\,{{\rm MHz}}}\right)^{\alpha (\hat{\boldsymbol{n}})} \;. \end{equation} \tag{ B2 }$$

Unfortunately this map lacks both the small scale angular fluctuations (because of the limited resolution of the Haslam map) and any spectral variations (because of the power law extrapolation) that would be present on the real sky. It is essential to include these to accurately test any foreground removal method.

To include these fluctuations we could add Gaussian realizations of Equation (A1) (with the galactic synchrotron parameters, see Table 1) to the base map, which contain frequency and angular fluctuations at arbitrary resolutions. However the Haslam map already constrains what the sky looks like on scales ≳ 1°, and the extrapolation with the spectral index map is a constraint on the sky at 1420 MHz (on scales larger than 51, the resolution of the Global Sky Model). Therefore, we would like the combined simulated map to be consistent with these observations. We can do this by constraining the realizations to ensure there are no fluctuations on the scales constrained. In practice, we do this by manipulating the amplitudes of the two highest valued eigenmodes of C_l(ν, ν') (from Equation A1) in each realization, to ensure that the 408 MHz and 1420 MHz slices are zero when smoothed on 1° and 51 scales, respectively.

A further problem is that we know the amplitude of small-scale fluctuations varies over the sky, however our realizations are statistically isotropic. This is clearly demonstrated in the analysis of La Porta et al. (2008), which shows that the amplitude of the angular power spectrum traces the galactic structure. To reproduce this we use the rms amplitude of fluctuations across the Haslam map in ∼4° patches (corresponding to Healpix pixels with NSIDE = 16, giving 3072 pixels covering the sky), to rescale the fluctuations. In particular this generates an angular power spectrum on the whole sky which is consistent with a single power-law even when crossing through the beam-scale of the Haslam map into the simulated fluctuations. We do not include variations of the power law index of the angular power spectrum as there appears to be no structure to the small variations found in La Porta et al. (2008).

B.2. Extra-galactic Point Sources

B.3. 21 cm Signal