IMAGING THE EPOCH OF REIONIZATION: LIMITATIONS FROM FOREGROUND CONFUSION AND IMAGING ALGORITHMS

The baselines are defined on a local plane whose origin is at the latitude and longitude of the array center. A two-dimensional vector x = [x  y] (in meter units) represents the projected baseline on the local plane. The corresponding baseline u = [u  v] (in wavelengths) at a frequency f is then given by u = x.(f/c), where c is the speed of light. A two-dimensional direction vector l toward any direction is represented by its direction cosines: l = [l  m]. An interferometer observation gives a visibility cube with two spatial axes representing the baseline and a frequency axis. An image cube on the other hand has two axes representing direction cosines on the sky and a third frequency axis. Assuming coplanarity of array elements, within the framework of Fourier synthesis imaging, u and l are Fourier conjugates.

The telescope sampling function is represented as $\bar{W}({\bf x}, f)$, which is the weight assigned to the baseline represented by the vector $\bf {x}$, and includes any visibility taper applied prior to imaging. In general, $\bar{W}({\bf x})$ depends on instrumental, observational, and processing parameters. Among instrumental parameters is the antenna layout that determines the set of available baselines. Observational parameters, which include pointing direction, duration of rotation synthesis, and mode of observation—drift scan or track mode, define the uv coverage. Observational parameters along with the antenna layout determines the natural weight at each x, which is proportional to the number of measured visibilities falling into the pixel centered at x. Among processing parameters are the visibility grid size and extent, convolution kernel used to grid the visibilities, and any visibility mask and/or taper used during imaging. The weighting function $\bar{W}({\bf x}, f)$ may be factored into a frequency-independent geometric component W(x), a frequency-dependent taper component T(x, f), and a frequency window function B(f):

$$\begin{equation} \bar{W}({\bf x}, f)=W({\bf x})\,T({\bf x},f)\,B(f). \end{equation} \tag{ 1 }$$

Instrumental and observational parameters are geometric in nature and hence frequency independent. They may be sufficiently represented by W(x). Processing parameters like visibility taper function and visibility masks are frequency dependent and are represented by T(x, f). B(f) is the window function in frequency assigned to each visibility.

An image made with a weighted baseline distribution $\bar{W}({\bf x}, f)$ has a frequency-dependent PSF $\tilde{A}\left({\bf l}, f\right)$, where l is the two-dimensional vector representing direction on the sky. Tilde notation is used for image domain functions as opposed to their Fourier domain counterparts. Consider the sidelobe contribution at zenith (l = [0  0]) from a source with a flux density of unity and spectral index of zero located at any l. This contribution $\tilde{S}({\bf l},f)$ varies with frequency due to a changing PSF. We wish to compute the energy in these spectral variations on different frequency scales in terms of the telescope weighting function $\bar{W}({\bf x}, f)$.

PSF $\tilde{A}({\bf l},f)$ is the Fourier Transform of the uv-weighting function A(u, f) and is given by

$$\begin{equation} \tilde{A}({\bf l},f)=\frac{\int _{-\infty }^{\infty }du\int _{-\infty }^{\infty }dv\, A({\bf u},f) \exp (-i2\pi {\bf u}.{\bf l})}{\int _{-\infty }^{\infty }du\int _{-\infty }^{\infty }dv\, A({\bf u},f)}, \end{equation} \tag{ 2 }$$

where the denominator is included to obtain $\tilde{A}({\bf l},f)$ as the normalized PSF for all frequencies. Any additional frequency weighting may be absorbed into the frequency window function B(f). We have assumed that A(u, f) includes the conjugate baselines, failing which $\tilde{A}({\bf l},f)$ may be defined as the real part of the integral in Equation (2). The function A is just a scaled version of the function $\bar{W}$ through the scaling relationship given by u = x.(f/c). We may represent this scaling relationship as

$$\begin{equation} A({\bf u},f) = \bar{W}\left({\bf u}\frac{c}{f},f\right) = \bar{W}\left({\bf x},f\right). \end{equation} \tag{ 3 }$$

The quantity uc/f has units of distance and is frequency invariant: it is simply the baseline distance in meters. We use this to make the substitution x = uc/f, y = vc/f and get

$$\begin{equation} \tilde{A}({\bf l},f) = \frac{1}{\mathcal {A_W}}\int _{-\infty }^{\infty } dx \int _{-\infty }^{\infty } dy \bar{W}({\bf x},f) \exp \left(-i2\pi {\bf x}.{\bf l}\frac{f}{c}\right), \end{equation} \tag{ 4 }$$

where $\mathcal {A_W}=\int _{-\infty }^{\infty } dx \int _{-\infty }^{\infty } dy \, \bar{W}({\bf x},f)$ is the area under $\bar{W}$. The importance of this substitution is that it expresses all baselines in units of meter as opposed to wavelength and makes them independent of frequency.

