Issue 
A&A
Volume 553, May 2013



Article Number  A105  
Number of page(s)  13  
Section  Astronomical instrumentation  
DOI  https://doi.org/10.1051/00046361/201220882  
Published online  17 May 2013 
Applying full polarization AProjection to very wide field of view instruments: An imager for LOFAR
^{1}
GEPI, Observatoire de Paris, CNRS, Université Paris Diderot,
5 place Jules Janssen,
92190
Meudon,
France
email: cyril.tasse@obspm.fr
^{2}
Department of Physics & Electronics, Rhodes
University, PO Box
94, 6140
Grahamstown, South
Africa
^{3}
SKA South Africa, 3rd Floor, The Park, Park Road,
7405
Pinelands, South
Africa
^{4}
Sterrewacht Leiden, PO Box 9513, 2300 RA
Leiden, The
Netherlands
^{5}
Netherlands Institute for Radio Astronomy (ASTRON),
Postbus 2,
7990 AA
Dwingeloo, The
Netherlands
^{6}
National Radio Astronomy Observatory, Socorro, NM
87801,
USA
Received:
11
December
2012
Accepted:
9
February
2013
The required high sensitivities and large fields of view of the new generation of radio interferometers impose high dynamic ranges, e.g., ~1:10^{6} to 1:10^{8} for the Square Kilometre Array (SKA). The main problem for achieving these high ranges is the calibration and correction of direction dependent effects (DDE) that can affect the electromagnetic field (antenna beams, ionosphere, Faraday rotation, etc.). It has already been shown that the AProjection is a fast and accurate algorithm that can potentially correct for any given DDE in the imaging step. With its very wide field of view, low operating frequency (~30–250 MHz), long baselines, and complex stationdependent beam patterns, the LOw Frequency ARray (LOFAR) is certainly the most complex SKA pathfinder instrument. In this paper we present a few implementations of the AProjection in LOFAR which can deal nondiagonal Mueller matrices. The algorithm is designed to correct for all DDE, including individual antennas, projection of the dipoles on the sky, beam forming, and ionospheric effects. We describe a few important algorithmic optimizations related to LOFAR’s architecture that allowed us to build a fast imager. Based on simulated datasets we show that AProjection can dramatically improve the dynamic range for both phased array beams and ionospheric effects. However, certain problems associated with the calibration of DDE remain (especially ionospheric effects), and the effect of the algorithm on real LOFAR survey data still needs to be demonstrated. We will be able to use this algorithm to construct the deepest extragalactic surveys, comprising hundreds of days of integration.
Key words: instrumentation: interferometers / techniques: interferometric / techniques: image processing
© ESO, 2013
1. Introduction: LOFAR and directiondependent effects
With the building or development of many large radio telescopes (LOFAR, EVLA, ASKAP, MeerKAT, MWA, SKA, eMerlin), radio astronomy is undergoing a period of rapid development. New issues arise with the development of new types of interferometers, and certain approximations applicable to imaging with the older generation of instruments are not valid anymore. Specifically, they have wide fields of view and will be seriously affected by direction dependent effects (DDE). Dealing with DDE represents serious challenges, in the theoretical, numerical, and technical aspects of calibration and imaging (see Bhatnagar 2009, for a detailed review).
Fig. 1 LOFAR stations are phased arrays characterized by a large field of view. X/Y polarimetric measurements made with it are therefore nontrivial compared to a radio telescope using dishes: we have to take into account the projection of the dipoles that are generally nonorthogonal on the sky and the fact that this angle varies across the field of view. This figure shows the Mueller matrix corresponding to baseline (01) in a given time and frequency slot, normalized by the Mueller matrix at the phase center (see text). Each pixel in the plot (i,j) shows the amplitude of the (i,j) Mueller matrix element in a certain direction s using a logarithmic scale. Even in this normalized version, the offdiagonal Mueller terms are as high as 10% and cannot be neglected. 
This is particularly true for the LOw Frequency ARray (LOFAR), an instrument that observes in a mostly unexplored frequency range (ν ≲ 240 MHz), and will be one of the largest radio telescopes ever built in terms of collecting area. LOFAR’s design is based on a combination of phased array and interferometer concepts (see de Vos et al. 2009, for a description of the LOFAR system). It presently consists of 40 stations in the Netherlands, and 8 international stations (5 in Germany, and 1 each in France, Sweden, and the United Kingdom). Each high band antenna station (110–240 MHz, HBA hereafter) comprises 24 to 96 tiles of 4 × 4 coherently summed antennas, while in a low band antenna station (10–80 MHz, LBA) they are clustered in groups of 96 elements. At the station level, the signals from the individual antennas or tiles (for LBA and HBA, respectively) are phased and summed by the beamformer. This step amounts to forming a virtual antenna pointing towards the targeted field location. The data are transported from the various stations of the LOFAR array to the correlator. The whole process and the pipeline architecture are described in more detail in Heald et al. (2010). LOFAR is affected by many complex baselinetimefrequency4 dependent DDE, consisting mainly of effects in the antenna/station beams and the ionosphere. We currently have models of both the highband and lowband station beams (HBA and LBA, respectively).
As shown in Bhatnagar et al. (2008), AProjection allows sky images to be estimated, taking all possible complicated effects associated with DDE into account (see also Bernardi et al. 2011; Mitchell et al. 2012; Sullivan et al. 2012, in the context of the Murchison Widefield Array and forward modeling). However, unlike dishbased interferometers, whose fields of view are typically 0.5 degree and where the beam shape and polarization angle are affected by pointing errors and rotated on the sky by the parallactic angle (depending on the dish mount), LOFAR is based on phased arrays that have very wide fields of view (up to ~12 degrees), and nontrivial and quickly varying beams, thereby driving complicated polarization effects. Technically speaking, very wide fieldofview instruments that aim to reach a high dynamic range have to deal with baselinedependent nondiagonal Mueller matrices (see Sect. 2 for a detailed discussion). For the VLA, due to its relatively small field of view, it was sufficient for AProjection to take only the diagonal terms of the Mueller matrices into account to demonstrate corrections for instrumental polarization. This is not possible for LOFAR, however, due to its heavily nondiagonal baselineassociated Mueller matrices, where all 4 × 4 Mueller terms have to be taken into account (see Sect. 3 for a detailed discussion).
We show in this paper that the scheme described in Bhatnagar et al. (2008) can indeed deal with the heavily nondiagonal Mueller matrices associated with the very wide field of view of phased arrays. Our imaging algorithm can take any model or calibration solution or ionosphere phase screen as input. In Sect. 2 we describe the issues related with the usage of phased arrays in interferometers, and focus on LOFARrelated issues, i.e., polarization aspects and baseline dependence of DDE. We describe a few important algorithmic optimizations related to LOFAR’s architecture, which allow us to build a fast imager^{1}. In Sect. 3 we describe the AProjection algorithm first presented in Bhatnagar et al. (2008), and the various implementations and optimizations we have found that make it reasonably fast in the case of LOFAR. We present the results in Sect. 4 and show that beam and ionosphere corrections can both be performed at high accuracy. We summarize and discuss the next developments in Sect. 5.
Fig. 2 Correction that can be applied to the image before the minor cycle. This is a firstorder correction for the complicated phased array beam, which depends on time, frequency, and baseline. 
2. Polarization effects associated with very wide fields of view interferometers
In this section we describe the polarizations effects associated with the complex structure of DDE that are inherent in the use of phased arrays with very wide fields of view and nondiagonal, baselinedependent Mueller matrices (or nonunitary Jones matrices in the case of interferometers having similar antennas, see Fig. 1 and the discussion in Sect. 2.1).
For an instrument with long baselines, large fractional bandwidth, very wide field of view, and stationdependent effective Jones matrices (due to beam, ionosphere, Faraday rotation), the Mueller matrices to be considered are not only nondiagonal, but also baselinedependent. To highlight some of the main complications associated with very wide fieldofview instruments, in Sect. 2.1 we describe in detail the structure of the linear operator introduced by Bhatnagar et al. (2008). In Sect. 2.2 we propose a method to approximately correct for the associated effects corrupting the image plane (see also Rau & Cornwell 2011, for other examples of the use of linear operators in the context of image synthesis and deconvolution).
2.1. Description of baselinedependence, DDE, and polarization effects
For convenience, in this section and throughout the paper, we do not show the sky term that usually divides the sky to account for the projection of the celestial sphere onto the plane, as this has no influence on the results. The DDE below are baselinedependent, since is the case for LOFAR, and the measurement equation formalism can properly model those effects (for extensive discussions on the validity and limitations of the measurement equation see Hamaker et al. 1996; Smirnov 2011, and see Appendix A for a short review of our requirements). If is the set of fourpolarization measurement (XX, XY, YX, YY) and G_{p} is the direction independent Jones matrix of antennap, then the corrected visibilities on baseline pq are , and we can write (1)where I is the fourpolarization sky, ⊗ the Kronecker product, Vec the operator12 that transforms a 2 × 2 matrix into a fourdimensional vector, and models the product of the effects of correlator, sky brightness and array geometry. The matrix is a 4 × 4 matrix, and throughout the text we refer to it as the Mueller matrix^{2}. We can also write Eq. (1) in terms of a series of linear transformations: (2)where are the visibility measurement points in the timefrequency block in which the DDE are assumed to be constant. If N_{pix} is the number of pixels in the sky model, the true sky image vector I has a size of 4N_{pix}, and it contains the full polarization information (XX_{x},XY_{x},YX_{x},XY_{x}) on the xth pixel at the 4x position. Then, contains the DDE, and is a (4N_{pix}) × (4N_{pix}) block diagonal matrix. On a given baseline (p,q), each of its 4 × 4 block is the Mueller matrix evaluated at the location of the xth pixel. F is the Fourier transform operator of (4N_{pix}) × (4N_{pix}). Each of its (4 × 4) blocks is a scalar matrix, the scalar being the kernel of the Fourier basis exp(−2iπφ(u,v,w,s)). The matrix is the uvplane sampling function for that visibility, of size , and is the diagonal matrix containing the weights associated with the visibilities.
We show in Fig. 1 that the Mueller matrix is nondiagonal for LOFAR. The various panels show the amplitude of the direction dependent Mueller matrix using our beam model alone (therefore it does not include Faraday rotation and ionosphere), for the baseline (p,q) = (0,1) in a given time and frequency slot. To minimize the offdiagonal elements, we computed in the direction s, where s_{0} is the phase center direction. Intuitively, the normalization of the Mueller matrix by makes the projected dipoles on the sky orthogonal at the phase centre (offdiagonal elements are zero there), while this becomes less true as we go further from the targeted location. The offdiagonal elements can be as high as ~10% at a few degrees from the pointed location and, unlike most interferometers, they cannot be neglected in the case of LOFAR.
We are generally interested in using the total set of visibilities over baselines, time, and frequencies, having 4N_{V} points. We can write (3)where ϵ is the noise, and a (4N_{V}) × (4N_{pix}) matrix made of all matrices on top of each other. Then has dimensions , where N_{block} is the total number of timefrequencybaseline blocks. We can also write , where has size 4N_{blocks}N_{pix} × N_{pix} and has all the effects on top of each other. ℱ is the blockdiagonal Fourier transform operator of size 4N_{blocks}N_{pix} × 4N_{blocks}N_{pix}, with all of its blocks 4N_{pix} × 4N_{pix} equal to the matrix F appearing in Eq. (2). and are the sampling and weight matrices of size (4N_{V}) × (4N_{blocks}N_{pix}) and (4N_{V}) × (4N_{V}) respectively. The transformation of Eq. (3) is exact. Estimating a sky from a sparsely sampled measured set of visibilities is less trivial, however.
2.2. Estimating a sky image: polarization effects
There are different ways to solve for I from the set of visibilities, and reversing Eq. (3) relies on the linearity of the sky term in the measurement equation. As mentioned above, our deconvolution scheme uses AProjection. This generalization of CSCLEAN is now better understood in the framework of compressed sensing, which includes other new techniques (see McEwen & Wiaux 2011, for a review). In this section we describe a method of approximately correcting the polarization effects in the image plane. To highlight the issues associated with polarization and the baseline dependence of the DDE, we simply write the sky least square solution as the pseudo inverse: (4)where the term is the 4N_{pix} dirty image, and is the (4N_{pix}) × (4N_{pix}) image plane deconvolution matrix. Its structure is rather complicated, and its size makes its estimation prohibitive.
In the simple case of a matrix being unity (no DDE), we can see that each 4 × 4 block number (x,y) of the matrix is the instrumental response to a source centered on the location of the xth pixel, evaluated at the yth pixel. Therefore, computing would involve estimating the point spread function (PSF) centered on the location of each pixel and inverting a 4N_{pix} × 4N_{pix} matrix. In the presence of nontrivial baselinedependent DDE involving nondiagonal Mueller matrices, the problem becomes more complex. However, we show below that under some assumptions, the operator (sometimes called the deconvolution matrix) affected by DDE is decomposable in a classical deconvolution matrix (containing information on the PSF), and a simple correction performed separately on each pixel.
Following the notation introduced above, we can write as a 4N_{pix} × 4N_{pix} matrix, with of size 4N_{blocks}N_{pix} × 4N_{blocks}N_{pix}. The latter matrix is block diagonal, and each of its 4N_{pix} × 4N_{pix} blocks describes the PSF of the instrument for a given baselinetimefrequency block. Their 4 × 4 xyblocks are scalar matrices, the scalar p_{t,ν,pq}(x,y) being the response of the instrument evaluated at the yth pixel to a source being centered at the xth pixel in the given baselinetimefrequency block. We then have (5)It is virtually impossible to compute this matrix, and this illustrates the difficulty of doing an image plane deconvolution in the presence of timefrequencybaseline dependent DDE. To apply a firstorder correction to the image plane, we assume that the DDE are constant enough across time, frequency, and baseline. Then we can write
where N_{t,ν,pq} is the number of baselinetimefrequency blocks, is a 4 × 4 matrix that is the average of over baselines, time, and frequency, and p(x,y) is the PSF stacked in baselines, time, and frequency. If the uvcoverage is good enough, then is block diagonal (p(x,y) = 0 for x ≠ y, p(x,y) = 1 otherwise). All the x ≠ y terms cancel out in the final product, and in the relevant part of the matrix the N_{pix} 4 × 4 blocks are on the diagonal. Applying to can then be done on each pixel separately, by computing the product , where I_{x} contains the full polarization information (XX_{x},XY_{x},YX_{x},XY_{x}) for the xth pixel. This provides a way to estimate an approximate leastsquare cleancomponent value from a flux density in the dirty image in a given direction s_{x}. The details of this imageplane normalization are discussed further in Sect. 3.1. Although a few assumptions underlie this normalization, from the simulations presented in Sect. 4.3 we argue that the complex 4 × 4 matrix normalization per pixel presented here brings clear improvement. However, it does not seem necessary for the sky solutions to convergence (see Sect. 7).
3. Implementation of AProjection for LOFAR
As explained above, the full polarization AProjection scheme has been described in Bhatnagar et al. (2008). However, for the VLA implementation only the diagonal Mueller matrix terms had to be taken into account. As explained in Sect. 2, LOFAR has very wide fields of view, and the baselinedependent Mueller matrices associated to the fourpolarization correlation products are nondiagonal. This basically means that each individual polarization cannot be treated independently of the others. In this section we describe in detail a few implementations of AProjection that allow correcting for the nondiagonal Mueller matrices. We propose optimizations in relation to the telescope architecture. We show in Sect. 4 that this algorithm can indeed deal with the heavily nondiagonal Mueller matrices associated with LOFAR’s wide fields of view.
Following Bhatnagar et al. (2008), we built our implementation of AProjection on top of a traditional Cotton & Schwab CLEAN algorithm (Schwab 1984). It is an iterative deconvolution algorithm that consists of two intricate steps. The CLEAN and AProjection algorithms perform a series of operations that can be described in the following way: (6)where the true sky estimate Î^{n + 1} at step n + 1 is built from its estimate Î^{n} at step n and Φ is a nonlinear operator. It basically estimates the deconvolved sky from the residual dirty image . The minor cycle performs operations that approximate as discussed in Sect. 2, and takes the estimated image vector to zero but at the strongest components above a certain threshold. For the sky solutions to converge, the predict step or degridding) has to be accurate and unbiased. AProjection is a fast way to apply or . If we separate the phase of the Fourier kernel in Eq. (1) as the sum of φ^{0}(u,v,s) = u.l + v.m and , then following Bhatnagar et al. (2008) invoking the convolution theorem Eq. (1) becomes (7)where ∗ is the convolution product, ℱ is a 2D Fourier transform, , and i and j index the polarization number (running over (XX, XY, YX, YY)). This method is efficient because the DDE are smooth on the sky, meaning the support of the corresponding convolution function can be small (no high frequency terms). Indeed, as shown in Fig. 1 the LOFAR station beam is very smooth, and depending on the field of view, the typical support is of the order of 5–11 pixels.
3.1. Naive implementation
The operation in Eq. (6) converts a sky model into visibilities. To apply this operator to a massive amount of data in an algorithmically efficient way, we apply the scheme outlined in Eq. (7). First a twodimensional, fast Fourier transform (FFT) is applied to the sky model image Î, independently on the 4 polarizations. Then for each baseline in a timefrequency block Δ(t,ν), where the DDE are assumed to be constant, the 16 convolution functions are computed as the Fourier transform of the image plane Kronecker product of the DDE. For the LOFAR beam on a given baseline this block is typically 10 min and 0.2 MHz wide. The residual visibilities in each polarization are interpolated from the sum in Eq. (7). In practice, to minimize the support of the Wterm and associated aliasing effects, we multiply the DDE in the image plane by a Prolate spheroidal function (see Appendix B for more details).
The predicted visibilities are removed from the measured visibilities by computing the residual visibilities . Then applies a correction to the residual visibilities, projects the result onto a grid (the gridding step) and Fourier transforms it. In practice this is done as (8)and (9)The resulting dirty image is still corrupted by the phasedarray beamrelated effects discussed in Sect. 2.2. Before the minor cycle we can either multiply each xth fourpolarization I_{dirty}(x) pixel by the 4 × 4 matrix , or simply normalize each polarization of I_{dirty} by , the diagonal elements of . As shown in Sect. 4.3, fully applying to each pixel in I_{dirty} lead to a minor improvement.
The computational cost of taking DDE into account using AProjection depends on (i) whether its baseline dependence; (ii) the angular size at which the effect needs to be sampled (thereby constraining the size support of the convolution function); and (iii) the amplitude of the terms of the 4 × 4 matrix. In the case of the full polarization AProjection, the data need to be corrected per baseline, time, and frequency slot. For each of those data chunks, to recover the corrected fourpolarization visibilities, one needs to take all 16 terms of the 4 × 4 Mueller matrix into account, along with the four visibilities built from the 2D Fourier transform of the sky model. Therefore in addition to the 16 convolution function estimates per baseline and timefrequency slot, in the gridding and degridding steps one needs to compute operations per fourpolarization visibility point, where N_{S} is the support of the convolution function. The algorithmic complexity is discussed further in Sect. 4.4, but this implementation is too slow to be practical.
3.2. Separating antenna beam pattern and beamformer effects
Depending on the assumptions and system architecture one can find powerful algorithmic optimizations. We show here that in the case of LOFAR, we can use the fact that although stations are rotated with respect to each other, the elementary antennas are parallel. The effective phasedarray beam B_{p,s} of station p is modeled as B_{p,s} = a_{p,s}.E_{p,s}, where a_{p,s} is the array factor, and E_{p,s} is the element beam pattern. The term a_{p,s} depends on the phased array geometry and on the delays applied to the signal of the individual antennas before the summation (by the beamformer of each individual LOFAR stations). The term E_{p,s} models both the individual element antenna sensitivity over the sky and its projection on the sky. We have (10)where a_{p,s} is scalar valued and E_{p,s} nondiagonal because (intuitively) the element beam projection on the sphere depends on the direction. Applying the convolution theorem to Eq. (10) we obtain (11)All LOFAR stations have different layouts on the ground, so the scalarvalued array factor a_{p,s} is station dependent. However, all individual antennas of all stations point in the same direction, and we can assume that the Mueller matrix is baseline independent, i.e. . Although this requires an additional correction step of the gridded uvplane visibilities, as we see below this is an interesting algorithmic shortcut because the elementbeam correction can be applied to the baselinestacked grids.
The element beam is very smooth over the sky, and in most cases it can be assumed to be constant on time scales of an hour, so that the polarization correction step does not need to be applied often. The degridding step goes as follows. (i) In each timefrequency slot Δ(t,ν)_{E} where the Mueller matrix of the element beam is assumed to be constant, the polarization correction is applied to the (XX, XY, YX, YY) grids as the sum of convolved grids by the terms of . We then loop^{3} over baseline (pq), and timefrequency range Δ(t,ν)_{a} where the array factor and wcoordinate are assumed to be constant within Δ(t,ν)_{E}. For each step in the loop; (ii) the oversampled convolution function for baseline (pq) is estimated in Δ(t,ν)_{a} for the term in the second line of Eq. (11); and (iii) it is used to interpolate the predicted visibilities at the given uvcoordinate, separately on each polarization.
As explained in Sect. 4.4, the computing time for estimating the convolution functions can be quite significant, and this scheme allows us to compute only one convolution function per baseline instead of 16, and 4 gridding/degridding steps instead of 16. We note, however, that the assumption of baseline independence of the Mueller matrix on which this optimization is based starts to become invalid in case of directiondependent differential Faraday rotation, or for the longer baselines where the curvature of the earth starts to be important (in that case the element beams are not parallel). As discussed in Sect. 4.4, the computing time of this implementation is dominated by the convolution function estimate.
3.3. Separating the Wterm: hybrid wstacking
The support of the Aterm is determined by the minimum angular scale to be sampled in the image plane. The beam or ionospheric effects are in general very smooth on the sky so that only small numbers of pixels are needed to fully describe the effects, thus corresponding to a small convolution function support size (typically 11 × 11 pixels). The highest spatial frequency in the image plane is the Wterm and its support can be as large as ~500 × 500 pixels for the long baselines in wide fields of view, when the target field is at low elevation. This forces us to (i) compute a convolution function with a large support; and (ii) grid each individual baseline using the large, Wtermdominated convolution function.
We note, however, that the Wterm is in itself baseline independent^{4}: two baselines characterized by different ionosphere and beams, but with the same wcoordinate, will have the same Wterm. We therefore slightly changed the piping of the algorithm at this point by taking the ATerm and the Wterm into account separately, as follows (12)We consecutively grid or degrid the data in wslices, i.e., which have similar wcoordinates. This algorithm is also known as Wstacking^{5}. In addition, we take into account that the points can lie either above or below the associated wplane central coordinate, using the term exp(−2iπφ^{1}(Δw,s)^{)}, where Δw = w − w_{plane}. This step is similar to the traditional wprojection algorithm. If we have enough wstacking planes, Δw is small, and the support of the baselinetimefrequency dependent convolution function remains small, leading to a dramatic decrease in the total convolution function estimation time. Conversely, given a convolution function support we can find the maximum usable Δw and derive the number of wstacking planes as a function of the observation’s maximum w coordinate (see Sect. 4.4.3 for more detailed discussion). In the case of LOFAR, choosing a convolution function support of ~21 pixel gives a number of wstacking planes of ~30.
This requires yet an additional step as compared to the implementation described in Sect. 3.2, the degridding step , which can be described as follows. First, following the notation introduced above, (i) in each timefrequency interval Δ(t,ν)_{E} we correct the fourpolarization grids from the element beam (including projection effects) using . Then (ii) we loop over the number of wplanes (ranging from −w_{max} to w_{max}, see Sect. C), and convolve the grid obtained in (i) by the associated wterm (appearing in the second line of 12). Finally in (iii) for each wplane obtained in each step of the loop (ii), we loop over the set of baselines (pq)_{w} and the timefrequency range Δ(t,ν)_{a,w} associated with the given wplane. We interpolate the predicted visibilities at the given uvcoordinate, separately for each polarization, based on the oversampled convolution function. As discussed in Sect. 4.4.3, this is the fastest implementation of AProjection we have obtained so far.
4. Simulations
Fig. 3 Essential part of the AProjection algorithm relies on the predict step, which transforms a 2D sky image (the projection of the celestial sphere on a plane) into a set of visibilities. We have simulated a dataset having one offaxis source, where the true visibilities (black dashed line) have been estimated using Eq. (1), taking the beam into account. This plot shows the comparison in all measured polarizations between the exact value of the visibility of a given baseline and the AProjection estimate (gray line). Contrary to a traditional predict step, the visibilities are modulated by the beam amplitude (dotted line), and we have timedependent polarization leakage. The overplotted graph shows a zoom in the small region shown in the topright panel. In the degridding step, we use a computationally efficient closestneighbor interpolation, creating steps in the predicted visibilities. 
Fig. 4 Dramatic differences between the estimated deconvolved sky when using three different modified CLEAN algorithms on a simulated dataset containing 100 sources. The visibilities have been generated by taking (i) the individual antennas; (ii) their projection on the sky and (iii) the beamforming step (the scalar array factor) into account. The topleft image shows the deconvolved sky as estimated with a traditional imager not taking time, frequency, baseline, directiondependent effects into account. The top right and lower left images have been generated by considering, respectively, the array factor only and both the array factor and the element beam. The lower right panel shows that the input flux densities are correctly recovered. 
To test the algorithmic framework described above, we performed a series of tests on LOFAR simulated datasets. In this section we summarize those results and discuss the computational costs of AProjection for LOFAR.
4.1. One offaxis strongly polarized source
As discussed above, by using AProjection we can compute the exact value of the fourpolarization visibilities on a given baseline at a given time and frequency, from (i) the sky model and (ii) the baselinetimefrequencydirection dependent effects. In a first step, we focus on testing the full polarization AProjection degridding (or predict) step (). The accuracy of this step is vital for the convergence of the CLEAN/AProjection algorithm. Our experience suggests that any small numerical systematic bias in this operation can lead to strong divergence of CLEAN.
To test this transformation, we simulated a dataset having only one strongly polarized offaxis pointlike source in a lowband 62 MHz LBA dataset. The fourpolarization visibilities were generated using BBS (BlackBoard Selfcal^{6}), which directly calculates the visibilities following Eq. (1). It takes the beams of the individual elementary antennas into account (and their projection on the sky), as well as the phasing of the individual antenna within a LOFAR station (the array factor). We located the simulated source a few degrees from the phase center, and its flux density is strongly polarized, in a nonphysical way (with Stokes parameters of I,Q,U,V = 100, 40, 20, 10 Jy). Figure 3 shows the real part of BBS and AProjection predicted visibilities on baseline (01). The residuals are centered at zero, and AProjection performs very well in predicting the fourpolarization visibilities, taking the complicated effects of the LOFAR phased array station beams into account. A traditional predict using simple Fourier transform, facets, or WProjection would suffer from much higher residuals, driving systematics in the deconvolution, and thereby limiting the dynamic range. Here the residual errors are dominated by the type of interpolation we use in the uv domain (closest neighborhood, see Appendix B for details).
4.2. Dataset with many sources
4.2.1. LOFAR station beam
To test our modified implementation of the entire CLEAN algorithm (involving gridding and degridding steps), we simulated a dataset containing 100 sources, with the source count slope following the 1.4 GHz NVSS source count (Condon et al. 1998). Most of these sources have flux densities between 10^{2} and 1 Jy, to which we added two much stronger sources of 100 and 10 Jy, respectively.
As in the dataset described above, the visibilities are generated using BBS6. We have considered the Jones matrices of the individual elementary antennas, as well as the arrayfactor (the beamforming step). As explained in Sect. 3.2 and shown in Fig. 3, because the LOFAR stations are rotated with respect to each other, all baselines will be affected differently by beam effects. We have applied firstorder corrections for the beam effect to the visibilities by computing , where D_{p}(s_{0}) is the Jones matrix of station p computed at the center of the field. This mostly compensates for the element beam effects, in particular the projection of the dipoles on the sky. However, as shown below, the LOFAR fields are large, and the projections of the dipoles vary across the field of view.
Figure 4 shows the restored images produced using three different modified CLEAN algorithms. Each of the images is ~3000 × 3000 pixels in size, and is ~6 degrees wide. We used 15 wstackings, 128 Δwplanes (see Sect. 3.3), a maximum baseline of 5 kλ, a Briggs^{7} weighting scheme, a cycle factor of 2.5 and 10 000 minor cycle iterations. The first map was generated using WProjection (Cornwell et al. 2008) as implemented in CASA^{8}. Strong artifacts are present around the brightest offaxis source, and the dynamic range is only 1:230. In the second image we used our implementation of AProjection taking the array factor only into account. This has reduced the effect of the residual visibility levels on each individual baselines and raised the dynamic range to ~1:3.400. In the third image we have taken into account all LOFAR beam effects: individual antenna sensitivity, spatially varying projections, and the array factor. The dynamic range increased to ~1:12.000 ^{9}.
The source flux densities in the restored maps were determined using the LOFAR source extraction software pyBDSM^{10}. As shown in Fig. 4, the input flux densities are very well recovered.
4.2.2. LOFAR station beam and the ionosphere
To test the ionospheric correction with AProjection, we simulated a dataset containing 100 sources, affected by a simulated ionosphere. The ionospheric effects essentially consist of a purely scalar, directiondependent phase, without including Faraday rotation. In addition to these purely scalar ionospheric phase effects, the visibilities are affected by the directiondependent LOFAR station beam effects discussed above.
Fig. 5 Dirty images at the location of a bright source in a simulated dataset. The images are 1deg in diameter. Without applying the AProjection correction for an ionospheric phase screen (left panel) the dirty image shows important distortions as a result of ionospheric effects. Application of this correction (right panel) clearly shows improvements. 
The ionosphere is modeled as an infinitesimally thin layer at a height of 200 km. The total electron content (TEC) values at a set of sample points^{11} are generated by a vector autoregressive random process. As described in van der Tol (2009) the spatial correlation is given by Kolmogorov turbulence, and the TEC values needed to construct the phase screen are determined using Kriging interpolation.
Fig. 6 This figure shows the deconvolved image synthesized from the simulated dataset described in Sect. 4.2.2. In the left panel, the ionospheric effects are not taken into account, and our deconvolution scheme naturally produces severe artifacts and high residuals in the reconstructed sky. The deconvolved image shown in the right panel has been estimated using our implementation of AProjection with the timedependent ionospheric phase screen. 
Figure 5 shows dirty images at the location of a bright source before and after AProjection correction of beam and ionosphere. This suggests that the dirty undistorted sky is properly recovered from the corrupted visibilities. We compare in Fig. 6 the cleaned image with and without the AProjection correction. These simulations demonstrate that AProjection can indeed deal with the complicated effects associated with the ionosphere.
Fig. 7 The evolution of the estimated flux density as a function of the major cycle iteration for one polarized offaxis source. From top to bottom, and left to right, the panels show results generated by using (i) WProjection only; (ii) AWProjection with the array factor only; (iii) AWProjection taking the full beam model into account; and (iv) image plane corrections as described in Sect. 2.2. 
4.3. Convergence
We have also studied the influence of the various corrections described in this paper on the convergence speed of the flux densities estimated through the major cycle loop. For this we used the dataset described in Sect. 4.1 containing one offaxis source having Stokes parameters (IQUV) = (100,40,20,10) Jy.
Figure 7 shows the evolution of the estimated flux density as a function of the major cycle iteration number for different algorithms. Using only the WProjection merely shows a weak convergence towards the observed flux density (i.e., to the “beammultiplied” sky). The situation improves somewhat when we apply a full polarization AProjection without the element beam (therefore assuming the dipole projections on the sky are constant across the field of view). We note that in the absence of polarization, the situation is not as bad. Taking the element beams into account, the algorithm makes the estimated flux densities converge to the true values to within one percent. In this version of the algorithm, the image plane correction is the same in all polarizations and just a scalar, the average primary beam normalization. The situation improves slightly in terms of convergence speed by applying the imageplane renormalization described in Sect. 2.2.
In any case, the accuracy of the recovered flux densities seems to be guaranteed by the accuracy of the degridding step. Our experience in implementing AProjection suggests that any small systematic error in the degridding step can lead to biases in the recovered flux densities or even to divergence in the more severe cases. Adding the imageplane polarization correction seems to improve the convergence speed, but does not appears to be necessary.
4.4. Computational and memoryrelated costs
Given the large amount of data that has to be processed for imaging an interferometric dataset, reducing the algorithmic complexity is of primary importance.
4.4.1. Memoryrelated issues
The Aterm is generally very smooth in the image plane, with corresponding small support convolution functions, and under those conditions the AProjection algorithm is virtually free interms of computing time as explained in Bhatnagar et al. (2008) (and in Cornwell et al. 2008, in the case of WProjection): the Aterm and the low wcoordinate convolution function support is less than, or comparable to, the spheroidal function support that needs to be applied anyway in a traditional gridder/degridder. However, as described in Sect. 3.2, all LOFAR stations are rotated with respect to each other, and the synthesized stations beams are all different in a given direction. This gives rise to a serious algorithmic complication because it makes convolution functions baselinedependent: the number of possible convolution functions are 16 × N_{times} × N_{Freqs} × N_{Stations} × (N_{Stations} − 1)/2. With ~800 to ~1500 baselines, even with only one convolution function every ten minutes, a typical observing run of eight hours, one frequency channel, and an average support size of ~30 pixels (taking the wterm into account), this requires ~100 Gbytes of storage. This is optimistic because in the case of other DDE such as the ionosphere, the Aterm will have to be evaluated every 30 s. Even if the storage is done at the minimum necessary resolution, those numbers are clearly prohibitive for storing the convolution into memory. The convolution functions therefore have to be computed on the fly, and algorithmically, this represents an additional cost.
4.4.2. Naive and elementbeamseparated AProjection
The LOFAR station beams are smooth on the sky, and the corresponding convolution function supports are small (typically 11–15 pixel complex images), while the Wterm needs up to ≳ 500 pixels depending on the wcoordinate (with an average of ~30 pixels). The computing time depends on the convolution function support N_{S}, and for the implementations described in Sects. 3.1 and 3.2, we have (13)where the subscripts S, A_{p}, A_{q} and W stand for the prolate spheroidal (~7–11 pixels), Aterms of antenna p and q, and Wterm, respectively. For a typical fieldofview, we have N_{S}(A) = 9−15 and , where D is the fieldofview diameter and w is the Wcoordinate of the given baseline in the given timeslot (see Appendix B). For that baseline, time and frequency slot we can write the total computing time as , where and are the estimated gridding and convolution function times, respectively. For LOFAR data, in most cases we have N_{S} ~ N_{S}(W) and . Our experience has shown that t_{CF} is dominated by the computing time of the Fourier transform of the zeropadded convolution function (see Fig. 8), whose size is , where is the oversampling parameter that controls the quality of the nearest neighborhood interpolation. If is the number of visibility buffers associated with the Wplane (ΔT_{win} and Δν_{win} are the time/frequency window intervals in which the DDE are assumed to be constant), then we have . The gridding time for a given wplane is simply , where is the number of visibilities associated with the wplane. The total computing time can then be written as (14)where a and b are constants, N_{el} the number of time/frequency blocks in which the element beam is assumed to be constant, and is the computing time needed to apply the element beam to the grids (c = 16 for full polarization imaging, and c = 8 for only Istokes). For the implementation described in Sect. 3.1, we have t_{el} = 0, but both a and b are 16 times higher (8 for only Istokes) than in the case of the algorithm described in Sect. 3.2. Figure 8 shows the estimated gridding and convolution function times as a function of the wcoordinate of the given baseline in the given time slots. From this figure it is clear that the Wterm is the most important limiting factor and that the estimate of the convolution function represents a major limitation of those implementations, especially in cases of a rapidly varying DDE, such as the ionosphere, where the convolution function calculation would largely dominate.
4.4.3. Hybrid Wstacking and AProjection
As explained in Sect. 3.3 for the Wstaking implementation, the baselinetimefrequencydependent oversampled convolution function is the Fourier transform of the zeropadded imageplane product of the spheroidal and Aterms, plus a Wterm accounting for the Δwdistance between the given visibilities and the corresponding wplanes. The larger the support of the convolution function, the fewer wplanes we need to fully correct for the wterm. As explained in Appendix B, in order to properly sample the wterm in the image plane, to a given wcoordinate and field of view corresponds a certain convolution function support. We can then obtain Δw = a.N_{S}/D^{2} (with ), and the necessary number of wplanes between −w_{max} and w_{max} is N_{W} = w_{max}D^{2}/(a.N_{S}). Because we compute the Wterm convolution in the image plane, the cost of this step is . The total computing time is then (15)where N_{el} is the number of time/frequency buffers in which the element beam is assumed to be constant, N_{buf} is the number of time/frequency buffers in which the A/Δwterms are assumed to be constant, and b and c are constants. In Table 1 we present the typical computing times for a major cycle with the implementation discussed here and presented in Sect. 3.3.
Fig. 8 This figure shows the computing time of both the gridding and convolution function steps for a 10 degree diameter field of view, a maximum wcoordinate of 10^{4} m and a timewindow of T_{window} = 1200 s. When the DDE vary rapidly (such as in the case of the ionosphere), a baselinedependent convolution function is often required, and the convolution function computing time can largely exceed the gridding computing time. Our hybrid Wstaking algorithm is not affected by this scaling law. 
5. Summary and future work
The new generation of lowfrequency interferometers (EVLA, MWA, LWA) and pathfinders and precursors of the SKA (LOFAR, ASKAP, MeerKAT) cannot be used properly and efficiently without the development of new calibration and imaging techniques that take the many DDE that influence the electromagnetic field into account. These effects mainly include complicated antenna/station beam effects, fast ionosphere phase, and amplitude variations, along with the associated Faraday rotation.
In this paper we have focused on the issues associated with the application of new techniques for the LOFAR elementary antenna and station beams, which involves using few levels of phased arrays (see Sect. 1). Using the measurement equation formalism, the associated high complexity (wide field of view, individual station rotations, projection of the elementary dipoles on the sky), the problem of imaging and deconvolution LOFAR calibrated datasets can be solved by applying the AProjection algorithm described in Bhatnagar et al. (2008). Because of its very large field of view (~5–10 degrees in diameter), a fullpolarization AProjection implementation dealing with nondiagonal Mueller matrices is needed for LOFAR. We showed that AProjection can indeed deal with the heavily baselinedependent, nondiagonal, directiondependent Mueller matrices associated with LOFAR baselines. We also demonstrated that efficient ionospheric corrections can be performed using AProjection. We proposed a series of implementations of AProjection for LOFAR, taking nondiagonal Mueller matrices into account, aiming at accuracy and computing efficiency.
5.1. Optimizations
As explained in Sect. 4.4 the DDE, although varying rapidly, are smooth on the sky, so the convolution functions have a small support. However, in the case of LOFAR, the wide field of view requires considering (i) the Wterm; and (ii) the offdiagonal terms of the Mueller matrices (due to the varying projection of the dipoles on the sky or to Faraday rotation). In addition (iii), the baseline dependence renders the storage of the convolution functions prohibitive, so they have to be computed on the fly.
In our first implementation (Sect. 3.1), the constraint (i) often makes the Wterm convolution function support dominant; while (ii) requires the implementation of a complicated gridding and degridding step, taking all polarizations into account to correct for polarization leakage. In the case of a quickly varying DDE such as the ionosphere; step (iii) can completely dominate the computing time (usually set by gridding/degridding operations).
Using the fact that LOFAR station elementary antennas (responsible for the complicated projection effects) are parallel, although the station layouts are rotated with respect to each other, we could separate the first implementation into two steps (Sect. 3.2). The first is a purely scalar gridding/degridding, which suffers from (i) and (iii), while the second is only affected by (ii). As the nondiagonality of the Mueller matrix is corrected in the baselinestacked uvplane, we win a computational factor of 10 to 16 as compared to the first implementation. It is important to note that this optimization breaks down when the nondiagonal Mueller matrix becomes baselinedependent, as in the case of very long baselines (due to the earth’s curvature) or due to the ionosphere’s differential Faraday rotation.
In the last implementation (Sect. 3.3), we still applied the directiondependent nondiagonal Mueller matrix of the baselineindependent elementbeam, but go further by separating the Wterm from the ATerm. The former is responsible for the large support sizes and is baseline independent (around a given wplane), while the latter has small support and is baseline dependent. In this implementation, we grid/degrid the data based on a small support baselinedependent, oversampled convolution function, and convolve the input and output grids (forward and backward steps, respectively) with the nonoversampled Wterm convolution function. We saved computing time by calculating the oversampled convolution function of a generally much smaller support, giving a net gain of two to ten over the second implementation (Sect. 3.2, and Table 1).
Such optimizations are vital for the feasibility of the LOFAR extragalactic surveys, given their huge total integration times of hundreds of days. Also for the SKA, it will be important to take such algorithmic optimizations into account, which are linked to the instrument and system architecture, because they can reduce the algorithmic complexity by orders of magnitude.
5.2. Application to LOFAR survey data and future work
We conducted many experiments on real LOFAR data in order to test our imager. Specifically, we repeatedly observed the same set of objects, each time shifting their pointing centers by a few degrees. Our implementation of AProjection gives coherent corrected flux densities at the level of 5–10% on the edge of the field. However, although we have shown in this paper that the algorithm gives excellent results in simulated datasets, it shows little or no dynamic range improvement in the resulting images, as compared to a traditional imager.
LOFAR is a very complex machine and feeding the imager with an incorrect directiondependent calibration solution will not lead to any improvement in the deconvolved sky, and can even decrease the dynamical range. To improve LOFAR’s dynamic range we will have to make progress in understanding the calibration aspects of DDE, especially those related to ionosphere and differential Faraday rotation. Much effort is being given to this, and DDE calibration algorithms of low complexity are under development or have already been achieved, such as SAGECal (Yatawatta 2008).
On the imager side, other ongoing developments include (i) implementation on GPU, of the gridding, degridding, or the convolution function calculations; (ii) compressed sensing in the image plane; (iii) uvplane interpolation techniques different from zeropadding FFT; and (iv) wideband AProjection Bhatnagar et al. (in prep.). Ionosphere and true beam calibration, in combination with peeling techniques, will hopefully allow us to use the framework presented in this paper to reach the high ~10^{5}−10^{6} dynamic range needed to construct the deepest extragalactic LOFAR surveys.
See for example Maxim Voronkov’s presentation at http://www.astron.nl/calim2010/presentations in the context of ASKAP.
See http://www.lofar.org/wiki/doku.php?id=public:documents:lofar_documents&s[]=bbs for a review of BBS functionalities.
With a robust parameter of 0, see http://www.aoc.nrao.edu/dissertations/dbriggs
It seems the accuracy of the output images of our imager is presently limited to ~10^{4} due to some numerical precision problem: throughout the code, singleprecision floatingpoint numbers are used, and the roundingoff of products and sums involved in the algorithm seem to generate limitations at this level. Detailed tests based on comparing of single and multithreading consistently show an instability at this level. For the LOFAR surveys we might not need higher accuracy however, since we plan to use directiondependent peeling to subtract the brightest sources.
See http://www.lofar.org/wiki/doku.php for more information.
For each station, the sample points correspond to the piercing points in the ionosphere of five directions located at the four corners and the center of the field of view. This way the ionospheric layer is sampled in the most relevant area and there are at least five sample points within the field of view of each station.
Acknowledgments
This work was done using the R&D branch of the CASA code base. We would like to thank the CASA Group and casacore team for the underlying libraries. We thank Sarod Yatawatta for useful discussions of the element beam modeling that lead to powerful optimizations. We are grateful to our anonymous referee and to Wim van Driel for his extensive and constructive look at the form and content of this paper.
References
 Bernardi, G., Mitchell, D. A., Ord, S. M., et al. 2011, MNRAS, 413, 411 [NASA ADS] [CrossRef] [Google Scholar]
 Bhatnagar, S. 2009, in The LowFrequency Radio Universe, eds. D. J. Saikia, D. A. Green, Y. Gupta, & T. Venturi, ASP Conf. Ser., 407, 375 [Google Scholar]
 Bhatnagar, S., Cornwell, T. J., Golap, K., & Uson, J. M. 2008, A&A, 487, 419 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Condon, J. J., Cotton, W. D., Greisen, E. W., et al. 1998, AJ, 115, 1693 [NASA ADS] [CrossRef] [Google Scholar]
 Cornwell, T. J., Golap, K., & Bhatnagar, S. 2008, IEEE J. Select. Top. Signal Proc., 2, 647 [NASA ADS] [CrossRef] [Google Scholar]
 de Vos, M., Gunst, A. W., & Nijboer, R. 2009, IEEE Proc., 97, 1431 [Google Scholar]
 Hamaker, J. P., Bregman, J. D., & Sault, R. J. 1996, A&AS, 117, 137 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Heald, G., McKean, J., Pizzo, R., et al. 2010, in ISKAF2010 Science Meeting [Google Scholar]
 McEwen, J. D., & Wiaux, Y. 2011 [arXiv:1110.6137] [Google Scholar]
 Mitchell, D. A., Wayth, R. B., Bernardi, G., Greenhill, L. J., & Ord, S. M. 2012, J. Astron. Instrument., 1, 50003 [NASA ADS] [Google Scholar]
 Rau, U., & Cornwell, T. J. 2011, A&A, 532, A71 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Schwab, F. R. 1984, AJ, 89, 1076 [NASA ADS] [CrossRef] [Google Scholar]
 Smirnov, O. M. 2011, A&A, 527, A106 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Sullivan, I. S., Morales, M. F., Hazelton, B. J., et al. 2012, ApJ, 759, 17 [Google Scholar]
 van der Tol, S. 2009, Ph.D. Thesis, TU Delft [Google Scholar]
 Yatawatta, S. 2008 [arXiv:0809.0208] [Google Scholar]
