3. DATA REDUCTIONS

3.1. Calibrations

All calibrations and data reductions were performed using standard routines in the Common Astronomy Software Applications package, CASA ¹, using versions 4.2 and earlier. Paper III describes in detail the process we have been following for calibrations and reductions. A summary follows here:

In order to reduce ringing across the band from strong RFI, it was necessary to Hanning smooth our data. Antenna-based delay calibrations have to be applied prior to this smoothing, which was done in an initial stage of the calibrations.

The absolute flux density scale of the data was determined by observing one of two standard primary calibrators; 3C286 or 33C48 (the latter for three of our early right ascension galaxies, NGC 891, NGC 660 and part of the NGC 2613 observations), and applying the Perley-Butler 2010 flux density scale. We used the primary calibrator to determine the bandpass corrections as well. See Table 1 for a list of the calibrators which were used for each galaxy in the sample.

3.1.1. Polarization calibration

For the polarization calibrations, the primary calibrator was used as a polarization angle calibrator. One of the NRAO standard zero polarization calibrators was used to calibrate instrumental polarization leakage. We corrected for cross-hand delays, leakage terms and polarization position angles. Where in some cases the polarization leakage calibrator could not be used (due to loss of scan during observation or being more heavily flagged than the other calibrators), we instead utilized the secondary calibrator where the range in parallactic angle was sufficient, i.e. more than 60 degrees. This was for example done for the galaxies in a C-band scheduling block (ID: 4809749), for which the polarization leakage calibrator scan was lost. Q and U fluxes were calculated in the standard way – see Paper II and references therein. Our tests have shown no difference in the result whether the secondary calibrator or polarization leakage calibrator was used.

3.2. Flagging

For all scans, calibrators and galaxy scans alike, we inspected the uv-data by eye. All bad data, whether from RFI or instrumental effects, were manually flagged. A number of tests were carried out using automatic CASA flagging routines as well as the CASA pipeline; however, we found significantly improved results by continuing with manual flagging.

Generally, fewer antennas could be used for our L-band observations than for C-band, since not all VLA antennas had yet been equipped with L-band receivers at the time of observation. Additionally, there is typically more RFI at L-band.

Calibrations and flagging were carried out iteratively, until a well-calibrated data set could be obtained. The final calibrated galaxy scans were then separated from the calibrator scans, in order to be imaged in Stokes I, Q and U.

3.3. Imaging

We have produced Stokes I, Q and U images, as well as spectral index maps (α, where I_ν∝ ν^α) and uncertainty maps, Δα (see Sect. 3.4) from the calibrated data.

The clean algorithm in CASA, which we use for our wide-field synthesis deconvolution, allows for multi-scale multi-frequency synthesis, ms-mfs (see Rau & Cornwell 2011, for a full description). We set the multi-scale feature to look for flux components at a variety of spatial scales, to account for the extended emission of our galaxies. Typically our scales ranged from zero (’classic clean’) up to approximately five times the synthesized beam, but sometimes required adjustments. In L-band, and for the uv-tapered weighting in particular, fewer scales were at times necessary in order to remove artifacts for point-like sources; for example a classic clean was set for NGC 2683, NGC 2992, NGC 4157 and NGC 4845.

Additionally, the clean algorithm enabled simultaneous fitting of a spectral index across the band width (’in-band spectral index’) via a simple power law. In order for this to take place, we set the number of Taylor coefficients used to model the sky frequency dependence,nterms, to 2 for all CHANG-ES galaxies. We defer the details of spectral information to Section 3.4, and will describe the spatial fitting, as implemented for the CHANG-ES project, in this section.

3.3.1. Weightings

For each galaxy, the Briggs robust weighting (Briggs 1995) was used for both C and L band data. Additionally, we imaged a uv-tapered version of this weighting, which we introduced in order to emphasize the broadscale structures ². We achieved an increase in beam size of approximately 80% in C-band and 40% in L-band, by typically applying a 6 kλ Gaussian taper in C-band and a 2.5 kλ Gaussian taper in L-band to the uv-distribution. Image products of both these weightings are included in Data Release 1 and are referred to as robust 0 and uv-tapered weightings, respectively.

