We present a GUI-based interactive Python program, VIMAP , which generates radio spectral index maps of active galactic nuclei (AGN) from Very Long Baseline Interferometry (VLBI) maps obtained at different frequencies. VIMAP is a handy tool for the spectral analysis of synchrotron emission from AGN jets, specifically of spectral index distributions, turn-over frequencies, and core-shifts. In general, the required accurate image alignment is difficult to achieve because of a loss of absolute spatial coordinate information during VLBI data reduction (self-calibration) and/or intrinsic variations of source structure as function of frequency. These issues are overcome by VIMAP which in turn is based on the two-dimensional cross-correlation algorithm of Croke & Gabuzda (2008). In this paper, we briefly review the problem of aligning VLBI AGN maps, describe the workflow of VIMAP , and present an analysis of archival VLBI maps of the active nucleus 3C 120. 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( S_v2 / S_v1 )/ log( v₂ / v₁ ), where Sv is the flux at frequency v and v ₁ and v ₂ are two different observing frequencies with v ₂ > v ₁ . 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. VIMAP is written in Python version 2.7.3 and uses the following well-known numerical and astronomical addon packages: These packages are available from the corresponding web sites. VIMAP itself, plus supplementary information material, is freely available on the Internet. 7 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 v ₂ > v ₁ needs to be convolved with the restoring beam of the map obtained at the frequency v ₁ , or has to be reconstructed from its δ function CLEAN components using the beam of the map at v ₁ . 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 v ₁ ; we recommend a size of 1/20 (or less) of the beam size. 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. 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 I_{v 1,2,ij} is the intensity at frequency v _1,2 and at spatial coordinate ( i, j ), is the mean intensity of a map at v _1,2 . The output of the cross-correlation analysis is a 2D map displaying the correlation coefficient r_xy 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 r_xy . If necessary, VIMAP users can set offsets manually. 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 S_v,ij is the root-meansquared thermal noise. RM S_v,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 ( α_{v 1,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 . 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 program 8 (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. 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 r_xy array ( Figure 2 ). The maximum correlation (≈0.99) is found for a relative shift of ≈0.16 mas (cf., Table 2 ). 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). 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).