The eMERLIN and EVN view of FR0 radio galaxies

We present the results from high-resolution observations carried out with the eMERLIN UK-array and the European VLBI network (EVN) for a sample of 15 FR0s, i.e. compact core-dominated radio sources associated with nearby early-type galaxies (ETGs) which represent the bulk of the local radio galaxy population. The 5-GHz eMERLIN observations available for 5 objects exhibit sub-mJy core components, and reveal pc-scale twin jets for 4 out of 5 FR0s once the eMERLIN and JVLA archival visibilities data are combined. The 1.66-GHz EVN observations available for 10 FR0s display one- and two-sided jetted morphologies and compact cores. The pc-scale core emission contributes, on average, to about one tenth of the total extended radio emission, although we note a increasing core contribution for flat/inverted-spectrum sources. We found an unprecedented linear correlation between the pc-scale core luminosity ($\sim$10$^{21.3}$-10$^{23.6}$ W Hz$^{-1}$) and [O III] line luminosity, generally considered as proxy of the accretion power, for a large sample of LINER-type radio-loud low-luminosity active nuclei, all hosted in massive ETGs, which include FR0s and FRIs. This result represents further evidence of a common jet-disc coupling in FR0s and FRIs, despite they differ in kpc-scale radio structure. For our objects and for other FR0 samples reported in the literature, we estimate the jet brightness sidedness ratios, which typically range between 1 and 3. This parameter roughly gauges the jet bulk Lorentz factor $\Gamma$, which turns out to range between 1 and 2.5 for most of the sample. This corroborates the scenario that FR0s are characterized by mildly-relativistic jets, possibly as a result of lower spinning black holes (BHs) than the highly-spinning BHs of relativistic-jetted radio galaxies, FRIs.


Introduction
The presence of accreting supermassive black holes (BHs) (Active Galactic Nuclei, AGN) at the centre of local massive early-type galaxies (ETGs) has been largely reported by large-area surveys (e.g. [22,44,55,109,127,142]) and supported by theoretical studies (e.g. [36,39,128]). The most massive ETGs appear to be typically associated with radio-loud AGN (RL AGN or radio galaxies, RG, in general) able to launch the most powerful jets in the Universe. In opposition to the latter which have typically activity period of 10 7 -10 8 yr, low-power ( 10 23 W Hz −1 ) RGs appear to be steadily active although at low regime [122]. From the most powerful to the weakest, RGs encompass a large range of jet properties, such as luminosities, morphologies, duty cycles and speeds, but all sharing a single type of evolved host [63,97]. This capability of massive ETGs to be frequently associated with radio AGN has been recently object of several observational and theoretical studies and contextualised in the framework of AGN feedback, since RGs, even at low powers, can continuously inject energy in the central kpc of galaxies (e.g. [18,39,58,59,141]). 1 In the optical taxonomy of RL AGN, RL LINERs are typically named as Low Excitation Radio Galaxies (LERGs), while Seyferts as High Excitation Radio Galaxies [24,65]. LINER and LERG are, thus, equivalent optical classifications. Local FR Is are almost totally classified as LINERs, whereas FR IIs are associated with either LINER or Seyfert-like nuclear spectra. For example, In the Third Cambridge (3C) catalogue of RGs [17], FR Is are mostly all LERGs, with only a few possible candidates of FR I HERGs, e.g. 3C 120. limiting the amount of extracted energy available to launch and power the jet, thus resulting in compact, weak jets. FR 0s could be early-phase FR Is and will eventually evolve into the latter over a period of ∼10 8 -10 9 years [6]. FR 0 jets could be highly intermittent, switching on/off frequently during the course of their evolution [4,6,95], resulting in the absence of prominent extended jets. One of the best ways to probe the very inner parts of the ejection mechanism in this mostlyunexplored class of RGs is to study the pc-scale radio emission with very-long baselined imaging observations (VLBI technique), with the intention of finding the smoking-gun feature of their inability of launching fully-fledged jet structure. At parsec scale, higher-resolution radio observations of a sample of FR 0s with world-wide VLBI, the American Very Long Baseline Array (VLBA) and European VLBI Network (EVN) show resolved jets of a few pc for ∼80% of the sample [31,32]. The VLBI multi-epoch data and the symmetry of the radio structures indicates that the jet bulk speeds are mildly relativistic (between 0.08c and 0.51c) with low bulk Lorentz factors (between 1.7 and 6) and large viewing angles. However, all previous VLBI-based studies focus on particularly bright FR 0s (flux densities higher than 50 mJy, a factor 10 higher than the typical FR 0 flux selection threshold [6]). Therefore, a comprehensive study of pc-scale radio emission in a significant and representative sample of FR 0s, moving to lower radio luminosities, is missing. Such a study would help to build an impartial view of the mechanisms of jet production in RGs, linking FR Is to FR 0s.
Together with VLBA and EVN, the eMERLIN also provides a range of baseline lengths that permit unique studies of faint radio sources to be made over a wide range of spatial scales. In this work we will analyse the pc-scale radio emission of a sample of 15 low-power FR 0s with eMERLIN and EVN arrays.
This paper is organized as follows. In Sections 2 and 3 we present the sample, observations and maps of the 5 FR 0s observed with the eMERLIN and 10 with EVN, respectively (Tab. 1). The main properties of the pc-scale emission of FR 0s are presented in Sect. 4 and discussed in Sect. 5 in the context of disc-jet coupling in FR 0s in relation to other classes of RGs. We summarise the results and draw our conclusions in Sect. 6. A supplementary section A provides the eMERLIN and EVN map parameters.

