A transition from parabolic to conical shape as a common effect in nearby AGN jets

Observational studies of collimation in jets in active galactic nuclei (AGN) are a key to understanding their formation and acceleration processes. We have performed an automated search for jet shape transitions in a sample of 367 AGN using VLBA data at 15 GHz and 1.4 GHz. This search has found ten out of 29 nearby jets at redshifts z < 0.07 with a transition from a parabolic to conical shape, while the full analyzed sample is dominated by distant AGN with a typical z ≈ 1. The ten AGN are UGC 00773, NGC 1052, 3C 111, 3C 120, TXS 0815−094, Mrk 180, PKS 1514+00, NGC 6251, 3C 371, and BL Lac. We conclude that the geometry transition may be a common effect in AGN jets. It can be observed only when sufficient linear resolution is obtained. Supplementing these results with previously reported shape breaks in the nearby AGN 1H 0323+342 and M87, we estimate that the break occurs at 10–10 gravitational radii from the nucleus. We suggest that the jet shape transition happens when the bulk plasma kinetic energy flux becomes equal to the Poynting energy flux, while the ambient medium pressure is assumed to be governed by Bondi accretion. In general, the break point may not coincide with the Bondi radius. The observational data supports our model predictions on the jet acceleration and properties of the break point.


INTRODUCTION
Understanding the physical processes that determine the formation, acceleration and collimation of relativistic jets in active galactic nuclei (AGN) continues to be among the most challenging problems of modern astrophysics. There are a wide variety of analytical and numerical models for jet acceleration and its confinement (e.g., Vlahakis & Königl 2003;Beskin & Nokhrina 2006;McKinney 2006;Komissarov et al. 2007;Tchekhovskoy et al. 2011;McKinney et al. 2012;Potter & Cotter 2015) that consider different solutions for jet shapes, such as cylindrical, conical and parabolic. General relativistic magnetohydrodynamic simulations (e.g., McKinney et al. 2012) predict that a jet starting from its apex has E-mail: yyk@asc.rssi.ru a parabolic streamline within the magnetically dominated acceleration zone. At other scales it transitions to a conical geometry associated with equipartition between energy densities of the magnetic field and the radiating particle populations. It has been shown for cold jets that acceleration should not occur in a conical jet. This requires something akin to a parabolic jet shape closer to the jet base to allow differential expansion (Vlahakis & Königl 2004;Komissarov 2012).
In order to investigate these theories it is important to collect observational data on jet profile shapes for a large enough sample of AGN whose properties are well understood. The first observational evidence for a transition from parabolic to conical jet shape was detected in M87 (Asada & Nakamura 2012) at a distance of about 900 mas near the feature HST-1, about 70 pc in projection, corresponding to 10 5 Schwarzschild radii. A few more studies of nearby AGNs to probe their innermost jet regions were performed recently: Mkn 501 (Giroletti et al. 2008), Centaurus A (Müller et al. 2014), Cygnus A (Boccardi et al. 2016;Nakahara et al. 2019), NGC 6251 (Tseng et al. 2016), 1H 0323+342 , 3C 273 (Akiyama et al. 2018), NGC 4261 (Nakahara et al. 2018, 3C 84 (Giovannini et al. 2018), 3C 264 (Boccardi et al. 2019), NGC 1052 (Nakahara et al. 2020). Hovatta et al. (2019) have indirectly addressed this question for the 3C 273 jet close to the central engine on the basis of a model analysis of ALMA rotation measure data. Larger survey studies (Pushkarev et al. 2009) have typically probed regions farther away from the central nucleus, although Algaba et al. (2017) have used apparent parsec-scale jet base parameters closer in.
In a previous work (Pushkarev et al. 2017), we analyzed parsec-scale radio VLBI images of jets in 362 active galaxies from the MOJAVE program (Lister et al. 2018). This sample is dominated by compact radio bright blazars with a jet at a small angle to the line of sight and a typical redshift z ≈ 1. However, some low luminosity nearby radio galaxies were also included. Pushkarev et al. (2017) show that while the majority of resolved jets have a shape close to conical, a significant fraction of the sample has observed deviations. A systematic change in jet width profile has been noted by Hervet et al. (2017), who explain it by using a stratified jet model with a fast spine and slow but relatively powerful outer layer. In this paper, we investigate if this outcome is partly affected by the typical finite angular resolution of VLBI observations. We probe a possible dependence of the jet shape on the distance r from the nucleus. Furthermore, we perform a systematic search for a possible transition from one jet shape to another on the basis of 15 GHz and 1.4 GHz VLBA images.
The observation of jets with a change from parabolic to conical shape may provide an instrument to probe the MHD acceleration mechanism models as well as the ambient medium conditions. The change in jet shape in M87 (Asada & Nakamura 2012) is coincident with the stationary bright feature HST-1, which can be associated with the change in ambient pressure profile and appearance of a recollimation shock due to pressure drop and abrupt expansion. This interpretation is supported by the measurements of external medium pressure by Russell et al. (2015) almost down to the Bondi radius rB = 2GM/c 2 s (sphere of influence), with an observed mass density profile ρ ∝ r −1 (here cs is a sound speed). The recently observed jet shape in 1H 0323+342 ) demonstrates a similar behavior. On the other hand, there are models predicting a jet shape transition for a single power law pressure profile. The analytical model by Lyubarsky (2009) predicts the transition from parabolic to conical form for certain regimes, as well as quasi-oscillations in jet shape in the conical domain. This solution has been applied to the reconstruction of the recollimation shock properties of M87 by Levinson (2017), with a predicted total jet power on the order of 10 43 erg/s. The recent semi-analytical results for the warm jet matching the ambient medium with a total electric current closed inside a jet by Beskin et al. (2017) predicts a change in a jet shape from parabolic to conical for the Bondi pressure profile P ∝ r −2 . In this work we follow the latter model and consider the results for a warm outflow in more detail.
The structure of the paper is the following: section 2 presents our results of a search for the jet profile change from parabolic to conical in a large sample of AGN jets, we suggest a model and interpret our findings in section 3, a discussion is presented in section 4. We summarize our work in section 5. Throughout this paper we will use the term "core" as the apparent origin of AGN jets that commonly appears as the brightest feature in VLBI images of blazars (e.g., Lobanov 1998;Marscher 2008). We adopt a cosmology with Ωm = 0.27, ΩΛ = 0.73 and H0 = 71 km s −1 Mpc −1 (Komatsu et al. 2009).

