© ESO, 2011

Introduction

The field around the bright radio source 3C 147 is a favourite showcase for dynamic range (DR) demonstrations. 3C 147 itself is very bright (22 Jy at 21 cm), which ensures a high SNR for selfcal solutions, while the surrounding field boasts a spectacular collection of mostly point-like fainter sources. The absence of significant extended emission has allowed very accurate sky models to be constructed. It is then not surprising that 3C 147 was the first field to break the 10⁶ dynamic range barrier (de Bruyn 2006; de Bruyn et al. 2011, in prep.), on a single 12-h synthesis. This spectacular result was achieved with regular selfcal implemented in the NEWSTAR package. A major contributing factor is the relatively low level of beam-related DDEs at the WSRT, as discussed in Paper II (Smirnov 2011, Sect. 2.1)¹. An image of the 3C 147 field at 1 600 000:1 dynamic range is shown in Fig. 1. The making of this image is the subject of Sect. 1.

The very same benign properties that allow the WSRT to achieve record DR also make it a perfect instrument for studying DDEs. The latter are prominent enough to clearly show up in high-DR maps (see e.g. left inset of Fig. 1), but not severe enough to hinder the building up of very deep sky models during normal selfcal². The 3C 147 observations by de Bruyn et al. (2011, in prep.) are a perfect example of this. I have therefore decided to reprocess these data using MeqTrees, to see if DDEs can be eliminated through the use of a suitable form of the RIME.

In the presence of a dominant source (in this case 3C 147 itself), selfcal solutions will tend to subsume all effects in the direction of that source. DDEs will then manifest themselves as artefacts around other sources, which need to be at a certain distance from the dominant source for the effect to become apparent. Since the dominant source is usually placed at or near the pointing centre (i.e. on-axis), DDE-related artefacts are also called off-axis artefacts. The artefacts themselves are quite clear (Fig. 1, left inset): the nature of the effects responsible for them, far less so. At least four possibilities have been postulated:

• Pointing errors (Paper II, Smirnov 2011, Sect. 2.1.4);

• Differential tropospheric refraction (ibid., Sect. 2.2.4);

• Errors in antenna positions, including non-coplanarity;

• Other systematic coordinate errors.

One of the objects of this study was to narrow down these possibilities. Section 2 therefore analyses the differential gain solutions obtained during this calibration, with a view to characterizing the DDEs.

1. Calibration approaches and results

1.1. Observations and NEWSTAR reduction

The observational data in question were obtained by de Bruyn in 2003. A single 12-h synthesis was taken, using 8 × 20 MHz bands (of 64 channels each) from 1300 to 1460 MHz, with 30 s integration time. Due to a back-end problem, one of the cross-correlations was corrupted, so only the XX and YY correlations were usable. De Bruyn then successfully reduced the observation using the NEWSTAR package, achieving a world record 1 600 000:1 dynamic range (de Bruyn et al. 2011, in prep.). Only regular selfcal was done and no peeling was attempted, so the resulting image showed DDE-related artefacts around off-axis sources (Fig. 1, left inset). One result of the reduction was a very deep NEWSTAR-format sky model for the field, containing about 300 point sources. This provided a fantastic platform from which to begin my DDE study with MeqTrees. I had a ready-made sky model that was known to be good enough to reach the thermal noise with this particular dataset, and I had intermediate images from de Bruyn’s reduction that could be used as checkpoints.

The same observation was repeated in 2006 with somewhat different correlator settings (see de Bruyn et al. 2011, in prep., for details), and reduced in a similar manner.

For his NEWSTAR reduction, de Bruyn self-calibrated each channel independently, rather than explicitly calibrating for a bandpass (see below). This procedure is described in Paper II (Smirnov 2011, Sect. 1.1). It uses the following implicit form of the RIME: (1)and consists of the following steps:

• 1.

  Find that minimizes  | X_pq − D_pq |  in a least-squares sense. Compute “corrected data” as (The solution interval here was one timeslot, one frequency channel. Only 56 channels in each band were usable; these were further averaged down using Hanning tapering, so in the end only 28 frequency points per band were used.)

• 2.

  Find that minimizes . Compute “corrected data” as , where “ ÷ ” is element-by-element division – the inverse of “ ∗ ”. (The solution interval here was the full 12 h, one band.)

• 3.

  Compute “residual data” as .

This procedure was repeated for each band. Residual visibilities were imaged and summed across all 8 bands, then deconvolved using Högbom CLEAN. The sky model was then added back in using a Gaussian restoring beam.

1.2. Calibration in MeqTrees: an overview

In broad terms, selfcal in MeqTrees also consists of a least-squares fit of an equation such as (1) to the data. However, the following features are different:

• The structure of the equation is not fixed: arbitrary forms of theRIME may be constructed. Crucially for my purposes, these mayinclude DDE terms.

