Properties of the jet in M87 revealed by its helical structure imaged with the VLBA at 8 and 15 GHz Open Access

ABSTRACT

We present full-track high-resolution radio observations of the jet of the galaxy M87 at 8 and 15 GHz. These observations were taken over three consecutive days in 2009 May using the Very Long Baseline Array (VLBA), one antenna of the Very Large Array (VLA), and the Effelsberg 100 m telescope. Our produced images have dynamic ranges exceeding 20 000:1 and resolve linear scales down to approximately 100 Schwarzschild radii, revealing a limb-brightened jet and a faint, steep spectrum counter-jet. We performed jet-to-counter-jet analysis, which helped estimate the physical parameters of the flow. The rich internal structure of the jet is dominated by three helical threads, likely produced by the Kelvin–Helmholtz (KH) instability developing in a supersonic flow with a Mach number of approximately 20 and an enthalpy ratio of around 0.3. We produce a clean imaging bias-corrected 8–15 GHz spectral index image, which shows spectrum flattening in regions of helical thread intersections. This further supports the KH origin of the observed internal structure of the jet. We detect polarized emission in the jet at distances of approximately 20 milliarcseconds from the core and find Faraday rotation which follows a transverse gradient across the jet. We apply Faraday rotation correction to the polarization position angle and find that the position angle changes as a function of distance from the jet axis, which suggests the presence of a helical magnetic field.

1 INTRODUCTION

M87 (Virgo A, NGC 4486, 3C 274) is a giant elliptical galaxy which has an Active Galactic Nucleus (AGN) and an extended relativistic jet powered by a supermassive black hole (SMBH). The combination of proximity and SMBH makes M87 a prime target to study the outflow nature. The distance and the redsift for M87 are D =16.7 Mpc (Mei et al. 2007) and z = 0.00436 (Smith et al. 2000) correspondingly, which correspond to the angular scale of 1 mas ≈0.08 pc. The supermassive black hole mass obtained from the Event Horizon Telescope image analysis is (6.5 ± 0.2|_stat ± 0.7|_sys) × 10⁹M_⊙ (Event Horizon Telescope Collaboration et al. 2019). Thus, the Schwarzschild radius R_s ≡ 2GM/c² ≈ 1.9 × 10¹⁵cm, which corresponds to 1 mas ≈130R_s.

The relativistic jet in M87 is frequently observed with the Very Long Baseline Interferometry (VLBI) technique. The images typically show an edge-brightened structure and a faint feature on the counter-jet side, south-east of the core. The full track 15 GHz observations which were done using the Very Large Baseline Array (VLBA) clearly show this feature (e.g. Ly, Walker & Junor 2007; Kovalev et al. 2007); see also the collection of MOJAVE images (Lister et al. 2018). This structure is believed to be the counter-jet deboosted by relativistic aberration effects. But there is no spectral information to support this idea. Spectral index maps could also help to better understand the nature of the edge-brightened (e.g. Kovalev 2008) and triple-ridge jet structure (e.g. Asada, Nakamura & Pu 2016; Hada 2017; Savolainen 2021). The high-fidelity Very Large Array (VLA) radio images of the kiloparsec (kpc) scale jet of M87 clearly show a helical pattern (Pasetto et al. 2021). This structure can also be seen in other jets, like 3C 273 (Lobanov & Zensus 2001) or S5 0836+710 (Vega-García, Perucho & Lobanov 2019). The oscillatory pattern of the kpc-scale jet of M87 was studied in Lobanov, Hardee & Eilek (2003), where the observed threads were identified with the Kelvin–Helmholtz (KH) instability happening in the jet. To our knowledge, no helical pattern is found and analysed in the inner, parsec scale jet, and high dynamic range images are needed to study plasma instability.

Kovalev et al. (2007) also showed an extended high linear polarization region but the absence of Faraday rotation maps did not allow the authors to reconstruct a magnetic field direction. Deep observations can significantly complement the polarization studies of the M87 jet and reveal new regions not visible on Zavala & Taylor (2002), Park et al. (2019) and Pushkarev et al. (2023) polarization maps.

These questions can be addressed using multifrequency high-resolution observations. In this paper, we report on dual-frequency full-track high-sensitivity VLBA observations toward the M87 jet which have happened over three consecutive days. The paper is structured as follows. Section 2 describes observations, data calibration and imaging processes. The analysis of the M87 jet structure is given in Section 3. The spectral index analysis is presented in Section 4. In Section 5, we present the linear polarization results. The KH instability and the intrinsic jet parameters analyses are given in Sections 3.2 and 4.3.

2 OBSERVATIONS AND IMAGING

Full Stokes, full track dual-frequency observations of M87 were made in 2009 May using the Very Large Baseline Array (VLBA) involving the Karl Jansky Very Large Array (VLA) single antenna (Y1) and the 100-m Effelsberg Radio Telescope (Eb). The project code from the NRAO archive is BK145. We used Y1 to obtain the shortest baseline possible with this array, thereby improving the extended flux sensitivity, which is crucial in the case of the M87 jet. The Effelsberg was used to obtain higher sensitivity to weak and compact structures. Observations were performed on three consecutive days. This included U band with a central frequency 15.4 GHz (later 15 GHz) on May 22 and 24 as well as X band with a central frequency 8.4 GHz (later 8 GHz) on May 23. The length of each day’s observation was 12 h. The data were recorded in 16 baseband channels (intermediate frequencies, IFs), each of 8 MHz bandwidths using 512 Mbps recording rate and 2-bit sampling. Both right- and left-hand polarizations were recorded simultaneously, giving a total observing bandwidth of 32 MHz in each polarization.

2.1 A priori data calibration

The data were edited and calibrated following traditional methods¹ in AIPS (Greisen 2003). The quasar 0923+392 was used as a fringe finder. Polarization leakage was calibrated with the AIPS task LPCAL using a 1308+326 calibrator. The electric vector polarization angle (EVPA) was calibrated by comparing the measured EVPAs of OJ 287, 0923+392 and 1308+326 with the values obtained from the UMRAO data base (H. Aller & M. Aller)² and the VLA polarization data base³ at 8 and 15 GHz. Due to a paucity of measurements, the values for the exact day of VLBA observations were found by linear interpolation of the UMRAO EVPAs from nearby dates. The final visibility data were averaged over 10 s. The AIPS task UVMOD was applied to the May 22 and 24 visibility data sets, creating a combined (u, v)-file at 15 GHz, to provide comparable to 8 GHz data sensitivity and, consequently, a higher quality of spectral index reconstruction.

2.2 Imaging

