Formation of Limb-brightened Radio Jets by Angle-dependent Energy Extraction from Rapidly Rotating Black Holes

Using general relativistic magnetohydrodynamic simulations, it has been suggested that the rotational energy of a rapidly rotating black hole (BH) is preferentially extracted along the magnetic field lines threading the event horizon in the middle and lower latitudes. Applying this angle-dependent Poynting flux to the jet downstream, we demonstrate that the jets exhibit limb-brightened structures at various viewing angles, as observed from Mrk 501, M87, and Cyg A between 5° and 75°, and that the limb brightening is enhanced when the jet is collimated strongly. It is also found that the jet width perpendicular to the propagation direction shrinks at the projected distance of the altitude where the jet collimates from a conical shape (near the BH) to a parabolic one (in the jet). Comparing with Very Long Baseline Interferometry observations, we show that this collimation takes place within the deprojected altitude of 100 Schwarzschild radii from the BH in the case of the M87 jet.


INTRODUCTION
In many active galactic nuclei, collimated jets are launched from their central regions, where supermassive black holes (BHs) presumably exist.It is widely believed that these relativistic flows are energized either by the extraction of the BH's rotational energy via the magnetic fields threading the event horizon (Blandford & Znajek 1977) or by the extraction of the accretion flow's rotational en-Corresponding author: Kouichi Hirotani, Hsien Shang hirotani,shang@asiaa.sinica.eduergy via the magnetic fields threading the accretion disk (Blandford & Payne 1982).In the recent two decades, general relativistic (GR) magnetohydrodynamic (MHD) simulations successfully demonstrated that a Poynting-dominated jet can be driven by the former process, the socalled "Blandford-Znajek (BZ) process" (Koide et al. 2002;McKinney & Gammie 2004;Komissarov 2005;Tchekhovskoy et al. 2011;Qian et al. 2018).The BZ process takes place when the horizon-penetrating magnetic field (B) is dragged in the rotational direction of the BH by space-time dragging, and when there is a merid-ional current (J ) flowing just above the horizon (as a part of the global pattern of magnetospheric currents).The resultant J × B Lorentz force exerts a counter torque on the horizon, extracting the BH's rotational energy electromagnetically, and carrying the energy and momentum to large distances in the form of torsional Alfven waves when the toroidally twisted magnetic field lines unwind.
Subsequently, applying the GR particle-incell (PIC) technique to magnetically dominated BH magnetospheres, it is confirmed that the BZ process does work even when the Ohm's law, and hence the MHD approximation breaks down (Parfrey et al. 2019;Chen & Yuan 2020;Kisaka et al. 2020;Crinquand et al. 2021;Bransgrove et al. 2021;Hirotani et al. 2023).However, with current computational capability, the GRPIC technique cannot be applied to supermassive BHs, if we consider the magnetic fields whose strengths are comparable to what inferred from observations, because the plasma skin depth and/or the gyration motion cannot be easily resolved in such cases.As a result, GR-PIC codes have been applied to 2D (axisymmetric) magnetospheres around either a supermassive BH with an extremely weak magnetic-field strength, or a stellar-mass BH with a realistic magnetic-field strength.Nevertheless, it is confirmed that the BH's rotational energy can be efficiently extracted via the BZ procss by both GRMHD and GRPIC techniques.
It is possible that such extracted electromagnetic energy is converted into the kinetic and internal energies of electron-positron pairs at some distance from the BH, and eventually dissipated as radiation via the synchrotron and inverse-Compton processes in the jet downstream.In the present paper, we thus semianalytically examine the jet emission properties, connecting the jet-launching and jetdownstream regions.
In § 2, we apply published GRMHD simulation results to the jet-launching region near the BH, and constrain the angular dependence of the BZ flux.We then describe how such electromagnetic energies are converted into the plasma's energies in the jet in § 3, utilizing another published MHD simulation results on jet propagation.Then in § 4, we depict the expected VLBI maps from low-accretion BH systems and demonstrate that the limb-brightened structures, which has been observed in the jets of Mar 501 (Giroletti et al. 2004), M87 (Kovalev et al. 2007), 3C84 (Nagai et al. 2014), and Cyg A (Boccardi et al. 2016), can be formed by virtue of the angle-dependent BZ flux.Finally in § 5, we compare the present analysis with an analytical work done by Takahashi et al. (2018).

