© ESO, 2011

1. Introduction

Instruments such as the EVLA (Perley et al. 2009), ASKAP (Deboer et al. 2009) and LOFAR (de Vos et al. 2009) are among a new generation of broad-band radio interferometers, currently being designed and built to provide high-dynamic-range imaging capabilities. The large instantaneous bandwidths offered by new front-end systems increase the raw continuum sensitivity of these instruments and allow us to measure the spectral structure of the incident radiation across large continuous frequency ranges. Until recently, the primary goal of wideband imaging has been to obtain a continuum image that makes use of the increased sensitivity and spatial-frequency coverage offered by combining multi-frequency measurements. So far, effects due to the spectral structure of the sky-brightness distribution have been considered mainly in the context of reducing errors in the continuum image, without paying attention to the accuracy of a spectral reconstruction. But now, the new bandwidths (100%) are large enough to allow the spectral structure to also be reconstructed to produce a meaningful astrophysical measurement. To do so, we need imaging algorithms that model and reconstruct both spatial and spectral structure simultaneously, and that are also sensitive to various effects of combining measurements from a large range of frequencies (namely varying ranges of sampled spatial scales and varying array-element response functions).

The simplest method of wide-band image reconstruction is to image each frequency channel separately and combine the results at the end. However, single-channel imaging is restricted to the narrow-band uv-coverage and sensitivity of the instrument, and source spectra can be studied only at the angular resolution allowed by the lowest frequency in the sampled range. Also, for complicated extended emission, the single-frequency uv-coverage may not be sufficient to produce a consistent solution across frequency. While such imaging may suffice for many science goals, it does not take full advantage of what an instantaneously wide-band instrument provides, namely the sensitivity and spatial-frequency coverage obtained by combining measurements from multiple receiver frequencies during image reconstruction.

Multi-frequency-synthesis (MFS) (Conway et al. 1990) is the technique of combining measurements at multiple discrete receiver frequencies during synthesis imaging. MFS was initially done to increase the aperture-plane coverage of sparse arrays by using narrow-band receivers and switching frequencies during the observations. Wide bandwidth systems (~10%) later presented the problem of bandwidth-smearing, which was eliminated by splitting the wide band into narrow-band channels and mapping them onto their correct spatial-frequencies during imaging. It was assumed that at the receiver sensitivities of the time, the measured sky brightness was constant across the observed bandwidth. The next step was to consider a frequency-dependent sky brightness distribution. Conway et al. (1990) describe a double-deconvolution algorithm based on the instrument’s responses to a series of spectral basis functions, in particular, the first two terms of a Taylor series. A map of the average spectral index is derived from the coefficient maps. Sault & Wieringa (1994) describe the MF-CLEAN algorithm which uses a formulation similar to double-deconvolution but calculates Taylor-coefficients via a least-squares solution. More recently, Likhachev (2005) re-derives the least-squares method used in MF-CLEAN using more than two series coefficients.

So far, these CLEAN-based multi-frequency deconvolution algorithms used point-source (zero-scale) flux components to model the sky emission, a choice not well suited for extended emission. Cornwell (2008) describes the MS-CLEAN algorithm which does matched-filtering using templates constructed from the instrument response to various large scale flux components. Greisen et al. (2009) describe a method simular to MS-CLEAN, and Bhatnagar & Cornwell (2004) describe the ASP-CLEAN algorithm that explicitly fits for the parameters of Gaussian flux components and uses scale size to aid the separation of signal from noise. We show in this paper that with the MF-CLEAN approach, deconvolution errors that occur with a point-source model are enhanced in the spectral-index and spectral-curvature images because of error propagation effects, and that the use of a multi-scale technique can minimize this.

For high dynamic range imaging across wide fields-of-view, direction-dependent instrumental effects need to be accounted for. Bhatnagar et al. (2008) describe an algorithm for the correction of time-variable and wide-field instrumental effects for narrow-band interferometric imaging. For wide-field wide-band imaging, these algorithms must be extended to include the frequency dependence of the instrument.

In this paper, we describe MS-MFS (multi-scale multi-frequency synthesis) as an algorithm that combines variants of the MF-CLEAN and MS-CLEAN approaches to simultaneously reconstruct both spatial and spectral structure of the sky-brightness distribution. Frequency-dependent primary-beam correction is considered as a post-deconvolution correction step¹. In Sect. 5, we show imaging examples using simulations, multi-frequency VLA data and wideband EVLA data, to illustrate the capabilities and limits of the MS-MFS algorithm.

1.1. Wide-band imaging

We begin with a discussion of how well we can reconstruct both spatial and spectral information from an incomplete set of visibilities sampled at multiple observing frequencies. An interferometer samples the visibility function of the sky brightness distribution at a discrete set of spatial-frequencies (called the uv-coverage). The spatial-frequencies sampled at each observing frequency ν are between and , where u is used here as a generic label for the uv-distance² and b represents the length of the baseline vector (in units of meters) projected onto the plane perpendicular to the direction of the source. The maximum spatial-frequency measured at each frequency defines the angular resolution of the instrument at the frequency (θ_ν = 1/u_max(ν)). The range of spatial-frequencies between u_min at ν_max and u_max at ν_min represents the region where the visibility function is sampled at all frequencies in the band, and there is sufficient information to reconstruct both spatial and spectral structure. The spatial-frequencies outside this region are sampled only by a fraction of the band and the accuracy of a broad-band reconstruction depends on how well the spectral and spatial structure are constrained by an appropriate choice of flux model.

Radio interferometric measurements of wideband continuum emission can be described as one or more of the following situations.

• 1.

  Flat spectrum sources: for a flat-spectrum source, measurements at multiple frequencies sample the same spatial structure, increasing the signal-to-noise of the measurements in regions of overlapping spatial-frequencies, and providing better overall uv-plane filling. The angular resolution of the instrument is given by u_max at ν_max. Standard deconvolution algorithms applied to measurements combined via MFS will suffice to reconstruct source structure across the full range of spatial scales measured across the band.

• 2.

  Unresolved sources with spectral structure: consider a compact, unresolved source (with spectral structure) that is measured as a point source at all frequencies. The visibility function of a point source is flat across the entire spatial-frequency plane. Therefore, even if u_max changes with frequency, the spectrum of the source is adequately sampled by the multi-frequency measurements. Using a flux model in which each source is a δ-function with a specific spectral model (for example, a smooth polynomial), it is possible to reconstruct the spectral structure of the source at the maximum possible angular resolution (given by u_max at ν_max).

• 3.

  Resolved sources with spectral structure: for resolved sources with spectral structure, the accuracy of the reconstruction across all spatial scales between u_min at ν_min and u_max at ν_max depends on an appropriate choice of flux model, and the constraints that it provides. For example, a source emitting broad-band synchrotron radiation can be described by a fixed brightness distribution at one frequency with a power-law spectrum associated with each location. Images can be made at the maximum angular resolution (given by u_max at ν_max) with the assumption that different observing frequencies probe the same spatial structure but measure different amplitudes (usually a valid assumption). This constraint is strong enough to correctly reconstruct even moderately resolved sources that are completely unresolved at the low end of the band but resolved at the higher end. On the other hand, a source whose structure itself changes by 100% in amplitude across the band would break the above assumption (band-limited signals).In this case, a complete reconstruction would be possible only in the region of overlapping spatial-frequencies (between u_min at ν_max and u_max at ν_min), unless the flux model includes constraints that bias the solution towards one appropriate for such sources.

• 4.

  Spectral structure at large spatial scales: the lower end of the spatial-frequency range presents a different problem. The size of the central hole in the uv-coverage of a typical interferometer increases with frequency, and spectra are not measured adequately at spatial scales corresponding to spatial-frequencies below u_min at ν_max. In the extreme case where most of the visibility function lies within this uv-hole, a flat-spectrum large-scale source (for example) can be indistinguishable from a relatively smaller source with a steep spectrum. Additional constraints in the form of short-spacing spectra may be required for an accurate reconstruction.

