© R. Angioni et al. 2019

1. Introduction

Radio-loud active galactic nuclei (AGN) host symmetric, highly relativistic and well-collimated jets, ejected from the vicinity of the central supermassive black hole (SMBH). These jets produce bright non-thermal radiation across the whole electromagnetic spectrum. AGN dominate the γ-ray sky, as revealed by the Fermi Large Area Telescope (LAT) (Acero et al. 2015). The largest class of identified γ-ray sources is blazars, that is, radio-loud AGN with jets oriented at small angles (θ <  10°) to the observer’s line of sight (Ackermann et al. 2015). The alignment and the relativistic bulk speed of the jet imply that its radiation is beamed in the direction of motion and boosted via relativistic Doppler effects, which make these sources extremely bright and shortens their variability time scales.

After decades of studies, there are still important unanswered questions on jet physics and their high-energy emission. For example, the angular resolution of the best GeV γ-ray instrument, the Fermi-LAT, can only reach down to θ ∼ 0.1°. Therefore it is rarely possible to associate the γ-ray emission to a specific morphological component of the AGN, such as for example, the radio core, jet, or diffuse extended lobes.

According to the most popular models, the unbeamed parent population of blazars is made by classic radio galaxies (see for instance, Urry & Padovani 1995; Tadhunter 2016). In these sources, the jet is oriented at larger angles to our line of sight, therefore its radiation is much less affected by Doppler boosting. For this reason, they are unfavored γ-ray emitters, making up only about 1–2% of identified sources in the Fermi-LAT catalogs, with a total of approximately 20 objects (Ackermann et al. 2015). The faint fluxes and small sample size make understanding γ-ray radio galaxies as a population challenging. Nonetheless, studying these sources provides us with a view of jets which is less masked by orientation-based relativistic effects. This allows us to better reconstruct the intrinsic properties (Abdo et al. 2010c).

The combination of γ-ray data with high-resolution radio observations can be a powerful tool in investigating these questions. Very Long Baseline Interferometry (VLBI) radio observations allow us to achieve an angular resolution which is orders of magnitude better than that of any γ-ray instrument, down to milliarcsecond scales. Since AGN are variable sources, it is possible to identify the γ-ray emission region on VLBI scales by looking for correlated variability in radio morphology or flux, and high-energy emission (see for example, Casadio et al. 2015). Additionally, multi-epoch VLBI observations provide the only direct measure of relativistic jet motion, allowing us to derive relevant jet parameters such as the apparent speed and jet orientation angle. This motivation led to the development of long-term radio VLBI programs monitoring a large number of AGN (see for instance, MOJAVE, Lister et al. 2009).

TANAMI (Tracking Active galactic Nuclei with Austral Milliarcsecond Interferometry) conducts the only large monitoring program of this kind for sources south of −30° declination, with a sample of about 130 AGN. First-epoch 8.4 GHz images were presented for an initial sample of 43 sources by Ojha et al. (2010). The sample was regularly expanded adding new γ-ray detected sources or otherwise interesting new objects. First-epoch images of the first-ever high-resolution observations, for most sources, were recently published (Müller et al. 2018). The radio VLBI monitoring is complemented by excellent multi-wavelength coverage, including near-infrared, optical, UV, X-ray and γ-ray data. These provide a quasi-simultaneous broadband view of the sources, as is required for detailed studies of variable sources such as AGN (see for example, Krauß et al. 2016). An overview of the multi-wavelength program and selected TANAMI results can be found in Kadler & Ojha (2015).

In this paper, we focus on the radio VLBI properties of the γ-ray detected TANAMI radio galaxies, that is Pictor A (a.k.a. 0518−458), PKS 0625−35 (a.k.a. 0625−354, 0625−35), Centaurus B (a.k.a. Cen B, 1343−601) and PKS 1718−649 (a.k.a. 1718−649, NGC 6328), and the peculiar AGN PKS 0521−36¹ (a.k.a. PKS 0521−365, 0521−365). We produce high-resolution images at 8.4 GHz for all available epochs between 2007 and 2013. We perform a kinematic analysis by fitting Gaussian components to the jet and tracking their motion across the epochs. To complement this, we produce γ-ray light curves and spectra from Fermi-LAT data, investigating a possible connection between the two bands.

In a subsequent paper, we will focus on the γ-ray faint TANAMI radio galaxies, and present the radio kinematic results, along with updated Fermi-LAT flux upper limits, to investigate if and how the two subsamples differ from each other.

The paper is organized as follows. In Sect. 2 we review the TANAMI sample definition, investigate the completeness of the TANAMI radio galaxy sample, and give an overview of the source sample. In Sect. 3 we illustrate the data reduction procedures for the radio VLBI (Sect. 3.1) and the γ-ray data (Sect. 3.2). In Sect. 4 we present the multi-frequency VLBI imaging results (Sect. 4.1), the VLBI kinematics results (Sect. 4.2), and the results of the Fermi-LAT analysis (Sect. 4.3). Finally, in Sect. 5 we discuss the scientific implications of the radio and γ-ray data, and draw our conclusions in Sect. 6. All the results are presented on a source-by-source basis. Throughout the paper we assume a cosmology with H₀ = 73 km s⁻¹ Mpc⁻¹, Ω_m = 0.27, Ω_Λ = 0.73 (Komatsu et al. 2011), and the radio spectral indices refer to the convention S ∝ ν^+α.

2. The TANAMI radio galaxy sample

The TANAMI sample was defined starting from two subsamples of sources south of δ = −30°, that is a radio-selected sample and a γ-ray selected sample. The radio-selected subsample was based on the catalog of Stickel et al. (1994)², with a flux density cut at S_5 GHz >  2 Jy and a spectral index cut at α >  −0.5 (S ∝ ν^+α) between 2.7 and 5 GHz. The γ-ray selected subsample included all identified EGRET blazars in the given declination range. The sample also includes a few additional sources of interest which did not satisfy this selection, and were added manually to the monitoring program.

We investigated whether the sample of radio galaxies in TANAMI, including sixteen sources is representative of a complete sample of radio galaxies in our declination range. To do this, we cross-matched the Véron–Véron 13th edition AGN catalog (Véron-Cetty & Véron 2010) with the Parkes radio catalog (Wright & Otrupcek 1990). While the Véron–Véron catalog is not complete in a statistical sense, it is a comprehensive compendium of known AGN. The Parkes catalog on the other hand provides extensive information on the radio properties of the sources.

Starting from this cross-matched catalog, we first applied the declination cut for TANAMI, that is, δ <  30°. Since we are only interested in radio galaxies for the purpose of this study, we had to clean the sample by excluding known BL Lacs and QSOs. Physically, the only criterion distinguishing radio galaxies from blazars should be the jet viewing angle. However, it can be challenging to obtain this information, and estimates often suffer from a large uncertainty. We therefore made this selection based on the optical spectral classification provided in the Véron–Véron catalog.

This left us with a total of 83 sources. Twelve out of sixteen TANAMI radio galaxies are included in this sample³. We looked for existing VLBI measurements for each of these 83 sources individually on the NASA/IPAC Extragalactic Database (NED), and found existing data for six non-TANAMI radio galaxies, from the Very Long Baseline Array (VLBA) at 8.6 GHz (Petrov et al. 2005) or the Long Baseline Array (LBA) at 8.4 GHz or 4.8 GHz (Fey et al. 2004; Hancock et al. 2009). Therefore, in total we found 22 southern radio galaxies with compact radio emission on VLBI scales. The TANAMI radio galaxy sample, with 16 sources, can be considered representative of the subpopulation of radio galaxies with a bright compact core.

Our sample includes well-known sources for which TANAMI provides the highest-resolution data available. The full list of TANAMI radio galaxies can be found in Table 1. Several radio galaxy subclasses are present, from classic FR I (such as Centaurus A) and FR II (such as Pictor A) to young radio sources (such as PKS 1718−649) and peculiar or misclassified AGN (such as PKS 0521−36). In this work we focus on γ-ray detected sources. Only half of the sample has been associated to a Fermi-LAT source so far. This is not surprising, given that radio galaxies are faint γ-ray emitters. Since TANAMI is, to date, the largest VLBI program monitoring southern-hemisphere sources at milliarcsecond resolution, including a considerable number of radio galaxies, it provides the first data set suitable for kinematic studies on these scales. The most notable member of the TANAMI radio galaxy sample is the closest radio-loud AGN, Centaurus A, which has been studied extensively using TANAMI data, revealing the complex dynamics of the jet’s inner parsec. A recent kinematic analysis indicates downstream jet acceleration (Müller et al. 2014b), and the spectral index map between 8.4 GHz and 22.3 GHz suggests a transverse structure that can be explained within the spine-sheath scenario (Ghisellini et al. 2005).

The peculiar source PMN J1603−4904 has also been extensively studied within the TANAMI Collaboration before (Müller et al. 2014a, 2015, 2016; Goldoni et al. 2016; Krauß et al. 2018). Therefore these two sources are not included in this study.

3. Observations and data reduction

3.1. Radio data

TANAMI monitors the sources in its sample with a cadence of about two observations per year since 2007, at 8.4 GHz and 22.3 GHz. The VLBI array consists of the Australian LBA, supported by associated antennas in South Africa, New Zealand, Antarctica and Chile (a full list of the participating antennas is reported in Table 2).

Details on the array configuration for each observation are provided in Tables A.1–A.5. The data were calibrated using the National Radio Astronomy Observatory’s Astronomical Image Processing System (AIPS; Greisen 1998) and imaged using the CLEAN algorithm (Högbom 1974) as implemented in the Difmap package (Shepherd et al. 1994) following the procedures described in Ojha et al. (2010).