• The elements of the Jones matrices are not necessarily simple solvable parameters (though they may be), but can be represented by arbitrary functions. For example, rather than solve for E-Jones (or Z-Jones) elements directly, MeqTrees can derive them from some model of the primary beam (or ionosphere) and solve for the parameters of the model. An example of this is given in Sect. 1.6.

• Different solvables may have different solution intervals, even in a simultaneous solution.

Because of the inherent flexibility of MeqTrees, calibration can be pursued in a great variety of ways. During this study, a specific methodology was narrowed down, implemented and tested. This became the basis of the “Calico” calibration framework that is now included with MeqTrees. The steps (and terminology) of calibration with Calico are as follows:

• A desired form of the RIME is constructed, by selecting a skymodel, and picking a series of Jones terms (plus, optionally,interferometer-based errors). For example, a form similar toNEWSTAR’s implicit RIME (Eq. (1)), would be: V_pq is usually called the “corrupted predict”.

• 2.

  “Corrupted predict” is fitted to “data”. That is, the MeqTrees solver is instructed to find values of solvable parameters that minimize  | D_pq − V_pq |  in a least-squares sense. This can be (and usually is) done in multiple stages, e.g. first, followed by , etc.

• Output visibilities are computed as either “corrected data” or as “corrupted residuals” R_pq = D_pq − V_pq, or as “corrected residuals”

Subsequent imaging steps are not considered part of Calico or MeqTrees per se, since they apply equally to data calibrated using any other means, and are accomplished via external tools such as the lwimager program (part of the casarest package). These steps may also differ from project to project. I routinely image the per-band corrected residuals , since these provide the best visual indicator of the quality of a calibration. For this particular study, I subsequently made 8-band residual images in MFS mode, using all data. The images were further deconvolved using Cotton-Schwab CLEAN (Schwab 1984), and the sky model was added back in using a Gaussian restoring beam.

While Calico can produce “corrected data” (for purposes of imaging, etc.) at any step of the reduction, it does not use it as input for subsequent stages like NEWSTAR and other 2GC packages do. The fitting at step 2 is always done using the original³ observed data D_pq, and calibration consists of building up the RIME until the “corrupted predict” fits the observations. Once DDEs enter the picture, “corrected data” (in the conventional understanding of data corrected for instrumental errors) do not really exist any more, since visibilities can only be corrected for the value of a DDE in a particular direction. Hence the Calico philosophy is to work with the original data at all times.

1.3. Bandpass selfcal

A per-channel, per-timeslot solution for G_p in Eq. (1) has 2N complex unknowns, and N(N − 1) complex measurements, where N is the number of stations (14 for the WSRT). While an acceptable ratio (this is, after all, why selfcal works in the first place!), it does not leave a lot of room for introducing DDE-related parameters. In general, we want to allow our solutions as few DoF’s as possible, and making the solution intervals (in time and/or frequency) larger is one way of ensuring this.

The WSRT has a reasonably stable bandpass, so an obvious way to reduce the parameter count is to separate G_p into a bandpass component to capture the frequency structure (with little to none variation in time), and a frequency-independent, rapidly varying complex gain. Per each station/receiver, this replaces N_chan × N_time parameters with only N_chan + N_time of them.

For the MeqTrees reduction, I therefore started with the following RIME: (2)Here, X_spq, E_s, G_p and M_pq have the same meaning as in the NEWSTAR equation above (including a similar smearing correction term in X_spq). The B_p term is a second diagonal Jones matrix representing the bandpass. Note that MeqTrees itself makes no special distinction between G and B. Both are generic diagonal complex Jones matrices, with solvable real and imaginary parts. It is only when we specify the solution intervals that these Jones terms acquire their intended meanings:

• The G_p solution interval is 30 s, all channels (thus, one independent solution every timeslot, per the entire band).

• The B_p solution interval was initially set to 30 min, and one channel (thus, one independent solution every channel, per 60 timeslots). Note that as the bandpass has significant structure, I did not attempt to fit it with any smooth function, but rather allowed each channel to be fitted as an independent parameter, with a timescale of 30 min. The latter interval was intended to accommodate slow drift in the bandpass.

The initial result of this calibration was profoundly unsatisfactory (Fig. 2, upper left). The residual image was dominated by spoke-like artefacts centred on 3C 147, at about 10 000:1 level (relative the flux of 3C 147 itself). These spokes correspond to edges of the 30-min solution intervals (being an E-W array, WSRT has a one-dimensional instantaneous PSF). The obvious explanation for the error is short-term bandpass instability. Decreasing the solution interval of B_p to 15 and 7.5 min reduced the artefacts to levels of 50 000:1 to 100 000:1, but did not eliminate them entirely. Smoothing the 7.5 min solution with a spline (along the time axis, per each channel independently) produced a marginal improvement (Fig. 2, lower right). Subtraction artefacts are still plainly visible in the map, although at a level not significantly above thermal noise.

