THE PARABOLIC JET STRUCTURE IN M87 AS A MAGNETOHYDRODYNAMIC NOZZLE

In order to inspect the property of the NRMHD jet acceleration in the trans-Alfvénic regime, let us follow the wave dynamics of TAWs propagating along a global poloidal magnetic field, which is originally argued by Shibata & Uchida (1985). We consider a poloidal magnetic flux tube given by

$$\begin{eqnarray} &&B_{p} \Sigma = \Psi = {\rm constant}, \end{eqnarray} \tag{ 4 }$$

where Σ is a cross section of the flux tube and Ψ is the magnetic flux. Σ is assumed to be Σ∝r². Under the WKB approximation, the conservation of the energy flux density of the TAW is considered through a circular (cross-sectional) area around the rotation axis in cylindrical coordinates (r, ϕ, z) as

$$\begin{eqnarray} &&\frac{\rho \delta V_{\phi }^2 V_{\rm A}}{B_{p}}={\rm constant}, \end{eqnarray} \tag{ 5 }$$

where ρ, δV_ϕ, and V_A ≡ B_p/(4πρ)^1/2 are the gas density, azimuthal component of the velocity field (a transverse displacement against the poloidal magnetic line of force drives a TAW), and Alfvén speed, respectively. From this equation, the following relation can be obtained:

$$\begin{eqnarray} &&\delta V_{\phi } \propto \rho ^{-1/4}. \end{eqnarray} \tag{ 6 }$$

Perturbed quantities in the transverse direction satisfy

$$\begin{eqnarray} &&\delta V_{\phi } = - \frac{\delta B_{\phi }}{\sqrt{4 \pi \rho }}, \end{eqnarray} \tag{ 7 }$$

which is another property of the shear Alfvén wave (Priest 1981). Thus, TAWs can be characterized as an equality of ratios between poloidal and toroidal components in the magnetic field,

$$\begin{eqnarray} &&\frac{\delta B_{\phi }}{B_{p}} = - \frac{\delta V_{\phi }}{V_{\rm A}}. \end{eqnarray} \tag{ 8 }$$

Note that the wave amplitudes (δB_ϕ and δV_ϕ) are not necessarily small, but can be arbitrarily finite. In a linear regime (|δB_ϕ| ≪ B_p), a perturbation of the TAW is incompressible (only magnetic tension force (B_p · ∇)B_ϕ induces the incompressible Alfvén mode), while if a wave amplitude is not negligibly small |δB_ϕ|/B_p ≳ O(1), the TAW also drives the compressional Alfvén mode. In a steady, axisymmetric MHD solution of BP82-type cold MHD outflow, |δB_ϕ| ≃ B_p is expected in the trans-Alfvénic flow V_p ≈ V_A (e.g., Königl & Pudritz 2000). On the other hand, in a non-cold MHD outflow, the initial bulk acceleration takes place by the thermal pressure gradient force and then magnetic acceleration processes take over beyond the MSP. In general, steady, axisymmetric solutions exhibit |δB_ϕ| ≳ B_p at around the Alfvénic point in both NRMHD (Kudoh & Shibata 1997; Vlahakis et al. 2000) and RMHD (Polko et al. 2013) regimes. Hereafter δB_ϕ is replaced by B_ϕ in the paper.

We therefore remark that large-amplitude, nonlinear TAWs play an important role in the trans-Alfvénic jet dynamics, which has been first examined in full two-dimensional MHD simulations (Uchida & Shibata 1985; Shibata & Uchida 1985, 1986). A magneto-centrifugal process of BP82 by quasi-linear TAWs is no longer effective, while the toroidal magnetic pressure gradient force in the poloidal direction $\nabla _{p} (B_{\phi }^2/8\pi)$ becomes dominant in the trans-Alfvénic regime. The uncoiling of B_ϕ produces a compression of the plasma in the poloidal direction, which drives fast magnetosonic waves, powering the bulk acceleration. On the other hand, the hoop stress, $-B_{\phi }^2/4 \pi r$, squeezes the plasma toward the central axis. The physical picture described above is called the "sweeping magnetic twist" mechanism by Uchida & Shibata (1986).

