Parsec-scale Properties of eight Fanaroff-Riley Type 0 Radio Galaxies

We report the high-resolution radio observations of eight Fanaroff-Riley type 0 radio galaxies (FR 0s), selected from the published FR 0 sample. These observations were carried out with the Very Long Baseline Array (VLBA) and European VLBI Network (EVN) at frequencies of 5 and 8 GHz with a highest resolution of 0.6 mas. All eight sources show compact structures on projected physical sizes of 0.3-10 parsec. Six sources show a two-sided structure and two sources show a one-sided jet structure. J1025+1022 shows an X-shaped jet structure, which could result from a reorientation of the jet axis due to a restart of the central engine or a projection of a highly curved inner jet, but more studies are needed to examine these scenarios. Proper motions for 22 jet components of the eight sources are determined to be between -0.08 c and 0.51 c. Although most of the sources exhibit flat spectra, other observed characteristics, such as, low-amplitude flux density variations, low jet proper motion speeds, and symmetric two-sided jet structures, tend to support that the pc-scale FR 0 jets are mildly relativistic with lower bulk Lorentz factors and larger viewing angles.


INTRODUCTION
The first observations of radio galaxies were made in the 1950s (Jennison & Das Gupta 1953), and their radio emission structure is often far beyond the physical extent of their host galaxies (Burbidge 1956). The development of the radio interferometers has allowed us to identify radio sources with radio galaxies and quasars (e.g. Edge et al. 1959) and to image the large-scale jet structures in detail (e.g. MacDonald et al. 1968). Radio galaxies were first classified into two distinct populations, namely Fanaroff-Riley (FR) type I and II, according to their luminosity and large-scale morphology (Fanaroff & Riley 1974). FR Is have relatively prominent inner jets and hotspots but diffuse emission at their edges (edge-darken) (e.g., 3C 31, Laing et al. 2008). On the contrary, FR II are more powerful and exhibit two symmetric edge-brightening lobes spanning a scale of several hundred kiloparsecs (Carilli & Barthel 1996). Powerful radio galaxies have two primary evolutionary tracks, corresponding to the formation of high-jet-power sources (FR IIs) and low-jet-power sources (FR Is). Each radio source starts as a Compact Symmetric Object (CSO), and grows into a Medium-sized Symmetric Object (MSO), and then a fraction of MSOs may successfully evolve into Large Symmetric Objects (LSOs) including FR II and I; the end ★ E-mail: xcheng@kasi.re.kr (KASI) † E-mail: antao@shao.ac.cn product, either FR II or I, depends on the initial jet power, the duration of nuclear activity, re-currency of the radio activity, and the interactions between the jet and the interstellar medium (Kunert-Bajraszewska et al. 2010;An & Baan 2012). This FR classification reveals a morphological sequence that is related to the ability of a jet to transport the momentum and energy to a large distance (An & Baan 2012).
Recent radio survey observations have found that most detected sources are associated with low-redshift early-type galaxies, and their typical observational characteristic in radio band is the lack of large-scale extended structure (Best & Heckman 2012;Sadler et al. 2014). In general, their Faint Images of the Radio Sky at Twenty-Centimeters (FIRST) made with the Very Large Array (VLA; Becker et al. 1995) show unresolved images with a resolution of 5 (corresponding to a projected linear size of 5 kpc) (Baldi et al. 2018b). The deficit of extended radio emission and high core dominance of these sources distinguish them from the classical FR I and FR II galaxies. These sources seem to a stand-alone class, and therefore named as a separate class, FR type 0 radio galaxies (FR 0s) (Ghisellini 2011). FR 0s have similar host galaxy and central engine properties with FR Is, but with lower radio power and higher core dominance. FR 0s are the dominant radio source population in the local Universe, and the number density of FR 0s is 5 times higher than that of FR Is at < 0.05. Therefore FR 0s are important for understanding the overall nature of extragalactic radio sources.
Currently, the study of FR 0s is very limited to the available observations in the radio band, and in particular, high resolution radio observations are lacking. The physical nature of the absence of extended emission in FR 0s and the relation to other AGN classes (e.g., FR classes, blazars, radio-quiet AGN) are still subjects of debate (Sadler et al. 2014;Baldi et al. 2015Baldi et al. , 2018b. Recently, a larger sample of 18,286 sources (Best & Heckman 2012) was cross-matched with the Sloan Digital Sky Survey (SDSS, York et al. 2000;Abazajian et al. 2009), NRAO VLA Sky Surve (NVSS, Condon et al. 1998), and FIRST (Becker et al. 1995) datasets with 104 FR 0s identified (FR0CAT, Baldi et al. 2018b). The comparison between the FR 0 and FR I galaxies suggests that FR 0 jets may have lower bulk Lorentz factors than FR I jets, and that FR 0s have shorter active periods that may be recurrent (Baldi et al. 2018b). However, FR 0s are not just a down-sized version of FR Is. The difference in the observational characteristics between these two classes might be dominated by an intrinsic physical mechanism that reduces the production ability of extended radio emission in FR 0s. The low luminosity and high number density suggest that FR 0s may represent short-lived AGN activity and the short duty cycles are not long enough for a galaxy to develop large-scale radio jets (Sadler 2016;Baldi et al. 2018bBaldi et al. , 2019a. Capetti et al. (2020) explored the recurrent behaviour of FR 0s based on the observations made by the Low-Frequency ARray (LOFAR, van Haarlem et al. 2013) at 150 MHz. They identified 66 FR 0s from the available LOFAR data, and investigated their low-frequency radio properties. Of the 66 FR 0s, 54 are unresolved with sizes of 3-6 kpc, while the remaining twelve sources are in the range of 15-40 kpc. The spectral slopes between 150 MHz and 1.4 GHz show that 20 per cent of their FR 0 sample sources have steep spectra ( 1 < −0.5). These observations suggest that FR 0s and FR Is may represent the two extremes of the jetted source population with a continuous distribution of size and radio luminosity, while FR 0s lie at the lower end of the size and radio power (Capetti et al. 2020).
In addition, Baldi et al. (2019b) randomly selected 18 FR 0s from the FR0CAT and made new observations with the Karl G. Jansky VLA at 1.5, 4.5 and 7.5 GHz with a highest resolution of ∼ 0.3 . Most of these sources remain unresolved, suggesting a compact radio structure on the sub-kpc scales and confirming the general lack of extended radio emission and high core dominance of FR 0s. The size distribution favors the possibility that FR 0s can not be all young radio galaxies that will eventually become FR I/II (Baldi et al. 2019b).
We have studied the radio properties of 14 FR 0s on pc scales from very long baseline interferometry (VLBI) imaging observations (Cheng & An 2018). The diverse distribution of low brightness temperatures, slow structure and flux density variations, compact structures, steep spectra and moderate jet speeds implies that the VLBI-detected FR 0s might be mixed with GHz-peaked spectrum (GPS) and compact steep-spectrum (CSS) sources (Cheng & An 2018). Sadler et al. (2014) also found that FR 0s is a heterogeneous population, in addition to the genuine FR 0s, this population is also mixed with GPS, CSS as well as a minority of beamed sources. Young radio galaxies, GPS and CSS radio galaxies (see recent review in O'Dea & Saikia (2021)), must be present in FR 0s (Sadler et al. 2014), but their proportion is not yet well constrained. GPS and CSS galaxies represent an early stage in the evolution of extended radio sources (An & Baan 2012) and carry important information about the onset of nuclear activity. The overall radio properties of the 1 Spectral index is defined as S ∝ .
VLBI-detected FR 0s (Cheng & An 2018) are consistent with what expected for young radio sources, but there are still some open questions to be addressed. Since young radio galaxies have extremely compact radio structures, only VLBI imaging observations are able to explore the pc-scale sources structure. In addition to revealing the pc-scale morphology, VLBI observations can determine the jet/hotspot kinematics and spectral properties providing additional important information for discerning the nature of FR 0s. Following our previous results (Cheng & An 2018), we selected and observed eight FR 0s, which have only one epoch VLBI data, with the Very Long Baseline Array (VLBA) and European VLBI Network (EVN) at 5 and 8 GHz, to study the detailed jet structure, spectral distribution and jet kinematic properties on parsec scales. Of particular interest is the determination of the advanced speed of the hotspot or jet knots, which can be diagnostic of whether the radio structure is still growing.
In this paper, we present radio images of eight FR 0s at 5 and 8 GHz and report the spectral index distribution and jet proper motions. Section 2 describes the sample selection and observations. The observational results are presented in Section 3 and discussed in Section 4. Section 5 gives a summary. Throughout the paper, we assume a cosmological model with H 0 = 73 km s −1 Mpc −1 , Ω M = 0.27, and Ω Λ = 0.73. At a redshift of = 0.01, an angular size of 1 mas corresponds to a projected linear size of 0.196 pc, and the conversion from angular velocity to projected linear speed is 1 mas yr −1 = 0.65 .

Sample selection
Cross-matching of the FR 0s catalogue (Baldi et al. 2018b) with the mJy Imaging VLBA Exploration at 20 cm (mJIVE-20 Deller & Middelberg 2014) catalogue and Astrogeo database 2 identified 14 FR 0s (Cheng & An 2018). Among the 14 sources: (i) three sources have only compact and unresolved cores at 8.4 GHz; (ii) three sources have multi-epoch data with extended jets of 1-6 mas in length, for which we have obtained the jet proper motions and source spectral indices in Cheng & An (2018); (iii) the remaining eight sources have either core-jet structure or core-double-hotspot structure, but only one epoch of VLBI data. Therefore, we selected the remaining eight sources, whose main properties are listed in Table 1, and observed with the VLBA and EVN, so that we have a total of 11 sources for our jet proper motion analysis. We should note that the FR 0 sources selected for this paper is actually the fraction with the highest radio flux density. As we have discussed in Cheng & An (2018), the remaining weaker sources in the FR0CAT may not have similar radio power or morphological characteristics to these 14 sources, and thus the conclusions of this study may not be universally applicable to the entire FR 0 sample. This would require a systematic deep VLBI imaging of a large sample to help clarify the physical nature of the radio emission of FR 0s on pc scales.
Six of the selected sources show a two-sided jet structure. Thus, another motivation of our observations is to identify the core position of the sample and to classify the radio source structure.

Observations
We obtained 10.5 h of observations with the EVN and VLBA, respectively. Based on the same motivation, the schedule of the EVN and VLBA observations are almost identical. We observed the eight sources in total flux density mode at 5 and 8 GHz. To ensure the success of the observations and to optimise the coverage of each source, we split these sources into two groups. The flux densities of these sources are not variable on weeks time scale, so the observations were scheduled to be performed within one week. Each source was observed three to four separate scans, with an effective on-source time of approximately 30 min.
The VLBA observations are summarized in Table 2 and were carried out from 2018 August 4 to August 12 (project codes of BC241A to BC241D; PI: X.-P. Cheng). All 10 VLBA antennas were used for BC241C. The other three sessions used 8 or 9 antennas due to the maintenance of some antennas. The observed data were recorded with two consecutive 128 MHz channels, sampled in 2 bits, at an aggregate data rate of 2 giga bits per second (Gbps). We used 3C 345 and 4C 39.25 as the finger searching calibrators. The raw data were correlated at the Array Operations Centre in Socorro (New Mexico, U.S.) using the DiFX software correlator (Deller et al. 2007(Deller et al. , 2011 with an averaging time of 2 s, 256 frequency channels per IF, and uniform weighting. The EVN observations are summarized in Table 2 and were carried out at 5 and 8 GHz frequency bands from 2019 February 24 to March 7 (project codes of EC063A to EC063D; PI: X.-P. Cheng). Fifteen antennas were used with 2-bit sampling and an aggregate data rate of 1 Gbps (dual-polarization, 8 consecutive 16 MHz frequency channels). Twelve antennas were used at 8 GHz with the same data rate. We used 3C 345 and 4C 39.25 as the calibration sources. The correlation was performed at the Joint Institute for VLBI ERIC (JIVE), Dwingeloo, the Netherlands, using the EVN software correlator (SFXC, Keimpema et al. 2015) .
The correlated VLBA and EVN data were transferred to the computer clusters at the China SKA Regional Centre , where the post-processing including calibration and imaging analysis was performed.

Data Reduction and Model Fitting
The data were calibrated with the U.S. National Radio Astronomy Observatory (NRAO) Astronomical Imaging Processing Software (AIPS) package (Greisen 2003) following the standard procedure. Firstly, we imported the data using the task FITLD and inspected the data quality. The voltage offsets of the samplers were corrected by auto-correlation. A-priori amplitude calibration was performed using the information from the gain curve (GC), system temperature (TY), and weather (WX) tables for all antennas to correct for the atmospheric opacity. The ionospheric delay was corrected using total electron content measurements monitored by the global positioning system (GPS). The phase contributions from the antenna parallactic angles were removed before the subsequent phase corrections were applied. The station Pie Town (PT) and Effelsberg (EF) were chosen as the reference antennas during the VLBA and EVN data calibration process. The fringe-fitting and bandpass calibration were performed with 3C 345 and 4C 39.25. Finally, the data in each sub-band were averaged and exported into single-source files. To remove residual amplitude and phase errors, we carried out the hydrib mapping (Pearson & Readhead 1984) in D software package (Shepherd 1997). Several iterations of self-calibration and imaging were carried out and the final images were made with nat-ural weighting. We found that the overall antennas gain correction factors determined were small, typically within 10 per cent.
A number of circular Gaussian components were fitted to model the brightness distribution of each source in the visibility domain in D . The fitted parameters of the gaussian models are presented in Table B1. Typical flux density uncertainties are less than 10 per cent and are mainly contributed by the visibility amplitude calibration errors and the antenna gain calibration errors. The uncertainty in the fitted component size is typically less than 15 per cent of the deconvolved size of the fitted Gaussian model. We followed Lister et al. (2009) to estimate the errors in the best-fit positions of the Gaussian jet components: i.e., ∼20 per cent of the component size convolved with the synthesised beam size. Since the analysis of the proper motions of the 8 sources contains data obtained from our own EVN and VLBA, as well as the archived VLBA data from the Astrogeo database, to reduce errors due to difference in data quality (in practice, the data quality difference is not much), we used the same model fitting uncertainties for all the data.

Radio Morphology on Parsec-scale
In the first and second columns of Fig. 1, we show the naturally weighted total intensity EVN images of the 8 sources at 5 and 8 GHz. The typical image size displayed is 15 mas × 15 mas and the dimension of the restoring beam in each image is given in Table 3. The VLBA images show similar morphology to the EVN images, so we only present them in Appendix A. The elliptical Gaussian restoring beam is indicated in the bottom-left corner of each image in Figures 1 and A1. The parameters associated with the images are listed in Table 3. The quasi-simultaneous dual-frequency VLBI observations allow estimation of the spectral index (listed in the last column of Table 3), which we calculated using the integrated flux densities on the VLBI images. We found that seven sources show a flat spectrum (−0.5 < < 0.5), with only J1037+4335 exhibiting a slightly steeper spectrum. J1213+5044 shows a change of the spectral index from = 0.19 in 2018 August to = −0.61 in 2019 February, which may be related to a variation in the opacity in the core region due to a change in the emission structure B).
Six sources (J0906+4124, J0909+1928, J1025+1022, J1205+2031, J1213+5044 and J1230+4700) display a two-sided structure on parsec scales, in agreement with the previous results (Cheng & An 2018). Two sources (J1037+4335 and J1246+1153) exhibit a faint one-sided jet structure. In our new observations of J1037+4335, we do not detect the SW component, but detect three new jet components (J1 -J3), which suggest a core-jet structure extending to the northwest (NW) direction. J1246+1153 shows a core-jet structure with the jet extending to the NW direction.
The identification of the core is crucial for kinematic studies, as it provides the stationary reference point for all epochs. In our sample, all sources have a prominent core, which greatly simplifies the identification of the reference component. With the high-resolution data at two frequencies, we can easily get the spectral index of individual components to help the core identification which corresponds to a flat-spectrum component. In most cases, the identification of jet features across epochs proves to be relatively straightforward due to the slow evolution of the source structure with time. We also used an additional check on the continuity of flux density and size to confirm the cross-identification. To avoid misidentification of the component and incorrect estimates of apparent jet speeds, we discuss the detailed characteristics of the radio morphology and components identification for each individual source in Appendix B.

Jet Proper Motion
In order to quantitatively study the jet kinematics, we used a nonacceleration, two-dimensional vector fit to the jet components position with time, using the core as a reference (Lister et al. 2016). Combined with our previous studies (Cheng & An 2018), data from three epochs over a maximum baseline of ∼ 14 years are available for the jet proper motion measurement. One sources (J1037+4335) was unresolved in the previous observations , so its proper motion is determined only based on two epochs of 2018 and 2019. The other seven sources have data from three epochs that allow a linear regression fitting of the jet proper motions. By model-fitting to our VLBI data, we were able to identify 22 jet components of 8 sources across the different epochs.
Column 4 in Fig. 1 shows the radial distance of the jet component versus the observing time. The slope of the linear regression fitting to the temporal change of jet distance gives the proper motion speeds. The corresponding values of the angular speed, apparent speed and observing frequencies are presented in Table B2.
We note that 13 jet components in six sources (J0906+4124, J1037+4335, J1205+2031, J1213+5044, J1230+4700 and J1246+1153) basically move along a constant radial direction. Since the jet components move away from the core, their flux densities tend to decrease in general and their sizes increase. In J0909+2918, we find that the north component (N1) moves along a constant direction, but the south jet components (S1 and S2) exhibit non-radial motion as they bend at 5 mas. During our observation period, S1 shows little motion and little change in position angle. The trajectory of S2 appears a slight eastward bending and a change in the position angle from 140 • to 155 • . In J1025+1022, W3 follows a straight trajectory toward the west direction to the core, but the position angles of W2 and E2 vary significantly from ±90 • to 130 • along the NW-SE direction. However, each of these three jet components seems to move in ballistic trajectories.
Combining with proper motion results of three FR 0s reported in Cheng & An (2018), we obtained jet proper motions of a total of 24 jet components in 11 FR 0 sources to study the jet kinematics. The derived jet speeds range from −0.08 to 0.51 , suggesting stationary motion to mildly relativistic speeds. The only exception of jet proper motion speed exceeding 1, is associated with the J1 component in J1037+4335, which has only two epochs data separated by <7 months. Follow-up VLBI monitoring of this source would place strict constraint on its jet proper motion.
Column 3 of Fig. 1 shows the spectral index maps for all sources between 5 and 8 GHz. We also created the spectral index maps using the VIMAP program (Kim & Trippe 2014) from the 5 and 8 GHz EVN data by aligning the maps on the core positions. The 8-GHz images were created with the same pixel size (0.05 mas) and restoring beam size as the 5-GHz maps. First, we excluded the core region with an elliptical mask approximately twice the size of the restoring beam. Then, VIMAP calculated the two-dimensional cross-correlation product to determine the shift between the two images. Due to the relatively closeness of the two frequencies, the typical shift is about 0.05 mas due to the frequency-dependent opacity effect. After correcting this relative offset, we computed the spectral index maps and superimposed them (in coloured scale) on the 5 GHz total intensity contours. The spectral index maps confirm the flat-spectrum core as the brightest component and reveal twin-jet structure for 6 sources and one-sided jet structure for 2 sources. The outer jet region has a steeper spectrum with a spectral index of < −0.5, typical for optically-thin AGN jets (Hovatta et al. 2014).

DISCUSSION
Several low frequency observations have confirmed the absence of extended radio emission from FR 0s, which is the main difference with other FR classification sources (Baldi et al. 2019b;Capetti et al. 2020). Therefore, the physical properties of FR 0s and their possible relation to large-scale FR Is depend on high-resolution observations on pc scales. Our previous study has shown that FR 0s detected with VLBI could not be a homogeneous population of radio sources (Cheng & An 2018). Thanks to the new EVN and VLBA data, we obtain the spectral index information for 12 FR 0s and jet kinematic properties for 11 FR 0 sources. The VLBI data allow us to study the physical properties of FR 0s at the onset phase of their long-term evolution. The cores other than J1230+4700 occupy more than 50 per cent of the integrated VLBI flux density, and the core dominance is not directly related to the radio structure: that is, sources with a onesided jet do not show a higher core flux density fraction than those with two-sided jets, although the sample in this paper is too small to draw statistically significant conclusions yet. Although most of the cores exhibit flat radio spectra and core dominance, it is still not sufficient to conclude that these cores have strong Doppler boosting: most jets exhibit two-sided symmetric structures, indicating rather large viewing angles; the variability of the cores is considerably low, below 30 per cent for all radio cores except for J1213+5044 which has a variability of 40 per cent at 8 GHz; as mentioned earlier, FR 0 sources have lower jet bulk Lorentz factors than FR Is. These characteristics distinguish these FR 0s from quasars which contain highly relativistic jets and prominent variability. The observed flat radio spectrum may be due to the fact that the nuclear region is still active, with a constant injection of fresh relativistic electrons. This also suggests that FR 0s (at least the brightest fraction in radio) are still growing. The number density of FR 0s in local Universe is much higher than that of FR Is (Baldi et al. 2018b), implying that most FR 0s cannot grow into extended sources with sizes of a few hundred kilo-parsec to Mega-parsec. The successful evolution of an FR 0 galaxy into an extended source depends not only on the intrinsic physical properties of the central engines (the duration and persistence/intermittency of the nuclear activity, initial jet power), but also on the environmental factors such as the density gradient of the ambient medium in the host galaxy, jet-ISM interactions along the jet path and within the terminal lobes (see An & Baan (2012) and references therein).
Interestingly, J1025+1022 shows an X-shaped jet on pc scales. This relatively large jet misalignment (∼ 40 • ) has been observed in other nearby radio galaxies, e.g., 3C 84 (Giovannini et al. 2018) and Mrk 231 (Wang et al. 2021). If the X-shaped jet of J1025+1022 is associated with re-started AGN activity, then the change in the jet direction could be the result of a change in the direction of the spin axis of the black hole (e.g., Ekers et al. 1978;Dennett-Thorpe et al. 2002). The continued evolution depends mainly on the continuity of the nuclear activity and the persistence of the jet (An & Baan 2012). Or the jet misalignment could simply be a projection of a helically twisted jet. Detailed follow-up studies can help to clarify these physical pictures.
Four sources (J0909+1928, J1205+2031, J1230+4700 and J1246+1153) show flat spectra. But their jet velocities, ranging from ∼ 0.04 to about 0.4 , are small compared with the typical flat-spectrum quasars (Lister et al. 2016). From our previous study, 4 other sources (J0910+1841, J0943+3614, J1604+1744, J1606+1814) also show flat spectra but small jet speed (Cheng & An 2018). This suggests that the bulk Lorentz factor Γ of the jet is low for most of the VLBI-detected FR 0s. It is important to note that the VLBI sample we used in this study is at the high flux density end of the overall FR 0 population; since the jets of these sources are found to be mildly relativistic (lower Lorentz factor), it is reasonable to believe that the jets of the majority of FR 0 sources have an even lower Lorentz factor. This is consistent with the conclusion obtained by other groups through other methods (e.g., Sadler et al. 2014;Baldi et al. 2018b;Capetti et al. 2020). These sources do not show any stationary hotspots at the ends of their jets, but the lack of terminal hotspots cannot be used to rule out the possibility that FR 0s represent short-lived episodes of AGN activity that do not last long enough for a galaxy to develop large-scale radio jets. If this scenario is correct, the non-detection of extended relic emission of past radio activity indicates a rapid radiative or adiabatic loss leading to a shorter lifetime of the large-scale lobes.
There is no obvious trend of increasing/decreasing apparent component speed with increasing core distance in these 8 sources, indicating no clear acceleration or deceleration in the pc-scale radio jets. The mildly relativistic jets and the lack of extended radio emission indicate relatively lower magnetic field strength on the pc scales and lower radiative efficiency in the vicinity of the SMBHs. The low jet speeds are more likely a consequence of a lower BH mass and spin than that in FR I/IIs (Baldi et al. 2018b;Garofalo & Singh 2019). Therefore, the extended jets should be dissipative and lie below the surface-brightness threshold of existing large-scale radio surveys (Sadler 2016). Follow-up studies with lower-frequency observations (e.g., Capetti et al. 2020) should be performed on larger samples and should also be expanded to lower radio power FR 0s.
We must point out that in most of our FR 0 sources, neither evidence of deceleration along the jet nor double-double morphology (manifesting re-current AGN activity) has been observed. The present results cannot distinguish between these two possibilities: short-lived recurrent sources and low jet-speed sources associated with low BH mass and spin. Obviously, we can draw a conclusion that the jets are accelerated in the sub-parsec regions from the central engine and remain mildly relativistic on pc scales. Since the value of the BH parameter in FR 0s is probably lower than that of FR I/IIs, the deceleration of the jet flow may occur at large distances (>10 pc), which is beyond the region that the present VLBI images can show. EVN + eMERLIN observations could provide intermediate resolutions to reveal the radio emission on tens of pc scales (Baldi et al. 2018a).

SUMMARY
We have observed 8 FR 0s with the VLBA and EVN at 5 and 8 GHz, and studied the inner jet morphology, spectral indices and proper motions. All eight sources show compact structures. We identified a total of 33 jet components from these sources. The focus of this paper is on jet kinematics and spectral property studies. Combining these results with our earlier VLBI measurements (Cheng & An 2018), we find that: 1. Six sources (J0906+4124, J0909+1928, J1025+1022, J1205+2031, J1213+5044 and J1230+4700) show two-sided structure and two sources (J1037+4335 and J1246+1153) exhibit a one-sided jet structure. We first detect three new jet components in J1037+4335, which suggests a core-jet structure. The X-shaped jetted source J1025+1022 is supposed to be a restarted source. 2. Proper motions of 22 jet components in eight sources are determined, between −0.08 and 0.51 . Most of the sources have a two-sided jet structure indicating a relatively large viewing angle of the jet; the low-level variability and low jet proper motion speeds exclude strong beaming effect and suggest sub-relativistic or mildly relativistic jet flow. These observations tend to support low bulk Lorentz factors for the FR 0 jets. 3. The spectral indices for these sources are distributed over a relatively wide range from −0.72 to 0.29. Most of these sources have flat spectra, indicating that the nuclear region is in an active state with a continuous injection of fresh relativistic particles. The spectrum of J1037+4335 is slightly steeper. The spectral index of J1213+5044 has changed significantly between 2018 and 2019. No correlation is found between the spectral index and the radio structure as well as the core dominance. The sample in this paper is relatively small and the analysis of these correlations requires more data. We created spectral index maps for all eight sources, revealing the spectral index distribution of the core and jet.
We expect to further study the jet kinematics at high spatial resolution in a larger FR 0 sample to enrich our understanding of the radio source nature of FR 0s and the possible evolutionary relation to other radio galaxy classes. Upcoming next-generation radio facilities, such as the Square Kilometre Array (SKA) and the next generation Very Large Array (ngVLA, Murphy et al. 2018), will provide a huge increase in sensitivity (Kapinska et al. 2015), up to several orders of magnitude, spanning the whole radio frequency band, and offer a promising opportunity to discover a large population of FR 0s. Although black-hole spins are difficult to measure observationally, the large-scale environment is associated with BH spins and host galaxy dynamics, which would reflect the intrinsic differences between FR 0s and FR Is and shed light on the low-power radio galaxy evolution. data presented in this publication are derived from the following EVN project code: EC063. The VLBA observations were sponsored by Shanghai Astronomical Observatory through an MoU with the NRAO (Project code: BC241). The Very Long Baseline Array is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc.

DATA AVAILABILITY
The correlation data used in this article are available in the EVN data archive (http://www.jive.nl/select-experiment) and VLBA data archive (https://archive.nrao.edu/archive/ archiveproject.jsp). The calibrated visibility data underlying this article can be requested from the corresponding authors.

APPENDIX A: VLBA IMAGES
The VLBA images of eight FR 0 sources are presented in Figure A1. The images parameters are referred to Tables 3.

APPENDIX B: DESCRIPTION ON INDIVIDUAL SOURCES
Here, we report the detailed image structure and component identification based on the properties of morphology, spectral index and kinematics obtained from the previous results and our observations.

B2 J0909+1928
The source J0909+1928 is associated with a galaxy at redshift = 0.028 3 . The source is unresolved at 150-MHz LOFAR observation and 1.4-GHz NVSS observation (Condon et al. 1998;Capetti et al. 2020). Our 5-and 8-GHz EVN images (second row of Figure 1) exhibit an asymmetric two-sided structure at mas resolution, with one jet extending to the north and the counter jet apparently bending eastward at ∼5 mas (2.5 pc), suggesting a local interaction between the jet and dense interstellar medium (ISM) in the narrow-line region of the host galaxy. Such jet-ISM interaction could deflect the jet flow and change its direction (Kino et al. 2018). The source has a complex radio spectrum: a flat spectrum with a spectral index of = −0.42 between 150 MHz and 1.4 GHz (Capetti et al. 2020), an inverted spectrum with = 0.34 between 1.4 GHz and 5 GHz (Capetti et al. 2019   Notes. Column description: (1) J2000 name; (2) observing date; (3) VLBI array; (4) observing frequency; (5) integrated flux density; (6) peak specific intensity; (7) off-source rms noise in the clean image; (8) major axis of the restoring beam (FWHM); (9) minor axis of the restoring beam (FWHM); (10) position angle of the major axis, measured from north through east; (11) spectral index.
Using the core as a reference, a linear regression fit gives the proper motion velocities of three components: N1 (along the jet axis) 0.219±0.017 mas yr −1 , S1 (almost no motion and position angle changing) −0.005±0.054 mas yr −1 , and S2 ( non-radial motion) 0.104±0.037 mas yr −1 . These angular velocities correspond to 0.39 , −0.08 , 0.18 for N1, S1, S2 respectively.  Table 3. The contours increase in steps of 2. The grey-colored ellipse in the bottom-left corner of each panel denotes the restoring beam.

B3 J1025+1022
The source J1025+1022 is at a redshift of 0.046 (Huchra et al. 2012). The source is unresolved at 150-MHz LOFAR observation and 1.4-GHz NVSS observation (Condon et al. 1998;Capetti et al. 2020). Our observation shows a rich two-sided structure with the jet extending roughly in the northwest-southeast direction on parsec scales. The inner jet follows a straight trajectory along the east-west direction ( ±90 • ) until 3 mas, but the outer jet has a position angle P.A.= ∼ −50 • (or 130 • ) and extends from northwest to southeast. The high-resolution VLA observations reveal a slightly convex radio spectrum: the flux density at 4.5 GHz is typical ∼20 per cent higher than that at 1.5 GHz and 15 per cent higher than that at 7.5 GHz (Baldi et al. 2019b). We calculated the total flux density spec-  (Cheng & An 2018). A new component E3 is identified due to the high EVN resolution in the east-west direction. J1025+1022 could be a restart of the central power source at a different jet axis direction, leading to an X-shaped source. The possible adiabatic expansion of the components E1 and W1 result in the eventual dissipation and disappearance at 8 GHz (An & Baan 2012). Alternatively, the apparent X-shaped morphology is a projection of a helical jet.

B4 J1037+4335
The source J1037+4335 is associated with a quasar at redshift = 0.025 (Falco et al. 1999). The 150-MHz LOFAR image and 1.4-GHz NVSS image only detect a compact core with a flux density of ∼ 260.9 mJy and ∼ 132.2 mJy (Condon et al. 1998;Capetti et al. 2020). J1037+4335 was resolved into two isolated components, separated by 19.7 pc in the previous study (Cheng & An 2018). Our high resolution observations first show a faint one-sided jet extending to the north-west (see Fig. 1). The core component shows the most compact size and the highest flux. Three new components (J1 -J3) are detected. The SW component is not detected in our EVN and VLBA observations at 5 and 8 GHz, although our sensitivity is higher (1 =0.12 mJy beam −1 at 8 GHz) compared to the previous image (Cheng & An 2018). Considering the large distance from the core and almost being perpendicular to the jet, we suggest that the SW component might be an artifact. The total flux density spectrum is steep ( 0.15−1.4 = −0.30, 1.4−5 = −0.56) at low frequencies (Capetti et al. 2020). The total flux density spectrum is also steep ( 4.87−8.37 VLBA = −0.58, 4.99−8.42 EVN = −0.72) between 5 and 8 GHz with the quasi-simultaneous two frequency observations, suggesting an emission peak at < 1.4 GHz. The flux density variation is less than 10 per cent from the three epochs of 8.4 GHz VLBI data.
The source is unresolved in the previous data, so we used the data from only two new epochs to determine the proper motion. Using the core as a reference, a linear regression fit yields the proper motion velocities of the three components: for J3 0.178±0.465 mas yr −1 , for J2 0.249±0.520 mas yr −1 , and for J1 0.783±0.991 mas yr −1 . These angular velocities correspond to 0.30 , 0.42 , and 1.32 for J3, J2, and J1 respectively. We note that we measured the separation velocities between two epochs on a small time baseline, so the uncertainties in the estimated values are high. A longer time span is needed to better constrain the proper motion.

B6 J1213+5044
The source J1213+5044 is associated with the galaxy NGC 4187 at redshift = 0.031 (Falco et al. 1999). The 150-MHz LOFAR image and 1.4-GHz NVSS image only detected a compact core with a flux density of ∼ 178.1 mJy and ∼ 96.5 mJy, respectively (Condon et al. 1998;Capetti et al. 2020). J1213+5044 shows a symmetric two-sided structure with jet extending in the west-east direction on parsec scales, in agreement with the previous results (Cheng & An 2018). The total flux density spectrum is steep (  Interestingly, the west component (W) has an optically thick spectral feature with of +1.01±0.36. One possibility is that W is strongly absorbed by the intervening ionised gas. Using the core as a reference, a linear regression fit gives the proper motion velocities of the two components: 0.118±0.033 mas yr −1 for W and 0.044±0.041 mas yr −1 for E. These angular velocities correspond to 0.24 and 0.18 for W and E, respectively.

B7 J1230+4700
The source J1230+4700 is associated with a galaxy at redshift = 0.039 (Mahdavi & Geller 2004). The source is unresolved in 150-MHz LOFAR observation and 1.4-GHz NVSS observation (Condon et al. 1998;Capetti et al. 2020). Our 5-and 8-GHz EVN images show a two-sided structure with jet extending in the northwestsoutheast direction, consistent with the previous results (Cheng & An 2018). The source has a flat spectrum with a spectral index of = −0.14 between 150 MHz and 1.4 GHz (Capetti et al. 2020 Using the core as a reference, the linear regression fitting gives the proper motion velocities of three components: 0.073±0.014 mas yr −1 for S2, −0.088±0.022 mas yr −1 for N, and 0.064±0.027 mas yr −1 for S1. These angular velocities correspond to 0.22 , 0.16 , and 0.21 for S2, N and S1, respectively. The low radial velocity obtained from proper motion measurements suggests an intrinsic low jet speed.

B8 J1246+1153
The source J1246+1153 is at redshift = 0.047 (Binggeli et al. 1985). J1246+1153 is unresolved in 150-MHz LOFAR observation and 1.4-GHz NVSS observation (Condon et al. 1998;Capetti et al. 2020). Our high resolution VLBI observations show a core-jet structure with the jet extending to the northwest direction, consistent with the previous results (Cheng & An 2018). The source has a flat spectrum with a spectral index of = 0.24 between 150 MHz and 1.4 GHz (Capetti et al. 2020 GHz VLBI data shows a very large variation. Using the core as a reference, the linear regression fitting gives the proper motion velocities 0.128±0.101 mas yr −1 for J1, corresponding to an angular velocity of 0.39 . This paper has been typeset from a T E X/L A T E X file prepared by the author. Notes: column description: (1) sources name; (2) VLBI array; (3) observing frequency; (4) component label; (5) integrated flux density of the fitted Gaussian model component; (6) peak intensity; (7) radial distance from the core; (8) position angle respect to the core, measured from north through the east; (9) component size