At this stage I had to conclude that the WSRT bandpass exhibits some very low-level, but extremely short-term jitter, precluding a separate bandpass selfcal at extreme dynamic ranges. On the other hand, this result also shows that where a single-band dynamic range of no more than 100 000 is expected (as is the case for many other observations), bandpass selfcal provides perfectly adequate results, and can cut down on the number of solvable parameters significantly. In the meantime, for the 3C 147 study I had to revert to the per-channel selfcal approach of de Bruyn.

1.4. Per-channel selfcal

The RIME for per-channel selfcal is just Eq. (2) without the bandpass term. It is, in fact, very similar to NEWSTAR’s implicit Eq. (1)⁴: (3)Per-channel selfcal is achieved by setting the solution interval of G_p to one channel and one timeslot. The resulting single-band residual images are dominated by closure errors at a level of  ~100 μJy (or 1:200 000 relative to 3C 147 itself); these go away after an M_pq solution (Fig. 3). These images are qualitatively very similar to those obtained by de Bruyn during his NEWSTAR calibration (which is to be expected, given the similarity of the equations).

The remaining artefacts in Fig. 3 are associated with the three next-brightest sources⁵ in the field: B (42 mJy), C (52 mJy), and D (22 mJy). The furthest of these (B) is about 20′ away from centre. The artefacts themselves are under 50   μJy (thermal noise in one band being  ~30 μJy), or at a level of 1:1000 relative to the associated sources. This is why de Bruyn (2006) talks about an “off-axis dynamic range” of a 1000: while 3C 147 itself (22 Jy) is subtracted without a trace (down to the thermal noise), the off-axis sources are only subtracted to a precision of about 1000. Some of the artefact structure is doubtlessly due to slightly under- or overestimated sky model fluxes (this produces regular rings matching the WSRT PSF), which can be taken care of during subsequent deconvolution. Most of it, however, is due to DDEs and does not deconvolve, producing artefacts in the final 8-band images such as the one shown in the left inset of Fig. 1 (the inset is, in fact, a close-up of source B).

1.5. Interferometer-based errors

It is not clear what causes closure errors at the WSRT. Common sense suggests the analogue part of the signal chain is to blame, but there is no hard evidence either way. What is evident is that high-DR images exhibit concentric ring-like artefacts such as those in the left image of Fig. 3, and that these go away once a solution for an interferometer-based multiplicative error – the M_pq term of Eq. (3) – is applied. A single solution per band, per the entire 12 h (per correlation and interferometer) is sufficient. The M_pq solutions are usually within 10^-3−10^-4 of unity (as is the case here), but can be much higher in some observations, for reasons that remain mysterious.

The latter fact suggests an intrinsic time variability, but solving for M_pq on short time intervals is very dangerous. Any solution for M_pq will also try to compensate for observed flux that is not present in the sky model. Unless the solution interval is sufficiently long, there will be unmodelled sources with a fringe rate slow enough that their vector average visibility over the solution interval will be significantly non-zero. These sources will then tend to be attenuated by the M_pq solutions. The 3C 147 observations provide a perfect example (Fig. 4). On the left is an 8-band residual image with a 12 h M_pq solution; on the right is one with a 30 min solution. Model sources are indicated by crosses. Attenuation of unmodelled sources towards phase centre is clearly visible in the right image.

This implies that closure errors cannot be reliably solved for on observations shorter than 12 h, unless a “complete” (i.e. down to the noise) sky model of the center of the field is available.

1.6. Pointing selfcal

Since the main cause of artefacts in WSRT images is commonly considered to be pointing error (see Paper II, Smirnov 2011, Sect. 2.1.4), I decided to implement a form of the pointing selfcal algorithm suggested by Bhatnagar et al. (2004). This proved to be a straightforward exercise in MeqTrees, since only a small modification of the RIME of Eq. (3) was required: (4)Instead of a per-source beam gain E_s = E(l_s,m_s) (where E(l,m) is the cos³ approximation given by Paper II, Smirnov 2011, Sect. 2.1.1), which I had been using in the previous equations, here I used a per-antenna, per-source beam gain E_sp, defined as follows: (5)Per-antenna pointing offsets δl_p,δm_p were then treated as solvable parameters.

The results of this proved inconclusive. Even though the solution converged to some physically-sensible offsets (on the order of arcmin), no tangible reduction of artefacts was observed in residual images. This could be due to a number of reasons:

• 1.

  The cos³ approximation is not good enough – unlikely, as it has been independently verified at least for the core of the main lobe, which the sources in question sources are well within.

• 2.

  With only a few sources sufficiently bright to exhibit DDEs, we simply don’t have enough constraints for a pointing solution on this field.