In C-band, images sufficiently clear of artifacts with uv-tapered weighting could not be achieved for NGC 4438 and NGC 4845 due to particularly strong contaminating sources in the field. In L-band, the uv-tapered weighting was either challenging to achieve (for the same reason) or unnecessary (the sources were imaged as point sources already at robust 0) for six of our galaxies; NGC 3448, NGC 4192, NGC 4388, NGC 4438, NGC 4845 and UGC 10288, and the attempt was discarded for these. Hereafter any images shown will be of robust 0 weighting, unless otherwise specified.

3.3.2. Stokes I imaging

Each imaging run was cleaned down to a flux density level of 2.5-3 σ. With few exceptions, we imaged the entire field (i.e. without specifying regions), since a lower rms could be obtained with this strategy. Additionally, the vast majority of our data do not suffer from significant artefact patterns with spurious sources outside the source region, which could jeopardize this strategy. When needed, one to two self calibrations were performed in order to deal with remaining artifacts.

Our self calibration strategy was based on iteratively refining the calibration table, by self calibrating at successively deeper thresholds and thus improving the input model, which in turn acts on the original visibilities. At each step, care was taken to a) check the model to be used for the self calibration that it did not include any artifacts and b) check that the peak intensity of the source would not decline (such a scenario was aided by opting for a phase-only self calibration at the first iteration, or by applying the self calibration at a shallower threshold).

In some cases, our cleaning strategy had to be adapted to the circumstances of the data (the quality of a dataset, strength of emission, contamination in the field, etc). This particularly applies to data with high dynamic ranges (HDR), due to galaxies with strong, often point-like centre sources or strong field sources (for example NGC 660, NGC 4594, Virgo galaxies). Our strategies to deal with these situations include one or more consecutive runs of self calibrations, specifying regions, attempts at peeling for strong field galaxies (see 3.3.6) and clean parameter adjustments (including multi-scale adjustments). Generally, the best one can achieve in terms of HDR, is a signal over noise of the order of 10000:1 (S. Bhatnagar, priv.comm.). Our dynamic ranges are listed in column 8 of Tables 4 and 5. Additionally, in cases of strong centre sources, choosing a smaller cell size will characterize the steep beam better and render a cleaner result (e.g. NGC 4594 and NGC 4845, both in C-band).

3.3.3. Polarization imaging

Stokes Q and U images were produced similarly to I. These maps were used to create polarization intensity and polarization angle maps. The polarization intensity maps were corrected for bias caused by the noise level ³ and the polarization angle maps are cut off at the 3σ level. We note that if the fraction of polarization intensity over Stokes I intensity is less than 0.5%, the polarization is not believable since the calibration cannot be guaranteed below that level (S. Myers, priv. comm.).

A uniform Faraday screen resulting from the ionosphere is corrected for in the standard polarization data reduction. However, if there is differential Faraday rotation, for example between the location of the primary calibrator and the source, then this is not corrected for ⁴. This effect is negligible at C-band. There could, however, be minor effects at L-band, which will be investigated in a subsequent polarization paper.

Our vector maps (panels (d), (e) and (f) of the figures in the appendix) are also ’apparent B-vectors’ which simply are E-vectors rotated by 90 degrees, because internal Faraday rotation in the galaxy has not been corrected for.

3.3.4. Primary beam corrections

Wideband primary beam corrections were carried out with the CASA task widebandpbcor after cleaning. This task accesses information about the beam from the visibility data, and calculates a known NRAO-supplied model of the primary beam (PB) at each frequency channel in each band. A PB cube is formed, and Taylor terms are found in the same way as was done for the data set when imaging. Each Taylor coefficient image (two in our case, see Sect. 3.3 and 3.4.1) is then corrected.