Many more antennas are available at 8.4 GHz (than at 22.3 GHz), yielding better angular resolution in this band, despite the lower frequency. Observations in this band are also more frequent. Hence, 8.4 GHz data are used for the kinematic analysis. We have used selected epochs at 22.3 GHz, simultaneous with the 8.4 GHz ones, to produce spectral index maps of our targets. The spectral information is crucial in order to identify the VLBI core component, which usually presents a flat spectral index, and as an input for estimating the jet viewing angle. Since the absolute position information is lost in the imaging process due to the phase self-calibration, it is necessary to properly align the maps before computing the spectral index. This is done through a 2D cross-correlation procedure (Fromm et al. 2013), referenced on an optically thin region of the jet, whose position should not vary with frequency. After aligning the maps and convolving them with the same beam (typically the one of the map with lowest resolution), the spectral index is calculated for each point as

(1)

A full estimate of the uncertainty of spectral index maps requires taking into account frequency-dependent flux density errors, pixel-by-pixel signal-to-noise ratios, and an analysis of the error in the shift between the two maps through Monte Carlo simulations (see Appendix D in Fromm et al. 2013). Since we did not perform a detailed study and in-depth interpretation of our spectral index maps, such an extensive error analysis is beyond the scope of this paper.

In order to study the evolution of the jet, we fit the self-calibrated visibilities with circular Gaussian components using the Modelfit task in Difmap. We then cross identified them in the different epochs by selecting a component in the first epoch, and searching for the closest component in the following epochs. This selection was then corroborated by visual inspection of the maps and the component properties, such as the evolution of their flux density. We then fit the motion of the components which were robustly detected in at least 5 epochs, separately in RA and Dec, to derive the two components of the velocity vector, following Lister et al. (2009):

(2)

where μ_x and μ_y are the angular speeds in RA and Dec, μ is the resulting vector modulus, and β_app is the apparent speed in units of speed of light, obtained using the luminosity distance D_L and the redshift z. We did not fit accelerations since second-order terms are difficult to constrain without a long enough monitoring. Lister et al. (2009) require a component to be detected in at least 10 epochs in order to fit an acceleration, and the maximum number of TANAMI epochs for the sources in our sample is 9. For the uncertainty on the component position, we adopted the approach described in Lico et al. (2012). The error is calculated as the component size divided by its signal-to-noise ratio (S/N). The latter is calculated as the ratio between the peak brightness of the component and the map noise. For bright components, this estimate can yield unrealistically small uncertainties. As a lower limit for this case, we take the largest value between half the component size and 1/5 of the beam major axis. The uncertainty on the flux is given by the calibration errors, which we estimate to be 15%.

Since the TANAMI array composition can vary significantly across epochs, the (u, v) coverage, and consequently the beam will also be inhomogeneous across the epochs for a single source. This complicates the identification of components across the epochs, since one component detected in a low-resolution epoch may be resolved into multiple components in another epoch. To overcome this problem, we re-imaged the sources applying a Gaussian taper to the highest-resolution maps in order to downweight the visibilities from the longest baselines, and approximately match the beam size of the epoch with lowest resolution along the jet direction.

Using the apparent speed measured from kinematics (β_app) together with the jet-to-counterjet flux density ratio R = S_jet/S_counterjet, it is possible to set some limits on the intrinsic jet speed β = v/c and viewing angle θ.

R and β_app can be expressed as a function of these two parameters:

(3)

where α is the jet spectral index, and the index p takes a value of 2 or 3 depending on whether R is calculated integrating the flux density over the jet or using a single component, respectively. We used the former method for cases where a counterjet is detected, while we used a single component in the sources with no detectable counterjet. In the cases where no counterjet emission is detected, we placed a lower limit on R. To do this, we took the ratio between the peak flux of the innermost jet component and the peak flux on the counterjet side. Out of the different measurements coming from the multiple available epochs, we used the most constraining one, given by the map with the best quality. These relations allowed us to set limits in the space of intrinsic parameters β and θ (see Sect. 4.2).

3.2. Fermi-LAT data

The LAT is a pair-conversion telescope, launched on 2008 June 11 as the main scientific instrument on board the Fermi Gamma-ray Space Telescope (Atwood et al. 2009). Its main energy range is 0.1–100 GeV, but its sensitivity extends down to 60 MeV and up to 2 TeV (Ajello et al. 2017). The telescope operates almost exclusively in survey mode, scanning the entire γ-ray sky approximately every three hours.

We used the Python package Fermipy (Wood et al. 2017) throughout the analysis, and assumed as starting model the third Fermi-LAT source catalog (3FGL, Acero et al. 2015). We considered a Region of Interest (ROI) of 10° around the target position, and included in the model all sources from the 3FGL within 15° from the ROI center, together with the models for the galactic and isotropic diffuse emission (gll_iem_v06.fits and iso_P8R2_SOURCE_V6_v06.txt, respectively). We performed a binned analysis with 10 bins per decade in energy and 0.1° binning in space. We used the source event class, and the Pass 8 response function P8R2_SOURCE_V6. We enabled an energy dispersion correction to take into account the degradation in energy resolution at low energies. We selected only times when the zenith angle of the incoming events with respect to the telescope was smaller than 100°, to avoid contamination by the Earth’s limb. We took advantage of the new Pass 8 characterization of events in different PSF quartiles, based on the quality of the direction reconstruction, from the worse quartile (PSF0) to the best (PSF3). We modeled each type separately and combined the resulting likelihood using the summed-likelihood method. Since the angular resolution of the LAT degrades at low energy, we progressively increased the low-energy cut for the worse PSF quartiles. The details of the selection of the different components are listed in Table 3.

We fit the ROI with the initial 3FGL model, freeing all the parameters of the target source and the normalization of all sources within 5° of the ROI center. Since our data set more than doubles the integration time with respect to the 3FGL catalog, we looked for new sources with an iterative procedure. We produced a map of excess Test Statistic (TS). The TS is defined as 2log(L/L₀) where L is the likelihood of the model with a point source at the target position, and L₀ is the likelihood without the source. A value of TS = 25 corresponds to a significance of 4.2σ (Mattox et al. 1996). A TS map is produced by inserting a test source at each map pixel and evaluating its significance over the current model. We looked for TS> 9 peaks in the TS map, with a minimum separation of 0.3°, and added a new point source to the model for each peak, assuming a power-law spectrum. We then fit again the ROI, and produced a new TS map. This process was iterated until all significant excesses were modeled out. We also performed a localization analysis on the target source and all new sources with TS > 25 found in the ROI.

As listed in Table 1, roughly half of the TANAMI radio galaxies are detected in the 3FGL catalog. For those that were not already extensively studied in previous TANAMI papers (see Sect. 2) we considered 103 months of LAT data (∼8.5 years). We produced 0.1–100 GeV Spectral Energy Distributions (SEDs) with a binning optimized for each source in order to balance spectral coverage with sufficient statistics in each bin. We also computed light curves with monthly bins and, when possible, a finer binning to reveal more detailed flaring activity. We also report the average source properties in the 1–100 GeV range in order to compare with those presented in the 3FGL.

For the variability analysis, we looked at the energy range 0.1–300 GeV. All the other parameters of the multi-component analysis are the same as those described above. We first performed a standard analysis over the whole time range, then we performed a full likelihood fit in each time bin. The number of free parameters in the latter fit depends on the statistics in each bin. We first attempted a fit with the normalization of all sources in the ROI left free to vary. If the statistics did not allow the fit to converge, we progressively reduced the radius including sources with free normalization, first to 3°, then to 1°. If all these methods failed to result in a successful fit, we freed only the normalization of the target source.

We consider the source detected in each bin if TS > 10 and the signal to noise ratio is higher than 2, otherwise we place a 95% confidence flux upper limit.

4. Results

4.1. Radio imaging results

In Fig. 1 we present first-epoch VLBI images of the γ-ray loud TANAMI radio galaxies, while the full set of multi-epoch images and related map parameters is presented in Appendix A. The imaging results for each individual source are discussed in this subsection.

Fig. 1. First-epoch 8.4 GHz contour maps of γ-ray detected TANAMI radio galaxies. The blue ellipse indicates the convolving beam, while the black bar shows the linear scale. Top left to bottom right: Pictor A, PKS 0521−36, PKS 0625−35, Centaurus B, PKS 1718−649 (the B1950 name is indicated in the top right corner of each image). The full set of multi-epoch images and the associated map parameter tables are presented in Appendix A.

0518–458 (Pictor A). Pictor A is a classical powerful FR II radio galaxy. In a previous kinematic study Tingay et al. (2000) characterized the pc-scale jet of the source with three components, with a fastest apparent motion of β_app = 1.1 ± 0.5. They did not detect any counterjet at this scale, and additionally they found an apparent bend in the jet at ∼10 mas from the core.

TANAMI monitoring has provided 5 epochs for this source. As shown in Fig. 1, the jet extends for ∼30 mas westward from the brightest component, and several features can be identified and tracked. We consistently detect emission eastward of the brightest component, which we assumed to be the VLBI core, in all epochs. This feature is not present in the first epoch map of Ojha et al. (2010), but its detection in multiple epochs in this work indicate that it is real. Additionally, this emission feature is also detected at 22.3 GHz.

To test if this emission feature can be considered as counterjet emission, we produce a spectral index map of the source between the quasi-simultaneous 8.4 GHz (epoch 2008 Nov. 27) and 22.3 GHz (epoch 2008 Nov. 29) images.

The resulting map is presented in Fig. 2. We see that the spectrum flattens around the brightest component, and becomes optically thin again in the eastward region. This indicates that this is counterjet emission, and should therefore be taken into account when computing jet-to-counterjet flux density ratios.

Fig. 2. Spectral index map of Pictor A (0518−458) between 8.4 GHz and 22.3 GHz for epoch 2008 Nov. 27. Black contours are from the 8.4 GHz image. The convolving beam is represented in gray in the lower-left corner.