Consider a point source with a flux density of unity and a spectral index of zero located at l. The response at zenith due to a sidelobe on this source is given by

$$\begin{equation} \tilde{S}({\bf l},f) = \tilde{A}({\bf l},f). \end{equation} \tag{ 5 }$$

We are interested in the energy in $\tilde{S}({\bf l},f)$ at different frequency scales. In other words, to evaluate PSF contamination, we have to compute the power spectrum $\tilde{S}({\bf l},\eta)$ by Fourier transforming $\tilde{S}({\bf l},f)$ with frequency offset Δf and delay η as Fourier conjugates. We use Δf instead of f because the spatial EoR fluctuations are evaluated about an origin r₀ corresponding to a frequency f₀ and redshift z₀.

$\tilde{S}(\eta,{\bf l})$ may be expressed as

$$\begin{equation} \tilde{S}({\bf l},\eta) = \int _{-\infty }^{\infty }d(\Delta f)\, \tilde{A}({\bf l},f) \exp (+i2\pi \eta \Delta f), \end{equation} \tag{ 6 }$$

where we have used a positive sign in the exponent as we are transforming a function of an image space variable (f) to a function of a Fourier space variable (η). Using Δf = f − f₀,  f₀ = (f_L + f_H)/2, where f_L and f_H are the lowest and highest frequencies in the instrument bandpass, and writing the conjugate variable as f instead of Δf we get

$$\begin{equation} \tilde{S}({\bf l},\eta) = \exp (-i2\pi \eta f_0)\int _{-\infty }^{\infty } df\nonumber \tilde{A}({\bf l},f) \exp (+i2\pi \eta f). \end{equation} \tag{ 7 }$$

Substituting from Equation (4) for $\tilde{A}({\bf l},f)$, from Equation (1) for $\bar{W}({\bf x},f)$, and interchanging the order of integrations we get

$$\begin{eqnarray} \tilde{S}({\bf l},\eta) &=& \exp (-i2\pi \eta f_0) \int _{-\infty }^{\infty }dx\int _{-\infty }^{\infty }dy\, W({\bf x})\nonumber \\ && \times\int _{-\infty }^{\infty }df\, \frac{T({\bf x},f)}{\mathcal {A_W}}B(f) \exp \left[-i2\pi f\left(\frac{{\bf x.l}}{c}-\eta \right)\right].\nonumber\\ \end{eqnarray} \tag{ 8 }$$

The inner integral, in general, may have to be evaluated numerically. However, to have an insight into the nature of $\tilde{S}({\bf l},\eta)$, we evaluate Equation (8) for the case where there is no taper or mask applied in the uv domain. In this case, T(x, f) ≡ 1, and $\mathcal {A_W}$ is frequency independent. The inner integral is now just a scaled and shifted version of the Fourier transform of the frequency window function. Consider a rectangular frequency window defined by

$$\begin{equation} B(f) = \left\lbrace \begin{array}{@{}rl}\dsty\frac{1}{\mathcal {B}} & \mbox{ if $f\in [f_L,f_H]$} \\\ms 0 &\mbox{ otherwise,} \end{array} \right. \end{equation} \tag{ 9 }$$

where $\mathcal {B}=f_H-f_L$ is the frequency bandwidth. The inner integral now evaluates to

$$\begin{equation} \tilde{B}\left(\frac{{\bf x.l}}{c}-\eta \right) = \frac{\sin \left[\pi \mathcal {B}\left(\frac{{\bf x.l}}{c}-\eta \right)\right]}{\pi \mathcal {B}\left(\frac{{\bf x.l}}{c}-\eta \right)} \exp \left[-i2\pi f_o\left(\frac{{\bf x.l}}{c}-\eta \right)\right]. \end{equation} \tag{ 10 }$$

This gives us the relationship we are looking for

$$\begin{equation} \tilde{S}({\bf l},\eta) = \frac{\exp (-i2\pi \eta f_o)}{\mathcal {A_W}}\int _{-\infty }^{\infty }dx\int _{-\infty }^{\infty }dy W({\bf x})\, \tilde{B}\left(\frac{{\bf x.l}}{c}-\eta \right). \end{equation} \tag{ 11 }$$

The power spectrum $\tilde{S}({\bf l},\eta)$ is a convolution of the weighting function W(x) and the convolving kernel $\tilde{B}$, which is a shifted version of the Fourier transform of the frequency window function B(f). Readers conversant with the Fourier synthesis imaging literature may recognize Equation (9) as the "delay beam" or "fringe washing function" for a rectangular bandpass window function (see Thompson et al. 2001 and Bridle & Schwab 1989 for details) with x.l/c as the geometric delay term and η as the instrumental delay term. The delay beam has traditional been used to evaluate bandwidth smearing in synthesis telescopes. Equation (11) essentially shows that the PSF contamination in delay space is given by the convolution between the telescope sampling function (in meter units) and the delay beam.