This task also corrects the spectral index image above a given threshold (we use 5σ), as well as corrects a map of the formal error of the spectral index (see Sec. 3.4.2). Subsequently, we also applied primary beam corrections to the polarization intensity maps, by dividing the image with the beam created in the widebandpbcor task.

Note that in this paper, any displayed images are not PB corrected, but any measurements are made from the PB-corrected versions of the maps, unless otherwise indicated. For total intensity images, use of a model for the PB rather than an observation-specific PB introduces negligible error over the scales of our galaxies. However, the effect on the spectral index maps must be investigated in more detail, which we do in Sec. 3.4.2.

3.3.5. Mosaicking two pointings

For the eight large galaxies which were observed in two pointings at C-band (see Table 2), each individual pointing was calibrated and imaged separately. Although it would be preferrable to mosaic the pointings onto the same grid while still in the uv-plane, this feature is not yet implemented in CASA at the time, for the clean settings that we required (i.e. nterms = 2). Instead, we combined the images from the two pointings, forming a weighted (by the primary beam) average in the overlapping region. The images were combined in the overlapping regions out to 10% of the primary beam. This rather inclusive choice of cut-off worked well for the D configuration data, but will likely be adjusted to 50% for the B and C configuration data in the future.

This strategy was applied to both the non-PB-corrected images, as well as the two primary beam pointings. The combined Stokes I image was divided by the combined primary beam in order to do the PB correction.

The two-pointing spectral index maps, which had already been individually PB-corrected, were similarly combined, as well as their accompanying error maps.

We assessed whether the two pointings had differences within the overlap regions, such as in noise levels. In such cases additional weightings, based on the noise level, would have to be considered. This was however not the case for any of the CHANG-ES galaxies, where separate pointings generally had been observed in turns in the same scheduling blocks, or in similar conditions.

3.3.6. Peeling of strong field sources

For a couple of cases, for example NGC 660 and NGC 4388, bright field sources interfere strongly in the cleaning process, causing a less than ideal imaging result fraught with severe artifacts (L-band only). An attempt was made at peeling these interfering sources in those examples, that is, to remove or decrease the intensity of the interfering source so that its impact is lessened in the final deconvolution. The process includes self calibrations with the interfering source centred, and subtracting the source from the uv-data using the CASA task uvsub, in several steps. For more information on peeling, see Pizzo & de Bruyn (2009), as well as the PhD thesis by Adebahr (2013). Table 3 lists the sources we peeled from NGC 660 and NGC 4438.

3.4. In-band Spectral Indices

3.4.1. Formation of Spectral Index Maps

As indicated in Sect. 3.3, the ms-mfs algorithm allows for spectral fitting across the wide bands that have been used in the CHANG-ES survey. Thus, a single observation over a single wave-band permits the determination of spectral indices, α, as a function of position within the galaxy. As described in Rau & Cornwell (2011) (see also Paper II), a specific intensity, I_ν, at some frequency, ν, within the band, can be fit with a function of the form I_ν = I_ν0 ν^{α +β log(ν/ν₀)}, where β is a curvature term and ν₀ is a reference frequency (our central frequency, ν₀, given in Tables 4 and 5). As implemented, such a fitting function is expanded in a Taylor series about ν₀, such that the first Taylor term (the map, TT0) corresponds to a map of specific intensities at ν₀ and the second Taylor term (the map TT1) contains information on α such that α = TT1 / TT0. The third Taylor term (TT2) allows for the recovery of the spectral curvature, β, where β = TT2 / TT0 − α(α − 1)/2.

In principle, any number of Taylor terms could enter into such an expansion with increasing numbers of terms improving the spectral fit (see Rau & Cornwell 2011, for examples). In practice, however, the effective number of terms is limited by the signal-to-noise (S/N) ratio and calibration precision of the data. For the CHANG-ES galaxies, tests that we have run gave poorer results with 3 terms than with 2; that is, the rms noise in the images is higher for such fits, when cleaning to the same threshold.

