Low frequency specific problems

The visibility measured by a properly calibrated interferometer is given by

$$V(u,v,w)=\iint I(l,m)P(l,m)e^{-2\pi\iota (ul+vm+w(\sqrt{1-l^2-m^2}-1))} {dldm \over \sqrt{1-l^2-m^2}}$$ (7.1)

where ($u,v,w$) are the antenna co-ordinates in the $uvw$ co-ordinate system, $l,m,$ and $n$ are the direction cosines in this system, $I(l,m)$ is the source brightness distribution (the image) and $P(l,m)$ is the far field antenna reception pattern (the primary beam). For further analysis we will assume $P=1$, and drop it from all equations (for typing convenience!), remembering all along that the effect of $P(l,m)$ is to limit the part of sky from which the antenna can receive radiation to $\sim\lambda/D$ radians, where $D$ is the diameter of the antennas.

For a small field of view ( $l^2+m^2 « 1$), the above equation can be approximated well by a 2D Fourier transform relation. The other case in which this is an exact 2D relation is when the antennas are perfect aligned along the East-West direction. Here, we discuss the problem of mapping with non-East-West arrays. The derivation presented here of results for devising algorithms used for imaging large fields of view presented here follow the treatment of Cornwell & Perley (1992) and Cornwell & Perley (1999).

Mapping with non co-planar arrays

Equation 4.1 also reduces to a 2D relation for a non-East-West array, if the integration time is sufficiently small (snapshot observations). However modern arrays are designed to maximize the uv-coverage with the antennas arranged in a 'Y' shaped configuration (non East-West arrays) (Mathur1969). Fields, such as the ones observed for this dissertation, with emission at all angular scales, require maximal uv-coverage. Telescopes such as the GMRT use the rotation of earth to improve the uv-coverage and observations of complex fields typically last for several hours. Hence, Equation 4.1 needs to be used to map the full primary beam of the antennas, particularly at low frequencies.

Let $n=\sqrt{1-l^2-m^2}$ be treated as an independent variable. A 3D Fourier transform of $V(u,v,w)$ can be written using ($u,v,w$) and ($l,m,n$) as a set of conjugate variables as

$$F(l,m,n) = \iiint{V(u,v,w) e^{2\pi\iota (ul+vm+wn)} du dv dw}$$ (7.2)

Substituting for $V(u,v,w)$ from Equation 4.1 we get

$$F(l,m,n) = \iint \Big\{ {I(l^\prime,m^\prime) \over \sqrt{1-l^{\prime^2}-m^{\prime^2}}}$$

$$\iiint{e^{-2\pi\iota (u(l^\prime-l)+v(m^\prime-m)+ w(\sqrt{1-l^{\prime^2}-m^{\prime^2}}-n-1))}} du dv dw$$

$$\Big\} dl^\prime dm^\prime$$ (7.3)

The $-1$ in the coefficient of $w$ in Equation 4.1 comes from fringe rotation. In the above equation, it corresponds to a shift by one unit along the $n-$axis. Since this is only a change of the origin of the co-ordinates system, it can be absorbed in the $n$ without loss of generality. Using the general result