We used the Caltech DIFMAP (Shepherd, Pearson & Taylor 1994) software to obtain Stokes I, Q, and U maps from the calibrated visibility data using hybrid mapping in combination with superuniform, uniform, and natural weighting with pixel sizes of 0.08 mas at 15 GHz and 0.16 mas at 8 GHz. The final imaging results yielded >20 000:1 dynamic range at each frequency. The noise level is estimated (Hovatta et al. 2014) as |$\sigma = (\sigma _{\textrm {rms}}^2 + (1.5\sigma _{\textrm {rms}})^2)^{1/2} \approx 1.8(\sigma _{\textrm {rms}})$|⁠, where the first term is connected with the off-source rms noise level of the image σ_rms ≃ 35 |$\mu \mathrm{Jy\, beam}^{-1}$|⁠, and the second term is connected with the uncertainties of the clean (Clark 1980) procedure.

We created two data sets: with and without data on baselines to Eb. All analysis in this paper was done with both data sets. For convenience, the final results in each section of the paper represent an average between results obtained with and without Eb baselines. In sections, where the difference between data sets is significant, the results of only one data set are shown. For example, the resulting spectral index maps with the Effelsberg data show higher uncertainties as compared to the VLBA + Y1 configuration. Thus, the images presented in this paper use the conservative VLBA  + Y1 configuration only as the ‘Atlantic gap’ gives poor (u, v)-coverage between the main array and the Eb, which makes the image reconstruction process more complicated, which can be seen in the spectral index map. In addition, this makes a Kelvin–Helmholtz thread analysis more complicated. The resulting Stokes I images at 8 and 15 GHz are shown in Fig. 1.

The lack of short baselines leads to a significant loss of flux from faint and extended features far from the core, especially at relatively high frequencies such as 8 and 15 GHz. Nevertheless, the data allowed reconstruction of the jet image until ≈ 450 mas and detect an HST-1 feature located at ≈ 850 mas from the core. To achieve a better dynamic range and to compare the results with the previous VLBA and EVN (European VLBI Network) HST-1 detections (Cheung, Harris & Stawarz 2007; Chang et al. 2010; Giroletti et al. 2012), the images were convolved with an 8 × 3 mas beam, corresponding to a VLBA beam at λ 18 cm (Fig. 2).

3 THE M87 JET STRUCTURE

3.1 Jet shape

Jet morphology provides crucial information for understanding the formation and propagation of relativistic outflows. The observed shape of an AGN jet depends on a number of parameters describing physical conditions in the jet plasma and in the ambient medium. So the knowledge of the shape and geometry of the flow can help us estimate intrinsic jet parameters (Pushkarev et al. 2017; Algaba et al. 2017; Kovalev et al. 2020; Nokhrina et al. 2019; Nokhrina, Kovalev & Pushkarev 2020).

In M87, the jet manifests a complex and asymmetric triple-ridge internal structure developing in a predominantly straight, expanding flow. To analyse this structure, we first determine the direction of the jet axis by estimating the overall jet position angle (PA). To avoid potential uncertainties and errors due to the local curvature of the jet near its origin, we estimated the general PA of the jet position of the knot HST-1 in the large-scale image of the jet obtained from our data and presented in Fig. 2. Using the intensity maximum of the HST-1, we obtained P.A. = 293|${_{.}^{\circ}}$|3 ± 0|${_{.}^{\circ}}$|5. Note the significant difference of several degrees from results by Plavin, Kovalev & Pushkarev 2022 due to the different utilized methods of P.A. measurements.

To simplify the further analysis, we rotate the jet images clockwise by ψ = 23|${_{.}^{\circ}}$|3 to coincide the jet axis and relative right ascension and obtain transverse profiles of brightness distribution in the jet. In order to highlight the triple-ridge structure, a stack profile is shown in Fig. 3. It was made by averaging all individual profiles with a step of 0.05 mas at 15–25 mas from the VLBI core. In the jet images produced with the nominal restoring beam, the jet structure up to ≈40 and ≈80 mas, in the 15 and 8 GHz image, respectively can be traced. The 8 GHz image can be further convolved with a larger beam ≈ 1/3 of the jet width in the outer parts (3 mas), which allows us to trace the jet structure with the transverse brightness distribution profile measured up to ≈200 mas separation from the jet origin.

Below, we apply the transverse profiles to measure the evolution of the jet width and quantify the properties of the internal structure of the flow. For both of these tasks, we fitted the obtained profiles with multiple Gaussian components defined as

where x is a transverse distance, the ith component is described by A_i peak amplitude, b_i peak location, and FWHM full width at half maximum, calculated from FWHM_i = 2(2ln 2)^1/2c_i.

The number of Gaussian components used in fitting of each profile was chosen using an χ² analysis. Each profile was fitted by a multi-Gaussian function with a different number of components, where the maximum quantity is 3. Afterwards, the reduced χ² with q degrees of freedom for all fits for a given intensity profile were compared according to Avni (1976). The number of degrees of freedom q = n − p, where n is the number of beams which can be fitted inside the transverse profile, and p is the number of the fitting parameters. Using a cumulative distribution function of χ², the theoretical value of Δ(q, α) ≥ χ² can be calculated, where the probability of this expression equals |$\alpha = 95{{\ \rm per\ cent}}$| (the confidence level). In this method, the q + 1 degrees of freedom model is preferred if the following ratio is satisfied

The width of the jet in a particular image profile was defined as a Full Width at a Quarter Maximum (FHQM) of the fit of the profile. The quarter maximum level was chosen because of the complex and asymmetric transverse structure of the M87 jet, where multiple components can have more than twice the brightness difference. In this case, the usage of the half maximum level can be misleading by measuring the width of one component instead of the whole jet.

After obtaining the jet width as described above, we deconvolved it with the beam projected onto a transverse profile. In our case, one-dimensional (1D) deconvolution procedure was done according to |$w = (\textrm {FWQM}_\textrm {jet~width}^2 - \textrm {FWQM}_\textrm {beam}^2)^{0.5}$|⁠, where w is the deconvolved jet width. In this case, a 1D deconvolved profile can be distorted by a two-dimensional (2D) elliptical beam which influences several profiles asymmetrically. In our data, the major axis of a beam is oriented roughly perpendicular to the jet direction, thus the distortions should be negligible. At the same time, in addition to restoring elliptical beams, we reconstructed images using equivalent area circular beams. Such maps should deliver even less biased results and be easier to analyse. Indeed no significant difference was found between the geometry obtained with circular and elliptical beams.