• 3.

  The model fluxes/positions for the sources in question are not sufficiently accurate.

• 4.

  The dominant DDE affecting this observation is not pointing error. This will be elaborated on further in Sect. 2.

• 5.

  There is something wrong with my implementation of pointing selfcal, especially since the figures in Sect. 2 suggest that mispointings are detectable.

Trying to get a better grip on the problem, I eventually settled on a “controlled experiment”: locating a field containing multiple bright off-axis sources, and observing it with deliberately exaggerated mispointings, to see if these can be more readily recovered from the data. This experiment became known as the “QMC Project” (in honour of the long-defunct Quality Monitoring Committee of WSRT), and was successfully carried out. The results of this are still being processed, and will be presented in a follow-up paper. In the meantime, I had to look for alternative approaches to DDEs in the 3C 147 field.

1.7. Differential gains: the “flyswatter”

In the spirit of “phenomenological” equations discussed in Paper II (Smirnov 2011, Sect. 1.3), I decided to introduce a differential gain Jones (ΔE-Jones) into my form of the RIME: (6)The ΔE_sp term is meant to subsume all DDEs associated with source s and antenna p (with the exception of the nominal primary beam gain, which is already represented by E_s). Solving for this term requires some caution, lest too many DoF’s be introduced into the equation. I approached this as follows:

• ΔE was fixed at unity for all sources except B, C, and D;

• For B, C and D, ΔE was set to a diagonal matrix with solvable complex elements;

• The solution intervals for the ΔE elements were set to one solution per 30 min, per entire band (and per source, antenna, receptor).

Figure 5 shows the effect of ΔE solutions. All artefacts associated with sources B, C and D disappeared completely, to the point where artefacts around fainter sources became just about visible. These were later eliminated once ΔE solutions were enabled for those sources. In fact, any source to which a solvable ΔE term was assigned promptly vanished from the residual maps, hence my name for the differential gains approach: the “flyswatter” algorithm.

Once 8-band residual images were created, the increased sensitivity made it apparent that four more sources were exhibiting a small amount of DDEs. A more in-depth look at the sky model also showed that all seven sources were in fact slightly extended; the NEWSTAR model represented each by a tight cluster of point sources. Since all “sources” in such a cluster are subject to the same DDE, it seemed sensible to make sure the same ΔE_sp term was applied to all of them. This was easily done using Tigger, a sky model manager/viewer tool included with MeqTrees. Tigger automatically detects source clusters and assigns source identifiers appropriately (e.g. B232, B232a, ... B232e) It was then a simple matter to tell the Calico framework to use the same ΔE_s0p term for all sources s associated with cluster s₀.

The seven source clusters for which differential gain solutions were eventually obtained are summarized in Fig. 6. Two of them are somewhat noteworthy. Source C270 is very close to centre, and therefore shouldn’t be affected by DDEs as much as the other sources. It is, however, a complicated and highly polarized source, so perhaps the artefacts it exhibits after regular selfcal are primarily due to sky model inaccuracies, which the ΔE solutions absorb (see discussion in Sect. 2.1). Source ae317 is almost the opposite: it is very faint, but far enough off-axis to be in a sidelobe of the primary beam, and so subject to especially severe DDEs.

1.8. The showcase result

The ultimate result of my calibration of the 2003 observation is shown in Fig. 1. This image is a true showcase for the differential gains approach. The precise steps leading to this image were as follows:

• 1.

  Each of the 8 bands was independently calibrated usingper-channel selfcal, interferometer-based errors, andΔE solutions on seven source clusters, as described above. Corrected residuals were generated.

• 2.

  The residuals for all 8 bands were imaged together (in MFS mode) to produce a single residual image. This revealed a large number of fainter sources not visible in the per-band maps.

• 3.

  The 8-band image was deconvolved using Cotton-Schwab CLEAN (Schwab 1984).

• 4.

  The sky model was added back into the deconvolved image, using a Gaussian restoring beam.

The resulting image is completely artefact-free. Presumably, all other sources in the field are too faint to exhibit any DDE-related artefacts. With a dynamic range of 1 600 000:1, this image is the deepest and cleanest single-synthesis radio map in the world to date.

1.9. Flyswatter limitations

While it can help produce spectacular images, the flyswatter has some serious caveats and drawbacks that need to be explored. First of all, it is a brute-force approach, in the sense that it squashes all effects into a single ΔE term. This includes inaccuracies in the sky model! Indeed, any missing source flux or error in source position can be accommodated with a suitable ΔE. Even unmodelled source structure can “leak” into differential gain solutions (Sect. 2.1). Thus, differential gains are good for subtracting sources, but at the cost of mashing up information on the source per se. (One does not use a flyswatter to probe a fly’s anatomy!)