$$\begin{split}\delta(l^\prime-l) = &\int e^{-2\pi\iota u(l^\pr... ...1 & l=l^\prime 0 & otherwise \end{array} \right. \end{split}$$ (7.4)

we get

$$\begin{split}F(l,m,n)&=\iint{ {I(l^\prime,m^\prime) \over \sq... ...delta(\sqrt{1-l^2-m^2}-n) \over {\sqrt{1-l^2-m^2}}} \end{split}$$ (7.5)

This equation then provides the connection between the 2D sky brightness distribution given by $I(l,m)$ and the result of a 3D Fourier inversion of $V(u,v,w)$ given by $F(l,m,n)$ referred to as the Image volume. Here after, I use $I(l,m,n)$ to refer to this Image volume.

The effect of including the fringe rotation term ($-2\pi w$) would be a shift of the Image volume by one unit in the conjugate axis ($n-$axis) (shift theorem of Fourier transforms; Bracewell1986and later eds). Hence, the effect of fringe stopping is to make the $lm-$plane coincide with the tangent plane at the phase center on the Celestial sphere (the point where the tangent plane touches the Celestial sphere) with the rest of the sphere completely contained inside the Image volume (Fig. 4.1).

Figure 4.1: Graphical representation of the geometry of the Image volume and the Celestial sphere. The point at which the Celestial sphere touches the first plane of the Image volume is the point around which the 2D image inversion approximation is valid. For wider fields, emission at points along the intersection of Celestial sphere and the various planes (labeled here as the Celestial Sphere) needs to be projected to the tangent plane to recover the undistorted 2D image. This is shown for 3 points on the Celestial sphere, projected on the tangent plane.

Figure 4.2: Graphical illustration to compute the distance between the tangent plane and a point in the sky at an angle of $\theta$.

Noting that the third variable $n$ of the Image volume is not an independent variable and is constrained to be $n=\sqrt{1-l^2-m^2}$, Equation 4.5 gives the physical interpretation of $I(l,m,n)$. Imagine the Celestial sphere defined by $l^2 + m^2 + n^2=1$ enclosed by the Image volume $I(l,m,n)$, with the top most plane being tangent to the Celestial sphere as shown in Fig. 4.1. Equation 4.5 then tells that only those parts of the Image volume correspond to the physical emission which lie on the surface of the Celestial sphere. However, the Image volume will be convolved by the telescope transfer function. The telescope transfer function is the Fourier transform of the sampling function $S$ in the $uvw$ frame (see Chapter 2, page ). The telescope transfer function, referred to as the dirty beam and defined as $B(l,m,n) = {\mathcal{FT}}[S(u,v,w)]$, also defines a volume in the image domain. The dirty image volume defined by the relation $I^d(l,m,n) = I(l,m,n) \star B(l,m,n)$ is a convolution of $I(l,m,n)$ with $B(l,m,n)$. Since the dirty beam is not constrained to be finite only on the Celestial sphere, $I^d$ will be finite away from the surface of the Celestial sphere corresponding to non-physical emission due to the side lobes of $B(l,m,n)$. A 3D deconvolution using the dirty image and the dirty beam volumes will produce a Clean image volume. An extra operation of projecting all points in the CLEAN image-volume along the Celestial sphere onto the 2D tangent plane to recover the 2D sky brightness distribution is therefore required. Graphical representation of the geometry for this is shown in Fig. 4.1.

3D imaging

The most straight forward method suggested by Equation 4.5 for recovering the sky brightness distribution, is to perform a 3D Fourier transform of $V(u,v,w)$. This requires that the $w$ axis be also sampled at the Nyquist rate (Bracewell1986and later eds; Brigham1988and later eds)). For most observations, it turns out that this is rarely satisfied and doing a FFT along the third axis would result into severe aliasing. Therefore, in practice, the Fourier transform on the third axis is usually performed using the direct Fourier transform (DFT) on the un-gridded data.

To perform the 3D FT (FFT along the $u-$ and $v-$ axis and DFT along the $w-$axis) one still needs to know the number of planes needed along the $n-$axis. This can be found using the geometry shown in Fig. 4.2. The size of the synthesized beam along the $n-$axis is comparable to that along the other two directions and is given by $\approx \lambda/B_{max}$ where, $B_{max}$ is the longest projected baseline length. The separation between the planes along $n$ should be $\le\lambda/2B_{max}$. The distance between the tangent plane and a point separated by $\theta$ from the phase center, for small values of $\theta$, is given by $1-\cos(\theta)\approx \theta^2/2$. For a field of view of angular size $\theta$, critical sampling would be ensured if the number of planes along the $n-$axis, $N_n$, is

$$\begin{split}N_n =&B_{max}\theta^2/\lambda =& \lambda B_{max}/D^2 (\mathrm{for \theta=Full  Primary  beam}) \end{split}$$ (7.6)

where $D$ is the diameter of the antenna ( $B_{max} \approx 25$ km for the GMRT at 325 MHz). Therefore, for mapping a $1^\circ$ field of view without distortions, one would require 8 planes along the $n$-axis. However, for mapping with Central Square alone ( $B_{max} \approx 1$ km), one plane is sufficient. At these frequencies, it becomes important to map most of the primary beam since the number and often the intensity of the sources in the field increase and the side-lobes due to these sources limit the dynamic range in the maps. Hence, even if the source of interest is small, it requires a full 3D inversion (and deconvolution).

Another reason why more than one plane would be required for very high dynamic range imaging is as follows. Strictly speaking, the only point which lies in the tangent plane is the point at which the tangent plane touches the Celestial sphere. All other points in the image, even close to the phase center, lie slightly below the tangent plane. Deconvolution of the tangent plane then results into distortions for the same reason as the distortions due to the deconvolution of a point source which lies between two pixels in the 2D case (Briggs1995). As in the 2D case, this problem can be minimized by over sampling the image which, in this case, implies having more than one plane along the $n$-axis, even if Equation 4.6 implies that one plane is sufficient.

Figure 4.3: Approximation of the Celestial sphere by multiple tangent planes (polyhedron imaging).

Polyhedron imaging

As mentioned above, emission from the phase center and from points close to it, lie approximately in the tangent plane. Polyhedron imaging relies on exploiting this by approximating the Celestial sphere by a number of tangent planes, referred to as facets, as shown in Fig. 4.3. The visibilities are recomputed to shift the phase center to the tangent points of each facet and a small region around each of the tangent points is then mapped using the 2D approximation.

The number of planes required to map an object of size $\theta$ can be found simply by requiring that the maximum separation between the tangent plane and the Celestial sphere be less than $\lambda/B_{max}$, the size of the synthesized beam. As shown earlier, this separation for a point $\theta$ degrees away from the tangent point is $\approx \theta^2/2$. Hence, for critical sampling, the number of planes required is equal to the solid angle subtended by the sky being mapped ( $\theta_f^2$) divided by $\theta^2/2$( $=\lambda/2B_{max}$)

$$\begin{split}N_{poly} =& 2\theta_f^2 B_{max}/\lambda =& 2B... ...da/D^2  (\mathrm{for \theta_f = Full primary beam}) \end{split}$$ (7.7)

Notice that the number of planes required is twice as many as required for 3D inversion. However, since a small portion around the tangent point of each plane is used, the size of each of these planes can be small, off-setting the increase in computations due to the increase in the number of planes required. Another approach which is often taken for very high dynamic range imaging, is to do a full 3D imaging on each of the planes. This would effectively increase the size of the field that can be imaged on each tangent plane, thereby reducing the number of planes required.

The polyhedron imaging scheme is implemented in the current version of the AIPS data reduction package and the 3D inversion (and deconvolution) is implemented in the (no longer supported) SDE package developed by T.J. Cornwell et al. Both these schemes, in their full glory, are available in the (recently released) AIPS++ package.

The GMRT 325-MHz data was imaged using the IMAGR task in AIPS. This program implements the polyhedron algorithm and requires the user to supply the number of facets to be used and a list of the locations of the centre of each facet with respect to the image centre and the size of each facet. This list of facets and their parameters were computed using the task SETFC in AIPS, typically resulting in a grid of $5\times 5$ facets, each of size $256 \times 256$ pixels.

Bandwidth Smearing

The effect of a finite bandwidth of observation as seen by the multiplier in the correlator, is to reduce the amplitude of the visibility by a factor given by $\frac{sin(\pi l \Delta \nu/\nu_o \theta)}{(\pi l \Delta\nu / \nu_o \theta)}$ where $\theta$ is the angular size of the synthesized beam, $\nu_o$ is the center of the observing band, $l$ is location of the point source relative to the field center and $\Delta \nu$ is the bandwidth of the signal being correlated.

The distortion in the map due to the finite bandwidth of observation can be understood as follows. For continuum observations, the visibility data integrated over the bandwidth $\Delta \nu$ is treated as if the the observations were made at a single frequency $\nu_o$ - the central frequency of the band. As a result the $u$ and $v$ co-ordinates and the value of visibilities are correct only for $\nu_o$. The true co-ordinates at other frequencies in the band are related to the recorded co-ordinates as

$$(u,v)=\left(\nu_o\frac{u_\nu}{\nu}, \nu_o\frac{v_\nu}{\nu}\right)$$ (7.8)

Since the total weights $W$, used while mapping, does not depend on the frequency, the relation between the brightness distribution and visibility at a frequency $\nu$ becomes

$$V(u,v)=V\left(u_\nu\frac{\nu_o}{\nu}, v_\nu\frac{\nu_o}{\nu}\righ... ...(\frac{\nu}{\nu_o}\right)^2 I\left(l\frac{\nu}{\nu_0},m\frac{\nu}{\nu_0}\right)$$ (7.9)

Hence the contribution of $V(u,v)$ to the brightness distribution get scaled by $(\nu/\nu_o)^2$ and the co-ordinates get scaled by $(\nu/\nu_o)$. The effect of the scaling of the co-ordinates, assuming a delta function for the Dirty Beam, is to smear a point source at position $(l,m)$ into a line of length $(\Delta \nu/\nu_o)\sqrt{l^2 + m^2}$ in the radial direction. This will get convolved with the Dirty Beam and the total effect can be found by integrating the brightness distribution over the bandwidth as given in Equation 4.9

$$I^d(l,m)=\left[ {\int\limits_0^\infty \vert H_{RF}(\nu)\vert^2 \l... ...u \over {\int\limits_0^\infty \vert H_{RF}(\nu)\vert^2 d\nu}} \right]*DB_o(l,m)$$ (7.10)

where $H_{RF}(\nu)$ is the band-shape function of the RF band and $DB_o$ is the Dirty Beam at frequency $\nu_o$. If the synthesized beam is represented by a Gaussian of standard deviation $\sigma_b=\theta_b/\sqrt{8 \mathrm{ln} 2}$ and the bandpass by a rectangular function of width $\Delta \nu$, the fractional reduction in the strength of a source located at a radial distance of $r=\sqrt{l^2+m^2}$ is given by

$$R_b=1.064{\theta_b\nu_o \over {r\Delta\nu}}\mathrm{erf}\left(0.833{r\Delta\nu \over {\theta_b\nu_o}}\right)$$ (7.11)

Equation 4.10 is equivalent to averaging a large number of maps made from quasi-monochromatic visibilities at $\nu$. Since each such map scales by a different factor, a source away from the center would move along the radial line from one map to another, producing the radial smearing mentioned above, convolved with the Dirty Beam.

The effect of bandwidth smearing can be reduced if the band is split into frequency channels with smaller channel widths. This effectively reduces the bandwidth $\Delta \nu$ as seen by the mapping procedure and while gridding the visibilities, $u$ and $v$ can be computed separately for each channel and assigned to the appropriate uv-cell. The FX correlator used in GMRT provides up to 128 frequency channels over the entire bandwidth of observation and the visibilities can be retained as multi-channel in the mapping process to reduce bandwidth smearing. Although purely from the point of view of bandwidth smearing, averaging $\sim10$ channels at 325 MHz would be acceptable, keeping the visibility database with all the 128 channels is usually recommended to allow identification and flagging of narrow bands RFI.