To investigate the jet geometry, we follow Kovalev et al. (2020) fitting the expansion profile by a power-law function w ∝ (r + r₀)^k, where r is the distance from the apparent VLBI core, r₀ is the distance from the core to the jet base, the expansion index k shows the jet geometry. In Fig. 4, the fitting results of the expansion profile are presented. For this, we used the jet images convolved with a circular beam of the FWHM-size 0.86 and 1.56 mas having an equivalent area of 15 and 8 GHz elliptical beam, respectively. In addition, the geometry of the 200 mas long jet image at 8 GHz convolved with 3 mas circular beam is also presented in Fig. 4. Jet geometry parameters were obtained for images with different array configurations (VLBA + Y1; VLBA + Y1 + Eb), frequencies (8 and 15 GHz), and for extended 200 mas jet at 8 GHz. Finally, all obtained results were averaged to obtain the final M87 jet geometry parameters: k = 0.532 ± 0.008, r₀ = 3.8 ± 1.1 mas. The r₀ estimate by this method has a large uncertainty and is higher than expected; compare with Hada et al. (2012), where the VLBI core separations are about 0.2 mas for 8 and 15 GHz. A possible r₀ overestimation may result from methodological complications. First, the r₀ value in the fitting procedure is sensitive to the absolute value of the jet width in contrast to k, which is sensitive to the relative width. The uncertainty of absolute values is mainly driven by the deconvolution procedure. Second, physical conditions in the jet can change geometry locally, leading to over or underestimations in r₀ estimates. In order to improve future measurements, a dedicated study of a possible deconvolution bias and geometry variations is needed on a large sample of AGN jets (e.g. Kovalev et al. 2020).

3.2 Kelvin–Helmholtz instability

Transverse oscillations in the M87 jet were analysed in previous studies. The 17-yr observations at 43 GHz by VLBA show a shift of the transverse position of the jet on a quasi-periodic 10-yr time-scale that is consistent with the Kelvin–Helmholtz instability (Walker et al. 2018). Recent studies of the M87 jet show a 1-yr period wiggles in multi-epoch 22 GHz KaVA VLBI observations (Ro et al. 2023a). It is still unclear what is the origin of these fast oscillations, but current-driven instability (CDI) was chosen as a preferable mechanism. Cui et al. (2023) analysed the periodicity of the P.A. of the jet using 22–43 GHz EAVN, VLBA, and EATING VLBI observations, where the phenomenon was interpreted as the jet nozzle precession with a period of ≈11 yr.

The data analysed in this article provide the most prominent, sensitive, and extended pc-scale jet structure. In Fig. 1, the brighter jet limb changes from the northern to the southern edge of the jet at a relative RA ≈ −10 mas and vice versa at a relative RA ≈ −25 mas), and this effect is observed in both the 8 GHz and the 15 GHz maps (Figs 1–3). Similar behaviour (changes of the brighter limb and multiple threads inside the jet) can be seen in the pc scales VLBA and RadioAstron observations (Walker et al. 2018; Savolainen 2021), the kpc scales in the VLA (Owen, Hardee & Cornwell 1989; Pasetto et al. 2021) and the HST (Sparks, Biretta & Macchetto 1996) images. The pattern was interpreted as resulting from the KH instability in the flow (Lobanov, Hardee & Eilek 2003).

At the distances within 15–25 mas from the core a triple-ridge structure is observed in 15 GHz map (Figs 1 and 3). The central filament has been observed before in several studies both in proximity to the radio core (Asada, Nakamura & Pu 2016; Walker et al. 2018; Kim et al. 2018) and downstream (Hada 2017) and associated with the spine in the ‘spine-sheath’ jet model (Mertens et al. 2016). The simulations done by Pashchenko et al. (2023) show that the central filament might be found in an edge-brightened model due to a clean imaging artefact. At 15 GHz, the jet width coincides with three beam sizes in that region. This can cause an overlay of beam sidelobes, creating the apparent central filament. This hypothesis is supported by the absence of a prominent central ridge at 8 GHz image, where the jet width is about two beams. For this purpose, we use 8 GHz maps for this analysis.

In Section 3.1, the jet transverse profiles were modelled by Gaussian components. The positions of the peaks of those components obtained for the inner and outer parts of the jet were superimposed together on the intensity image (bottom plots of Fig. 5), thus allowing us to trace the development of the jet structure on scales of up to ∼200 mas from the observed jet origin. The revealed jet pattern suggests the presence of up to three intertwining helical threads inside the jet, which can be related to the development of KH instability inferred for the kpc-scale jet.

Identifications of each of the individual threads were made by requiring a continuous and smooth evolution of thread parameters, such as position, intensity, and width along the jet. Using this approach, a self-consistent picture of the evolving thread-like patterns inside the jet was reconstructed.

In addition to the rich evolution of the internal structure, the jet exhibits a slight bend ≈80 mas which is clearly visible in Fig. 2. To quantify the position and the magnitude of this bent, the 400 mas jet image with 18 cm VLBA beam was used. The image was rotated 23|${_{.}^{\circ}}$|3 clockwise, and all transverse intensity profiles were fitted by a single Gaussian, obtaining the ridgeline (the blue line in the top of Fig. 5).

Both, the observed bent and the pronounced internal structure observed in the jet can result from the development of the KH instability in the flow (Hardee 2003; Lobanov, Hardee & Eilek 2003). The changes in the ridge line and the evolution of the individual threads identified inside the flow can be well described by a three-dimensional (3D) helix in a Cartesian coordinate system (x, y, z):

where z is the distance from the jet origin, ϕ is the phase, and the amplitude and the wavelength depend on the jet radius A ∝ R_jet(z), λ ∝ R_jet(z) (Hardee 2000).

First, we fitted the global jet curvature by a 3D helix projected onto the sky plane, for which the jet angle to the line of sight θ = 17|${_{.}^{\circ}}$|2 (Mertens et al. 2016) was used. The resultant mode can be associated with the helical surface mode (H_s) of the KH instability as it apparently leads to oscillations of the entire jet around its average propagation direction. In order to measure more accurately the parameters of the modes affecting the internal structure of the flow, this mode was subtracted from the component positions before the fitting procedure. A robust identification of the other modes can be obtained using the characteristic wavelength λ* = λ_nm(n + 2m + 0.5), where λ* is the characteristic wavelength, λ_i is the observed wavelength, n is the azimuthal wavenumber, and m is the order of the mode. The pinch (n = 0), helical (n = 1), and elliptical (n = 2) modes are expected to be most prominent in relativistic jets. The order of the mode m determines whether the corresponding perturbation affects the surface (m  = 0) or the interior (m > 0) of the jet. The characteristic wavelength depends only on the physical conditions in the jet; thus, should have a similar value for different modes (Lobanov & Zensus 2001). The final fitting results are presented in Figs 5 and 6 by semitransparent thick pink, red, orange, and blue curves. A zoom into the innermost 60 mas section of the jet is shown in the bottom panel of Fig. 5. The identification criteria discussed earlier were applied to the extracted individual threads resulting in one plausible combination of KH modes shown in Table 1. The red and blue threads have nearly equal wavelengths and ≈180° phase difference, which is expected to result from the two regions of increased pressure and density produced in the jet by the elliptical surface (E_s) mode of instability (Hardee 2000). The orange thread has a much shorter wavelength with a smaller amplitude, which is expected in the case of the body mode. In the case of a first-order elliptical body mode (E_b1), one should observe the accompanied thread. For illustration purposes, the thread with the same parameters, but with the 180° phase difference was plotted as a dashed curve of the same colour in Fig. 5. Due to the lack of high-enough resolution, the accompanied thread was not identified robustly and extracted from the data since it is blended with the other modes. However, it is observed that the dashed curve plotted in Fig. 5 complements the KH model well, describing the intensity images better.