Automated search of candidates with a change in jet geometry
For the purposes of our study, we made use of data at 15 GHz from the MOJAVE program, the 2 cm VLBA Survey, and the National Radio Astronomy Observatory (NRAO) data archive (Lister et al. 2018) for those sources that have at least five VLBA observing epochs at 15 GHz between 1994 August 31 and 2016 December 26 inclusive. We used the 15 GHz VLBA total intensity MOJAVE stacked epoch images supplemented by single epoch 1.4 GHz VLBA images to derive apparent jet widths, d, as a function of projected distance r from the jet core, and determined jet shapes similar to Pushkarev et al. (2017). In that work we fitted the d-r dependence with a single power law d ∝ r k . The index is expected to be k ≈ 0.5 for a quasi-parabolic shape and 1.0 for a conical jet. We note that even single-epoch observations at 1.4 GHz adequately reproduce source morphology, i.e., effectively fill jet cross-section due to a steep spectrum of synchrotron emission of the outflow, with a typical spectral index −0.7 measured between 2 and 8 GHz (Pushkarev & Kovalev 2012) and −1.0 between 8 and 15 GHz (Hovatta et al. 2014), making the low-frequency observations sensitive enough to probe jet morphology at larger scales. In our analysis we use the jet width measurements made at 1.4 GHz only on large scales, not covered by the 15 GHz data. These scales are typically beyond 10 mas. This allows us to neglect the core shift effect (e.g., Kovalev et al. 2008;Sokolovsky et al. 2011;Plavin et al. 2019a), which is expected to be about 1 mas between 1.4 and 15 GHz. Its value can not be easily derived since it requires simultaneous observations at different frequencies. As a result, the jet widths estimated at 15 GHz smoothly transition to those at 1.4 GHz. We have carried out a similar analysis allowing for a change in the jet shape. Using all available data (15 GHz only or combined data set at 15 GHz and 1.4 GHz) for each source, we performed a double power law fit of the jet width as a function of distance, dividing the jet path length in a logarithmic scale by two parts in proportion of 1:1, 1:2 and 1:3 to search for cases when the fitted k-index at inner scales was 0.5 ± 0.2, while at outer scales it was 1.0 ± 0.2. After such cases were identified automatically, we tuned the fits by setting the distance of the transition region by eye.
We ended up dropping 36 AGN jets from the original samle of 367 objects as having unsatisfactory fits caused by either (i) non-optimal ridge line reconstruction for jets with  Figure 1. Jet profiles with an indication of transition from parabolic to conical shape in ten well resolved nearby active galaxies. The dependence of the jet width on projected distance from the apparent jet base is shown. The cyan and orange dots show measurements at 15 GHz and 1.4 GHz, respectively. The red and black stripes represent Monte Carlo fits for jet regions before and beyond the jet shape transition region, respectively. The projected distance is shown in pc for targets with known redshift and in mas for 0815−094 which has no redshift information. General properties of these AGN are presented in Table 1, parameters of the fits -in Table 2, parameters of the shape transition region -in Table 4. Table 2. Derived best-fit parameters of the two fitted dependencies d = a 1 (r + r 0 ) k 1 and d = a 2 (r + r 1 ) k 2 before and after the jet break, respectively ( Figure 1). We used the VLBA data at 15 GHz only (band 'U') or 15 GHz and 1.4 GHz (band 'UL') between r min and rmax distance from the apparent core. Note that all the values of r are projected on the plane of the sky.

Source
Band  Table 3. Derived best-fit parameters of a single fit dependence d = a(r + r 0 ) k for 319 AGN. Their k values are presented in Figure 2. We used the VLBA data at 15 GHz only (band 'U') or 15 GHz and 1.4 GHz (band 'UL') between r min and rmax distance from the apparent core. Full strong bending, (ii) numerous large gaps in jet emission, (iii) too short a jet length (iv) low intensity regions not captured well by our jet width fitting. This resulted in a sample comprising 331 AGN jets. As a result of this analysis, we found a shape transition in ten jets ( Figure 1, Table 1) out of 367 analyzed. We emphasize that all the AGNs with detected transition of the jet shape turned out to have low redshifts z < 0.07, i.e., have a high linear resolution of 15 GHz VLBA observations -better than 1 pc. This is highly unlikely to occur by chance and provides additional strong evidence that this result is not an observational artifact but a real effect. See discussion of the rest of analyzed low redshift AGN in the sample in subsection 2.4. Among the ten sources, there is one, the radio galaxy 0238−084 (NGC 1052), that shows a two-sided jet morphology. For this object, we analyzed the approaching, brighter outflow propagating to north-east direction, determining the position of a virtual VLBI core using a kinematic-based minimization method described in Vermeulen et al. (2003).
Following our discovery of the shape transition preferentially occurring in nearby AGN, we supplemented our initial AGN sample of 362 targets from Pushkarev et al. (2017) with stacked images of five more low-z AGN which had five or more 15 GHz VLBA observing epochs after the Pushkarev et al. analysis was finished. These were: 0615−172, 1133+704, 1200+608, 1216+061, 1741+196. All the stacked images are available from the MOJAVE database 1 .

Rigorous fitting of the jet shape
For each of the 10 sources found to have a jet geometry transition, we fit the data with the following dependencies: d = a1(r + r0) k 1 and d = a2(r + r1) k 2 , describing a jet shape before and after the break. Here r0 is understood as the separation of the 15 GHz apparent core from the true jet origin due to the synchrotron opacity (e.g., Lobanov 1998;Pushkarev et al. 2012a), while r1 shows how much one underestimates the jet length if it is derived from the data only beyond the geometrical transition of the jet. We note that this approach is more accurate but more computationally intensive than that used by Pushkarev et al. (2017) and applied in the original selection of jet break candidates. It is needed in order to better fit for jet shape close to the apex.
We fit these dependencies with Bayesian modeling using the NUTS Markov Chain Monte Carlo sampler based on the gradient of the log posterior density. It was implemented in PYMC3 (Salvatier et al. 2016), which automatically accounts for uncertainties of all the parameters in further inferences. The best fit parameters are listed in Table 2, showing that initially the jets are quasi-parabolic with k1 close to 0.5, while beyond the break point region the outflow manifests a streamline close to conical, with k2 ≈ 1. The location of the jet shape break given in Table 4 is estimated as the intersection point of these two d − r dependencies. Note that Table 4 includes results on the jet shape transition region for two more sources, 1H 0323+342 and M87 taken from Hada et al. (2018) and Nokhrina et al. (2019), respectively. We also note that the shown error of the deprojected position of the break is propagated from the fitting procedure, it does not include uncertainties on the viewing angle and black hole mass.
For other sources without a detected shape break, we fit a single power-law d = a(r + r0) k for consistency. We excluded objects with unreliable ridge line detection or patchy structure in images (15 sources) and those with nonphysical d − r dependence (24 sources) after visual inspection. They constitute only about one tenth of the dataset and thus the exclusion should not bias our estimates. To account for increased uncertainties of jet width measurements further from the core, the power law model is complemented as following: Here all of a, r0, k, R, σ1, b, σ2 are treated as unknown parameters and inferred simultaneously using a Nested Sampling algorithm as implemented in PolyChord (Handley et al. 2015). As expected, σ2 is typically significantly larger than σ1. We find that this model generally captures the d − r dependence and its uncertainty well. Fitting results are given in Table 3 and the source distribution of exponents k is shown in Figure 2. Even though the estimates for individual sources have a large spread, the median exponent is very close to 1. This indicates a conical average outflow shape, and agrees with previous results using slightly different estimation method (Pushkarev et al. 2012b). We note the peak in the histogram bin at k = 0.5 which corresponds to the parabolic jet shape; the number of objects with k ≈ 0.5 is not high enough in the sample to make it significant (Table 3).