eMERLIN 2.1. Sample, observations and maps
A sample of several thousands RGs selected by cross-correlating the optical SDSS and 1.4-GHz NVSS and FIRST datasets (hereafter the SDSS/NVSS sample) has been recently created [19], which mostly (∼80% of the sample) consists of compact unresolved objects at 5"-scale with FIRST [3]. This sample includes compact RGs with FIRST flux densities > 5 mJy up to z ∼0.3, covering the range of radio luminosity ∼ 10 22 -10 26 W Hz −1 . Based on a JVLA follow-up of 12 sources in L and C band, [4] classified seven FR 0s based on their optical and radio properties at 1.4, 4.5 and 7.5 GHz. We proposed eMERLIN observation of this sample, but only five FR 0s (named eMERLIN FR 0 sample) have been observed in C band (CY4213 project), see Table 1 (upper part).
The five sources are randomly selected from the seven FR 0s studied with the JVLA [4]. Two sources appear extended on a few kpc at sub-arcsec resolution with JVLA: J0020-0028 show a twin 1" structure, probably associated with the dusty disk of an edge-on lenticular galaxy and J0101-0024 displays an 0.8" elongated bent structure. The five FR 0s have FIRST flux densities in the range 5-25 mJy and are characterised by a steep radio spectra at 1.4-4.5 GHz, but a spectral flattening emerges at higher frequencies (towards 7.5 GHz). They are all hosted in red massive ETGs with BH masses >10 7.5 M .
The eMERLIN observations of the five FR 0s were carried in December 2016 with 6 telescopes (without the Lovell telescope) spread across the UK (having a maximum baseline ); (6)-(9) JVLA core luminosity (erg s −1 ), morphology (SR slightly resolved or extended) and size (kpc) and radio spectral type (S steep, F flat, I inverted) with the 1.4-4.5 GHz JVLA spectral slope (S ν ∝ ν α ). length of 217 km), reaching a nominal resolution of ∼40 mas at the C band. The observations were centered at 5 GHz with a bandwidth of 512 MHz (4 spectral windows) and correlated at Jodrell Bank Observatory. The observations were divided in two observing blocks, where two or three targets are stored together with their phase and flux calibrators. We scaled the flux density using ∼15 min scans of 3C 286 and ∼10 min scan of the bandpass calibrator OQ 208 at the beginning of each run. Then we tracked the complex antenna gains using regular ∼2.5 min scans of the bright phase reference source (see list in Tab 1) which we interleaved between 6.5 min scans on the target field, with a total target time-on-source of ∼2 hours. We calibrated the dataset with the eMERLIN CASA data reduction pipeline [96], which is intended to automate the procedures required in processing and calibrating radio astronomy data from the eMERLIN correlator. After removing the low-level RFIs with the AOFlagger software [112], we also performed a standard calibration procedure by fitting the offsets in delays using the CASA task 'fringefit', by calculating the bandpass and phase solutions 2 . The complex gains (phase and amplitude) were improved with several runs of self-calibration on the gain calibrators and then applied to all the targets. After a further inspection of the target data to check the quality and remove possible bad scans, the calibrated datasets were imaged (Stokes I) within the CASA environment using the task 'tclean' with the deconvolver mode 'mtmfs', which performs a multiscale, multi-frequency synthesis algorithm [116], with a cell size of 0.013" and natural weighting. For the targets with flux densities higher than 5 mJy, we carried out a few rounds of self-calibration in phase and a final one in phase and amplitude, using 1-2 min integration times and using a 3σ minimum threshold for valid solution.

Name
To analyse the source parameters in the radio maps, we used 'imfit', part of the CASA 'viewer', which fits two-dimensional Gaussians to an intensity distribution on a region selected interactively on the map, providing the position, the deconvolved size, peak flux density, integrated flux density, and position angle (PA) of the compact source (all listed in Tab. A1). Whilst this procedure is valid for compact components, to estimate the total brightness of an extended component is by interactively marking the region around the irregular shape of the source with CASA 'viewer'. Figure 1 displays the the 5-GHz maps of the five FR 0s with angular resolution (restoring beam) of ∼40 mas and rms between 40 -80 µJy beam −1 .
The five sources are detected with flux densities ranging from 0.4 and 2.7 mJy beam −1 at 5 GHz (see Tab 2). The sources appear all slightly resolved, except J2336+0004 which shows a weak second component with PA ∼ 7 • . The angular lengths, estimated from the convolved-beam major axis, considered as the largest angular scale of the source, are of ∼50 mas, apart from J2336+0004 which is 0.15": the linear physical scales correspond to ∼70-240 pc. All the targets, even the two JVLA resolved sources (J0020-0028 and J0101-0024), do not show significant evidences of jetted extended morphologies at parsec scale.