Considering the trans-Alfvénic NRMHD bulk flow powered by nonlinear TAWs, order-of-magnitude estimations of the one-dimensional (z) equation of motion give

$$\begin{eqnarray} &&V_{z} \sim \frac{B_{\phi }^2}{4 \pi \rho } V_{\rm A}^{-1}. \end{eqnarray} \tag{ 9 }$$

By combining Equations (6), (8), and (9), Shibata & Uchida (1985) derive

$$\begin{eqnarray} &&V_{z} \sim B_{p}^{-1} = {\rm constant}. \end{eqnarray} \tag{ 10 }$$

Thus, the terminal velocity of the TAW-driven MHD jet model with an initially uniform axial field (B_p = constant) is a constant with an order of trans-Alfvénic speed (e.g., Shibata & Uchida 1985; Kudoh et al. 1998; Kigure & Shibata 2005).

On the other hand, the terminal velocity of a TAW-driven MHD jet with a non-uniform axial field (B_p ≠ constant) can be super-Alfvénic (even super-fast magnetosonic) in the downstream. We thus suggest

$$\begin{eqnarray} &&V_{z} \propto z^{2/a}. \end{eqnarray} \tag{ 11 }$$

Again, a denotes the power-law index of the parabolic streamline z∝r^a(1 < a ⩽ 2) of a magnetic flux tube. In the specific case of a purely parabolic stream with a = 2, we note the trans-Alfvénic jet may exhibit that the flow speed increases linearly as a function of the distance V_z∝z, which can be observed in NRMHD simulations (e.g., Contopoulos 1995; Ouyed & Pudritz 1997).

In Figure 4, we show the bulk speed β(= V/c) corrected by Equation (1) (again, adopting θ_v = 14°) under an assumption of Equation (2) as a function of the de-projected distance, showing a feature of the acceleration from β ≃ 10⁻² to 0.5 in the range of z = 10²–10⁴ r_s. A distribution of β exhibits a linear increase as a function of z as β∝z around z ≃ 10²–10³ r_s, suggesting a trans-Alfvénic flow. In order to fit the data points and derive a solid picture, more sampling from this region is desired. Therefore, higher angular resolution and higher sensitivity observation with detailed analysis to identify more components and trace the proper motions are important. Beyond z ≳ 10³ r_s, it seems that the bulk acceleration is weakened slightly, indicating that the jet may transit from trans-Alfvénic to trans-fast magnetosonic, but the underlying flow is still accelerated toward a relativistic regime β ≲ 1.

Zoom In Zoom Out Reset image size

Of our particular interest is the trans-(fast) magnetosonic flow in the parabolic streamline. While the fraction of magnetic energy converted into bulk kinetic energy may still be small at the CFP, a further conversion from magnetic to kinetic energy and therefore a continued bulk acceleration would take place beyond the CFP via the so-called magnetic nozzle effect, which was originally pointed out by Camenzind (1989); poloidal magnetic flux is required to diverge sufficiently rapidly so that most of the Poynting flux in the initially magnetically dominated outflow can be converted into the kinetic energy flux beyond the CFP. This fundamentally important hypothesis has been confirmed in steady, axisymmetric solutions in both RMHD (Li et al. 1992; Li 1993; Vlahakis 2004) and NRMHD (Kudoh & Shibata 1997) regimes.

We examine a steady (∂/∂t = 0) and axisymmetric (∂/∂ϕ = 0), NRMHD outflow along the poloidal streamline in cylindrical coordinates (z, ϕ, z). The system of ideal MHD equations consists of the continuity, momentum (including the gravity), entropy, and Maxwell equations and Ohm's law (e.g., Lovelace et al. 1987):