Checking consistency of the fits and analyzing for possible biases
By setting r = 0 we can estimate the apparent core size  Table 5. Two sources, 0111+021 and 0415+379, show a good agreement between d MC c and d uv c , while for the other seven objects d MC c is somewhat larger than d uv c . This is likely due to a non-ideal determination of the core position throughout the epochs, which is used to align single-epoch maps to produce stacked images. The radio galaxy 0430+052 is the only source having d MC c < d uv c for reasons that are unclear.
A bias related to this effect might affect the results. Statistically analyzing jet shapes for the whole sample of 331 sources with stacked VLBA images by introducing different ridge line path length limits we have found the following. A near-parabolic streamline for quasars and BL Lacs can be derived if the innermost jet, only up to ∼1 mas from the apparent core, is considered. This is not a real effect. The bias is found to be the most pronounced for curved jets or jets with features emerging at different position angles over time (Lister et al. 2013). This is confirmed by an apparent artificial correlation of median jet width with the number of epochs in a stacked image for such AGN. Uncertainties in the core position also contribute to this effect due to the imperfect alignment of images while performing the stacking. Variability of opacity conditions and apparent position of the core (Plavin et al. 2019a) affect this partially even though the alignment of the stacked single epoch images is done on the core position. Together, it causes an additional artificial widening of the jet near the core region up to distances r ≈ 0.3 mas. The effect quickly vanishes at larger scales. Thus, if we exclude jet width measurements at distances 0.4 mas, the effect becomes much weaker and disappears completely if we rule out the measurements within 0.5 mas from the core. We also note that radio galaxies, being at Jet width (pc) 2200+420 Figure 3. Dependence of the jet width on projected distance to apparent jet base for 0430+052 and 2200+420 in which 43 GHz data are used instead of 15 GHz. The green and orange dots show measurements which are used in the fits at 43 and 1.4 GHz, respectively. The 15 GHz measurements are not included in the fitting, they are shown by the grey color. The red and black stripes represent Monte Carlo fits for jet regions before and beyond the jet shape transition region, respectively. Parameters of the fits are as follows. For 0430+052: a 1 = 0.193 ± 0.007 pc 1−k 1 , r 0 = 0.009 ± 0.042 pc, k 1 = 0.586 ± 0.047, a 2 = 0.241 ± 0.026 pc 1−k 2 , r 1 = −0.593 ± 0.178 pc, k 2 = 1.116 ± 0.029. For 2200+420: a 1 = 0.187 ± 0.025 pc 1−k 1 , r 0 = 0.043 ± 0.074 pc, k 1 = 0.571 ± 0.097, a 2 = 0.430 ± 0.018 pc 1−k 2 , r 1 = −1.219 ± 0.082 pc, k 2 = 1.126 ± 0.011. Note that the derived k-values agree with 15 GHz results presented in Table 2.
low redshift and thus having apparently wider outflows, are much less subject to this effect. The same is true for the sources with a jet shape break shown in Figure 1, as these are low-redshift objects. Only for BL Lac, as the most remote source among them and also having a bright quasi-stationary component near the core (Cohen et al. 2014), we put a conservative limit of 0.9 mas. For the other sources we used the non-cut intervals listed in Table 2, because dropping measurements at r < 0.5 mas did not significantly change the fit parameters. For the remaining sources, we have dropped all measurement for r < 0.5 mas while analyzing the data (Table 3).
Another possible problem might be related to cases where the jet width is completely unresolved. Indeed, this was found for some AGN targets at some epochs from the visibility model fitting of the core (e.g., Kovalev et al. 2005;Lister et al. 2019). We have addressed this issue by dropping all measurement for r < 0.5 mas. Interestingly, the rest of the measured deconvolved jet width values are always positive. If we assume that this is some sort of a positive bias overestimating the width, it should not depend on r for unresolved jets and will result in k values close to zero. This behavior was not seen in our fitting results.
We have also compared the fitted parameter r0 with the core offset from the jet base estimated from the core shift measured between 15 GHz and 8 GHz (Pushkarev et al. 2012a) assuming an inverse frequency dependence r ∝ ν −1 . These quantities, also listed in Table 5, agree well within the errors in four out of six sources having measured core shifts. The large discrepancy for two sources can be explained by the recently recently established phenomenon of significant core shift variability (Plavin et al. 2019a) or the difference between the true jet shape derived by us and the assumed conical jet shape in (Pushkarev et al. 2012a). We note that this result opens a new way to estimate the distance to the true jet origin which does not require an assumption regarding the jet geometry.
We checked and complemented our analysis using 43 GHz data from the Boston University (BU) AGN group 2 2 https://www.bu.edu/blazars/VLBAproject.html for 0430+052 and 2200+420 (Figure 1), which are present in both the MOJAVE and BU samples. For each of these sources we (i) produced stacked total intensity 43 GHz maps, aligning single epoch-images by the position of the VLBA core derived from structure modelfitting of the visibility data, (ii) determined the reconstructed jet ridge line, and (iii) fitted the transverse jet width as a function of distance from the core ( Figure 3). It resulted in the same k-values before and after the break as in our original analysis within the errors (compare Figure 3 and Table 2). The jet shape transition region is found at core separations comparable to those from the 15 GHz data fits but has shifted slightly. We note however that 7 mm jet width estimates are systematically lower than those found from the 15 GHz data due to the weak high frequency synchrotron emission coming from the jet edges. Robust estimates of jet geometry and particularly of the jet width require high dynamic range images which are better sampled at intermediate radio frequencies.
A good agreement between 15 GHz and 1.4 GHz width measurements increases the robustness of our results.
We warn readers about deriving jet shapes from structure model fitting of single-epoch data (e.g., Hervet et al. 2017), as the jet may appear quasi-parabolic (k < 1) up to a certain (typically short) distance from the core and then change its shape to conical (k ≈ 1). This effect occurs in the sources that show variations in their inner jet position angle. Lister et al. (2013) established this as a common, decadetimescale phenomenon for the most heavily monitored AGNs in the MOJAVE sample. Thus, single-epoch VLBI maps may not reveal the whole jet cross-section, but rather a portion of it, especially in the inner jet regions where images are dynamic range limited. Therefore, the conclusions regarding jet geometry based strictly on a modelfit approach should be treated with caution.