Secondly, solvable differential gains can lead to a proliferation of DoF’s. Per-channel selfcal has N_ant unknowns per measurements, or unknowns per measurement; differential gains add unknowns per measurement, where N_time and N_freq are the sizes of the solution interval for ΔE. This ratio remains favourable for small N_src and large N_time and/or N_freq (as is the case for my 3C 147 reduction), but one must be careful.

The third caveat is processing cost. While usually not as expensive in terms of I/O or CPU as peeling (which, in addition to the solutions themselves, requires repeated subtraction and phase shifting steps), the flyswatter is not free. Every source with a differential gain solution adds 4N_ant unknowns (assuming a diagonal complex ΔE_sp term, hence 4 real values per matrix) to the equations. As the number of unknowns (N_unk = 4N_antN_src) grows, inversion of the normal matrix within the least-squares solver becomes a CPU bottleneck, since it scales as . This makes it impractical to solve for ΔE’s for more than a handful⁶ of sources at a time.

One way to mitigate the solver bottleneck is to decompose the ΔE_sp’s into nearly-orthogonal sets of unknowns. For example, we can treat the set of ΔE_sp’s associated with one source s as independent from all other sources. The (4N_antN_src)² normal matrix inside the solver then becomes block-diagonal, composed of N_src blocks of size (4N_ant)². Inversion of this matrix then scales as O. This scheme was tested in MeqTrees, and it was found that the trade-off is slower convergence, requiring more iterations. For large numbers of sources, however, this becomes very favourable.

2. Analysis of differential gain solutions

It is time to see whether any useful information can be gleaned from the differential gains solutions themselves. As a result of the reduction, I had obtained: per each source direction (7 of these: see Fig. 6), per each antenna (14), per each band (8), per 30-min interval (24 of these in a 12-h synthesis), two complex numbers representing the apparent differential gain of the X and Y dipole (“differential” being relative to the gain in the direction of 3C 147 – almost at the centre of the field – which had been taken care of by regular selfcal).

I then adopted the following approach. Given the relatively low fractional bandwidth, I didn’t expect much variation with frequency in ΔE. I therefore treated each set of 8 per-band solutions as independent samples of the same variable. The mean of the 8 samples was used as an estimator of variable, and the standard deviation of the 8 samples as an estimator of the error (i.e. the error bar).

Since it quickly became apparent that the ΔE solutions were exhibiting some very interesting behaviour, I applied exactly the same procedure to the 2006 observations, so that comparisons could be made. The plots below show the results from both observations.

Figure 7 is a summary of the differential gain-amplitudes per source, per antenna. The precise quantity plotted here was computed as follows. First, I computed the norm of the ΔE matrix (diagonal by construction) as

I then normalized (divided) this value by the mean value per source (that is, the mean across all time intervals, bands, and antennas), which was meant to take out the effect of incorrect model fluxes (see below). The resulting “normalized norm” was then plotted as a function of time, per source, per antenna.

Figure 8 is a similar plot of the differential gain-phases, computed as

The most striking feature of Figs. 7 and 8 is the high SNR. They show a high degree of temporal continuity in the solutions, and statistically significant structure. This strongly suggests that the solutions represent real physical or numerical effects. As to the nature of these effects, I still do not have satisfactory answers, though it is hoped that the “QMC Project” mentioned earlier will shed some more light on the issues. The rest of this section discusses some of the more prominent questions, and proposes some rather speculative explanations.

2.1. Absorbing errors in the sky model

As already mentioned, a major caveat of the flyswatter approach is that ΔE solutions will tend to absorb inaccuracies in the source model. By analogy (and for exactly the same reasons), classic selfcal alone cannot solve for absolute positions or fluxes. Indeed, if the true position of a source l = (l,m) is offset from the model position l^(mod) = l + δl = (l + δl,m + δm), while the true brightness differs from the model brightness by a multiplicative matrix factor: (the latter being a straightforward generalization of a scalar factor a²), then the coherency term of the RIME for the model source may be written out as If a solvable differential gain is then assigned to the source, the model can be made to fit the data by absorbing the K_p(δl)A factor into the ΔE_ps solutions.

Even more insidiously, differential gains can absorb some source structure. Consider a source that is slightly extended in one direction, enough to be resolved on the longest baselines. An E-W array like the WSRT has a one-dimensional instantaneous fan beam. It will “see” the source as a point source when the fan beam is aligned with the source orientation, and start resolving it when the fan beam becomes perpendicular to the source. In other words, the apparent flux of the source will remain constant in time on short baselines, and vary in time on the long baselines as the source resolves. If such a source is represented by a point source in the sky model, the model flux will be constant on all baselines. Now, if some antennas are predominantly involved in long baselines (RTC and RTD, in the case of WSRT), ΔE solutions can compensate for some of the flux discrepancy by changing the gain-amplitudes of these antennas. I would expect to see a variation of  || ΔE ||  with a 12-h period. Since most of the baselines to RTC and RTD are mutually redundant (0-C equals 1-D, etc.), their variation in  || ΔE ||  should be very similar.

