A new approach to multifrequency synthesis
in radio interferometry
^{1}^{1}institutetext: MaxPlanck Institut für Astrophysik (KarlSchwarzschildStr. 1, D85748 Garching, Germany) ^{2}^{2}institutetext: LudwigMaximiliansUniversität München (GeschwisterSchollPlatz 1, D80539 München, Germany)
We present a new approach to multifrequency synthesis in radio astronomy. Using Bayesian inference techniques, the new technique estimates the sky brightness and the spectral index simultaneously. In principle, the bandwidth of a wideband observation can be fully exploited for sensitivity and resolution, currently only limited by higher order effects like spectral curvature. Employing this new approach, we further present a multifrequency extension to the imaging algorithm resolve. In simulations, this new algorithm outperforms current multifrequency imaging techniques like MSMFCLEAN.
1 Introduction
In radio astronomy, multifrequency observations are widely used for many different purposes. Examples include the investigation of spectral line emission, analyzing continuous synchrotron spectra of radio sources, the determination of Faraday rotation in polarization imaging, correlating brightness structures at largely different wavelengths or improving the sensitivity and coverage of an interferometer without introducing more antennas into the array (for a review see e.g. Taylor et al., 1999). Most of these only become possible with observations that span over many frequencies and data analysis methods that can handle this kind of observations.
For the remainder of this work, we focus entirely on multiwavelength radio continuum studies in total intensity. Neither spectral line observations, polarization imaging, nor wavelengths studies beyond pure radio observations will be a direct topic. In Sec. 4, we will comment on possible extensions of the presented work into other domains of multifrequency radio astronomy.
Historically in radio interferometry, the term multifrequency synthesis is mainly used to denote techniques that focus on the direct combination of singleinstrument observations at different frequencies to improve the resolution and sensitivity of a radio interferometer (see Conway et al. (1990) and references therein). For this purpose, a number of methods have been devised, most notably double deconvolution (Conway et al., 1990) and multifrequency CLEAN (Sault & Wieringa, 1994). These methods usually work by Taylorexpanding a spectral model function around a reference frequency. We will comment on these methods in more detail in Sec. 2.1.
Conversely, investigating the spectral behavior of a radio synchrotron source on its own usually is achieved through totally independent standard imaging of the surface brightness in the sky at a number of frequencies. Resolution is kept uniform over all frequencies to correctly recombine the images at individual frequencies. Subsequently, the spectral parameters are determined by fitting a function through all the single frequency images, usually assuming a powerlaw shaped spectral evolution.
Only recently, implementations of MFCLEAN also make it possible to constrain the spectral properties (see Reid & CASA Team, 2010), effectively starting to merge both multifrequency high resolution imaging and spectral analysis into one algorithm. It is this combined notion, how we like to understand and use the term multifrequency synthesis.
Further development of multifrequency imaging techniques is paramount for being able to fully exploit the data from the new generation of radio telescopes, such as the upgraded VLA, LOFAR, the SKA pathfinder missions or ultimately the SKA itself (see e.g. Garrett, 2012). Their unprecedented, broadband frequency coverages of many GHz of possible bandwidths and previously unknown frequency regimes offer many advances in astrophysical and cosmological sciences. But at the same time, they are a challenge for current data analysis methods (see Sec. 2.1). One part of the challenge rises from the fact that current algorithms might not meet the expected sensitivity or fidelity because previously acceptable models and approximations break down with the quality and quantity of data now available. The other side is that the huge quantities in which new data sets come, represent a huge computational challenge on their own, regardless of the algorithms used. It should be emphasized that this study focuses entirely on the first part of the challenge (for this topic, see also Sec. 4).
In this paper, we present a new approach that uses Bayesian statistical inference techniques to combine a simultaneous statistical estimation of the sky brightness and the spectral behavior with the benefits in resolution and sensitivity of classical multifrequency synthesis. Conceptually, this new approach has a number of advantages over standard methods. In principle, the full bandwidth of the observation is used for maximum theoretical sensitivity. No subsequent imaging at all single frequencies is employed and, thereby, less reconstruction artifacts are introduced and no artificial downgrading of resolution is necessary. Our model for the spectral behavior is not approximated through Taylorexpansions around a frequency, and in the moment only is limited by possible higher order effects like spectral curvature. In this way, our approach is conceptually very similar to the new method of Faraday Synthesis (Bell & Enßlin, 2012) in polarization imaging (see Sec. 2.2 for more details). Furthermore, the spatial correlation of the parameters describing the spectral behavior is used to constrain and improve the estimation of spectral properties, and can further be viewed as a new scientific result on its own. Finally, it is possible to approximate the statistical uncertainty of the spectral index estimate in addition to that of the total intensity.
To demonstrate the viability of our new approach, we present a multifrequency extension to the radio extended emission imager resolve (Junklewitz et al., 2013) as an alternative to multiscalemultifrequency CLEAN (Rau & Cornwell, 2011), the standard method for wideband multifrequency synthesis of extended radio sources.
For our derivations, we will often refer to the work presented in Junklewitz et al. (2013), henceforth Res2013.
2 Theory & Algorithm
In this section, we first briefly review the current status of multifrequency synthesis techniques (Sec. 2.1), then develop the new approach (Sec. 2.2) based on information field theory (Enßlin et al., 2009) and the inference framework presented in Res2013. We further present a multifrequencyextension to the imaging algorithm resolve (see Res2013) using the developed framework (Sec. 2.2).
2.1 Multifrequency aperture synthesis and spectral index reconstruction
Aperture synthesis is the technique of connecting an array of telescopes in such a way that we can effectively synthesize a combined instrument with a much larger aperture and therefore resolution (Ryle & Hewish, 1960; Thompson et al., 1986). It can be shown that for observations at a single frequency , and under the assumption of measuring the sky as flat in a plane tangent to the phase center of the observation, such a radio interferometer approximatively takes incomplete samples of the Fourier transformed brightness distribution in the sky (Thompson et al., 1986):
(1) 
For our purposes, working under this assumption suffices for the analysis undertaken in this paper. For a detailed derivation of (1) consider (Thompson et al., 1986). The quantity is called the visibility function. The coordinates and are vector components describing the distance between a pair of antennas in an interferometric array, where this distance is usually referred to as a baseline. They are given in numbers of wavelengths, with and usually parallel to geographic eastwest and northsouth, respectively. The coordinates and are a measure of the angular distance from the phase center along axes parallel to and , respectively. is a sampling function defined by the layout of the interferometric array. It is zero throughout most of the space, apart from where measurements have been made where it is taken to be unity.
The visibility function is what our instrument measures, but we are actually interested in the brightness distribution of the source in the sky. Unfortunately, an inversion of (1) gives us not the true brightness distribution, but its convolution with the inverse Fourier transform of the sampling function, better known as the dirty beam :
(2) 
Here, we have introduced a symbolic Fourier operator with and , the common notation , dirty image, for the simple Fourier inversion of the visibilities, and the symbol to denote a convolution operation.
For imaging at a single frequency, one usually proceeds using a deconvolution algorithm aiming to solve (2) approximatively. The most common algorithm is CLEAN (Högbom, 1974), or one of its many variants (Clark, 1980; Schwab, 1984; Sault & Wieringa, 1994; Cornwell, 2008; Rau & Cornwell, 2011), which basically model the sky brightness to be a collection of delta peak point sources that are iteratively assigned by searching for the peak values in the dirty image. But alternatives exist, especially for imaging extended emission, like the Maximum Entropy Method (Cornwell & Evans, 1985), Adaptive Scale Pixel decomposition (Bhatnagar & Cornwell, 2004) or the recently published Bayesian extended emission imager resolve (Res2013), which will be used in this paper in an expanded multifrequency version (see Secs. 2.2 and 3).
Originally, radio telescopes observed with relatively narrow bandwidths and at only a few different frequencies. The standard approach for imaging the spectral properties of a source since then is to take single frequency observations, or observations averaged over close channels for higher sensitivity, image them separately, and then fit a spectral model to the observations. Since continuous synchrotron emission is known to show a powerlaw spectrum (Rybicki & Lightman, 1985), for many astrophysical purposes such a model is sufficient:
(3) 
Sometimes, higher order spectral deviations are modeled with a second term in the exponential of (3):
(4) 
effectively enhancing the linear function in space (3) to a quadratic polynomial. This model is referred to as spectral curvature.
Simultaneous multifrequency imaging at different frequencies has been introduced to radio astronomy upon the realization that observations at many frequencies, if stacked together appropriately, can improve the sampling in the plane (see Conway et al. (1990) and references therein). That way, the sensitivity of a radio interferometric observation can be enhanced considerably. Because the coordinates are measured in numbers of wavelengths, the same interferometer samples different parts of the Fourier space at different frequencies. An interferometer includes baselines and therefore points, where is the number of antennas. In theory, observing at a number of frequencies enhances the measured baselines to , which is the equivalent effect of introducing roughly extra antennas (Conway et al., 1990).
This approach was first developed by Conway et al. (1990) for midsized fractional bandwidths of around , together with the method of double deconvolution to mitigate spectral errors when using CLEAN for a spectrally combined data set. Later, double deconvolution was developed into the more advanced multifrequency CLEAN (MFCLEAN) (Sault & Wieringa, 1994) and multifrequency multiscale CLEAN (MFMSCLEAN) (Rau & Cornwell, 2011).
All these methods assume the spectral dependence to be a powerlaw, like Eqs. (3) or (4) and propose to approximate it during the CLEANlikedeconvolution process with a Taylorexpansion. To our knowledge, current implementations use a few terms of a Taylorexpansion in or around a reference frequency . In Rau & Cornwell (2011) it is discussed that a direct decomposition into one or two terms of a polynomial in would actually be the most accurate representation^{1}^{1}1Actually, for the spectral model (4), this would be no approximation at all since it simply is a quadratic polynomial in space., but for current implementations this is discarded by the authors because of numerical instabilities.
During the CLEAN deconvolution process, the iteratively updated sky model is used to also update the coefficients of the spectral Taylor expansion. In this way, spectral index or even spectral curvature can be constrained by these coefficients. Care must be taken, since this expansion is usually stopped after a few terms and might not be valid over large bandwidths. Because of this, the only implementation of MFMS CLEAN known to the authors to date, within the radio astronomical software package CASA (Reid & CASA Team, 2010), is considered to be experimental and on a sharedrisk basis with regard to the spectral reconstructions with a higher number of terms.
For modern wideband data sets from the new generation of instruments, spanning several GHz of bandwidths, this means that probably even more emphasis must be put onto the handling of the spectral effects, either by invoking higher order terms in the expansions (see Sault & Wieringa, 1994; Rau & Cornwell, 2011) or by shifting to a different approach, where the sky brightness and its spectral properties are fully considered simultaneously, and estimated to fit the entire data using some global minimization function. This last approach was actually mentioned by Conway et al. (1990), but found to be unnecessarily complicated and computationally expensive for the typically modest fractional bandwidths at the time.
We will now consider the latter approach and present a statistical solution.
2.2 A multifrequency extension to the algorithm resolve
In the course of this section, we often refer to the detailed derivations layed out in Sec. 2 and App. 1 of Res2013 that form the basis from which we derive our multifrequency algorithm.
(5) 
where again and , the term introduces measurement noise, the data has been introduced, which is basically the visibility function with measurement noise, and the spectral dependence of the source is kept general for the moment.
In order to simplify notation and to identify (5) as a general inverse inference problem as analyzed in Res2013, we henceforth drop all explicit dependence on and and combine all known instrumental effects into a response function , leading to
(6) 
There are many instrumental effects beyond the sampling in the plane given by like an antenna sensitivity pattern or an direction dependent, variable sampling. For this work, we stay with the basic definition and refer the reader to Res2013 concerning the possibility of including other effects within this framework. We just emphasize that, without loss of generality, most of these effects can be in principle included^{2}^{2}2In principle this approach can be extended into a full RIME (Radio Interferometer Measurement Equation), considering all Stokes parameters and instrumental gains (Smirnov, 2011a, b).. We also assume the instrument to be fully calibrated, and thus the response to be known.
In accordance with standard radio interferometric literature (Thompson et al., 1986), we assume Gaussian noise statistics, mainly induced by the antenna electronics and independent between measurements at different frequencies and time steps of the observation (Thompson et al., 1986). Henceforth, the noise will be assumed to be drawn from a multivariate, zero mean Gaussian distribution of dimension :
(7) 
Solving (6) exactly for is not possible, since all information is lost on the Fourier modes not sampled by . This is just a different way of stating that a direct Fourier inversion of (1) yields the dirty image and not an exact representation of the sky brightness.
Instead, the aim is to find a statistical estimate for the most probable sky brightness signal, given all observational and noise constraints. In Res2013, it was shown in detail how this can be done in a Bayesian statistical framework for a radio astronomical data model like (1). We briefly repeat the main points, and otherwise refer the reader to Res2013.
To find an optimal statistical estimate for the sky brightness signal , we regard it as a random field with certain a priori statistical properties expressed in the prior distribution , but fully constrained by the data through the statistics of the likelihood distribution . The likelihood distribution summarizes how the data are obtained with a measurement of the true sky brightness signal, and for our problem can be expressed as
(8) 
which is a Gaussian over with the covariance structure of the uncorrelated noise .
Prior and likelihood statistics can be combined into the posterior distribution that holds the important information of how much the sky brightness signal is statistically constrained by the data. From there, an estimate for the signal can be obtained by calculating a suitable statistic of the posterior, most prominently its mean or its mode, corresponding to the minimization of different error norm measures between the signal and its estimate (see Res2013; or Jaynes (2003), Caticha (2008) and Lemm (1999), Enßlin et al. (2009) for a comprehensive review of Bayesian statistics or inference on fields respectively).
The exact choice of an appropriate inference algorithm at this stage largely depends on the complexity of the problem (i.e. the posterior). For the problem at hand, since the likelihood is already known (8), this comes down to the question of the prior statistics of .
The most general way conceivable would be not to explicitly model the spectral dependence of at all, as for instance in (3) or (4). Instead, can be interpreted as a three dimensional continuous field and should be inferred as a whole from the entire data. In general, let us assume not only pointlike but also extended emission in the sky, and some kind of extended structure in spectral space as well. In such a setting, could be set a priori as a statistical field with an unknown and probably nonisotropic crosscorrelation structure in the combined sky and spectral space. Such a complex, unknown crosscorrelation structure at the one hand complicates the problem enormously, but could also be used to guide the reconstruction, if correctly estimated with the field itself.
A number of statistical methods have already been developed to handle simpler problems of signal reconstructions with unknown but isotropic correlation structure, many of them solving the problem for Gaussian fields using information field theory (Enßlin & Frommert, 2011; Enßlin & Weig, 2010) or GibbsSampling MonteCarlo methods (Sutter et al., 2012; Karakci et al., 2013), or recently also for lognormal fields (Enßlin & Frommert, 2011; Enßlin & Weig, 2010; Oppermann et al., 2013; Greiner, 2013; Selig & Enßlin, 2013). Most notably this method was also used to create resolve (Res2013). A full combined spatial and spectral reconstrcution as outlined above, would require substantial further development and is outside of the scope of this work.
We can choose a more direct strategy, still residing within our approach of statistical inference, but fix a spectral model and infer instead the spectral index (or curvature) as a field on its own. For the rest of the paper, we will choose the model (3) for the simplicity of the approach, and for a functional similarity with the algorithm resolve that makes it very natural to include into a combined method.
resolve works under the assumption that the extended surface brightness at a single frequency is a priori assumed as a random field drawn from lognormal statistics (Res2013). For our multifrequency problem, this basically turns (6) into
(9) 
where is a Gaussian random field (such that the logarithm of is a Gaussian random field again), and is a constant to e.g. normalize the system to the right units. Although, the frequency dependence of was not explicitly shown in Res2013, the derivation of resolve implicitly assumed the algorithm to work for a single frequency in the way presented here. resolve assumes the spatial correlation of an extended source in the sky to be reflected by the covariance of the Gaussian random field , which is unknown a priori, and thus estimated from the data itself together with the sky brightness. The covariance of a Gaussian field is equivalent to its twopoint correlation function and is handled as such by resolve in form of the power spectrum , which is the Fourier transformation of the correlation function^{3}^{3}3We actually assume the spatial correlation to be a priori rotationally and translationally invariant, and thus the power spectrum to be diagonal . However, this does not imply that the correlation structure must be invariant under any transformation a posteriori as well. A more detailed discussion can be found in Res2013.. A deeper analysis of this can be found in Res2013.
We now make the central assumption that the spectral index can be modeled a priori as a Gaussian random field with its own spatial correlation structure in the sky. At least for an extended source, we have every reason to assume that the spectral index should be a field with spatial extension itself. Observational constraints strongly imply that typical extended radio structures show as well extended and smooth spectral index structures, for instance radio halos and relics of galaxy clusters (Feretti et al., 2012), radio galaxy lobes (Kassim et al., 2005), or supernova remnants (Green, 2001).
In Res2013, we argued extensively that a Gaussian random field would be the ideal choice for a signal prior of an extended field with a priori unknown correlation structure, as long as the field is not assumed to vary strongly on orders of magnitude and not necessarily needs to be positive definite. Both constraints apply very well to known spectral index maps, where variations usually do not reach even one order of magnitude, and nothing prevents the spectral index in principle to change the sign. For details of these arguments see Res2013, we now proceed under this assumption.
If we rewrite (3) only slightly,
(10) 
it reveals that, if has a Gaussian prior, (10) naturally turns into a model for a lognormal prior, only different in shape (and more complicated) from (9) because of the term .
It is important to note that we have not specified yet. Thus, at this point, our inference approach to multifrequency synthesis is in principle compatible with any method that reconstructs and deconvolves the surface brightness at a single reference frequency . At least as long as it seems consistent with the source of interest to assume that the spectral index is an extended and spatially correlated field. In an extreme case, like single, unresolved point sources, this method probably will not yield optimal results.
For this paper, we take the choice to combine (9) and (10) into one single method, where we assume our double lognormal measurement model to be
(11) 
with again a constant, from now on set to one w.l.o.g., and the signal fields and having Gaussian prior distributions and :
(12)  
(13) 
We now write down a posterior distribution for each signal field, while the other field (and its covariance) is assumed to be known, held constant and regarded as part of two distinct versions of the response operator , henceforth called and (where the symbol denotes where the field needs to be inserted that the operator acts on):
(14)  
(15) 
As in Res2013, it is not possible to calculate the posterior mean for either of the two signal fields or without invoking complicated or expensive perturbation or sampling methods (see Res2013). We therefore continue to use the procedure already presented there, and calculate the posterior maximum to estimate both signals
(16) 
In signal inference, this procedure is called Maximum A Posteriori (MAP) (Jaynes, 2003). The resulting fixpoint equations from Eqs. (16) need to be solved numerically using a nonlinear optimization scheme. For this, we resort to the same implementations as in Res2013 (see App. 2 therein).
With this choice, we basically extend resolve to a multifrequency algorithm by integrating a second complete resolve step for the spectral index into the method, and iterating between the statistical estimation of with its covariance , and with its covariance . As outlined above, we always hold one of the fields (and their respective covariances) as constant and regard them as part of the response during this process.
It should be emphasized again that the resolve step for the spectral index could in principle be combined with any other method to reconstruct at a single frequency .
3 Tests
We have integrated multifrequency capability into the existing implementation of resolve and tested the algorithm using simulated data^{4}^{4}4To get access to the preliminary code prior to its envisaged public release, please contact or .. The code is written in Python using the signal inference library NIFTy (Selig et al., 2013), for details of the implementation we refer the reader to App. 2 of Res2013.
As in Res2013, we constructed simulated observations with the tool makems^{5}^{5}5See http://www.lofar.org/wiki/lib/exe/fetch.php
?media=software:makems.pdf. using a realistic coverage from a VLA observation in its AConfiguration. The VLA samples the plane nonuniformly at irregular intervals, and the response includes thereby a convolutional gridding and degridding operator using a KaiserBessel kernel (for details see App. 2 in Res2013). We simulated observations over a range of 2 GHz, with 20 separate frequency channels. The observations are short snapshots of approximatively 20 minutes per frequency, with a total of 42 120 visibility measurements at each frequency channel (see Fig. 1). This setting leads to an especially sparse sampling of the plane at a single frequency, but to a much better coverage for the combined multifrequency data (see Fig. 1).
For all tests, the signal is drawn from a Gaussian distribution, finally entering the formalism as an exact lognormal field^{6}^{6}6In Res2013 it was demonstrated that resolve also works beyond that on realistic signals drawn from real CLEAN maps.
In order to use a spectral index signal with some correlation to the brightness signal, we model the spatial source dependence of the spectral index in an ad hoc fashion.
For the spectral index signal , we use a sum of two Gaussian fields (see Sec. 3.1), one being independently drawn, while the second is the surface brightness signal itself, downweighted with a suitable factor. This is done in order to introduce a spectral index signal with some correlation to the surface brightness signal, at least in an adhoc fashion. Typically, extended sources in radio astrophysics show a spatial crosscorrelation between both, which is a result of different physical processes underlying the structure and formation of these sources. A notable example are radio halos and relics in galaxy clusters (Feretti et al., 2012). Of course, our model is only an adhoc approximation for the proof of concept undertaken in this work. We also emphasize that such a crosscorrelation between and is not exploited by resolve (see Sec. 2.2 for an outlook).
The complex, Gaussian input noise variance in space is equal for all visibilities and frequencies. As with singlefrequency resolve, the algorithm does not require equal noise variances and can in principle handle varying variances. The noise variance was set to a low^{7}^{7}7In comparison to the signal strength in Fourier space, the chosen value ensures a high signal to noise ratio. The unit Jy was used here for convenience. Effectively, it stands for whatever units the simulated signal is interpreted to be given in. value of .
In continuation with Res2013, we use a relative  norm measure of the difference in the signal to the reconstruction to measure the accuracy of the estimate of both brightness and spectral reconstructions:
(17)  
(18) 
where the sums are taken over all pixels of the reconstruction. For the motivation behind this choice see Sec. 3 in Res2013.
In Sec. 3.1 we focus exclusively on the reconstruction of the two signals and . Then, in Sec. 3.2, we compare the results from multifrequency resolve with standard imaging procedures. The reconstruction of the signal power spectra is discussed separately in Sec. 3.3.
3.1 Main test results
We start by showing the reconstruction of the presented simulated observations using multifrequency resolve. In Fig. 3, a surface brightness and a spectral index signal are shown, together with the respective reconstructions obtained with resolve and absolute difference maps of both signals to their reconstructions. We show the spectral index maps in full, and overlayed with a mask that focuses on the part of the observed field that contains the brightest part of the surface brightness signal. Later, in Fig. 4, we show the reconstruction for different masks (along with a comparison to other methods, see Sec. 3.2). We choose the three different masks mainly qualitatively by visual comparison to show only the parts of the sky with a surface brightness signal above a certain threshold. The actual threshold values for the surface brightness were Jy/px, Jy/px, and Jy/px. In Fig. 3, the second conservative mask is used. The error measures are for the surface brightness, and for the spectral index.
It can be seen that resolve recovers very accurately the original surface brightness. For this particular choice of relative low noise, the algorithm succeeds in reconstructing even small scale features of the signal. The effects of the instrumental point spread function are successfully deconvolved. These findings are in perfect agreement with the results presented on singlefrequency resolve in Res2013.
The general structure of the spectral index is well reconstructed, even for outer regions where the brightness signal is very weak. Overall, small scale features are much better recovered in the inner regions, where the main brightness sources are located, as can be seen by comparison between the full and the masked images in Fig. 3. It is expected that the quality of the spectral index reconstruction depends on the strength of the observed surface brightness, as is illustrated by the poor performance of standard methods in recovering any structures outside the strong source regions, presented in the next section (see Fig. 4). In most real applications, the outer parts of the observed fields should therefore usually be not a focus of the investigation. We thus choose the relatively conservative mask in Fig. 3 to highlight this important part of the reconstruction. Nevertheless, it should be noted that – at least for this low noise example with relatively high sensitivity due to the broad bandwidth – resolve is able to reliably extrapolate its estimation also into the weaker regions around the main sources.
In addition, with resolve, an uncertainty of the spectral index reconstruction can be estimated. As for the single frequency reconstructions presented in Res2013, the second derivative of the posterior is used to approximate its covariance (for details, see App. A and Res2013). This way, the full estimate for the spectral index signal becomes
(19) 
For the given set of simulated data, roughly of the original signal values within the unmasked regions in Fig. 3 lie within a interval , but only roughly lie within a interval for the full spectral index image. Due to the nonlinear nature of the inference problem (16), it is expected that the uncertainty estimate is not exactly what would be expected from a pure Gaussian covariance (i.e. in the interval), the MAP estimate is not guaranteed to lie very close to the real posterior covariance (see Res2013 for this problem). It is no surprise that the estimate worsens for all the outer regions with only very weak brightness structures present, which explains the poor performance on the whole spectral index map. Furthermore, as discussed in Res2013, the calculation of an uncertainty estimate is computationally costly. Due to a lack of accessible computer power only for testing purposes, we stopped the calculations at some point and smoothed the outcome. This should not be a great problem, since the uncertainty is expected to be smooth and we are mainly interested in a proof of concept at this point.
3.2 Comparison to standard methods
In this section, we compare the results of multifrequency resolve with two standard methods: A straight forward powerlaw fit for (3) using singlefrequency CLEAN images, and a MFMSCLEAN reconstruction for the surface brightness and the spectral index. The reconstructions were performed on the same simulated observation as in Sec. 3.1. All results were obtained using the radioastronomical software package CASA (Reid & CASA Team, 2010). The MFMSCLEAN reconstructions were obtained using 1500 iterations, a small gain factor of , uniform weighting and ten different multiscales ranging from a single pixel to moderately large structures. We further used two terms for the Taylor expansion in (see Sec. 2.1) and only half the available bandwidth, since a larger frequency range seems to lie outside the convergence radius of the Taylorexpansion (see Sec. 2.1 used in the implementation of CASA.
In Fig. 4, a comparison is shown between resolve, a powerlaw fit, and MFCLEAN results for the spectral index and differently strong masks. The full field is not shown, because neither the powerlaw fit nor the CLEAN reconstruction could recover any structures in the remaining regions. The error measures are listed in Tab. 1. It can be seen that for this example resolve overall outperforms the other methods. The advantage of resolve is more pronounced in the outer regions, where the surface brightness signal is weak, but resolve still gives the best result even for the most central parts of the reconstruction (see also the full reconstruction in Fig. 3).
Fig. 5 shows a comparison of multifrequency surface brightness reconstructions with resolve and MSMFCLEAN for the simulated observation of Sec. 3.1. For this simulation, resolve is more successful in reconstructing the overall structure, and especially recovers more of the small scales. It should be noted that this particular CLEAN image was achieved using the full bandwidth and coverage (other than for the spectral index images, as stated at the beginning of the section).
Algorithm \Mask  liberal mask  medium mask  conservative mask 

resolve  0.34  0.33  0.32 
powerlaw fit  132.79  0.94  0.48 
MSMF CLEAN  3.07  1.89  1.84 
3.3 Power spectrum reconstructions
The power spectrum of the spectral index is reconstructed together with the field itself. It represents the spatial correlation of the spectral index over the observed sky. As explained in Sec. 2.2, we expect the typical spectral index structure of an extended (i.e. spatially correlated) radio source to be spatially correlated to itself. For all the details on power spectrum reconstructions we refer the reader to Res2013, since almost all derivations for standard singlefrequency resolve power spectrum reconstructions are valid as well for the spectral index.
In this section we only discuss the spectral index power spectrum, since the reconstructions for the brightness power spectra are identical to the ones already presented in Res2013. We simply note that, in principle, a multifrequency reconstruction recovers structure on smaller scales, and thus, it is expected that multifrequency resolve should be able to map spatial power spectra up to higher Fourier modes (i.e smaller correlation structures).
A typical result from multifrequency resolve is illustrated in Fig. 2. It shows the original spectral index power spectrum of the Gaussian signal field, used in the simulated observations of Sec. 3.1, and its reconstruction that belongs to the same iteration step as the presented signal reconstructions earlier. The power spectrum is reconstructed relatively well, no power is lost on high modes (i.e. small scales), which is just a consequence of the fact that the simulated observation was conducted with low noise.
As for the brightness reconstructions in Res2013, we emphasize that the power spectrum should not only be viewed as a byproduct of the algorithm in order to accurately estimate the spectral index signal. For instance, some astronomical objects show very distinct spectral index structures and can even be classified after this criterion. A prominent example might be radio halos and relics of galaxy clusters, both of which typically show steep spectral indices that evolve spatially roughly like the source itself (Feretti et al., 2012). Measuring the spatial power spectra of these objects might lead to a more exact and quantitative classification scheme. Another application might lie in the investigation of different physical processes within a single source that lead to spectrally very different regions. Estimating the spectral correlation structure over such regions offers a new way of quantitative analysis of the interplay of these processes.
4 Conclusion
We presented a multifrequency extension to the imaging algorithm resolve (Junklewitz et al., 2013). The combined algorithm is optimal for multifrequency imaging of extended radio sources. It simultaneously estimates the surface brightness at a reference frequency, and the spectral index of the source. Within the assumption of a spectral index model, no further expansions or parameterdependent modeling is used in the reconstruction. Multifrequency resolve is thus capable of exploiting the full bandwidth of a modern radio observation for maximum sensitivity and resolution, only limited by higher order spectral effects like spectral curvature.
Multifrequency resolve has been tested successfully using simulated observations of the VLA in its Aconfiguration. For the presented tests, the algorithm can outperform standard imaging methods in both surface brightness and spectral index estimations.
The algorithm uses a Gaussian prior and an effective lognormal model for the spectral index. This approach is not necessarily restricted to a combination with singlefrequency resolve, and can in principle be combined with any other imaging or deconvolution method for the surface brightness reconstruction.
For the sake of feasibility, many details of radio interferometric observations have been left out of the analysis. The response might realistically contain a number of additional effects, for instance an instrumental primary beam or widefield and direction dependent effects. We refer the reader to Junklewitz et al. (2013), where many of these problems already have been discussed in the outlook.
In its current form, the presented algorithm has relatively high computational costs and numerical demands (for an analysis of algorithmic efficiency, see Junklewitz et al., 2013). In general, this is true for most Bayesian statistical inference algorithms. This could pose an obstacle for daytoday applicability, since especially modern broad band data sets tend to be very large, and therefore are already a numerically challenge on their own. Since this study was centered on fundamental algorithmic development, numerical efficiency was not a focus. Future work might be needed to obtain a more efficient implementation of the algorithm.
For the future, it seems to be desirable to refrain completely from an explicit spectral model, and try to infer a full, nonparametric, three dimensional spectral intensity . Possibly, such a development could benefit greatly from reconstructing a full crosscorrelation structure between the sky and spectral space. We leave this more complete approach for a future publication.
References
 Bell & Enßlin (2012) Bell, M. R. & Enßlin, T. A. 2012, A&A, 540, A80
 Bhatnagar & Cornwell (2004) Bhatnagar, S. & Cornwell, T. J. 2004, A&A, 426, 747
 Caticha (2008) Caticha, A. 2008, ArXiv eprints
 Clark (1980) Clark, B. G. 1980, A&A, 89, 377
 Conway et al. (1990) Conway, J. E., Cornwell, T. J., & Wilkinson, P. N. 1990, MNRAS, 246, 490
 Cornwell (2008) Cornwell, T. J. 2008, IEEE Journal of Selected Topics in Signal Processing, 2, 793
 Cornwell & Evans (1985) Cornwell, T. J. & Evans, K. F. 1985, A&A, 143, 77
 Enßlin & Frommert (2011) Enßlin, T. A. & Frommert, M. 2011, Phys. Rev. D, 83, 105014
 Enßlin et al. (2009) Enßlin, T. A., Frommert, M., & Kitaura, F. S. 2009, Phys. Rev. D, 80, 105005
 Enßlin & Weig (2010) Enßlin, T. A. & Weig, C. 2010, Phys. Rev. E, 82, 051112
 Feretti et al. (2012) Feretti, L., Giovannini, G., Govoni, F., & Murgia, M. 2012, A&A Rev., 20, 54
 Garrett (2012) Garrett, M. A. 2012, in From Antikythera to the Square Kilometre Array: Lessons from the Ancients
 Green (2001) Green, D. A. 2001, in American Institute of Physics Conference Series, Vol. 558, American Institute of Physics Conference Series, ed. F. A. Aharonian & H. J. Völk, 59–70
 Greiner (2013) Greiner, M. 2013, The Galactic Free Electron Density: A Bayesian Reconstruction, Master Thesis
 Högbom (1974) Högbom, J. A. 1974, A&AS, 15, 417
 Jaynes (2003) Jaynes, E. T. 2003, Probability Theory: The Logic of Science
 Junklewitz et al. (2013) Junklewitz, H., Bell, M. R., Selig, M., & Enßlin, T. A. 2013, ArXiv eprints (Res2013)
 Karakci et al. (2013) Karakci, A., Sutter, P. M., Zhang, L., et al. 2013, ApJS, 204, 10
 Kassim et al. (2005) Kassim, N. E., Lazio, T. J. W., Lane, W. M., et al. 2005, in XRay and Radio Connections, ed. L. O. Sjouwerman & K. K. Dyer
 Lemm (1999) Lemm, J. C. 1999, ArXiv Physics eprints
 Oppermann et al. (2013) Oppermann, N., Selig, M., Bell, M. R., & Enßlin, T. A. 2013, Phys. Rev. E, 87, 032136
 Rau & Cornwell (2011) Rau, U. & Cornwell, T. J. 2011, A&A, 532, A71
 Reid & CASA Team (2010) Reid, R. I. & CASA Team. 2010, in Bulletin of the American Astronomical Society, Vol. 42, American Astronomical Society Meeting Abstracts #215, 479.04
 Rybicki & Lightman (1985) Rybicki, G. B. & Lightman, A. P. 1985, Radiative processes in astrophysics.
 Ryle & Hewish (1960) Ryle, M. & Hewish, A. 1960, MNRAS, 120, 220
 Sault & Wieringa (1994) Sault, R. J. & Wieringa, M. H. 1994, A&AS, 108, 585
 Schwab (1984) Schwab, F. R. 1984, AJ, 89, 1076
 Selig et al. (2013) Selig, M., Bell, M. R., Junklewitz, H., et al. 2013, A&A, 554, A26
 Selig & Enßlin (2013) Selig, M. & Enßlin, T. 2013, ArXiv eprints
 Smirnov (2011a) Smirnov, O. M. 2011a, A&A, 527, A106
 Smirnov (2011b) Smirnov, O. M. 2011b, A&A, 527, A107
 Sutter et al. (2012) Sutter, P. M., Wandelt, B. D., & Malu, S. S. 2012, ApJS, 202, 9
 Taylor et al. (1999) Taylor, G. B., Carilli, C. L., & Perley, R. A., eds. 1999, Astronomical Society of the Pacific Conference Series, Vol. 180, Synthesis Imaging in Radio Astronomy II
 Thompson et al. (1986) Thompson, A. R., Moran, J. M., & Swenson, G. W. 1986, Interferometry and synthesis in radio astronomy
Appendix A Details of the algorithm
In this appendix, we repeat some of the basic derivations from Res2013, and derive the details of how reconstructing the spectral index using resolve differs from the standard procedure presented there. In principle, many of the original derivations are still valid for multifrequency resolve, estimating and . The derivations for power spectrum reconstructions and uncertainty calculations stay especially close to the equations already presented in Res2013, and we do not repeat them explicitly here.
a.1 Reconstruction of the sky brightness signal field
The estimation of the sky brightness at a reference frequency can mostly be conducted with the standard singlefrequency resolve.
The only difference is the more complex response operator . Dropping the explicit signal index for a moment and writing the specific measurement model (11) for the signal with explicit operator and field indices
(20) 
reveals that the operator actually spans the singlefrequency Fourier transform over all observed frequencies into a large, manyfrequency data space^{8}^{8}8As a reminder: In (20) we use an implicit convention to sum or taking the integral over repeated discrete or continuous indices, see Res2013 for details.. Conversely, the adjoint response now includes a sum over all frequencies to collapse everything back into the two dimensional sky at reference frequency. Since both operations are used in resolve, it is by this procedure that the multifrequency resolve effectively uses information at all frequencies to constrain the estimation of . Of course, the quality of this reconstruction depends on the accuracy of the current estimation for , used in the reconstruction step for to define .
a.2 Reconstruction of the spectral index
The lognormal model for the spectral index, , only differs in the term from the model used for the sky brightness signal. This slightly complicates the calculations for the signal estimation with MAP and the power spectrum reconstruction and eventual uncertainty calculation. For both, derivatives of the posterior (15) with respect to are needed (see Res2013), and this will add extra terms into the equations. The results of this are summarized in Sec. A.3.
As for the sky brightness signal (see App. A.1), the response operator for the reconstructions is more complex. Effectively, the term in the basic lognormal model would prevent us from just defining our signal space to be the twodimensional sky since acts on . In the actual implementation of the multifrequency resolve algorithm, we circumvent this problem by assuming that acts on and all other terms, including the exponential operation of the lognormal model, are part of the operator^{9}^{9}9This actually renders a nonlinear operator.. In this way, and can be understood to act in the same way as their respective counterparts for do (i.e. they also span up and collapse into the full frequency data space).
a.3 Combined algorithm
Using our findings in Res2013 and in the previous subsections, multifrequency resolve comes down to solving iteratively two only slightly different, subsequent sets of equations:
(21)  
(22)  
(23)  
(24)  
(25)  
(26) 
A rigorous derivation for all equations can be found in Res2013. The two sets of equations only differ in form by the  terms that show up in the spectral index reconstruction because of the derivatives used in order to calculate the MAP estimate. The quantities , and , are defined as
(27)  
(28)  
(29)  
(30) 
, or are projection operators onto a band of Fourier modes denoted by the index , while are parameters to model the unknown power spectrum into a number of such bands , or (see Res2013 for details). The quantities , and are parameters of a power spectrum prior for the signal or the spectral index, and is an operator, which enforces a smooth solution of the power spectrum .
Eqs. (21) and (24) are the fix point equations that need to be solved numerically to find a MAP signal estimate or for the current iteration. The second equations (22) and (25) result from calculating the second derivative of the respective posteriors for the signal estimates or , their inverses serve as an approximation to the signal uncertainty or at each iteration step. The last equations (23) and (26) represent an estimate for the signal power spectra (and therefore their autocorrelation functions), using the signal uncertainties or to correct for missing signal power in the current estimates or . The iteration is stopped after a suitable convergence criterion is met (see App. 2 in Res2013).
