To verify the robustness of the fitting procedure, we generated a synthetic data set reproducing the internal structure with the parameters of all fitted modes. For this, we first created a model image using the parabolic jet model threaded by the KH modes obtained from the fits to the 8 GHz image. Then, we created artificial uv-data using the same uv-coverage and thermal noise, as in the original 8 GHz data set. The same algorithms of the image reconstruction and the KH modes fitting were applied to the artificial data set. The derived parameters of the KH modes were found to be consistent with the ground truth values within 25 per cent. First, this demonstrates that our approach can reconstruct the parameters with the accuracy required to distinguish between two modes (surface and body). Second, the obtained accuracy of the method shows that two threads of the elliptical surface mode are consistent with having the same observed wavelength and the 180° phase offset.

4 SPECTRAL INDEX MAP

The main mechanism of AGN radio jet emission is synchrotron radiation. Assuming a power-law particle energy distribution N(E) dE ∝ E^−pdE, the jet spectrum for the optically thin regions is then described by a distribution I_ν ∝ ν^α with a spectral index α = (1 − p)/2. Therefore, the spectral index distribution provides important information about the physical conditions of different jet regions. Due to phase self-calibration, the information about absolute celestial coordinates is lost. Thus, images at different frequencies need to be well aligned to determine spectral index distributions. For this reason, the normalized 2D cross-correlation (2DCC) method was used (Lewis 1995; Walker et al. 2000; Croke & Gabuzda 2008; Fromm et al. 2013; Hovatta et al. 2014). The method uses optically thin parts of a jet as a reference. Due to their transparency, these parts can be assumed to be located in the same place at both frequencies. So choosing the jet image regions far away from the core and applying a normalized 2DCC, we can align the images at different frequencies.

The spectral index of the M87 jet obtained between 8 and 15 GHz is shown at the top of Fig. 7. Note that, the two used images have comparable sensitivity. Due to the sparsity of the visibility plane coverage, a bias in the resulting spectral index image might be significant; see analysis in Pashchenko et al. (2023). Thus, we estimated the effect of the bias and corrected for it, see details in Section 4.2.4. The final unbiased spectral index maps are shown in Figs 7 and 10. The spectral index profile along the jet axis is presented in Fig. 8.

Quasi-periodic oscillations of width along the jet are observed in the geometry profile Fig. 4. Additionally, Fig. 8 shows deviations of a measured jet width at 8 and 15 GHz from the fitted 15 GHz power-law model. The oscillations exhibit quasi-periodic deviations from the parabolic shape with an amplitude up to 0.7 mas. This result cannot be an imaging effect since the oscillations are seen at both observing frequencies with equal periodicity and at the same distances. Moreover, the contraction near 5 mas is also observed at different epochs (Mertens et al. 2016; Walker et al. 2018, and Kravchenko et al. in preparation). We suggest that this is caused by stationary recollimation shocks. The spectral index map profile in Fig. 8 supports this assumption. There is evidence that recollimation shocks are seen in VLBI observations. In addition, a simulated VLBI total intensity map obtained by computing the radio continuum synchrotron emission using the relativistic magnetohydrodynamic (RMHD) model shows similar periodic contractions (Gómez et al. 2016). Periodic oscillations are also seen in other RMHD simulations (Fuentes et al. 2018; Mizuno et al. 2015). Our results are qualitatively consistent with other observations and simulations, but individual simulations are needed for a detailed comparison.

The M87 jet spectral index analysis results at 8–15 GHz obtained by Hovatta et al. (2014) show a moderate steepening of the spectra within 10 mas down to α ≈ −1.5. The recent results, observed at 22–43 GHZ show a rapid steepening down to α ≈ −2.5 within 10 mas from the VLBI core (Ro et al. 2023b). However, in this paper, there is no significant steepening in the same region. The interesting part of all these data is they all have a similar resolution, thus it is not a consequence of a resolution effect. In order to explain this discrepancy, several reasons are proposed. If the effect is intrinsic, then it can be caused by synchrotron radiation losses (Kardashev 1962; Pacholczyk 1970). This interpretation can be used in the case of NGC 315, where a similar effect was observed in two papers Park et al. (2021) and Ricci et al. (2022). For this source in the same region, the spectral index steepens with the frequency. The other reasonable interpretation is the temporal variations of the spectral index in combination with the multilayer jet structure. The last can produce several electron plasma populations with different energy distributions that will create a broken power-law spectrum. If the effect is instrumental, it can be caused by high-frequency flux losses due to clean bias (Pashchenko et al. 2023). In this case, the residuals of the high-frequency image are convolved with a high-frequency beam, that is smaller than the lower-frequency beam. Thus, if clean is not deep enough, this will reduce the flux value in the extended regions of the high-frequency image. In this paper, to avoid this effect, deep clean was performed. In addition, the unique feature of the data in this paper is twice the difference in exposure time between 8 and 15 GHz. It was done to reduce the noise in 15 GHz data to make it correspond to the 8 GHz data. Thus, the spectral index map produced in this paper may be less affected by this effect, which is observed as the absence of the significant spectral index steepening along the jet.

4.1 Core shift

Due to synchrotron self-absorption, the location of the VLBI core depends on the frequency of observations as |$\vec{\boldsymbol{ r}}_{\mathrm{ c}}(\nu) \propto \nu ^{-1}$| (the so-called core shift effect), under the assumption that the jet is expanding freely, and there is equipartition between the particle kinetic and the magnetic field energy (Lobanov 1998). The observations show that the core shift in the many AGN jets indeed follows dependence |$\vec{\boldsymbol{ r}}_{\mathrm{ c}}(\nu) \propto \nu ^{-1}$| (Sokolovsky et al. 2011). Moreover, using the phase-referencing multifrequency observations of M87, Hada et al. (2011) show that the core shift is described by |$\vec{\boldsymbol{ r}}_{\mathrm{ c}}(\nu) \propto \nu ^{-0.94 \pm 0.09}$|⁠. Assuming |$\vec{\boldsymbol{ r}}_{\mathrm{ c}}(\nu) \propto \nu ^{-1}$|⁠, we can estimate the distance from the 15 GHz core ν₁₅ to the jet origin as