This is exactly what we’re seeing in Fig. 7! The plots very strongly suggest that the top three sources (B, C and D) are indeed slightly extended (more so than in the model, that is). If this is the case, then the dominant contribution to  ||ΔE||  on antennas RTC and RTD is due to source structure rather than any actual DDE.

2.2. Amplitude behaviour

If the behaviour of differential amplitude on RTC and RTD is due to source structure, this still leaves effects on the other antennas unexplained. First let us consider what a pointing error would look like. With the exception of source ae317 (for which the solutions are much too noisy anyway), the other six sources are well within the main lobe of the primary beam, where the beam gain can be expected to decrease smoothly with distance from pointing centre. We should therefore expect to see  ||ΔE_sp||  > 1 if antenna p mispoints towards source s (the source appears brighter on antenna p), and  ||ΔE_sp||  < 1 if it mispoints away. WSRT dishes are equatorially mounted, so mispointings due to mechanical or electronic errors in either axis drive would be stationary with respect to the sky, and thus cause constant  ||ΔE||  offsets. Mispointings due to wind pressure, thermal or gravitational deformation, on the other hand, would be intrinsically time-variable (but would perhaps correlate between adjacent antennas).

Quite a few plots in Fig. 7 do show (mostly) static offsets. It can be illuminating to present  ||ΔE||  in a format I call the “rogues gallery” (Figs. 9 and 10). This shows, for each of the 14 antennas, a 12-h average  ||ΔE||  per source, using circles of varying size placed at the position of the source. The magnitude of ( || ΔE||  −1) is indicated by circle size, and the sign by colour. A static mispointing in, e.g., a Northern direction would show up as blue circles in the top half of the plot (i.e. sources appearing brighter), and red circles in the bottom half.

The galleries show exactly such a pattern for antennas RT5 through RT8 (and perhaps RT4), for both the 2003 and, to a lesser extent, the 2006 observations. It is a little bit strange that 4 (or even 5) adjacent antennas would so consistently mispoint North, and do the same three years later. Perhaps this is another, poorly understood consequence of unmodelled source structure. Some pointers to this are that antennas RT4–8, being in the middle of the array, form up predominantly shorter baselines, and that the long-baseline antennas RTC and RTD exhibit the opposite behaviour. (What hinders such an analysis is the unfortunate fact that the three brightest off-axis sources all exhibit some structure, and all three lie in the bottom half of the field.) Another puzzling feature is the consistently low  ||ΔE||  for sources F, H and K on antenna RTC in 2003 (and to a far lesser extent in 2006). If due to source structure, why does it not repeat on RTD? Perhaps RTC is mispointing to the South?

Antennas RT0–2, RT9 and RTA, on the other hand, show completely different patterns, with little to no similarity between 2003 and 2006. Some of these are consistent with a static mispointing. Some antennas (RT8 and RT9, and RTB especially) also show a hint of time variability in  ||ΔE|| .

In any case, it is clear that the complicated interaction between source structure and differential gain-amplitudes makes the latter extremely difficult to interpret. Note also that my (or rather de Bruyn’s) source model was built by NEWSTAR based on regular selfcal, so there’s bound to be some contamination from DDE-related artefacts in the source parameters. Truly robust methods for disentangling source structure from DDEs have yet to be developed. It is also clear that an approach that parametrizes the DDEs in a “global” way, such as pointing selfcal, is the way forward – what is not yet clear is how much the global solution itself can be affected by unmodelled structure in the brighter sources, and what to do about it.

2.3. Phase behaviour

A rather prominent feature of the phase plots in Fig. 8 is their continuity from antenna to antenna, and the fact that the phases at the two ends of the array exhibit opposite temporal trends. This is suggestive of an evolving phase slope over the array. Fitting such a slope (per source) produces some very striking results (Fig. 11). The dominant phase effect is clearly global rather than antenna-based, and is extremely consistent across both observations.

A phase slope over the array can be interpreted as an apparent position offset. It appears that the slope behaviour in Fig. 11 can be fitted quite well by constant position offsets. The best-fitting position offsets are indicated in Fig. 12, and the corresponding slope curves are plotted in green in Fig. 11.

Figure 12 immediately suggests a field rotation. And indeed, the entire collection of phase slopes (for both the 2003 and 2006 observation), is, to first order, consistent with a rotation of 45″ around the phase centre. The corresponding slope curves are plotted in red. While there are some significant differences in the brighter sources, it seems clear that the dominant effect is not an instrumental DDE at all, but a systematic rotation of the sky model. The model positions are derived by NEWSTAR from direct fits to the visibilities, and de Bruyn (priv. comm.) has independently cross-checked the positions of distant sources against the NVSS, which seems to preclude a rotation in the model itself. Note that a 45″ rotation can also be introduced by a clock error of about 2.9 s, or a corresponding rotational error in conversion of uvw coordinates from apparent to J2000. Since NEWSTAR and MeqTrees use completely different tool chains and visibility data formats, I cannot exclude a coordinate conversion error somewhere along the line. This needs to be urgently investigated. If indeed the entire sky model is slightly rotated, then perhaps the image of Fig. 1 can even be improved upon!