We investigated the significance of a putative jet bending by plotting the position angle of the Gaussian components versus the radial distance. This is shown in Fig. 3. The error on the position angle ψ has been calculated adapting Eq. 9 from Schinzel et al. (2012), as:

(4)

Fig. 3. Position angle of Gaussian jet components for Pictor A (0518−458) versus radial distance from the core. The color coding indicates different components, and is the same as the upper left panel of Fig. 7. The sign of the counter-jet component CJ1 has been reversed to improve visualization.

where a is the component size, r is the component radial distance and the S/N is defined as in Sect. 3.1. There is no obvious break in the distribution around ∼10 mas, in contrast with the findings of Tingay et al. (2000).

0521–365. This is a nearby AGN with an uncertain classification. Leon et al. (2016) derive limits on the jet viewing angle, speed, and Doppler factor using the Atacama Large Millimeter Array (ALMA). Their results suggest a jet viewing angle in the range 16 ° ≤θ ≤ 38°, from the jet-to-counterjet ratio. Their detection of a large scale double structure already suggests that Doppler boosting effects in this source are not dominant on the kpc scale. Giroletti et al. (2004) studied PKS 0521−36 as part of a sample study on low-redshift BL Lacs. Using 5 GHz VLBA observations, the authors constrain the viewing angle of the source to lie in the range 21 ° ≤θ ≤ 27°. D’Ammando et al. (2015) constrain the same parameters using SED modeling including γ-ray data, obtaining a more aligned jet viewing angle of 6 ° ≤θ ≤ 15°. These results point to an intermediate jet viewing angle between a blazar and a steep spectrum radio quasar (SSRQ) or Broad Line Radio Galaxy (BLRG). However, due to its uncertain nature, PKS 0521−36 was not included in the large study of γ-ray properties of misaligned AGN by Abdo et al. (2010c).

Previous VLBI observations performed with the VLBA and with the Southern Hemisphere VLBI Experiment (SHEVE) at 4.9 GHz and 8.4 GHz provided an upper limit on the apparent speed of jet components β_app <  1.2 (Tingay & Edwards 2002). This is also consistent with the hypothesis that the jet of PKS 0521−36 is not strongly beamed.

Due to its known variability at γ-ray energies (D’Ammando et al. 2015), PKS 0521−36 is one of the more densely monitored sources in the sample, with nine TANAMI epochs, which provide an excellent data set for a kinematic analysis. The full resolution maps and the corresponding image parameters are presented in Appendix A.

Our TANAMI images show a faint jet flow extending out to ∼60 mas from the core to the north-west (see Fig. 1). There is a drop in the jet brightness around 10-15 mas from the core, which is seen consistently in all epochs. The jet structure is remarkably consistent across the epochs, hinting at a slow speed. The extended jet is not detected in the last epoch due to the lack of the shortest baseline ATCA-Mopra which provides the necessary sensitivity to the larger scale emission.

We produced a spectral index map between the quasi-simultaneous 8.4 GHz (epoch 2008 Mar. 28) and 22.3 GHz (epoch 2008 Mar. 26) images, which is presented in Fig. 4.

Fig. 4. Spectral index map of PKS 0521−36 (0521−365) between 8.4 GHz and 22.3 GHz for epoch 2008 Mar. 28. Black contours are from the 8.4 GHz image. The convolving beam is represented in gray in the lower-left corner.

0625–354. This is an FR I radio galaxy, but shows an optical spectrum similar to a BL Lac object. Its γ-ray properties (see Sect. 4.3) also suggest a moderately aligned jet, similar (but less extreme) to the case of IC 310 (Aleksić et al. 2014). PKS 0625−35 is also one of the very few radio galaxies to have been detected at TeV energies by Cherenkov telescopes (Abdalla et al. 2018).

TANAMI is the first multi-epoch VLBI data set for this source, providing 9 epochs. A previous single-epoch observation with the LBA at 2.29 GHz by Venturi et al. (2000) provided constraints on the jet angle to the line of sight and intrinsic jet speed using a lower limit on the jet-to-counterjet ratio and an estimate of the core dominance. The latter method gives the most constraining estimates of θ ≤ 43° and β ≥ 0.74.

Our TANAMI images show a one sided core-jet structure with a faint extended jet in the south-east direction (see Fig. 1). We produce a spectral index map between 8.4 GHz (epoch 2008 Nov. 27) and 22.3 GHz (epoch 2008 Nov. 29), which is presented in Fig. 5. The index is relatively steep even in the brightest region.

Fig. 5. Spectral index map of PKS 0625−35 (0625−354) between 8.4 GHz and 22.3 GHz for epoch 2008 Nov. 27. The black contours are from the 8.4 GHz image. The convolving beam is represented in gray in the lower-left corner.

1343–601 (Centaurus B). This classic FR I radio galaxy was added to the TANAMI sample after being detected by Fermi-LAT in the second source catalog (Nolan et al. 2012). TANAMI imaging reveals a smooth, one-sided jet extending in the south-west direction (see Fig. 1). Only two epochs are currently available, which is insufficient to perform a robust kinematic analysis. Using one well-defined jet component, we can obtain a rough upper limit on the jet apparent speed, which is lower than 0.89c.

1718–649. This is one of the most classic examples of the Compact Symmetric Object (CSO) and Gigahertz-Peaked Spectrum (GPS) source classes (see for example, O’Dea 1998; Giroletti & Polatidis 2009), that is, a young radio galaxy. These sources are typically small (linear size < 1 kpc), and the radio emission is dominated by symmetric mini-lobes which can show hot-spots similar to the large scale lobes of FR II radio galaxies. The advance speed of these hot-spots provides a kinematical age estimate for these sources, which is typically t_age <  10⁵ yr (Tingay et al. 1997).

TANAMI provides the first multi-epoch data set for this well-studied source, and therefore the first opportunity for a measurement of its kinematical age. The TANAMI images shows a clear structure of two components separated by ∼8 mas (see Fig. 1), which is consistent across the epochs.

We present a spectral index map between 8.4 GHz and 22.3 GHz in Fig. 6. The spectral morphology suggests that the core of the young radio source is strongly absorbed at these frequencies, and is located between the two emission components, corresponding to the region of highly inverted spectral index.

Fig. 6. Spectral index map of PKS 1718−649 between 8.4 GHz and 22.3 GHz for epoch 2008 Feb. 07. The black contours are from the 8.4 GHz image. The convolving beam is represented in gray in the lower-left corner.

4.2. Radio kinematic analysis results

The main results of the kinematic analysis for Pictor A, PKS 0521−36, PKS 0625−35 and PKS 1718−649, following the procedure described in Sect. 3.1, are summarized by the plots of the identified Gaussian component’s core separation versus time (Fig. 7 and 8). The resulting values for the component speed and estimated ejection date are listed in Tables 4 through 7. Only components with at least five epochs are listed, and ejection dates are given only for sources with speed not consistent with zero. Finally, the parameter space of intrinsic jet parameters β_app and θ allowed by our results is illustrated in Fig. 10, for sources with significant component motion. Additional information illustrating the kinematic analysis results is provided in Appendix B.

Fig. 7. Jet kinematics of our radio galaxies: core distance of jet features as a function of time. The solid lines represent a least squares fit to their positions (the slope is the apparent speed). Top left to bottom right: Pictor A, PKS 0521−36, PKS 0625−35, PKS 1718−649. The B1950 name is given above each image.

Fig. 8. Jet kinematics of PKS 0521−36 (0521−365): core distance of jet features as a function of time, including the previous VLBI dataset of Tingay & Edwards (2002). The solid lines represent a least squares fit to their positions (the slope is the apparent speed). Compare with top-right panel of Fig. 7.

0518–458 (Pictor A). We re-imaged the source applying a Gaussian taper in order to match the resolution of epoch 2010 Jul. 24 (see Fig. A.1).

We find at least mildly relativistic apparent speeds with a minimum significant value of 0.1c (J1) and a maximum of 2.6c (J2), as allowed within the 1σ errors. The lower limits are consistent with the estimates in Tingay et al. (2000), while the increased number of epochs allows us to reveal one component which is not consistent with subluminal motion (currently at a significance level of ∼1.3σ) in the jet of Pictor A.

The upper left panel of Fig. 10 shows the resulting limits on the intrinsic jet speed and viewing angle for Pictor A obtained via the combination of the apparent speed information with the jet-to-counterjet ratio. The estimates for R represent the minimum (dashed blue line), mean (continuous blue line) and maximum (dot-dashed blue line) values. The central estimate for β_app is the one for the fastest component, while the minimum and maximum values represent its error. The relatively small value of the jet ratio combined with the mildly superluminal apparent speed results in a tight constrain on the viewing angle, which lies in the range 76 ° < θ <  80°. To account for the observed apparent speed with such a large jet angle, the intrinsic jet speed should be β >  0.96.

0521–365. We re-imaged all epochs applying a Gaussian taper, in order to approximately match the resolution of epoch 2010 Mar. 12 (see Fig. A.2).

Since we obtain slow jet speeds, we attempted to fit our TANAMI jet model together with the previous 8.4 GHz VLBI dataset from Tingay & Edwards (2002). It is possible to cross-identify and fit four components between the two datasets. This is shown in Fig. 8. The updated values of μ and β_app for these components are reported in Table 6. We note that the measured apparent speed values are significantly larger when using a longer time range for the kinematic analysis. This indicates that the limited time coverage of our TANAMI observations may lead us to underestimate the apparent speed in the case of slow apparent motions, which typically require a larger time range to be adequately constrained (see for instance, Piner & Edwards 2018).