Consider a single point source at direction cosine l_o entering the zenith pixel through the PSF sidelobe at l_o. This is depicted in Figure 2. Due to the lf invariance of the PSF, the frequency structure of the PSF in a bandwidth δf is identical to the sidelobe structure of the PSF (evaluated at f_o) in the interval [l_o − δfl_o/(2f_o), l_o + δfl_o/(2f_o)]. This truncated version of the PSF can be extracted by multiplying the PSF with a window centered at l_o and having width δf l_o/f_o. A Fourier transform of the product of these two functions is the convolution of their individual Fourier transforms which are the baseline weighting function W(x) and the delay beam (or convolution kernel) from Equation (11). The sinc-function envelope of the delay beam has a main lobe width equal to 2c/(l_oΔf), which in general is far narrower than the visibility range 2x_max. The rapid variation in the delay beam is a complex sinusoid in x with frequency l_o f_o/c. At any η, the convolution in Equation (11) is sensitive to a small range in W(x) that is centered at x = cη/l_o with a width c/(l_oΔf). This explains the sharply peaked PSF contamination in η centered at η = x_max l_o/c (see Figure 1), where x_max l_o/c is indeed the geometric delay for a source at l_o. Sources from zenith to horizon (l_o = 0 to |l_o| = 1) suffer geometric delays from η = 0 to η = x_max/c respectively, and hence contribute to foreground contamination at different fairly localized values of η. This is shown in Figure 3 for sources at five different locations. The relative magnitude of these impulses is determined by the relative sidelobe levels at the respective directions. For the case depicted in Figure 3, the rectangular shape of W(x) gives a 1/l roll-off in the PSF sidelobe levels which corresponds to a 1/η roll-off in the contamination at different values of η.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We have used the phrase sharply peaked to describe the PSF contamination for a source entering the zenith pixel from the direction l_o. The exact profile is described by the sinc envelope of the delay beam (see Figure 3). Since the maximum possible geometric delay (source at horizon) is x_max/c, the delay space may be divided into two regions. Region $\mathcal {R}_1$: η ∈ [0, x_max/c] defines a region of relatively high foreground contamination. Region $\mathcal {R}_2$: η > x_max/c, however, has significantly lower but finite PSF contamination due to the sidelobes of the delay beam. The contamination in $\mathcal {R}_2$ falls off as 1/η since the sidelobes of a sinc-function envelope of the delay beam have a 1/x-form roll-off. The level of contamination for η > x_max/c may be reduced further by a judicious choice of a bandpass window B(f) that gives a delay beam with low sidelobe levels. $\mathcal {R}_2$ defines an EoR window where we may focus our efforts to detect the EoR. Figure 3 demonstrates a reduction in delay beam sidelobe contamination in $\mathcal {R}_2$ by over three orders of magnitude with a Blackman–Nuttall frequency window function (Nuttall 1981). This reduction in contamination comes at the cost of reduced resolution in η space. The performance of such a frequency window function will rapidly deteriorate as interference ridden frequencies are flagged. The improvement in $\mathcal {R}_2$ contamination falls below an order of magnitude even if 1% of randomly selected channels are flagged. In any case, $\mathcal {R}_2$ defines a region of significantly lower PSF contamination—an EoR window where we may focus our effort to detect the EoR.

The instrumental k space is an effective space to represent regions and extent of foreground contamination. The instrumental k space is two dimensional with the comoving LOS wavenumber k_|| and comoving transverse wavenumber k_⊥ as axes. If k_x, k_y, and k_z are comoving wavenumbers along three spatial axes with origin at some point about which the EoR fluctuations will be estimated, then Morales & Hewitt (2004) have shown that

$$\begin{eqnarray} k_{\perp } &=& \sqrt{k_x^2+k_y^2},\,\,\,\, k_x = \frac{2\pi u}{D_M(z)},\,\,\,\,k_y = \frac{2\pi v}{D_M(z)},\nonumber\\ k_z &=& k_{||}\approx \eta \frac{2\pi H_of_{\rm H\,\mathsc{i}}E(z)}{c(1+z)^2}, \end{eqnarray} \tag{ 12 }$$

where the approximation holds for a small bandwidth (small redshift range). Here, $f_{{\rm H\,\mathsc{i}}}=1420$ MHz is the rest frequency of 21 cm H i line emission, H₀ = 70 km sec⁻¹ Mpc⁻¹ is the present value of the Hubble constant, $E(z)=(\Omega _M(1+z_o)^3+\Omega _k(1+z_o)^2+\Omega _{\Lambda })^\frac{1}{2}$, and D_M(z) is the transverse comoving distance at redshift z. Figure 4 shows a plot of the measurement space (in k_⊥ − k_|| domain) for a Fourier synthesis array with uniform visibility coverage. A Fourier synthesis array is sensitive to measurements in a region in the k space bounded by some maximum and minimum values of k_⊥ and k_||. The highest k_⊥ mode sampled by the instrument is related to the uv extent of the array and is given by