Jet shape transition: a common effect in AGN jets, its consequences and prospects
We have found evidence for geometry transition in many jets for which sufficient linear resolution was achieved. This means that a change in jet shape is a common phenomenon Table 4. Derived parameters of the jet shape break for 10 AGN with addition of 0321+340 adopted from Hada et al. (2018) and 1228+126 from Nokhrina et al. (2019). Columns are as follows: (1) source name (B1950); (2) jet width at the break in mas; (3) same as (2) but in pc; (4) projected distance of the break from the apparent core along the jet in mas; (5) projected distance of the break from the BH along the jet in mas; (6) same as (5) but in pc; (7) deprojected distance of the break from the BH along the jet in pc, the parameter uses the estimated viewing angle; (8) (2012). This refers to the well-known HST-1 feature (Chang et al. 2010) which is located too far downstream to be sampled by typical MOJAVE images. Table 5. Angular size of the VLBA core at 15 GHz, d MC c , and its offset from the true jet origin, r MC 0 , derived from our Monte Carlo modeling of the jet width compared with independent MOJAVE core size measurements in the visibility plane, d uv c (Lister et al. 2019), and the core offset, r cs 0 , estimated from the multi-frequency core shift measurements (Pushkarev et al. 2012a which has significant consequences for many high angular resolution astrophysical and astrometric studies. It is difficult to conclude if the geometry transition with measured properties is specific to only nearby radio galaxies and BL Lacs, or can be extended to the AGN class in general. The radio luminosities of the nearby (z < 0.07) jets are much lower than the rest of the sample and this might affect the geometry and transition zone. We note that Figure 4 presents a consistent picture of the power index dependence on the downstream distance for nearby and distant jets.
In total, indications of the transition from parabolic to conical shape are found in 10 out of 29 nearby (z < 0.07) jets observed as part of the MOJAVE program or by other investigators. VLBA archival data from the latter were processed by the MOJAVE team. The reasons for non-detection of a geometry transition in nearby AGN jets are varied. Some jets, e.g., 0007+106 and 1959+650, have too compact structure to study their shapes. Some others, e.g., 0241+622, 0316+413, 1216+061, show purely parabolic streamlines (Table 3), and their transition regions are expected at larger angular scales than those probed by our observations. E-MERLIN or low-frequency VLBA observations are needed. For example, the nearby radio galaxy 1216+061 (z = 0.0075, scale factor 0.15 pc mas −1 ; not shown in Figure 4) has a parabolic streamline with k = 0.64 ± 0.05 out to 7 mas at 15 GHz, corresponding to a deprojected distance of only ≈ 1 pc. We are studying the remaining 12 low-redshift jets that show no sign of a profile break in a followup approved VLBA program.
The other jets in the sample (Table 3), namely 97 %, do not show a clear significant change in jet geometry. We explain this by (i) a large scale factor of the order of 8 pc mas −1 Figure 4. Best fit k-index values plotted against deprojected distance from the 15 GHz VLBA core ( Table 2, Table 3) for the sources listed in Table 1 with measured redshift and viewing angle. Filled circles show fits at 15 GHz only, while empty circles denote results from analyzing measurements at 15 GHz and 1.4 GHz. Horizontal lines denote the scale over which the k-index was measured for every target. The symbols are placed at the median core distance of the analyzed jet portion. Eleven AGN with detected jet shape transition are shown in blue: 0111+021, 0238−084, 0321+340, 0415+379, 0430+052, 1133+704, 1514+004, 1637+826, 1807+698, 2200+420, and M87. The data for 0321+340 and M87 are taken from Hada et al. (2018) and Nokhrina et al. (2019), respectively.
for a typical source in the sample at a redshift of z ∼ 1 and (ii) a small viewing angle typically about several degrees (Pushkarev et al. 2017). Jet power may also play a role, since the MOJAVE sample is flux-density limited and the AGN with z > 0.1 typically have jet luminosities ∼ 2 orders of magnitude higher than the lower-redshift ones. The jets with a detected shape change have an average scaling factor of 0.7 pc mas −1 and, on average, larger viewing angle since 6 out of 12 are radio galaxies. Thus, if a transition region is located at a distance of a few tens of pc, it corresponds to a projected angular separation of < ∼ 1 mas from the apparent jet base at 15 GHz, which is comparable to the typical interferometric restoring beam size. VLBI observations at higher frequencies may be more effective in registering the jet shape transition, since they provide a better angular resolution and are less subject to opacity effects. This would probe scales closer to the jet apex and possible dependencies between acceleration zone extension and the maximum bulk Lorentz factor or jet power, as predicted by Potter & Cotter (2015). On the other hand, the steep spectrum of the optically thin jet emission hinders the tracing of the jet for long distances. The small viewing angles of the bright AGN jets set another limit on any jet shape investigation in the innermost parts. The streamline of an outflow can be studied down to distances at which the jet half-opening angle is still smaller than viewing angle. As shown by Pushkarev et al. (2017), the intrinsic jet opening angle reaches values of a few degrees at scales of the order of 10 pc. This suggests that the jet shape transition phenomenon might be more effectively studied for nearby AGNs that are oriented at larger angles to the line of sight. After considering all the points discussed above, we have begun a dedicated VLBA program in 2019 to search for geometry transitions in 61 AGN jets with z < 0.07 from observations at 15 GHz and 1.4 GHz.
It is a challenging problem to estimate the consequences of this result on astrometry and astrophysics of AGN. VLBI astrometry delivers the position of the true jet apex only if the opacity driven core shift is proportional to the frequency as r ∝ ν −1 (Porcas 2009). However, this is expected only for conical jets and synchrotron opacity (Lobanov 1998). A non-conical jet base results in an extension of the true jet length between the apex and the observed opaque core. This also produces somewhat larger VLBI-Gaia offsets for AGN positions Plavin et al. 2019b) than predicted by Kovalev et al. (2008).