POYNTING FLUX NEAR THE BLACK HOLE
Let us start with constraining the meridional dependence of the BZ flux in the BH vicinity.In a stationary and axisymmetric BH magnetosphere, the radial component of the BZ flux (or equivalently the Poynting flux) is given by where Ω F ≡ F tθ /F θϕ designates the angular frequency of the rotating magnetic field, B µ ≡ * F µ t does the µ component of the magnetic field; F and * F denote the Faraday and Maxwell tensors, respectively.The strength of the toroidal component of the magnetic field is expressed up to a numerical factor by the enclosed poloidal current I(A φ ) by where A φ denotes the magnetic flux function.In a stationary and axisymmetric magnetosphere, a magnetic field line resides in a constant A φ surface.
In this paper, we model the jet flow lines in the poloidal plane with a parabolic-like geometry (Broderick & Loeb 2009;Takahashi et al. 2018), where r H = M + √ M 2 − a 2 denotes the horizon radius in geometrized unit (i.e., c = G = 1), and θ does the colatitude; c and G refer to the speed of light and the gravitational constant, respectively.Equation (3) gives a pure parabolic flowline if q = 1 and a conical one if q = 0. We adopt this flow-line model, because the M87 jet is known to be parabolic-like up to the Bondi radius (Asada & Nakamura 2012).
Differentiating equation ( 3) with respect to θ, we obtain the radial component of the magnetic field, (4) where a denotes the BH's spin parameter, ∆ ≡ r 2 − 2M r + a 2 vanishes as the horizon, and Σ ≡ r 2 + a 2 cos 2 θ.At the horizon, r = r H , We thus obtain the BZ flux, where ω H ≡ a/(2M r H ) denotes the angular frequency of the spinning BH.
It is reasonable to assume that the magnetic field lines are predominantly radial in the BH vicinity even in a magnetically dominated magnetosphere (in which plasma's rest-mass energy density is negligible compared to that of the magnetic field), because of the causality at the horizon and the particle's inertia (e.g., Hirotani et al. 1992;McKinney et al. 2012;Bransgrove et al. 2021).Then Tchekhovskoy et al. (2010) performed a GRMHD simulation and demonstrated that Ω F ≈ 0.5ω H and I(A φ )/(A max ω H ) ≈ 3 sin((π/2)(A φ /A max ) hold for a radial field line.We thus obtain where S r 0,H ≡ a 2 ω 2 H A 2 max /(8πΣ 2 ) denotes the typical, angle-averaged BZ flux at the horizon, and w ≡ Ω F /(0.5ω H ).An additional factor f = f (θ) is of an order of unity, and is set to be f = 1 in most cases.However, to test the effect of the vanishment of B φ on the equator, we will consider the case of a vanishing f (θ) at θ = π/2 in § 4.2.Nevertheless, the main conclusion of this paper, that is, the formation of a limb-brightened jet, is little affected by the deviation of f from unity.
We plot S r H /S r 0,H as a function of θ in figure 1.The red solid curve shows the BZ flux when f = 1.The black dashed one does the BZ flux when f = x 2 /(x 2 + 1), where x = 5(θ − π/2).In both cases, the BZ flux is proportional to θ 4 at |θ| ≪ 1.It is, therefore, difficult to contrive a Poynting-dominated jet along θ = 0 axis within the force-free or MHD framework.We will use the red solid case as the main scenario in § 4, and use the black dashed case only for comparison in § 4.2

JET MODEL
In this paper, we assume that the jet is composed of electron-positron pair plasmas, and consider both thermal and non-thermal components.

Kinetic flux inferred by the BZ flux
When a jet is composed of a pure pair plasma, its kinetic flux can be computed by where n th * and n nt * denote the co-moving number density of thermal and non-thermal electrons, respectively; βc refers to the fluid velocity, Γ ≡ 1/ 1 − β 2 the bulk Lorentz factor, Θ ≡ kT e /m e c 2 the dimensionless temperature of thermal electrons, ⟨γ⟩ the mean Lorentz factor of randomly moving non-thermal electrons, and m e c 2 the electron's rest-mass energy.
Neglecting the energy dissipation in the jet (Celotti & Fabian 1993), we find that the summation of electromagnetic and kinetic energies is conserved along each flux tube.We thus obtain the kinetic-energy flux at position (r,θ) in the jet, where σ denotes the magnetization parameter, B p (r, θ) the poloidal magnetic field strength at distance r from the BH and at colatitude θ in the jet.F H /B p (r H , θ).In addition, B p (r, θ) can be computed by equation (4), and S r H /S r 0,H is known ( § 2).We constrain the normalization factor S r 0,H ∝ A 2 max by adjusting the integrated flux density (at 86 GHz in the preset paper) with the observed value.Therefore, the right-hand side of equation ( 8) is obtained, once σ = σ(r, θ) is specified.

Emission and absorption coefficients
Let us consider how to constrain the comoving densities, n th * and n nt * , combining equations ( 7) and ( 8).For thermal pairs, the emission coefficient is given by Leung et al. (2011), and the absorption coefficient can be computed by the Kirchhoff's law.Their distribution function obeys the Maxwell-Jüttner distribution, which is characterized by Θ.For non-thermal pairs, emission and absorption coefficients are given in e.g., ?Le Roux (1961); Ginzburg & Syrovatskii (1965).We present the explicit expressions of these coefficients in the Appendix.
To compute the emission and absorption coefficients, it is essential to constrain both n th * and n nt * .Let us parameterize the former by n th * = ζn nt * , which gives the total proper density n * = (1+ζ)n nt * .Ito et al. (2008) investigated the kinetic powers and ages of four Fanaroff-Riley II radio galaxies, and derived 4 < 1 + ζ < 310 under a minimum-energy condition (e.g., Miley 1980;Kellermann & Pauliny-Toth 1981).

Assumptions
For thermal pairs, we assume that they have a semi-relativistic temperature, Θ = 1.0, and dominate the non-thermal ones by a factor ζ = (1 − 0.007)/0.007≈ 142; that is, 0.7 % of the pairs are non-thermal.It is noteworthy that ζ merely changes the normalization of the jet luminosity.Thus, when ζ has a greater value, for instance, we can obtain a very similar jet solution by adopting an appropriately smaller value of A max .In another word, the main conclusions of the present paper -limb brightening and constriction of jets -are little affected by the choice of such parameters.
For the non-thermal pairs, we set the powerlaw index p = 3.0, which gives the spectral index of α = −1.0 in the optically thin frequency regime.The lower and upper cut-off Lorentz factors are γ min = 1 and γ max = 10 5 .
We assume an energy equipartition between the energy density of (thermal+non-thermal) pairs and that of the random magnetic field, B 2 * /8π.
For the jet geometry, we model the jet flow lines so that they may coincide with the magnetic-flux surface of a parabolic-like geometry (eq.[3]).In a super-fast-magnetosonic region, magnetic field lines are substantially wound up in the trailing direction (Camenzind 1986); accordingly, the motion of an MHD flow becomes predominantly poloidal far outside the outer light surface (Takahashi et al. 1990;Toma & Takahara 2013).Thus, we neglect the rotational motion of fluids around the jet axis.
In this paper, we assume that q (eq.[3]) increases linearly with r inside the collimation altitude r = r c , and parameterise it by Accordingly, the flow lines change from a conical shape (i.e., q = 0) near the BH at r − r H ≪ r c to a parabolic one (i.e., q = q ∞ ) beyond the altitude r = r c .We adopt an observational value, q ∞ = 0.75 (Nakamura et al. 2018).
As an example, we apply the present method to the M87 jet, and adopt a BH mass 6.4 × 10 9 M ⊙ (Gebhardt et al. 2011;Event Horizon Telescope Collaboration et al. 2019;Liepold et al. 2023) and the luminosity distance 16.7 Mpc.An angular scale of 1 mas corresponds to a spatial scale 140R S , where R S ≡ 2GM c −2 = 2M denotes the Schwarzschild radius.We assume a rapidly rotating BH and adopt a = 0.9M (Li et al. 2009;Dokuchaev 2023;Feng & Wu 2017).

Constraining the jet kinematics
Analyzing multi-epoch VLBI data with the wavelet-based image segmentation and evaluation method, Mertens et al. (2016) investigated the kinematics of the M87 jet.They revealed a stratification of the jet flow with a slow, subluminal and fast, superluminal (i.e., relativistic) velocity components.Since the fast velocity component is present in all the three emission ridges, they argue that it can be associated with the bulk fluid motion in the sheath, rather than the spine of the jet.On the other hand, the slow component may be interpreted as a pattern speed associated with an instability.
On these grounds, we consider that the fast, relativistic velocity component represents the plasma motion in the jet.

Evolution of the bulk Lorentz factor
Using the observed apparent velocity as a function of distance, Mertens et al. (2016) revealed that the bulk Lorentz factor increases from Γ ≈ 2 at distance ϖ ≈ 0.4 mas from the jet axis (i.e., at de-projected distance r ≈ 1 mas from the BH) to Γ ≈ 4.5 at ϖ ≈ 1.5 mas (i.e., at r ≈ 20 mas).It should be noted that r in Mertens et al. (2016) designates the distance ϖ from the jet axis in the present notation.We thus approximate the Lorentz factor evolution by in the present paper.At larger scales, r ≫ 1 mas, the bulk Lorentz factor is found to be Γ ∼ 6 by HST observations (Biretta et al. 1999).We thus set an upper limit of 6 in the right-hand side.Nevertheless, the saturation of Γ at 6 does not play a role in the present analysis, because we consider only the inner jet, r ≪ 100 mas, in the de-projected distance from the BH.

Evolution of the magnetization parameter
In relativistic MHD jet models, most of the conversion of energy from Poynting to kinetic takes place when highly compressed magnetic field lines (in the upstream) expands in the super-magnetosonic region (Li et al. 1992;Chiueh et al. 1998;Vlahakis & Königl 2003).Accordingly, the increased plasma's mass winds the field lines into the counter-rotational direction to collimate the field lines into a parabolic geometry due to the hoop stress (Tomimatsu & Takahashi 2003).The fast-magnetosonic point is located typically at several light cylinder radii, ϖ LC (Li et al. 1992;Chiueh et al. 1998), whereas ϖ LC typically becomes a few Schwarzschild radii (Znajek 1977;Takahashi et al. 1990).For instance, assuming a quasiparabolic jet flow line, and applying a cold, ideal MHD model to the M87 jet, Mertens et al. (2016) finds that the magnetization parameter σ, the ratio between the Poynting and kinetic energy fluxes, is much greater than 1.0 near the BH, but decreases with distance from the BH to become 2.0 at fast point, whose projected distance is z ≡ r sin θ v ≈ 1 mas along the jet from the BH, where θ v denotes the observer's viewing angle of the jet.In what follows, z designates the ordinate of the expected VLBI maps (figs.2, 3 7, and 7), whereas ϖ does the abscissa of figures 2-6.
On these grounds, we assume that σ evolves with the de-projected distance r from the BH such that which mimics the MHD simulation of Mertens et al. (2016).Outside the fast point, σ continuously decreases and attain σ ≪ 1, which means that the jet becomes kinetic-dominated asymptotically (e.g., Chiueh et al. 1998).The terminal Lorentz factor tends to σ 0 1/3 , where σ 0 (≫ 1) denotes the σ at the jet launching point.
It is noteworthy that an efficient gamma-ray emission, as observed during the VHE flares in year 2008 and 2010 (Acciari et al. 2009;Aliu et al. 2012;Abramowski et al. 2012), takes place as a result of synchrotron self-Compton process only when the jet becomes particle-dominated (Blandford & Levinson 1995).In addition, the VHE gamma-rays are appeared to be emitted within 20R S from the black hole (Hada et al. 2013).In this context, the rapid conversion of energy from Poynting to kineitic within the central mas scales (eq.[11]) is consistent with the VHE observations.

FORMATION OF LIMB-BRIGHTENED JETS
In this section, we examine if a jet appears limb-brightened in radio frequencies.We adopt the photon frequency ν = 86 GHz in our frame of reference, and compute expected VLBI maps.

Jet collimation altitude
We follow the method described in § 3 to constrain the emission and absorption coefficients in the jet co-moving frame ( § 3.2).Then we apply their Lorentz invariant relations to convert these quantities in our frame of reference.Using the obtained emission and absorption coefficients, we integrate the radiative transfer equation to compute the specific intensity I ν along each line of sight.Finally, the surface brightness is obtained by multiplying I ν with the solid angle subtended by each area on the map.
In this subsection ( § 4.1), we constrain the altitude r c below which the jet is collimated, comparing the results with VLBI observations.For this purpose, we fix the observer's viewing angle at θ v = 0.30 rad (i.e., 17 • ), thermal lepton's temperature at Θ = kT e /m e c 2 = 1.0, and the non-thermal lepton's power-law index at p = 3.0.
In figure 2, we compare how the jet appears constricted (perpendicularly to the propagation direction) when a collimation (from conical to parabolic-like flow line) takes place at various de-projected altitude r = r c from the BH.In panels (a), (b), and (c), we adopt r c = 50R S , r c = 100R S , and r c = 200R S , respectively.It follows that the jet width perpendicular to the propagation direction shrinks.This is because the jet is assumed to change its shape from conical to parabolic one inside the de-projected distance r c (eq. 9), and because this change of flowline geometry leads to a shrink of transverse width in the celestial plane.
The constricted structure appears at z = 0.10 mas, 0.22 mas, and 0.48 mas when r c /R S = 50, 100, and 200, respectively.On the other hand, VLBI observations of M87 at 86 GHz reveals that a constricted structure appears at z = 0.2-0.3mas from the VLBI core (Hada et al. 2016).Thus, we consider that the M87 jet finishes its collimation within a de-projected altitude r ≈ 100R S .

Limb-brightening of jets
Figure 2 also shows that the jet exhibits a limb-brightened structure for various values of the collimation height r c .This is because the BH's rotational energy is preferentially extracted along the magnetic field lines threading the event horizon in the middle or lower latitudes (fig.1).In another word, the energy carried along each magnetic flux tube increases with the distance ϖ from the jet axis, resulting in an efficient mass loading and emission along the field lines near the jet limb.In this paper, we adjust the normalization constant A max so that the integrated flux density of figure 2b may be consistent with observed values (∼ 1 Jy at 86 GHz) (e.g., Kim et al. 2018;Hada et al. 2016).We accordingly obtain B p = 30 G.
Let us next consider how the limb-brightened jet changes its appearance as a function of the observer's viewing angle, θ v .In figure 3, we adopt θ v = 0.10 rad, 0.60 rad, and 1.00 rad, as denoted in each panel, keeping other parameter unchanged from figure 2b.Comparing figure 3 with figure 2b, we find that the jet increases its brightness with decreasing θ v as a result of the relativistic beaming effect.It also follows that the jet perpendicular thickness decreases with increasing θ v as a result of a projection effect.Namely, the apparent half opening angle, χ app = atan(tan χ int / sin θ v ), increases with decreasing θ v for a fixed intrinsic half opening angle, χ int .

Dependence on the magnetic-field geometry
Figure 2. Surface brightness distribution of synchrotron emission from the M87 jet at 86 GHz.For all the panels, we set θ v = 0.30 rad, Θ ≡ kT e /m e c 2 = 1.0, and p = 3.0.In panels (a), (b), and (c), we set r c = 50R S , 100R S , and 200R S , respectively.Color coding is common to all the panels and defined in the color bar on the right.The contour level increases by 10 0.5 toward the VLBI core, which is optically thick for synchrotron self-absorption.We put f = 1 in equation ( 6) in all panels.
Let us examine how the limb-brightened structure depends on the large-scale magnetic field configurations.

Dependence on poloidal field geometry
We investigate how the jet structure depends on the poloidal componenents of the magnetic field.First, we consider its dependence on the collimation altitude, r c .In figure 4, we depict the brightness profile (or equivalently, the transverse intensity profile) as a function of ϖ (i.e., the distance across the jet axis), adopt-ing four discrete collimation altitudes, r c = 100R S , 200R S , 300R S , and 400R S .The black solid, red dashed, green dash-dotted, blue dotted, and cyan dash-dot-dot-dotted curves show the transverse intensity at projected distance z = 1.0 mas, 2.0 mas, 4.0 mas, 8.0 mas, and 16.0 mas from the BH.We fix other parameters than r c , adopting θ v = 0.30 rad, Θ ≡ kT e /m e c 2 = 1.0, p = 3.0, and q ∞ = 0.75.It follows from the top two panels of figure 4 that the transverse intensity profile does not change as long as r c < 200R S .From this reason, we con- 2, but the viewing angle is changed.We set θ v = 0.1, 0.6 and 1.0 for panels (a), (b) and (c), respectively.Other parameters are commonly set as r c = 100R S , Θ = 1.0, and p = 3.0.Color coding is common to all the panels and defined in the color bar on the right.We put f = 1 for all panels.All panels should be compared with fig.2b.strain r c , using the observed constricted structure of the M87 jet (Hada et al. 2016), instead of the observed transverse intensity profile (e.g., Kim et al. 2018).It follows from the bottom two panels of figure 4 that the limb-brightening is enhanced when r c > 300R S .However, if we adopt such a large (de-projected) collimation altitude, a constricted structure appears at a much greater projected distance z than observed (Hada et al. 2016).Thus, we do not consider such a large r c in what follows.
Second, we consider its dependence on the collimation degree q (eq.9).In figure 5, we depict the transverse intensity profile, adopting two representative values of q ∞ .The left panel corresponds to the result when the jet geometry tends to a quasi-parabolic one with q ∞ = 0.75 (Asada et al. 2016;Nakamura et al. 2018), whereas the right panel shows the case of a pure-parabolic geometry (i.e., q ∞ = 1).As expected, jet perpendicular width shrinks when the jet is collimated more strongly (i.e., when q is greater).
It is worth comparing the obtained transverse intensity profile with observations.In figure 6, we plot the results at projected distance z = 0.6 mas (left) and z = 0.8 mas (right), for q ∞ = 0.5 (blue dash-dotted), q ∞ = 0.75 (red solid), and q ∞ = 1.0 (black dashed).The thin black dotted curves in each panel denote what obtained from the M87 jet at 86 GHz (Kim et al. 2018).It shows that the limb brightening is more enhanced when q ∞ increases, or equivalently, when the jet is more strongly collimated.We also find that the value q ∞ = 0.75, which was inferred by observations at much larger scales, gives a consistent width of the brightened limbs.However, the peak intensity of the two limbs attains only half of what is observed.
We also plot the expected VLBI map in figure 7. Panels (a) and (b) represent the results for q ∞ = 0.5 and q ∞ = 1.0, respectively, and should be compared with the q ∞ = 0.75 case (fig.2b).The limb-brightening structure, as well as the constricted structure, commonly appears for a wide range of q.

Dependence on toroidal field geometry
We investigate how the jet structure depends on the toroidal componenents of the magnetic field.Particularly, we are interested in the case when B φ , and hence the BZ flux vanishes on the equator.Thus, to constrain the BZ flux at the jet base, we adopt the black dashed curve, instead of the red solid one, in fig. 1, or equivalently, we suppress f (θ) near the equator as described at the end of § 2. We present the resultant VLBI map in figure 7c.It shows that the jet appearance changes only mildly when B φ vanishes near the equator.
On these grounds, we consider that main conclusions of the present paper -limb brightening and constriction of a jet -are robust for a reasonable change of the large-scale magneticfield geometry.

Dependence on lepton energy distribution
We next consider how the surface brightness distribution depends on the energy distribution of thermal and non-thermal leptons.In fig-ure 8a, we present the result when the thermal leptons have higher temperature, kT e = 3.0m e c 2 .Comparing with figure 2b, we find that the limb-brightening feature is not affected by the lepton temperature.The difference of temperature appears only at the jet base, within the central a few mas.
Let us also examine the difference when the non-thermal lepton's power-law index p varies.In figures 8b and 8c, we show the cases when p = 2.5 and p = 3.5.Comaparing with figure 2b (i.e., the case of p = 3.0), we find that the limb-brightened structure is not affected by the changle of p, as long as it is between 2.5 and 3.5.

Counter jet
The dim structure appearing in z < 0 in figures 2, 3, 7 and 8 shows the counterjet.The counter jet also exhibits a limbbrightened structure.Because of the relativistic beaming, the jet is brighter than the counterjet.Let us consider the fiducial case of r c = 100R S , Θ = 1.0, p = 3.0, and θ v = 0.30 rad (i.e., fig.2b).For the jet, the peak brightness attains 35.2Jy mas −2 , 32.4Jy mas −2 , 17.5Jy mas −2 , and 9.4Jy mas −2 at θ v = 0.10, 0.30, 0.60, and 1.00, respectively.For the counterjet, it attains 0.57Jy mas −2 , 0.75Jy mas −2 , 0.88Jy mas −2 , 1.4Jy mas −2 , at the same θ v 's.The integrated flux density within the central 2 and 5 mas becomes 0.11 Jy and 0.23 Jy for the jet, whereas it is 2.2 mJy and 2.7 mJy for the coutnerjet.Accordingly, we find that the jet is brighter than the counterjet by a factor R, which lays between 50 and 85 in the core region.

DISCUSSION
To constrain the angle-dependent energy extraction from a rotating BH, we use the existing MHD simulations of BH magnetospheres (Tchekhovskoy et al. 2010)  kinetic energy flux, we use another MHD simulation of Mertens et al. (2016) for the M87 jet.Then we compute the lepton density by the method described in § 2, and demonstrate that the M87 jet naturally exhibits a limb-brightened structure as a result of this angle-dependent energy extraction from the BH.We also showed that the jet exhibits a constricted structure because of a collimation from conical to paraboliclike shape.

Jet versus counter jet brightness ratio
In § 4.5, we find that the jet-counterjet brightness ratio becomes 50 < R < 85 in the core region.However, these values are a few times greater than the value R = 27.1 ± 9.1 obtained within the central 0.2-0.5 mas in the M87 jet ( § 3.6 of Kim et al. 2018).We could in fact, obtain a consistent R with such VLBI observations, if we adopted a slower evolution of Γ or σ (than eqs. 10 or 11) as a function of r.However, such a fine-tuning of functions or parameters are out of the scope of the present paper, because the main conclusions -the limb brightening and constriction of a jet -are not affected by such details.

Formation of kinetic-dominated jets
Let us briefly comment on the evolution of the magnetization parameter σ as described by (pure-parabolic) Figure 5. Transverse jet intensity profiles for q ∞ = 0.75 (left) and q ∞ = 1.00 (right).The black solid, red dashed, green dash-dotted, blue dotted, and cyan dash-dot-dot-dotted curves show the brightness at altitude z = 1, 2, 4, 8, and 16 mas in the projected distance from the BH along the jet.The abscissa denotes the distance (mas) measured perpendicular to the jet.
equation ( 11).We adopt this equation as an assumption, because it mimics the 2D, cold, ideal MHD simulation of Mertens et al. (2016), specifically their models 1 and 2, which best reproduce the measured acceleration of jet components (their § 4).However, to the best of the authors' knowledge, no 3D MHD simulations that selfconsistently solve the system from the BH vicinity to the jet downstream, have successfully demonstrated such a rapid conversion of the electromagnetic energy into the kinetic one, e.g., within the initial 100R S in the case of the M87 jet.On the other hand, it is an observational fact that the jet becomes already kineticdominated (Celotti & Fabian 1993;Blandford & Levinson 1995) in the high-energy emission re-gion, which is possibly located within the inner 20R S for the M87 jet (Hada et al. 2013).
To settle down this energy conversion argument, namely how such a rapid energy conversion could take place in the jet-launching region, it may be necessary to quantitatively investigate the pair-production cascade in a jet (Levinson & Blandford 1995, 1996), and the resultant massloading issue.

Comparison with Takahashi et al. (2018)
To demonstrate the formation of a limbbrightened jet, we used the angle-dependent energy extraction from the BH, and constrained the plasma density at each position in the jet.On the other hand, Takahashi et al. (2018) assumed that the density peaks within a ringlike geometry away from the jet axis.In the present paper, we obtained a symmetric bright-Distance across the jet (mas) Surface brightness ( Jy / mas 2 ) r = 0.6 mas r = 0.8 mas Distance across the jet (mas) Figure 6.
Transverse jet intensity profiles at altitude z = 0.6 mas (left) and 0.8 mas (right).The blue dash-dotted, red solid, black dashed curves denote the brightness computed by the present model for q ∞ = 0.50 (weakly parabolic), q ∞ = 0.75 (quasi parabolic), and q ∞ = 1.00 (purely parabolic), respectively.The thin dotted curve does the observed values reported by Kim et al. (2018).
ness distribution with respect to the jet axis, because we neglect the rotational motion of fluids around the jet axis.This assumption, a ballistic motion of jet fluids, will be justified if the jet is kinetic-dominated.Since σ ∼ 1 is considered to be achieved at the projected distance of z = 1.5 mas (Mertens et al. 2016), the jet will become kinetic-dominated at z ≫ 2 mas.On the other hand, Takahashi et al. (2018) took the fluid's rotational motion into account in the intermediate region where the electromagnetic energy is being converted into the kinetic energy, and showed the formation of asymmetric jets due to relativistic beaming of radiation.

Stronger emissivity along the limb
As discussed in § 4.3.1,we failed to obtain a strong enough intensity at the limbs, comparing with what are obtained in the M87 jet (fig.6).
This might be because we adopted a stationary, ideal MHD model in the present paper.Nevertheless, it is likely that the jet carrying strong energy flux near the limb (fig. 1) causes fluid instability or magnetic reconnection at the layer with the disk wind, and accelerate nonthermal particles at the jet limb (Matsumoto & Masada 2013;Parfrey et al. 2015).Alternatively, pair-production cascade taking place in a BH magnetosphere (Beskin et al. 1992;Hirotani & Okamoto 1998;Neronov & Aharonian 2007;Rieger & Aharonian 2008;Levinson & Rieger 2011;Broderick & Tchekhovskoy 2015;Hirotani & Pu 2016) could result in a angledependent mass loading in the jet launching region.This effect works most efficiently when the mass accretion rate is so small that the gap longtitudinal (radial) thickness becomes compa- rable to the transverse width in the meridional direction (Hirotani et al. 2016).In this case, the voltage drop along magnetic field lines becomes a good fraction of the electro-motive force exerted on the BH surface.Accordingly, the sum of the kinetic energy density of the accelerated charges and the radiation energy density of the emitted photons becomes a good fraction the typical electromagnetic energy density near the BH.If it happens along the magnetic field lines threading horizon in the middle latitude, stronger emission may be expected along such flow lines, resulting in a stronger emission from the limbs.

Formation of the central-ridge emission
VLBI observations have revealed that the M87 jet has a triple-ridge structure from the jet base whose distance from the central ring structure is approximately 30R S (Lu et al. 2023) to larger scales whose distance is greater than 100R S (Kim et al. 2018).Near the jet base, the in- tensity of the the central-ridge structure (near the jet axis) attains about 60 % of that of the outer-ridge (i.e., the limb).At larger scales, it attains about 45 % of that of the outer ridge.
On these grounds, as long as the M87 jet is concerned, it is not enough to discuss only the two brightened limbs.Nevertheless, we consider that the central ridge emission is formed when the polar funnel of the BH magnetosphere is failed to be described within the MHD framework, which we have based on in the present paper.For instance, a strong BH gap is formed in the polar region when the accretion rate is much small compared to the Eddignton rate, and when the BH is extremely rotating (Song et al. 2017).Therefore, taking account of such non-force-free or non-MHD effects, we will investigate if the emission is enhanced along the jet axis in our subsequent papers.

Figure 1 .
Figure 1.Stationary and axisymmetric BZ flux as a function of the colatitude.The red solid curve denotes the BZ flux when f = 1 (see text), whereas the black dashed one does the BZ flux when f , and hence B φ vanishes on the equator.
the B p at the same point.At the jet base, we obtain

Figure 3 .
Figure 3. Similar figure as fig.2, but the viewing angle is changed.We set θ v = 0.1, 0.6 and 1.0 for panels (a), (b) and (c), respectively.Other parameters are commonly set as r c = 100R S , Θ = 1.0, and p = 3.0.Color coding is common to all the panels and defined in the color bar on the right.We put f = 1 for all panels.All panels should be compared with fig.2b.
. To quantify how such a Poynting flux is converted into the jet's Distance across the jet (mas) Surface brightness

Figure 4 .
Figure 4. Transverse jet intensity profiles for r c = 100R S (top left), 200R S (top right), 300R S (bottom left), and 400R S (bottom right).The black solid, red dashed, green dash-dotted, blue dotted, and cyan dash-dot-dot-dotted curves denote the brightness at z = 1, 2, 4, 8, and 16 mas from BH.No difference can be seen between the two top panels.

Figure 7 .
Figure 7. Similar figure as fig.2, but the global magnetic field configuration is changed.Panels (a) and (b): We change the poloidal field configuration, adopting q ∞ = 0.5 (weakly parabolic) and q ∞ = 1.0 (pure parabolic) flow lines in panel (a) and (b), respectively.Panels (c): we change the toroidal field configuration, suppressing f below 1 near the equator (black dashed curve in fig.1; see also the end of § 2).All panels should be compared with fig.2b.In the present paper, panel (c) of this figure shows the only case in which f deviates from unity in equation (6).

Figure 8 .
Figure 8. Similar figure as fig.2, but the electron energy distribution is changed.Panel (a): we consider higher temperature for thermal leptons, and put Θ = 3.0, while other parameters are common with fig 2a.Panels (b) and (c): we consider different power-law index for nonthermal leptons, and put p = 2.5 in panel (b) and p = 3.5 in panel (c).Color coding is common.We put f = 1 in all panels, which should be compared with fig.2b.