where |$\Delta \vec{\boldsymbol{ r}}$| is the core shift vector. To obtain the core shift vector, we measured the core positions and the alignment shift. To find the core positions, we used the MODELFIT function in difmap: the M87 image was fitted with a model consisting of 2D-Gaussian components. The coordinates of the model components, which have the highest brightness, were used as the core position relative to the phase centre (⁠|$\vec{\boldsymbol{ r}}_{8}$| for 8 GHz and |$\vec{\boldsymbol{ r}}_{15}$| for 15 GHz). The shift between images |$\vec{\boldsymbol{ r}}_{8 \xrightarrow {} 15}$| was obtained using 2DCC at the beginning of Section 4, where images were aligned for constructing the spectral index map. Finally, we calculate the core shift vector as |$\Delta \vec{\boldsymbol{ r}}_{8 \xrightarrow {} 15} = \vec{\boldsymbol{ r}}_{8} - \vec{\boldsymbol{ r}}_{15} - \vec{\boldsymbol{ r}}_{8 \xrightarrow {} 15}$|⁠.

Averaging the results of the measurements for all restored images with different beams and antenna configurations, we estimated |$|\vec{\boldsymbol{ r}}_{\mathrm{ c}}(\nu _{15})| = 0.2 \pm 0.1$| mas. This result is consistent with the phase-referencing core shift measurements Hada et al. (2012). Neglecting the distance from the central SMBH to the jet origin (which is predicted to be about 2.5–4 R_g; Doeleman et al. 2012), |$|\vec{\boldsymbol{ r}}_{\mathrm{ c}}(\nu _{15})|$| represents the separation of the 15 GHz VLBI core from the central SMBH. Finally, using θ = 17|${_{.}^{\circ}}$|2 ± 3|${_{.}^{\circ}}$|3 (Mertens et al. 2016) (Section 4.3), we obtained the deprojected distance from the 15 GHz VLBI core to the black hole: r_deprojected = 0.7 ± 0.3 mas or (5 ± 2) × 10⁻² pc.

4.2 UNCERTAINTIES

4.2.1 Statistical error of the map

The lower the intensity, the lower the signal-to-noise ratio and the larger the errors in the resulting spectral index map. The errors can also be influenced by uncertainties of image shift correction. To estimate the accuracy of the spectral index measurement, a spectral index error map was made. The error consists of the error caused by the uncertainty of the intensity measurement σ_I and the uncertainty of the alignment of the image σ_cs, thus the spectral index random error is |$\sigma _{\alpha } = (\sigma ^{2}_{\textrm {I}} + \sigma ^{2}_{\textrm {cs}})^{1/2}$|⁠.

For the errors caused by the intensity measurement uncertainties

where σ_X, σ_U are the intensity maps errors, and I_X, I_U are the intensity maps at 8 and 15 GHz, respectively. The errors caused by the uncertainty of the image alignment σ_cs were obtained with the help of the algorithm

• Define the uncertainty of alignment. In this case, the error is 2 pixels or 0.025 mas.

• Shift the map in the alignment uncertainty value in four different directions and obtain spectral index maps.

• Subtract modulo for each of the four obtained maps with the original spectral index map.

• Average all four maps.

The bottom image of Fig. 7 displays the final spectral index error map. The error rises to the edges discussed earlier. Also, one can see that the error grows dramatically near the VLBI core. This is caused by the large gradient in intensity images near the core region, so the image alignment uncertainties produce significant errors. If the difference between neighbouring pixels along the jet axis in intensity images is high, then even a 1-pixel shift can dramatically change the resulting spectral index value. The rapid decrease of intensity in images southeast of the core is caused by the synchrotron self-absorption of a jet base and the transition from Doppler boosted to deboosted emission of the jet and the counter-jet. The region of a relatively large error ends at ≈0.5 mas south–east of the core, and the counter-jet area errors are significantly lower. This is another indication that the faint feature south–east to the core is counter jet. The errors in the regions further than ≈50 mas are high, so we decided to show a spectral index map within 50 mas.

4.2.2 Effects of the sparsity of visibility plane

Since the (u, v)-coverage is sparse, the restored beam takes a complicated shape with notable side lobes. The clean algorithm is unable to completely remove them because the flux is subtracted not only from the regions of maximum intensity but elsewhere. This creates a bias on a residual image, and artefacts can appear. After that, a clean model is convolved with a clean beam or a 2D-Gaussian function approximation of a dirty beam. Thus, the side lobes are not taken into account, and convolution errors (bias) can appear. A spectral index map is especially sensitive to this effect.

The coverage of the visibility plane depends on frequency. Even in the case of the same antennae configuration in dual-frequency observations, the pattern of the obtained (u, v) coverage will be the same, but the size will be different. So the coverage at 8 GHz will look like a stretched version of the 15 GHz (u, v) data. Thus, initially, the (u, v)-ranges do not correspond to each other, and it can cause errors or imaging artefacts which can become significant in spectral index maps, such as a steepening of the spectral index along a jet (Hovatta et al. 2014).

To test how inconsistency of the (u, v)-ranges affects the spectral index map, the (u, v)-coverage was clipped to be the same in both bands. Surprisingly, a comparison of original and (u, v)-clipped data made an insignificant difference. Thus, uv-clipping was not applied to obtain the final spectral index map, shown in the paper. The sparsity effects are much more difficult to check since it is based on the fundamental problem of deconvolution. However, it is possible to estimate the influence of this effect on the spectral index map. Here are the steps of bias checking applied in the paper:

• Create an artificial model out of the clean components whose structure is similar to the real jet in its form and intensity.

• Create visibility functions with the same model but different (u, v)-coverage (at 8 and 15 GHz), using the UVMOD task in aips.

• Do imaging in difmap.

• Obtain a spectral index map.

The obtained spectral index map should have zero values in a range of all maps, in the absence of a bias. In the real-intensity image, we can see the edge-brightening effect. To describe the jet and create a synthetic model, we assumed that the jet is hollow and used simple geometrical shapes. Since the edge-brightened jet is transparent and the core is opaque, we can define the core as a circle and the limbs as two expanding rails. The overall jet model was also inclined in a corresponding P.A. This model was chosen to bring the simulation closer to the observations.

Analysing the results with different beams, we conclude that the spectral index map convolved with the 8 GHz elliptical beam has minimal influence on the convolution effect with the error Δα ⪅ 0.13 in the central ridge region. Since, in the case of the VLBA + Y1 + Eb data set, the long-baseline coverage is poor, the problem of the convolution effect is more significant compared to the data with flagged Eb. It is seen that with increasing the resolution, the bias amplitude increases too. But its pattern is not truly aligned to the flattened spectrum in a ridge which was noticed in the spectral index map (Fig. 7). The results of modelling also showed that the amplitude of the bias is less than the flat spectrum region values. Thereby, according to the analysis done above, the convolution effects do not affect the spectral index map.

