1. INTRODUCTION

Active galactic nuclei (AGN) emit strong radio continuum synchrotron radiation and show spatial structure at scales ranging from sub-parsecs to kiloparsecs. Typical features are a core and outflows, especially jets, of usually complex structure. Multiple theoretical and observational studies indicate that the accretion of gas onto rotating supermassive black holes and acceleration of the accreted matter by magnetohydrodynamical processes is responsible for launching relativistic outflows and powering the observed highenergy synchrotron radiation (see, e.g., Boettcher et al. 2012 for a recent review).

Multi-frequency Very Long Baseline Interferometry (VLBI) observations are crucial for studying the physics and evolution of the outflows of AGN. Synchrotron radiation is characterized by its spectral index, α ≡ log(Sv2/Sv1)/ log(v2/v1), where Sv is the flux at frequency v and v1 and v2 are two different observing frequencies with v2 > v1. Measurements of the spatial distribution of α provide valuable physical information. Walker et al. (2000) revealed the spatial geometry of ionized gas associated with the accretion disk of 3C 84 by analyzing its continuum absorption spectra. Fromm et al. (2013) interpreted a radio flare in CTA 102 as interaction of a traveling shock with a stationary structure (potentially a recollimation shock). O’Sullivan & Gabuzda (2009) measured particle densities and magnetic field strengths in the parsec-scale jets of six active galaxies in order to test the classical model of AGN jets proposed by Blandford & Königl (1979). Analyses of the structural and spectral variability of radio jets complement single-dish studies of temporal AGN variability (e.g., Park & Trippe 2012, 2014; Kim & Trippe 2013) and of the interplay between accretion and jet formation (e.g., Allen et al. 2006; Trippe 2014).

In practice, obtaining spatially resolved spectral index information for AGN outflows is challenging mainly for two reasons. Firstly, an AGN core corresponds to a radio photosphere, e.g., an optical depth τv = 1 surface, for which the observed position is a function of frequency (a phenomenon known as core-shift effect; Lobanov 1998). Secondly, phase self-calibration, which is an important step in VLBI data reduction, removes absolute coordinate information (Kameno et al. 2003; Kadler et al. 2004). Both effects imply the need for a careful spatial alignment of radio maps obtained at different frequencies. Traditionally, two methods for aligning multi-frequency VLBI images have been used. The first one uses optically thin compact components, which are assumed to not change their positions as function of frequency, as reference points in a map (e.g., Kadler et al. 2004). The second method employs spatial correlations of optically thin extended jet structure for achieving alignment and/or for measuring the coreshift (Walker et al. 2000; Croke & Gabuzda 2008).

Currently, standard radio astronomical software packages, especially AIPS, CASA,1 and Difmap (Shepherd 1997) lack specialized tasks dedicated to image alignment. Motivated by this, we developed a highly interactive graphical user interface (GUI) based Python program, VIMAP, which adopts the two-dimensional cross-correlation scheme of Croke & Gabuzda (2008). VIMAP is open to the public and is capable of perform-ing map alignment and the generation of spectral index maps in a straightforward manner. In Section 2 we outline how VIMAP works. In Section 3 we present an analysis of archival VLBI maps of the AGN 3C 120 for demonstration.

2. THE VIMAP WORKFLOW

2.1. Prerequisites

VIMAP is written in Python version 2.7.3 and uses the following well-known numerical and astronomical addon packages: Numpy version 1.8.1 or higher2 Scipy version 0.14.0 or higher3 Matplotlib version 1.3.1 or higher4 PyFITS version 3.2.3 or higher5 wxPython version 2.8.12.1 or higher6

These packages are available from the corresponding web sites. VIMAP itself, plus supplementary information material, is freely available on the Internet.7

2.2. Preparation of Radio Maps

In order to begin the analysis, the sizes (in pixels) and pixel scales (angular extension of a pixel) of the two maps need to be made equal. Once this is achieved, the map obtained at frequency v2 > v1 needs to be convolved with the restoring beam of the map obtained at the frequency v1, or has to be reconstructed from its δ function CLEAN components using the beam of the map at v1. This is easily achieved using the standard AIPS and/or Difmap data reduction packages. The pixel scale should be much smaller than the synthesized beam size at v1; we recommend a size of 1/20 (or less) of the beam size.

2.3. Masking of Cores

In general, AGN cores need to be excluded from image alignment calculations because their positions may vary as function of frequency. In VIMAP, users can place an elliptical mask onto a core with variable (i) center position, (ii) semimajor/minor axis lengths, and (iii) orientation. Depending on how well the core can be identified and separated from the outflows in a given map, it may be necessary to iterate the alignment several times with different choices of masks.

2.4. Cross-Correlation

The relative shift between two images is estimated via a two-dimensional cross-correlation. The definition of the correlation coefficient rxy implemented in VIMAP is

(Croke & Gabuzda 2008) where Iv1,2,ij is the intensity at frequency v1,2 and at spatial coordinate (i, j), is the mean intensity of a map at v1,2.

The output of the cross-correlation analysis is a 2D map displaying the correlation coefficient rxy as function of relative offsets in right ascension (RA) and declination (DEC). The spatial shift between the images is given by the location of the maximum value of rxy. If necessary, VIMAP users can set offsets manually.

2.5. Spectral Index Maps

Eventually, VIMAP generates a spectral index map showing the index α (as defined in the introduction) as function of position, as well as a spectral index error map.

The spectral index error map is based on the assumption that the uncertainty of the flux density in each of the two radio maps used is a combination of (i) a systematic error coming from imperfect visibility amplitude calibration and (ii) random thermal noise from the observed source, sky, and instruments. This implies a total intensity error σv,ij , as function of frequency v and image location (i, j), of