• 5.

  Frequency dependence of the instrument: array-element responses usually vary with frequency, direction and time. Standard calibration accounts for the frequency and time dependence for the direction in which the instrument is pointing. Away from the pointing-direction, the frequency-dependent shape of the primary-beam is the dominant remaining instrumental effect, and this results in artificial spectral structure in the images. To recover both spatial and spectral structure of the sky brightness across a large field of view, the frequency dependence of the primary beam must be modeled and removed before or during multi-frequency synthesis imaging.

To summarize, just as standard interferometric image reconstruction uses a priori information about the spatial structure of the sky to estimate the visibility function in unmeasured regions of the uv-plane, multi-frequency image reconstruction algorithms need to use a priori information about the spectral structure of the sky brightness. By combining a suitable model with the known frequency-dependence of the spatial-frequency coverage and element response function, it is possible to reconstruct the broad-band sky brightness distribution from incomplete spectral and spatial-frequency sampling.

2. Multi-scale multi-frequency deconvolution

The MS-MFS algorithm described here is based on the iterative image-reconstruction framework described in Rau et al. (2009) and summarized in Appendix A. Sections 2.1 to 2.7 formulate the algorithm and summarize its implementation in the CASA package. Differences between the multi-scale and multi-frequency parts of MS-MFS with the original MF-CLEAN and MS-CLEAN approaches are highlighted in Sects. 3.1 and 3.2.

2.1. Parameterization of spatial structure

An image with multi-scale structure is written as a linear combination of images at different spatial scales (Cornwell 2008). (1)where I^m is a multi-scale model image³, and is a collection of δ-functions that describe the locations and integrated amplitudes of flux components of scale s in the image. N_s is the number of discrete spatial scales used to represent the image and is a tapered truncated parabola of width proportional to s. The symbol  ⋆  denotes convolution.

2.2. Parameterization of spectral structure

The spectrum of each flux component is modeled by a polynomial in frequency (a Taylor series expansion about ν₀). (2)where represents a multi-scale Taylor coefficient image, and N_t is the order of the Taylor series expansion.

These Taylor coefficients are interpreted by choosing an astrophysically appropriate spectral model and performing a Taylor expansion to derive an expression that each coefficient maps to. One practical choice is a power law with a varying index, represented by a second-order polynomial in space. (3)Here, represents an average spectral-index, and represents spectral-curvature. The motivation behind this choice of interpretation is the fact that continuum synchrotron emission is usually modeled (and observed) as a power law distribution with frequency. Across the wide frequency ranges that new receivers are now sensitive to, spectral breaks, steepening and turnovers need to be factored into models, and the simplest way to include them and ensure smoothness, is spectral curvature⁴.

A Taylor expansion of Eq. (3) yields the following expressions for the first three coefficients from which the spectral index  and curvature images can be computed algebraically. (4)Note that with this choice of parameterization, we are using a polynomial to model a power-law, and N_t rapidly increases with bandwidth. A power-series expansion about and will yield a logarithmic expansion (i.e. I vs. log ν) which requires fewer coefficients to represent the same spectrum⁵.

2.3. Multi-scale multi-frequency model

A wideband model of the sky brightness distribution is constructed from Eqs. (1) and (2). A wideband flux component is a spatial basis function (, Gaussian or parabola) whose integrated amplitude follows a Taylor polynomial in frequency. A region of emission in which the spectrum varies with position will be modeled as a sum of these wide-band flux components. The image-reconstruction process simultaneously solves for spatial and spectral coefficients of these flux components.

The image at each frequency can be modeled as a linear combination of Taylor-coefficient images at different spatial scales. (5)Here, N_s is the number of discrete spatial scales used to represent the image and N_t is the order of the series expansion of the spectrum. represents a collection of δ-functions that describe the locations and integrated amplitudes of flux components of scale s in the image of the tth series coefficient.

2.4. Measurement equations

The measurement equations⁶ for a sky brightness distribution parameterized by Eq. (5) are (6) is a vector of n × 1 visibilities measured at frequency ν. w_ν are Taylor-weights (shown in Eq. (5)).  [S_ν] is an n × m projection operator that represents the spatial-frequency sampling function for frequency ν. The image-domain convolution of model with is written as a spatial-frequency-domain multiplication. where is a spatial-frequency taper function. All images are vectors of shape m × 1.

Equation (6) can be re-written to include all frequencies by stacking  [S_ν] for multiple frequencies. Let N_c be the number of frequencies. (7)where V^obs is a vector of nN_c × 1 visibilities,  [S] is an nN_c × m sampling matrix representing the multi-frequency uv-coverage of the synthesis array. is a diagonal nN_c × nN_c matrix of weights, and consists of N_c diagonal blocks each of size n × n and containing .

If the summations over t and s are written as a block-matrix dot-product, the full measurement matrix has the shape nN_c × mN_sN_t. When multiplied by the set of N_sN_t model sky vectors each of shape m × 1, it produces nN_c visibilities.

For N_t = 3,N_s = 2 the measurement equations can be written as follows, in block matrix form. The subscript p denotes the pth spatial scale and the subscript q denotes the qth Taylor coefficient of the spectrum polynomial. (8)

2.5. Normal equations

The normal equations for the system described in Eq. (6) can be written in block matrix form, with each block-row (for scale size s, and Taylor term t) given by (9)Here, each is an m × m block of the Hessian matrix, and  is one of N_sN_t dirty images.  [W^im]is a diagonal matrix of data-weights (and imaging-weights) and is a diagonal matrix containing Taylor-weights . is the Hessian matrix formed using only one frequency channel, and is a convolution operator containing a shifted version of the single-frequency point-spread-function in each row (see Appendix A.2 for details).  [F^†T_sF] and  [F^†T_pF] are also convolution operators with and as their kernels. The process of convolution is associative and commutative, and therefore, is also a convolution operator whose kernel is given by (15)The dirty images on the RHS of Eq. (9) can be written as follows. where is the dirty image formed by direct Fourier inversion of weighted visibilities from one frequency channel.

When all scales and Taylor terms are combined, the full Hessian matrix contains N_tN_s × N_tN_s blocks each of size m × m, and N_t Taylor coefficient images each of size m × 1, for all N_s spatial scales.

The normal equations in block matrix form for the example in Eq. (8) for N_t = 3,N_s = 2 is shown in Eq. (21). The Hessian matrix consists of N_s × N_s = 2 × 2 blocks (the four quandrants of the matrix), each for one pair of spatial scale s,p (the upper indices). Within each quadrant, the N_t × N_t = 3 × 3 matrices correspond to various pairs of t,q (Taylor coefficient indices; the lower indices). This layout shows how the multi-scale and multi-frequency aspects of this imaging problem are combined and illustrates the dependencies between the spatial and spectral basis functions. (21)This is the system of equations to be solved. The spatial-frequency sampling of a real interferometer is always incomplete ( [S] is rank-deficient). Therefore, each Hessian block, and the entire Hessian matrix is singular, and an exact inverse does not exist⁷. An accurate reconstruction can be obtained only via successive approximation (iterative numerical optimization).

2.6. Principal solution

The principal solution⁸ of the normal equations is an approximate solution that can be computed via diagonal approximations of all Hessian blocks. This solution is then used to pick flux components within the minor-cycle of iterative deconvolution.

Each Hessian block is a convolution operator with a shifted version of a point-spread-function in each row (their centers are aligned on the diagonal). A diagonal approximation represents the assumption that the PSFs are δ-functions, or that the amplitudes in the dirty image at the location of a source reflect the true flux of the source. Also, with the assumption of spatially invariant PSFs, all elements on the diagonal within each Hessian block are the same. Therefore, we can reduce each to one number. The full Hessian reduces to an N_tN_s × N_tN_s element matrix  [H^peak], The principal solution can be obtained by inverting  [H^peak] once, and applying it to the dirty image vectors, one pixel at a time. Such a solution will be correct only at the locations of the centers of isolated flux-components and must be augmented with an iterative optimization approach to ensure accuracy. In the case of perfect sampling (where the Hessian blocks are truly diagonal and PSFs are δ-functions), the principal solution will directly give correct images of Taylor-series coefficients.

2.6.1. Properties of  [H^peak]