$$\begin{eqnarray} &&\nabla \cdot (\rho {\boldsymbol V})=0, \end{eqnarray} \tag{ 12 }$$

$$\begin{eqnarray} &&\rho ({\boldsymbol V}\cdot \nabla){\boldsymbol V}=-\nabla p + \frac{1}{c}{\boldsymbol J}\times {\boldsymbol B}- \rho \nabla \phi _{\rm g}, \end{eqnarray} \tag{ 13 }$$

$$\begin{eqnarray} &&{\boldsymbol V}\cdot \nabla (p/\rho ^{\Gamma })=0, \end{eqnarray} \tag{ 14 }$$

$$\begin{eqnarray} &&\nabla \times {\boldsymbol B}=\frac{4 \pi }{c}{\boldsymbol J}, \ \nabla \times {\boldsymbol E}= 0, \ \nabla \cdot {\boldsymbol B}= 0, \end{eqnarray} \tag{ 15 }$$

$$\begin{eqnarray} &&{\boldsymbol E}+\frac{{\boldsymbol V}}{c} \times {\boldsymbol B}=0, \end{eqnarray} \tag{ 16 }$$

where Γ, ϕ_g = −GM_•/(r² + z²)^1/2, J, and E are the specific heat ratio, gravitational potential, current density, and electric field, respectively. By using Equation (12) and the solenoidal constraint in Equation (15), the velocity field can be split into poloidal and toroidal components as

$$\begin{eqnarray} &&{\boldsymbol V}= {\boldsymbol V}_{p}+V_{\phi } \hat{{\boldsymbol e}}_{\phi }, \end{eqnarray} \tag{ 17 }$$

with

$$\begin{eqnarray} &&{\boldsymbol V}_{p} = V_{r} \hat{{\boldsymbol e}}_{r} + V_{z} \hat{{\boldsymbol e}}_{z} =\frac{1}{4 \pi \rho r}\left(-\frac{\partial \Phi }{\partial z} \hat{{\boldsymbol e}}_{r} + \frac{\partial \Phi }{\partial r} \hat{{\boldsymbol e}}_{z}\right), \end{eqnarray} \tag{ 18 }$$

where Φ(r, z) is the Stokes' poloidal stream function and constant along a streamline.

We consider that the super-Alfvénic jet carries a substantial energy in the form of Poynting flux along the axial (z) direction, especially for a limit of B_ϕ ≫ B_p (in other words, the system can be described as a current-carrying jet; e.g., Contopoulos 1995). As is mentioned above, the divergence of poloidal magnetic flux is one of the key factors to attain an efficient conversion of the Poynting flux into the kinetic energy flux. Kudoh & Shibata (1997) investigate two cross-sections of a poloidal magnetic flux Σ in steady, axisymmetric jets; the solution for Σ∝r² almost fails to convert the Poynting flux into the kinetic energy flux beyond the CFP so that no further bulk acceleration occurs, with the Poynting flux remaining dominant (two-thirds of the total energy) at z/r₀ ∼ 10⁶, where r₀ denotes the radius at the footpoint of the jet. On the other hand, if the poloidal magnetic flux expands more than the purely parabolic case beyond the Alfvén radius (Σ∝r² + r³/r_A), the bulk acceleration continues until the Poynting flux decreases to zero, well beyond the CFP at around z/r₀ ≃ 10⁵.