The curvature maps also show large and sometimes unrealistic variations, point-to-point. While globally averaged values of β can still be useful (see Paper II for an example), we have adopted two Taylor terms for the CHANG-ES project as indicated earlier. A two-term Taylor fit then takes the simple form I_ν = I_ν0 + I_ν0 α ((ν − ν₀) / ν₀). A map of random errors describing the accuracy of the fit, Δα is also formed ⁵.

The beauty of the spectral fitting is that it is carried out by flux component which means that a uniform resolution is achieved over the band, with a restoring beam applied to the TT0 and TT1 maps at the end of the fitting process. Also, the S/N of the entire band applies to each flux component. If one were to form a spectral index in the classical way within a band, one would need to break up the band into its various channels (or groups of channels), smooth each channel to the poorest resolution of the lowest band frequency, and accept the very high noise of these channels as a cut-off when forming α maps.

Clearly, the method, as implemented in CASA, is far superior for any individual band ⁶. Nevertheless, there are limitations, as we outline below.

3.4.2. Post-Imaging Corrections to the Spectral Indices

Because the primary beam (PB) varies with frequency, it imposes its own spectral index, α_PB, onto the α maps. This effect is significant; for example, at the 50% level at L-band, α_PB = −1.6. We use the CASA task widebandpbcor, described in Section 3.3.4, to correct for α_PB, and at this point, we cut off the α maps at 5σ, where σ is the rms of the image map (i.e. the TT0, or I_ν0 map).

Once PB-corrected α maps are formed, they can show large variations around their peripheries and the corresponding error maps, Δα, show correspondingly large errors. In addition to the 5σ cut originally applied, we also cut off all α maps wherever the corresponding error map exceeds a value of 1.0. The choice of this cutoff is arbitrary; however, various tests showed that more aggressive cuts (e.g. where Δα > 0.75) resulted in undesirable ‘holes’ in some α maps. Note that this step does not change the PB correction of the α maps. It simply discards data points that are known to have large errors.

Both α and Δα maps are computed for each map pixel (cell). However, these pixels are not independent. For example, at L-band, there may typically be 64 cells/beam and at C-band, 50 cells/beam (both uv-tapered). Consequently, the ‘per pixel’ errors are much greater than the ‘per beam’ errors and the α maps, as produced, show variations pixel-to-pixel which are larger than beam-averaged variations would be. In addition, maps of Δα show artifacts, i.e. regions of low Δα that have a ‘thread-like’ appearance through the map. These artifacts depend on the choice of multi-scales that are used during imaging and shift position if the map is remade with different scales. If the adopted scales are changed, these ’threads’ also shift position.

To ameliorate the above issues, we have ‘smoothed/ averaged’ our maps of α and Δα over the size of a beam. This is accomplished by artificially inserting into the map header a Gaussian beam size whose area (Gaussian-weighted) is equivalent to the pixel area and then convolving the map to the correct clean beam size ⁷. The result does not change the spatial resolution of the images but recognizes the non-independence of the pixels and strongly minimizes pixel-to-pixel variations. For example, the rms of an α map could decrease by as much as a factor of 2 after such a convolution. These are the final maps that are shown in the appendix, panels h), i), k) and l).

In summary, our final spectral index maps and all calculations performed on them apply to α after

1. PB-correction,
2. a 5σ cutoff,
3. a cutoff wherever Δα > 1.0, and
4. convolving with a Gaussian to smooth over pixel-to-pixel variations within any clean beam.

3.4.2. Uncertainties in Spectral Index Maps

As indicated above, α is determined via fitting flux components over the band during the imaging process. Fitting errors should, for the most part, be accounted for by the Δα maps, but those maps consider random errors only. Therefore, a variety of tests have been carried out on the in-band spectral index maps to help us understand realistic uncertainties, and we outline them here. Largely qualitative descriptions are provided, with a quantitative summary at the end. In the following, we refer to the final α maps, as described above, unless otherwise indicated.