Appendix A: Measurement equation formalism
To model the complex directiondependent effects (DDE – station beams, ionosphere, Faraday rotation, etc.), we make extensive use of the measurement equation formalism developed by Hamaker et al. (1996), which provides a model of a generic interferometer. Each of the physical phenomena that transform or convert the electric field before the correlation is computed is modeled by linear transformations (2 × 2 matrix). If is a sky direction, and ^{H} stands for the Hermitian transpose operator, then the correlation matrix V_{pq} between antennas p and q can be written as (A.1)where D_{p,s} is the product of directiondependent Jones matrices corresponding to antenna p (e.g., beam, ionosphere phase screen and Faraday rotation), G_{p} is the product of directionindependent Jones matrices for antenna p (like electronic gain and clock errors). The matrix K_{p,s} describes the effect of the array geometry and correlator on the observed phase shift of a coherent wavefront between antenna p and a reference antenna. This effect is purely scalar, so it is reducible to the product of a scalar and the unity matrix and we can write , where (u,v,w) is the baseline vector between antennas p and q in wavelength units and 1 is the unity matrix. We then have , where the −1 term models the correlator effect when phasing the signals in the direction of w. Finally, the product is the sky contribution in the direction s and is the true underlying source coherency matrix [[XX,XY], [YX, YY]].
This elegant formalism enables us to model the full polarization of the visibility as a function of the true underlying electric field correlation. In a simple and consistent way it takes the direction dependent and direction independent effects into account. Indeed, most of the Jones matrices in the measurement equations have a fairly simple formulation, and radio calibration problems amount to finding the various components of G and E, given a sky model .
The measurement equation introduced above can be written in a more extended and continuous form that is better suited to imaging by applying the Vec operator^{12} to . We obtain (A.2)where ⊗ is the Kronecker product, Vec is the operator that transforms a 2 × 2 matrix into a dimension 4 vector.
Appendix B: Further algorithm details
In this section we describe some important details of the various implementations of AProjection for LOFAR. In particular, we make extensive use of the Prolate spheroidal function for resolution adaptation and zeropadding for uvplane interpolation.
As explained in Sect. 3, since LOFAR stations are characterized by different primary beams, the gridding and degridding steps are baseline dependent. Therefore the convolution functions cannot be computed once and then kept stored in memory as is done for Wprojection. Instead, they have to be computed on the fly (see Sect. 4.4). However, because the station’s beam model is fairly complex and costly to evaluate (coordinates transformation, estimate of the Element beam Jones matrix), we only store the 4polarization image plane beams in memory (their Jones matrices) at the minimal resolution. The necessary resolution is simply estimated as λ/(2.D_{station}), where D_{station} is the given station’s diameter.
The Wterm is also estimated only once and stored in memory at the minimal resolution. This amounts to finding the maximum frequency to be sampled in the image plane, and the corresponding number of pixels corresponding to the minimum support required for the Wterm convolution function. If the image has an angular diameter D_{im}, the resolution needed to properly sample to Wterm is the inverse of the highest spatial frequency, located in one of its corners. The support of the Wterm is then .
To interpolate the visibilities on the grid in the gridding step (or conversely in the degridding step), or to adapt the resolution of the A and Wterms we use a zeropadding interpolation. Since this scheme can produce artifacts due to the presence of sharp edges and aliasing problems, we have to make extensive use of a Prolate spheroidal function. It is computed at the maximum resolution in the image plane. We then Fouriertransform it, find its support , “cut” it down to that size, and store it in the uvplane (hereafter ).
For the various algorithms presented in this paper, we have to compute the products of various DDE in the image plane. For example for the algorithm described in Sect. 3.1, we have to adapt the A and Wterm resolutions before multiplying them in the image plane. We first have to find the support N_{S} of the convolution function as in Eq. (13). We first compute the image plane spheroidal at the resolutions of the A and Wterms ( and respectively), as follows (B.1)where ℱ is the Fourier transform, and the zeropadding operator that puts the input into a grid of size n. To estimate the Aterm interpolated on an N_{S} × N_{S} pixel image: (B.2)We obtain the image plane effects at the same resolution and multiply them according to the specific needs of the various implementations described in Sect. 3 and obtain the image plane product P_{im}. We then compute the oversampled convolution function as (B.3)where the resulting interpolated convolution function has pixels, with as the oversampling parameter. If the final image is N_{pix} × N_{pix} pixels, because we have used the spheroidal function, after applying the inverse Fourier transform, we have to normalize the dirty image by the spheroidal function .
As outlined above, all LOFAR stations are different, and the convolution functions are baseline dependent. The for loops described in Sect. 3 are therefore parallelized on baseline. For the degridding step from a common readonly grid, the residual visibilities are estimated independently using different threads. The code has been created in the LOFAR package and depends on the casacore and casarest packages.
Appendix C: Limiting the maximum wcoordinate
Fig. C.1 The Wterm increases for lower elevations of the target field. Using the WProjection algorithm, the computing time increases with w^{2}. This figure shows the normalized cumulative distribution of the wcoordinate (left panel) for a typical LOFAR observation, and the corresponding normalized cumulative computing time (right panel). We can see that rejecting the ~5% of the points with w > 10^{4} saves ~70% of the computing time. 
For the traditional interferometers operating at relatively high frequencies compared to LOFAR, in general the field of view is small enough so that the Wterm can always be neglected (). However, for the new wide fields of view, longbaseline interferometers, the Wterm is very important, and not taking it into account produces artifacts and considerably reduces the dynamic range. For a given field of view
or a given angular distance to the phase center, the importance of the Wterm increases as the wcoordinate value, i.e., when the targeted field is at low elevation. It is therefore important to stress that wide fields of view or long baselines do not directly mean that the Wterm will be important: irrespective of its baseline or field of view, a planar array that would observe at zenith would always have a null wcoordinate.
Algorithmically, for AProjection, the Wterm is one of the main limiting factors. Using the WProjection algorithm (Cornwell et al. 2008) and assuming the Wterm support is higher than the Prolate spheroidal support, the gridding time evolves as t_{grid,w} ∝ w^{2}.D^{4}, because the highest spatial frequency in the image plane has to be properly sampled. We found that in a typical LOFAR dataset this nonlinear behavior generally makes the ≲5% of the points with the highest wcoordinates responsible for ≳ 70% of the computing time (as in Fig. C.1). This small amount of data does not necessarily add to the overall sensitivity or resolution. Setting a wmax value above which the visibilities are not gridded significantly increases the computational efficiency, without losing sensitivity or resolution. Not that our hybrid Wstacking (Sect. 3.3) algorithm is not affected by these limitations.
All Tables
All Figures
Fig. 1 LOFAR stations are phased arrays characterized by a large field of view. X/Y polarimetric measurements made with it are therefore nontrivial compared to a radio telescope using dishes: we have to take into account the projection of the dipoles that are generally nonorthogonal on the sky and the fact that this angle varies across the field of view. This figure shows the Mueller matrix corresponding to baseline (01) in a given time and frequency slot, normalized by the Mueller matrix at the phase center (see text). Each pixel in the plot (i,j) shows the amplitude of the (i,j) Mueller matrix element in a certain direction s using a logarithmic scale. Even in this normalized version, the offdiagonal Mueller terms are as high as 10% and cannot be neglected. 