4.2.3 Spectral index within the inner 6 mas

To examine in detail the spectral index distribution in the inner 6 mas, we used a circular beam obtained from an average of 8 and 15 GHz elliptical beams equivalent area. The resultant image revealed a double structure. The flattening of the spectral index in the image coincides with the peaks of total intensity at 8 and 15 GHz, while the spectral index steepens toward the jet edges.

We also made use of recent VLBA observations performed on 2018 April 28 simultaneously at 24 and 43 GHz, which are presented in Kravchenko et al. (2020). The full-intensity images at two frequencies were convolved with a common 24 GHz equal-area circular beam and were aligned using a 2DCC procedure. Besides, the (u, v)-coverage was matched for this pair of frequencies. The resultant spectral index map revealed a two-humped structure which corresponds well to the result at 8–15 GHz observations.

These two independent maps reveal a similar structure, see for details fig.1 in Pashchenko et al. (2023), but is it real? A detailed analysis was performed to check its significance and concluded that the two-humped structure is actually a product of an imaging bias. Due to its importance, the analysis was presented in a dedicated paper Pashchenko et al. (2023).

4.2.4 Bias correction of the spectral index images

Pashchenko et al. (2023) employed a series of simulations with various jet brightness models and dual-frequency VLBI data sets. They found that the spectral index maps of the M87 parsec scale jet are heavily affected by systematical effects. For the data set analysed in this paper, these effects flatten the spectrum in a series of stripes along the jet. This is similar to the observed image: two stripes near the core turn into the central and two outer stripes of the spectral flattening further out at r ≈ 10 mas. The outer stripes are shown at the top of Fig. 7 by red horizontal stripes at 3 and −3 mas from the central stripe. The simulations of the data set reveal that the systematic spectral index effects trace the bias of the low frequency (8 GHz) Stokes I image. At the same time, the bias of the 15 GHz Stokes I image is down-weighted by convolving with a lower frequency clean beam. This made it possible to successfully compensate the spectral index bias in simulations by recreating such effects in the high frequency (15 GHz) Stokes I image. Indeed, assuming that the Stokes I bias is small and does not heavily depend on the brightness distribution, which follows from the simulations, it cancels out in the expression for the spectral index. The procedure of the spectral index bias correction consists of three steps. First is an interpolation of the original clean model at the 15 GHz on the uv-points of the 8 GHz data set, creating visibility data of the 15 GHz model with 8 GHz (u, v)-coverage. Next, the resulting data set is imaged in the same way as the original 8 GHz data set. In the final step, a bias-compensated spectral index map is produced (Fig. 7). These results heavily contradict previous conclusions. It is evident, that simple intensity models, as they are used in Section 4.2.2, are insufficient, and a more sophisticated approach is required to analyse the convolution effects and correct the bias in the spectral index map.

4.3 Jet to counter-jet flux ratio

The jet-to-counter-jet flux density ratio is a direct measurable value useful for estimating the basic physical parameters of relativistic jets such as viewing angle θ and Lorenz factor Γ. We measured the flux density of the counter-jet in the following way. The jet structure was modelled with several 2D Gaussian components, the brightest and closest to the phase centre was identified as the core. The core component only was convolved with the restoring beam. Then the obtained core image was subtracted from the original one, giving us well-separated jet and counter-jet structures. To get the flux ratio, we obtained fluxes from geometrically corresponding regions. For this, the length of the detected counter-jet was measured using the estimated in Section 3.1 distance to the central SMBH. Thus, the boundary between the jet and the counter-jet was defined.

To measure viewing angle θ and Lorenz factor Γ, we used

where δ_jet, δ_cjet are the Doppler factors of the jet and the counter-jet, β_app is the apparent jet speed (Boettcher, Harris & Krawczynski 2012). In the obtained images, continuous emission is observed, thus |$\delta _\textrm {jet}/\delta _\textrm {cjet} = ({F_{\textrm {jet}}}/{F_{\textrm {cjet}}}) ^{\frac{1}{2-\alpha }}$|⁠, where F_jet and F_cjet are the fluxes of the jet and the counter-jet (Scheuer & Readhead 1979). The estimations require a spectral index and apparent velocities. Since the length of observed counter-jet is roughly 4 mas, the corresponding 4 mas downstream the jet were considered. We take the spectral index of the jet from this region (α = −0.5 ± 0.2), which is the median value of the distribution in the map (Section 4). The jet speed was chosen as β_app = 0.34 ± 0.35, which was estimated as an average of the apparent velocities of the jet features within 4 mas from the core (Mertens et al. 2016). Finally, we obtained results with different array configurations (VLBA + Y1; VLBA + Y1 + Eb) and clean beams (elliptical; round). The averaged outcome is summarized here

5 LINEAR POLARIzATION AND ROTATION MEASURE

Polarized emission moving through magnetized plasma can be affected by Faraday rotation. The effect can depolarize the light and change its EVPA. A relativistic jet is surrounded by a slow cocoon, optically active matter located in the magnetic field (Savolainen et al. 2021), thus, in order to investigate the magnetic field in a jet, the effect should be taken into account (Burn 1966; Hovatta et al. 2012). The rotation measures RM of extragalactic sources and can be used as a probe of the intergalactic medium to correct EVPA for Faraday rotation. The RM has been determined by the linear fit

where ϕ_obs is the observed polarization angle, ϕ₀ is the intrinsic EVPA and RM is the rotation measure. The nπ ambiguity of the ϕ_obs measurements leads to uncertainties in RM estimations. Thus,