r-dependence of Σ defines a divergence of the poloidal magnetic flux, which is physically related to the degree of the magnetic collimation of poloidal magnetic flux by the hoop stress. Three-dimensional NRMHD simulations of current-carrying jets exhibit that the distribution of the axial component of the magnetic field becomes steeper, B_z∝r⁻³ in the super-Alfvénic regime (Nakamura & Meier 2004). In general, the inwardly directed hoop stress (magnetic tension) becomes dominant in the Lorentz forces at the inner core (spine), while the outwardly directed magnetic pressure, which is dominant in the Lorentz forces at the outer edge (sheath), will be balanced with the external pressure so that the boundary shape of jets may not be magnetically self-determined (the hoop stress plays a minor role), but be controlled by the ambient ISM (Nakamura et al. 2006), as we will also discuss in Section 7.1. In late years, RMHD simulations confirm the intrinsic correlation between an efficiency of bulk accelerations and non-uniform lateral expansions of poloidal magnetic flux; it is referred to as "differential bunching" of field lines (Tchekhovskoy et al. 2009) or "differential collimation" (Komissarov 2011), but they are essentially the same as the magnetic nozzle effect (Camenzind 1989).

An edge-brightened feature is widely confirmed in the M87 jet from VLBI to Very Large Array scales (Reid et al. 1989; Owen et al. 1989; Junor et al. 1999; Kovalev et al. 2007; Ly et al. 2007; Hada et al. 2011; Asada & Nakamura 2012), indicating a hollow structure. A schematic illustration of our MHD nozzle model is shown in Figure 5. The jet expands in a parabolic shape z∝r^a(1 < a ⩽ 2) with a dominant axial component of the velocity field V_z. Our domain of examination is a hollow tube with a variable thickness δr(z), which increases as a function of z. We assume that the unit normal to the area surface is parallel to the z-axis. A concept of the differential collimation means a gradual concentration of poloidal magnetic flux toward the central part of the flow; a faster collimation of the inner flux surfaces than of the outer ones. This is mathematically expressed (e.g., Komissarov 2011) as

$$\begin{eqnarray} &&\frac{d}{d z} \left(\frac{\delta r}{r}\right) > 0, \end{eqnarray} \tag{ 19 }$$

which is equivalent to

$$\begin{eqnarray} &&\frac{1}{\delta r}\frac{d \delta r}{d z} > \frac{1}{r}\frac{d r}{d z}. \end{eqnarray} \tag{ 20 }$$

Following Contopoulos (1995; see also Liffman & Siora 1997), we apply Gauss' divergence theorem to Equation (12) with a cross-sectional area rδr:

$$\begin{eqnarray} &&\rho V_{z} r \delta r = {\rm constant} = \kappa (\Phi). \end{eqnarray} \tag{ 21 }$$

Faraday's law in Equation (15) together with Equation (16) implies the "frozen-in" field equation ∇ × (V × B) = 0. Using Stokes' theorem, we derive

$$\begin{eqnarray} &&V_{z} B_{\phi } \delta r = {\rm constant} = \chi (\Phi). \end{eqnarray} \tag{ 22 }$$

The entropy equation (14) gives the adiabat

$$\begin{eqnarray} &&p/\rho ^{\Gamma } = {\rm constant} \equiv Q (\Phi), \end{eqnarray} \tag{ 23 }$$

where above functions κ, χ, and Q of Φ are conserved quantities along a flow streamline.

Zoom In Zoom Out Reset image size

Considering Equation (13) in the one-dimensional axial (z) direction at z ≫ r, we can use the following equation:

$$\begin{eqnarray} &&\rho V_{z} \frac{d V_{z}}{d z}=-\frac{d p}{d z}-\frac{d}{d z} \left(\frac{B_{\phi }^2}{8 \pi }\right)-\rho \frac{G M_{\bullet }}{z^2}. \end{eqnarray} \tag{ 24 }$$

By differentiating Equations (21)–(23) and substituting them into Equation (24), we obtain

$$\begin{eqnarray} &&\left(V_{z}^2-V_{f}^2 \right)\frac{1}{V_z} \frac{d V_z}{d z}=\frac{V_{f}^2}{\delta r}\frac{d \delta r}{d z}+\frac{C_{s}^2}{r}\frac{d r}{d z}-\frac{G M_{\bullet }}{z^2}, \end{eqnarray} \tag{ 25 }$$