Some properties of  [H^peak] are worth noting, to understand the numerical stability of this approach and its dependence on the choice of spectral and spatial basis functions (sky model), and spectral and spatial-frequency sampling functions (data and instrument).

• 1.

  There are N_sN_t elements on the diagonal of  [H^peak]. Each is a measure of theinstrument’s sensitivity to a flux component of unit total fluxwhose shape and spectrum are described by one pair of spatial andspectral basis functions ( and ). Each diagonal element is given by(22)Note that the instrument’s sampling function and data-weights are included in this expression.

• 2.

  The off-diagonal elements measure the orthogonality⁹ between the various basis functions, for the given uv-coverage and weighting scheme. They measure the amount of overlap between basis functions in the measurement domain. Smaller values indicate a more orthogonal set of basis functions that the instrument is better able to distinguish between.

• 3.

  The condition number of this matrix (or of blocks within this matrix) will indicate if the chosen set of basis functions and spatial-frequency coverage provide enough constraints to provide a stable solution, and can be used as a metric to choose a suitable basis set. For a simple example, if a 3-term solution is attempted with data from only two distinct frequencies,  [H^peak] will be singular. Or, for some choice of multi-frequency uv-coverage, the visibilities measured by the instrument for two different spatial scales may become hard to distinguish. Then, the cross-term element of  [H^peak] corresponding to this combination could have a higher value, indicating that the two parameters are highly coupled, and there is insufficient information in the data and sampling pattern to distinguish between the scales. A similar situation can arise to create ambiguity between spatial and spectral structure (an extreme example is multi-frequency measurements from only one baseline).

• 4.

  In general,  [H^peak] will be a positive-definite symmetric matrix whose inverse can be easily computed via a Cholesky decomposition¹⁰. The value of N_s is usually  < 10, making the inversion of  [H^peak] tractable as a one-time operation.

• 5.

  Some further approximations can be made about the structure of  [H^peak] to simplify its inversion, and it is important to understand the numerical implications of these trade-offs. One is a block-diagonal approximation of  [H^peak] (i.e. using only those blocks of the Hessian in Eq. (21) for which s = p; top-left and bottom-right quadrants). This approximation treats each spatial scale separately and assumes that the scale basis functions are orthogonal. For each scale, cross-terms between Taylor functions are preserved, and a multi-frequency principal solution can be obtained separately for each spatial scale. Note that a set of tapered truncated paraboloids is never orthogonal, but this separate-scale approximation works because of the iterative χ²-minimization process. This approximation makes the Hessian inversion easier, but to preserve accuracy, the update step of the iterative deconvolution still needs to evaluate the full LHS of the normal equations while subtracting out a flux component.

2.7. MS-MFS algorithm

This section describes an iterative process that solves the normal equations (Eq. (9)) and produces a set of N_t Taylor-series coefficient images at N_s different spatial scales. Appendix B lists the algorithm-steps in pseudo-code format, reflecting the implementation of MS-MFS in the CASA¹¹ software package.

Pre-compute Hessian:

each block of the Hessian is a convolution operator, consisting of a shifted version of the same convolution kernel in each row. Therefore, it suffices to compute and store one kernel per Hessian block. Convolution kernels for all distinct blocks in the N_sN_t × N_sN_t Hessian are evaluated via Eq. (15). All kernels are normalized by the sum-of-weights such that the peak of is unity, and the relative weights between Hessian blocks is preserved. A block-diagonal approximation of  [H^peak] is done, and a set of N_s matrices each of shape N_t × N_t and denoted as are constructed. Their inverses are computed and stored in .

Initialization:

all N_sN_t model images are initialized to zero (or an a priori model to start from).

Major and minor cycles:

iterative image reconstruction in radio interferometry is usually split into major and minor cycles (see Appendix A.3). The major cycles compute the RHS of the normal equations, and the minor cycle inverts the Hessian and gets an estimate of the model. This model is used in the next major cycle to compute new RHS vectors from the residuals, and the process repeats itself until convergence is achieved. Steps 1 and 5 (of the list of steps given below) form one major cycle, and repetitions of Steps through 4 form the minor cycle.

• 1.

  Compute residual images: the RHS vectors (residual or dirtyimages) -1 }  of the normal equations are computed viaEq. (20) by first computing the multi-frequency dirtyimages and then convolving them by the scale basis functions.

• Find a flux component: the principal solution is computed for all pixels, one scale at a time. (23)Here, is a list of N_t Taylor-coefficients for the pixel-location pix and scale s. is the sth block (of size N_t × N_t) in the list of diagonal-blocks of  [H^peak], and is the N_t × 1 vector constructed from -1 }  for one pixel pix. This step is performed on all pixels, separately for all scales s, resulting in N_s sets of N_t Taylor-coefficient images. The most appropriate set of Taylor coefficients must now be chosen. The search is performed across pixels and scales. Many heuristics can be used here. For example, in iteration i, choose the N_t element solution-set with the dominant q = 0 component across all scales and pixel locations. Or, pick the set of components that makes the largest impact on the value of χ². Or, choose the location from the peak in the t = 0 residual image, and compute the principal solution only for that pixel. The result of this step is a set of N_t model images, each containing one δ-function that marks the location of the center of a flux component of shape (p represents the scale of the chosen component, out of all possible values of s). The amplitudes of these N_tδ-functions are the Taylor coefficients that model the spectrum of the integrated flux of this component. Let these N_t model images from iteration i be denoted as .

• 3.

  Update model images: a set of N_t multi-scale model images are accumulated. (24)where g is a loop-gain that takes on values between 0 and 1 and controls the step size for each iteration in the χ²-minimization process.

• 4.

  Update RHS: the RHS residual images¹² are updated by evaluating and subtracting out the entire LHS of the normal equations. However, since the chosen flux component corresponds to just one scale, this becomes a summation over only Taylor terms (and not scale; there is only one scale p in the current model). (25)Repeat from Step 2 until the minor-cycle flux limit is reached.

• Predict: model visibilities are computed from the set of Taylor-coefficient images via Eq. (6). Residual visibilities are computed as . Repeat from Step 1 until a global convergence criterion is satisfied. (After the first iteration, is used as the new in Step 1.)

Restoration:

after convergence, the model Taylor-coefficient images can be interpreted in different ways.

• 1.

  The most obvious data products are the Taylor-coefficientimages themselves, which are directly smoothed by the restoringbeam. Residual images are added back in after computing theprincipal solution from the residuals obtained in the last instanceof Step 1, to ensure that any undeconvolved flux has the correctflux values¹³.

• 2.

  For the study of broad-band radio emission, the spectral coefficients can be interpreted in terms of a power law in frequency with varying index (as described in Sect. 2.2). The data products are images of the reference-frequency flux , the spectral-index and the spectral curvature . Spectral index and curvature images are calculated only in regions where the values in are above a chosen threshold.

• 3.

  An image cube can be constructed by evaluating the spectral polynomial via Eq. (2) for each frequency. This data product is useful for sources whose emission is not well modeled by a power law, but is a smooth polynomial in frequency. Band-limited signals that taper off smoothly in frequency are one example.

• 4.

  When multiple sources along a given line-of-sight have different spectra, the Taylor-coefficients will represent the combined spectrum. To compute spectral index and curvature maps for foreground sources, a polynomial background-subtraction must be done on the Taylor-coefficient images before Eqs. (26) to (28) are evaluated.

Primary-beam correction:

for wide-field imaging, the spatial and spectral structure of the primary beam (array-element response function) contributes to the measured signal. If this instrumental effect is not accounted for, the output Taylor-coefficient images approximately represent the product of the sky and the primary-beam. (29) P_υ0 is the primary beam at the reference frequency, and P_α and Pβ are spectral index and curvature due to the frequency dependence of the primary beam.

A correction for the average primary-beam and its frequency dependence can be done as a post-deconvolution step. The primary-beams are first evaluated or measured as a function of frequency, and the frequency-dependence per pixel modeled by a power-law or a polynomial (perferably the same spectral polynomial used for the image reconstruction). Primary-beam correction can then be done as follows. Note that if a polynomial is fit for the frequency-dependence of the primary beam, and P_υ0,Pα,Pβ computed from it, the above operation is numerically identical to doing a polynomial division in terms of two sets of coefficients (for N_t ≤ 3). A brute-force polynomial-division using more series coefficients will yield a more accurate solution. Note however, that such a correction will not be accurate if there are time-dependent variations in the primary-beam, and will require integration with the AW-Projection algorithm discussed in Bhatnagar et al. (2008).