To solve the nπ problem successfully, at least three frequencies are needed (Simard-Normandin, Kronberg & Button 1981). Unfortunately, in our case, there are only two of them, so the problem of the nπ ambiguity is especially important. In Park et al. (2019), it was shown that the rotation measure of the jet within 15 mas is changing in time though slowly, as compared to the measurement errors. So one can assume the stability of RM over a long period of time, at least 20 yr. Since the only unknown parameter in equation (8) is n, it can be selected to match RM obtained using our data with the average value of RM, according to Zavala & Taylor (2002) and Park et al. (2019). The average RM value for the period of 1995–2015 is ≈−4500 rad m⁻². In addition, we calculated the RM value, which provides full depolarization with the condition of the π EVPA rotation. This makes an upper bound for |n|. Finally, the average RM and the maximum |n| were used to solve the nπ ambiguity. For the case of n = 0, the resulting mean RM ≈ −4000 rad m⁻² is consistent with previous measurements (Zavala & Taylor 2002; Park et al. 2019). We assess the significance of selecting the integer value of n by comparing the shift ΔRM associated with a change of Δn = 1 with the errors |$\sigma _{\textrm {RM}}^{\textrm {old}} \approx 1000~\mathrm{rad\,m}^{-2}$| in the average RM ≈ −4500 rad m⁻² value reported in the literature, which we use as a reference point. The errors |$\sigma _{\textrm {RM}}^{\textrm {old}}$| should be less than the shift ΔRM, otherwise, multiple integer solutions for n would exist, leading to ambiguity. In our particular case, a change of Δn = 1 corresponds to a shift of ΔRM ≈ 2600 rad m⁻², which leads to a change of EVPA by Δϕ ≈ 60°. This supports the choice of n, since the errors of RM measurement in Zavala & Taylor (2002) and Park et al. (2019) are two times smaller. The galactic RM correction was not applied to the M87 RM map, since it is negligible for 8 and 15 GHz (Kravchenko, Kovalev & Sokolovsky 2017). We present the reconstructed deep M 87 dual-frequency polarization images and the Rotation Measure map in Fig. 9. Assuming that the linear polarization angle lies in the emitting particle orbit plane, i.e. perpendicular to the magnetic field lines, it is possible to reconstruct the magnetic field lines direction. In Fig. 9, we observe that EVPA is oriented perpendicular to the jet direction near the edges, so the magnetic field lines are parallel to the jet axis. The opposite can be seen in the central jet regions. These results are consistent with the observations made by Zavala & Taylor (2002) and Park et al. (2019), where the perpendicular orientation of EVPAs in the northern edge of the M87 jet is observed.

The Rotation Measure map (Fig. 9) shows the difference between the northern and the southern limb of the jet. The value changes from approximately −5000 to −2000 rad m⁻². Having the mean value of −4070 ± 1030 rad m⁻², the individual pixel errors are smaller ∼300 rad m⁻². The significance of the RM gradient should be carried out carefully. Previous studies proposed several criteria for RM gradient establishment (Taylor & Zavala 2010; Kravchenko, Kovalev & Sokolovsky 2017). Thus, the gradient of RM in the map can be considered significant, since several criteria are satisfied:

• The M87 jet has a rich resolved transverse structure: more than three clean beams across the jet.

• A change in RM is more than three times greater than the typical error.

• The gradient is located in an optically thin region.

• The RM change is monotonic and smooth.

Although the results show a smooth distribution of RM, it is important to note that the conclusion of the significance of the RM gradient was made under the assumption of a uniform value of n across the whole jet due to the limited frequency coverage.

6 DISCUSSION

6.1 Jet morphology, spectral index distribution, and KH instability

Morphological and spectral properties of the jet can be compared by overlaying the KH model threads from Section 3.2 onto the spectral index map from Fig. 7. This overlay is shown in Fig. 10 where one can see that the locations of the KH threads generally coincide with the regions of the flatter spectrum. This positional coincidence is in agreement with the theoretical expectation (Hardee 2000; Perucho et al. 2006) for the elliptical surface mode and different body modes of KH instability to produce higher pressure regions near the jet boundary and in its interior. In these regions, the combination of intrinsic heating (Lobanov, Hardee & Eilek 2003; Hardee & Eilek 2011) and decreasing optical depth at higher frequencies should manifest itself with an apparent flattening of the synchrotron spectrum.

Basic physical parameters of the jet and ambient medium can be derived from the mode identification described in Table 1 by applying linear analysis of the KH instability (Hardee 2000). This allows us to estimate the jet Mach number, M_jet, the ratio, η = h_jet/h_ex of specific enthalpy in the jet and the external medium, and the respective speeds of sound in the external medium and the jet, a_ex, a_jet. In these calculations, we use the estimates of the jet viewing angle, θ = 17|${_{.}^{\circ}}$|2 ± 3|${_{.}^{\circ}}$|3, apparent speed v_app = (2.31 ± 0.14)c and the pattern speed of instability v_w/c = 0.34 ± 0.21 (Mertens et al. 2016) and obtain M_jet  = 20 ± 17, η = 0.3 ± 0.5, a_jet = 0.05 ± 0.03, a_ex = 0.03 ± 0.01. The resulting estimated Mach number and the enthalpy ratio are higher than what is typically expected in relativistic jets (Rossi et al. 2008), for example in a relativistic jet simulation for 3C 273 and 3C 31 show M_jet ∼ 3, η ∼ 0.02 (Perucho, Lobanov & Martí 2005; Perucho & Martí 2007). Conversely, the respective sound speed in the jet plasma is lower than expected from those models. This apparent discrepancy may be explained by the effect of the dynamically important magnetic field which may affect the instability pattern (Hardee 2007) but is not included in the plain KH models. Alternatively, it can result from underestimating the true jet speed from the apparent measured speed. At the viewing angle of the jet in M 87, the detection of plasma motion at Lorentz factors ≳ 3.38 may be hindered by the differential Doppler boosting. Allowing for this effect, we can adopt the highest detected speed in M87 (v_app ≈ 6c; Biretta, Sparks & Macchetto 1999) for our estimates, which results in M_jet ≈ 5, η ≈ 0.014, and a_jet ≈ 0.24. These values are in good agreement with the estimates obtained for the kiloparsec-scale jet in M87 (Lobanov, Hardee & Eilek 2003).

The edge brightening of the jet which can be seen in Fig. 1 and also has been reported previously (Asada, Nakamura & Pu 2016; Walker et al. 2018; Kim et al. 2018; Janssen et al. 2019) can be explained by the stratification of a jet in velocity, density, and internal energy (Walker et al. 2018; Bruni et al. 2021) or by the presence of dynamically important magnetic field (Janssen et al. 2021; Kramer & MacDonald 2021). Each of these effects should also provide a flattening of the spectrum toward the geometrical axis of the jet. In the first case, an interaction between layers in a velocity-stratified jet can heat the plasma, increasing the internal energy and the spectral index in the spine region. In the second case, the spectral index can increase due to opacity variations induced by the magnetic field (Clausen-Brown, Lyutikov & Kharb 2011). A combined action of these factors is also plausible (Kim et al. 2018), but our present results do not warrant a more detailed assessment of the physical conditions causing the observed edge brightening of the jet.