where $V_{f}^2=C_{s}^2+V_{\rm A \phi }^2$ denotes the classical fast magnetosonic speed (C_s ≡ Γp/ρ and $V_{\rm A \phi } \equiv B_{\phi }^{2}/4 \pi \rho$ are the sonic and toroidal Alfvén speeds). Note that we eliminate the second term in the right-hand side of Equation (25) because of Equation (20) and our assumption, V_f ≃ V_Aϕ in the flow:

$$\begin{eqnarray} &&\left(V_{z}^2-V_{f}^2 \right)\frac{1}{V_z}\frac{d V_z}{d z}=\frac{V_{f}^2}{\delta r}\frac{d \delta r}{d z}-\frac{G M_{\bullet }}{z^2}. \end{eqnarray} \tag{ 26 }$$

Now supposing that z∝δr^a(1 < a ⩽ 2), then our approximate Hugoniot equation for the parabolic MHD nozzle can be described as

$$\begin{eqnarray} &&\frac{z}{V_{z}}\frac{d V_{z}}{dz} \approx \frac{V_{f}^2/a-G M_{\bullet }/z}{V_{z}^2-V_{f}^2}. \end{eqnarray} \tag{ 27 }$$

Thus, for a specific case a = 2, we speculate that the poloidal jet speed exceeds the local escape speed V_esc around the trans-fast magnetosonic regime:

$$\begin{eqnarray} &&V_{z} \approx V_{f} \approx V_{\rm esc}, \end{eqnarray} \tag{ 28 }$$

where $V_{\rm esc}=\sqrt{2GM_{\bullet }/z}$. As is discussed in Kudoh & Shibata (1997), when a cross-section expands more rapidly than Σ∝r², the classical fast magnetosonic point approaches near the Alfvén point by an enhancement of bulk accelerations in the diverging magnetic nozzle effect, but the fast magnetosonic point is fairly bound to the location around $1/2 V_{p}^2 \simeq \phi _{ \rm g}$ in both cases (Σ∝r² and Σ_∞∝r³).

Note Equation (27) represents an identical formula (when a = 2) with the equation of the hydrodynamic solar wind (Parker 1958) by replacing C_s with V_f (the thermal pressure is substituted by the toroidal magnetic pressure). Thus, the nature of the trans-fast magnetosonic flow in a parabolic MHD nozzle, which is powered by TAWs, is similar to the one of transonic solution of Parker's solar wind. The hoop stress differentially bunches the poloidal magnetic flux toward the central axis so that a rapid expansion of the cross-sectional area of the poloidal magnetic flux operates on the hollow parabolic MHD stream; we speculate a de Laval-like (Dessler 1967) diverging nozzle effect in the super-fast magnetosonic regime (see also Figure 5).

The distribution of V_esc along the parabolic streamline of the M87 jet is drawn in Figure 4, showing the location where the jet is trans-escape as well as trans-fast magnetosonic. It is notable that the jet obtains a speed of ∼0.04 c in order to escape under control of the gravity of the SMBH at around 0.2–0.3 pc (≃ 6 × 10² r_s). Thus, we speculate that, as a consequence of the energy conversion from Poynting to kinetic energy fluxes, the MHD bulk acceleration takes place beyond a super-fast magnetosonic regime; the magnetic nozzle effect (Camenzind 1989) is activated by the sweeping magnetic twist mechanism (Uchida & Shibata 1986) so that a further energy conversion from Poynting to kinetic energy fluxes can be expected in the parabolic MHD nozzle as examined above.

If the jet is highly magnetized ($V_{\rm A \phi }^2 \gg C_{s}^2$) and/or a cross-section of the poloidal magnetic flux expands more rapidly than Σ∝r², then a position of the classical fast magnetosonic point approaches near the Alfvén point as mentioned above. This could be the case in the M87 jet. Thus, we may not distinguish clearly the trans-Alfvénic region from the trans-magnetosonic region in Figure 4. However, based on an intrinsic nature of MHD outflows, the accelerated flow first becomes trans-Alfvénic, and then trans-magnetosonic in the dynamical evolution. Note that the distribution of the bulk speed in Figure 4 may indicate a transition from trans-Alfvénic to trans-magnetosonic regime at ∼10³ r_s, where the speed becomes comparable to the local escape speed. A sequential MHD bulk acceleration, which can be also powered by the toroidal magnetic pressure gradient, maintains even after weakening of the acceleration at super-Alfvénic regime as we examine the magnetic nozzle effect above.