3. Relation to MF-CLEAN and MS-CLEAN

The MS-MFS algorithm is a combination of the general ideas used in MS-CLEAN and MF-CLEAN, but there are some subtle differences. The next two sections briefly discuss these differences and their numerical implications.

3.1. Relation to MF-CLEAN

The MF-CLEAN algorithm (Sault & Wieringa 1994) models the sky as a collection of point sources with a Taylor polynomial spectrum. A point-source version of MS-MFS can be derived by setting N_s = 1 and using -function in the derivations in Sects. 2.4 and 2.5. The normal equations can be written in block matrix form (for example, for N_t = 3). (33)Each block  [H_t,q] is a convolution operator with as its kernel. On the other hand, the MF-CLEAN algorithm described in Sault & Wieringa (1994) follows a matched-filtering approach using functions called spectral-PSFs, which are equivalent to the convolution kernels from the first row of Hessian blocks (q = 0) in Eq. (34). In MF-CLEAN, the Hessian elements and RHS vectors are calculated by convolving spectral-PSFs with themselves and the residual images. Formally, this matched filtering approach is exactly equal to the calculations shown in Eqs. (34) and (35) only under the conditions that there is no overlap on the spatial-frequency plane between measurements from different observing frequencies, and all measurements are weighted equally across the spatial-frequency plane (uniform weighting). In practice, MF-CLEAN incurs errors for arrays with dense spatial-frequency coverage where tracks from different baselines and frequency-channels intersect. When applied to simulated EVLA data, numerical instabilities limited the fidelity of the final image, especially with extended emission, and this instability was eliminated by changing the computations from Eqs. (36) and (37) to Eqs. (34) and (35). Also, the MF-CLEAN formulation uses a two-term Taylor-polynomial, which can be shown to result in a dynamic-range limit of 10⁴ for a bandwidth ratio of 2:1 and source spectral index of –1.0.

3.2. Relation to MS-CLEAN

The MS-CLEAN algorithm (Cornwell 2008; Greisen et al. 2009) models the sky as a combination of multiscale flux components with no spectral structure.

A narrow-band (or flat-spectrum) version of MS-MFS can be derived by setting N_t = 1 in the derivations in Sects. 2.4 and 2.5. The normal equations can be written in block matrix form (for example, for N_s = 2). The peaks of the convolution kernels from the diagonal blocks of the Hessian are a measure of the sensitivity of the instrument to a particular spatial scale. (38)In MS-MFS, a diagonal approximation of this Hessian is used to compute the principal solution. This is equivalent to normalizing the residual images (RHS vectors) by the sum of weights for each spatial scale, before searching for peaks.

In both existing forms of MS-CLEAN, this normalization is replaced by a scale bias, an empirical term that de-emphasises large spatial scales. The scale bias b_s = 1−0.6   s/s_max used by Cornwell (2008) (where s_max is the width of the largest scale basis function) is a linear approximation of how the inverse of the area under each scale function changes with scale size¹⁴. The algorithm described by Greisen et al. (2009) uses b_s ≈ 1.0/s^2x where x ∈  {0.2,0.7} , to approximate a normalization by the area under a Gaussian, for the case when images are smoothed by applying a uv-taper that tends to unity for the zero spatial-frequency. Both these normalization schemes are approximations of using a diagonal approximation of  [H^peak].

Once we have this understanding, we can see that the full Hessian  [H^peak] (and not just a diagonal approximation) can be inverted to get the normalization exactly right, giving an accurate estimate of the total flux at each scale. This becomes useful for sources that contain overlapping flux components of different spatial scales. However, this solution gives correct values only at the locations of the centers of the flux-components, and introduces large errors in the PSF sidelobes. Therefore, for reasons of stability, a diagonal approximation is still a more appropriate choice (demonstrated on simulated EVLA data).

Another difference lies in the minor-cycle updates. The update steps in MS-MFS and the Cornwell MS-CLEAN evaluate the full LHS of the normal equations to update the smoothed residual images and subtract out flux components within the image domain. This allows each minor cycle iteration to search for flux components across all spatial scales. The update step in the Greisen MS-CLEAN ignores the cross-terms, and performs a full set of minor cycle iterations on one scale at a time. A choice between all these methods will depend on trade-offs between the accuracy within each minor cycle iteration, the computational cost per step, and optimized global convergence patterns to control the total number of iterations.

4. Hybrids of narrow-band and continuum techniques

The preceding sections discussed multi-frequency solutions that used data from all measured frequencies together, to take advantage of the combined spatial-frequency coverage. However, there are some situations where single-channel methods used in combination with multi-frequency-synthesis (and no built-in spectral model) will be able to deliver scientifically useful wide-band reconstructions at the continuum sensitivity.

The basic idea of a hybrid wide-band method is to combine the advantages of single-channel imaging (simplicity and non-dependence on any spectral model) with those of continuum imaging (deconvolution with full continuum sensitivity).

• 1.

  Deconvolve each channel separately upto the single-channelsensitivity limit σ_chan. Only sources brighter than σ_chan will be detected anddeconvolved.

• 2.

  Remove the contribution of bright (spectrally varying) sources by subtracting out visibilities predicted from the model image cube. At this stage, the peak residual brightness is at the level of the single-channel noise limit σ_chan.

• 3.

  Perform MFS imaging (flat-spectrum assumption) on the continuum residuals to extract flux that lies between σ_chan and σ_cont. According to Conway et al. (1990) (and as discussed later in Sect. 6), errors due to this flat-spectrum assumption become visible only above a dynamic range of  ~1000 (for α =  −0.7 and a 2:1 bandwidth ratio). Therefore, as long as the sensitivity improvement between a single-channel and the full band is less than  ~1000, this second step will incur no errors even if the remaining flux has spectral structure. This requirement translates to N_chan < 10⁶, which is usually satisfied¹⁵.

• 4.

  Add model images from both steps, and restore the results. For unresolved sources, it may be appropriate to use a clean-beam fitted for the highest frequency, but in general, to not bias spectral information, all channels should be restored using a clean beam fitted to the PSF at the lowest frequency in the range.

The advantages of this approach are its simplicity, and that it can handle wide-band reconstructions with band-limited signals and spectral-lines. The disadvantages are that the angular resolution of the images and spectral information will be restricted to that given by the lowest frequency (a factor of two worse than what is possible for the 2:1 bandwidth available with the EVLA L-band). Also, the single-frequency spatial-frequency coverage may not be sufficient to unambiguously reconstruct all the spatial structure of interest at all frequencies, which in turn will affect spectra measured from these single-channel images. In some cases, additional constraints can be used. Bong et al. (2006) describe spatio-spectral MEM, an entropy based method in which single-channel imaging is done along with a smoothness constraint applied across frequency.

In general, single-channel methods can be used for wide-band imaging mainly to construct an image of the continuum flux. Only if there is sufficient single-frequency uv-coverage to reconstruct an accurate model of the source structure (for example, fields of isolated point sources), may reliable spectral information also be derived from such an approach.

5. MS-MFS imaging results

This section contains three imaging examples that demonstrate the MS-MFS algorithm’s basic capabilities. Section 6 discusses various sources of error and how they manifest themselves, and Sect. 7 demonstrates the limits to which the algorithm can be reasonably pushed.

5.1. EVLA simulation

Data:

wide-band EVLA observations were simulated for a sky brightness distribution consisting of one point source with spectral index of  − 2.0 and two overlapping Gaussians with spectral indices of  − 1.0 and  + 1.0. The spectral index across the resulting extended source varies smoothly between  − 1.0 and  + 1.0, with a spectral turnover in the region of overlap of the two Gaussians, corresponding to a spectral curvature of approximately +0.5. Data were simulated for the EVLA in C-configuration between 1–2 GHz, for an 8-h observation with measurements 5 min apart. The goal of this test was to assess the ability of the MS-MFS algorithm to reconstruct both spatial and spectral information accurately enough for astrophysical use.

