On the M87 jet structure near the central engine

At present, there is no doubt that relativistic jets observed in active galactic nuclei pass from highly magnetized to weakly magnetized stage, which is observed as a break in the dependence on their width d jet ( z ) on the distance z to the central engine. In this paper, we discuss the possibility of observing another break, which should be located at shorter distances. The position of this break can be associated with the region of formation of the dense central core near the jet axis which was predicted both analytically and numerically more than a decade ago, but has not yet received sufficient attention. In this case, the observed width should be determined by the dense core, and not by the total transverse size of the jet. The calculations carried out in this paper, which took into account both the transverse electromagnetic structure of the jet and the change in the spectrum of emitting particles along its axis, indeed showed such behaviour. We also found the evidence of the predicted break in the jet expansion profile using stacked 15 GHz VLBA image of M87 radio jet and constrain the light cylinder radius.


INTRODUCTION
A real breakthrough in the study of the internal structure of relativistic jets emanating from active galactic nuclei (AGN), achieved over the past decade thanks to the high angular resolution obtained using very long baseline interferometry (VLBI) observations (Kovalev et al. 2007;Lister et al. 2009;Homan et al. 2015;Lister et al. 2016;Mertens et al. 2016) allows us to investigate directly their internal structure.In particular, it give us the direct information about the dependence of the jet width djet(z) on the distance z to the central engine.It finally became possible to compare the results of observations with the predictions of a theory, the basic elements of which have long been constructed (Blandford 1976;Lovelace 1976;Ardavan 1976;Blandford & Znajek 1977;Begelman et al. 1984;Camenzind 1986;Heyvaerts & Norman 1989;Camenzind 1990;Takahashi et al. 1990;Pelletier & Pudritz 1992;Beskin & Pariev 1993;Heyvaerts 1996).
We recall that the nature of relativistic jets from AGN is believed to be associated with highly magnetized magnetohydrodynamical (MHD) outflows generated by fastly rotating supermassive black holes (Krolik 1999;Camenzind 2007;Beskin 2010;Meier 2012).Within this approach, the key role is played by the poloidal magnetic field generated in the accretion disc as within this model both plasma outflow and the energy flux propagate along magnetic field lines from the central engine to active regions.As was shown by Blandford & Znajek (1977), a strongly magnetized flow (in which the main role in energy transfer is played by the electromagnetic field, i.e. the Poynting vector flux) is really capable of taking the energy away from a rotating black hole and transmitting it to infinity.As a result, thanks to numerous works devoted to both analytical and numerical analysis of MHD equations (Chiueh et al. 1991;Appl & Camenzind 1992;Eichler 1993;Bogovalov 1995;Ustyugova et al. 1995;Beskin 1997;Lery et al. 1999;Beskin & Malyshkin 2000;Vlahakis & Königl 2003;McKinney 2006;Beskin & Nokhrina 2006;Komissarov et al. 2007;Beskin & Nokhrina 2009;Romanova et al. 2009;Lyubarsky 2009;Beskin 2010;Tchekhovskoy et al. 2011;Porth et al. 2011;McKinney et al. 2012;Potter & Cotter 2015;Beskin et al. 2017), a consistent model was constructed, which is now generally accepted.
Let us note from the very beginning the following key points characterizing the MHD model of relativistic jets; the assumptions of the model we use will be described in detail in Section 2. First, it should be stressed the role of the ambient pressure Pext which actually determines the transverse size of the jet d(z) and, hence, their inner structure (Beskin 1997;Lyubarsky 2009).In particular, it is the pressure dependence Pext(z) on the distance z that is to determine the dependence of transverse size of the jet on the distance from the central engine.
Second, as was found both analytycally (Beskin & Nokhrina 2006;Lyubarsky 2009) and numerically (McKinney 2006;Komissarov et al. 2007;Tchekhovskoy et al. 2011;Porth et al. 2011), for a sufficiently high external pressure Pext (i.e.not so far from the central engine), the longitudinal poloidal magnetic field within the jet is to be homogeneous.However, as the external pressure decreases, a denser core begins to form in the very centre of the jet, while the poloidal magnetic field Bp and the number density ne of the outflowing plasma begin to decrease significantly with distance ϖ from the jet axis.But the flow still remains highly magnetized.Recently this result was comfirmed by Beskin et al. (2023b) for relativistic jets with zero velocity along the axis.
Third, as in the magnetically dominated flow the energy of particles increases with the distance ϖ from the jet axis, at large enough distances z corresponding to lower external pressures (where the jet width djet(z) also becomes large enough), inevitably, saturation should occur, when almost all the electromagnetic energy flux (which, as was already stressed, is dominated near the black hole) will be transferred to the energy of outflowing plasma (see Sect. 2 for more detail).
Currently, there are direct indications that saturation occurs at the distances of several parsecs (10 5 -10 6 of gravitational radii Rg = GM/c 2 ) from the central engine (Nokhrina et al. 2020).Such a transition from highly magnetized to weakly magnetized regime was predicted to be observed as a break in the dependence on the width of the jet djet(z) on the distance z (Beskin et al. 2017;Nokhrina et al. 2019).This results in the observed parabolic streamline in the upstream region, and the jet shape becomes conical downstream a few 10 5 Rg.This structural transition of a jet shape was first discovered in M87 (Asada & Nakamura 2012) and the similar transition was lately found in a number nearby sources (Tseng et al. 2016;Hada et al. 2018;Nakahara et al. 2018;Akiyama et al. 2018;Kovalev et al. 2020;Nakahara et al. 2020).
Remember that for M87 such a nature of the break has also been discussed by Mertens et al. (2016).But they talked about a feature located at a distance about 10 mas, i.e., hundreds of times smaller.Therefore, it is of interest to discuss the question of whether the feature at the very base of the jet in the M87 is not connected with the formation of the central core, i.e. with the transition from the jet with homogeneous poloidal magnetic field to the jet having a dense central core.An indication of such a break can be found in the work of Hada et al. (2016).However, this issue has not yet been discussed in detail.
Thus, in this paper, we investigate the internal structure of relativistic jets at small distances z from the central engine.
Having decided on the model of the internal structure of electromagnetic fields (which, within the framework of the MHD approach, allows us to determine the plasma number density as well), we then find the observed width of the jet, the radiation of which is associated with the synchrotron radiation of outflowing particles.Our goal is to find out whether, in addition to external break in the dependence of the width of the relativistic jet djet(z) on the distance z, at which the transition from magnetically dominated to particle dominated flow occurs, there is also an internal break in which the central core near the jet axis begins to form.
Let us recall that the source of radiating leptons are currently unknown.Since the seminal paper by Blandford & Znajek (1977) it has been believed that they can be formed as a result of the collision of thermal gamma rays emitted by the highly heated internal regions of the accretion disk.Nu-merical simulations for pair production from the thermal seed radiation of plasma in a vicinity of a black hole event horizon by Mościbrodzka et al. (2011) provide the number of electronpositron pairs enough to explain the jet particle number density implied by a core-shift method (see e.g.Nokhrina et al. 2015 and references therein) further downstream.It is assumed (see e.g.Blandford & Königl 1979;Konigl 1981) that the jet emission in radio is due to synchrotron self-absorbed radiation.In optically thin regime far from a radio core the observed power-law spectrum (Hovatta et al. 2014;Kutkin et al. 2014) can be explained by the power-law energy distribution of emitting particles.The standard approach usually separates the MHD jet modelling and the mechanisms of emitting plasma acceleration and formation of an expected power-law spectrum (see e.g.Kirk et al. 1994;Sironi et al. 2013;Ostrowski 1990Ostrowski , 1998)).This is due to the complexity of the full problem of merging MHD with physical kinetics (particle-in-cell method, see Kagan et al. 2015, for a review) in one scheme.Sironi et al. (2016) superimposed the development of blobs of highly relativistic plasma, forming due to reconnection, on a jet-like structure.But modelling the spectral flux maps to explain the observational data is achieved at the moment by imposing the assumed non-thermal plasma distribution onto the jet structure (Porth et al. 2011;Potter & Cotter 2012, 2013;Pashchenko & Plavin 2019;Fromm et al. 2019;Ogihara et al. 2021;Kramer & Mac-Donald 2021;Frolova et al. 2023).In this paper, we follow the latter approach.
The paper is organized as follows.In Section 2 the basic relations of the MHD theory are formulated, allowing a detailed description of the internal structure of relativistic jets.Based on this theory, in Section 3 we determine the transverse size of relativistic jets near their base, which allows us to find the position of the internal break.Finally, in Section 4 we employ the observational data on the break and discuss the imposed constraints.