2.4. Feeding differential gains back into the sky model

The results above suggested that I could improve my sky model by feeding back in some information extracted from the ΔE solutions. In the previous section, I obtained a correction to the model positions of the seven sources⁷. Following the discussion of Sect. 2.1, I could also provide corrections for the I and Q fluxes by applying the per-source average ΔE amplitudes: where t_i and ν_j represent the time and frequency solution intervals of ΔE. In terms of the I and Q fluxes, the correction becomes: I therefore applied these corrections for I, Q, and position to my sky models (independently for the 2003 and 2006 observations), and repeated the calibration procedure. An improvement in single-band residuals was immediately apparent (Fig. 13) – after G_p and M_pq solutions, the seven off-axis sources subtracted noticeably better. Residuals after ΔE_sp solutions, on the other hand, looked pretty much the same (this is not surprising, since differential gains had already taken care of the visible off-axis errors in the original reduction), with a very slight improvement around 3C 147 itself, which can be explained by improved G_p solutions due to the more accurate sky model.

2.5. Phase behaviour II

Presumably, the remaining residual structures in Fig. 13 are more representative of the instrumental DDEs per se, since inaccuracies in the sky model have been significantly reduced. We should also expect the differential gain solutions to be more indicative of the actual DDEs (apart from the issue of resolved sources affecting  ||ΔE||  on antennas RTC and RTD, which the improved sky models do not address at all). Of particular interest is the effect that the improved model has on the differential gain-phases. As for the gain-amplitudes, we would expect them to differ by only an overall per-source scaling factor. Indeed, making the same  ||ΔE||  plot as in Fig. 7 confirms this – it is, to all intents and purposes, identical (and omitted here to save space), since the plotted amplitudes are renormalized by the per-source average  ||ΔE|| .

The phases, on the other hand, show a marked difference, since the formerly dominant effect – that of position offsets – has been taken out. The argΔE solutions themselves are shown in Fig. 14. Phase slopes are still very much in evidence, as can be seen in Fig. 15. Somewhat surprisingly, these slopes indicate that some residual position offsets remain, at a level of 15% to 20% of the original offsets (Fig. 16). This suggests that my procedure of fitting phase slopes to argΔE solutions, followed by fitting position offsets to the slopes, systematically underestimates the true position offsets. This is possibly an effect of the complex averaging implicit in having one ΔE solution per a relatively large solution interval (20 MHz by 30 min). If so, this could perhaps be incorporated as a multiplicative correction factor in the model update procedure. Further work is required to fully understand the effect.

The brighter sources B, C, and (to a lesser extent) D show clear second-order phase effects, both in the phase slopes, and in the phase solutions themselves. The temporal continuity in the phase slopes can be interpreted as a time-variable position offset. I can speculatively offer two explanations for such an offset:

• Unmodelled source structure (again!) For any given hour angle,an E-W array only sees an integrated cross-section through thesource in a given direction. If the source is slightly resolved withan asymmetric “hotspot”, the zero-order moment of each suchcross-section will be slightly different.

• Differential tropospheric or large-scale ionospheric refraction, including perhaps apparent change of baseline caused by refraction (the Anderson effect).

Another puzzling feature of the argΔE solutions in Fig. 14 are the significant and (to first degree) constant phase offsets of some sources (e.g. H, K, ae) on some antennas. The offsets are mostly (though not completely) consistent between the 2003 and 2006 observations. None of the explanations offered above are consistent with a constant phase offset! Could this be the phase component of the primary beam? There are too few sources in this reduction to infer any sort of directional dependence, but perhaps the “QMC Project” can provide more insights on this effect.

2.6. The lurking errors

The two calibrations (with the original and the improved sky models) described above have produced what appear to be identical final maps. This shows that the “flyswatter” can accommodate for significant errors in the sky model. On the other hand, the detailed structures in Fig. 15 suggest that (even in the very benign case of 21 cm WSRT observations!) moderately bright off-axis sources still require some form of DDE correction even if the model is perfect. If this is the case, then a legitimate question is: why worry about getting the sky model right, if we need to do ΔE solutions anyway, which will absorb any imperfections? (Besides the obvious caveats of the “flyswatter” discussed in Sect. 1.9, that is.)