Edge-effects:

In spite of the cutoffs and smoothing outlined in the previous section, values in the spectral index maps (as well as their errors) tend to become extreme at the edges (for example, see the L-band panels (i) and (l) of NGC 3735 in the appendix).

Primary Beam Correction:

The NRAO-supplied PB model cannot take into account variations in the PB that are specific to any given observation. The true PB will change with flexure in the antennas and therefore with position on the sky. It will also rotate on the sky as the source is tracked and therefore any non-axisymmetry in the PB will also introduce error ⁸. As a result, PB correction errors, which increase with distance from the pointing centre, must be understood. Bhatnagar et al. (2013) have shown that errors in the α maps are negligible out to the half-power point but increase significantly beyond that distance (see their Fig. 4, right). For L-band, our largest α map (NGC 4565) extends to a radius of 6.8 arcmin which is well within the HPBW of the L-band PB (radius of 15 arcmin, or FWHM of 30 arcmin). At C-band, however, the PB radius is only 3.75 arcmin and consequently α maps with radii larger than this will show such errors in their outer regions.

We have investigated such potential PB errors in two ways:

The first is to examine α maps using cases in which the observations of a given galaxy were carried out in two different observing sessions, such as NGC 660 (see Table 1). As such, it is possible to form two different α maps for the same galaxy and same sky pointing but observed on different days. This has allowed us to compare the resulting α maps and to compare the variations between those days with the random errors as given in the Δα maps. We tested both low and high dynamic range cases and compared the two results quantitatively.

The second is to examine our C-band observations of large galaxies where two offset pointings were observed (see Sec. 2.2.1). A comparison of α made for the individual pointings allowed us to investigate α errors at larger distances from the PB centre.

As an example, Fig. 1 (panel a) shows a map of the absolute pixel-by-pixel difference between the PB-corrected spectral index measurements in the two C-band pointings of NGC 891. An increase of the differences with distance from the map centre is clearly apparent. We formed the spectral index error (shown in panel b) by computing the PB-weighted deviation from the average spectral index of the two pointings ⁹.

Figure 1. a) Absolute difference between the PB-corrected spectral index measurements in the two C-band pointings of NGC 891. There is an increase in difference with distance from the midpoint between the pointing centres, even in regions with high S/N. b) Error of the PB-weighted average of the PB-corrected spectral index measurements in the two pointings. The red contour is placed at an error of 0.1. c) Spectral index error map determined by the ms-mfs clean algorithm (same as panel (k) of the appendix figure of NGC 891, but without the cutoff at Δα > 1.0). Contours represent an error of 0.1 in this map (white) and in the map shown in panel b) (red). d) Ratio of the two different error maps (panel b) divided by panel c)). The ratio is mostly close to unity, except for those parts of the disk lying outside the 0.7 PB level of either pointing. In each panel, the two PB circular contours for each pointing are placed at the 0.5 and 0.7 level.

This error map is for the most part comparable to the error calculated by the ms-mfs algorithm (panel c)), as the map of the ratio of these two errors (panel d)) illustrates. A major exception to this is the disk of the galaxy, but only in those parts that are located outside the 70% PB level (i.e. where the PB gain is less than 0.7) of either pointing. The average error ratio in these narrow regions around the mid-plane is ∼ 5 (with individual pixel values up to ∼ 15), i.e. here the error originating from the pointing differences is on average ∼ 400% higher than the error resulting from the ms-mfs spectral fitting (for comparison, the average error ratio of NGC 4565 is 3.3). Such high ratios are primarily a consequence of the small ms-mfs-based errors in regions of high signal-to-noise, yet the increase of the pointing differences with distance from the pointing centres is significant. In particular, beyond the half-power point the mid-plane errors in panel b) of Fig. 1 increase up to ∼ 0.5, whereas in panel c) they remain below 0.1 throughout the disk, as shown by the displayed contours. While the errors in panel b) are in rough agreement with Bhatnagar et al. (2013) in the sense that they do not exceed 0.1 out to approximately the half-power point (not considering the above-mentioned edge effects), the error ratio map suggests that the simple PB model used by the widebandpbcor task already shows significant inaccuracies at the 70% level.