Combining JVLA and eMERLIN data
Since we own the calibrated C-band JVLA observations of these five sources (carried out in 2012-2013), we combined the calibrated eMERLIN visibilities with those from JVLA in the same band to perform a hybrid map, obtained with long and short baselines. To match the overlapping bandwidth of two datasets, we exported the first two spectral windows of the eMERLIN visibilities (centered at 4816.5 and 4944.5 MHz) and the last two spectral windows of the JVLA visibilities (centered at 4819 and 4947 MHz). First, we performed a first Table 2: Column Description: (1) name; (2)-(5) core peak flux density (mJy beam −1 ), core brightness temperature (units of 10 4 K), total core flux density (mJy) and its core luminosity (erg s −1 ) from eMERLIN maps at 5 GHz; (6)-(7) radio morphology (TS two-sided, OS one-sided, SR slightly resolved, UR unresolved) and size (pc); (8) core peak flux density (mJy beam −1 ) from the JVLA maps [4]; (9)-(10) core peak flux density (mJy beam −1 ) and total emission from the combined JVLA-eMERLIN maps; (11)-(13) PA (degree), radio morphology and size (kpc); (14) jet sidedness. shallow re-weight of each individual visibility by its rms value with the CASA task 'statwt' to contribute equally in the combined dataset given that two radio arrays both have similar image noises and antennas sizes. This process is equivalent to combining the data sets prior to Fourier transformation into the sky-plane. Both the CASA measurement sets files were imaged simultaneously with task 'tclean', using multiscale parameters [0,5,10] to give the best compromise between diffuse structure and unresolved sources. A further issue to solve is that the combined JVLA-eMERLIN beam is not a regular Gaussian due to the presence of wide shoulders owing to the short spacings provided by the JVLA. Therefore, following the procedure discussed by [117], we reduced the possible flux density offsets between the combined maps and the original ones by deconvolving the combined dirty maps to deeper thresholds than those reached in the individual JVLA and eMERLIN maps (see Table for the typically lower rms obtained), and re-weighing the data so that the resultant beam more closely represents a Gaussian. Systematic effects and complications, caused by an extreme deconvolution, can be mitigated with continuous prudent iterations [118]. Figure 2 presents the maps of the combined JVLA-eMERLIN datasets for the 5 FR 0s centered at 4.9 GHz with an angular resolutions less than 0.1" (see Tab A1). Their rms (a few tens of µJy) are lower than that of the eMERLIN maps because deep self-calibration allowed to reduce the noise levels and to reduce the problem of the irregular combined JVLA-eMERLIN beam 3 . All the targets appear extended, apart from J0020-0028 where only a core with scattered components along the radio axis (PA ∼ 45 • ) has been detected 4 . Twosided jetted morphologies with extents ranging (0.2-0.8") ∼0.3-1.5 kpc have been resolved, with core components typically of 1 mJy beam −1 , apart from J2346+0059 which outstands for its brightness, 7.3 mJy beam −1 . The twin jets detected in J2336+0004 have a PA of 46 • , significantly different from what has been measured from the marginally one-sided pc-scale extent (PA ∼7 • ) detected with the eMERLIN array singularly. J0101-0024 displays two-sided structures with lobes, with a slightly convex bending of the jets (PA moving from 38 • to 29 • from the north to the south jets), consistent with the elongated gentle curvature noted in the JVLA map [4].
eMERLIN resolves the core component, probably related to the jet basis, for the five FR 0s with a core dominance F core /F NVSS of typically ∼6 %. The apparent absence of pc-scale jets in the eMERLIN maps is probably due to the intrinsic weakness of their radio brightness. The kpc-scale jetted structures detected in the combined maps consist of 40-80% of the total emission and the pc-scale core components detected in full-resolution with the eMERLIN consist of 10-50 % of the total emission detected in the combined maps.

