Matter density distribution of general relativistic highly magnetized jets driven by black holes
DDraft version February 17, 2021
Typeset using L A TEX twocolumn style in AASTeX63
Matter density distribution of general relativistic highly magnetized jets driven by black holes
Taiki Ogihara, Takumi Ogawa,
2, 1 and Kenji Toma
3, 1 Astronomical Institute, Graduate School of Science, Tohoku University, Sendai, Miyagi, 980-8578, Japan Center for Computational Sciences, University of Tsukuba, Tsukuba, Ibaraki, 305-8577, Japan Frontier Research Institute for Interdisciplinary Sciences, Tohoku University, Sendai, Miyagi, 980-8578, Japan
ABSTRACTHigh-resolution very long baseline interferometry (VLBI) radio observations have resolved the de-tailed emission structures of active galactic nucleus jets. General relativistic magnetohydrodynamic(GRMHD) simulations have improved the understanding of jet production physics, although theo-retical studies still have difficulties in constraining the origin and distribution of jetted matter. Weconstruct a new steady, axisymmetric GRMHD jet model to obtain approximate solutions of black hole(BH) magnetospheres, and examine the matter density distribution of jets. By assuming fixed poloidalmagnetic field shapes that mimic force-free analytic solutions and GRMHD simulation results andassuming constant poloidal velocity at the separation surface, which divides the inflow and outflow, wenumerically solve the force-balance between the field lines at the separation surface and analyticallysolve the distributions of matter velocity and density along the field lines. We find that the densitiesat the separation surface in our parabolic field models roughly follow ∝ r − in the far zone from theBH, where r ss is the radius of the separation surface. When the BH spin is larger or the velocity at theseparation surface is smaller, the density at the separation surface becomes concentrated more near thejet edge. Our semi-analytic model, combined with radiative transfer calculations, may help interpretthe high-resolution VLBI observations and understand the origin of jetted matter. Keywords:
Active galactic nuclei (16), Black hole physics (159), General relativity (641), Relativisticjets (1390), Magnetohydrodynamics (1964) INTRODUCTIONRelativistic collimated outflows (or jets) are observedin active galactic nuclei (AGNs). They originate fromthe system composed of a supermassive black hole (BH),an accretion disk, and a surrounding gas at the centerof the galaxy. Very long baseline interferometry (VLBI)radio observations have revealed detailed emission struc-tures of AGN jets. Limb-brightened structures havebeen observed at mm-cm wavelengths in jets of M87(Kovalev et al. 2007; Walker et al. 2008; Hada et al. 2011,2016; Mertens et al. 2016; Kim et al. 2018; Walker et al.2018), Mrk 501 (Piner et al. 2008), Mrk 421 (Piner et al.2010), Cyg A (Boccardi et al. 2016), and 3C84 (Nagaiet al. 2014; Giovannini et al. 2018), and a triple-ridgestructure composed of the limb-brightened componentsand an additional central bright component has been
Corresponding author: Taiki [email protected] observed in the M87 jet (Asada et al. 2016; Hada et al.2017; Walker et al. 2018). The Event Horizon Telescope(EHT) has revealed the ring-like emission structure justaround the BH horizon of M87 (Event Horizon Tele-scope Collaboration et al. 2019a,b,c,d,e,f). These high-resolution observations have provided us hints for testingtheoretical jet models.The plausible driving mechanism of relativistic jets isthe Blandford-Znajek (BZ) process (Blandford & Znajek1977). The frame-dragging effect on a rotating BH mag-netosphere twists the poloidal magnetic field lines thatthread the ergosphere and creates the outward Poyntingflux extracting the energy and angular momentum fromthe BH (Koide et al. 2002; Komissarov 2004; Dermer &Menon 2009; Toma & Takahara 2014, 2016; Kinoshita &Igata 2018). General relativistic magnetohydrodynamic(GRMHD) simulations show that a highly magnetizedregion emerges around the rotational axis, where the BZprocess works (e.g. McKinney & Gammie 2004; McKin-ney & Blandford 2009; Tchekhovskoy et al. 2011; Taka-hashi et al. 2016; Nakamura et al. 2018; Porth et al. a r X i v : . [ a s t r o - ph . H E ] F e b Ogihara, Ogawa & Toma e + e − pair-creation gap formation (e.g.Blandford & Znajek 1977; Beskin et al. 1992; Broderick& Tchekhovskoy 2015; Hirotani et al. 2016; Levinson &Segev 2017; Kisaka et al. 2020), and/or the annihilationof high-energy photons from the accretion disk (Levin-son & Rieger 2011; Mościbrodzka et al. 2011; Kimura& Toma 2020). However, these non-thermal effects arenot considered in GRMHD simulations because of un-certainties in the physics of mass loading. Instead, theyset the artificial density floor values to prevent the nu-merical difficulty (McKinney & Gammie 2004; Riordanet al. 2018).Solving the steady equations do not need the densityfloor. The mass density distribution of the solutionsmay constrain the mass-loading mechanisms. The an-alytical solution of the steady axisymmetric GRMHDflow along a prescribed field line has been developedin the literature (Bekenstein & Oron 1978; Camenzind1986; Takahashi et al. 1990; Globus & Levinson 2014).Once the field line configuration and the Bernoulli pa-rameters are given, one can solve the Bernoulli equation(or wind equation) for each field line. Pu et al. (2015)showed an analytic solution for the single field line whichmimics the GRMHD simulation results (Tchekhovskoyet al. 2010; Nakamura et al. 2018). Pu & Takahashi(2020) presented a wind solution, which successfully passthe fast points, for multiple field lines with the pre-scribed toroidal field, but they do not consider the force-balance between the field lines (or the Grad-Shafranov(GS) equation). To solve the GS equation is computa-tionally demanding (Nitta et al. 1991; Fendt 1997; Be-skin 2010; Nathanail & Contopoulos 2014; Pan et al.2017; Mahlmann et al. 2018). Recently, Huang et al.(2019) and Huang et al. (2020) succeeded in construct-ing steady axisymmetric numerical solutions for themonopole and parabolic field line configurations by it-eratively solving the wind equation and the GS equa- tion. In Huang et al. (2020), they introduced the ‘load-ing zone’. The inner boundary of the loading zone isthe null-charge surface, and the outer boundary is thesurface where the Bernoulli equation’s solution of theoutflow becomes u p = 0 . The outflow and inflow startat these surfaces. They assumed the fluid number fluxper magnetic flux η as constant for different field lines,and did not try to constrain the mass-loading mecha-nism.One can also investigate the distribution of mass den-sity as well as other physical quantities by comparingobservations with the synthetic images derived by ra-diative transfer calculations in the jet models (Broderick& Loeb 2009; Porth et al. 2011; Dexter et al. 2012; Luet al. 2014; Mościbrodzka et al. 2014, 2016, 2017; Dexter2016; Jiménez-Rosales & Dexter 2018; Kawashima et al.2019; Chael et al. 2018, 2019; Davelaar et al. 2019; Chat-terjee et al. 2020; Jeter et al. 2020). Based on specialrelativistic force-free steady jet models, Takahashi et al.(2018) showed that a rapidly rotating BH accelerates theflow efficiently with remaining small toroidal velocity,leading to symmetric limb-brightened structure, which iscompatible with the large-scale observations of M87 jet,and Ogihara et al. (2019) showed that the velocity fieldstructure naturally produces the observed central ridgeemission of M87 jet by relativistic beaming effect (seealso Chernoglazov et al. 2019). They also demonstratedthat the synthetic images of large-scale jets strongly de-pend on the density distribution of jets near the BHs.Kawashima et al. (2020) performed GR radiative trans-fer calculations to show that the emission at the bot-tom of the separation surface can reproduce the ring-likeimage just around the M87 BH. Future EHT observa-tions of M87 will unveil the emission structure betweenthe limb-brightened one and the ring-like image (Hada2019). Parametric studies of steady GRMHD jet modelsand their comparison with the observations will extractthe mass density distribution and other specific condi-tions around the separation surface.In this paper, we introduce a new way to constructa steady axisymmetric GRMHD approximate solutionto examine the jet’s density distribution without be-ing suffered from the density floor problem. We fix thefield line configuration which mimics the ones of force-free or highly-magnetized GRMHD simulation results(Tchekhovskoy et al. 2008, 2010; Nakamura et al. 2018;Porth et al. 2019) with an additional term for obtainingtrans-fast-magnetosonic solutions (c.f. Beskin & Nokh-rina 2006; Pu et al. 2015). We numerically solve thetransverse force-balance between the field lines at theseparation surface to determine the distribution of theBernoulli parameters, which include η , and analytically atter density distributions of jets MODELTo examine the density distribution inside the jet, westudy the general relativistic equations for steady, ax-isymmetric, cold ideal MHD flows. We analytically solvethe dynamics parallel to the prescribed poloidal fieldlines and numerically keep the transverse force balanceat the separation surface.The spacetime geometry is given by the Kerr metricin the Boyer-Lindquist coordinates, ds = − (cid:18) − r Σ (cid:19) dt − ar sin θ Σ dtdφ + Σ∆ dr + Σ dθ + [( r + a ) − ∆ a sin θ ] sin θ Σ dφ , (1)where a is the dimensionless spin parameter, Σ = r + a cos θ , and ∆ = r − r + a . Hereafter we use the unitof c = G = 1 and the BH mass M = 1 . The radius ofthe event horizon is r H = 1 + √ − a , and the angularvelocity of the BH is ω H = a/ r H .2.1. Bernoulli equation
The basic equations are the energy-momentum equa-tion, the Maxwell equations, and the mass flux conser-vation law with the steady, axisymmetric, ideal MHDconditions. They lead to the four integral constantsalong a field line (i.e. the Bernoulli constants), whichare the total energy flux per particle E , the total angu-lar momentum flux per particle L , the number densityflux per magnetic flux η , and the so-called the ‘angularvelocity of the magnetic field’ Ω F (Bekenstein & Oron1978). They are given as follows: ˆ E (Ψ) = − u − Ω F B πµη , (2) ˆ L (Ψ) = u − B πµη , (3) η (Ψ) = − nu B G t = − nu B G t , (4) Ω F (Ψ) = F F = F F . (5)Here we have introduced the magnetic flux function Ψ( r, θ ) and the electromagnetic tensor F µν . For theaxisymmetric field, F A3 = ∂ A Ψ( A = 1 , . Thenthe magnetic field is given by B µ = 1 / (cid:15) νµλσ ξ ν F λσ ,where ξ µ = (1 , , , is the time-like Killing vector, (cid:15) µνλσ = √− g [ µνλσ ] , [ µνλσ ] is the permutation symbol,and g ≡ det( g µν ) . F A0 = E A is the electric field. n is the fluid-frame number density of the matter, u µ isthe four-velocity, and µ is the specific enthalpy. In thecold limit, µ is the rest energy of the plasma particle. ˆ E = E/µ , ˆ L = L/µ , and G t = g + Ω F g .Combination of Equations (2), (3), (4), and (5) isreduced to the Bernoulli equation, i.e., a fourth-orderequation of the poloidal velocity u p = √ u u + u u , (cid:88) i =0 A i u i p = 0 , (6)where A = 1 ,A = k B p πµηG t ,A = 1 + ˆ E k + (cid:18) k B p πµηG t (cid:19) ,A = B p ( k − ˆ E k )2 πµηG t ,A = k ( k − ˆ E k ) (cid:18) B p πµηG t (cid:19) . (7)Here, k = − ( g +2Ω F g +Ω F g ) , k = (1 − ˆ L Ω F / ˆ E ) , k = [ g + 2 g ˆ L/ ˆ E + g ( ˆ L/ ˆ E ) ] / ( g − g g ) , and B p = √ B B + B B is the poloidal magnetic field.Equation (6) is identical with Equation (34) in Pu et al.(2015) and Equation (23) in Huang et al. (2019).Given the integral constants { ˆ E (Ψ) , ˆ L (Ψ) , η (Ψ) , Ω F (Ψ) } and the field line configuration Ψ( r, θ ) , one can solvethe Bernoulli equation. Then one can derive the num-ber density distribution along the field line from thedefinition of η (Equation 4). The toroidal field can becalculated by B = − πµη G φ ˆ E + G t ˆ LM − k , (8)with the Alfven Mach number M = πµη /n and G φ = g + Ω F g . Ogihara, Ogawa & Toma model ν (cid:15) a Ω F (Ψ = 1) u p , ss ˆ E M0 . . H − P1 − . . H − − P2 − . . H − − P3 − .
95 0 . H − − P4 − . . H × − − P5 − . . H . × − − Table 1.
Parameter values used in our models.
Flux function model
We assume that Ψ( r, θ ) has a form of Ψ( r, θ ) = C (cid:20)(cid:18) rr H (cid:19) ν (1 − cos θ ) + 14 (cid:15)r sin θ (cid:21) , (9)where C is the normalization factor setting Ψ( r H , π/
2) =1 . ν and (cid:15) are the model parameters controlling thepoloidal field line shape. For (cid:15) = 0 , the field line con-figurations of ν = 0 and ν = 1 have the monopole andparabolic shape in the far zone, respectively. Ψ( r, θ ) with ν = 0 and (cid:15) = 0 is the exact solution of theforce-free magnetosphere in a Schwartzschild spacetime(Blandford & Znajek 1977), and Ψ( r, θ ) with ν = 1 and (cid:15) = 0 is the dominant term of the exact solu-tion of the force-free magnetosphere in a Schwarzschildspacetime (Blandford & Znajek 1977; Lee & Park 2004).The jet-disk boundary in GRMHD simulations definedas the magnetic-to-matter energy flux ratio σ = | − Ω F B / (4 πµηu ) | = 1 matches the parabolic configura-tions with a constant ν in a large computational domain(Nakamura et al. 2018; Porth et al. 2019). (cid:15) representsa small disturbance from the force-free field lines, whichmakes the outflow accelerate by converting the Poyntingflux to the kinetic energy flux (c.f. Beskin & Nokhrina2006; Pu et al. 2015). When (cid:15) = 0 , B p r sin θ becomesconstant along the field line in a far zone (Pu & Taka-hashi 2020). For the outflow to pass through the fastmagnetosonic point, B p r sin θ needs to decrease withthe radius around the fast magnetosonic point (Begel-man & Li 1994; Beskin 2010; Toma & Takahara 2013).When (cid:15) > , the field lines are collimated and B p r sin θ decreases. Then, the outflow can pass through the fastmagnetosonic point.In this paper, we consider two models, a split-monopole configuration model ( ν = 0 , (cid:15) = 0 ), and a per-turbed parabolic configuration model ( ν = 1 , (cid:15) = 10 − ).The aim of performing the split-monopole model isto check if our method can approximately reproducethe split-monopole force-free solution of slowly rotat-ing Kerr BH magnetosphere, in which Ω F ≈ . H and ˆ E (Ψ) ∝ sin θ (Blandford & Znajek 1977). We set a = 0 . because such a small BH spin does not change . . . . . . . . . Ω F / Ω H Ψ . . . . . .
03 0 0 . . . . L i n Ω F / E i n Ψ Figure 1. Ω F (Ψ) and L in Ω F /E in (Ψ) of the M0 model. Ourresult satisfies the conditions of the force-free split-monopolemagnetosphere, Ω F = 0 . H and L Ω F /E = 1 , within 1%accuracy in Ψ > . . the field line structure significantly from the monopoleforce-free solution near the horizon (Tanabe & Nagataki2008; Tchekhovskoy et al. 2010). After confirming thatour model can reproduce the split-monopole analyticsolution well, we perform a calculation of a perturbedparabolic magnetosphere around a rapidly rotating BH( a = 0 . ) to investigate the density distribution in thejet.We list the parameter values used in our models in Ta-ble 1. The M0 model is the split-monopole configurationmodel, and the results are shown in Section 3.1. The P1model is our fiducial model in the parabolic configura-tion models. We analyze the details of the P1 model inSection 3.2. The other parabolic configuration models(P2, P3, P4, and P5) are also analyzed to investigatethe parameter dependence. atter density distributions of jets − − − u p r − r H Ψ = 0 .
1Ψ = 0 .
5Ψ = 0 . . . . − − B r − r H Ψ = 0 .
1Ψ = 0 .
5Ψ = 0 . Figure 2. u p , B along the field lines Ψ = 0 . , . , and . of the split-monopole configuration model. We note that u p is the absolute value of the poloidal velocity. The constant B along the field line represents that the energy conversionfrom the Poynting flux to the kinetic energy flux is inefficientin this model. Determining the Bernoulli parameters
Solving the wind equation (Equation 6) requires theBernoulli parameters { ˆ E (Ψ) , ˆ L (Ψ) , Ω F (Ψ) , η (Ψ) } . Wedetermine ˆ L (Ψ) , η (Ψ) , and Ω F (Ψ) in the same methodfor both of the field line configuration models, and de-termine ˆ E (Ψ) differently. ˆ L (Ψ) and η (Ψ) are determined simultaneously fromthe assumed poloidal velocity at the separation surfaceand the Znajek condition (Znajek 1977), u p ( r ss , θ ) = u p , ss = const . (10) B ( r H , θ ) = − (cid:34)(cid:18) g g (cid:19) / ( ω H − Ω F ) ∂ θ Ψ (cid:35) r=r H , (11) where r ss denotes the radius of the separation surface.Combining u p , ss and the Znajek condition with Equation(6) at the separation surface and Equation (8) at thehorizon, one can obtain the two equations for solving ˆ L and η . We set u p , ss = 10 − for the M0, P1, P2, and P3model. The P2 and P3 models are used to investigatethe BH spin dependence, while the P4 and P5 models areused to investigate the u p , ss dependence. The separationsurface is defined by ∂ p k = 0 , where ∂ p ≡ B ∂ + B ∂ is the differentiation in the direction along the field line. Ω F (Ψ) is determined to satisfy the force balance be-tween field lines at the separation surface. The compo-nent of the equation of motion perpendicular to the fieldline is e (n)A [ ρu ν u A ; ν − F Aν I ν ] = 0 , (12)where e (n)A is the unit vector in the direction perpen-dicular to the field line in the poloidal plane e (n)A = E A / (cid:112) E B E B ( A, B = 1 , ). We assume Ω F (Ψ = 1) =0 . H for the M0 model, and Ω F (Ψ = 1) = 0 . H for the parabolic configuration models. Then we choose Ω F (Ψ) by minimizing the left-hand side of Equation(12). We do this only at the separation surface to ob-tain the approximate solutions. We confirm that thismethod works well for the M0 model and then apply itto the parabolic configuration models. We also evaluatehow well our approximate solutions keep the transversebalance at the regions other than the separation surface.Regarding of ˆ E (Ψ) , we take different treatments forthe split-monopole and parabolic cases. (i): For the M0 model, we set ˆ E = ˆ E sin θ and ˆ E =10 , following the dependence of the Poynting fluxin the force-free monopole magnetosphere (Bland-ford & Znajek 1977). The small value of u p , ss istaken for the flow to be Poynting-flux dominated. (ii): For the parabolic configuration models, we deter-mine ˆ E for the outflow to accelerate smoothlypassing through the fast magnetosonic point. If ˆ E is too low, u p diverges after passing the Alfvenpoint, while if ˆ E is too large, u p does not di-verge, but M does not reach M ≡ k +( G B ) / ( ρ B ) , where M = M should be sat-isfied at the fast magnetosonic point (Takahashiet al. 1990). We adjust ˆ E so that u p does not di-verge and there is a point where | − M /M | < − .We iteratively adjust ˆ E and Ω F to obtain theirvalues satisfying both of the condition for smoothacceleration and the force balance at an expectedprecision. Ogihara, Ogawa & Toma
From Equation (4), η has the opposite sign for the in-flow and outflow. The sign of the second term of Equa-tion (2) and (3) change between the inflow and outflow.On the other hand, the sign of the first term u doesnot change, and u is almost − at the separation sur-face. Due to this, ˆ E and ˆ L must have different valuesfor the inflow and outflow to keep the electromagneticfield smooth and continuous at the separation surface.Thus, we set ˆ E in = − ˆ E out + 2 , where ˆ E in and ˆ E out arethe values for inflow and outflow, respectively. RESULTSIn Section 3.1, we perform the split-monopole con-figuration model and confirm that our model, whichtakes the trans-field balance only at the separation sur-face, can make the electromagnetic field consistent withthe force-free solution. In Section 3.2, we perform theparabolic configuration model and show the density dis-tribution inside the jet.3.1.
Split-monopole configuration model
Figure 1 shows Ω F (Ψ) and ˆ L out Ω F / ˆ E out (Ψ) at the sep-aration surface. Ω F / Ω H = 0 . and the force-free condi-tion ˆ L out Ω F / ˆ E out = 1 are satisfied within 1% accuracyin Ψ > . , which mean that our approximate solu-tions are consistent with the force-free monopole solu-tion (Blandford & Znajek 1977). The deviation fromthe monopole force-free solution decreases, as either theBH spin is smaller, ˆ E is larger, u p , ss is smaller, or Ω F (Ψ = 1) is closer to . H .Figure 2 shows u p and B along the field lines of Ψ = 0 . , . , . . The outflow does not pass through thefast magnetosonic point, unlike in the parabolic fieldconfiguration case, as discussed in Camenzind (1986). B of each flow is almost constant along the field lineunless it diverges. This means that the conversion fromthe Poynting flux to the fluid energy flux is inefficientin the monopole field configuration and that the electro-magnetic field is almost force-free in the whole region.3.2. Parabolic configuration model
P1 model
In this subsection, we focus on the results of the cal-culation of the P1 model. We show the two dimensionaldistribution of u p , n , and σ in Figure 3. Here, we nor-malize the number density by n norm ≡ (cid:20) B B + B B + B B πµ (cid:21) (r=r ss , Ψ=1) . (13)The inflow and outflow smoothly accelerate from theseparation surface to relativistic speeds, and the densitydecreases with the distance from the separation surface. We note that the density does not diverge at the sepa-ration surface since u p , ss is not zero. At the separationsurface, σ has the almost same value as ˆ E (Ψ) because | u | ≈ (see the top left panel of Figure 4). Along theseparation surface, it increases from ≈ at the jet edgeto the maximum value ≈ , and then decreases to ≈ near the axis. Along the field line, σ decreases as theradius increases for both the inflow and outflow.Figure 4 shows η (Ψ) E (Ψ) , ˆ E (Ψ) , Ω F (Ψ) , and n/n norm at the separation surface. ˆ E (Ψ) has a peak at Ψ ∼ . . The Poynting flux becomes zero at the axis, whichmeans ˆ E (Ψ = 0) = − u ≈ . Ω F increases toward Ω F = 0 . H from the edge to the axis, but it decreasesnear the axis at Ψ ≈ . . ˆ E ≈ ˆ L Ω F is satisfied within1% accuracy. This means that the flow is Poynting fluxdominated at the separation surface. ηE roughly follows ∝ sin θ H , while this dependence is that of B on θ H at the horizon for a (cid:28) (Equation 11). The numberdensity at the separation surface has the peak at the jetedge and decreases to nearly zero toward the jet axis.Figure 5 shows M , M , and M along the fieldlines of Ψ = 0 . , . , and . , where M ≡ k . Theintersections of M and M are the Alfven points,and the ones of M and M are the fast magne-tosonic points. Figure 6 shows u p along the field lines Ψ = 0 . , . , and . . Both the inflow and outflow startfrom the separation surface and accelerate to relativisticspeeds.We evaluate the trans-field force-balance by introduc-ing χ ≡ | ( f + − | f − | ) / ( f + + | f − | ) | , where we gather thepositive components of Equation (12) to f + and the neg-ative ones to f − . χ ranges from 0 to 1, and χ = 0 means the complete force balance. The results showthat χ is less than − for all the field lines at theseparation surface. For the outflow above the separa-tion surface, χ increases rapidly to ∼ − and thenturns to decrease, while for the inflow, it increases upto ∼ and then decreases. The Lorentz force di-rects toward the axis e (n)A F Aν I ν > and the inertialforce directs in the opposite way − e (n)A ρu A u ν ; ν < .These forces are well balanced at the separation sur-face e (n)A ρu A u ν ; ν ≈ e (n)A F Aν I ν . In other radii than theseparation surface, the Lorentz force is larger than theinertial one, and the net force is directing for the flowto collimate. 3.2.2. Parameter dependences
We calculate the parabolic configuration models withdifferent parameter values listed in Table 1 to investigatethe dependencies of the density distribution on the BHspin a and u p , ss . atter density distributions of jets Figure 3.
Two dimensional distribution of u p , n/n norm , and σ of the P1 model. We use a = 0 . and . for the P2 and P3 mod-els, respectively, to investigate the BH spin dependence. ˆ E, Ω F , ηE, n/n norm are shown in Figure 4. Interestingly,as a becomes larger, the density gets larger near the jetedge and smaller near the axis, while ˆ E changes in theopposite way. ηE of the P2, P3 models also roughlyfollow ∝ sin θ H . ηE becomes larger in all the field lineswith a because of the increase of B ( r = r H ) and thePoynting flux. Ω F also increases with a . The force-freeness − L Ω F /E is smaller than 0.2 for all threemodels. The minimum value − L Ω F /E ≈ . real-izes where ˆ E is maximum.We perform calculations with different u p , ss . We use u p , ss = 6 × − and . × − for the P4 and P5models, respectively. The results are shown in Figure7. When u p , ss is smaller, n/n norm changes in a similarfashion as a gets larger. n/n norm at the jet edge changesproportional to u − , ss . The P1, P4, and P5 models showthat ηE does not significantly depend on u p , ss . As u p , ss decreases, ˆ E increases and η decreases for all the fieldlines. DISCUSSION4.1.
Density on the separation surface
The matter density distribution will constrain themass-loading mechanism. In Figure 8, we show n/n norm of the P1, P2, and P3 models as a function of r ss . n/n norm is largest at the jet edge and decreases as Ψ get smaller as shown in Figure 4. r ss decreases as a getslarger. At the far zone, the normalized density roughlyfollows n/n norm ∝ r − in all the models. Figure 9 shows n/n norm as a function of r ss for the different values of u p , ss . We also have the dependence n/n norm ∝ r − in the far zone in these models. For r ss < , the r ss de-pendence of n/n norm is steeper.The annihilation of high-energy photons from the ac-cretion disk is one of the proposed mass-loading mech-anisms (Levinson & Rieger 2011; Mościbrodzka et al.2011; Kimura & Toma 2020). This process leads to the e + e − density distribution n ∝ r − for the case in which γ -ray emitting region is compact near the BH (Mości-brodzka et al. 2011), and then the particles are injectedmainly at the base for the outflow, i.e., at the sepa-ration surface, which is similar to the situation of ourmodel. However, the dependence r − appears too steepcompared with the results in our model. If γ -ray emit-ting region is extended, say around r ∼ , in the ac-cretion flow (Kimura & Toma 2020), the r dependenceof the e + e − density can be much shallower and mightbe consistent with our results. However, it should benoted that this e + e − injection model can provide par-ticle number density sufficient for screening the sparkgap (i.e., larger than the Goldreich-Julian number den-sity), but not sufficient for the radio synchrotron fluxof M87 jet (Kimura & Toma 2020). Other injectionmechanisms such as magnetic reconnection (e.g. Parfreyet al. 2015; Mahlmann et al. 2020) and/or fluid insta-bility (e.g. Globus & Levinson 2016; Nakamura et al.2018; Chatterjee et al. 2019; Sironi et al. 2020) could beefficient for injecting electrons that produce the limb-brightened radio emission.The results of our model are compatible with the ob-served limb-brightened emission structure of jets (seealso Figure 3), although it is uncertain what fraction ofthe matter contribute to the non-thermal emission. Therelative amount of the density near the axis compared tothe one near the jet edge is also important because the Ogihara, Ogawa & Toma . . . . − ˆ E i n Ψ a = 0 . a = 0 . a = 0 . . . . . . . . . . Ω F / Ω H Ψ a = 0 . a = 0 . a = 0 . . . . . . . . . . . . η i n E i n Ψ a = 0 . a = 0 . a = 0 . ∝ sin θ H . . . . . n / n n o r m Ψ a = 0 . a = 0 . a = 0 . Figure 4. ˆ E in , Ω F , η in E in , and n/n norm at the separation surface. The black, blue, and red lines are the results of P1, P2, andP3 models, respectively. The number density at the separation surface is concentrated more near the jet edge when the BH spin a is larger. emission near the axis will be Doppler-boosted, formingthe central ridge of the observed triple-ridge emissionstructure of M87 jet (Ogihara et al. 2019).4.2. Comparison to other studies
The distribution of the Bernoulli parameters can becompared to the results of other studies. Huang et al.(2020) numerically solved the GS equation for the wholeregion. They set the loading zone between the sepa-ration surface and the null-charge surface, although itis unclear whether such a large loading zone is nec-essary for obtaining the solutions. ˆ E (Ψ) has a peaknear the jet edge in their result, while our models havethe one relatively closer to the axis. They showedthat Ω F monotonically increases toward the axis, and Ω F (Ψ = 0) = 0 . H , which is set as a boundary condi-tion, while in our model, Ω F (Ψ) decreases rapidly nearthe axis. η of their model is assumed by a given mag-netization parameter and the poloidal magnetic field atthe null-charge surface. It is larger near the axis, whichis the opposite trend from our results. Ω F (Ψ) distribution is also shown in Beskin & Zhel-toukhov (2013), in which they solve the GS equation ofa cylindrical jet in a special relativistic regime. In theirresult, Ω F (Ψ) decreases near the axis like our results.This trend is also seen in the GRMHD simulation ofMcKinney et al. (2012).4.3. Magnetic bending profile atter density distributions of jets − − − − − − M , M A l f , M f a s t r − r H M (Ψ = 0 . M (Ψ = 0 . M (Ψ = 0 . − − − − − − M , M A l f , M f a s t r − r H M (Ψ = 0 . M (Ψ = 0 . M (Ψ = 0 . − − − − − − M , M A l f , M f a s t r − r H M (Ψ = 0 . M (Ψ = 0 . M (Ψ = 0 . Figure 5. M , M , and M along the field lines of Ψ =0 . , . , and . in the P1 model. The intersections of M and M are the Alfven points and the ones of M and M are the fast magnetosonic points. − − − u p r − r H Ψ = 0 .
1Ψ = 0 .
5Ψ = 0 . Figure 6. u p along the field lines Ψ = 0 . , . , and . inthe P1 model. We note that u p is the absolute value of thepoloidal velocity. Pu & Takahashi (2020) introduced the reasonableshape of the function of E p , ZAMO /B T , ZAMO , which canbe rewritten by the bending angle of the field line,and derived wind solutions with the prescribed func-tion (see also Tomimatsu & Takahashi 2003; Takahashi& Tomimatsu 2008). E p , ZAMO = | G φ B p / ( G t √ g ) | is the poloidal electric field strength and B T , ZAMO = | B / ( α √ g ) | is the toroidal magnetic field strength inthe zero angular momentum observer frame.There are some constraints to this function. First, E p , ZAMO /B T , ZAMO = 1 at the horizon (i.e. theZnajek condition). Second, E p , ZAMO /B T , ZAMO = 0 at the null-charge surface, where the field line coro-tate with the spacetime ( − g /g = Ω F ). Finally, ( E p , ZAMO /B T , ZAMO ) < − / ˆ E for the outflow inorder to prevent u p from diverging. For the outflow, E p , ZAMO /B T , ZAMO needs to increase for the flow to ac-celerate. Additionally, E p , ZAMO /B T , ZAMO should besmooth and continuous.The previous studies mentioned above prescribed ( E p , ZAMO /B T , ZAMO ) as a constant value for out-flow and derived u p using it. The assumed ( E p , ZAMO /B T , ZAMO ) and the derived one using Equa-tion (8) was not self-consistent. We solved Equation(6) which does not explicitly depend on B , and de-rive ( E p , ZAMO /B T , ZAMO ) afterwards. We show theself-consistent profile of E p , ZAMO /B T , ZAMO of Ψ =0 . , . ,and . in Figure 10. ( E p , ZAMO /B T , ZAMO ) inour result follows all the conditions listed above.4.4. Model limitations Ogihara, Ogawa & Toma . . . . − ˆ E i n Ψ u p,ss = 6 × − u p,ss = 1 × − u p,ss = 1 . × − . . . . . . . . . Ω F / Ω H Ψ u p,ss = 6 × − u p,ss = 1 × − u p,ss = 1 . × − . . . . . . . . . η i n E i n Ψ u p,ss = 6 × − u p,ss = 1 × − u p,ss = 1 . × − ∝ sin θ H . . . . . n / n n o r m Ψ u p,ss = 6 × − u p,ss = 1 × − u p,ss = 1 . × − Figure 7. ˆ E in , Ω F , η in E in , and n/n norm at the separation surface. The black, blue, and red lines are the results of P1, P4,and P5 models, respectively. The number density at the separation surface is concentrated more near the jet edge when u p , ss issmaller. The inflows do not pass through the fast magnetosonicpoint in our results. The inflow diverges at the regionvery close to the horizon. Additionally, the indicator ofthe force-balance between the field lines χ is also largein the inflow region. We try to find different values of ˆ E in from those satisfying the trans-field force-balancefor the inflows to pass through the fast magnetosonicpoints, and find that | ˆ E in | needs to be smaller only bya few in the case of the P1 model. The difference of E in is smaller for larger Ψ . Adjustment of Ψ( r, θ ) should beconsidered in future work.We note that if more particles are injected at regionscloser to the BH, like in the case of the annihilationof high-energy photons, the inflows are highly affectedby mass loading and should be modeled with varying η along a field line (cf. Globus & Levinson 2014), whilethe current model of outflows is still applicable.The distribution of u p , ss (Ψ) may differ between theinjection mechanisms. In our model, we assume that u p , ss (Ψ) is constant for all the field line. If we changethe distribution, it may change the density profile alongthe separation surface. This will be discussed in a sep-arate work. Injection via instabilities and the magneticreconnection may occur by the interaction between thejet edge and the disk wind (Globus & Levinson 2016;Nakamura et al. 2018; Chatterjee et al. 2019; Sironi et al.2020). This may change the field line configuration andthe density distribution near the jet edge.Although our current approximate solutions ofGRMHD jets have the limitations mentioned above, atter density distributions of jets − −
10 100 n / n n o r m r ss a = 0 . a = 0 . a = 0 . r − Figure 8. n/n norm as a function of r ss . The black, blue,and red lines are the results of the P1, P2, and P3 models,respectively. When the BH spin is larger, r ss decreases. Atthe far zone, n/n norm roughly follows ∝ r − . − −
10 100 n / n n o r m r ss u p,ss = 6 × − u p,ss = 1 × − u p,ss = 1 . × − r − Figure 9. n/n norm as a function of r ss . The black, blue,and red lines are the results of the P1, P4, and P5 models,respectively. When the BH spin is larger, r ss decreases. Atthe far zone, n/n norm roughly follows ∝ r − . they are useful not only for constraining the mass-loading mechanism but also for obtaining key ingredi-ents from the observed complex emission structures. In-deed, by comparing the approximate special relativisticforce-free steady jet models with observations of far-zonejets, it is found that the special relativistic beaming ef-fect is essential for understanding the symmetry of theobserved limb-brightened structure and the central ridgeemission and that the emission image strongly dependson the matter density distribution at the base of outflow − − − − ( E p , Z A M O / B T , Z A M O ) r − r H Ψ = 0 .
9Ψ = 0 .
5Ψ = 0 . Figure 10.
The radial profiles of ( E p , ZAMO /B T , ZAMO ) of Ψ = 0 . , . , and . in the P1 model. They are smoothand continuous and satisfy ( E p , ZAMO /B T , ZAMO ) = 1 at thehorizon from the Znajek condition, ( E p , ZAMO /B T , ZAMO ) =0 at the null-charge surface, and ( E p , ZAMO /B T , ZAMO ) < − / ˆ E in the outflow region. (Takahashi et al. 2018; Ogihara et al. 2019). This hasmotivated us to develop GRMHD steady jet models. SUMMARY & PROSPECTSWe constructed a steady, axisymmetric GRMHD jetmodel to examine the density distribution inside thejet. We fixed the poloidal field line configuration, whichmimics the ones of force-free or GRMHD simulation re-sults with an additional term for obtaining trans-fast-magnetosonic outflows, and assumed the poloidal veloc-ity of the flow at the separation surface u p , ss as con-stant for different field lines. We numerically solved thetrans-field force-balance (Equation 12) at the separationsurface to determine the Bernoulli parameters and ana-lytically solved the fourth-order equation of the poloidalvelocity along the field lines (Equation 6).We consider the two field line configurations; (i)the split-monopole configuration model and (ii) theparabolic configuration model. (i) In the split-monopoleconfiguration model, we confirmed that our methodcould reproduce the BZ solution well except aroundthe axis. (ii) In the parabolic configuration model,we obtained approximate solutions of highly-magnetizedGRMHD flows. The force-balance between field lines atthe separation surface is satisfied with high accuracy.The outflow successfully passes the fast magnetosonicpoint and satisfies a good force-balance. We found thatthe number density distribution at the separation sur-face roughly scales as n ( r = r ss , Ψ) ∝ r − in the farzone. We examined the parabolic configuration models2 Ogihara, Ogawa & Toma with the different BH spins a and u p , ss . As the BH spinincreases, the density at the separation surface is con-centrated more at the jet edge. We obtained the similartrend when u p , ss is smaller.Our semi-analytic model can be utilized to examinethe density distribution of the highly magnetized jet re-gion, while current GRMHD simulations cannot sincethey use the artificial density floor, and can survey alarger parameter space in detail because of the smallercomputational budget. It will be interesting to apply ra-diative transfer calculations to our model and comparethem to the observed emission structure, as done withthe special relativistic force-free jet models (Takahashiet al. 2018; Ogihara et al. 2019).Observations of core shifts, collimation profiles, propermotions of blobs, and Faraday rotation maps in variousfrequencies and scales (Kovalev et al. 2007; Asada et al.2013; Hada et al. 2011, 2013, 2018; Nakamura & Asada2013; Mertens et al. 2016; Park et al. 2019; Park et al. 2019; Kim et al. 2020; Akiyama et al. 2018; Hiura et al.2018; Giovannini et al. 2018; Kravchenko et al. 2020) willalso be useful for testing jet models (Kino et al. 2014,2015; Chael et al. 2019; Jeter et al. 2020). Future EHTand EATING-VLBI observations will provide crucial in-formation on the jet launching region of M87 (Hada &Eavn/Eating VLBI Collaboration 2020).ACKNOWLEDGMENTSWe thank the anonymous referee for useful comments.We thank Shigeo Kimura, Hung-Yi Pu, Fumio Taka-hara, and Masaaki Takahashi for fruitful discussions.T.O. acknowledges financial support from the GraduateProgram on Physics for the Universe of Tohoku Univer-sity. This work is partly supported by Grant-in-Aid forJSPS Fellowship 20J14537 (T.O.) and KAKENHI No.18H01245 (K.T.).REFERENCES Akiyama, K., Asada, K., Fish, V., et al. 2018, Galaxies, 6,15, doi: 10.3390/galaxies6010015Asada, K., Nakamura, M., Doi, A., Nagai, H., & Inoue, M.2013, The Astrophysical Journal, 781, L2,doi: 10.1088/2041-8205/781/1/l2Asada, K., Nakamura, M., & Pu, H.-Y. 2016, TheAstrophysical Journal, 833, 56,doi: 10.3847/1538-4357/833/1/56Begelman, M. C., & Li, Z.-Y. 1994, ApJ, 426, 269,doi: 10.1086/174061Bekenstein, J. D., & Oron, E. 1978, Physical Review D, 18,1809, doi: 10.1103/physrevd.18.1809Beskin, V. S. 2010, MHD Flows in Compact AstrophysicalObjects, doi: 10.1007/978-3-642-01290-7Beskin, V. S., Istomin, Y. N., & Parev, V. I. 1992, SovietAstronomy, 36, 642Beskin, V. S., & Nokhrina, E. E. 2006, Monthly Notices ofthe Royal Astronomical Society, 367, 375,doi: 10.1111/j.1365-2966.2006.09957.xBeskin, V. S., & Zheltoukhov, A. A. 2013, AstronomyLetters, 39, 215, doi: 10.1134/s1063773713040014Blandford, R. D., & Znajek, R. L. 1977, Monthly Notices ofthe Royal Astronomical Society, 179, 433,doi: 10.1093/mnras/179.3.433Boccardi, B., Krichbaum, T. P., Bach, U., Bremer, M., &Zensus, J. A. 2016, Astronomy & Astrophysics, 588, L9,doi: 10.1051/0004-6361/201628412Broderick, A. E., & Loeb, A. 2009, The AstrophysicalJournal, 697, 1164, doi: 10.1088/0004-637x/697/2/1164 Broderick, A. E., & Tchekhovskoy, A. 2015, TheAstrophysical Journal, 809, 97,doi: 10.1088/0004-637x/809/1/97Camenzind, M. 1986, A&A, 162, 32Chael, A., Narayan, R., & Johnson, M. D. 2019, MonthlyNotices of the Royal Astronomical Society, 486, 2873,doi: 10.1093/mnras/stz988Chael, A., Rowan, M. E., Narayan, R., Johnson, M. D., &Sironi, L. 2018, doi: 10.1093/mnras/sty1261Chatterjee, K., Liska, M., Tchekhovskoy, A., & Markoff,S. B. 2019, MNRAS, 490, 2200,doi: 10.1093/mnras/stz2626Chatterjee, K., Younsi, Z., Liska, M., et al. 2020, MNRAS,499, 362, doi: 10.1093/mnras/staa2718Chernoglazov, A. V., Beskin, V. S., & Pariev, V. I. 2019,Monthly Notices of the Royal Astronomical Society, 488,224, doi: 10.1093/mnras/stz1683Davelaar, J., Olivares, H., Porth, O., et al. 2019,Astronomy & Astrophysics, 632, A2,doi: 10.1051/0004-6361/201936150Dermer, C. D., & Menon, G. 2009, High Energy Radiationfrom Black Holes: Gamma Rays, Cosmic Rays, andNeutrinosDexter, J. 2016, Monthly Notices of the Royal AstronomicalSociety, 462, 115, doi: 10.1093/mnras/stw1526Dexter, J., McKinney, J. C., & Agol, E. 2012, MonthlyNotices of the Royal Astronomical Society, 421, 1517,doi: 10.1111/j.1365-2966.2012.20409.x atter density distributions of jets Event Horizon Telescope Collaboration, Akiyama, K.,Alberdi, A., et al. 2019a, ApJL, 875, L1,doi: 10.3847/2041-8213/ab0ec7—. 2019b, ApJL, 875, L2, doi: 10.3847/2041-8213/ab0c96—. 2019c, ApJL, 875, L3, doi: 10.3847/2041-8213/ab0c57—. 2019d, ApJL, 875, L4, doi: 10.3847/2041-8213/ab0e85—. 2019e, ApJL, 875, L5, doi: 10.3847/2041-8213/ab0f43—. 2019f, ApJL, 875, L6, doi: 10.3847/2041-8213/ab1141Fendt, C. 1997, A&A, 319, 1025Giovannini, G., Savolainen, T., Orienti, M., et al. 2018,Nature Astronomy, 2, 472,doi: 10.1038/s41550-018-0431-2Globus, N., & Levinson, A. 2014, The AstrophysicalJournal, 796, 26, doi: 10.1088/0004-637x/796/1/26Globus, N., & Levinson, A. 2016, MNRAS, 461, 2605,doi: 10.1093/mnras/stw1474Hada, K. 2019, Galaxies, 8, 1, doi: 10.3390/galaxies8010001Hada, K., Doi, A., Kino, M., et al. 2011, Nature, 477, 185,doi: 10.1038/nature10387Hada, K., & Eavn/Eating VLBI Collaboration. 2020, inPerseus in Sicily: From Black Hole to Cluster Outskirts,ed. K. Asada, E. de Gouveia Dal Pino, M. Giroletti,H. Nagai, & R. Nemmen, Vol. 342, 73–76,doi: 10.1017/S1743921318005173Hada, K., Kino, M., Doi, A., et al. 2013, The AstrophysicalJournal, 775, 70, doi: 10.1088/0004-637x/775/1/70—. 2016, The Astrophysical Journal, 817, 131,doi: 10.3847/0004-637x/817/2/131Hada, K., Park, J. H., Kino, M., et al. 2017, Publications ofthe Astronomical Society of Japan, 69,doi: 10.1093/pasj/psx054Hada, K., Doi, A., Wajima, K., et al. 2018, TheAstrophysical Journal, 860, 141,doi: 10.3847/1538-4357/aac49fHirotani, K., Pu, H.-Y., Lin, L. C.-C., et al. 2016, TheAstrophysical Journal, 833, 142,doi: 10.3847/1538-4357/833/2/142Hiura, K., Nagai, H., Kino, M., et al. 2018, Publications ofthe Astronomical Society of Japan, 70,doi: 10.1093/pasj/psy078Huang, L., Pan, Z., & Yu, C. 2019, The AstrophysicalJournal, 880, 93, doi: 10.3847/1538-4357/ab2909—. 2020, The Astrophysical Journal, 894, 45,doi: 10.3847/1538-4357/ab86a3Jeter, B., Broderick, A. E., & Gold, R. 2020, MNRAS, 493,5606, doi: 10.1093/mnras/staa679Jiménez-Rosales, A., & Dexter, J. 2018, Monthly Notices ofthe Royal Astronomical Society, 478, 1875,doi: 10.1093/mnras/sty1210 Kawashima, T., Kino, M., & Akiyama, K. 2019, TheAstrophysical Journal, 878, 27,doi: 10.3847/1538-4357/ab19c0Kawashima, T., Toma, K., Kino, M., et al. 2020, arXive-prints, arXiv:2009.08641.https://arxiv.org/abs/2009.08641Kim, J.-Y., Krichbaum, T. P., Lu, R.-S., et al. 2018,Astronomy & Astrophysics, 616, A188,doi: 10.1051/0004-6361/201832921Kim, J.-Y., Krichbaum, T. P., Broderick, A. E., et al. 2020,A&A, 640, A69, doi: 10.1051/0004-6361/202037493Kimura, S. S., & Toma, K. 2020, arXiv e-prints,arXiv:2003.13173. https://arxiv.org/abs/2003.13173Kino, M., Takahara, F., Hada, K., et al. 2015, TheAstrophysical Journal, 803, 30,doi: 10.1088/0004-637x/803/1/30Kino, M., Takahara, F., Hada, K., & Doi, A. 2014, TheAstrophysical Journal, 786, 5,doi: 10.1088/0004-637x/786/1/5Kinoshita, S., & Igata, T. 2018, Progress of Theoretical andExperimental Physics, 2018, 033E02,doi: 10.1093/ptep/pty024Kisaka, S., Levinson, A., & Toma, K. 2020, TheAstrophysical Journal, 902, 80,doi: 10.3847/1538-4357/abb46cKoide, S., Shibata, K., Kudoh, T., & Meier, D. L. 2002,Science, 295, 1688, doi: 10.1126/science.1068240Komissarov, S. S. 2004, Monthly Notices of the RoyalAstronomical Society, 350, 427,doi: 10.1111/j.1365-2966.2004.07598.xKomissarov, S. S., Barkov, M. V., Vlahakis, N., & Königl,A. 2007, Monthly Notices of the Royal AstronomicalSociety, 380, 51, doi: 10.1111/j.1365-2966.2007.12050.xKomissarov, S. S., Vlahakis, N., & Königl, A. 2010,MNRAS, 407, 17, doi: 10.1111/j.1365-2966.2010.16779.xKovalev, Y. Y., Lister, M. L., Homan, D. C., &Kellermann, K. I. 2007, The Astrophysical Journal, 668,L27, doi: 10.1086/522603Kravchenko, E., Giroletti, M., Hada, K., et al. 2020, A&A,637, L6, doi: 10.1051/0004-6361/201937315Lee, H., & Park, J. 2004, Physical Review D, 70,doi: 10.1103/physrevd.70.063001Levinson, A., & Rieger, F. 2011, The AstrophysicalJournal, 730, 123, doi: 10.1088/0004-637x/730/2/123Levinson, A., & Segev, N. 2017, Physical Review D, 96,doi: 10.1103/physrevd.96.123006Lu, R.-S., Broderick, A. E., Baron, F., et al. 2014, TheAstrophysical Journal, 788, 120,doi: 10.1088/0004-637x/788/2/120 Ogihara, Ogawa & Toma
Lyubarsky, Y. E. 2009, Monthly Notices of the RoyalAstronomical Society, 402, 353,doi: 10.1111/j.1365-2966.2009.15877.xMahlmann, J. F., Cerdá-Durán, P., & Aloy, M. A. 2018,Monthly Notices of the Royal Astronomical Society, 477,3927, doi: 10.1093/mnras/sty858Mahlmann, J. F., Levinson, A., & Aloy, M. A. 2020,MNRAS, 494, 4203, doi: 10.1093/mnras/staa943McKinney, J. C., & Blandford, R. D. 2009, MonthlyNotices of the Royal Astronomical Society: Letters, 394,L126, doi: 10.1111/j.1745-3933.2009.00625.xMcKinney, J. C., & Gammie, C. F. 2004, The AstrophysicalJournal, 611, 977, doi: 10.1086/422244McKinney, J. C., Tchekhovskoy, A., & Blandford, R. D.2012, Monthly Notices of the Royal AstronomicalSociety, 423, 3083, doi: 10.1111/j.1365-2966.2012.21074.xMertens, F., Lobanov, A. P., Walker, R. C., & Hardee,P. E. 2016, Astronomy & Astrophysics, 595, A54,doi: 10.1051/0004-6361/201628829Mościbrodzka, M., Dexter, J., Davelaar, J., & Falcke, H.2017, Monthly Notices of the Royal AstronomicalSociety, 468, 2214, doi: 10.1093/mnras/stx587Mościbrodzka, M., Falcke, H., & Shiokawa, H. 2016,Astronomy & Astrophysics, 586, A38,doi: 10.1051/0004-6361/201526630Mościbrodzka, M., Falcke, H., Shiokawa, H., & Gammie,C. F. 2014, Astronomy & Astrophysics, 570, A7,doi: 10.1051/0004-6361/201424358Mościbrodzka, M., Gammie, C. F., Dolence, J. C., &Shiokawa, H. 2011, The Astrophysical Journal, 735, 9,doi: 10.1088/0004-637x/735/1/9Nagai, H., Haga, T., Giovannini, G., et al. 2014, TheAstrophysical Journal, 785, 53,doi: 10.1088/0004-637x/785/1/53Nakamura, M., & Asada, K. 2013, The AstrophysicalJournal, 775, 118, doi: 10.1088/0004-637x/775/2/118Nakamura, M., Asada, K., Hada, K., et al. 2018, TheAstrophysical Journal, 868, 146,doi: 10.3847/1538-4357/aaeb2dNathanail, A., & Contopoulos, I. 2014, The AstrophysicalJournal, 788, 186, doi: 10.1088/0004-637x/788/2/186Nitta, S.-y., Takahashi, M., & Tomimatsu, A. 1991, PhysicalReview D, 44, 2295, doi: 10.1103/physrevd.44.2295Ogihara, T., Takahashi, K., & Toma, K. 2019, TheAstrophysical Journal, 877, 19,doi: 10.3847/1538-4357/ab1909Pan, Z., Yu, C., & Huang, L. 2017, The AstrophysicalJournal, 836, 193, doi: 10.3847/1538-4357/aa5c36 Parfrey, K., Giannios, D., & Beloborodov, A. M. 2015,Monthly Notices of the Royal Astronomical Society:Letters, 446, L61, doi: 10.1093/mnrasl/slu162Park, J., Hada, K., Kino, M., et al. 2019, The AstrophysicalJournal, 871, 257, doi: 10.3847/1538-4357/aaf9a9Park, J., Hada, K., Kino, M., et al. 2019, ApJ, 887, 147,doi: 10.3847/1538-4357/ab5584Piner, B. G., Pant, N., & Edwards, P. G. 2010, TheAstrophysical Journal, 723, 1150,doi: 10.1088/0004-637x/723/2/1150Piner, B. G., Pant, N., Edwards, P. G., & Wiik, K. 2008,The Astrophysical Journal, 690, L31,doi: 10.1088/0004-637x/690/1/l31Porth, O., Fendt, C., Meliani, Z., & Vaidya, B. 2011, TheAstrophysical Journal, 737, 42,doi: 10.1088/0004-637x/737/1/42Porth, O., Chatterjee, K., Narayan, R., et al. 2019, ApJS,243, 26, doi: 10.3847/1538-4365/ab29fdPu, H.-Y., Nakamura, M., Hirotani, K., et al. 2015, TheAstrophysical Journal, 801, 56,doi: 10.1088/0004-637x/801/1/56Pu, H.-Y., & Takahashi, M. 2020, The AstrophysicalJournal, 892, 37, doi: 10.3847/1538-4357/ab77abRiordan, M. O., Pe’er, A., & McKinney, J. C. 2018, TheAstrophysical Journal, 853, 44,doi: 10.3847/1538-4357/aaa0c4Sironi, L., Rowan, M. E., & Narayan, R. 2020, arXive-prints, arXiv:2009.11877.https://arxiv.org/abs/2009.11877Takahashi, H. R., Ohsuga, K., Kawashima, T., & Sekiguchi,Y. 2016, The Astrophysical Journal, 826, 23,doi: 10.3847/0004-637x/826/1/23Takahashi, K., Toma, K., Kino, M., Nakamura, M., &Hada, K. 2018, The Astrophysical Journal, 868, 82,doi: 10.3847/1538-4357/aae832Takahashi, M., Nitta, S., Tatematsu, Y., & Tomimatsu, A.1990, The Astrophysical Journal, 363, 206,doi: 10.1086/169331Takahashi, M., & Tomimatsu, A. 2008, Physical Review D,78, doi: 10.1103/physrevd.78.023012Tanabe, K., & Nagataki, S. 2008, Physical Review D, 78,doi: 10.1103/physrevd.78.024004Tanaka, S. J., & Toma, K. 2020, Monthly Notices of theRoyal Astronomical Society, 494, 338,doi: 10.1093/mnras/staa728Tchekhovskoy, A., McKinney, J. C., & Narayan, R. 2008,Monthly Notices of the Royal Astronomical Society, 388,551, doi: 10.1111/j.1365-2966.2008.13425.x—. 2009, The Astrophysical Journal, 699, 1789,doi: 10.1088/0004-637x/699/2/1789 atter density distributions of jets15