It is widely believed that AGN jets are highly supersonic from their earlier evolution stage nearby the central engine. Because the jet interior is out of causal contact with its outer edges by sonic waves in the lateral direction, a freely expanding ballistic (conical) flow is considered (Blandford & Königl). We, however, suggested that the M87 jet is not a case of the conical jet, but the parabolic jet up to over 10⁵ r_s (AN12). Furthermore, if we consider MHD jets in general, which is one of the popular theoretical frameworks to organize the jet dynamics, it may be the norm that AGN jets are still in causal contact in the lateral expansion through a communication by magnetosonic waves.

MHD jets can be a parabolic z∝r^a (1 < a ⩽ 2) if the ISM pressure is decreasing as p_ism∝z^−b with b ≃ 2 (e.g., Tchekhovskoy et al. 2008; Zakamska et al. 2008; Komissarov et al. 2009; Lyubarsky 2010). On the other hand, the jet boundary does not adjust to the ambient ISM pressure profile, but instead asymptotes to a conical shape for b > 2 (e.g., Zakamska et al. 2008; Komissarov et al. 2009; Lyubarsky 2010; Kohler & Begelman 2012). This indicates that the jet lateral expansion time becomes shorter than the signal propagation time of magnetosonic waves across the bulk flow, leading to a loss of causal contact to form a freely expanding ballistic jet. Non-conical, parabolic streamlines require a transverse pressure support by the extremal ISM (Blandford & Rees 1974; Blandford & Königl 1979). Thus, we suggest that the M87 jet may be magnetically accelerated, but thermally confined by the ISM (rather than self-confined in the MHD jet dynamics) as we will discuss below.

A power-law-like distribution of the ISM pressure with b ∼ 2 would be expected toward the ISM environment of M87 inside the Bondi accretion radius R ≲ R_B (in a spherical geometry) as a consequence of hot RIAFs, which are generally considered in low-luminosity AGNs. During the past several decades, various types of RIAF models have been examined, and we can now categorize them into three main models: the advection-dominated accretion flow (ADAF; e.g., Ichimaru 1977; Narayan & Yi 1994), the convection-dominated accretion flow (CDAF; e.g., Narayan et al. 2000; Quataert & Gruzinov 2000), and the advection-dominated inflow-outflow solution (ADIOS; e.g., Blandford & Begelman 1999, 2004).

Given the self-similar analytical treatment among these RIAF models, the gas density distribution can be described as a power-law scaling ρ_ism∝R^−k, where k = 3/2 for Bondi (1952) accretion flow (BAF) and ADAF, in which the gas undergoes a nearly free fall. On the other hand, k = 1/2 in CDAF is found because the radial convection transports the angular momentum inward, causing a reduction of the mass accretion rate and a much flatter radial density profile. The ADIOS model is known to have intermediate values 1/2 < k < 3/2 by taking into account a wind carrying away the mass, angular momentum, and energy. Assuming the polytropic equation of state, p_ism∝z^−b = R^{−(k + 1)} (we here denote z = R). Thus, 1.5 ⩽ b ⩽ 2.5 are allowed in RIAF models.