$$\begin{equation} k_{\perp {\rm max}}=\frac{2\pi x_{{\rm max}}f_{{\rm H\,\mathsc{i}}}}{c(1+z)D_M}, \end{equation} \tag{ 13 }$$

while the smallest k_⊥ probed by the instrument is related to the physical extent of the primary antenna element d and is given by

$$\begin{equation} k_{\perp {\rm min}}=\frac{2\pi df_{\rm H\,\mathsc{i}}}{c(1+z)D_M}. \end{equation} \tag{ 14 }$$

The highest value of k_|| sampled by the instrument is related to the frequency resolution and is given by

$$\begin{equation} k_{|| {\rm max}}\approx \eta _{{\rm max}}\,\frac{2\pi H_of_{\rm H\,\mathsc{i}}E(z)}{c(1+z)^2}\,\,,\,\, \eta _{{\rm max}}=\frac{1}{\Delta f}. \end{equation} \tag{ 15 }$$

The lowest value of k_|| sampled by the instrument in a snapshot observation is related to the frequency bandwidth and is given by

$$\begin{equation} k_{|| {\rm min}}\approx \eta _{{\rm min}}\,\frac{2\pi H_of_{\rm H\,\mathsc{i}}E(z)}{c(1+z)^2}\,\,,\,\, \eta _{{\rm min}}=\frac{1}{\mathcal {B}}. \end{equation} \tag{ 16 }$$

Zoom In Zoom Out Reset image size

The final comment in this section is related to the implications of lf invariance of the PSF to our understanding of the location and level of contamination in the instrumental k space. Consider the foreground contamination in the neighborhood (δl ≪ 1) of some image pixel $\mathcal {P}$ due to a point source at some l_o removed from $\mathcal {P}$. Due to the lf invariance, the transverse structure of contamination around $\mathcal {P}$ is identical to the LOS contamination at $\mathcal {P}$. Since δl/l_o = δf/f_o, the energy in transverse fluctuations around $\mathcal {P}$ at any wavenumber k_⊥ will be the same as the energy in LOS fluctuations at $\mathcal {P}$ at wavenumber η = (D_Ml_o)/(2πf_o) k_⊥. Using Equation (12), the (k_⊥,  k_||) pairs with identical contamination in them are defined by $k_{||}=k_{\perp }\,\,\frac{H_o E(z)}{c(1+z)}\,l_o$.

For a uniformly weighted visibility distribution with uniform coverage, there is significant energy only in the transverse mode k_⊥o = 2πx_max f_o/(cD_M). This gives us a corresponding $\eta _o=l_o/f_o\,k_{\perp _o}=x_{{\rm max}}l_o/c$: a result we established in Section 3.1. If $\mathcal {P}$ is the zenith pixel, η_o reaches a maximum value of η_m = x_max/c for a source at the horizon (|l_o| = 1). The corresponding values of k_||o and k_||m may then be obtained from Equation (12). The PSF contamination in this case will lie on the line joining (k_⊥max, k_||m) and (k_⊥max, k_||min) in the instrumental k space. If the visibility plane is not uniformly filled there will be contamination in various values of k_⊥ < k_⊥o. The PSF contamination however will always lie within the shaded triangle in Figure 4. The instrumental k space sampled by the telescope outside of this shaded triangle then defines the EoR window in k_⊥ − k_|| space where we may focus our efforts to detect the EoR.

Fourier synthesis arrays have traditionally formed image cubes in lmf space by (1) computing the baselines (in wavelength units: u = x.f/c) generated by a given antenna distribution at every frequency channel, (2) weighting the visibilities to get desirable PSF characteristics and to modify the spatial filtering, (3) performing gridding convolution on the visibilities to re-grid them on a regular uv grid at every frequency channel, and (4) using an efficient FFT algorithm to compute a dirty sky map at every frequency channel. Among these signal processing steps, gridding convolution has not been considered in the analytical treatment of Section 2. We now describe the influence of gridding convolution on foreground contamination in η space.