Imaging results:

two MS-MFS imaging runs were done and the results compared. Figure 1 illustrates the algorithm’s performance by comparing several truth images describing the true sky brightness (left column) and reconstructions from a multi-scale (middle-column) and a point-source (right-column) version of MS-MFS. The multi-scale version used a flux model in which N_t = 3 and N_s = 4 with scale sizes defined by widths of 0,6,18,24 pixels, and the point-source version used N_t = 3 and N_s = 1 with one scale function given by the δ-function (to emulate the MF-CLEAN algorithm). From top to bottom, the rows correspond to the intensity image at the reference frequency I₀ = I_υ0, the first-order Taylor-coefficient I₁ = I_αI_υ0, the second-order Taylor-coefficient I₂ = (I _α(Iα − 1)/2 + I_β)I_υ0, the spectral index I_α = I₁/I₀ , and the spectral curvature I_β = (I₂/I₀) − I_α(I_α − 1)/2.

• 1.

  With a multi-scale multi-frequency flux model (MS-MFS) thespectral index across the extended source was reconstructed to anaccuracy of δα < 0.02 with the errors rising at the edges of the source wherethe signal-to-noise ratio decreases. The spectral curvature acrossthe extended source was estimated to an accuracy of δβ < 0.05 in the centralregion with the maximum error of δβ ≈ 0.2 in the regions where thecurvature signal goes to zero and the source surface brightness isalso minimum (the outer edges of the source).

• 2.

  With a multi-frequency point-source model (MF-CLEAN) the accuracy of the spectral index and curvature maps was limited to δα ≈ 0.2,δβ ≈ 0.6. This is because the use of a point source model will break any extended emission into components the size of the resolution element and this leads to deconvolution errors well above the off-source noise level. Error propagation during the computation of spectral index and curvature as ratios of these noisy Taylor-coefficient images leads to high error levels in the result (see Sect. 6.3). This clearly shows the importance of using a multiscale flux model when there is extended emission.

• 3.

  The Gaussian on the left has α =  −1.0 and β = 0.0, and N_t = 3 is not sufficient to model the power-law accurately, leading to a value of β ≈  + 0.1 in the truth image as well as in the reconstruction. However, the Gaussian on the right has α =  + 1.0 and β = 0.0 (a straight line), and N_t = 3 is more than sufficient to model it, leading to an accurate value of β ≈ 0 in the truth image and MS-MFS reconstruction. The point-source has α =  −2.0, and the error on β is proportionally higher. The use of N_t > 3 reduces these errors.

5.2. Multi-frequency VLA observations of Cygnus-A

Objective:

multi-frequency VLA observations of the bright radio galaxy Cygnus A were used to test the MS-MFS algorithm on real data and to test standard calibration methods on wide-band data. Cygnus A is an extremely bright (1000 Jy) radio galaxy with a pair of bright compact hotspots (α =  −0.5) located about 1 arcmin away from each other on either side of a very compact flat-spectrum core, and extended radio lobes associated with the hotspots with broad-band synchrotron emission at multiple spatial scales (α =  −0.6 ~  − 1.0) (Carilli & Barthel 1996). The best existing images of Cygnus-A and its spectral structure have been from large amounts of multi-configuration narrow-band VLA data (Carilli et al. 1991) designed to measure the spatial structure as completely as possible at two widely separated frequencies (1.4 and 4.8 GHz).

The goal of this test was to use multi-frequency snapshot observations of Cygnus A to evaluate how well the MS-MFS algorithm is able to reconstruct both spatial and spectral structure from measurements in which the single-frequency uv-coverage is insufficient to accurately reconstruct all the spatial structure at that frequency.

Data were taken in April 2009 when the VLA was transitioning to the EVLA. 15 out of 27 antennas had new wideband EVLA feeds, but the correlator was that of the VLA (narrow-band). Wideband data were taken as a series of narrow-band snapshots spread across 8 h and 9 distinct frequencies across the L-Band (30 min per frequency tuning). Flux calibration at each frequency was done via observations of 3C 286. Phase-only calibration was done using an existing narrow-band image of Cygnus A at 1.4 GHz (Carilli et al. 1991) as a model, and further self-calibration was done with the wide-band flux model derived from the MS-MFS algorithm. The final dataset used for imaging consisted of 9 spectral windows each of a width of about 4 MHz and separated by about 100 MHz.

These data were imaged using two methods, the MS-MFS algorithm with N_t = 3 and N_s = 10, and a hybrid of single-channel imaging followed by MFS on the residuals. Note that these observations do not have dense single-frequency uv-coverage, and the purpose of applying the hybrid method was to emphasize the errors that can occur if this method is used inappropriately. Due to the small angular size of Cygnus-A, the effect of the L-band primary-beam was ignored in both runs.

Results:

Figure 2 shows the reconstructed total-intensity images (top), the spectral-index map obtained via MS-MFS (second row) and the spectral-index map constructed from single-subband maps (third row). For comparison, the image at the bottom is a spectral-index map constructed from existing narrow-band images at 1.4 and 4.8 GHz, each constructed from a combination of VLA A, B, C and D configuration data (Carilli et al. 1991).

• 1.

  The total intensity image (top row) has a peak brightness of77 mJy/beam at the hotspots and a peak brightnessof about 400 mJy/beam for the fainter extendedparts of the halo. Both the methods (MS-MFS and hybrid) gavevery similar total-intensity images and residual images. Theon-source and off-source residuals for the MS-MFS algorithmare 30 mJy and 25 mJy and with thehybrid algorithm are 50 mJy and30 mJy respectively.

• 2.

  The spatial structure seen in the MS-MFS spectral index image (second row) is very similar to that seen in the two-point (1.4–4.8 GHz) spectral index image (bottom). This shows that despite having a comparatively small amount of data (20 VLA B-configuration snapshots at 9 frequencies) the use of an algorithm that models the sky brightness distribution appropriately is able to extract the same astrophysical information as traditional methods applied to large amounts of data (full synthesis runs in multiple VLA configurations at two frequencies). The estimated errors on the spectral index map are  < 0.1 for the brighter regions of the source (near the hotspots) and  ≥ 0.2 for the fainter parts of the lobes and the core. In contrast, the Hybrid spectral index map (third row) clearly shows errors arising due to non-unique solutions at each separate frequency (due to insufficient narrow-band spatial-frequency coverage) as well as smoothing to the angular resolution at the lowest frequency.

• 3.

  The MS-MFS spectral curvature map (not shown here) contains values corresponding to  △ α ≈ 0.4 for the brightest regions on the source. Such a change in α appears high, and we did not have sufficient single-frequency data at other bands to verify these values (the next section contains an example where we could verify the curvature maps with other data). Also, the values of curvature changed between imaging-runs with N_t = 3 and N_t = 4, suggesting over-fitting errors due to the low signal-to-noise ratio of the curvature measurement (which could arise from inaccurate bandpass calibration).

5.3. Multi-frequency observations of M 87

A similar observation of M 87 was done with the goal of measuring the 1–2 GHz spectral index in different parts of the inner core-jet-lobe system and outer halo filaments.

M 87 is a bright (200 Jy) radio galaxy located at the center of the Virgo cluster. The spatial distribution of broad-band synchrotron emission from this source consists of a bright central region (spanning a few arcmin) containing a flat-spectrum core, a jet (with known spectral index of  −0.55) and two radio lobes with steeper spectra (−0.5 > α >  −0.8) (Rottmann et al. 1996; Owen et al. 2000). This central region is surrounded by a large diffuse radio halo (7 to 14 arcmin) with many bright narrow filaments (≈10′′ × 3′).

Multi-frequency VLA data were taken similar to the observations of Cygnus-A, with a series of 10 snapshots at 16 different frequencies within the sensitivity range of the EVLA L-band receivers. These data were imaged using MS-MFS with N_t = 3 and N₂ = 12. The top row of images in Fig. 3 shows the intensity, spectral index and spectral curvature maps of the bright central region at an angular resolution of 3 arcsec (C+B-configuration).