In the text 
Fig. 2 Correction that can be applied to the image before the minor cycle. This is a firstorder correction for the complicated phased array beam, which depends on time, frequency, and baseline. 

In the text 
Fig. 3 Essential part of the AProjection algorithm relies on the predict step, which transforms a 2D sky image (the projection of the celestial sphere on a plane) into a set of visibilities. We have simulated a dataset having one offaxis source, where the true visibilities (black dashed line) have been estimated using Eq. (1), taking the beam into account. This plot shows the comparison in all measured polarizations between the exact value of the visibility of a given baseline and the AProjection estimate (gray line). Contrary to a traditional predict step, the visibilities are modulated by the beam amplitude (dotted line), and we have timedependent polarization leakage. The overplotted graph shows a zoom in the small region shown in the topright panel. In the degridding step, we use a computationally efficient closestneighbor interpolation, creating steps in the predicted visibilities. 

In the text 
Fig. 4 Dramatic differences between the estimated deconvolved sky when using three different modified CLEAN algorithms on a simulated dataset containing 100 sources. The visibilities have been generated by taking (i) the individual antennas; (ii) their projection on the sky and (iii) the beamforming step (the scalar array factor) into account. The topleft image shows the deconvolved sky as estimated with a traditional imager not taking time, frequency, baseline, directiondependent effects into account. The top right and lower left images have been generated by considering, respectively, the array factor only and both the array factor and the element beam. The lower right panel shows that the input flux densities are correctly recovered. 