A baseline vector x (in meter units) will suffer a displacement of δu = xδf/c (in wavelength units) between adjacent frequency channels in the uv plane. Here, δf is the inter-channel frequency spacing. We call this phenomenon "baseline migration" and depict it in Figure 5. The figures show a one-dimensional visibility axis at two adjacent channels. While grid pixel number 1 and 8 go from being filled to being unfilled between channel n and n + 1, visibility pixels 5 and 12 go from being unfilled to being filled. Pixels such as 2 and 10 continue to be filled but experience a step change in their weighted visibility value. Such inter-channel variations result in step changes in gridded visibility values of pixels within the support of the convolving kernel. Many such step changes in visibilities due to the displacement of many baseline vectors results in a stochastic jitter in the frequency structure of the PSF. PSF jitter due to baseline migration may induce channel to channel fluctuations in flux versus frequency on the image plane, thereby significantly contaminating the EoR window $\mathcal {R}_2$. This PSF jitter due to baseline migration is expected to have high contribution (1) from longer baseline vectors due to the relatively higher quantum of baseline displacement δu and (2) from regions of low uv density. In regions of higher uv density a given uv pixel, during gridding convolution, receives contributions from several neighboring baselines. The visibility contribution to a given uv pixel then depends on the relative displacement of the neighboring baselines, rather than their actual displacement. On the other hand, consider the extreme case of very low uv density as depicted in Figure 5. Here the actual inter-channel displacement of the baseline vectors induces relatively high steps in the post-gridding convolution visibilities. Natural weighting results in significantly lower weight assigned to regions with (1) longer baseline vectors and (2) low uv density and is expected to result in low PSF jitter in frequency. Uniform weighting, on the other hand, increases the weighting for longer baselines and is expected to result in substantially high PSF jitter.

Zoom In Zoom Out Reset image size

The angular power spectrum of Ectic sources has been observed to have a power-law-like form (Blake & Wall 2002) with substantially greater power on small angular scales; the universe is homogeneous and isotropic on large scales. Consequently, short baselines are not expected to have substantially greater response to extragalactic continuum sources. However, the mean sky brightness will have a response at zero spacing (u = 0, v = 0). Energy in various spatial Fourier modes is seldom perfectly isolated and the zero spacing component will leak into short baselines. Consequently, channel to channel changes in the visibility distribution close to u = 0, v = 0 may result in step changes in flux versus frequency at image pixels thereby contaminating higher values of η in $\mathcal {R}_2$. While extragalactic sources are expected to have a smooth distribution on large angular scales, flux from the Galactic synchrotron emission is expected to follow a power law with significant spatial distribution power at large angular scales. Consequently, short baselines close to u = 0, v = 0 will have significant response to our Galaxy and baseline migrations in this part of the uv plane will create significant PSF jitter that may further contaminate $\mathcal {R}_2$.

Fourier synthesis arrays usually have a central uv void as the antenna elements have a minimum spacing due to constraints imposed by their own physical size. Antenna cross-talk considerations often result in array configurations with a central uv void that is several times the antenna size. In such cases, the central uv void spans several uv pixels. We mentioned in Section 4.1 that baseline migrations in regions of high uv density do not contribute significantly to inter-channel PSF jitter. Though short spacing baselines lie in a densely populated part of the uv plane, with increasing frequency, most baseline vectors are displaced radially away from any uv grid point close to u = 0, v = 0, as opposed to other uv grid points where baseline vectors also get displaced toward them. This asymmetric neighborhood near short spacing visibilities results in relatively higher inter-channel step changes in the visibilities on short baseline vectors. To alleviate this problem, a frequency-independent uv mask that is larger than the central uv void may be employed to reduce such step due to asymmetric migrations in short baselines.

We have seen how channel to channel changes in the visibility weighting W(u, f) contaminates the LOS wavenumber dimension. Here, W(u, f) includes the effect of all visibility weighting and tapers used, and W(u, f) determines the PSF via the imaging transform. If the visibility distribution has uniform coverage, maintaining invariance of W(u, f) with frequency is possible through application of a frequency-independent visibility taper function. In such a case, W(u, f) is same at all frequencies and there is no frequency dependence in the PSF. In practice, there are many unfilled uv pixels and Jelić et al. (2008) have suggested that a visibility mask may be used that rejects the uv pixels that are not filled at all (or most) frequencies prior to image formation. At one end, we may exploit the complete instrument data and sensitivity at the cost of increased contamination in $\mathcal {R}_2$, and at the other end we may have a frequency-independent PSF with some loss in information and sensitivity. However, there exists a middle path that assures that no data are rejected while restricting the foreground contamination in $\mathcal {R}_2$—image formation using the CZT (Rabiner et al. 1969). We first present a brief description of CZT.

The Z transform of any finite discrete sequence x_n,  n = 0, 1, 2....N − 1 is defined as

$$\begin{equation} X(z)=\sum _{n=0}^{n=N-1}x_nz^{-n}. \end{equation} \tag{ 17 }$$

The Z transform is a function of the complex variable z. When evaluated at discrete points along the unit circle defined by z_k = exp (i2πk/N), where k = 0, 1, 2....N − 1, Z transform of the finite discrete sequence x_n reduces to its Discrete Fourier Transform (DFT). The CZT is another special case of the Z transform that is evaluated at discrete points on the locus given by z_k = AW^−k where A and W are complex constants. Such a locus describes a spiral in the complex z plane. Note how for A = 1 and W = exp (− i2π/N) the CZT again reduces to a DFT.