INTERNAL STRUCTURE OF RELATIVISTIC JETS
First of all, let us dwell in detail on our model of the internal structure of relativistic jets.It is based on the results obtained by Beskin & Nokhrina (2006), according to which, for strongly collimated jets, one can consider their internal structure as a sequence of cylindrical flows depending on a jet radius ϖ only.For such cylindrical configurations, the Grad-Shafranov (GS) equations describing internal structure of ideal MHD flows (see Heyvaerts & Norman 1989;Camenzind 1990;Pelletier & Pudritz 1992;Beskin & Pariev 1993;Heyvaerts 1996 for more detail) can be easily integrated (Beskin 1997;Beskin & Malyshkin 2000;Lyubarsky 2009) as it reduces to two ordinary differential equations of the first order.As a result, by specifying four integrals of motion at the base of a magnetically dominated flow (see below), it is possible to obtain transverse profiles of electromagnetic fields (and, within the framework of this approach, hydrodynamic speed and particle number density as well) at an arbitrary distance from the central engine.
A detailed discussion of the model used can be found in our previous works (Beskin & Nokhrina 2009;Nokhrina et al. 2015;Beskin et al. 2017;Chernoglazov et al. 2019).Here we recall only two significant circumstances.First, as already noted, one of the important features of the model is that at sufficiently large distances from the central engine, a denser core begins to form in the very center of the jet, and the poloidal magnetic field Bp and the number density ne of the outflowing plasma begin to decrease significantly with distance ϖ from the jet axis.This structure, shown in Figure 1, was also confirmed by numerical simulation (Komissarov et al. 2007;Tchekhovskoy et al. 2011;Porth et al. 2011).
Second, despite the cumbersome form of the GS equation, rather simple asymptotic solutions have been formulated, which will be enough for us to state the main points of this paper.As was already stressed, the structure of the poloidal magnetic field within relativistic jets outflowing from active galactic nuclei substantially depends on the ambient pressure Pext which, in turn, depends on the distance z from the central engine.Therefore, in what follows, we consider all quantities as functions of the distance z.
As a result, with a good accuracy, poloidal magnetic field Bp can be represented as , Bp(0) < Bcr. (2) Here is the core radius where the integration constant γin ∼ 1 corresponds to the Lorentz-factor of a flow along the axis, and RL = c/Ω0 is the radius of the light cylinder (the constant Ω0 will be defined below).
In other words, at small distances from the central engine, poloidal magnetic field remains actually constant (so that B2 0 /8π ≈ Pext)1 , while at large distances it decreases with distance from the axis ϖ.As was already stressed, this result has been repeatedly confirmed by numerical calculations (Komissarov et al. 2007;Tchekhovskoy et al. 2008;Porth et al. 2011).2) (black dashed lines) with the structure of a poloidal magnetic field (blue solid lines) obtained by solving a system of two ordinary differential equations describing internal structure of a cylindrical jet.In this figure, we terminated the solid lines of numerical solution at smaller radii so that the analytical approximation may be better seen.We plot three jet crosscuts for a conical jet with intrinsic opening angle 10 • .
Moreover, this structure is realized for zero hydrodynamical velocity along the jet axis (γin = 1) as well (Beskin et al. 2023b).As shown in Figure 2, analytical approximation ( 1)-( 2) perfectly reproduces the structure of a poloidal magnetic field obtained by solving a system of two ordinary differential equations describing internal structure of cylindrical jet (see Beskin & Nokhrina 2009;Lyubarsky 2009;Chernoglazov et al. 2019 for more detail).Wherein, the characteristic magnetic field Bcr, at which the formation of the central core begins, can be presented as where is so-called Michel (1969) magnetization parameter.Here µ ≈ mec 2 is the relativistic enthalpy, and η = nup/Bp is the particle-to-magnetic flux ratio 2 .Due to numerous evaluations (see, e.g.Nokhrina et al. 2015), for most jets from active galactic nuclei the Michel magnetization parameter σM ∼ 10-50.For M87 we assume σM ∼ 20.In turn, the magnetic field BL is determined from the relation where Ψtot is the total magnetic flux within the jet.In what follows we assume BL = 100 G, which corresponds to RL = 10Rg.Further, the exponent α changes from 0 to 1 as the ambient pressure decreases to values Peq = B 2 eq /8π, where Beq = σ −2 M BL.In this region, the magnetic field on the axis B0 weakly depends on Pext.In all this magnetically dominated region, the asymptotic behaviour γ(x) = x = Ωϖ/c is fulfilled.At an even lower ambient pressure Pext < B 2 eq /8π, the flow becomes partically dominated.In this domain α ≈ 1.
It should be emphasized that in order to determine the explicit form of dependences B0(z) and α(z), it is necessary to specify either the dependence of ambient pressure Pext(z) on z or the shape of the jet boundary.However, since the ambient pressure is not known with sufficient accuracy, below we use the second approach.It is based on the fact that the parameter of the problem considered by Beskin & Nokhrina (2006) was not the distance from the central engine z, but the transverse radius of a jet rjet.Therefore, by specifying the shape of the jet, one can also determine the dependence of the parameters on the coordinate z.Below we use the conical model (see Figure 1), which seems to be the most reasonable at the very base of the jet.As we see below, this assumption does not contradict the observed parabolic shape of the jets.
Using now the results of the work of Beskin & Nokhrina (2006), we obtain the following approximation for the functions α (z) and B0 (z) defining the poloidal magnetic field (1)-( 2) in the case of a conical structure for z > zcr α (z) ≈ 0.05 ln (z/zcr) , ( 7) In ( 8) B0 (z) is determined from the conservation of the total magnetic flux.
From here it can be seen that at short distances, before the formation of the central core, we have B0 (z) ∝ z −2 .It should be noted that the approximate function ( 7) at distances z ≫ zcr becomes more than one, which contradicts the condition in the region when the flow of electromagnetic energy almost completely passes into the plasma energy flow.Nevertheless, within the framework of this paper, such distances will not be considered, and this question concerns an approximate function such that the inequality 0 < α (z) < 1 is satisfied.
The advantage of the approach considered here (in contrast to works in which it was necessary to solve a system of two ordinary differential equations) is that the expression (2) allows one to write down the magnetic flux Ψ = 2π Bpϖ dϖ explicitly In turn, it allows us to find all the integrals of motion in any point, as they in standard GS approach (Heyvaerts & Norman 1989;Pelletier & Pudritz 1992) depend on magnetic flux Ψ only.Indeed, the key role in the MHD theory of strongly magnetized jets play the integrals (Beskin 2010) where E(Ψ) (integral Bernoulli) corresponds to the energy flux and L(Ψ) to the angular momentum flux.Here ΩF(Ψ) is the angular velocity which is also constant on magnetic surfaces (Ferraro's isorotation law), and I is the total current inside the magnetic surface.The value Ω0 we introduced earlier is equal to ΩF(0) by definition.Besides, in the black hole magnetosphere both the electric current I and the angular veloc-ity ΩF are determined from the critical conditions on the singular surfaces.Wherein, the current I entering the relations (10)-( 11) turns out to be close to the so-called Goldreich-Julian current (appropriate current density jGJ = ΩFB0/2π), and the angular velocity ΩF is close to half the angular velocity of the black hole rotation ΩH.Thus, with good accuracy, the Bernoulli integral can be written in the form and besides, in the region of strongly magnetized flow, the contribution of particles can be neglected.On the other side, as we know, for both conical (Michel 1973) and parabolic (Blandford 1976) flows near the axis, we can set In what follows, however, we use the integrals (Nokhrina et al. 2019;Chernoglazov et al. 2019) which allow us to include in consideration the reverse electric current flowing within the jet.As a result, knowing the potential Ψ(ϖ, z) and explicit expressions for the integrals of motion ΩF(Ψ), E(Ψ), and L(Ψ), we can determine all components of electric and magnetic fields within the jet Here, however, we add the additional factor (1 + ε) into the expression for the toroidal magnetic field, where ε = ε(Ψ) is a small parameter, ε ≪ 1.This increase of the toroidal magnetic field with respect to exact force-free solution makes it possible to simulate the absence of particle acceleration at large distances, when the electromagnetic energy flux has already been transferred to the plasma flow.Moreover, such a correction does not violate the basic relationship E + (ΩF × r/c) × B = 0. Indeed, using the fundamental theoretical result, according to which, outside the light cylinder, the hydrodynamic velocity V becomes almost equal to the electrical drift velocity (Tchekhovskoy et al. 2008;Beskin 2010) so the hydrodynamic Lorentz-factor of plasma can be written as Γ = Γ dr where Here x = ϖ/RL is the dimensionless distance to the jet axis, and ω = ΩF(x)/Ω0.As a result, at small distances from the central engine we get for ω = 1 i.e., well-known asymptotic solution for collimated magnetized jets (Beskin 1997;Lyubarsky 2009;Tchekhovskoy et al. 2011).As already emphasized, this asymptotic behaviour remains valid in the entire region Bext > Beq.On the other hand, at large distances, i.e. for x > (2ε) −1/2 , we have Γ dr ≈ (2ε) −1/2 ≈ const.This asymptotic behaviour models the saturation region, when the entire energy flow is already concentrated in the hydrodynamic particle flow and Lorentzfactor of the flow reaches its maximal value Γmax ≈ (2ε) −1/2 .As a result, the function ε(Ψ) can be easily found from Bernoulli equation (10) which just determines the maximum Lorentz-factor on a given magnetic surface: Using now relations ( 10)-( 11), we obtain where now due to ( 13)-( 14) and the value of σM (5) has already been defined above.Below we use the velocity field obtained to deternime the Doppler factor.

OBSERVED JET TRANSVERSE DIMENSION
Let us now discuss the observed transverse dimension of a jet djet(z) = 2rjet(z) at an arbitrary distance z from the central engine (see Figure 1).Due to relation (4), the appearance of the central core occurs at the moment when the transverse radius of the jet becomes of order rcr ∼ (γinσM) 1/2 RL.
Accordingly, the distance zcr to the central engine can be evaluated as For reference, we present here the corresponding estimates for parabolic flow as well As for the value rcr, it does not depend on the geometry of the jet.This important result will be used significantly below.
It is clear that as long as the longitudinal magnetic field (and hence the number density) remains uniform within the jet, observed transverse radius of the conical jet rjet(z) increase linearly with distance z as the total transverse radius Here Θjet is the total angular half-width of the conical jet.On the other hand, observed width of the jet djet(z), due to the decrease in the integrals of motion ( 10)-( 11) and the number density towards its edge at large distances z > zcr, is to be smaller than 2Θjetz.Nevertheless, below we use relations ( 24)-( 25) for the estimate of the observed jet width for small distances up to z = zcr.
To more accurately determine the observed transverse size of the jet at any distance from the central engine, we use standard relations (Beresnyak et al. 2003;Pariev et al. 2003;Lyutikov et al. 2003;Lyutikov et al. 2005), allowing to construct the brightness temperature maps.It is based on the hypothesis of the power-law spectrum of radiating particles in the comoving reference frame dN = n ′ γ f (γ ′ ) dγ ′ d 3 r ′ dΩ/4π, when the energy spectrum f (γ ′ ) can be presented in the form Here γ0 < γ ′ < γmax, n ′ γ is the number density of radiating particles in the comoving reference frame, p is the spectrum index, dΩ is the solid angle, and Ke is the normalising constant resulting from the condition f (γ ′ )dγ ′ = 1.
Here two important clarifications must be made.First, below we take into account the evolution of the spectrum of radiating particles due to conservation of the first adiabatic invariant (Beskin et al. 2023a) where is a magnetic field in the comoving reference frame; for magnetically dominated flow considered here one can put h ≈ Bp.
Since the condition p ′ ⊥ ≈ mecγ ′ is satisfied for ultrarelativistic particles, we finally obtain Assuming now that γ ′ max ≫ γ0 and γ0 ≫ 1 we obtain for the normalising constant Ke Ke Thus, when moving away from the central engine, due to the preservation of the first adiabatic invariant (30), the entire spectrum of radiating particles will shift towards low energies without changing its shape.Accordingly the normalising constant Ke will have the form It is easy to check that expression (34) remains the same if the lower limit corresponds to non-relativistic velocities (γ0 = 1).This is due to a decrease in the number of ultrarelativistic particles.
Second, despite the fact that the conservation of the first adiabatic invariant formally leads to a noticeable angular anisotropy of the angular distribution of radiating particles (Beskin et al. 2023a), turbulence and plasma instabilities should lead to isotropization of the distribution function.Therefore, in what follows we do not take into account the effects of anisotropy of the angular distribution of radiating particles.
On the other hand, below we are not limited to the approximation of an optically thin plasma, although this approximation is also valid for sufficiently high frequencies.The point is that we are interested in regions close enough to the central engine, where the effects of synchrotron self-absorption can play a significant role (Frolova et al. 2023).
In what follows, it is convenient to write down the number density of high-energy particles nγ(r) entering in (41) (which now corresponds to the laboratory reference frame) as where is the so-called Goldreich & Julian (1969) number density (i.e., the lowest particle density required to screen the longitudinal electric field), and the constant λγ is the multiplicity factor for radiating particles.The convenience of relation ( 43) is that for highly collimated flows (when the poloidal magnetic field Bp is actually directed along the rotation axis) it can be used at any distance from the central engine.Wherein, the entire dependence of number density nγ(r) on coordinates is contained in Bp(r).
It is clear that the number density of emitting particles λγnGJ should not exceed the number density of the outflowing plasma λnGJ, for which, as estimates show (see, e.g.Mościbrodzka et al. 2011;Nokhrina et al. 2015), λ ∼ 10 11 -10 13 .Therefore, below we assume that the number density of emitting particles is about one percent of the density of the runaway plasma, so that λγ = 10 10 .As a result, we finally get Here h0 is the magnetic field at the light cylinder.Figure 3 shows an example of predicted cross sections of the brightness temperature T br for characteristic jet values (magnetization parameter σM = 20, jet angular half-width Θjet = 10 • , spectrum index p = 2.5) for z = 320 RL (which is larger than zcr ≈ 30 RL) from the central engine assuming the presence and the absence of a dense central core.Everywhere in the calculations it is assumed ν = 15 GHz.As we see, the presence of a central core does result in a significant reduction in the observed transverse size of the jet.Recall that the observed asymmetry is associated with the presence of helical magnetic field (Clausen-Brown et al. 2011) and toroidal plasma velocity, due to which the maximum Doppler factor is achieved at a certain distance from the jet axis (Ogihara et al. 2021;Frolova et al. 2023;Beskin et al. 2023a).Also note that our chosen value of λγ = 10 10 results in reasonable brightness temperature T br .
In addition, Figure 4 shows the dependence of the observed

OBSERVATIONS AND ANALYSIS
According to Section 3, the appearance of the core is accompanied by the break in the jet geometrical profile.Jet in the radio galaxy M87 provides the best opportunity to search for the predicted phenomenon due to its proximity.The mass of the black hole in M87 MBH ≈ 6.5 × 10 9 M⊙ and distance is 16.8 Mpc (EHT Collaboration et al. 2019) provide the angular to spatial scale ratio of 1 milliarcsecond (mas) ≈ 0.08 pc ≈ 260 Rg.Thus, the expected position of the break (24) corresponds to its apparent position from tens of µas for large spins to few mas for moderate spin values a ≈ 0.1.

Searching for the break
For the analysis of the internal jet structure, we considered publicly available 4 15 GHz VLBA observations obtained between 1995 July 28 and 2022 September 29 in frame of the MOJAVE program (Lister et al. 2019), with inclusion of data 4 https://www.cv.nrao.edu/MOJAVE/sourcepages/1228+126.shtml from NRAO archive.Single epoch images often poorly represent large-scale structure and stacking of the single epoch images obtained within sufficiently large time interval could reveal the whole jet cross-section (Pushkarev et al. 2017).Jet of M87 is known to be a hard target for the CLEAN (Högbom 1974) imaging procedure.Pashchenko et al. (2023) showed that CLEAN produces a spurious central brightening in the intrinsically edge-brightened jets that could be suppressed by convolving CLEAN model with a larger beam.Also, M87 parsec-scale jet appearance at 15 GHz is dominated by the helical threads of the Kelvin-Helmholtz instability (Nikonov et al. 2023).These introduce significant problems in the analysis of the jet structure.Thus, we convolved each single epoch CLEAN model with a common circular beam of size 1.4 mas which corresponds to typical size of the naturally weighted beam major axis.This reduces CLEAN imaging artifacts, while decreases the resolution, thus, making our results more conservative.We also filtered out 8 single epochs with a significant image noise to increase the final stacking image fidelity (Pushkarev et al. 2023).
To extract the collimation profile we considered the distribution of the total intensity in the image domain along the slices, taken transverse to the jet direction, which was chosen to be at position angle P A = 17 • with respect to the Right Ascension axis5 .The slices were taken starting from 0.5 mas from the core in steps of 0.1 mas, and for the each cut we fitted two Gaussians to the brightness profile and taking the jet width D as profile width at 10% of the maximal fitted profile value.For the complex resolved asymmetric transverse structure of the M87 jet, where multiple components can have significant brightness difference, the usage of the half maximum level (FWHM) can be misleading by measuring the width of one component instead of the whole jet (Nikonov et al. 2023).As an additional check, we employed simulations with a known parabolic jet model to ensure that usage of 10% level does not bias the geometry estimate.The deconvolved jet cross section therefore was computed as djet = D 2 − θ 2 10% , where θ 10% is the full width of the restoring beam at 10% of the maximum.
In the deriving the jet expansion profile from the images one has to consider the dependence of the nearby pixels.This could affect the uncertainty estimation if the data points are assumed independent and identically distributed or even bias the parameters of the fit.To overcome this we considered a correlated noise model, where the noise is modelled via Gaussian Process (hereafter GP, Rasmussen & Williams 2006).We employed the Rational Quadratic kernel that contains a mixture of scales.
To search the possible jet break, we employed the following relations for the jet width (Kovalev et al. 2020): Here z break is the apparent distance of the break from the jet origin and z0 > 0 corresponds to a separation of the 15 GHz core component from the true jet origin.That is necessary in order to better fit the jet shape near the jet apex (Kovalev et al. 2020).The coefficient a2 is chosen for two power-laws to join each other at z break .We also forced the transition region to be smooth by cubic interpolation model predictions on both sides of the break at distance 0.5 mas.Both single and double slope (i.e. with a break) dependencies were fitted to the data using the Diffusive Nested Sampling algorithm (Brewer et al. 2011) implemented in the DNest4 package (Brewer & Foreman-Mackey 2018) for sampling the posterior distribution of the model parameters.We made use the Bayesian factors for the model selection (Trotta 2008).
Accounting for the correlated noise led to an increase in the uncertainty interval width of the break position by a factor of two, and for the errors of the slope -by four to five times.The resulting fits are presented in Figure 5 where the red colour represents the power-law model ( 46) and the blue colour corresponds to the GP component.The residual GP component significantly deviates from zero, especially in z obs > 15 mas region.This could be a real ef-fect or the result of the limited uv-coverage.The break at z obs = 5.80 +0.65  −0.62 mas is significant with logarithm of the Bayes factor ∆ log Z = 14.9,where Z is the evidence or the marginal likelihood of the model (Trotta 2008).The power-law before the break is consistent with almost conical expansion, while it becomes significantly flatter after the break.Although the break seems to be a robust feature (e.g. it is clearly seen in independent high-sensitivity 15 GHz image of Nikonov et al. 2023, see their Figure 4), the exact values of the exponents of the expansion profile ( 46) should be treated with caution.Pashchenko & Plavin (2019) showed that at least for a relativistic jet model of Blandford & Königl (1979) the jet shape measured from the stacked VLBI images could deviate systematically from the true unobserved value.
We assessed the robustness of our break detection procedure with several tests, including simulations with a known model brightness distribution.For this we created artificial multi-epoch uv-data set using the edge-brightned jet model without the break from Pashchenko et al. (2023) and the same uv-coverages and noise as for the real data from the MOJAVE archive.After imaging the artificial data and stacking single epoch images we employed the same procedure of the break estimation as described above for the real observed data.Below we summarize our tests: • We fitted the synthetic data based on the smooth parabolic jet model with a single power-law and different number of Gaussians (single, double and triple) describing the transverse profile.The resulting fits are generally consistent with the intrinsic parabolic expansion profile k = 0.5 with single Gaussian fit providing the worst estimate k 1G = 0.42 +0.02 −0.02 and double and triple Gaussian fit are consistent with the true value (k 2G = 0.46 +0.03 −0.02 and k 3G = 0.47 +0.01 −0.03 ).• We found that the residual oscillations of the jet width in the synthetic stacked images of the intrinsically smooth model significantly deviate from zero at z obs < 2 mas and z obs > 15 mas.Thus, imaging artifacts could also contribute to the observed oscillations of the M87 jet width as well as helical threads of Kelvin-Helmholtz instability observed in e.g.Nikonov et al. (2023) or oscillations of the expanding magnetized jet (Lyubarsky 2009;Komissarov et al. 2015;Mertens et al. 2016).
• We also fitted the observed transverse profiles with a single and three Gaussians instead of two.The position of the break is stable against different number of Gaussians.The single Gaussian fit results in insignificantly larger break distance z 1G break = 6.24 +0.39 −0.36 mas than a triple Gaussian fit z 3G break = 5.81 +0.56 −0.60 mas.• The fit of the synthetic data based on the smooth jet model revealed a significant -with even larger Bayes factor ∆ log Z = 22.4 than for the real data -break from the conical to the parabolic jet shape with z 2G break = 2.20 +0.50 −0.49mas.Interesting that this type of imaging artefact is opposite to the systematics described in Kovalev et al. (2020) who observed that jet profiles at z obs < 0.5 mas could be artificially flattened.The fits with single and triple Gaussians also showed a break: z 1G break = 3.26 +0.72 −0.58 with ∆ log Z = 21.9 and much less significant z 3G break = 2.24 +0.65 −0.77 mas with modest ∆ log Z = 2.6.• We also checked more the influence of the convolving beam size on the break appearance and employed circular beam size 0.85 mas, that is the typical equivalent area circular beam size from the analysis of Pushkarev et al. (2017).
We conclude that the detected change of the jet shape from conical to parabolic could be affected by the imaging systematics.Thus, the constraints from the break position in the MOJAVE 15 GHz data obtained in the following section should be treated with caution, but could be regarded as a firm limits.

Constrains from the break position
Adopting EHT Collaboration et al. ( 2019) estimation for M87 black hole mass MBH = (6.5 ± 0.9) × 10 9 M⊙ and the distance (16.8 ± 0.8) Mpc, we can use the observed jet width at the observed break position Equation 24 and its uncertainty to estimate the width of the allowed region in (σM, RL) space.As emphasized above, this parameter does not depend on the jet geometry and therefore can be used for both conical and parabolic flows.
Further constrains could be obtained by considering the observed apparent motion at the position of HST-1 (Biretta et al. 1999) and assuming that it corresponds to the maximal Lorentz-factor Γmax in the jet observed at the viewing angle Θ = 17.2 ± 3.3 • (Mertens et al. 2016).Equation B13from Nokhrina et al. (2022) connects Γmax and the initial magnetization σM for the considered MHD model.Using the observed jet width d break = 3.8 ± 0.3 mas at the position of the break in the MOJAVE data estimated in subsection 4.1 and assuming γin ≈ 1, we obtain the constraints presented in Figure 6.Here we also plot the constrained region that corresponds to the end of the conical expansion profile at z obs ≈ 0.6 mas observed in Hada et al. (2016).

DISCUSSION AND CONCLUSION
Thus, it was shown that the formation of a central core near the base of relativistic jets, which was first predicted theoretically and then confirmed by numerical modeling many years ago, should inevitably lead to a change in the dependence of the observed jet width d(z) on the distance z to the central engine.This is due to a sharp decrease in the poloidal magnetic field Bp with distance r from the jet axis at large enough distances z > zcr and, as a consequence, with an rapid decrease in the density of radiating particles ne ∝ Bp (43).As a result, the observed width of the jet djet(z), which at small distances z < zcr can be estimated as djet ≈ 2 rjet, inevitably becomes noticeably smaller than the geometric size of the jet at larger distances.
Thus, the theory predicts another important property of relativistic jets, which in principle makes it possible to carry out additional diagnostics of their physical parameters.The fact that this issue has not yet been discussed in detail up to now is due to the obvious difficulty of studying the internal structure of the jet at short distances z ∼ zcr from the central engine.Note that despite the different literal expressions for conical (25) and parabolic (27) flows, the numerical values of zcr for reasonable parameters (σM ∼ 10-50, γin ≈ 1, and the corresponding dimensionless spin a of M87 SMBH, obtained from the observed radius of the jet at the core appearance position, apparent speeds in HST-1 region and accounting for the uncertainties in the black hole mass, distance and viewing angle estimates.Shade of the color represents 1 and 2 σ intervals.The blue color represents the rcr observed in this paper subsection 4.1.The green color corresponds to the end of the conical expansion profile observed by (Hada et al. 2016, see their fig.9).θjet ∼ 10 • ) differ little from each other.As a result, the internal break should be located at a distance of the order of 10-100 RL, i.e.only 100-1000 gravitational radii.
On the other hand, as was already stressed, some indications of the existence of an internal break have already been noted for the closest source M87, for which the observed position of the break corresponds to a fairly large angular distance ∼ 1 mas.For example, Figure 15 given by Mertens et al. (2016) shows the break in the dependence of the hydrodynamic Lorentz-factors on the distance z to the central machine.Since in the region of magnetically dominated flow we are considering, the well-known connection Γ = r/RL must be satisfied (Beskin 2010), it can also be interpreted as a break in the dependence of the width of the jet djet(z).Accordingly, the break in the dependence of the observed jet width θjet on distance z is clearly visible in Figure 9 presented by Hada et al. (2016) at z obs = 0.6 mas.
We also searched for the predicted break in the jet geometry further downstream the jet using stacked image of 37 single-epoch 15 GHz VLBA observations of M87 radio jet from the MOJAVE archive.The transition from almost conical shape at z obs ≈ 6 mas was found.However, our tests with a synthetic data revealed the presence of the imaging systematical effects that could artificially steepen the expansion profile near the jet origin.
Although the presence of instabilities is generally expected in axisymmetric MHD jets, and the observations of M87 are consistent with development of Kelvin-Helmholtz instability modes (Lobanov et al. 2003;Nikonov et al. 2023), they seem not to disrupt the global jet structure.Porth & Komissarov (2015) have explored the stability of a central core, which we connect to the existence of the inner break.It turned out that we should expect the core stability for the ambient pressure Pext ∝ z b with b > 2, which is consistent with our results for

Figure 1 .
Figure 1.The structure of the magnetic field in the conical model under consideration.At a distances z > zcr from the central engine, when the transverse radius of the jet reaches the scale ∼ rcr, in a conical flow (in which the plasma density and the magnetic field weakly depend on the distance from the axis), a denser central core begins to form.The light cylinder ϖ = c/Ω 0 is shown by a dashed line.

Figure 2 .
Figure 2. Comparison of analytical approximation (1)-(2) (black dashed lines) with the structure of a poloidal magnetic field (blue solid lines) obtained by solving a system of two ordinary differential equations describing internal structure of a cylindrical jet.In this figure, we terminated the solid lines of numerical solution at smaller radii so that the analytical approximation may be better seen.We plot three jet crosscuts for a conical jet with intrinsic opening angle 10 • .

Figure 3 .
Figure 3. Transverse brightness temperature profiles (spectrum index p = 2.5, magnetization parameter σ M = 20, and the jet angular half-width Θ jet = 10 • ) at z = 320 R L from the central engine assuming the presence and the absence of a dense central core.

Figure 4 .
Figure 4. Dependence of the observed jet width d jet (z) at the levels of 10% and 50% of the maximum brightness temperature for two half-widths Θ jet of conical jets.Parabolic asymptotic behaviour d jet ≈ 0.06 z 0.57 pc (Nokhrina et al. 2019) shown by the bold dashed line.Thin dashed lines indicate asymptotic behaviour d jet ∝ z for small z.

Figure 5 .
Figure 5. Jet width versus apparent (projected) distance along the jet z obs fitted by a single (Upper ) and broken (Lower ) powerlaw models.Black points show data points.Lines represents 500 samples from the posterior distribution.Red colour corresponds to the power-law model, while blue colour shows Gaussian Process component.

Figure 6 .
Figure6.Constrains on the light cylinder R L (in units of rg) and the corresponding dimensionless spin a of M87 SMBH, obtained from the observed radius of the jet at the core appearance position, apparent speeds in HST-1 region and accounting for the uncertainties in the black hole mass, distance and viewing angle estimates.Shade of the color represents 1 and 2 σ intervals.The blue color represents the rcr observed in this paper subsection 4.1.The green color corresponds to the end of the conical expansion profile observed by(Hada et al. 2016, see their fig.9).