The upper right panel of Fig. 10 shows the resulting limits on the intrinsic jet speed and viewing angle for PKS 0521−36. Since there is no counterjet, we place a lower limit on R as described in Sect. 3.1. For β_app the values adopted are the minimum and maximum observed values, considering the four components that are cross-identified with the Tingay & Edwards (2002) dataset. We obtain a range of β >  0.67 and θ <  26°, respectively.

0625–354. The outer jet of this source is too faint to be modeled reliably with Gaussian components, therefore we re-imaged only the inner ∼20 mas for the kinematic analysis.

Müller et al. (2012) presented a preliminary kinematic analysis of the TANAMI data on PKS 0625−35, using the first six epochs. They found a highest apparent speed of β_app = 3.0 ± 0.5, which is consistent with our result (component J3), including three more epochs. Our results provide a robust confirmation of superluminal component motion in the pc-scale jet of PKS 0625−35.

The lower panel of Fig. 10 shows the limits on the intrinsic jet speed and viewing angle for PKS 0625−35 resulting from our observations. In this case, again, we do not detect a counterjet, therefore we place a lower limit on R. The estimates for β_app are given by the fastest observed speed and its uncertainty. Our observations limit the intrinsic jet parameters to β >  0.89 and θ <  53°.

1343–601 (Cen B). As mentioned above, there are only two available TANAMI epochs for Cen B, separated by about nine months. Although it is not possible to perform a detailed kinematic analysis, as for the other sources, the data can give some indications regarding the presence or absence of motions in the source. We tentatively identify a local maximum in the brightness distribution of the Cen B VLBI jet at a distance of ∼25 mas downstream of the core in both our images, which were taken approximately nine months apart (see Fig. 9). If this association is correct, which will be tested by forthcoming TANAMI epochs, the apparent jet speed is likely subluminal or at most mildly superluminal.

Fig. 9. Multi-epoch images of Centaurus B (1343−601). The colored crossed circles represent the circular Gaussian components that have been fit to the clean maps.

Fig. 10. Parameter space of intrinsic jet speed β and viewing angle θ allowed by our observations. The blue shaded area is the one allowed by the measurement of R (θ as function of β given R, that is, θR), while the red shaded area is the one allowed by the observed βapp (θ as function of β given βapp, that is, θβapp). For each source we provide a minimum, maximum and (except for PKS 0521−36) a central estimate of R and βapp. The top-right legend reports the resulting limits on θ and β. The dashed colored lines and boxes indicate constraints from previous works, namely Hardcastle et al. (2016), Tingay et al. (2000) for Pictor A (0518−458), Pian et al. (1996), Giroletti et al. (2004), D’Ammando et al. (2015) for PKS 0521−36 (0521−365), and Venturi et al. (2000) for PKS 0625−35 (0625−354).

1718–649. In this case the source is more complex than the classic core-jet morphology seen in the other radio galaxies presented here. We reference the kinematic analysis at the position of the brightest component (C1), even though it is not the core of the radio source in this case. The resulting kinematic values represent the evolution of the distance between the two components consistently detected in all epochs. We have preferred such a simple two-component modeling of the brightness distribution, although finer structures are clearly seen in certain epochs. However, due to differences in image noise and dynamic range across the epochs, an attempt to identify and track such faint substructures in a kinematic analysis would yield unreliable results. We therefore chose a more simplistic modeling approach.

In the lower-right panel of Fig. 7 we plot the distance between the two components (referenced to the brightest one) as a function of time, and the corresponding linear fit. The angular separation speed is μ = (0.13 ± 0.06) mas yr⁻¹, and the corresponding apparent linear speed is β_app = 0.13 ± 0.06. We are therefore able to estimate when the young radio source first ejected its two symmetric component, and find a zero-separation epoch of 1963 ± 22.

4.3. Fermi-LAT results

0518–458. Pictor A was reported as a Fermi-LAT detection for the first time by Brown & Adams (2012) based on three years of data, and confirmed in the 3FGL. It is the faintest source among the 3FGL TANAMI radio galaxies. While it has been established that diffuse structures such as lobes and hot spots may give a significant or even dominant contribution to γ-ray emission in radio galaxies (Abdo et al. 2010b; Ackermann et al. 2016), Brown & Adams (2012) found that this is not the case for Pictor A. Through SED modeling of the western hot-spot, the authors found that Synchrotron Self-Compton (SSC) and Inverse Compton with CMB photons (IC/CMB) models failed to reproduce the X-ray and γ-ray data at the same time, and therefore concluded that the γ-ray emission must be dominated by the AGN jet. This is also supported by the indications of variability observed in the LAT data.

Our analysis over ∼8.5 years results in a higher significance w.r.t. the 3FGL, and the source properties are consistent with the catalog results within the errors (see Table 8).

We produced a monthly light curve over the whole ∼8.5 years time period, as described in Sect. 3.2 (see top panel of Fig. 11). Pictor A is a relatively faint γ-ray source, therefore it is not detected on monthly timescales for most of this time range. The high state in September 2010 provided the statistics to allow a detection integrating over the first three years of Fermi-LAT data, with a monthly significance of TS ∼ 40. In order to better characterize the variability of Pictor A, we also produced 3-month and 6-month-binned light curves, shown in the center and bottom panel of Fig. 11, respectively. Integrating over longer timescales allows us to detect the source for a larger fraction of the time range (27% and 47%, respectively).

Fig. 11. Light curve of Pictor A (0518−458) between 0.1 and 300 GeV over 103 months of Fermi-LAT data. Blue points are detections, red arrows are upper limits. Top to bottom: 1-month, 3-months and 6-months binning, respectively.

0521–365. This source is a very bright γ-ray emitter, showing significant variability. Monthly and weekly light curves over ∼8.5 years of data are shown in Fig. 12. The flaring activity in 2010-2011 has been studied by D’Ammando et al. (2015) down to 12-h time scales. We find a second period of activity of comparable magnitude at the end of 2012. We investigated the short variability in this period by producing daily and 6-h time scale light curves, which are presented in Fig. 13. The source is especially variable on these time scales in 2012 Dec. The maximum peak on 6-h time scales is observed on 2012 Dec 12, with a flux of (4.0 ± 1.0)×10⁻⁶ photons cm⁻² s⁻¹. The sharpest flux change is the decrease from this bin to the following 6-h bin, with a variation of a factor ∼3.

Fig. 12. Light curve of PKS 0521−36 (0521−365) between 0.1 and 300 GeV over 103 months of Fermi-LAT data, with weekly (blue points) and monthly (red points) binning. Upper limits are indicated by arrows of the respective colors. The left y-axis reports the weekly flux values, the right one reports the monthly flux values.

Fig. 13. Light curve of PKS 0521−36 (0521−365) between 0.1 and 300 GeV over 103 months of Fermi-LAT data, with 6-h (blue points) and daily (red points) binning. Upper limits are indicated by arrows of the respective colors. The left y-axis reports the 6-h flux values, the right one reports the daily flux values.

0625–354. This source has been detected by Fermi-LAT since the first source catalog (Abdo et al. 2010a). It is a recent addition to the small sample of six radio galaxies which have been detected by Cherenkov Telescopes in the TeV range (Abdalla et al. 2018). This is mostly due to its notably hard spectrum (see Table 9). The source properties in the 3FGL energy range are consistent with the catalog values, and the significance is increased.

A monthly-binned light curve of the source over 8 years is presented in Fig. 14. No significant activity is detected.

Fig. 14. Light curve of PKS 0625−35 (0625−354) between 0.1 and 300 GeV over 103 months of Fermi-LAT data, with monthly binning. Blue points are detections, red arrows are upper limits.

1343–601. Centaurus B was detected by Fermi-LAT in the 2FGL (Nolan et al. 2012). Our analysis over ∼8.5 years in the 3FGL range gives results consistent with the catalog values within the error. It is interesting to note that in spite of its relatively steep spectral index, Cen B has been indicated as a good candidate for a TeV detection (Angioni et al. 2017). This is probably favored by its very low redshift (see Table 1).

The analysis of this source is particularly challenging due to its location behind the Galactic plane (b = 1.73°), which contains rich diffuse emission structures, complicating point-source analysis. This ROI required the addition of many sources in excess of the 3FGL, since we are using more than double the data with respect to the catalog, some of which may actually be due to contribution from improperly modeled Galactic diffuse emission.

A monthly light curve of Cen B is presented in Fig. 15. The source is undetected for most of the time range, and does not show any significant variability.

Fig. 15. Light curve of Centaurus B between 0.1 and 300 GeV over 103 months of Fermi-LAT data, with monthly binning. Blue points are detections, red arrows are upper limits.

1718–649. This source has recently become the first γ-ray detected young radio galaxy. Migliori et al. (2016) detected the source using 7 years of LAT data, with TS ∼ 36, confirming it as the counterpart of the unidentified catalog source 3FGL J1728.0−6446. We double-checked this result using ∼8.5 years of data. In order to avoid any bias toward either of the two hypotheses, that is, an unidentified γ-ray source at the position of 3FGL J1728.0−6446, or a γ-ray source associated to PKS 1718−649, we removed the catalog source from the model and evaluated the excess statistics in the ROI. Fig. 16 shows a TS excess map after subtracting 3FGL J1728.0−6446. The excess is nicely coincident with the position of PKS 1718−649. A new source is found by the source-finding algorithm (dubbed PS J1724.2−6459, for “Point Source”), and after localization we found that the source coincides with the position of the target within the errors.

Fig. 16. Fermi-LAT map of excess TS in the inner region of the ROI centered on PKS 1718−649, after removing the unidentified catalog source 3FGL J1728.0−6446. The ellipses represent the 68% and 95% confidence positional errors on the catalog source, while the black circles represent the same errors for the new source PS J1724.2−6459. The map size is 1.6 ° ×1.6°. Each pixel corresponds to 0.1°.