Consider a case where |A| = 1 and W = exp (− i2π/(Nα)), where α is a real constant. The locus of points on the complex z plane where the transform is evaluated now lies on a unit circle like in the case of a DFT. However, angular spacing between the points on the locus is now 2π/(nα) as opposed to 2π/N in the case of a DFT. The first point on the locus is also shifted by the argument of the complex constant A. The loci for the two cases are shown in Figure 6. Though |A| = 1, W = exp (− i2π/(Nα)) is a special case of the generalized CZT, we refer to it hereinafter as CZT with a chirp factor α. The only difference between DFT and CZT is the presence of the chirp factor α and an offset due to the presence of A. α simply stretches the locus on which the transform is evaluated and the CZT of a sequence is essentially a DFT of a stretched version of the sequence. We now describe how visibilities can be gridded in meter units and how a chirp factor can be used to compensate for the stretching of baseline vectors with frequency.

Zoom In Zoom Out Reset image size

As described in Section 4.1, the principal cause of channel to channel PSF variations in conventional gridding and imaging processes is baseline migration due to independent visibility gridding at every channel. However, the baseline vector generated by a given antenna pair moves across the visibility grid with change in frequency in a very predictable manner. The baselines in wavelength units at any frequency f are merely the baselines in meter units stretched by a factor α, where α = f/c. In principle, we may grid the visibilities only once in meter units and use the CZT with chirp factor α to form the image at every frequency channel. Given an image extent, the grid size in meter units is chosen to avoid aliasing even at the highest frequency in the visibility cube. Gridding the visibilities only once in meter units eliminates the problem of baseline migrations and consequently ameliorates the problem of contamination in $\mathcal {R}_2$. In summary, a given antenna pair will fall on exactly the same visibility grid point at all frequencies. We then compute the image on an lm grid that is constant over frequency and use the chirp factor α in a CZT to compensate for baseline stretching. Such a CZT-based imaging transform is given by

$$\begin{equation} S({\bf l}_{mn}) = \sum _{i=0}^{i=N-1}\sum _{j=0}^{j=N-1}V_1({\bf x}_{ij},f) \exp \left({-\frac{i2\pi }{N\alpha (f)}\,{\bf x}_{ij}{\bf l}_{mn}}\right), \end{equation} \tag{ 18 }$$

where i, j are uv (meter units) grid pixel indices, m, n are image grid pixel indices, function V₁(x_ij, f) represents the visibilities at frequency f gridded on to the common grid, and α(f) is the chirp factor given by α(f) = c/f. The transform in Equation (18) is similar to a DFT and may be evaluated by the same efficient FFT algorithms with a minor modification to the exponential terms. The CZT-based image formation results in reduced contamination in $\mathcal {R}_2$. In addition, it reduces the computational cost of gridding convolution because the gridding convolution at all frequency channels is done to a single common grid and the same baseline-timeslot to pixel-weight map is used for all frequencies. On the other hand, due to the introduction of α(f) into the exponential (see Equation (18)), the multiplicative exponential factors in the Radix-2 butterflies need to be recomputed at every frequency channel, thereby increasing the computational cost. Hence, there are competing effects which decide the efficacy of image formation using the CZT.

The threshold flux density S_C above which we expect to find one source per beam of full width at half-maximum θ_arcmin at frequency f_GHz is given by Subrahmanyan & Ekers (2002)

$$\begin{equation} S_C = 40\,\theta _{{\rm arcmin}}^{1.67}f^{-0.67}_{{\rm GHz}} \, \mu {\rm Jy}. \end{equation} \tag{ 19 }$$

The stochastic response to sources weaker than S_C contributes to sky pixel rms that is close to S_C. S_C for MWA evaluates to 2.3 mJy and is about three orders of magnitude higher than the expected thermal noise in 8 s snapshot images: MWA snapshot images will be confusion limited. We expect to find one source above 5S_C about in 10 beamwidths. It is common practice to use a cutoff flux corresponding to one source in 25 beamwidths (Hazard & Walsh 1959) if we need to distinguish point sources from a combination of weaker unresolved faint sources for purposes of making a source catalog. For purposes of peeling, this requirement may be relaxed, and we assume successful peeling of all sources above 5S_C = 10 mJy. The expected residual rms confusion noise after successful subtraction of sources above 10 mJy is about 3.6 mJy. Due to the dominance of confusion over thermal noise, we will neglect the effects of thermal noise hereinafter. However, it must be noted that post foreground modeling and successful subtraction, the residual image that will be used for EoR power spectrum estimation is expected to be dominated by thermal noise.