Sample, observations and maps
A catalogue of 104 bona-fide FR 0s galaxies from the SDSS/NVSS sample, named FR0CAT, have been recently set [6] by including compact sources with: 1) a limit on their deconvolved angular sizes of 4", corresponding to a linear size ∼5 kpc, based on the FIRST images, 2) optical spectrum characteristic of LINERs, 3) red massive ETGs and 4) z<0.05. Eighteen objects were randomly extracted from the FR0CAT sample and observed with the JVLA at L and C band in A-array configuration [8]. To study the pc-scale structure and then increase the chances to resolve the jets, ten out of these eighteen sources (see Tab 1, lower panel) with JVLA flux densities at 1.4 GHz larger than 20 mJy were selected to be observed with the EVN array. These sources have BH masses (10 7.5 -10 9 M ) and JVLA radio core luminosities (∼10 39 -10 41 erg s −1 ) similar to those of the eMERLIN FR 0 sub-sample. Steep (but flattening at 7.5 GHz), flat and inverted spectra in the range 1.4-4.5 GHz characterise the 10 sources [8]. Three sources reveal kpc-scale emission in the JVLA maps: J0907+3257, similar to J0020-0028, exhibit radio emission from an edge-on disc component; and J1213+5044 and J1703+2410 reveal genuine two-sided jets of a few kpc. (2)-(6) core peak flux density (mJy beam −1 ), core brightness temperature (K), total core flux density (mJy) and total emission from the entire structure (mJy), and its core luminosity (erg s −1 ) from EVN maps at 1.66 GHz; (7)-(9) PA (degree), radio morphology (TS two-sided, OS one-sided, SR slightly resolved, UR unresolved) and size (pc); (10) jet sidedness. For J1040+0910 the core location within the structure is ambiguous and we mark with * its core flux densities and luminosity, estimated assuming its most south brightest component as core. The observations of these 10 objects with EVN (EVN FR 0 sample) at 18 cm (1.66 GHz) occurred in May 2019 in phase-reference mode: target source and a nearby phase calibrator (see list in Tab. 1) were observed with a ∼ 5 minutes cycle. Data were correlated at the JIVE correlator 5 in Dwingeloo and the standard EVN pipeline [120] was used for a first calibration. These data were copied in Bologna where a final calibration and data reduction was done using AIPS. UV-data of calibrators and target sources were examined and bad data were flagged. A standard calibration procedure was followed: a gain correlation solutions were applied to the target sources from the phase and flux calibrators. A first image was obtained with the 'imagr' task in AIPS and sources with a correlated flux density of a few mJy were self-calibrated to correct the phase only. Analogously to the method explained in Sect 2.1 for the eMERLIN FR 0 sub-sample, flux densities and sizes for point-like or slightly extended sources are obtained with a Gaussian fit (CASA task 'imfit'), whereas for extended sources the flux density are measured by a flux integration over a region manually selected. Figure 3 depicts the final maps of the 10 FR 0s observed with EVN at 1.66 GHz with a typical angular resolution of 20 mas and rms noise level of several tens of µJy (see Tab.A1). All the targets are detected with total flux densities of several mJy (7.5-144 mJy), where the core component appears to be brightest element of the structure (from one fourth to almost the entire flux density of the total emission). In J0907+3257, which exhibit a disc-related extended emission in the JVLA maps [8], at higher resolution we resolve a triple morphology 6 with a PA, roughly perpendicular to the kpc-scale extended emission. Most of the sources (7/10) appear resolved in elongated structures where additional components are detected in some cases: one-sided jetted (3/7) or double or triple sources (3/7 two-sided). We abstain from the morphological classification of J1040+0910 because it shows a bent triple structure (bending from -12 • to 24 • ), where the core position is ambiguous: in the case of the core as the brightest component of its triple structure, it can be classified as one-sided jetted source, otherwise if the core is the central component, it resembles a two-sided structure. Here we assume that its putative core is the brightest south component. The physical scales of the 10 sources vary between ∼5.4 to 133 pc (∼30 mas to 0.15"). Four of these sources have been already observed with VLBI, VLBA and EVN at 2.3, 5 and 8.4 GHz [31,32] with an angular resolution (up to a factor 10) higher than our observations. J0943+3614 displays a single variable unresolved component with flux densities between 130-270 mJy beam −1 with the VLBI [31], consistent with our EVN detection. While we detect a single core, observations with VLBA and EVN of J1025+1022 from [32] resolved a two-sided jetted shape with a core component of 60 mJy beam −1 , a factor ∼0.7 lower than our detection. J1213+5044 reveals a VLBA and EVN east-west twin structure of a few pc [32], included within the (51 pc long) one-sided (PA ∼-25) morphology we detect. J1230+4700 has been observed with the VLBI, VLBA and EVN [31,32] and shows a compact North-South two-sided morphology of a few pc, roughly aligned with (128 pc long) triple structure we detect.
EVN resolves the pc-scale core component of several mJy for the all 10 FR 0s with a core dominance F core /F NVSS ranging on a large interval of values, 0.03-4, indicative of a wide variety (e.g. aligned, variable sources, see Sect. 5) of radio sources involved in this sample. For the resolved sources, the jet emission consists of 10-60% of the total emission, values which strongly depend on the morphology (one or two-sided).