5. Discussion

We now discuss the results presented in the previous sections, focusing first on the individual sources, and then on the properties of γ-ray radio galaxies as a sample.

5.1. Individual sources

0518−458. Based on their VLBI data, Tingay et al. (2000) place an upper limit on the jet viewing angle of θ <  51°. This is consistent with the allowed range of parameters from our TANAMI observations, that is 40 ° < θ <  78° and β >  0.78.

Hardcastle et al. (2016) observed the jet of Pictor A in X-rays using Chandra observations, combined with ATCA observations of the large scale radio structure. They detect an X-ray counterjet, and place an upper limit on the viewing angle of θ <  70° and a lower limit on the intrinsic jet speed of β >  0.3 from the jet sidedness. The authors noted that at the time there was no direct evidence of bulk relativistic motion in the jet of Pictor A, and pointed out the value of future VLBI studies on the source. Our results provide the first robust measurement of component motion in the pc-scale jet of Pictor A. Our lower limit on the intrinsic jet speed does not constrain the jet to be highly relativistic, but significantly reduces the parameter space with respect to the estimate by Hardcastle et al. (2016). Our lower limit on the viewing angle confirms that the X-ray emission mechanism cannot be Inverse Compton on CMB photons (IC/CMB), since this would require θ smaller than a few degrees and β ∼ 1.

The γ-ray light curve of Pictor A shows that the source is often undetected on monthly time scales. In a variability study of 12 radio galaxies using four years of Fermi-LAT data Grandi et al. (2013) found that FR II sources were detected only during short periods of activity (see for example, 3C 111, Grandi et al. 2012), while FR I sources were detected for a larger fraction of the investigated time range. This seems to be consistent with our results on Pictor A, a classic FR II radio galaxy, which is detected for a fraction of ∼5% of the time over 103 months.

Our TANAMI data, combined with the γ-ray variability allows us to shed some light on the nature of the high-energy emission in Pictor A. We note that a new VLBI jet feature emerged from the radio core of Pictor A between 2010 Jul. 24 and 2011 Aug. 13. Since the elevated γ-ray state that allowed the first Fermi-LAT detection was within this time range, it is plausible that it was caused by the passage of a new shock through the radio core, corresponding to the appearance of this new radio feature. This hypothesis is in agreement with the results of Brown & Adams (2012), who modeled the SED of the western hot-spot of Pictor A, and found that it cannot reproduce the observed X-ray and γ-ray emission at the same time. Therefore the LAT detection is probably associated with emission from the innermost part of the jet.

One can speculate whether the association with pc-scale jet activity is a defining feature of γ-ray emitting FR II radio galaxies. If our inference is correct, Pictor A would be yet another example of this kind of behavior, after 3C 111 (Grandi et al. 2012) and 3C 120⁴ (Casadio et al. 2015).

This would nicely fit with the findings of Grandi et al. (2013), if we assume that FR II sources can only be detected when there is significant activity in their innermost jet, and/or when the inner jet is temporarily well-aligned with our line-of-sight (see Casadio et al. 2015; Janiak et al. 2016, for the case of 3C 120). The presence of jet spine “wiggling” and precession is suggested by numerical simulations such as for instance, Liska et al. (2018). This would naturally lead to a low duty cycle for γ-ray activity, and consequently to fewer detected FR IIs, exactly as observed (see for example Grandi et al. 2016, and references therein).

0521–365. Our VLBI results complete the multi-wavelength picture provided by Tingay & Edwards (2002), D’Ammando et al. (2015), and Leon et al. (2016), confirming that the jet of PKS 0521−36 is not highly beamed, with viewing angles larger than 10° still compatible with our observations.

Pian et al. (1996) used multi-wavelength data to model the broadband SED of PKS 0521−36, and found that the source is likely not Doppler-boosting dominated, with an estimated viewing angle of θ = 30 ° ±6° and a Doppler factor δ ∼ 1.5. This range of parameters is completely ruled out by our kinematic analysis, as can be seen in Fig. 10⁵.

Fast γ-ray variability is typical of blazar jets, where the time scales are strongly reduced by the large Doppler factors. However, as shown by the VLBI results presented here, and by previous multi-wavelength studies, it is unlikely that the jet of PKS 0521−36 is strongly affected by Doppler boosting. Therefore the jet morphology does not appear to change in response to the strong γ-ray flaring activity. On the other hand, we observe a doubling of the VLBI core flux density during the first γ-ray flaring periods of 2010-2011 (see Fig. B.5), which suggests that the γ-ray emission region is located inside the radio core, and not in the jet.

This combination of slow pc-scale jet and fast γ-ray variability bears some resemblance to the case of IC 310, a transitional FR I-BL Lac object which shows minute-timescale variability at VHE (Aleksić et al. 2014), but no fast jet motion in VLBI images (Schulz 2016). For this source, Aleksić et al. (2014) favor a model where the fast γ-ray variability is produced by charge depletion in the supermassive black hole magnetosphere due to low accretion rate phases.

0625–354. The estimated range of viewing angles (θ <  38°) from our VLBI data does not allow us to settle the uncertain classification of PKS 0625−35, as it is consistent with jet orientations typical of radio galaxies (see for instance, Kim et al. 2018, for the case of M 87), but also allows for smaller viewing angles, more typical of blazar sources. The transitional nature of this source is also supported by the indications of a hard X-ray nuclear component and by its location in the parameter space of radio core luminosity at 5 GHz and X-ray non-thermal luminosity, which places it exactly in the region between FR Is and BL Lacs (Trussoni et al. 1999). Finally, a small viewing angle would partially explain its TeV detection (Abdalla et al. 2018). In this respect, our finding of superluminal jet motions in PKS 0625−35 appears in contrast with what is observed for the bulk of the TeV blazars population, which typically shows slow jet speeds. However, a recent study by Piner & Edwards (2018) has shown for the first time a superluminal tail in the jet speed distribution of TeV sources. PKS 0625−35 thus appears to belong to this elusive class of TeV-bright AGN jets with fast parsec-scale motion.

1718–649. Our kinematic analysis allows us to obtain a rough estimate on the age of this young radio source, albeit with a large relative error, as t_age = d_max/μ = (70 ± 30) years. This limit is consistent with a quite young age for the source, compared to usual estimates for young sources (t_age ∼ 10⁵ yr, Tingay et al. 1997, and references therein). This is mostly due to the notably small linear size of ∼2.5 pc, compared with the bulk of the young radio sources population (see for example Orienti & Dallacasa 2014). Our estimates are also in broad agreement with the one reported in Giroletti & Polatidis (2009). According to young radio source evolutionary models (see for instance, Snellen et al. 2000; An & Baan 2012; Orienti & Dallacasa 2014), at the scales relevant to PKS 1718−649 the flux density should decrease going back in time in the life of the source. If our estimate is correct, this could potentially be compared to the earliest flux density measurements for this object.

5.2. Sample properties

In this subsection, we investigate the correlated radio and γ-ray properties of the LAT-detected radio galaxies in the TANAMI program. To increase the sample size, we added all radio galaxies with published VLBI results from the MOJAVE survey, and performed the same LAT analysis described in Sect. 3.2. The resulting γ-ray properties of the sources with a significant γ-ray detection are listed in Table 10. We detect (TS > 25) three sources not previously published in any Fermi-LAT catalog, that is, NGC 315, NRAO 128, and PKS 1514+00, in addition to the well-known γ-ray source 3C 120. PKS 1128−047 is not classified as radio galaxy in the Fermi-LAT catalogs, but it is according to the NASA/IPAC Extragalactic Database (NED) and MOJAVE. We exclude NGC 1052 due to its lack of a clear VLBI core in the MOJAVE data, although we are aware that at higher frequencies the absorption by the surrounding torus is not present and a clear core is detected at 86 GHz (Baczko et al. 2016).

The resulting subsample includes sixteen radio galaxies with measured VLBI kinematics and LAT properties. This is the largest sample of γ-ray detected radio galaxies studied with VLBI techniques so far. For the MOJAVE sources, the kinematics results were taken from Lister et al. (2013, 2016). The exact reference for each source is given in Table 10.

In order to combine the VLBI data from the TANAMI and MOJAVE samples, we had to account for the different frequency in the two surveys, that is 8.4 GHz and 15 GHz, respectively. In order to do so, we have back-extrapolated the MOJAVE core fluxes down to 8.4 GHz using the spectral index values provided in Hovatta et al. (2014).

In order to investigate possible correlations between the radio and γ-ray properties of the sources in this sample, we have used the Kendall’s correlation coefficient (τ). The correlation coefficient is equal to zero in the case of uncorrelated data, one in case of maximum correlation, and minus one in case of maximum anti-correlation. The resulting correlation coefficients with errors and the relative p-values are listed in Table 11.

In Fig. 17 we show the maximum measured apparent speed versus the total VLBI luminosity of this subsample. It can be seen that there is a lack of low-luminosity sources (≲10²⁴ W Hz⁻¹) with high apparent speed (β_app  ≳  0.5). This was already found to be true for blazars (see for example Cohen et al. 2007; Lister et al. 2013), but this is the first time that it is observed in radio galaxies. On the other hand, while for blazars a population of high-luminosity, low-speed sources is observed (Cohen et al. 2007), this is not the case in our radio galaxy sample. Cohen et al. (2007) interpreted this subpopulation, mostly composed of BL Lac blazars, as sources where the measured VLBI apparent speed at 15 GHz is not representative of the jet bulk Lorentz factor. This is now a well-established fact for high-energy peaked BL Lacs (HBLs, see Piner & Edwards 2018, and references therein). We do not find clear evidence for a misaligned parent population of these intrinsically luminous, low-speed blazars.