Deprojected position of the jet break
We chose the MOJAVE-1 sample of 135 AGN ) to perform a direct comparison with the 12 jets showing the breaks. Our reasoning is as follows. Most of MOJAVE-1 targets were observed by VLBA not only at multiple 15 GHz epochs but also in a single epoch at 1.4 GHz, which increases the jet distance probed by our analysis. In addition, VLBI measurements of the apparent kinematics βapp (Lister et al. 2019) and variability Doppler factor estimates δ (Hovatta et al. 2009;Liodakis et al. 2017) are available for a large fraction of the sample. We need this information to derive deprojected distance values. These requirements result in a sample of 65 sources (Table 1) described in Pushkarev et al. (2017).
We derived viewing angle estimates through the relation θ = arctan 2βapp β 2 app + δ 2 var − 1 to convert the jet distance from angular projected to linear deprojected. Note that this assumes the same beaming parameters for the flux density variability and jet kinematics. For βapp we used the fastest non-accelerating apparent jet speeds from the MOJAVE kinematic analysis. For 1H 0323+342 we use θ = 6.3 • , based on the observed superluminal motion (Lister et al. 2016) assuming θ = (1 + βapp) −0.5 = γ −1 , which minimizes the required bulk Lorentz factor γ. The other possible viewing angle value for this target θ = 4 • is based on the variability time scale . For the BL Lac objects 0111+021 and 1133+704 we assumed a viewing angle of 5 • , typical for this class of AGN (Hovatta et al. 2009;Savolainen et al. 2010;Pushkarev et al. 2017;Liodakis et al. 2017). For the radio galaxy 1514+004 we assumed a viewing angle of 15 • which is typical for this class of AGN in our sample.
In Figure 4, we plot the corresponding single powerlaw k-index values derived from the 15 GHz and 1.4 GHz VLBA data (Pushkarev et al. 2017) versus deprojected distance from the 15 GHz VLBA core for 62 sources. There are eleven sources with known deprojected linear jet distance that have a jet shape transition (Figure 1, Table 1). They are shown by a pair of points each from the double power-law fits. The BL Lac object 0815−094 is not shown in Fig. 4, as it does not have a measured spectroscopic redshift. Our results on jet shape transition (Table 2, Table 4, Figure 4) are supplemented by multi-frequency data for M87 from Nakamura et al. (2018), with k1 = 0.57, k2 = 0.90, and break point position obtained by Nokhrina et al. (2019). For M87 we adopt  Table 1. When available, we use the mass estimates based on the velocity dispersion method, otherwise -those from reverberation technique. Note, the rightmost source with the detected transition from parabolic to conical shape is 1H 0323+342. Its mass estimate is based on reverberation mapping and might be strongly underestimated as argued by León Tavares et al. (2014) and Hada et al. (2018). θ = 14 • (Wang & Zhou 2009), consistent with more recent results by Mertens et al. (2016). For NLSy1 1H 0323+342 we use 1.4-2.3 GHz measurements from VLBA observations , with k1 = 0.6 and k2 = 1.41, for which the jet shape break point position is estimated.
Horizontal lines represent the scales at which k-indices were derived, starting from several tens of mas distance from the 15 GHz VLBA core (see subsection 2.3) and up to distances limited by the sensitivity of our observations. The nearby jets, for which we are probing closer to the central engine, have low k values and show a transition from quasiparabolic values at small scales to quasi-conical at larger scales ( Figure 4). It is possible that at scales greater than ∼ 100 kpc, where jets become diffuse and disruptive, their geometry further changes from conical to hyperbolic, characterized by more rapid expansion (Owen et al. 2000).
In order to plot the observed k-index values as a function of the deprojected distance along jets in gravitational radius rg = GM/c 2 units, we use the black hole masses estimated assuming virialized broad lines region (BLR) motion and correlation between BLR size and UV/optical luminosity (Torrealba et al. 2012;McLure & Jarvis 2002;Vestergaard & Peterson 2006;Landt et al. 2017;Palma et al. 2011;Shaw et al. 2012;Liu et al. 2006). We also use mass values inferred by stellar or gas kinematics methods (e.g., Woo & Urry 2002) for the closest sources. The mass values and references can be found in Table 1. We plot the data in Figure 5. It turns out that the sources with BH masses obtained by stellar velocity dispersion method or stellar/gas kinematics measurements are the subset of the sources with the detected jet shape break (i.e., the closest ones).
Since estimating the black hole mass is a complicated and strongly model-dependent method, some of the values might be significantly in error. By dropping the highest and lowest values as possible outliers of the derived jet break position r break measured in rg we are able to bound its values in the narrower range r break ∈ (10 5 , 10 6 )rg. This is an important result, especially when taken together with our finding that the jet shape transition may be a common phenomenon in nearby or even most of the AGN.
We note the following. The black hole mass of 1H 0323+342 is suspected to be underestimated (León Tavares et al. 2014;Hada et al. 2018). If we use for this source the mass M = 10 8.6 M , obtained using the relation between black hole mass and bulge luminosity (León Tavares et al. 2014), 1H 0323+342 yields r break = 5.6 × 10 6 rg, falling much closer to the discussed above range of r break /rg distances. This may provide an additional argument favoring a higher black hole mass for this source.
We have compared our results for the radio galaxy NGC 6251 with those obtained earlier for this source by Tseng et al. (2016). We have found that the jet shape transition region in this source is at (1.6 ± 0.2) × 10 5 rg, assuming viewing angle of 18 • and black hole mass of 6 × 10 8 M (see Table 1, Table 4). This is slightly smaller compared to (1 − 2) × 10 5 Schwarzschild radius estimated by Tseng et al. (2016), who assumed the same black hole mass and a viewing angle of 19 • . The small difference might be caused by different techniques used to derive it. First, we measured transverse jet widths from the stacked image of the source, using 14 epochs at 15 GHz from the MOJAVE program and archival VLBA data. Second, we have taken into account the synchrotron opacity of the jet base by introducing the parameter r0 that reflects an offset of the apparent 15 GHz core from the true jet apex.
Of 12 sources with observed change in a jet boundary shape 6 are FR I type, 2 are FR II type, and 4 have uncertain classification based on published radio images. This may mean that different environments expected in these two different types of sources on large scales are either the same on the smaller scales, or affect the jet shape in the same way up to 10 6 rg.