The observed spectral index distribution presented in Fig. 7 indicates that the emission in the compact region at the base of the jet is optically thick at frequencies below 15 GHz, most likely due to synchrotron self-absorption. The estimated errors of the spectral index are high in this region, making it difficult to locate the transition from optically thick to thin emission. With these errors taken into account, the transition may be occurring in the jet anywhere between ≈−1.1 mas and ≈+1.5 mas axial separation from the coordinate origin. On the counter-jet side, the spectral index is reliably measured at separations exceeding ≈1.5 mas, allowing for comparison of the emission properties in the jet and the counter jet. The median values of the spectral index distribution for the jet and the counter-jet regions are α_j = −0.5 ± 0.2, α_cj = −0.8 ± 0.2. The spectral index can also be estimated for the HST-1 feature, using the measured integral flux densities F_8GHz ≈ 19 mJy and F_15GHz ≈ 8 mJy, which yields α_HST-1 ≈ −1.4. Although this estimate should formally indicate pronounced synchrotron ageing of the emission, it should be viewed with extra caution because, at the location of HST-1, the 15 GHz flux density may be underestimated which would cause spurious steepening of the measured spectral index.

6.2 Magnetic field structure

Similarly to the total intensity maps, the images of fractional linear polarization in Fig. 9 are also edge brightened. For the optically thin emission, this effect can result from apparent or true depolarization in the central regions of the jet caused by a helical magnetic field (Gabuzda 2021) or a turbulent plasma flow (Marscher 2014). The possibility for such a structure to result from a clean imaging artefact due to the residual uncleaned polarized flux (Pushkarev et al. 2023) is not likely, as we employed a deep clean in this work Section 4.

The variations of the EVPA shown in Fig. 9 can be reconciled with in a hollow, edge brightened jet (Frolova, Nokhrina & Pashchenko 2023) seen at a viewing angle ∼1/Γ and threaded by a large scale helical magnetic field with the pitch angle >45° in the plasma frame (Lyutikov, Pariev & Gabuzda 2005). Similar EVPA morphology can be obtained for the force-free reverse field pinch (Clausen-Brown, Lyutikov & Kharb 2011) and pure helical magnetic field model, both analytically (Murphy, Cawthorne & Gabuzda 2013; Butuzova & Pushkarev 2023) and numerically (Kramer & MacDonald 2021).

Some indications in favour of the presence of a substantial regular magnetic field in the jet may also be found in the transverse gradient of the Faraday rotation measure observed in the jet. Although transverse gradients of the rotation measure may as well be produced by the differences in the density of a Faraday screen of thermal electrons, which surrounds the jet, the polarization images of the large jet in M 87 (Pasetto et al. 2021) also suggest the presence of a helical magnetic field. The observed anticorrelation between the side of the jet with higher degrees of polarization and the side with higher RM magnitude Fig. 9 could also point at the helical magnetic field (Gabuzda 2021).

7 SUMMARY

This paper presents an investigation of the physical properties of the parsec-scale jet in M 87 obtained from imaging the jet with augmented VLBA at 8 and 15 GHz. The main observational results of this work are:

• Images of total and polarized intensity are obtained at each frequency. The total intensity images, reaching a record dynamic range >20 000:1 show edge-brightening, faint counter-jet, HST-1 knot and reveal helical threads in the jet.

• clean bias, which strongly affects a spectral index map obtained from the individual total intensity images, is identified and corrected for. The bias-corrected spectral index map demonstrates a complex pattern.

• The linear polarization maps uncover the change of magnetic field lines from edges to the jet’s centre.

• The rotation measure map shows a significant gradient perpendicular to the jet direction.

A detailed analysis of the observational information obtained from the VLBA data yields the following conclusions:

• The helical threads observed in the jet can be explained by the Kelvin–Helmholtz instability in the jet. This interpretation is also supported by the spectral index map, where the flattening of the spectra traces well the observed helical threads.

• The edge brightening observed in the total intensity and fractional polarization can be interpreted either by transverse velocity stratification of relativistic plasma or by a large-scale helical magnetic field.

• The faint structure south–east to the core is confirmed as the counter-jet. According to the spectral index map analysis, this feature is optically thin α_cj = −0.8 ± 0.2 and hence this region cannot be the jet origin.

• Intrinsic physical jet parameters are estimated from modelling the observed jet structure: the jet viewing angle θ = 17° ± 8°, Lorenz-factor Γ = 1.2 ± 0.1, expansion index k = 0.532 ± 0.008, Mach number M_jet  = 20 ± 17, jet to ambient medium density ratio η = 0.3 ± 0.5 and the deprojected distance from the VLBI core to the SMBH r_c = (5 ± 2) × 10⁻² pc.

Our results reveal that Kelvin–Helmholtz instability starts to develop in the regions relatively close to the central engine of the jet (∼10²–10⁴R_g). This specific region corresponds to the jet’s formation and collimation zone (e.g. Nakamura et al. 2018; Kovalev et al. 2020), emphasizing the significant role of plasma instability in the jet morphology and evolution. Consequently, it is crucial to consider plasma instability when studying jet properties. The next step for the investigation will be the study of the temporal evolution of the helical pattern in the jet. However, due to the limited sensitivity of modern telescopes, it is challenging to track the helical pattern in survey mode. Thus, only rare full-track observations can be utilized for the analysis. The situation for M87-type jet studies will be greatly improved by ngVLA (Selina et al. 2018; Murphy et al. 2018). In addition to the sensitivity, short baselines will provide an opportunity to close the gap between parsec and kilo-parsec scale jets.

Facilities: VLBA, VLA, Effelsberg telescope.

Software:astropy (Astropy Collaboration et al. 2013, 2018, 2022), matplotlib (Hunter 2007), numpy (Harris et al. 2020), python (Van Rossum & Drake 2009), scipy (Virtanen et al. 2020).

ACKNOWLEDGEMENTS

We thank an anonymous referee for careful reading and valuable feedback, that have greatly enhanced the clarity of this paper. We thank M. L. Lister, D. C. Homan, K. I. Kellermann, E. E. Nokhrina, M. Janssen, J. Livingston, J. I. Nikonova, and J. Röder for useful discussions and comments as well as Elena Bazanova for language editing. This research was supported in part by the Russian Science Foundation (project 16-12-10481). This work is part of the M2FINDERS project which has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 Research and Innovation Programme (grant agreement no. 101018682). ASN received financial support for this research from the International Max Planck Research School (IMPRS) for Astronomy and Astrophysics at the Universities of Bonn and Cologne. The National Radio Astronomy Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc. This work is partly based on the observations with the 100-m telescope of the MPIfR (Max-Planck-Institut für Radioastronomie) at Effelsberg. This research has made use of the data from the University of Michigan Radio Astronomy Observatory which has been supported by the University of Michigan and by a series of grants from the National Science Foundation, most recently AST-0607523.

DATA AVAILABILITY

The correlated data underlying this article are accessible from the public NRAO archive. The 8 and 15 GHz Stokes I and polarization intensity images, the original and the corrected for bias spectral index maps are available from the article and as FITS files in its online supplementary material. The maps are also available from the public MOJAVE archive.⁴

Footnotes

References

Author notes