Here δv 0.1 is the systematic amplitude calibration error on amplitude and RM Sv,ij is the root-meansquared thermal noise. RM Sv,ij can be measured in empty regions of an image. The factor δv may be known a priori for a given interferometer; if not, it can be estimated by measuring the flux from a compact flux calibrator and cross-comparison to results found by other observatories. Eventually, the uncertainty of the spec-tral index αij – i.e., the final spectral index error map Err(αv1,2,ij ) – is found via standard error propagation:

When generating the final spectral index map, VIMAP employs two cutoffs that exclude trivial or unphysical spectral index values: (i) a marginal intensity boundary, outside of which the source flux is assumed to be too weak for extracting spectral index information; and (ii) an upper limit on the error on αij .

3. APPLICATION

In the following, we provide an analysis of a VLBI data set with VIMAP. We used two Very Long Baseline Array (VLBA) maps of the AGN 3C 120 obtained in 2006 at 8.4 and 12.1 GHz, respectively. We obtained raw data from the archive of the MOJAVE program8 (Lister et al. 2009) and converted them into images using the modelfit task of Difmap. Observation details and intensity maps are provided in Table 1 and Figure 1, respectively.

Table 1Notes. Target: 3C 120. Observing date: May 24, 2006. Both maps have a size of 1024×1024 pixels and use a pixel scale of 0.05 mas/pixel. a: The 12-GHz map is convolved with the beam of the 8-GHz map before alignment. b: r.m.s. noise estimated from source-free regions in the CLEAN map.

Figure 1.VLBI maps of the active galaxy 3C 120 at 8 and 12 GHz. Contour levels decrease by factors of from the peak. Grey ellipses indicate beam sizes. Top: 8-GHz map. The peak intensity is 0.70 Jy/beam with the lowest contour corresponding to 0.9% of the peak value. Center: 12-GHz map. The peak intensity is 0.96 Jy/beam with the lowest contour level corresponding to 0.6% of the peak value. Bottom: Same as center panel, but using the restoring beam of the 8-GHz map. The peak intensity is then increased to 1.09 Jy/beam.

In each map, we covered the core (cf., Section 2.3) with an elliptical mask about three times the size of the restoring beam in order to ensure that optically thick regions do not affect the correlation. After this, VIMAP calculated the normalized 2D cross-correlation for the two input maps (cf., Equation 1) and saved the resultant rxy array (Figure 2). The maximum correlation (≈0.99) is found for a relative shift of ≈0.16 mas (cf., Table 2).

Figure 2.2D cross-correlation map; the color scales indicates the value of the correlation coefficient. Offsets are in units of pixels. The 12-GHz map is shifted relative to the 8-GHz map. Black contours correspond to correlation values of 0.7, 0.73, 0.76, 0.79, etc.

Table 2Shift of the 12-GHz map relative to the 8-GHz map.

After correcting for the relative offset, we computed the spectral index map (cf., Section 2.5). We adopted an intensity limit of 0.9% of the peak value of the 8-GHz map and an upper limit on the spectral index error of 0.5. Spectral index maps before and after image alignment are shown in Figure 3; the impact of the offset between the input maps is evident. Figure 4 provides the spectral index error map (top panel) and the evolution of the spectral index along a one-dimensional cut following the jet (center and bottom panels). At least three features should be noted: (i) Up to 2 mas from the core (≈ beam size; first red star from the core in the center panel of Figure 4), the α profile obtained after image alignment rapidly decreases, whereas the profile obtained without alignment is rather flat. (ii) About 3 mas from the core (second red star), the spectral index obtained without image alignment jumps up to ≈ 1. Without image alignment, this feature might be incorrectly interpreted as an overdensity in the jet. (iii) About the same holds for the spectral index from 7 mas to around 10 mas (third to fourth red star).

Figure 3.Maps of spectral index α (color scale) of 3C 120 before (top) and after (bottom) correcting for the offset between the maps. The impact of the offset is evident.

Figure 4.Top: Spectral index error map, assuming δv = 0.1 at both frequencies; the map was cut off at 0.5. Center: Same as the lower panel of Figure 3. The magenta line indicates a 1D cut through the α map along which the evolution of the index is analyzed. Red stars on the magenta line indicate locations of notable features shown in the bottom panel. Bottom: 1D spectral index profile along the cut indicated in the center panel, before and after correcting for the offset between the 8-GHz and 12-GHz maps.

4. SUMMARY

We report the development of a GUI-based interactive Python program, VIMAP, for generating radio spectral index maps of active galactic nuclei from VLBI radio maps obtained at two different frequencies. For the crucial task of accurate spatial alignment of the input maps we employ the 2D cross-correlation algorithm of Croke & Gabuzda (2008).9 As yet, such a tool has not been available in the frame of standard radio interferometric data processing packages like AIPS or CASA. Our approach complements the method of image alignment via cross-identification of jet model components (Kameno et al. 2003; Kadler et al. 2004). As noted explicitly by Kadler et al. (2004) (in their Section 4), the model component approach depends on a reliable crossidentification of optically thin jet components which is not always possible. VIMAP offers a higher degree of flexibility because it does not rely on such crossidentifications and, in fact, does not even require modelling but may be applied to deconvolved (CLEANed) maps immediately. A detailed comparison of the two approaches is provided by Fromm et al. (2013) (their Section 2.1).

A key feature of VIMAP is its interactive GUI-based approach depending only on few standard, freely available, Python packages. Accordingly, VIMAP provides users with a new tool capable of handling large amounts of multi-frequency VLBI data in a straightforwardmanner, especially in view of users of the Korean VLBI Network (KVN) (Lee et al. 2014) which is capable of si-multaneous multi-frequency observations, and the KVN and VERA Array (KaVA) (Niinuma et al. 2014).