Fig. 17. Maximum measured apparent speed βapp as a function of median total VLBI luminosity for our TANAMI γ-ray radio galaxies (blue circles) and MOJAVE γ-ray radio galaxies (red squares). The downward arrows indicate upper limits on the apparent speed, for sources with speed consistent with zero within the uncertainties.

In Fig. 18 we show the LAT flux above 100 MeV versus the median VLBI core flux density. We test for a linear correlation between the log core flux density and the log of the γ-ray flux, and find that a Kendall’s test yields a correlation coefficient τ = 0.55 and a p-value of 0.3%. This result suggests that in radio galaxies, the observed γ-ray flux is related to the radio flux density of the innermost pc-scale jet. Such a correlation has already been highlighted for large, blazar-dominated AGN samples (see for example, Kovalev et al. 2009; Ackermann et al. 2011; Lico et al. 2017). In a study of the radio and γ-ray properties of the full TANAMI sample using the first year of Fermi-LAT data, Böck et al. (2016) found a similar correlation between radio and Fermi-LAT fluxes. This sample was dominated by blazars, while our (smaller) sample is focused on radio galaxies alone. We can therefore confirm that a similar dependence of γ-ray flux on the radio core flux density holds for misaligned AGN jets as well. The idea that γ-ray emission in radio galaxies with compact VLBI emission is related to the innermost jet is supported by the recent detection of diffuse γ-ray emission from the large-scale lobes of Fornax A (Ackermann et al. 2016). In this radio galaxy, the γ-ray emission has been determined to be spatially extended at > 5σ confidence level. The contribution of the extended lobes to the high-energy flux is dominant, while the radio core contribution is no more than 14%. Interestingly, Fornax A has not been detected with VLBI observations, suggesting that it does not show a bright, compact radio jet on parsec scales. Therefore, the subdominant contribution of the radio core to the γ-ray emission in Fornax A is to be expected, if the radio brightness of the innermost jet is indeed related to its γ-ray flux, as our results indicate.

Fig. 18. Fermi-LAT flux above 100 MeV as a function of median radio VLBI core flux density for our TANAMI γ-ray radio galaxies (blue circles) and MOJAVE γ-ray radio galaxies (red squares).

In Fig. 19 we show the LAT flux versus the VLBI jet flux, that is, the flux obtained subtracting the core modeled flux density from the total imaged flux density. A Kendall’s correlation test for these two quantities yields a coefficient τ = 0.27 with a p-value of 2%. Since this exceeds the 1% threshold, we conclude that the γ-ray flux is not significantly related to the pc-scale jet flux density.

Fig. 19. Fermi-LAT flux above 100 MeV as a function of median radio VLBI jet flux density for our TANAMI γ-ray radio galaxies (blue circles) and MOJAVE γ-ray radio galaxies (red squares).

In Fig. 20 we show the LAT luminosity versus the median VLBI core luminosity. The clear linear correlation observed between the log quantities is consistent with a 1:1 correlation, which suggests that it is induced by the common dependence on redshift. The linear correlation results in τ = 0.78 and a p-value = 2 × 10⁻⁵. When accounting for the redshift as a third parameter, the partial correlation analysis (see Sect. 2.1 of Akritas & Siebert 1996) yields τ = 0.39 and a p-value = 0.05. One might assume that controlling for the effect of redshift when correlating luminosity with luminosity should yield the same result as the flux-flux density correlation. However, there are two effects that can explain why this is not the case. First, the flux correlation and the luminosity correlation would only yield exactly the same result if all the sources were the same, without any difference with redshift. Second, while the VLBI core luminosity has a simple dependence in the form , the γ-ray luminosity is integrated over several orders of magnitude in energy, and therefore depends not only on flux and luminosity distance, but also on photon index.

Fig. 20. Fermi-LAT luminosity as a function of median radio VLBI core luminosity for our TANAMI γ-ray radio galaxies (blue circles) and MOJAVE γ-ray radio galaxies (red squares).

In Fig. 21 we show the median VLBI core brightness temperature versus γ-ray luminosity. The brightness temperature is calculated as

(5)

Fig. 21. Median radio VLBI core brightness temperature as a function of Fermi-LAT luminosity above 100 MeV for our TANAMI γ-ray radio galaxies (blue circles) and MOJAVE γ-ray radio galaxies (red squares). A correction factor has been applied to the MOJAVE data to account for the difference in frequency, but this might still imply a bias.

k_B is Boltzmann’s constant, S is the radio flux density, λ is the observing wavelength, z is the source redshift, and θ is the FWHM of the circular Gaussian model fit to the component (in our case, the core). Since the component sizes resulting from modeling can be much smaller than the actual resolution, we have calculated an upper limit on the component size based on its S/N, following Kovalev et al. (2005):

(6)

where b is the beam size, and the S/N was calculated as the ratio between the peak flux density of the residual map after subtracting the given component, in an area defined by the component size, and the post-fit rms in the same residual map area. When the fit component size was θ <  θ_lim, we have used θ_lim as an upper limit, and therefore the obtained T_b is a lower limit estimate. We did not compute the resolution limit for the MOJAVE sources, since we only made use of publicly available data.

A correction factor given by the ratio of the two wavelengths has been applied to the MOJAVE data to account for the different observing wavelength. Böck et al. (2016) found, for the full TANAMI sample, an indication of increasing core brightness temperature with increasing γ-ray luminosity. Interestingly, our results focused on γ-ray radio galaxies show that the γ-ray luminosity appears to be completely uncorrelated with the core brightness temperature for misaligned jets. This is confirmed by the correlation analysis (see Table 11). Since the core brightness temperature is often considered as an indicator of Doppler boosting (see for instance, Kovalev et al. 2005), this indicates that while the core flux density does reflect higher γ-ray fluxes, high-energy emission in radio galaxies is not Doppler boosting dominated.

Another indication supporting this inference comes from the VLBI core dominance (CD), another indicator of Doppler boosting. We define this parameter as a ratio between VLBI core flux density and total flux density, that is, . As can be seen in Fig. 22, the observed γ-ray luminosity is not correlated with the CD (see also Table 11), supporting the finding that it is not driven by Doppler boosting. In the first systematic study of the γ-ray properties of misaligned AGN, Abdo et al. (2010c) noted how γ-ray detected radio galaxies from the revised Third Cambridge catalog of radio sources (3CRR, Laing et al. 1983) showed preferentially higher CD than their undetected counterparts. This might indicate that radio galaxies with more aligned jets are preferentially detected by the Fermi-LAT. However, the definition of CD in Abdo et al. (2010c) and the one we adopted are quite distinct. Abdo et al. (2010c) define it as CD=log[S_core/(S_tot − S_core)], where S_core is the core flux density at 5 GHz measured with the Very Large Array (VLA), that is, at arcsecond scale, and S_tot is the total flux at 178 MHz on arcminute scale. It is clear how our definition of CD on pc-scales differs significantly from the one in Abdo et al. (2010c), and that therefore our results do not necessarily contradict the ones presented in this previous work.

Fig. 22. Median radio VLBI core dominance as a function of Fermi-LAT luminosity above 100 MeV for our TANAMI γ-ray radio galaxies (blue circles) and MOJAVE γ-ray radio galaxies (red squares).

Overall, the comparison of several VLBI and Fermi-LAT properties of our sample of γ-ray radio galaxies consistently suggests that high-energy emission in misaligned jet sources is not driven by orientation-dependent Doppler boosting effects, in agreement with the unified model of jetted AGN.

6. Conclusions

In this paper, we presented the first part of a study investigating the radio and γ-ray properties of radio galaxies in the TANAMI monitoring program. We presented the pc-scale jet kinematic properties of four well-known γ-ray emitting radio galaxies. High-resolution VLBI results were indicated by studies in different wavelengths as an important piece of the puzzle in order to complete our understanding of this subpopulation of AGN (see for example Tingay & Edwards 2002; Hardcastle et al. 2016). Our results fit ideally in this context. The TANAMI program provides the best multi-epoch milliarcsecond resolution images of these radio galaxies, allowing us to robustly study variations in the pc-scale jet structure for the first time. We complemented the VLBI analysis with Fermi-LAT γ-ray light curves, and investigated whether there are physical connections between the emission observed at the two ends of the spectrum. Our main results on individual sources can be summarized as follows:

• Pictor A: we find that the first γ-ray detection of this source was coincident with the passage of a new VLBI component through the compact core, an association that appears to be a defining feature of γ-ray FR II radio galaxies. Additionally, we detect a pc-scale counterjet for the first time in this source, which allows us to obtain a much better constraint on the viewing angle, which should lie approximately between 76° and 79°.

• PKS 0521–36: our VLBI results show subluminal motions on a ∼20 years time range, while the Fermi-LAT light curve shows fast variability down to 6-h time scales. Such a combination of fast high-energy flaring activity and slow jet motions bears some resemblance to the case of the hybrid FR I/BL Lac object IC 310 (Aleksić et al. 2014).

• PKS 0625–35: our TANAMI monitoring confirms the presence of superluminal motion in this jet, up to β_app ∼ 3. Our results constrain the jet viewing angle to values θ <  53°, leaving open the possibility of a transitional jet orientation between the radio galaxy and blazar classes. A small viewing angle would be consistent with the γ-ray properties of the source, that is, a hard Fermi-LAT spectrum and the observed TeV emission.

• Centaurus B: The small number of epochs does not allow a full kinematic analysis. However, using the two available images, we tentatively constrain the jet speed to be likely subluminal or at most mildly superluminal.

• PKS 1718–649: TANAMI provides the first multi-epoch pc-scale maps of the source. Our kinematic analysis yields a rough estimate of the age of this young radio object, which is on the order of 70 years. This is the first kinematic age measurement of this archetypal CSO. Moreover, a spectral index map between 8.4 GHz and 22.3 GHz suggests that the core of this young radio source is strongly absorbed at these frequencies.