• 1.

  The peak brightness at the center of the final restored inten-sity image was 15 Jy with an off-source rms of1.8 mJy and an on-source rms of between 3 and10 mJy. The spectral index map¹⁶ of the bright central region(at 3 arcsec resolution) shows a near flat-spectrumcore with α_LL =  −0.25, a jet with α_LL =  −0.52 and lobes with  −0.6 > α_LL >  −0.7. This bright central region hadsufficient (>100) signal-to-noise to be able to detect spectral curva-ture, and no obvious deconvolution errors. However the error barson the spectral curvature are at the same level as the measurementitself, and a reliable estimate can only be obtained as an averageover this entire bright region. The average curvature is measuredto be β_LL =  −0.5 which corresponds to a change in α across L-band by .

• 2.

  To verify consistency of this spectral-curvature value with the data, each of the 16 spectral-windows was imaged separately, and restored with the same clean-beam. The plot at the bottom of Fig. 3 shows this integrated flux spectrum (log (I) vs. log (ν/ν₀)) as round dots, along with the average spectrum calculated by MS-MFS (curved line), and straight dashed and dotted lines that correspond to constant spectral indices of  − 0.43,  − 0.53 and  − 0.63. These plots show that a change in α of about 0.2 is consistent with the data. Note that the scatter seen in the single-spectral-window data points is at the 1% level of the source flux. This illustrates the accuracy at which bandpass calibration must be done in order to measure a physically plausible spectral-curvature signal across a 2:1 bandwidth.

• 3.

  These numbers were further compared with two-point spectral indices computed between 327 MHz (P-band), 1.4 GHz (L-band), and 4.8 GHz (C-band) from existing images (Owen et al. 2000; Owen F., priv. comm.). Across the bright central region, two-point spectral indices are  −0.36 > α_PL >  −0.45 and  −0.5 > α_LC >  −0.7. The measured values from our experiment (−0.4 > α_LL >  −0.7 and  △ α ≈ 0.2 across L-band) are consistent with these independent calculations.

6. Sources of error

There are various sources of error that can affect the imaging process, and leave artifacts both on-source and off-source. As with any image-reconstruction algorithm, signs of these errors must be looked-for in the output images before astrophysical interpretation.

6.1. On-source polynomial-fit errors

The errors on the polynomial coefficients and quantities derived from them will depend on the number of measurements of the spectrum, their distribution across a frequency range, the signal-to-noise ratio of the pattern being fitted, and the order of the polynomial used in the fit. These errors will affect regions in the image both on-source and off-source, and the resulting error patterns and their magnitudes will depend on the available uv-coverage, and the choice of reference frequency. Although the physical parameters I_υ0,I_α and I_β can be obtained from the first three coefficients of a Taylor expansion of a power-law with varying index (Eqs. (26) to (28)), a higher order polynomial may be required during the fitting process to improve the accuracy of the first three coefficients¹⁷.

Figure 4 summarizes the errors obtained when the order of the polynomial chosen for imaging is not sufficient to model the power-law spectrum of the source. Data for 8 h synthesis runa with EVLA uv-coverages were simulated for 5 different frequency ranges around 2.0 GHz. The sky brightness distribution used for the simulation was one point source whose flux is 1.0 Jy and spectral index is –1.0 with no spectral curvature. The bandwidth ratios¹⁸ for these 5 datasets were 100%(3:1), 66%(2:1), 50%(1.67:1), 25%(1.28:1), 10%(1.1:1).

MS-MFS was repeated on all these datasets with N_t = 1 to N_t = 7, and one scale N_s = 1, a δ-function. All these datasets were imaged using a maximum of 10 iterations, a loop-gain of 1.0, natural weighting and a flux threshold of 1.0 µJy. No noise was added to these simulations (in order to isolate and measure numerical errors due to the spectral fits). The top left panel in Fig. 4 shows the peak residuals, measured over the entire 0th order residual image. The other three panels show on-source errors for I_υ0 (top right), I_α (bottom left) and I_α (bottom right). Errors on I_υ0,I_α,I_β were computed at the location of the point source by taking differences with the ideal values of I_υ0 = 1.0,α =  −1.0,β = 0.0.

One noticeable trend from these plots is that with sufficient signal-to-noise in the measurements, all errors decrease exponentially (linearly in log-space) as a function of increasing order of the polynomial, and as a function of decreasing total bandwidth. As expected, for very narrow bandwidths, the use of high-order polynomials increases the error. Also, when the order of the polynomial used is too low, the peak residuals are much smaller than the on-source error incurred on I_υ0, I_α and I_β. These trends are based on one simple example, and represent the best-case scenario in which all sources can be described as point sources. For extended emission, there can be additional errors due to deconvolution artifacts. In the case of very noisy spectra, or at a late stage of image reconstruction when the signal-to-noise ratio in the residual images is low, errors can arise from attempting to use too many terms in the polynomial fit.

6.2. Off-source errors and dynamic-range limits

Consider the errors on the continuum image when spectral structure is ignored during MFS imaging (N_t = 0; flat-spectrum assumption). Spectral structure will masquerade as spurious spatial structure, leading to error patterns that resemble the instrument’s response to the first-order term in the Taylor-series expansion. A rough rule of thumb for EVLA uv-coverages is that for a point source with spectral index α =  −1.0 measured between 1 and 2 GHz, the peak error obtained if the spectrum is ignored is at a dynamic range of  < 10³. In other words, if the dynamic-range allowed by the data (peak brightness/thermal noise) is 800 (for example), there will be no visible artifacts if the spectral structure upto α =  −1.0 is ignored. These numbers are consistent with error estimates derived in Conway et al. (1990) that predict errors at the level of Iα/X where X = O(500) is a factor that depends on the uv-coverage of the instrument and the choice of reference frequency. Therefore, if the only goal is to obtain a continuum image over a narrow field-of-view, it may be possible to achieve the maximum-possible dynamic range by dividing out an average spectral index (one single number over the entire sky) from the visibilities before imaging them via MFS with a flat-spectrum assumption (Conway et al. 1990).

At high dynamic-ranges, off-source errors trace the spectral response functions for higher-order Taylor terms. As long as they are visible above the noise, they can be eliminated with a higher-order polynomial fit. Figure 5 shows a set of four images of the 3C 286 field made from wideband EVLA data taken in October 2010. These data consist of four EVLA snapshots spread across 90 min, and cover a frequency range of 1.02 GHz to 2.1 GHz (800 MHz usable bandwidth after accounting for radio-frequency-interference). 3C 286 is 14.4 Jy/beam at 1.5 GHz, with a spectral index of –0.47. The four panels correspond to MS-MFS imaging runs with N_t = 1 (top left), N_t = 2 (top right), N_t = 3 (bottom left), N_t = 4 (bottom left), all with N_s = 1 (a δ-function). All images are displayed with the same grey-scale levels, and show the levels at which the error patterns appear. The dynamic range (calculated as the peak brightness to the peak residual measured near the brightest peak) ranges from 1.6 × 10³ when spectral structure is ignored (N_t = 1), to 1.1 × 10⁵ with N_t = 4.

6.3. Propagation of multi-scale errors

Deconvolution errors contribute to the on-source error in the Taylor coefficient images, and these errors propagate to the spectral index map which is computed as a ratio of two coefficient images (Iα = I₁/I₀). Errors that result when a point-source flux model is used for extended emission can increase the error bars on the spectral index and curvature by an order of magnitude (as demonstrated by the example shown in Fig. 1). These errors are approximately given by (39)where  △ I₀ and  △ I₁ are the absolute errors measured in the first two Taylor-coefficient images. For the example shown in Fig. 1, the errors on the coefficient images were measured as the rms of the deviation from the truth images within the region of high signal-to-noise. These errors result in a prediction of  △ Iα ≈ 0.05 for the multi-scale version, and  △ Iα ≈ 0.7 for the point-source version, which approximately matches the rms of the deviation of the output α image from the truth α image.

6.4. Frequency dependence of the primary beam