The threshold to which continuum sources may be identified and subtracted also depends on our ability to distinguish spillover flux due to sidelobes on the sky. We now roughly estimate the expected rms fluctuation at any image pixel due to the sidelobes on all the confusing sources in the sky—a quantity we call "sidelobe confusion." If we assume (1) the sidelobes of the PSF to be a zero mean random variable, and (2) the confusion flux in each pixel to be a random variable that is uncorrelated with the sidelobe level, then the approximate rms sidelobe confusion can be estimated fairly easily: C_S = ∑_iP_iS_i where P_i is the sidelobe level at pixel i and S_i is the primary beam weighted true sky flux due to confusing sources at pixel i. Since 〈P_i〉 = 0, the variance in sidelobe confusion is 〈C²_S〉 = ∑_i〈P²_iS_i²〉 = ∑_i〈P²_i〉 〈S²_i〉. The far sidelobes of MWA are expected to have an rms variation of about 0.05% giving 〈P²_i〉 = 0.25 × 10⁻⁶, and the rms confusion is about 3.6 mJy giving 〈S²_i〉 = 13 mJy². The sidelobe level assumed here is true for sidelobes far from the phase center at which sidelobe confusion is being computed. In reality sidelobe levels vary from sky pixel to pixel and a constant value has been assumed here for simplicity. Far sidelobe level has been chosen since far sidelobes result in foreground contamination at higher values of η where we expect to detect the EoR with minimal contamination, and this calculation will present the worst case. While P_i is assumed to be independent of i, S_i is assumed to be the confusion noise weighted by the primary beam gain, and the summation is performed using the expected primary beam of the MWA primary antenna element (full width to first null of 60°) to evaluate sidelobe confusion. The resulting sidelobe confusion from the entire sky is 〈C_S〉 ≈ 1.9 mJy. The similarity in the values of blending confusion and sidelobe confusion is an interesting consequence of the similarity between the number of resolution elements in the image (weighted by the primary beam) and inverse of the sidelobe level of the synthesized beam. This implies that after subtracting bright point sources from MWA snapshot images, the resulting residual image has statistical fluctuations with roughly equal contributions from the stochastic response to the extragalactic sources and the sidelobe levels.

The EoR brightness fluctuations are expected to have an rms of about 25 mK on scales comparable to the MWA 5 arcmin resolution. This corresponds to an rms flux of about 36 μJy. Hence, apart from confusing foreground sources, sidelobe confusion can confuse any attempt to detect the EoR signal through a simple measurement of image plane variance. That said, flux from confusing sources and PSF sidelobes have structure in the frequency domain which can be exploited to detect the EoR fluctuations despite a high rms fluctuation in the image plane. While spectral structure of confusion extragalactic sources follows relatively simple power laws, spectral structure of instrumental sidelobe confusion is more complex as described in Sections 2, 3, and 4. We now quantitatively estimate the foreground contamination in frequency space (and eventually in η space) for the case of snapshot imaging with the MWA.

This subsection describes a simulation to estimate the foreground contamination in η space for the case of snapshot imaging of extragalactic point sources using the MWA. The simulation assumes a frequency-invariant sky model and assumes that all sources brighter than 10 mJy have been successfully subtracted. The sky model is then interpolated into an lm grid with grid size at least twice as fine as the MWA instrument resolution. The grid extent covers the interval [ − 0.5  0.5] in direction cosines that corresponds to the MWA primary beam full width at first null, which is 60° at 150 MHz. Gridding convolution of visibilities is performed using a triangular kernel like the one shown in Figure 5. A simple kernel was chosen purely to illustrate the nature of contamination in η space. At every frequency channel, three separate PSFs are generated using natural weighting, uniform weighting, and uniform weighting with a Gaussian taper (hereinafter referred to simply as Gaussian weighting). The simulated sky in lm is then convolved with the PSF at zenith and the result is summed. This accumulated value at every frequency channel gives the dirty image flux density at the zenith pixel. This zenith pixel flux versus frequency is used to evaluate the spillover of foreground contamination by Fourier transforming to η domain. The simulation was performed for three cases each involving a different gridding and imaging algorithm; we discuss these below.