A recent analytical model for a giant ADAF (Narayan & Fabian 2011) shows b ∼ 2.3 from beyond R_B ≃ 5 × 10⁵ r_g down to the SMBH, and it is universally derived for slow rotation of the ISM (≲ 30% of the Keplerian speed) with the wider range of the α-viscosity parameter (0.001–0.3; Shakura & Sunyaev 1973). A spherically symmetric accretion with MHD turbulence (Shcherbakov 2008) produces b ∼ 2.25 without a strong dependence on the parameters and the equation of state. On the other hand, due to comprehensive studies of numerical simulations over the past decade by many authors for exploring nonlinear dynamics on RIAFs, we learn that the density profile flattens compared to the original BAF/ADAF analytical solutions, but the system does not approach a full CDAF, staying at the intermediate ADIOS, ρ_ism∝R^−k with k ≲ 1, corresponding to p_ism∝z^−b with b ≲ 2 (e.g., Yuan et al. 2012a, 2012b for recent progresses on numerical simulations of RIAFs and useful references therein).

Spatial resolutions toward M87 in current X-ray instruments may not be adequate to resolve the region (≲ 1'') inside the Bondi accretion radius of 3.8 × 10⁵ r_s (∼250 pc) although Chandra X-ray observations confirmed ISM properties, such as the King core radius r_c ≃ 1.4 kpc and the Bondi radius R_B ∼ 250 pc (see, e.g., Young et al. 2002; Di Matteo et al. 2003). We, however, remark that Chandra X-ray observations revealed the RIAF-like feature in a nearby low-luminosity AGN, NGC 3115 (Wong et al. 2011). Stellar kinematics of the bulge of NGC 3115 give the SMBH mass estimation as M_• ∼ (1–2) × 10⁹ M_☉, identifying this source to be the nearest >10⁹ M_☉ SMBH because of its proximity D ∼ 9.7 Mpc (Wong et al. 2011, and references therein). The Bondi radius (R_B ∼ 4''–5'') is spatially resolved for the first time in AGNs; the radial density profile within 4'' is described as $\rho _{\rm ism} \propto R^{-1.03^{+0.23}_{-0.21}}$, giving a nice agreement with the RIAF model. Thus, we speculate that p_ism∝z^−b with b ∼ 2 may be moderately reasonable.

We first discuss the hot accretion disk as an origin of the MHD outflow in M87. In the ADAF (RIAF) model (C_s ∼ V_K, where the Keplerian speed V_K is equal to $c/\sqrt{2 r/r_{\rm s}}$), the disk (ion) temperature is considered as nearly virial, T_i = 5.4 × 10¹²(r/r_s)⁻¹ K (Narayan & Yi 1994). Given r ∼ 10 r_s in Figure 3 (where VLBI cores are located), the sound speed C_s ≃ 0.28 c is also comparable to the local escape speed $V_{\rm esc} = c/\sqrt{r/r_{\rm s}} \sim 0.3\, c$. Therefore, we would seek much faster motions in the downstream than measured speeds 0.01 c–0.1 c at 10²–10³ r_s based on more than 10 yr of VLBA observations at 15 GHz (Kovalev et al. 2007; see also Figure 4), which are considered as the bulk speed in Section 2. If the M87 jet indeed originated from the surface of RIAFs (B_ϕ ∼ B_p; e.g., Meier 2001), as attributed to Polko et al. (2013) in Section 3, we may not expect a standard magneto-centrifugal acceleration process (BP82) in the sub-Alfvénic region, but favor an initially trans-Alfvénic flow (see also Contopoulos 1995). Thus, this may give some constraint on the required MHD jet model with an initial jet speed with an order of ∼0.1 c¹ at the jet footpoint.