Qualitative consideration
Both analytical (see below) and phenomenological (Potter & Cotter 2013 considerations as well as numerical simulations (Komissarov et al. 2009;Tchekhovskoy et al. 2009;Porth et al. 2011) show that for moderate initial magnetization of a jet σM ∼ 10-10 2 , where is the Michel magnetization parameter, the flow transits from a magnetically dominated regime at small distances r from the origin to a particle dominated regime at larger distances. Here Ψ0 and Ω0 are the total magnetic flux and characteristic angular velocity of the "central engine" respectively. Accordingly, µ = mpc 2 + mpw is the relativistic enthalpy, where w is the non relativistic enthalpy, and mp is a particle mass. Here we assume a leptonic jet, so mp is the electron mass. Below for simplicity we consider not so large temperatures, so that w c 2 . Finally, η is the particle-tomagnetic flux ratio.
Indeed, the physical meaning of the Michel magnetization parameter is the maximum Lorentz factor γ of the hydrodynamical flow when all the electromagnetic energy flux is transferred to particles. On the other hand, for quasicylindrical jets the following asymptotic solution for magnetically dominated flow exists (see e.g., Beskin 2009) where RL = c/Ω0 is the light cylinder radius, and r ⊥ is the distance from the jet axis. For the black hole spin a * = 0.5, RL ≈ 14.9 rg ≈ 2.2 × 10 15 (MBH/10 9 M ) cm ≈ 7.1 × 10 −4 (MBH/10 9 M ) pc. Here and below we use the maximum BH energy extraction rate condition ΩF = ΩH/2 (Blandford & Znajek 1977). For observed pc scale jets, the jet width d at the jet shape break point reaches 1 pc. This means that at the transition point d/2RL > σM, and the flow cannot be still magnetically dominated. As was shown by Nokhrina et al. (2015) who have analysed about 100 AGN jets, σM ∼ 10 − 50 is a reasonable value constrained by the observations. The observed median value of 1.02 for the kindex also clearly points to a ballistic plasma motion. This suggests that the jet is dominated by the plasma bulk motion kinetic energy at the deprojected distance longer than ∼ 100 pc or ∼ 10 7 rg rather than by the Poynting flux, as expected close to the launching region. For this reason we aim to explain the break in the d(r) dependence as a consequence of a transition from the magnetically dominated to the particle dominated regime. Below we present the main results of our semi-analytical consideration. Our goal is in evaluating the dependence of the jet width d on an ambient pressure profile Pext(r). The results for the cold jet are presented in Beskin et al. (2017), while here we consider the semi-analytical results for a warm outflow.

Semi-analytical model
Basic equations describing the internal structure of relativistic and non relativistic jets within the Grad-Shafranov approach are now well-established (Heyvaerts & Norman 1989;Pelletier & Pudritz 1992;Lery et al. 1998;Beskin & Malyshkin 2000;Beskin 2009;Lyubarsky 2009). This approach allows us to formulate the problem of finding a stationary axisymmetric magnetohydrodynamic outflow structure (a jet solution) using a set of two differential equations on a magnetic flux function Ψ and an Alfvénic Mach number M. These equations are Bernoulli equation and Grad-Shafranov equation of a force balance perpendicular to magnetic surfaces. The approach allows us to determine the internal structure of axisymmetric stationary jets knowing in general case five "integrals of motion", i.e., energy E(Ψ) and angular momentum L(Ψ) flux, electric potential which connects with angular velocity ΩF(Ψ), entropy s(Ψ), and the particle-to-magnetic flux ratio η(Ψ). All these values are to be constant along magnetic surfaces Ψ = const. Once the Grad-Shafranov and Bernoulli equations are solved for the given integrals, all the other flow properties, such as particle number density, four-velocity, electric current, and Lorentz factor, can be determined from algebraic equations (e.g., Beskin 2009). In particular, it was shown that a jet with total zero electric current can exist only in the presence of an ex-ternal medium with non-negligible pressure Pext. Thus, it is the ambient pressure Pext that is expected to determine the transverse dimension of astrophysical jets. In general, it is a complicated problem to solve the set of Bernoulli and Grad-Shafranov equations. An additional complication is connected with the change of a system type from elliptical to hyperbolic. So, to tackle the problem different simplifications are introduced. Here we simplify the problem, assuming the flow is highly collimated and can be described within the cylindrical geometry, in which case it can be solved numerically (Beskin & Malyshkin 2000).
On the other hand, careful matching of a solution inside the jet with the external medium has not been achieved up to now. The difficulty arises with having a very low energy density of the external medium in comparison with the energy density inside the relativistic jet. For this reason, in most cases an infinitely thin current sheet was introduced. Moreover, an ambient pressure was often modelled by homogeneous magnetic field B 2 ext /8π = Pext. Below we use the approach developed by Beskin et al. (2017). This paper is later referred to as B17. We propose a flow with an electric current closing fully inside a jet. This is achieved by a natural assumption that the integrals L and ΩF vanish at the jet boundary. The second assumption of the model is a vanishing flow velocity at the jet boundary, which leads to vanishing of a poloidal magnetic field component along with a toroidal due to current closure. As a consequence, only a thermal pressure, defined by a sound velocity cjet and particle number density njet, is left at the jet boundary to balance the ambient medium pressure without a current sheet. We solve Grad-Shafranov and Bernoulli equations for the flux function Ψ(r ⊥ ) and the square of an Alfvénic Mach number M 2 (r ⊥ ). The local non-relativistic enthalpy w for a polytropic equation of state with politropic index Γ = 5/3 can be written as where the local particle number density n is obtained from the equation (4) We solve the system of MHD equations (B17) for the boundary conditions Ψ(0) = 0 and We should note that due to vanishing of the integrals L(Ψ) and ΩF(Ψ) at the jet boundary, the thickness of the final current closure domain tends to zero and in B17 it was is not resolved. However, as it was shown, that the total pressure in this region is strictly conserved: This means that the solution we obtain up to the boundary does contain the residual current and, thus, the toroidal magnetic field Bϕ.
The main difference between the result presented here and the result by B17 is in more accurate account for the thermal terms, which can be seen in Equation 4. To obtain the solution, we employ the following iterative procedure. For each fixed fast magnetosonic Mach number at the axis M 2 0 we initially set Pext at the jet boundary. It defines the particle number density at the boundary njet, and together with M 2 0 -the particle number density at the axis n0. Having set the latter, we solve MHD equations across a jet from the axis outwards and calculate the jet pressure at the boundary provided by the solution P (solution) . By iterations we find self-consistently such Pext that is equal to one, provided by the solution: P (solution) = Pext. Thus, we obtain the dependence of a jet pressure at the boundary as a function of a local jet width d.
This procedure fully determines the solution of our problem. For each magnitude of the external pressure the obtained solution is a crosscut at r = const. Piling of these different crosscuts is a solution for an outflow in which one may neglect by the derivatives over r in comparison with the derivatives over r ⊥ in the two-dimensional Grad-Shafranov and Bernoulli equations. This can be done for highly collimated, at least as a parabola, outflows (Nokhrina et al. 2015) and flows with small opening angles (Tchekhovskoy et al. 2009).
We find that for the chosen sound velocity at the boundary c 2 0 = 0.001c 2 the thermal effects may be neglected in the outflow volume, playing an important role only at the outflow boundary. It turns out that the resultant dependence of pressure at the jet boundary as a function of jet radius obtained by B17 and here start to differ somewhat only for large M 2 0 (this value is of an order of 10, but depends on the initial magnetization), affecting the flow boundary shape downstream of the equipartition transition, and the effect on k2-index is of the order of a few per cent. We will address the particular effects of higher temperature in the future work. The proposed jet model with an electric current enclosed inside the jet has a natural sheath structure, observed, for example, in the M87 jet . Due to choice of integrals, the outer parts of a jet have slower velocities, tending to non-relativistic with γ(d/2) = 1. Such a sheath may be produced by different mechanisms: it may be a slower disk wind or an outer jet disturbed and slowed down by the pinch instability (Chatterjee et al. 2019). In our model it appears naturally as a consequence of a jet transiting into the ambient medium with the hydrodynamical discontinuity only (B17).

Transition from magnetically dominated to particle dominated flow
It is necessary to stress that this system of equations can describe both magnetically and particle dominated flow, with the physical answer (including the jet boundary radius d/2) depending on one external parameter only, namely, on the ambient pressure Pext. In Figure 6 we show the dependence of the dimensionless ambient pressurep on a dimensionless jet widthd obtained by solving numerically the system of Grad-Shafranov and Bernoulli equations B17. The pressure is plotted in units of so that Pext =p p0, and the jet width in units of light cylinder radius is d =d RL. We observe (see Figure 6) that the pressure has a different power law dependence on the jet radius for small and large d. For each magnetization σM, this behavior holds, with the change between two profiles occurring at different jet widths. For σM = 50 the pressure changes its dependence on d from closer to the jet base to further downstream. The particular exponents of the power laws depend weakly on σM.
We assume the equilibrium between jet and ambient medium pressure. In order to model a jet shape break position along the jet, we need to introduce the exerted pressure dependence on r, which we choose in the power law form Such a pressure profile is consistent with Bondi flow (Quataert & Narayan 2000;Shcherbakov 2008; Narayan & Fabian 2011) having b ∈ (1.5; 2.5) for different models, with the limiting value 2.5 for classical supersonic Bondi flow. This power law with b ≈ 2.0 allows us to reproduce well both the parabolic jet form upstream the break and conical downstream. Using power laws Equation 8, Equation 9, and Equation 10, we obtain for small distances r (magnetically dominated regime) Accordingly, for large distances (saturation regime) d ∝ r 0.83 .
As we see, qualitatively, the power indices are in good agreement with the observational data. Thus, we are able to reproduce the jet boundary shape behaviour without introducing two different pressure profiles, as was done in (Asada & Nakamura 2012). Having the reasonable pressure dependence on a distance, we reproduce both power laws in a jet shape. For example, for a central mass M = 10 9 M and black hole spin a * = 0.5 the light cylinder radius is RL ≈ 7 × 10 −4 pc. We also set the total magnetic flux in an outflow Ψ0 = 10 32 G cm 2 (Zamaninasab et al. 2014;Zdziarski et al. 2015;Nokhrina 2017), which gives the value B(rg) ≈ 1400 G. Thus, for these test parameters the jet width at the break, designated by a star in Figure 6, has typical values 0.2 − 1.0 pc in agreement with the observational results in Table 4. In dimensionless units the point of transition from one power law for pressure as a function of a jet width to the other is defined by one parameter only: the jet initial magnetization. In the equipartition regime the jet bulk Lorentz factor is γ = σM/2. The observed kinematics in parsec-scale jets constrains the initial magnetization to a value 100 (Lister et al. 2016), while estimates for σM based on coreshift effect measurements provide the preferred value 20 (Nokhrina et al. 2015). In dimensional units the jet width at the break depends also on BH mass and spin. The distance to a shape transition along the jet is determined by the total magnetic flux in a jet and the ambient medium pressure. We address the question of bounding these parameters in the next paper (in preparation).

Magnetization
In this subsection we check whether the break in a jet shape corresponds to the transition from the magneticallydominated into the equipartition regime. The jet magnetization is defined as the ratio of Poynting flux to particle kinetic energy flux where n is particle number density in the jet proper frame. Using the standard expressions for ideal MHD velocities and electric and magnetic fields, one obtains the following expression for the magnetization: Using the definitions of bulk Lorentz factor γ and total current I, we rewrite it as In order to check σ along the jet, we calculate the maximal magnetization across the jet for each given distance r. The magnetization is always much less than the unity at the jet axis and at the jet boundary. The first holds everywhere, since the Poynting flux behaves at the jet axis as if the current density j has no singular behavior at r ⊥ = 0. Thus, σ → 0 at the axis. The same holds for the boundary in a case of the full electric current closure. Due to specific choice of integrals E(Ψ), L(Ψ), and ΩF(Ψ) (B17), the Poynting flux together with the magnetization reach their maximum values at Ψ = Ψ0/2. It is at this magnetic field  Figure 7. An example of a jet boundary shape (blue solid line) for σ M = 50 and P 0 = 10 −6 dyn/cm 2 at r 0 = 10 pc. The jet magnetization at a given distance from its base is plotted by a red solid line, with black vertical line marking σmax = 1. The transition from one power law to the other (green dashed lines) for the jet boundary roughly coincides with the point where the outflow transits from the magnetically dominated to particle dominated (equipartition) regime.
line the flow attains its highest Lorentz factor across the jet for the given distance from the central source. Thus, we choose the maximal magnetization reaching approximately unity as a criteria of a flow attaining the ideal MHD equipartition regime. In Figure 7 we present the maximal magnetization and the break in a jet form. We plot the modelled jet boundary shape for σM = 50, BH and jet parameters the same as in subsection 3.3. The position of a jet shape break along a jet depends on an ambient pressure profile (Equation 7 and Equation 10), and we use here, as an example, P0 = 10 −6 dyn/cm 2 at r0 = 10 pc. We see that the break in jet shape occurs roughly at the distance from the BH, where the flow magnetization becomes equal to unity. For the higher initial magnetization it takes the larger transverse jet dimension in RL to accelerate the flow up to equipartition, according to Equation 2.

Role of a Bondi sphere
In this paper we propose that the jet form change, observed in a dozen of nearby sources, may be explained by an internal flow transition from magnetically dominated to particle dominated regime with the smooth external pressure profile P ∝ r −2 . There are indications, however, that the ambient pressure may have different profiles at different scales. The measurements of particle number density in ISM by Russell et al. (2015) suggest ρ ∝ r −1 from about 400 pc down to expected Bondi radius rB ∼ 100 − 250 pc. The temperature profile on scales 100 − 1000 pc is roughly constant. This means that just outside, or even inside, the Bondi radius, pressure profile is P ∝ r −1 , with no information on it inside a sphere ∼ 150 pc. The position of a sphere of influence is expected to be at a distance 10 5 − 10 6 rg (Blandford et al. 2019). The position of a transition point r break from magnetically dominated to particle dominated regimes predicted by our model for reasonable parameters lay in general in the r break r B shock ACZ P r 2 P r 1 d r 0. 5 d r 0 .8 Figure 8. A schematic jet boundary shape for an ambient pressure with different profiles, changing at the Bondi radius r B . The jet accelerates while sustaining its boundary as a parabola (acceleration and collimation zone, ACZ). After reaching σ = 1 at r break the jet form becomes almost conical up to the Bondi radius.
same interval or inside rB. For example, in the case of M87 we observe r break ≈ 40 pc (Nokhrina et al. 2019) smaller than rB. The same phenomenon has been noted by Nakahara et al. (2018) for NGC 4261, where the structural transition lies well inside the expected sphere of influence. In Figure 8 we present a cartoon for a jet shape with different ambient pressure profile. Inside the Bondi sphere the jet is accelerating effectively up to the distance r break , with predicted parabolic boundary shape described by Equation 11. This is the acceleration and collimation zone (ACZ) discussed by Blandford et al. (2019). For r break < r < rB the jet assumes a close to conical form Equation 12. Up to rB the jet stays in equilibrium with the ambient pressure Pjet = Pext. If for r > rB the ambient pressure has a more shallow profile, the conical particle-dominated jet may become overpressured with a possible appearance of a standing shock. Thus, we predict the presence of a standing bright feature, associated with a shock, outside the Bondi sphere and downstream the break in jet shape. At this shock we may expect plasma heating, with the flow continuing a conical expansion (Blandford et al. 2019). The position of HST-1 in M87 jet in a close vicinity of expected rB and downstream the r break supports this picture.

Additional observational evidence of the break point and predicted evolution of plasma acceleration
For each of the 10 sources with a jet geometry transition detected (Table 4), we checked for slow pattern (βapp < 0.2c) jet features in Lister et al. (2019). We examined if their median locations with respect to the core are positionally associated, i.e., they match within the errors with the position of the derived jet shape break. We found that five sources have a quasi-stationary feature in the region where jet changes its shape, as expected (see discussion in subsection 4.1). This is a factor of 1.5 larger compared to a ratio from the overall statistics of jet kinematics analysis performed at 15 GHz, which reveals a fraction of quasi- stationary jet features to be about 30% (Lister et al. 2019), applying the criterion βapp < 0.2c. We underline that the MOJAVE kinematic analysis uses conservative criteria in cross-identifying components between epochs and selecting robust ones (Lister et al. 2019). This means that the 50% fraction of sources which show a standing feature in the break point region should be considered as a lower limit. This analysis is also conservative because of the requirement of the feature to be coincident with the detected break point. As discussed above, the shock may be located downstream the jet in the vicinity of rB, which position is usually not known. We note that two sources included from other studies, 1H 0323+342 and M87, have the jet shape transition at distances larger than maximum angular scales probed by the MOJAVE 15 GHz observations. We plot in Figure 9 the maximum Lorentz factor of a bulk plasma motion along a jet, which we obtain within our semi-analytical model. The predicted pattern of a bulk Lorentz factor acceleration in magnetically dominated domain is γ ∝ r ⊥ , which provides for a parabolic jet γ ∝ r 0.5 . After the flow reaches equipartition, the acceleration continues slower than any power-law (logarithmically slow) (e.g., Beskin & Nokhrina 2006). There is also a transitional zone between the two regimes. Thus, we would expect for the sources with the detected jet shape break and superluminal motion the following kinematics pattern: efficient Lorentz factor growth before the break point, and cessation of it in the conical region. This expected Lorentz factor behaviour was reported by Hada et al. (2018). The observed in radio band velocity map in M87  shows the acceleration saturation much earlier than the jet shape break. However, observations in the optical band (Biretta et al. 1999) support the acceleration of plasma continuing further, with reported γ = 6 at HST-1, situated downstream the jet shape break. This may point to non detection of fast components in radio.
This prediction is consistent with observations by the MOJAVE program that acceleration is a common property of jet features (e.g., Homan et al. 2015;Lister et al. 2019), reflecting a tendency for increasing Lorentz factors near the base of the jet, with decreasing or constant speeds being more common at projected distances 10 − 20 parsecs (Homan et al. 2015). While decreasing speeds are not a prediction of this model for a change in jet shape, they could naturally occur if the reduction in positive acceleration is also accompanied by entrainment of external material into the jet. Pushkarev et al. (2017) studied AGN jet shapes by measuring the power low index k assuming a d ∝ r k dependence of the observed deconvolved jet width d on the apparent distance from its core r. Most of the jets exhibited k values in the range from 0.5 to 1.5. As it was clearly demonstrated by Pushkarev et al. (2017), high-quality, high-dynamic-range stacked images are needed for an analysis of this kind in order to trace the full jet channel. In view of a few recent exciting reports on jet shape transitions from parabolic to conical (e.g., Asada & Nakamura 2012;Giroletti et al. 2008;Tseng et al. 2016;Hervet et al. 2017;Hada et al. 2018;Akiyama et al. 2018;Nakahara et al. 2018Nakahara et al. , 2019, we have performed a systematic search of such transition using MOJAVE 15 GHz stacked images, supplementing some of them with available single epoch 1.4 GHz VLBA images to trace larger scales.

SUMMARY
Using an automated analysis approach, we have found 10 jets with such transition out of 367 analyzed: 0111+021, 0238−084, 0415+379, 0430+052, 0815−094, 1133+704, 1514+004, 1637+826, 1807+698, 2200+420. Their redshifts lie in the range z < 0.07 except for 0815−094, whose redshift is unknown. For the full analyzed sample the redshift values cover the range from 0.004 to 3.6 with the typical value being about 1. This low-z coincidence is unlikely to have occurred by chance. Taken together with an analysis of possible biases, we conclude that a genuine effect is present in the data for which VLBA reaches the linear resolution better than 1 pc. We would also predict that the BL Lac object 0815−094 is a nearby AGN.
This finding leads to the following important conclusion. A transition from parabolic to conical shape may be a general property of AGN jets. At the same time, we note that AGN observed at higher redshifts typically have higher luminosities and kinetic power, which can affect the collimation properties. This conclusion has important implications for jet models, astrophysics and astrometry of AGN. Measuring this phenomenon requires a search within nearby AGN which is the subject of our current followup study, or increasing the resolution by using Space VLBI (e.g., Giovannini et al. 2018) or high dynamic range high frequency VLBI imaging.
The deprojected distance r break from the nucleus to the break zone is found to be typically 10 pc. Even more interesting due to its relation to jet formation and acceleration models is this value measured in gravitational radius units. We find the range to be r break ∈ (10 5 , 10 6 )rg which corresponds to the typical Bondi radius.
We have developed the following model to explain the observed jet shape break. The accurate matching of a jet outflow with an ambient medium B17 predicts a change in jet shape from parabolic to conical if the ambient medium pressure is assumed to be governed by Bondi accretion. Within the model, a smaller external pressure is needed to support a jet than in earlier models. The transition of predicted jet shape from parabolic to conical occurs in the domain where the bulk plasma kinetic energy flux becomes equal to the Poynting energy flux, i.e., where the bulk flow acceleration reaches saturation (Beskin & Nokhrina 2006). From studying the break properties we can estimate black hole spin and/or mass, jet total magnetic flux, and ambient medium properties as discussed by Nokhrina et al. (in prep.).
The following two model predictions are supported observationally. The break point, where jets start to be plasma dominated energetically, might be a preferable domain for shocks. We detect standing jet features in this region from MOJAVE analysis (Lister et al. 2019) in at least a half of the AGN targets. The plasma acceleration is predicted to decrease significantly at the transition region, which is consistent with MOJAVE acceleration results (Homan et al. 2015;Lister et al. 2019).
Our finding also implies the following (see also discussion in Algaba et al. 2017). The well-known effect of the apparent shift of the core position with frequency due to synchrotron self-absorption does not follow the rcore ∝ ν −1 law all the way up to the true jet base, since a −1 power low index is expected only for a conical jet (Blandford & Königl 1979;Lobanov 1998). Geometrical and physical estimates made on the basis of core shift measurements will need to take this into account while VLBI and VLBI-Gaia astrometry applications will need to correct for it (Porcas 2009) in cases where very high accuracy is required. Table 1: Properties for 12 sources with a detected jet shape break from this study (Figure 1) as well as Hada et al. (2018, 1H 0323+342) and Asada & Nakamura (2012, M 87). They are supplemented by the MOJAVE-1 sources for which redshift values, Doppler factor estimates, and robust jet shape fits (Table 3) are available. Columns are as follows: (1) B1950 name; (2) alias; (3) optical class, where Q = quasar, B = BL Lac, G = radio galaxy, N = Narrow Line Seyfert 1 (NLSy1); (4) redshift; (5) literature reference for the data in column (4); (6) maximum apparent radial speed from Lister et al. (2019), (7) variability Doppler factor from Hovatta et al. (2009); (8) viewing angle; (9) black hole mass estimated basing on assumption of virialized broad lines region (BLR) movement and correlation between the size of BLR and UV/optical luminosity; (10) literature reference for the data in column (9); (11) black hole mass estimated by a stellar velocity dispersion method and associated fundamental plane method (for 2200+420); (12) literature reference for the data in column (11). Names of the sources with the shape break are highlighted by the boldface font.

Source
Alias  Liodakis et al. (2017). b Assumed θ value as typical for BL Lacs. c Assumed θ value as typical for radio galaxies in the list which do not show a strong counter-jet.