Case 1 employed independent gridding convolution at every frequency channel, and the dirty image was formed using an efficient FFT algorithm. Figure 7 shows the flux density versus frequency and contamination power in η for this case. We have shown earlier that for η < x_max/c in $\mathcal {R}_1$, the exact structure of spillover depends on the nature of relative weightings given to visibilities. Weighting here includes the effect of uv plane distribution of baselines, user defined weighting, and any tapers applied. For the case of uniformly weighted visibilities with complete uv coverage, the η spillover from a point source is expected to be localized since the convolution in Equation (11) has a significant contribution only due to the discontinuity in the weighting function at |x| = x_max. In the case of MWA snapshot visibilities, the weighting function is not smooth for |x| < x_max and the convolution in Equation (11) has a significant contribution from all |x| < x_max. Though this results in a distributed spillover, a significant part of the spillover energy is still restricted to $\mathcal {R}_1$, and $\mathcal {R}_2$ will, in the absence of baseline migrations, be contaminated only due to the sidelobes of the convolving kernel from Equation (11) as discussed earlier. This is the PSF contamination that is expected to persist in $\mathcal {R}_2$. Any additional channel to channel variations in the PSF will contaminate $\mathcal {R}_2$ beyond the PSF contamination. Such channel to channel variations are evident in the upper panel of Figure 7 in the form of jitter, especially for the cases of uniform weighting and Gaussian weighting. This jitter was found to be due to inter-channel baseline migrations as discussed in Section 4. To quantitatively estimate the amount of gridding contamination due to baseline migration, the Blackman–Nuttall window with very low sidelobes was used as a frequency window function B(f). As described in Section 3.2, such a windowing function will significantly suppress PSF contamination in $\mathcal {R}_2$ (ηc > 500 m), and the dominant residual contamination in $\mathcal {R}_2$ will be due to the inter-channel jitter arising from the gridding process.

Zoom In Zoom Out Reset image size

It is interesting to note that baseline migrations have induced relatively higher jitter for the cases of uniform weighting and Gaussian weighting as opposed to natural weighting. As discussed in Section 4.1 this is due to two reasons. First, for an inter-channel frequency spacing of δf, a baseline vector x (in meter units) suffers an inter-channel migration on the uv plane of length xδf/c (in wavelength units). Second, and more importantly, a relatively higher amount of inter-channel jitter due to baseline migrations arises from regions of lower uv density (longer baseline vectors) compared to regions of higher uv density (smaller baseline vectors). Since uniform weighting, and in this case Gaussian weighting, places higher emphasis on longer baseline vectors as compared to natural weighting, the amount of jitter in Figure 7 goes in decreasing order of magnitude from uniform weighting to Gaussian weighting to natural weighting. Earth rotation synthesis and drift-scan synthesis will increase the uv density and is expected to ameliorate the jitter.

Case 2 involves image formation with the CZT algorithm with a similar image resolution at every frequency. As described in Section 4.4, image formation with the CZT algorithm is expected to lower gridding contamination in $\mathcal {R}_2$ arising from inter-channel jitter due to baseline migrations. Figure 8 shows the LOS contamination for the case of image formation using the CZT algorithm. The jitter in flux density versus frequency (upper plot) has significantly reduced due to the absence of baseline migrations inherent to independent gridding at every frequency channel—an expected outcome of gridding the visibilities in meter units. Similar image resolution at every frequency channel was achieved by scaling the Gaussian taper to have the same form in wavelength units, and restricting the longest baseline vector that was used in the gridding and imaging process to the same value (in wavelength units) at every frequency channel. This resulted in inter-channel steps in the weighting function at the far end of the uv plane that resulted in a small but significant amount of jitter and hence a contamination in $\mathcal {R}_2$. Nevertheless, the contamination in $\mathcal {R}_2$ for this case is reduced by a factor of 1.7, 6.4, and 21.7 for natural weighting, uniform weighting, and Gaussian weighting, respectively. Relative contamination for natural weighting has not improved significantly as regions in the uv plane that contribute most to jitter are down-weighted by natural weighting.

Zoom In Zoom Out Reset image size

Case 3 involves image formation with the CZT algorithm with a frequency-dependent image resolution. In this case, the set of baselines that are used in the gridding and imaging routine are identical at all frequencies. It may be noted that despite a frequency-dependent image resolution, the CZT algorithm computes the sky flux density on the same lm grid at all frequencies. The results of this simulation are shown in Figure 9. The contamination in $\mathcal {R}_2$ has dropped considerably owing to the elimination of channel to channel jitter in the PSF by this method. The extremely low sidelobes of the Blackman–Nuttall window and the absence of baseline migration in the CZT-based imaging algorithm we have proposed herein result in extremely low levels of contamination in $\mathcal {R}_2$. Quantitatively, the contamination in $\mathcal {R}_2$ for this simulation as compared to Case 1 has reduced by more than three orders of magnitude for all types of weighting schemes.

Zoom In Zoom Out Reset image size

It is important to note that the amount of inter-channel jitter is largely dependent on the visibility distribution and in particular on the filling fraction in different parts of the uv plane. Therefore, the relative merits of the three different weighting schemes discussed will depend on the array configuration and observing strategy. It is however expected that the contamination in $\mathcal {R}_2$ will be orders of magnitude lower if gridding and imaging are done using the proposed CZT algorithm along with a good frequency window. The simulations and results in this section are not exhaustive, and are included as a means to appreciate the various factors that influence the frequency dependence of the PSF and illustrate the potential reduction in foreground contamination.