Next, we discuss the Kerr black hole as an origin of the MHD outflow in M87. As we introduced in Section 3, the outermost sheath of the ergosphere-driven jet exhibits non-relativistic motion γ ≃ 1 at ∼5 × 10³ r_g as seen in GRMHD simulations (McKinney 2006); from a polar angle cross-sectional distribution of the gas density, we find that the outermost sheath has much higher mass loading (about four orders of magnitude) than the innermost jet spine. This may cause an inefficient bulk acceleration at least up to this spatial scale. A matter piling up at the outer surface of the jet may also be responsible for observed edge-brightened features in M87. The location of the CFP is about a few hundred r_s near the surface layer of the ergosphere-driven jet where the jet is still Poynting-flux dominated. Thus, we may speculate a behavior of the jet at sub- to trans-Alfvénic region z ≲ 100 r_s as a magneto-centrifugal driving process (by not a disk rotation, but an ergosphere) similar to BP82, and then the jet seems to be powered by TAWs (the toroidal magnetic pressure gradient force) at trans-Alfvénic to trans-fast magnetosonic regions. Note that a magneto-centrifugal driving allows the flow to start at a very slow speed V₀ compared to the rotational motion of the magnetic field (Contopoulos 1995). In the case of non-relativistic Keplerian disk, V_K/V₀ ≃ 300 can be obtained in steady, axisymmetric NRMHD outflow solutions (Kudoh & Shibata 1997; Vlahakis et al. 2000). Thus, even for relativistic field rotations with 0.1 c–c near the black hole, we may speculate an initial jet speed with V₀ ≪ 0.1 c at the jet footpoint.

One of the crucial tests to constrain the nature of the upstream flow within 100 r_s is to image extended jet emissions from the VLBI core. We are conducting new GMVA observations at 86 GHz in order to resolve the innermost jet structure at a few 10^{1 − 2} r_s scales.

In Section 5, we show that three additional data points of VLBI cores seem to fairly follow a parabolic streamline derived by AN12, implying that the jet maintains a single non-conical structure with nearly five orders of magnitude in distance. We here discuss about limitations of a frequency depending VLBI core shift measurement to seek the true origin of the M87 jet. The median values of the VLBI core shifts in AGN jets are reported as 0.080 mas (Sokolovsky et al. 2011) and 0.128 (Pushkarev 2012) between 15.4–8.1 GHz. A scale of ∼0.1 mas corresponds to 1–10 pc in de-projection for typical blazars with a small viewing angle θ_v ∼ 5°, which is equivalent to 10^{3 − 4} r_s with M_• = 10⁹ M_☉. Therefore, if one applies this method to blazars at even higher frequencies, an axial offset between the position of the SMBH (the origin) and/or accretion disk (the equatorial plane) and the jet footpoint may be ignorable.

On the other hand, in the case of M87, the value of the VLBI core shift between 15.4 and 8.4 GHz is similar (∼0.1 mas; Hada et al. 2011), but it corresponds to ∼10 r_s, approaching a limit to apply the method of frequency-dependent VLBI core shift measurements. By using Equation (3), we derive Δz (230 GHz) = 4.3 ± 2.2 r_s, which may be comparable to the so-called stagnation surface (SS) seen in general relativistic MHD simulations (e.g., McKinney 2006). If we consider the MHD jet launching from the rotating SMBH magnetosphere (Blandford & Znajek 1977), the SS (where the poloidal velocity component changes the sign −/+: the inflow onto the SMBH inside the SS, while the outflow as a jet outside the SS) has an offset from the SMBH. The location of the SS depends on a spin of the black hole: the balance between the gravity and the centrifugal forces along a poloidal magnetic field in the co-rotating frame with the line of force.

McKinney (2006) prescribes a black hole spin of j = 0.9375, and the SS is located at around 2–3 r_s. A lower spin will put the radius of the SS further downstream >Δz (230 GHz). Even if we consider that the MHD jet originates from an RIAF-type accretion disk H ∼ r, r ≳ 2.5 r_s (the lensed ISCO radius for the prograde spinning j = 1; Doeleman et al. 2012), then the disk surface is presumably the SS so that the situation will be unchanged. Thus, it may be challenging to explore the jet footpoint in M87 by using VLBI cores at higher frequencies. If we take an offset of the SS into consideration, three VLBI cores may not be located on a parabolic streamline in Figure 3. Thus, future sub-millimeter VLBI imaging of the jet footpoint together with the SMBH shadow is crucial.