In order to study the connection between pc-scale properties and high-energy emission in radio galaxies in a more general fashion, we included public results from the MOJAVE sources, which resulted in a total sample of fifteen sources with VLBI monitoring and detected by Fermi-LAT, the largest sample of γ-ray detected radio galaxies studied with VLBI techniques so far. We find that the VLBI core flux density correlates with the observed γ-ray flux, as observed in blazars, while the γ-ray luminosity does not correlate with typical Doppler boosting indicators such as core brightness temperature and core dominance. This indicates that while the compact pc-scale emission does drive the observed high-energy flux, the observed γ-ray luminosity is not driven by Doppler boosting effects, as is observed in blazars. This difference reinforces the orientation-based unified model of jetted AGN, since radio galaxies are not expected to be Doppler boosting-dominated, having misaligned jets w.r.t. our line of sight.

The TANAMI program is still ongoing, as new observing epochs are added and analyzed. This will allow us to obtain more robust kinematic results, including a first measurement of (or limit to) the separation speed between the two lobes of PKS 1718−649, and further investigation of the new VLBI component possibly associated with the LAT detection of Pictor A. Fermi-LAT monitoring will also continue to provide γ-ray data in the coming years.

In the second paper of this series, we will report on the VLBI properties and updated γ-ray flux upper limits of LAT-undetected TANAMI radio galaxies. We will compare the two subpopulations, and also investigate in a more general fashion the possible connections between VLBI and high-energy properties in radio galaxies.

Acknowledgments