In the text 
Fig. 5 Dirty images at the location of a bright source in a simulated dataset. The images are 1deg in diameter. Without applying the AProjection correction for an ionospheric phase screen (left panel) the dirty image shows important distortions as a result of ionospheric effects. Application of this correction (right panel) clearly shows improvements. 

In the text 
Fig. 6 This figure shows the deconvolved image synthesized from the simulated dataset described in Sect. 4.2.2. In the left panel, the ionospheric effects are not taken into account, and our deconvolution scheme naturally produces severe artifacts and high residuals in the reconstructed sky. The deconvolved image shown in the right panel has been estimated using our implementation of AProjection with the timedependent ionospheric phase screen. 

In the text 
Fig. 7 The evolution of the estimated flux density as a function of the major cycle iteration for one polarized offaxis source. From top to bottom, and left to right, the panels show results generated by using (i) WProjection only; (ii) AWProjection with the array factor only; (iii) AWProjection taking the full beam model into account; and (iv) image plane corrections as described in Sect. 2.2. 

In the text 
Fig. 8 This figure shows the computing time of both the gridding and convolution function steps for a 10 degree diameter field of view, a maximum wcoordinate of 10^{4} m and a timewindow of T_{window} = 1200 s. When the DDE vary rapidly (such as in the case of the ionosphere), a baselinedependent convolution function is often required, and the convolution function computing time can largely exceed the gridding computing time. Our hybrid Wstaking algorithm is not affected by this scaling law. 

In the text 
Fig. C.1 The Wterm increases for lower elevations of the target field. Using the WProjection algorithm, the computing time increases with w^{2}. This figure shows the normalized cumulative distribution of the wcoordinate (left panel) for a typical LOFAR observation, and the corresponding normalized cumulative computing time (right panel). We can see that rejecting the ~5% of the points with w > 10^{4} saves ~70% of the computing time. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.