Results
Starting from their sub-arcsec JVLA compact morphologies, the 5-GHz eMERLIN and 1.66-GHz EVN observations of 15 FR 0s (at z<0.1) have resolved their core component within an extended structure for 11 sources. The central core, which pinpoints the location of the putative BH, represents the jet base, where we can explore the mechanisms of jet launching and propagation at parsec scale. For our FR 0 sample, the core flux densities range from sub-mJy to several tens of mJy beam −1 . More precisely, the EVN FR 0 sample has higher flux densities than the eMERLIN sample due to the selection criteria of the EVN observations. The fraction of the pc-scale core flux density with respect to the JVLA component gradually increases with the radio frequencies 7 , on average, from around ∼10-20% at 1.4 GHz, to 20-30% at 4.5 GHz, to 50-70% at 7.5 GHz. This suggests that a flat-spectrum core emerges from our high-resolution observation, even in the steep-spectrum sources. As exception, the EVN core flux density of J1530+2705 is higher than those detected previously with the JVLA, a hint for further processes involved (e.g. Doppler boosting, accretion/ejection variability).
The core brightness temperature of the eMERLIN sample is a few 10 4 K (Tab. 2), lower than the typical value 10 7 K, generally expected from relativistic particles emitting synchrotron emission (e.g. from jets or AGN), due to their low peak flux densities (sub-mJy beam −1 level). The much brighter cores and a smaller beam of the EVN sample at 1.66 GHz result in higher brightness temperatures (Tab. 3), i.e. 10 7 -10 8 K, but lower than the nominal equipartition brightness temperature ∼10 10 K [119] expected from a relativistic beaming of the jet.
The most remarkable result of this study is the unprecedented concatenation of the eMERLIN and JVLA visibilities for five FR 0s. This procedure turned out to be a positive test to detect sub-kpc scale jets in FR 0s when clearly resolved jetted structure both at arcsec and mas scale resolution do not emerge. Four out of five FR 0s reveal the presence of two-sided jets with radio extents ∼0.3-1.5 kpc.
Since FR 0s have been defined as compact radio sources which lack substantial extended emission, the detection of jetted structures for most of the sources from tens of pc to kpc is surprising. The pc-scale jetted shapes are typically aligned with the JVLA jets or those detected at higher resolution by [31,32]. The radio morphologies of the pc-scale extended jets are predominantly two sided (two-sided:one-sided=7:3) showing lobes and plumes 8 . Among the two-sided jet group, core brightened morphologies are more frequent (corebrightened:edge-brightened=4:3) than the edge-brightened morphologies (triple and lobed), despite of an evident difficulty of classification due to the low brightness of the jets.
Particularly interesting are the two cases of RGs which emit kpc-scale disc emission probably from edge-on lenticular galaxies. At higher resolution, J0907+3257 reveals the presence of a triple structure of 133 pc with PA roughly perpendicular to the extended disc emission, whereas J0020-0028 shows scattered radio emission along the JVLA radio structure with a north-west marginally elongated core towards the JVLA-detected twin shape.
Considering the radio spectra of the sample, we note a relation between the core dominance of the source (i.e ratio between the pc-scale core measured from our maps and the total radio emission measured from NVSS) and the spectral slope (S ν ∝ ν α ) measured between 1.4 and 4.5 GHz from JVLA data taken from [4,8] (see Fig. 4). The contribution from the core component of the total extended emission increases with the spectral flattening. This is in line with a tendency of the flat/inverted-spectrum sources of being unresolved or one-sided jetted. In addition, it is noteworthy to point out that three sources (J0943+3614, J1025+1022, J1530+2705), which are characterised by an inverted/flat spectrum and by a substantial flux increase of their JVLA detections with respect to their NVSS-FIRST measurements by a factor 1.6-3.2 [4,8], show a (one-sided) compact core, possibly indicating of a marginal alignment of the source with the line of sight.
The pc-scale (peak) core luminosities L core of the eMERLIN and EVN FR 0 sub-samples range between 10 21.3 -10 23.6 W Hz −1 (integrated luminosities νL core ∼ 10 37.5 -10 39.8 erg s −1 ), a factor 10-100 weaker that the FR 0s studied by [31,32] (10 23 -10 24 W Hz −1 ). The pc-scale cores detected by eMERLIN and EVN offer a more robust estimate of the jet base energetics than the sub-arcsec cores detected by JVLA. In this interpretation, Figure 5 represents a diagnostic plot, radio vs [O III] luminosity, i.e. a proxy of the jet energetics vs a proxy of the AGN bolometric luminosities 9 . The two quantities are expected to correlate following the radio-line relationship valid for FR 0s and FR Is [4,8], suggesting that radio emission efficiency, i.e. the fraction of the radio emission produced with respect to the AGN accretion power, is constant between FR 0s and FR Is. As we can note from  red dots for the eMERLIN FR 0 sub-sample, blue filled squares for the EVN FR 0 sub-sample, black stars for FR 0s from [31,32], green empty circle for RL Palomar LINERs, upwards orange triangles for 3C/FR Is, downward pink triangle for Core Galaxies (see Sect. 4 for details). The dotted line indicates the best linear correlation found for all the sources. core and line distributions could be the result of a sharper dependence on possible Doppler flux boosting or larger core variability on parsec scales than those from a sub-arcsec view of the JVLA cores, apart from a possible flattening of the slope of the radio-line correlation. If we include the FR 0s studied with the VLBI technique by [31,32] (black stars in Fig. 5), the data points cluster in limited intervals of radio and [O III] luminosities 10 , preventing from a statistical study of the presence of a L core -L [O III] correlation valid for FR 0s.
Since an affinity between the nuclear properties of FR 0s and FR Is has been clearly assessed [4,8,136] and low-power RGs have been found to be the scaled-down version of FR Is (e.g., [7,11,13,14,102]), it is plausible to search for a general connection between the jet and accretion strength of a population of low and high power RGs with comparable properties: hosted in massive ETGs and characterised by a LINER spectrum 11 . With the purpose of  [24,67]), see Fig 5, for: LINERs from the Palomar sample of nearby galaxies [66] (green circles, from [1,77,83,102]), which appear as unresolved cores or core-jet structures at JVLA scales in the radio band and are classified as RL by [11]; 3C/FR Is (all LINERs, open upwards triangles, from [52,72,83,139]); and Core Galaxies (open downward triangles, from [102]), which have LINER spectra and are known to be miniature FR Is [13]. With the inclusion of these RL AGN in the L core -L [O III] plot, it is evident that FR 0 data points belong to a broad sequence between the two luminosities correlated on ∼4 order of magnitudes. We fit the data points present in this sequence with a linear (in a log-log plot, hereafter) relation and we find a correlation in the form L [O III] ∝ L 0.72±0.10 core with a Pearson correlation coefficient (r-value) of 0.782 which indicates that the two quantities do not correlate with a probability of 3×10 −11 . This statistically-robust relation corroborates the idea that nuclei of FR 0s are similar to those of LINER-like RL AGN. The large scatter of the correlation, ∼0.27 dex, could be caused by the Doppler boosting, the nuclear variability and the different data frequencies (particularly sensitive in the case of steep spectrum). Furthermore, the slope of the established linear correlation is consistent with that found for nearby low-luminosity RL LINERs studied with eMERLIN [11]. However, a study on a large sample of FR 0s at lower radio luminosities is still needed to confirm the existence of this linear regression.
With the purpose of studying the ejection mechanism in FR 0s, for the extended sources of our eMERLIN and EVN sub-samples, we estimate the jet/counter-jet ratio measuring the brightness ratio in two symmetric regions as near as possible to the nuclear emission by avoiding the core and considering the surface brightness at similar distance from the core. The jet-to-counter-jet sidedness ratios measured for the JVLA+eMERLIN and EVN sub-samples (Tab 2 and 3) are between 1.1-1.8 and 1-3, respectively. For each source also resolved by the JVLA, the mas-scale jet ratio is higher than the corresponding value measured with JVLA at sub-arcsec scale (in the range, 1-2 [8]). Figure 6 (upper panel) displays the pc-scale jet sidedness distribution measured from our data and from VLBI observations taken from [31,32] 12 . The distribution of jet-to-counter-jet brightness ratio for our FR 0s and the sample from [31] goes between 1 and 3 with one exceptional value at 6.5 (J0907+3257). The shaded histogram in Fig. 6 (upper panel) represents the lower limits for the one-sided jetted sources.
From this parameter we estimate the bulk Lorentz Γ factor of the jet following the procedure discussed by [16]. Since the jet parameters are unknown, this method assumes the value for angle θ m that maximizes β, i.e. cos θ m = β, where β = v/c is the intrinsic jet speed, in units of c and θ is the angle of the jet to the line of sight. With this simplification, the maximum Γ bulk is where R is the jet sidedness and m = k − α, where k is the parameter that accounts for the geometry of the ejecta, with k = 2 for a continuous jet and α is the spectral index of the emission (S ν ∝ ν α , taken from Tab. 2 and Tab. 3 and from [31,32]). With these assumptions, the range of Γ bulk measured for our sample and the FR 0s from [31,32] is between ∼1 and 2.5 with two values at ∼3.5 and 12 for the targets J0933+1009 and J0907+3257, respectively (Fig.  6, lower panel). We plan to re-observe these last two FR 0s with VLBI array to check their particular extreme jet properties with respect to the rest of the sample.

