Gamma-ray emission in radio galaxies under the VLBI scope -- I. Parsec-scale kinematics and high-energy properties of $\gamma$-ray detected TANAMI radio galaxies

In the framework of the TANAMI multi-wavelength and VLBI monitoring, we study the evolution of the parsec-scale radio emission in radio galaxies in the southern hemisphere and their relationship to the $\gamma$-ray properties. In this first paper, we focus on Fermi-LAT-detected sources. We perform a kinematic analysis for five $\gamma$-ray detected radio galaxies using multi-epoch 8.4 GHz VLBI images, deriving limits on intrinsic jet parameters. We analyzed Fermi-LAT data in order to study possible connections between the $\gamma$-ray properties and the pc-scale jets of Fermi-LAT-detected radio galaxies, both in terms of variability and average properties. We discuss the individual source results and draw preliminary conclusions on sample properties including published VLBI results from the MOJAVE survey, with a total of fifteen sources. We find that the first $\gamma$-ray detection of Pictor A might be associated with the passage of a new VLBI component through the radio core. For the peculiar AGN PKS 0521-36, we detect subluminal parsec-scale jet motions, and we confirm the presence of fast $\gamma$-ray variability in the source down to timescales of 6 hours. We robustly confirm the presence of significant superluminal motion, up to $\beta_{app}\sim$3, in the jet of the TeV radio galaxy PKS 0625-35. Finally, we place a lower limit on the age of the Compact Symmetric Object (CSO) PKS 1718-649. We draw some preliminary conclusions on the relationship between pc-scale jets and $\gamma$-ray emission in radio galaxies. We find that the VLBI core flux density correlates with the $\gamma$-ray flux, as seen in blazars. On the other hand, the $\gamma$-ray luminosity does not show any dependence on the core brightness temperature and core dominance, two indicators of Doppler boosting, suggesting that $\gamma$-ray emission in radio galaxies is not driven by orientation-dependent effects.


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, radioloud AGN with jets oriented at small angles (θ < 10 • ) to the observer's line of sight . 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. 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 longterm 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 ). 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 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 0 = 73 km s −1 Mpc −1 , Ω m = 0.27, Ω Λ = 0.73 (Komatsu et al. 2011), and the radio spectral indices refer to the convention S ∝ ν +α .

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) 2 , 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 3 . 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  (Acero et al. 2015), unless otherwise indicated. (c) Originally misclassified as BL Lac, this source has been classified as a young radio galaxy based on multi-wavelength studies (Müller et al. 2014a(Müller et al. , 2015. (d) Goldoni et al. (2016). (e) First γ-ray detection reported by Migliori et al. (2016).
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 radioloud 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(Müller et al. , 2015Goldoni et al. 2016;Krauß et al. 2018). Therefore these two sources are not included in this study.

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 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 selfcalibrated 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 A148, page 3 of 35 A&A 627, A148 (2019) RA and Dec, to derive the two components of the velocity vector, following Lister et al. (2009): 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: 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).

Fermi-LAT data
The LAT is a pair-conversion telescope, launched on 2008 June 11 as the main scientific instrument on board the Fermi Gammaray 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 A148, page 4 of 35 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 2 log(L/L 0 ) where L is the likelihood of the model with a point source at the target position, and L 0 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 powerlaw 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.

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.

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) and22.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.
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: 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 • .  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 quasisimultaneous 8.4 GHz (epoch 2008 Mar. 28) and 22.3 GHz (epoch 2008 Mar. 26) images, which is presented in Fig. 4.

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) and22.3 GHz (epoch 2008 Nov. 29), which is presented in Fig. 5. The index is relatively steep even in the brightest region.
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 5 yr (Tingay et al. 1997).
TANAMI provides the first multi-epoch data set for this wellstudied 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.

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.

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 jetto-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 crossidentify 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.
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 −1 , 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.

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).   0.5 ± 0.2 1.9 ± 0.9 1950 ± 23 8 J6 0.0665 ± 0.0009 0.242 ± 0.003 * 8 0.25 ± 0.07 0.9 ± 0.3 12 ∼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 −6 photons cm −2 s −1 . The sharpest flux change is the decrease from this bin to the following 6-h bin, with a variation of a factor ∼3.

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.
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).  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)   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.  Notes. 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.

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.

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 4 (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 lineof-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 multiwavelength 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   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 . 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 luminos-ity, 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 5 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

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 galax-ies 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. (2013Lister et al. ( , 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 24 W Hz −1 ) 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 highluminosity, 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 A148, page 15 of 35   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  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.

TANAMI MOJAVE
γ-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.
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.
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 −5 . When accounting for the redshift as a  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).  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). 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 L VLBI core ∝ D 2 L S VLBI core , the γ-ray luminosity is integrated over several orders of magni-  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).  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. tude in energy, and therefore depends not only on flux and luminosity distance, but also on photon index.

TANAMI MOJAVE
In Fig. 21 we show the median VLBI core brightness temperature versus γ-ray luminosity. The brightness temperature is calculated as  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).

TANAMI MOJAVE
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): 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, CD = S VLBI core /S VLBI tot . 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.
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.

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 pcscale 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 LATundetected 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.

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.    Obs. date Array configuration (a) S total   Obs. date Array configuration (a) S total  Obs. date Array configuration (a) S total

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.