The rather striking image of Fig. 17 shows that the final maps are not in fact identical, although the difference is buried in the noise. This image was produced by subtracting the original-model 8-band residual map from the improved-model map (2003 observation). Since the noise term in both maps is the same, subtraction reveals very faint structures that would normally be hidden in the noise. We’re beginning to see more limitations of the “flyswatter”. In the original reduction, apparent position offsets of the off-axis sources caused phase gradients in t,ν-space in the differential gain-phases. These were approximated by a stepwise ΔE solution (since I solved for only one ΔE term per 30 min, per entire band), which proved to be good enough to drive off-axis errors to a level below the thermal noise. Improving the model positions has effectively “flattened out” these gradients, reducing the error made by a stepwise approximation even further. Figure 17 demonstrates the improvement. The radial spokes correspond to “jumps” at the boundaries of the solution intervals, but the other structures (especially the half-circles) are rather more difficult to explain, and will have to be addressed in follow-up work.

The implications of this result is that any errors in the sky model (or uncorrected DDEs) will propagate into the selfcal solutions, and result in faint but highly coherent structures in the residual maps. We may think we are reaching the thermal noise, but in the process, we are producing “submerged” calibration artefacts at levels below the noise, where we can’t even see that something is still going wrong! This is of particular concern to ongoing work on detection of the Epoch of Reionization (EOR) signature, which relies on statistical analysis of residual images to find sub-noise artefacts of astrophysical origin (Harker et al. 2010; Morales et al. 2006). Such analysis will have to reckon with these lurking selfcal artefacts.

3. Conclusions

One of the biggest selling points of the RIME formalism is the flexibility it offers for describing observational effects. Unfortunately, to date only three software packages have exploited the power of the RIME (CASA, MeqTrees, and the LOFAR BBS system). Of these, only MeqTrees allows for truly arbitrary forms of the RIME. This paper has explored some practical applications of one such form of the RIME: a form that includes differential gain terms. I have demonstrated that the differential gain approach (the “flyswatter”) can be a powerful way of dealing with DDEs on a source-by-source basis. This has been used with WSRT data to produce artefact-free maps of 3C 147 at record dynamic ranges of well over a million-to-one. While the differential gain solutions themselves absorb inaccuracies in the sky model as well as the DDEs themselves, I have demonstrated that at least flux and positions corrections can be recovered, so iterative improvements to the sky model are possible.

The latter may also prove be necessary: I have demonstrated that even a perfect-looking map produced using differential gains contains a large number of selfcal artefacts hidden in the thermal noise, which can be significantly reduced by improving the sky models. These “invisible” artefacts have hitherto been ignored, but they should be of particular concern to projects relying on statistical signal extraction, such as the ongoing search for the EOR signature.

The nature of the remaining DDEs (as seen in the differential gain solutions) has not yet been adequately explained. Some of the amplitude effects are consistent with pointing error. The phase behaviour is even more difficult to understand, but may be due to unmodelled source structure. Further work is required on the subject.

I have shown that differential gain-phase solutions can be used to detect position shifts to within small fractions of the synthesized beam size. Offsets of less than 0.05″ (well under 0.01 of the PSF size!) have been reliably detected. There is a very clear indication of a systematic rotational offset of  ~45″ in the sky model generated by NEWSTAR, when interpreted using MeqTrees. This is may be due to a coordinate conversion error somewhere in the visibility data processing tool chain, and needs to be investigated further.

Finally, I should consider some wider implications of my results. All currently mooted schemes of DDE calibration for LOFAR (Nijboer & Noordam 2007), the MWA (Mitchell et al. 2008) and the ionosphere in general (Intema et al. 2009; Cotton et al. 2004) revolve around the use of “beacon sources” to probe the ionosphere and/or the primary beam. It is rather difficult to envisage a closed-loop scheme without beacons (how else would one sample a DDE?), so future telescopes such as the SKA will most likely need to use something very similar. Any such scheme predicates on there being a sufficient number of sufficiently bright in-beam beacons for any direction on the sky. This is not a problem at the LOFAR and MWA end of the spectrum, since the low-frequency sky is so much brighter, but it has been a bit of a worry for the higher frequencies, where FoVs are narrower and sources are fainter.

My 3C 147 results suggest that calibration beacons can be a lot fainter than previously thought. What has been established

is that for this particular configuration of the WSRT, sources as faint as 2 mJy can provide meaningful DDE solutions. This result can be scaled to future telescope designs by comparing their expected sensitivity with that of the 3C 147 observation.

---

Acknowledgments

I would like to thank a succession of managers for putting up with me all these years, and Jan Noordam for making this process considerably easier (especially for the managers), and for many other things besides. Ger de Bruyn has been more than generous with data, models and wisdom. Johan Hamaker started it all, and Wim Brouw has provided an avalanche of insights. Last but certainly not least, the rest of the MeqTrees team has been instrumental in making everything work.

References

All Figures