Discussion
Our eMERLIN and EVN observations of 15 FR 0s confirm the ability of this class of RGs to launch pc-scale jets, at lower luminosities ( 10 23 W Hz −1 ) than previous VLBI studies of FR 0s. Radio structures with extents from a few tens to hundreds of pc are present for ∼86 % of the FR 0s studied so far with VLBI technique (this work, [31,32]), despite 4 FR 0s reveal their jetted structures only after combining JVLA and eMERLIN datasets. This result eventually attests the idea that FR 0s are genuine RGs with extended jets, although of limited size. The resolved jetted structures typically host a steep-spectrum core, while the flat/inverted spectrum FR 0s appear as one-sided jets or individual core components. This morphologyspectrum connection is the result of a core dominance with respect to the total emission, increasing with the flattening/inverting of the radio spectra. This highlights the power of the VLBI-based observations in resolving the pc-scale flat-spectrum core in 'compact' sources like FR 0s.
The effective capability of launching pc-scale of FR 0 jets adds up to a further similarity with classical FR Is, which generally exhibit core-brightened 13 radio morphologies at parsec scale [43,52,138]. At kpc-scale, about half of the FR Is shows large jet asymmetries with a jet/counter-jet flux ratio larger than 2 [23,115]. At smaller scales, the jet-to-counter-jet sidedness ratio of FR Is [49][50][51]138] is generally larger than that measured with the JVLA [140]. 12 We measure the jet components for the one-sided and two-sided jets at L and C band (or X band if the previous two bands are not available) from VLBI, VLBA and EVN data from [31,32]. 13 However, it is correct to point out that the observed pc-scale core-brightening morphology of FR Is could be the effect of a selection bias because only FR Is with bright cores have been observed with the VLBI [83].
This morphological jet variation of FR Is at different scales is interpreted with a change of the jet bulk speed, initially from relativistic, Γ >3, to sub-relativistic speeds on kpc scales by decelerating for entrainment [20,75]. In comparison to FR Is, for the sample of 22 FR 0s studied at parsec scale (our work, [31,32]), the jet-to-counter-jet ratio is less prominent: about one third of FR 0s has jet sidedness larger than 2. With inevitable assumptions on jet properties, the estimated FR 0 jet bulk Γ factor, maximised to reproduce the observed distribution of jet sidedness, ranges between 1-2.5 for ∼90 % of the FR 0 sample. This is in agreement with results from [31,32] that FR 0 jets are mildly relativistic with low jet proper motions. However, a proper systematic analysis on larger samples of FR 0s matching in redshifts and luminosities with FR Is is needed to draw a final conclusion on the jet bulk speed difference between the two populations.
Another outstanding result of this study, which strengthens the connection between FR 0s and FR Is, is the presence of a broad correlation between the pc-scale core luminosities and the [O III] line luminosities, used as proxy of the AGN strength, established for low-power LINER-like RGs, which encompass FR 0s, FR Is, and low-luminosity RL LINERs in general. This result suggests that RL LINERs, other than sharing the same type of host (massive ETGs with BH mass > 10 7.5 M ), also shares a single type of central engine, able to launch (from low to high power) core-dominated jets [7,10,27] with (from mildly to highly) relativistic speeds. A radiatively inefficient accretion disc is the most plausible engine to account for the low Eddington ratios and low bolometric luminosities of FR 0s and FR Is [14,15,136]. This hypothesis is supported by the idea of an ADAF disc coupled with a jet, typically attributed to the RL LINER population in general, in agreement with previous studies (e.g. [11,33,41,108]). In fact, RIAF discs are supposed to be efficient at emanating jets, even at low energies, as suggested by theoretical studies including analytical works (e.g., [91,103,107] and numerical simulations (e.g., [90,135]).
The best scenario, which would simultaneously predict low-brightness pc-scale jets of FR 0s and a disc-jet coupling and host conditions similar to those of more powerful FR Is, but differing in large-scale environment, is a specific cosmological evolution of RGs based on BH spin and clustering. In an accretion-driven scenario, FR 0s host prograde low-spinning BHs, which will eventually evolve into highly spinning FR Is when the accreting material reaches ∼30% of the initial BH mass [46]. The poorer environment of FR 0s, as compared to the one of FR Is, would favor the longer stage of FR 0s due to the limited gas availability and thus would account for the larger space density of compact RGs than the extended ones [30]. Assuming a BH spin-Γ jet dependence (e.g., [86,134]), the lower bulk Lorentz factor of FR 0s would be caused by their lower BH spin, limiting so their extracted energy to launch mildly-relativistic jets. Instead, in a merger-driven scenario, assuming that the BH spin is the main result of a BH-BH coalescence event, in a poorer environment major mergers of equal mass galaxies are less frequent [71,131], resulting in low probabilities of reproducing highly-spinning BHs, a condition required to create extended relativistic jets. Therefore, the large-scale environment generates different merger histories and different BH spin distributions between FR 0s and FR Is, resulting eventually in a large and small population of RGs with compact and fullyfledged morphologies, respectively.
While the majority of FR 0 population conforms with the idea that they are compact core-dominated sources with mildly-relativistic core-brightened jets, a small fraction exhibits different radio features from this picture. In fact, different studies at low and high radio frequencies suggested that FR 0s certainly are a mixed population of low-power radio sources which includes genuine young RGs, such as CSS and GPS [94,111,125]. In our sample, we must mention that one source (J0943+3614) has an inverted spectrum, similar to what is typically seen in GHz-peaked sources [111], but previous VLBI and EVN observations of this target resolve only a core component. In addition, two sources (J0907+3257 and J1230+4700) show a mas-scale double-lobed structures consistent with a classical picture of a new-born RG. The three variable radio sources (J0943+3614, J1025+1022, J1530+2705), instead, could reconcile with the picture of the aligned counterparts of FR 0s, as the BL Lacs population has been claimed as the parental population of FR 0s [88]. All these possible exceptions to the general FR 0 rule contribute to an approximately estimated ∼10% of contamination present in the FR 0 population [6], where young RGs, blazars, radio-quiet AGN and low-luminosity compact sources 14 [4,94,111,125] would erroneously fall in the FR 0 category, if their radio-optical properties are not deeply analysed. Nevertheless, it is possible that a genuine population of FR I progenitors or aligned weak BL Lacs whose properties have a large overlap with those of misaligned FR 0s, cannot be fully disentangled from the bona-fide FR 0 population.

Summary and Conclusions
We present eMERLIN and EVN observations of a sample of 15 FR 0s, respectively, 5 sources at 5 GHz and 10 sources at 1.66 GHz. Parsec scale one-and two-sided jets have been detected for 11 sources of the sample and cores for the remaining objects with angular resolution of a few mas and sensitivity of several tens of µJy.
Since we also own calibrated JVLA C-band data for the five FR 0s observed with eMERLIN, for the first time we combined the visibility datasets of the two arrays in this band for this population of apparently jet-less compact RGs. This method aims at probing the intermediate scales of the jet formation and propagation between the kpc-scale emission lacking in FR 0s and the sub-kpc core resolved with JVLA. This procedure turns out to be successful in detecting pc-scale jets, which were missing in the two original datasets. This result means that the jets appear unresolved and resolved out, respectively, in the JVLA and eMERLIN maps. We can thus conclude that FR 0s, although apparently lacking of extended emission, are effectively able to emanate pc-scale jets, whose low brightness make them hard to isolate and detect. The FR 0s, which do not show extended jets, are typically variable and/or one-sided or having unresolved core, suggesting a possible jet orientation effect. In any case, the combination of long and short baseline represents a powerful tool to study the jet properties of the FR 0 population.
The compact component detected with mas-scale angular resolution, at the centre of the resolved jetted structures of FR 0s studied with the VLBI technique, delimits the parsec-scale section of the jet along its extension. This radio-emitting region of a few pc represents a closer look to the jet launching region and thus a better proxy of the jet energetics than previous radio studies. The core identification of FR 0s enabled us to carry on an unbiased comparison with the rest of widely-studied population of RGs. In fact, it has been found that FR 0s, together with a homogeneous population of low-power RGs, all characterised by a LINER optical spectrum, an early-type host and large BH masses (>10 7.5 M , satisfying the radio-loudness criteria [34]) appear to follow a broad correlation between the pc-scale core luminosities, considered as good proxy of the jet kinetic power, and the [O III] line luminosities, indicative of the accretion power. This result sets a disc-jet connection similar to FR Is, where a RIAF disc supports the launch of a core-brightened radio structure, valid for FR 0s and generally for lowpower LINER-type RGs. In conclusion, a particular set of nuclear and host properties defines a homogeneous class of RGs with a common accretion-ejection mode, from low-luminosity RL AGN to powerful FR Is. However, although they show large differences in radio behavior from the point of view of their sizes, luminosities, and morphologies, LINER-like low-luminosity AGN and FR 0s represent the low-end part of this continuous population distribution, which ends with fully-fledged FR Is. An increasing radio luminosity, possibly due to higher BH spins or other higher values of the internal parameters (BH mass, magnetic field, etc.), would smoothly connect the low-power RGs and FR 0s to the FR Is.
Low and high frequency radio, optical and X-ray observations of FR 0s have shown that FR 0s differ from FR Is for the large space density, the large-scale radio structure and Mpc-scale environment, unlike their nuclear and host properties which appear indistinguishable. These similarities and differences have been interpreted in the framework of RG evolution: FR 0 state consists of a long phase of the life-time of RGs, which inhabit moderately poor environment and are powered by low spinning BHs, which then struggle to spin up for the low accretion rate. This scenario predicts that FR 0s are characterised by low jet bulk speed, assuming a linear connection between the relativistic flow of the jets and the BH spin, from where the jet energy is extracted. This interpretation is supported by our eMERLIN and EVN study: the bulk jet speeds of FR 0s, estimated from the pc-scale sidedness, are mostly in the range 1-2.5, lower than typical values >3 found for FR Is. In a two-zone jet model with an inner fast spine and a slower layer [47], a weak/short (mildly) relativistic spine, which then possibly decelerates and disrupts within the galaxy via gas entrainment [20,80,81], could reconcile with the high-energy and radio properties of FR 0s [9]. Nevertheless, more dedicated studies on a statistically-complete sample of FR 0s with VLBI-technique observations are needed to confirm their slower jets than those of FR Is and to eventually understand the physical mechanism of jet launching/propagation in this still-marginally explored class of compact RGs.
The fast approaching high-resolution and high-sensitivity radio facilities, such as the SKA [73] and the ngVLA [110], offer promising opportunities to place more stringent constraints on the jet speed and unearth large population of FR 0s in the local Universe. The observations with these up-coming observatories will shed light on the mechanisms of accretion-ejection coupling in this abundant population of compact RGs within a general picture of evolutionary and orientation-dependent RL AGN unification scheme.
Author Contributions: R.D. Baldi had analysed the eMERLIN data and combined the eMERLIN and JVLA data. G. Giovannini has analysed the EVN dataset. All the authors contribute to the discussion of the result of this work. R.D. Baldi has written the article and the co-authors gave comments and suggestions to improve the quality of the manuscript.

Appendix A
We present Table A1, which provides the radio contours, rms and the properties of the restoring beams of the eMERLIN, combined JVLA-eMERLIN and EVN maps of our sample.