We thank Laura Vega García for the development of the Python GUI-based code that was used for the kinematic analysis and the spectral index maps. We thank the anonymous journal referee, the Fermi-LAT internal referee, Bindu Rani, the MPIfR internal referee, Rocco Lico, as well as Marcello Giroletti, David Thompson and Nicola Omodei for useful comments which improved the manuscript. This research has made use of the NASA/IPAC Extragalactic Database (NED), which is operated by the Jet Propulsion Laboratory, California Institute of Technology, under contract with the National Aeronautics and Space Administration; data from the MOJAVE database that is maintained by the MOJAVE team (Lister et al. 2018); APLpy, an open-source plotting package for Python (Robitaille & Bressert 2012); Astropy, (http://www.astropy.org) a community-developed core Python package for Astronomy (Astropy Collaboration 2013, 2018) C.M. acknowledges support from the ERC Synergy Grant “BlackHoleCam: Imaging the Event Horizon of Black Holes” (Grant 610058). R.S. gratefully acknowledge support from the European Research Council under the European Union’s Seventh Framework Programme (FP/2007–2013)/ERC Advanced Grant RADIOLIFE-320745. The Fermi-LAT Collaboration acknowledges generous ongoing support from a number of agencies and institutes that have supported both the development and the operation of the LAT, as well as scientific data analysis. These include the National Aeronautics and Space Administration and the Department of Energy in the United States; the Commissariat á l’Energie Atomique and the Centre National de la Recherche Scientifique/Institut National de Physique Nucléaire et de Physique des Particules in France; the Agenzia Spaziale Italiana and the Istituto Nazionale di Fisica Nucleare in Italy; the Ministry of Education, Culture, Sports, Science and Technology (MEXT), High Energy Accelerator Research Organization (KEK), and Japan Aerospace Exploration Agency (JAXA) in Japan; and the K. A. Wallenberg Foundation, the Swedish Research Council, and the Swedish National Space Board in Sweden. Additional support for science analysis during the operations phase is gratefully acknowledged from the Istituto Nazionale di Astrofisica in Italy and the Centre National d’Etudes Spatiales in France. The Long Baseline Array is part of the Australia Telescope National Facility which is funded by the Commonwealth of Australia for operation as a National Facility managed by CSIRO. This study made use of data collected through the AuScope initiative. AuScope Ltd is funded under the National Collaborative Research Infrastructure Strategy (NCRIS), an Australian Commonwealth Government Programme. This work made use of the Swinburne University of Technology software correlator, developed as part of the Australian Major National Research Facilities Programme. This work was supported by resources provided by the Pawsey Supercomputing Centre with funding from the Australian Government and the Government of Western Australia. Hartebeesthoek Radio Astronomy Observatory (HartRAO) is a facility of the National Research Foundation (NRF) of South Africa. This research was funded in part by NASA through Fermi Guest Investigator grants NNH10ZDA001N, NNH12ZDA001N, and NNH13ZDA001N-FERMI (proposal numbers 41213, 61089, and 71326, respectively). This research was supported by an appointment to the NASA Postdoctoral Program at the Goddard Space Flight Center, administered by Universities Space Research Association through a contract with NASA. This work was performed in part under DOE Contract DE-AC02-76SF00515. We would like to thank the Institute of Radio Astronomy and Space Research, AUT University, New Zealand for the use and operational support of their radio telescopes in the collection of data for this work.

References

Appendix A: Full resolution VLBI maps and image parameters

Here we present the full-resolution VLBI images for the sources studied in this paper, along with the tables including the details of each observation, for each source.

Fig. A.1. Full-resolution images of Pictor A. The map parameters for each epoch can be found in Table A.1. The blue ellipse represents the beam size, while the black line indicates the linear scale at the source’s redshift. Contours increase in steps of two starting from 1.8, 3.0, 3.5, 5.5, 4.0, 2.0 times the noise level in each map, from top left to bottom right, respectively.

Fig. A.2. Full-resolution images of PKS 0521−36. The map parameters for each epoch can be found in Table A.2. The blue ellipse represents the beam size, while the black line indicates the linear scale at the source’s redshift. Contours increase in steps of two starting from 4.0, 3.0, 3.0, 3.0, 3.0, 5.0 times the noise level in each map, from top left to bottom right, respectively.

Fig. A.2. continued. Contours increase in steps of two starting from 10.0, 4.0, 10.0, 1.0 times the noise level in each map, from top left to bottom right, respectively.

Fig. A.3. Full-resolution images of PKS 0625−35. The map parameters for each epoch can be found in Table A.3. The blue ellipse represents the beam size, while the black line indicates the linear scale at the source’s redshift. Contours increase in steps of two starting from 2.0, 2.0, 2.0, 4.0, 3.0, 4.0 times the noise level in each map, from top left to bottom right, respectively.

Fig. A.3. continued. Contours increase in steps of two starting from 4.0, 2.0, 4.0, 2.0 times the noise level in each map, from top left to bottom right, respectively.

Fig. A.4. Full-resolution images of Centaurus B. The map parameters for each epoch can be found in Table A.4. The blue ellipse represents the beam size, while the black line indicates the linear scale at the source’s redshift. Contours increase in steps of two starting from 9.0 times the noise level in each map.

Fig. A.5. Full-resolution images of PKS 1718−649. The map parameters for each epoch can be found in Table A.5. The blue ellipse represents the beam size, while the black line indicates the linear scale at the source’s redshift. Contours increase in steps of two starting from 5.0, 5.0, 3.0, 10.0, 10.0, 5.0 times the noise level in each map, from top left to bottom right, respectively.

Fig. A.5. continued. Contours increase in steps of two starting from 20.0 and 3.0 times the noise level.

Appendix B: Extended kinematic analysis information

Here we include additional information illustrating the results of our kinematic analysis of the multi-epoch TANAMI data. Figures B.1 through B.4 show the multi-epoch images of TANAMI radio galaxies, with the corresponding component identification and tracking. We note that, for ease of representation, the distance between the images at different epochs is not to scale in these figures. Moreover, the colored lines are not fits to the displayed component positions, but simple interpolations meant to clarify the identification and tracking of the different components.

Tables B.1 through B.4 list the flux density, radial distance, position angle and size for each circular Gaussian component identified in each source during the kinematic analysis. All components have been shifted so that the core position is always at (0, 0) coordinates in all epochs. The position angle is given in the range (π, −π), with the zero in the N–S direction (in image coordinates) and positive values in the counter-clockwise direction. Components that were not identified (blue crossed circles in Figs. B.1 through B.4) are not listed. An apparent speed was fit only for components detected in at least five epochs.

Fig. B.1. Multi-epoch tapered images of Pictor A. The colored crossed circles represent the circular Gaussian components that have been fit to the clean maps. The distance between the images at different epochs is not to scale. The colored lines are not fits to the displayed component positions, but simple interpolations meant to guide the eye.

Fig. B.2. Multi-epoch tapered images of PKS 0521−36. The colored crossed circles represent the circular Gaussian components that have been fit to the clean maps. The distance between the images at different epochs is not to scale. The colored lines are not fits to the displayed component positions, but simple interpolations meant to guide the eye.

Fig. B.3. Multi-epoch images of PKS 0625−35. The colored crossed circles represent the circular Gaussian components that have been fit to the clean maps. The distance between the images at different epochs is not to scale. The colored lines are not fits to the displayed component positions, but simple interpolations meant to guide the eye.

Fig. B.4. Multi-epoch images of PKS 1718−649. The colored crossed circles represent the circular Gaussian components that have been fit to the clean maps. The distance between the images at different epochs is not to scale. The colored lines are not fits to the displayed component positions, but simple interpolations meant to guide the eye.

Fig. B.5. Flux density evolution of the modeled jet features with time. Top left to bottom right: Pictor A, PKS 0521−36, PKS 0625−35, PKS 1718−649.

All Tables

All Figures

Fig. 1. First-epoch 8.4 GHz contour maps of γ-ray detected TANAMI radio galaxies. The blue ellipse indicates the convolving beam, while the black bar shows the linear scale. Top left to bottom right: Pictor A, PKS 0521−36, PKS 0625−35, Centaurus B, PKS 1718−649 (the B1950 name is indicated in the top right corner of each image). The full set of multi-epoch images and the associated map parameter tables are presented in Appendix A.

In the text

	Fig. 2. Spectral index map of Pictor A (0518−458) between 8.4 GHz and 22.3 GHz for epoch 2008 Nov. 27. Black contours are from the 8.4 GHz image. The convolving beam is represented in gray in the lower-left corner.
In the text

	Fig. 3. Position angle of Gaussian jet components for Pictor A (0518−458) versus radial distance from the core. The color coding indicates different components, and is the same as the upper left panel of Fig. 7. The sign of the counter-jet component CJ1 has been reversed to improve visualization.
In the text

	Fig. 4. Spectral index map of PKS 0521−36 (0521−365) between 8.4 GHz and 22.3 GHz for epoch 2008 Mar. 28. Black contours are from the 8.4 GHz image. The convolving beam is represented in gray in the lower-left corner.
In the text

	Fig. 5. Spectral index map of PKS 0625−35 (0625−354) between 8.4 GHz and 22.3 GHz for epoch 2008 Nov. 27. The black contours are from the 8.4 GHz image. The convolving beam is represented in gray in the lower-left corner.
In the text

	Fig. 6. Spectral index map of PKS 1718−649 between 8.4 GHz and 22.3 GHz for epoch 2008 Feb. 07. The black contours are from the 8.4 GHz image. The convolving beam is represented in gray in the lower-left corner.
In the text

	Fig. 7. Jet kinematics of our radio galaxies: core distance of jet features as a function of time. The solid lines represent a least squares fit to their positions (the slope is the apparent speed). Top left to bottom right: Pictor A, PKS 0521−36, PKS 0625−35, PKS 1718−649. The B1950 name is given above each image.
In the text

	Fig. 8. Jet kinematics of PKS 0521−36 (0521−365): core distance of jet features as a function of time, including the previous VLBI dataset of Tingay & Edwards (2002). The solid lines represent a least squares fit to their positions (the slope is the apparent speed). Compare with top-right panel of Fig. 7.
In the text

	Fig. 9. Multi-epoch images of Centaurus B (1343−601). The colored crossed circles represent the circular Gaussian components that have been fit to the clean maps.
In the text

Fig. 10. Parameter space of intrinsic jet speed β and viewing angle θ allowed by our observations. The blue shaded area is the one allowed by the measurement of R (θ as function of β given R, that is, θR), while the red shaded area is the one allowed by the observed βapp (θ as function of β given βapp, that is, θβapp). For each source we provide a minimum, maximum and (except for PKS 0521−36) a central estimate of R and βapp. The top-right legend reports the resulting limits on θ and β. The dashed colored lines and boxes indicate constraints from previous works, namely Hardcastle et al. (2016), Tingay et al. (2000) for Pictor A (0518−458), Pian et al. (1996), Giroletti et al. (2004), D’Ammando et al. (2015) for PKS 0521−36 (0521−365), and Venturi et al. (2000) for PKS 0625−35 (0625−354).

In the text

	Fig. 11. Light curve of Pictor A (0518−458) between 0.1 and 300 GeV over 103 months of Fermi-LAT data. Blue points are detections, red arrows are upper limits. Top to bottom: 1-month, 3-months and 6-months binning, respectively.
In the text

	Fig. 12. Light curve of PKS 0521−36 (0521−365) between 0.1 and 300 GeV over 103 months of Fermi-LAT data, with weekly (blue points) and monthly (red points) binning. Upper limits are indicated by arrows of the respective colors. The left y-axis reports the weekly flux values, the right one reports the monthly flux values.
In the text

	Fig. 13. Light curve of PKS 0521−36 (0521−365) between 0.1 and 300 GeV over 103 months of Fermi-LAT data, with 6-h (blue points) and daily (red points) binning. Upper limits are indicated by arrows of the respective colors. The left y-axis reports the 6-h flux values, the right one reports the daily flux values.
In the text

	Fig. 14. Light curve of PKS 0625−35 (0625−354) between 0.1 and 300 GeV over 103 months of Fermi-LAT data, with monthly binning. Blue points are detections, red arrows are upper limits.
In the text

	Fig. 15. Light curve of Centaurus B between 0.1 and 300 GeV over 103 months of Fermi-LAT data, with monthly binning. Blue points are detections, red arrows are upper limits.
In the text

	Fig. 16. Fermi-LAT map of excess TS in the inner region of the ROI centered on PKS 1718−649, after removing the unidentified catalog source 3FGL J1728.0−6446. The ellipses represent the 68% and 95% confidence positional errors on the catalog source, while the black circles represent the same errors for the new source PS J1724.2−6459. The map size is 1.6 ° ×1.6°. Each pixel corresponds to 0.1°.
In the text

	Fig. 17. Maximum measured apparent speed βapp as a function of median total VLBI luminosity for our TANAMI γ-ray radio galaxies (blue circles) and MOJAVE γ-ray radio galaxies (red squares). The downward arrows indicate upper limits on the apparent speed, for sources with speed consistent with zero within the uncertainties.
In the text

	Fig. 18. Fermi-LAT flux above 100 MeV as a function of median radio VLBI core flux density for our TANAMI γ-ray radio galaxies (blue circles) and MOJAVE γ-ray radio galaxies (red squares).
In the text

	Fig. 19. Fermi-LAT flux above 100 MeV as a function of median radio VLBI jet flux density for our TANAMI γ-ray radio galaxies (blue circles) and MOJAVE γ-ray radio galaxies (red squares).
In the text

	Fig. 20. Fermi-LAT luminosity as a function of median radio VLBI core luminosity for our TANAMI γ-ray radio galaxies (blue circles) and MOJAVE γ-ray radio galaxies (red squares).
In the text

	Fig. 21. Median radio VLBI core brightness temperature as a function of Fermi-LAT luminosity above 100 MeV for our TANAMI γ-ray radio galaxies (blue circles) and MOJAVE γ-ray radio galaxies (red squares). A correction factor has been applied to the MOJAVE data to account for the difference in frequency, but this might still imply a bias.
In the text

	Fig. 22. Median radio VLBI core dominance as a function of Fermi-LAT luminosity above 100 MeV for our TANAMI γ-ray radio galaxies (blue circles) and MOJAVE γ-ray radio galaxies (red squares).
In the text

	Fig. A.1. Full-resolution images of Pictor A. The map parameters for each epoch can be found in Table A.1. The blue ellipse represents the beam size, while the black line indicates the linear scale at the source’s redshift. Contours increase in steps of two starting from 1.8, 3.0, 3.5, 5.5, 4.0, 2.0 times the noise level in each map, from top left to bottom right, respectively.
In the text

	Fig. A.2. Full-resolution images of PKS 0521−36. The map parameters for each epoch can be found in Table A.2. The blue ellipse represents the beam size, while the black line indicates the linear scale at the source’s redshift. Contours increase in steps of two starting from 4.0, 3.0, 3.0, 3.0, 3.0, 5.0 times the noise level in each map, from top left to bottom right, respectively.
In the text

	Fig. A.2. continued. Contours increase in steps of two starting from 10.0, 4.0, 10.0, 1.0 times the noise level in each map, from top left to bottom right, respectively.
In the text

	Fig. A.3. Full-resolution images of PKS 0625−35. The map parameters for each epoch can be found in Table A.3. The blue ellipse represents the beam size, while the black line indicates the linear scale at the source’s redshift. Contours increase in steps of two starting from 2.0, 2.0, 2.0, 4.0, 3.0, 4.0 times the noise level in each map, from top left to bottom right, respectively.
In the text

	Fig. A.3. continued. Contours increase in steps of two starting from 4.0, 2.0, 4.0, 2.0 times the noise level in each map, from top left to bottom right, respectively.
In the text

	Fig. A.4. Full-resolution images of Centaurus B. The map parameters for each epoch can be found in Table A.4. The blue ellipse represents the beam size, while the black line indicates the linear scale at the source’s redshift. Contours increase in steps of two starting from 9.0 times the noise level in each map.
In the text

	Fig. A.5. Full-resolution images of PKS 1718−649. The map parameters for each epoch can be found in Table A.5. The blue ellipse represents the beam size, while the black line indicates the linear scale at the source’s redshift. Contours increase in steps of two starting from 5.0, 5.0, 3.0, 10.0, 10.0, 5.0 times the noise level in each map, from top left to bottom right, respectively.
In the text

	Fig. A.5. continued. Contours increase in steps of two starting from 20.0 and 3.0 times the noise level.
In the text

	Fig. B.1. Multi-epoch tapered images of Pictor A. The colored crossed circles represent the circular Gaussian components that have been fit to the clean maps. The distance between the images at different epochs is not to scale. The colored lines are not fits to the displayed component positions, but simple interpolations meant to guide the eye.
In the text

	Fig. B.2. Multi-epoch tapered images of PKS 0521−36. The colored crossed circles represent the circular Gaussian components that have been fit to the clean maps. The distance between the images at different epochs is not to scale. The colored lines are not fits to the displayed component positions, but simple interpolations meant to guide the eye.
In the text

	Fig. B.3. Multi-epoch images of PKS 0625−35. The colored crossed circles represent the circular Gaussian components that have been fit to the clean maps. The distance between the images at different epochs is not to scale. The colored lines are not fits to the displayed component positions, but simple interpolations meant to guide the eye.
In the text

	Fig. B.4. Multi-epoch images of PKS 1718−649. The colored crossed circles represent the circular Gaussian components that have been fit to the clean maps. The distance between the images at different epochs is not to scale. The colored lines are not fits to the displayed component positions, but simple interpolations meant to guide the eye.
In the text

	Fig. B.5. Flux density evolution of the modeled jet features with time. Top left to bottom right: Pictor A, PKS 0521−36, PKS 0625−35, PKS 1718−649.
In the text