When imaging wide fields-of-view, sources away from the pointing center will be attenuated by the value of the primary beam at each frequency. The EVLA primary-beams across a 2:1 bandwidth contribute an extra spectral-index of –1.4 at the half-power point (measured from simulated beams, as well as a Gaussian approximation of the main lobe of the primary beam (Sault & Wieringa 1994)). Therefore, even if a source has a flat spectrum, this artificial spectral index can result in imaging artifacts at the levels described in Sect. 6 (i.e. a dynamic range limit of  ~10³ for a flat spectrum source at the HPBW, due only to the spectral variation of the primary beam). This effect can be corrected as a post-deconvolution step, or by using wide-field imaging algorithms along with MS-MFS. So far, accurate post-deconvolution corrections have been demonstrated out to the 10% level of the highest-frequency primary-beam.

7. Imaging performance (non-standard conditions)

This section describes a set of simulations that test the limits of MS-MFS for different types of source structure. Three cases are studied; sources unresolved at the lowest sampled frequency and resolved at the upper end, sources whose visibility functions lie mostly within the central unsampled region of the uv-plane near the origin, and band-limited signals.

7.1. Moderately resolved sources

Consider a source with broad-band continuum emission and spatial structure that is unresolved at the low-frequency end of the band and resolved at the high-frequency end. Traditionally, spectral information would be available only at the angular resolution offered by the lowest observed frequency. With MS-MFS, the intensity distribution as well as the spectral index of such emission can be imaged at the angular resolution allowed by the highest frequency in the band. This is because compact emission has a signature all across the spatial-frequency plane and its spectrum is well sampled by the measurements. The highest frequencies constrain the spatial structure, and the flux model (in which a spectrum is associated with each flux component) naturally fits a spectrum at the angular resolution at which the spatial structure is modeled. Such a reconstruction is model-dependent and will be accurate only if the spectrum at the smallest measured scales can really be represented as a polynomial.

EVLA simulation:

wide-band EVLA data were simulated for the D-configuration across a frequency range of 3.0 GHz with 6 frequency channels between 1 and 4 GHz (600 MHz apart). This wide frequency range was chosen to emphasize the difference in angular resolution at the two ends of the band (60 arcsec at 1 GHz, and 15 arcsec at 4.0 GHz). The sky brightness chosen for this test consists of a pair of point sources separated by a distance of 18 arcsec (about one resolution element at the highest frequency), making this a moderately resolved source. These point sources were given different spectral indices ( + 1.0 for the top source and  − 1.0 for the bottom one). Figure 6 shows the six single-channel images of this source. At the low frequency end, the source is almost indistinguishable from a single flux component centered at the location of the bottom source whose flux peaks at the low-frequency end. The source structure becomes apparent only in the higher frequencies where the top source (with a positive spectral index) is brighter.

MSMFS Imaging results:

these data were imaged using the MS-MFS algorithm with N_t = 3 and N_s = 1 with only one spatial scale (a δ-function), and Fig. 7 shows the results. The intensity distribution, spectral index and curvature of this source were recovered at the angular resolution allowed by the 4.0 GHz samples (15 arcsec). This demonstrates that an appropriate flux model can constrain the solution accurately at the angular resolution given by the highest sampled frequency.

7.2. Emission at large spatial scales

Consider a very large (extended) flat-spectrum source whose visibility function falls mainly within the central hole in the uv-coverage. The size of this central hole increases with observing frequency. The minimum spatial-frequency sampled per channel will measure a decreasing peak flux level as frequency increases. Since the reconstruction below the minimum spatial-frequency involves an extrapolation of the measurements and is un-constrained by the data, these decreasing peak visibility levels can be mistakenly interpreted as a source whose amplitude itself is decreasing with frequency (a less-extended source with a steep spectrum). Usually, a physically-realistic flux model suffices as a constraint, but with the MS-MFS model (polynomial spectra associated with 2D Gaussian-like components), a large flat-spectrum source and a smaller steep-spectrum source are both allowed and considered equally probable. This creates an ambiguity between the reconstructed scale and spectrum that cannot always be resolved directly from the data, and requires additional information (a low-frequency narrow-band image at the reference frequency to constrain the spatial structure, or low-resolution spectral information).

EVLA simulation:

wide-band EVLA data were simulated for the D-configuration across a frequency range of 3.0 GHz centred at 2.5 GHz. (6 frequency channels located 600 MHz apart between 1.0 and 4.0 GHz). The size of the central hole in the uv-coverage was increased by flagging all baselines shorter than 100 m and the wide frequency range was chosen to emphasize the difference between the largest spatial scale measured at each frequency (0.3 kλ or 10.3 arcmin at 1.0 GHz, and 1.3 kλ or 2.5 arcmin at 4.0 GHz). The sky brightness chosen for this test consists of one large flat-spectrum (α = 0.0) 2D Gaussian whose FWHM is 2.0 arcmin (corresponding to 1.6 kλ, at the reference frequency of 2.5 GHz), and one steep spectrum point-source (α =  −1.0) located on top of this extended source at 30 arcsec away from its peak.

MS-MFS imaging results:

these data were imaged using the MS-MFS algorithm with N_t = 3 and N_s = 3 with scale sizes given by [0, 10, 30] pixels. Two imaging runs were performed with these parameters, without and with short-spacing information. Figure 8 shows images of the intensity (left), spectral index (middle) and residuals (right) for these two runs, and Fig. 9 shows the visibility amplitudes present in the simulated data (left column) as well as in the reconstructed model (right column) at each of the 6 frequencies. In both figures, the top and bottom rows correspond to runs without and with short-spacing information respectively and each pair (in Fig. 8) are displayed at the same intensity scale.

• 1.

  The first imaging run used only baselines longer than100 m to simulate a large central hole in theuv-coverage. The spectrum of the point-source was correctlyestimated as  − 1.0, but the extended source acquired a falsesteep-spectrum (α ≈  −0.8). The algorithm was able to reconstruct thecorrect flux and spectrum of the extended source only if themulti-scale basis functions were carefully chosen to match theknown scale size (i.e. a stronger a prioriconstraint). This shows that without additional constraints, it isnot always possible to distinguish large flat-spectrum sourcefrom a slightly less-extended steep spectrum source.

• 2.

  A second imaging run was performed including additional information in the form of short-spacing constraints, added in by retaining a small number of very short-baseline measurements at each frequency (baselines shorter than 25 m). The visibility plots and imaging results with this dataset show that the short-spacing flux estimates were sufficient to bias the solution towards the correct solution (flat-spectrum extended source).

• 3.

  The image residuals are at the same level in both runs, demonstrating that in the absence of additional short-spacing information, both flux models are equally poorly constrained by the data themselves. One way to avoid this problem altogether (but lose some information) is to flag all spatial-frequencies smaller than u_min at ν_max and not attempt to reconstruct any spatial scales larger than what ν_max allows.

7.3. Band-limited signals

Consider a source of radiation where frequency traces different physical structures in the source (as opposed to a fixed structure with each point emitting broad-band radiation). For such sources, emission may be detected in only part of the sampled frequency range and is for all practical purposes a band-limited signal.

Since the MS-MFS algorithm uses a polynomial to model the spectrum of the source (and is not restricted to a power-law spectrum) it should be able to reconstruct such structure as long as it varies smoothly. However, for a band-limited signal, the angular resolution at which the structure can be mapped will be limited to the resolution of the highest frequency at which the signal is detected (and not the highest resolution allowed by the measurements).

EVLA D-configuration data were simulated for the band-limited radiation observed with synchrotron emission from solar prominences where different frequencies probe different depths in the solar atmosphere. The structures are generally arch-like with lower frequencies sampling the top of the loop and higher frequencies sampling the legs. The MS-MFS algorithm was run on these simulated data, using N_t = 5 to fit a 4th-order polynomial to the source spectrum (to accomodate its nearly band-limited nature) and N_s = 3 with scales given by [0, 10, 30] pixels. Iterations were terminated after 200 iterations.