Resolved out structures:

Although the method of calculating the spectral index ensures that the resolution is common across the band, if there are structures that are resolved out at one end of the band compared to the other, this is not accounted for in the spectral index maps. For example, there may be structures resolved out at the high frequency end of the band, but not at the low frequency end, which would steepen the spectral index.

Comparison with Classical Spectral Index:

A comparison was made between the in-band spectral index and a classically formed spectral index from a given band. The classically formed α maps were made from the spectral index maps after the PB-correction was applied but not after applying further processing (i.e. PB-correction and a 5σ cutoff were applied, but not a cutoff based on the error map, nor smoothing to the size of the beam, as described in Section 3.4.2); they were formed by imaging the lowest end of the band and the highest end of the band, smoothing to a common spatial resolution and then combining in the classical way. The results were in agreement within errors.

Two-pointing PB cutoffs: For ‘two-pointing’ galaxies (large galaxies at C-band) the maps, as indicated in Section 3.3.5, were combined over regions in which the value of the primary beam exceeded 0.1 (10% of the peak). When mosaicking is carried out, each point is weighted by the PB such that points that are farther from the beam centre are weighted down and therefore the difference in total intensity maps, whether one cuts off points where the PB is > 10% or where the PB is > 50%, is entirely negligible. However, since α maps are known to increase quite strongly below the 50% level, we remade our mosaicked α maps with a 50% cutoff to compare it to the 10% cut-off values used in this paper. For these results, we found only minor differences (for example, if the two maps are subtracted, the rms in the regions between the 50% and 10% PB cutoffs is < 0.1).

In summary, our final α maps should present realistic results with the following cautions:

a) Extreme values around the edges are artifacts and should be ignored.
b) Measurements of α should not be quoted ‘per pixel’ but rather averaged over a beam when carrying out evaluations or comparisons.
c) The calculated Δα maps typically underestimate the true errors by ≈ 20% for galaxies of small angular size ¹⁰.
d) The largest source of uncertainty in Δα relates to the PB model correction away from the pointing centre, which most strongly affects our largest galaxies in C-band. Spectral index errors increase significantly with distance from the pointing centre, such that beyond the 70% PB level of either pointing (i.e. at distances greater than 2' 25" from the respective pointing centre at C-band) the Δα maps typically underestimate the true errors in the disk by a factor ∼5 (see Fig 1).

¹ http://casa.nrao.edu Back.

² For an example of a different weighting, see Paper II for which a naturally weighted map of NGC 4631 was included. Back.

³ Note that the maps shown in panels (e) and (f) of the images in the appendix are the noise biased maps. Peak polarization values have however been measured from the corrected maps. Back.

⁴ The software is not yet available in CASA. Back.

⁵ See Eqn. 39 of Rau & Cornwell (2011). Back.

⁶ For an example showing a result when more than one band is used, see Paper III. Back.

⁷ The CASA task imsmooth is used. Back.

⁸ Recently, it has become possible to correct for narrow-band time-varying PB effects during the imaging process using the ‘A-Projection’ algorithm (Bhatnagar et al. 2008, Bhatnagar et al. 2013). However, an algorithm applicable to CHANG-ES data (i.e. for wide bands) is not yet available or practical. Back.

⁹ √(w₁(α₁ − α_avg)² + w₂(α₂ − α_avg)²) / (w₁ + w₂), where w₁ = beam₁ / beam₂ and w₂ = beam₂ / beam₁ Back.

¹⁰ Estimated from results where galaxies have been observed in two separate observations. Back.