Figure 10 shows a comparison of the true and reconstructed structure at three different frequencies (1.2 GHz (left), 1.8 GHz (middle) and 2.6 GHz (right) ). The top row shows the true structure, and the bottom row is the reconstruction. The 3D structure is mostly recovered, with the largest errors being in the central region where the signal spans the shortest bandwidth. The point source on the right edge of the loop was reconstructed at an angular resolution slightly lower than that of the highest sampled frequency and corresponds to the highest frequency at which this spot is brighter than the background emission. Figure 11 compares the true and reconstructed spectra for three positions in the source (left leg (left), arch of the loop (middle), and right leg (right)). These spectra show the accuracy to which a 4th-order polynomial with MS-MFS was able to reconstruct the structure.

Spectral-lines are an extreme case of a band-limited signal, and the use of MS-MFS imaging is restricted to obtaining a wide-band model of the continuum flux from line-free channels, to be subtracted out of the data before spectral-line strengths are studied.

8. Discussion

The introduction of broad-band receivers into radio interferometry has opened up new opportunities for the study of wide-band continuum emission from a vast range of astrophysical objects. With imaging algorithms that account for the frequency dependence of the incoming radiation as well as of the instrument, we can minimize imaging artifacts, achieve continuum sensitivity and reconstruct spatial and spectral structure at the angular-resolution allowed by the highest observed frequency.

The MS-MFS algorithm models the wide-band sky brightness distribution as a collection of multi-scale flux components whose amplitudes follow a Taylor-polynomial in frequency. The data products are a set of Taylor-coefficient images, which can either be interpreted in terms of continuum intensity, spectral index and curvature, or used to evaluate a spectral cube, or serve as a wide-band model for self-calibration or continuum-subtraction. For wide-field imaging, multiple pointing-centers (mosaicing) and the w-term are accounted for during image-formation, and the effect of the primary beam and its frequency-dependence is approximately corrected in a post-deconvolution step.

For point sources for which standard imaging has a dynamic-range limit of a few thousand, this algorithm achieves dynamic ranges  >10⁵ on test observations with the EVLA, and  >10⁶ on noise-free simulations. For sources with smooth continuum spectra, it is able to reconstruct spectral information at the angular resolution allowed by the combined multi-frequency uv-coverage. However, if the visibility functions of very large-scale emission lie mostly within the central unsampled region of the uv-plane, the algorithm requires a priori knowledge about the spectral structure at those scales. For high SNR data, error-bars on the spectral-index images range from  <0.02 when multi-scale basis functions are chosen appropriately, up to  >0.2 in the extreme case of modeling smooth extended emission with a set of δ-functions. In practice, on-source errors of  △ α ≈ 0.1 have been achieved. For low SNR extended structure (<10), errors due to overfitting can dominate when a high-order polynomial is used.

There are several directions in which wide-band imaging techniques need to be extended and improved. The imaging errors in MS-MFS are currently dominated by the multi-scale aspect of the algorithm, and methods that do adaptive fits (Bhatnagar & Cornwell 2004) and use more a priori information about large-scale spectra may be more appropriate. Methods that adaptively find the optimal number of parameters to operate with would help in error-control. Implementations of algorithms such as MS-CLEAN and MS-MFS are inefficient in memory-use, and other approaches may be required for large image sizes. For wide-field imaging, the time variation of the antenna primary beams must be taken into account during imaging, and wide-band methods combined with algorithms for direction-dependent corrections (Bhatnagar et al. 2008, or Smirnov 2011). For full-Stokes wide-band imaging, where a Taylor-polynomial in frequency is not the most appropriate basis function to model Stokes Q,U,V emission, wide-band imaging with other flux models must be tried.

---

Acknowledgments

The authors would like to thank the National Radio Astronomy Observatory, the New Mexico Institute of Mining and Technology, and the Australia Telescope National Facility for support during the Ph.D. Thesis project that resulted in this algorithm and its implementation within CASA. The authors would like to thank J.A. Eilek in particular, for extremely helpful comments on the presentation of the thesis material that formed the basis of this paper. The authors would also like to thank S.Bhatnagar, K. Golap, R. Nityananda, F. N. Owen, R. J. Sault, and M. A. Voronkov, among others, for useful discussions pertaining to this work and its software implementation. This project used data from the (E)VLA telescope (test observations and project AR664) operated by the National Radio Astronomy Observatory, a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc.

References

Appendix A: Matrix notation and framework

The matrix notation used in this paper is explained here, in the context of an iterative χ² minimization process used to solve a system of linear equations. The basic idea is as follows. Let Ax = b be the system of equations to be solved (measurement equations). The goal is to find a set of parameters x that minimizes χ² = (Ax − b)^†W(Ax − b). Setting grad(χ²) = 0 to minimize χ² leads to a new system of equations called the normal equations  [A^†WA]x =  [A^†]Wb. The matrix on the LHS is called the Hessian  [H] =  [A^†WA]. Iterations begin with an initial guess for the parameters x. These parameters are updated in iteration i as x_i + 1 ← x_i + g [H^-1]A^†W(b − Ax_i), where g controls the step-size. Iterations continue until a convergence criterion is satisfied. The basic iterative framework used in most imaging and deconvolution algorithms in radio interferometry can be described using this matrix notation (Rau et al. 2009; Rau 2010).

A.1. Measurement equations

The measurement equation of an imaging instrument describes its transfer function (the effect of the measurement process on the input signal). For an ideal interferometer (a perfect spatial-frequency filter, with no instrumental gains), it can be written as follows in matrix notation as follows. Let be a pixelated image of the sky and let be a vector of n visibilities. Let S_n × m be a projection operator that describes the instrument’s sampling function (uv-coverage) as a mapping of m discrete spatial frequencies (pixels on a grid) to n visibility samples (usually n > m). Let F_m × m be the Fourier transform operator. Then, the measurement equations become .

A.2. Normal equations

The normal equations are the linear system of equations whose solution gives a weighted least-squares estimate of a set of parameters in a model (χ² minimization). For an ideal interferometer, it is given by where W_n × n is a diagonal matrix of weights and S^† denotes the mapping of measured visibilities onto a spatial-frequency grid. We can write the Normal equations as where the Hessian (matrix on the LHS) is by construction a convolution operator¹⁹ with a shifted version of the point-spread function in each row. The dirty image on the RHS is produced by direct Fourier inversion of weighted visibilities. The normal equations therefore state that the dirty image is the result of the convolution of with . The solution of these normal equations represents a deconvolution.

A.3. Iterative solution

Most existing iterative image reconstruction algorithms in radio interferometry consist of major and minor cycles. Major cycles compute the RHS of the normal equations, and minor cycles perform approximate (implicit) Hessian inversions to calculate updates for the sky model parameters I^m. The first major cycle starts by transforming the observed visibilities into an image as I^dirty =  [F^†S^†W]V^obs, and initializing the sky-model I^m. Minor cycle steps do a deconvolution to calculate updates to the model I^m. After several iterations, this updated model is fed into the next major cycle. This and all subsequent major cycles calculate model visibilities from the current model V^m =  [SF]I^m, calculate residual visibilities V^res = V^obs − V^m, and transform these residuals into images I^res =  [F^†S^†W]V^res. These residual images replace the initial dirty image, and a new set of minor-cycle iterations are done. This process continues until a convergence criterion is reached. Usually, convergence is defined as I^res being noise-like with no signal left to be extracted in the minor-cycle.

Appendix B: MS-MFS as implemented in CASA

The MS-MFS algorithm described in Sect. 2.7 has been implemented and released via the CASA software package (version 3.1 onwards). More recently, it has been implemented in the ASKAPsoft²⁰ package. Algorithm 1 contains a pseudo-code listing.

The main parameters that control the algorithm are

• 1.

  ν₀: a reference frequency chosen near the middle of the sampled frequency range, about which the Taylor expansion is performed,

• 2.

  N_t: the number of coefficients of the Taylor polynomial to solve for, and

• 3.

  N_s and : a set of scale sizes in units of image pixels to use for the multi-scale representation of the image. In order to always allow for the modeling of unresolved sources, the first scale function is forced to be a δ-function.

The data products are N_t Taylor-coefficient images, a spectral index image, and a curvature image (if N_t > 2). This wide-band image model can be used within a standard self-calibration loop.

All Figures