Parabolic Jets from the Spinning Black Hole in M87

According to a steady self-similar solution of the axisymmetric FFE jet (Narayan et al. 2007; Tchekhovskoy et al. 2008, hereafter NMF07, TMN08), we consider here an approximate formula of the magnetic stream function Ψ (r, θ) in polar (r, θ) coordinates in the Boyer–Lindquist frame (Tchekhovskoy et al. 2010):

$$\begin{eqnarray}&&{\rm{\Psi }}(r,\theta )={\left(\displaystyle \frac{r}{{r}_{{\rm{H}}}}\right)}^{\kappa }(1-\cos \theta ),\end{eqnarray} \tag{ 2 }$$

where ${r}_{{\rm{H}}}={r}_{{\rm{g}}}(1+\sqrt{1-{a}^{2}})$ is the radius of the BH horizon, and the dimensionless Kerr parameter a = J/J_max describes the BH spin. J is the BH angular momentum, and its maximum value is given by J_max = r_gMc = GM²/c, where G is the gravitational constant. The ranges 0 ≤ κ ≤ 1.25 and 0 ≤ θ ≤ π/2 are adopted in TMN08. Ψ is conserved along each field (stream)line in a steady solution of the axisymmetric MHD outflow.²³ If κ = 0 is chosen, the asymptotic streamline has a split-monopole (radial) shape z ∝ R (where, $R=r\,\sin \theta$ and $z=r\,\cos \theta$), whereas if κ = 1 is chosen, the streamline has the (genuine) parabolic shape $z\propto {R}^{2}$ at R ≫ r_g (BZ77). κ = 0.75 is the case of the parabolic shape $z\propto {R}^{1.6}$ (BP82), which is important in this paper.

In magnetized RIAF simulations about a nonspinning BH (Igumenshchev 2008), the poloidal magnetic field distribution takes an "hourglass" shape and has an insignificant vertical component on the equatorial plane outside the BH. This feature becomes more prominent in the case of spinning BHs as is shown in GRMHD simulations of the ergospheric disk with a vertical magnetic flux (the Wald vacuum solution; Wald 1974); as the poloidal magnetic flux and mass accrete onto the BH, all magnetic lines threading the ergospheric disk develop a turning point in the equatorial plane, resulting in an azimuthal current sheet (Komissarov 2005).

Due to strong inertial frame dragging inside the ergosphere, all plasma entering this region is forced to rotate in the same sense as the BH. Thus, the poloidal field lines around the equatorial plane are strongly twisted along the azimuthal direction. The equatorial current sheet develops further owing to the vertical compression of the poloidal field lines caused by the Lorentz force acting toward the equatorial plane at both upper and lower (z ≷ 0) directions. Magnetic reconnection (although numerical diffusion is responsible for activating the event in an ideal-MHD simulation) will change the field topology; all poloidal field lines entering the ergosphere penetrate the event horizon. A similar result is obtained in GRFFE simulations (Komissarov & McKinney 2007).

We speculate that the situation is qualitatively unchanged even if the weakly magnetized RIAF exists in the system. Strong poloidal fields in the ergosphere compress the innermost BH accretion flow vertically and reduce the disk thickness down to H/R ≃ 0.05, while H/R ≳ 0.3 (H: the vertical scale height) is the reference value of the disk body outside the plunging region (e.g., Tchekhovskoy 2015).

Based on the physical picture above, we assume that no poloidal magnetic flux penetrates the equatorial plane at R > r_H. Therefore, the outermost field line, which is anchored to the event horizon with the maximum colatitude angle at the footpoint θ_fp = π/2, can be defined as

$$\begin{eqnarray}&&{\rm{\Psi }}({r}_{{\rm{H}}},\pi /2)=1\end{eqnarray} \tag{ 3 }$$

in Equation ((2); e.g., Tchekhovskoy et al. 2008). Figure 1 shows the outermost streamlines of Ψ (r, θ) = 1 with different κ (κ = 1: BZ77 or κ = 0.75: BP82) with a fiducial BH spin (a = 0.9375; McKinney & Gammie 2004; McKinney 2006). Let us compare the outermost streamline of the funnel jet in GRMHD simulations with Equation (3) at a quasi-steady state.

Zoom In Zoom Out Reset image size

The public version of the 2D axisymmetric GRMHD code HARM²⁴ (Gammie et al. 2003; Noble et al. 2006) is used in our examinations. The code adopts dimensionless units GM = c = 1. We, however, occasionally reintroduce factors of c for clarity. Lengths and times are given in units of ${r}_{{\rm{g}}}\equiv {GM}/{c}^{2}$ and ${t}_{{\rm{g}}}\equiv {GM}/{c}^{3}$, respectively. We absorb a factor of $\sqrt{4\pi }$ in our definition of the magnetic field. HARM implements so-called modified Kerr–Schild coordinates: x₀, x₁, x₂, x₃, where x₀ = t, x₃ = ϕ are the same as in Kerr–Schild coordinates, but the radial r(x₁) and colatitude θ(x₂) coordinates are modified (McKinney & Gammie 2004). The computational domain is axisymmetric, expanding in the r-direction from r_in = 0.98r_H to r_out = 40r_g and in the θ-direction from θ = 0 to θ = π. An "outflow" boundary condition is imposed at r = r_out; all primitive variables are projected into the ghost zones. The inner boundary condition is identical at r = r_in < r_H (no back flow into the computational domain). A reflection boundary condition is used at the poles (θ = 0, π).

Typical 2D axisymmetric GRMHD simulations (e.g., Gammie et al. 2003; McKinney & Gammie 2004; McKinney 2006) adopt a dense "Polish Doughnut"-type torus (Fishbone & Moncrief 1976; Abramowicz et al. 1978), which is in a hydrodynamic equilibrium supported by the centrifugal and gas pressure (p) gradient forces. The torus is surrounded by an insubstantial, but dynamic, accreting spherical atmosphere (the rest-mass density ρ and the internal energy density u are prescribed in power-law forms as ρ_min = 10⁻⁴(r/r_in)^−3/2 and u_min = 10⁻⁶(r/r_in)^−5/2) that interacts with the torus. This is the so-called "floor model" that forces a minimum on these quantities in the computational domain to avoid a vacuum. The initial rest-mass density ρ₀ in the system is normalized by the maximum value of the initial torus ρ_0,max on the equator.

A poloidal magnetic loop, which is described by the toroidal component of the vector potential as a function of the density ${A}_{\phi }\propto \max ({\rho }_{0}/{\rho }_{0,\max }-0.2,0)$, is embedded in the torus. The field strength is normalized with the ratio of gas to magnetic pressure (the so-called plasma-β, hereafter ${\beta }_{{\rm{p}}}\equiv 2(\gamma -1)u/{b}^{2}$, where ${b}^{2}/2={b}^{\mu }{b}_{\mu }/2$ is the magnetic pressure measured in the fluid frame). The inner edge of the torus is fixed at (r, θ) = (6r_g, π/2) and the pressure maximum is located at (r, θ) = (r_max,π/2), where r_max = 12r_g is adopted. β_p0,min = 100, where β_p0,min denotes the minimum plasma-β at t = 0, is chosen²⁵ for our fiducial run. An ideal gas equation of state p = (γ − 1)u is used, and the ratio of specific heats γ is assumed to be 4/3. If not otherwise specified, a value of a = 0.9375 is used (McKinney 2006). For further computational details, readers can refer to Gammie et al. (2003) and McKinney & Gammie (2004), which adopt default parameters in HARM.

Given small perturbations in the velocity field, the initial state of a weakly magnetized torus with a minimum value of β_p0,min ≳ 50 is unstable against the magnetorotational instability (MRI; Balbus & Hawley 1991) so that a transport of angular momentum by the MRI causes magnetized material to plunge from the inner edge of the torus into the BH. The turbulent region in the RIAF's body and corona/wind gradually expands outward. The inner edge of the torus forms a relatively thin disk with the "Keplerian" profile on the equator in both cases of nonspinning (Igumenshchev 2008) and spinning (McKinney & Gammie 2004) BHs. The turbulent inflow of the RIAF body becomes laminar at the plunging region inside the innermost stable circular orbit (ISCO; e.g., Krolik & Hawley 2002; De Villiers et al. 2003; McKinney & Gammie 2004; Reynolds & Fabian 2008), whereas there is no other specific signature in the flow dynamics across the ISCO (Tchekhovskoy & McKinney 2012). During the time evolution of the system, PFD (highly magnetized, low-density) funnel jets are formed around both polar axes (θ = 0 and π).

Our goal in using GRMHD simulations is to examine the quasi-steady structure of the boundary between the low-density funnel interior (PFD jet) and the high-density funnel exterior (corona/wind outside the RIAF body). It can be directly compared with an outermost streamline of the semianalytical FFE jet solution, which is anchored to the horizon (Figure 1), when the funnel enters a long, quasi-steady phase.

We first examine a high-resolution (fiducial) model in McKinney & Gammie (2004); default parameters in HARM (as introduced above) are adopted with a fine grid assignment ${N}_{{x}_{1}}\times {N}_{{x}_{2}}=456\times 456$. The simulation is terminated at t/t_g = 2000, or about 7.6 orbital periods at the initial pressure maximum on the equator. The MRI turbulence in the RIAF body is not sustained at t/t_g ≳ 3000 owing to our assumption of axisymmetry as explained by the anti-dynamo theorem (Cowling 1934). However, the decay of the turbulence does not affect the evolution of the PFD funnel jet (McKinney 2006). Thus, we extended the time evolution up to t/t_g = 9000 in order to examine whether the quasi-steady state of the PFD funnel jet is obtained or not. This enables us to perform a direct comparison between the steady axisymmetric FFE solution and axisymmetric GRMHD simulations regarding the funnel shape. Constraining physical quantities at the funnel edge is important for further understanding the structure of relativistic jets in M87 and others in general.

Figure 2 shows the time sequence of the distribution of the relative densities of magnetic and rest-mass energy b²/ρ. This figure provides a quantitative sense for the spatial distribution of the PFD funnel jet, wind/corona, and RIAF body (see also Figures 5 and 7 for details). Following some previous work (e.g., McKinney 2006; Dexter et al. 2012), we confirmed that the funnel jet-wind/corona boundary can be identified to be where b²/ρ ≃ 1. At the stage t/t_g = 3000 (top panel), the MRI is well developed, and thus the magnetized material in the RIAF body is swallowed by the BH. A certain amount of the poloidal flux, which falls into the ergosphere, is twisted along the azimuthal direction, and powerful PFD jets are formed toward the polar directions. Consequently, the funnel expands laterally by the magnetic pressure gradient, and its outer boundary (b²/ρ ≃ 1; orange between red and yellow on the contour map) shapes a nonconical geometry, which is quantitatively similar to the parabolic outermost streamline (thick solid line) $z\propto {R}^{1.6}$ (BP82; see Figure 1).

Zoom In Zoom Out Reset image size

After this phase, the MRI-driven turbulence in the wind/corona region above the RIAF body decays gradually. The bottom panel shows the final stage ($t/{t}_{{\rm{g}}}=9000$) of the system; the turbulent structure is saturated in the wind/corona region, but it survives near the equator in the RIAF body. It is notable that the BP82-type parabolic structure of the funnel jet-corona/wind boundary is still sustained until this phase (unchanged at the quantitative level during t/t_g ≃ 3000–9000), suggesting that the PFD funnel jet is entering a quasi-steady state, while the outside region (wind/corona and RIAF body) still evolves dynamically. The distribution of b²/ρ can be divided into the following three regions: (i) the funnel (≳1), (ii) the wind/corona (≃10⁻³–10⁻¹), and (iii) the RIAF body (≲10⁻³), respectively. Thus, we clearly identified that the PFD funnel jet (orange to red) is outlined with the BP82-type parabolic shape, rather than the genuine parabolic shape (BZ77: $z\propto {R}^{2}$). It is also notable that the boundary of the funnel follows b²/ρ ≃ 1 during the whole time of the quasi-steady state at t/t_g ≳ 3000.

Figure 3 displays the poloidal magnetic field line distribution at the same times chosen for Figure 2. We can see that the ordered, large-scale poloidal magnetic flux only exists inside the PFD funnel jet region where b²/ρ ≳ 1 (Figure 2) during the quasi-steady state t/t_g ≳ 3000. There seems to be no such coherent poloidal magnetic flux penetrating the equatorial plane at R > r_H. This is also examined in Komissarov (2005) and Komissarov & McKinney (2007). At the stage t/t_g = 3000, a lateral alignment of the poloidal magnetic flux ends around the outermost parabolic streamline $z\propto {R}^{1.6}$. This holds until the final stage of t/t_g = 9000. Thus, the distribution of poloidal magnetic field lines also indicates that the funnel interior reaches a quasi-steady state with insignificant deviation when t/t_g ≳ 3000.

Zoom In Zoom Out Reset image size

The density of contours in Figure 3 directly represents the poloidal field strength (it may be a quantitatively reasonable interpretation, at least in the funnel area). At the interior of the funnel jet, the lateral spacing of each field line decreases ($R/{r}_{{\rm{g}}}\to 0$) around the event horizon (r/r_g ≲a few), suggesting an accumulation of the poloidal flux around the polar axis (z). This is caused by the enhanced magnetic hoop stress by the toroidal field component, and it is prominent if the BH spin becomes large ($a\to 1;$ Tchekhovskoy et al. 2010). On the other hand, we may also see this effect in the downstream region (r/r_g ≲ 20) along the polar axis, but a concentration of the poloidal flux is rather smooth and weak compared with the innermost region around the event horizon. It indicates that the toroidal field does not yet fully dominate the poloidal one on this scale, and no effective bunching of the poloidal flux takes place (see Figure 9 for a visible inhomogeneity further downstream). Figure 4 confirms this quantitatively; a concentration of the poloidal magnetic field becomes clear as r increases (at r/r_g ≳ 20 and θ ≲ 20°). We can also identify the gradual decrease of $| {B}^{r}|$ as a function of θ, implying a differential bunching of the poloidal flux. Further examinations are presented in Section 3.2.2 (behaviors at the downstream in our large domain computations with different BH spins).

Zoom In Zoom Out Reset image size

No visible (but very weak) bunching of poloidal magnetic flux in the range of a few ≲ r/r_g ≲ 20 (see middle and bottom panels in Figure 3) indicates that the local poloidal field can be approximately treated as a force-free system (Narayan et al. 2009). It is worth noting that the radial (r) distribution of b²/ρ inside the funnel jet has a weak dependence on the colatitude angle (θ), as is shown in Figure 2. This also implies no effective bunching of the poloidal flux, as well as no concentration of the mass density toward the polar axis. We consider that a quasi-uniform stratification of b²/ρ inside the funnel (near the event horizon) plays a critical role in determining the dynamics of the GRMHD jet (the lateral stratification of the bulk acceleration), as is treated in Sections 3.2.2, 3.2.3, and 5.2.

In order to confirm the boundary between the PFD funnel jet and the wind/corona region at a quantitative level, we provide Figure 5. Contours of b²/ρ = 1 and β_p = 1 are shown, at the final stage of t/t_g = 9000 during the long-term, quasi-steady state (t/t_g ≳ 3000) in our fiducial run. It is similar to Figure 2 in McKinney & Gammie (2004) and Figure 3 in McKinney (2006). Note that their snapshots correspond to t/t_g ≤ 2000, which is before when we find that the funnel reaches a quasi-steady state. Notably, the equipartition of these quantities (b²/ρ = 1 and β_p = 1) is maintained at the funnel edge along the outermost BP82-type parabolic streamline. There are other contours of β_p = 1, which are distributed outside the funnel and above the equatorial plane. Quantitatively, the latter corresponds to the boundary between the wind/corona and the RIAF body (β_p = 3; e.g., McKinney & Gammie 2004; McKinney 2006).

Zoom In Zoom Out Reset image size

We also identify the boundary between the PFD funnel jet and the external area (wind/corona and RIAF body) regarding the unbound/bound state of the fluid to the BH by means of the Bernoulli parameter (e.g., De Villiers et al. 2003; McKinney & Gammie 2004): $\mathrm{Be}=-{{hu}}_{t}-1$, where $h=1+(u+p)/\rho$ is the relativistic specific enthalpy and u_t is the covariant time component of the plasma 4-velocity. In the fiducial run (and other runs as well in the next section), we confirm p/ρ ≲ 10⁻¹ throughout the computational domain. Thus, we can approximately use h ≈ 1 and $\mathrm{Be}\approx -{u}_{t}-1$, which are adopted throughout this paper. A fluid is gravitationally unbound and can escape to infinity for Be > 0, and vice versa. The contour of −u_t = 1 (Be ≈ 0) is shown in Figure 5, forming "V"-shaped geometry (originally found in McKinney & Gammie 2004). We can clearly find the bound region Be < 0 inside the PFD funnel (close to the BH) and the whole outside (throughout the computational domain r/r_g ≤ 40). It is a well-known issue, though we emphasize here that the outgoing mass flux in the PFD funnel jet does not come from the event horizon, whereas the Poynting flux is extracted from there via the Blandford–Znajek process (Blandford & Znajek 1977) as is widely examined (e.g., McKinney & Gammie 2004; McKinney 2006; Globus & Levinson 2013; Pu et al. 2015).

As investigated above, our fiducial run provides the boundary condition of the funnel edge (at a quantitative level):

$$\begin{eqnarray}&&{b}^{2}/\rho \simeq 1,\ {\beta }_{{\rm{p}}}\simeq 1,\ \mathrm{and}\ -{u}_{t}\simeq 1\ (\mathrm{Be}\approx 0)\end{eqnarray} \tag{ 4 }$$

along the outermost BP82-type parabolic ($z\propto {r}^{1.6}$) streamline (the ordered, large-scale poloidal magnetic field line), which is anchored to the event horizon with the maximum colatitude angle θ_fp ≃ π/2 (at the footpoint). There is a discrepant distribution of these quantities in previous results (McKinney & Gammie 2004; McKinney 2006). We speculate that this may be just because the funnel does not reach a quasi-steady state. Note that an alignment of the funnel edge along the specific streamline depends weakly on the initial condition (such as ${\beta }_{{\rm{p}}0,\min }$); formation of a BP82-type parabolic funnel is confirmed for a reasonable range of $50\lesssim {\beta }_{{\rm{p}}0,\min }\lt 500$ under the fixed Kerr parameter (a = 0.9) up to r_out/r_g = 100 (see Appendix A). We also note that a qualitatively similar structure of the low-density funnel edge of b²/ρ ≃ 1 , which is anchored to the event horizon with a high inclination angle θ > 80°, is identified in 3D GRMHD simulations with a = 0.5 (e.g., Ressler et al. 2017). Thus, we suggest that our finding may not depend on the dimensionality (2D or 3D).

We further examine the distribution of the total (gas + magnetic) pressure, which is underlaid in Figure 5. It gives a good sense of the jet confinement by the ambient medium. External coronal pressure outside the funnel, which consists of both gas and magnetic contributions (β_p ≃ 1), is surely overpressured with respect to the magnetic-pressure-dominated region (β_p ≪ 1) inside the funnel. The situation seems unchanged even after the MRI-driven MHD turbulence saturates as is shown in Figure 6. In the vicinity of the BH (r/r_g ≲ 5), the total pressure of the RIAF body is dominated by the gas pressure (t/t_g ≲ 3000), while it is in an equipartition (β_p ∼ 1; t/t_g ≳ 3000). The important point of Figure 6 is that the funnel pressure (magnetically dominated) and the coronal pressure (β_p ≃ 1) are almost unchanged during t/t_g ≃ 3000–9000, suggesting that the quasi-steady parabolic shape of the funnel is sustained. Some convenient terminology is provided in Figure 5 for dividing the domain into the funnel (b²/ρ ≫ 1, β_p ≪ 1), corona (${b}^{2}/\rho \simeq {10}^{-2}$, β_p ≃ 1), and RIAF body (b²/ρ ≪ 10⁻², β_p ≳ 1), following the literature (De Villiers et al. 2003; McKinney & Gammie 2004; Hawley & Krolik 2006; Sadowski et al. 2013).

Zoom In Zoom Out Reset image size

We realize, however, that the ram pressure of the accreting gas becomes even higher than the total pressure in the RIAF body on the equatorial plane (by almost an order of magnitude). This is conceptually similar to the so-called "magnetically arrested disk" (MAD; e.g., Narayan et al. 2003; Igumenshchev 2008; Tchekhovskoy et al. 2011), although the accretion flow in our HARM fiducial run is identified as the "standard and normal evolution" (SANE; Narayan et al. 2012) without having an arrested poloidal magnetic flux on the equatorial plane at $R\gt {r}_{{\rm{H}}}$ (see Figure 3). Anyhow, as a consequence of this combination of effects, we could expect the PFD jet to be deformed into a nonconical geometry. Note that the funnel structure becomes radial (i.e., conical) if the magnetic pressure in the funnel is in equilibrium with the total pressure in the corona and the RIAF body in 3D GRMHD simulations (e.g., Hirose et al. 2004). We also remark that the low-density funnel edge with b²/ρ ≃ 1 is established even in the MAD state with 3D runs (Ressler et al. 2017).

Finally, Figure 7 provides further examination of outflows (and inflows as well) in the system. The outgoing radial mass flux density exists both inside and outside the funnel, but that in the corona is significantly higher than the funnel with ∼1–3 orders of magnitude, which is quantitatively consistent with other 3D simulations in the SANE state (e.g., Sadowski et al. 2013; Yuan et al. 2015). Aside from the terminology in McKinney & Gammie (2004), we consider outflows in the funnel as jets and those in the corona as winds (e.g., Sadowski et al. 2013), as is labeled. The former is highly magnetized with considerable poloidal magnetic flux and powered by the spinning BH, but the latter is not (see Figures 2 and 3). A division boundary of these outflows lies on the Be ≈ 0 contour along the outermost BP82-type parabolic streamline. That is, the BH-driven PFD funnel jet is unbound, but the RIAF-driven coronal wind is bound (at least in our simulation up to r_out/r_g = 100).

Zoom In Zoom Out Reset image size

In the coronal wind (bound), a considerable mass is supplied from the RIAF, and thus the outflow does not possess sufficient energy to overcome the gravitational potential. On the other hand, the funnel jet (unbound), which carries very little mass, becomes relativistic quite easily (due to the magnetic acceleration along the poloidal magnetic field) and could escape to infinity. However, the bound wind may not exist permanently at large distance; Be is presumably not a constant (in the nonsteady flow), and the sign can be positive with a 3D turbulent environment (Yuan et al. 2015). The jet stagnation surface in the PFD funnel, at which the contravariant radial component of the plasma 4-velocity becomes zero (${u}^{r}=0$; McKinney 2006; Broderick & Tchekhovskoy 2015; Pu et al. 2015), is clearly identified.

In Figure 7, we can see that the jet stagnation surface does reasonably coincide with the bound/unbound boundary: $\mathrm{Be}\approx 0\ (-{u}_{t}=1)$, except for the polar region (θ ∼ 0) of the V-shaped geometry. Outside the PFD funnel jet, the coronal wind can be identified with the outgoing radial mass flux density, but there is an accreting gas (identified by the ingoing radial mass flux density) around θ ≃ π/4. On the other hand, there is the outgoing gas in the RIAF body at π/3 ≲ θ ≲ 2/3π (see Figure 5). This is not a wind, but is associated with turbulent motions. Thus, both inflows and outflows are mixed up in the corona and RIAF body, suggesting the adiabatic inflow-outflow solution (ADIOS; Blandford & Begelman 1999, 2004). Note that u^r does not vanish along the jet/wind boundary, and thus the wind plays a dynamical role in confining the PFD jet (see also Figure 11).

The origin of the wind is beyond the scope of our consideration here, but the magnetocentrifugal mechanism (BP82) would be unlikely to operate; this is because of the absence of a coherent poloidal magnet field outside the PFD funnel (see Figure 3), where the toroidal magnetic field is dominant and the plasma is not highly magnetized (b²/ρ ≪ 1 and β_p ≳ 1; see Figures 2 and 5). Note that a dominant toroidal magnetic field may also be true in the SANE state with 3D runs (e.g., Hirose et al. 2004). The funnel-wall jet (Hawley & Krolik 2006), which could be driven by a high total pressure corona squeezing material against an inner centrifugal wall, does not seem to appear in our fiducial simulation. We do not find any significant evidence of it: there is no pileup of the mass flux and/or the total pressure at the funnel edge along the BP82-type parabolic outermost streamline. There is a key difference between the coronal wind and the funnel-wall jet: the former is bound, while the later is unbound at least in the vicinity of the BH (De Villiers et al. 2003). The outer boundary of the matter-dominated coronal wind (b²/ρ < 0.1; see Figures 2 and 7) is somewhat indistinct as is indicated by Hawley & Krolik (2006).

Based on our fiducial run, we further examine the BP82-type parabolic structure ($z\propto {R}^{1.6}$) of the PFD funnel jet against the varying BH spin. Different Kerr parameters are examined with the same value of β_p0,min = 100 in the extended computational domain r_out/r_g = 100 (with a grid assignment ${N}_{{x}_{1}}\times {N}_{{x}_{2}}=512\times 512$). We prescribe r_max/r_g = 12.95, 12.45, 12.05, and 11.95 for a = 0.5, 0.7, 0.9, and 0.99, respectively.

3.2.1. Funnel Structure

Figure 8 exhibits b²/ρ at the final stage t/t_g = 12,000 for various BH spins. First of all, we confirmed that the overall structure outside the funnel seems to be unchanged with different spins: b²/ρ ≃ 10⁻³–10⁻¹ is obtained in the wind/corona region, while the RIAF body sustains b²/ρ ≲ 10⁻³. The MRI-driven turbulence has decayed in both the wind/corona and RIAF body by this phase. The funnel boundary fairly follows the outermost BP82-type parabolic streamline, which is anchored to r = r_H, but shifts inward with increasing a (the funnel becomes "slimmer"). We find that b²/ρ ≃ 1 is maintained along the funnel edge for sufficiently high spins (a ≥ 0.9), while the value is somewhat smaller—around b²/ρ ≃ 0.5—for lower spins (a ≤ 0.7) in the quasi-steady phase; b²/ρ ≤ 5 is obtained at the jet stagnation surface (see Figure 14 in the next section), so that the downstream of the funnel jet does not hold a highly PFD state.

Zoom In Zoom Out Reset image size

We report one obvious but notable dependence on the BH spin a during the time evolution of the system in Figure 8; b²/ρ in the funnel grows larger when a increases. This implies that a large spin would be suitable for obtaining a large value of b²/ρ at the jet-launching region, which determines the maximum Lorentz factor (e.g., McKinney 2006; see also Figure 14 for reference values). Thus, we may consider the BH spin dependence on the asymptotic Lorentz factor at large distance, although another factor (i.e., the magnetic nozzle effect) may also affect this. Beneath the jet stagnation surface, there is an inflow u^r < 0 (Be < 0) along each streamline because the fluid is strongly bound by the hole's gravity. Thus, a shift of the jet stagnation can be expected, depending on the black hope spin; the fluid can escape from a deeper gravitational potential well of the hole owing to an enhanced magnetocentrifugal effect in a cold GRMHD flow (Takahashi et al. 1990; Pu et al. 2015) if a becomes large. See Section 3.2.3 (with Figure 14) for more details.

On the other hand, we confirmed that the shape of the funnel exhibits a weak dependence on the BH spin (a ≥ 0.5). As a increases, the angular frequency Ω of a poloidal magnetic field line (Ferraro 1937) increases where frame dragging is so large inside the ergosphere that it generates a considerably larger toroidal magnetic field. Consequently, a magnetic tension force due to hoop stresses ("magnetic pinch") would be expected to act more effectively on collimating the funnel jet. This, however, is not the case at the funnel edge; hoop stresses nearly cancel centrifugal forces, suggesting that a self-collimation is negligible (McKinney & Narayan 2007a, 2007b). Note that a self-collimation will be effective at the funnel interior as a increases (see below).

As is introduced in Section 2.1, the outermost BP82-type parabolic streamline ($z\propto {R}^{1.6}$) of the steady axisymmetric FFE jet solution (NMF07, TMN08), which is anchored to the event horizon with the maximum colatitude angle θ_fp = π/2, can suitably represent the boundary (Equation (3)) between the PFD funnel jet and the wind/corona region. The funnel edge can be approximately defined as the location where b²/ρ ≃ 0.5–1 along the specific curvilinear shape. This is, at least in our GRMHD simulations, valid on the scale of r/r_g ≲ 100 with a range of Kerr parameters a = 0.5–0.99 (due to the limited space, we omit the case of a = 0.998, but we confirmed similar structure of the PFD funnel jet as is shown in Figure 8). Again, we report a formation of the BP82-type parabolic funnel with 50 ≲ ${\beta }_{{\rm{p}}0,\min }\lt 500$ (a = 0.9) up to r_out/r_g = 100 in our GRMHD runs (see Appendix A).

Figure 9 displays the poloidal magnetic field lines, corresponding to each panel of Figure 8. We obtained qualitatively similar results to Figure 3 at the final stage (t/t_g = 9000), and the overall feature is unchanged with varying BH spin as is similar to Figure 8. Again, there is no coherent structure of the poloidal magnetic field lines in both the corona and RIAF body for all cases: a = 0.5–0.99. It is likely that the MHD turbulence (due to the MRI) saturates in our axisymmetric simulations, whereas there are some features of poloidal field tangling. As is widely examined in a 3D simulation with similar initial weak poloidal field loops lying inside the torus, the evolved field is primarily toroidal in the RIAF body and corona (e.g., Hirose et al. 2004); the magnetic energy of the coronal field is, on average, in equipartition with the thermal energy (β_p ≃ 1), and the corona does not become a magnetically dominated force-free state, which is quantitatively consistent with our results.

Zoom In Zoom Out Reset image size

By contrast, in the funnel region (b²/ρ ≳ 1) the field is essentially "helical" with a dominant radial component outside the ergosphere and a toroidal component that becomes larger with increasing BH spin (see also Hirose et al. 2004). In their result, the magnetic pressure in the funnel is in equilibrium with the total pressure in the corona and the inner part of the RIAF body. This is, however, different from our results; as is shown in our fiducial run (Figure 5), the total pressure in the PFD funnel (i.e., the dominant magnetic pressure) is smaller than the outer total pressure (especially at r/r_g ≲ 20–30, where the curvature of the funnel edge is large). We understand that this is one of the reasons why our PFD funnel jet is deformed into a paraboloidal configuration, rather than maintaining the radial shape.

Figure 10 shows that an excess of the external coronal pressure (compared with the funnel pressure) gradually decreases as r increases (as the curvature of the parabolic funnel becomes small). They are almost comparable at r/r_g = 100, indicating that the pressure balance is established, where the magnetic pressure inside the funnel decays by about two orders of magnitude compared with that at r/r_g = 10. By combining with Figure 11, it is implied that a coronal wind plays a dynamical role in shaping the jet into a parabolic geometry.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

An inhomogeneous spacing of poloidal magnetic field lines in the lateral direction (R) at z/r_g ≳ 50 (the tendency becomes strong with a ≥ 0.7) is confirmed in Figure 9 (the density of contours is high around the polar axis, while it becomes low near the funnel edge). This is widely seen in the literature and interpreted as the self-collimation by the magnetic hoop stress, which collimates the inner part of streamlines relative to the outer part (e.g., Nakamura et al. 2006; Komissarov et al. 2007, 2009; Tchekhovskoy et al. 2008, 2009, 2010). This is named the "differential bunching/collimation" of the poloidal magnetic flux (Komissarov et al. 2009; Tchekhovskoy et al. 2009), leading to a sufficient opening of neighboring streamlines for an effective bulk acceleration via the magnetic nozzle effect (e.g., Toma & Takahara 2013, for recent progress and references therein). In the next section, further quantitative examination of the differential bunching of the poloidal magnetic flux is presented (see also Figure 13), which is associated with the jet bulk acceleration.

3.2.2. Outflows

Next, we examine the nature of the coronal wind in Figure 11, corresponding to each panel of Figure 8. Again, we obtain qualitatively similar results to Figure 7 (the magnitude of the outgoing radial mass flux density of the coronal wind is higher than that of the funnel jet by ∼1–3 orders of magnitude) even as the BH spin varies. This is consistent with the SANE state of 3D GRMHD simulations, in which the magnitude of the outgoing radial mass flux density of the coronal wind does not exhibit a noticeable dependence on the BH spin (Sadowski et al. 2013). We also confirm that the distribution of $\mathrm{Be}\approx 0\ (-{u}_{t}=1)$ forms a V-shaped geometry as mentioned in Section 3.1. The right side of the V-shaped distribution of Be ≈ 0 reasonably follows the outermost BP82-type parabolic streamline, which is anchored to the event horizon, up to r_out/r_g = 100 for a = 0.5–0.99 (although there is some deviation near the outer radial boundary in a = 0.9), suggesting that the coronal wind is bound up to r_out/r_g ∼ 100. It is notable that the jet stagnation surface (${u}^{r}=0$) inside the funnel shifts toward the BH as a increases (e.g., Takahashi et al. 1990; Pu et al. 2015), but it coincides with the left side of the V-shaped distribution of Be ≈ 0 except for the polar region (θ ∼ 0; see also Figure 7).

Finally, we examine the velocity of the funnel jet. We evaluate the Lorentz factor Γ, which is measured in Boyer–Lindquist coordinates (e.g., McKinney 2006; Pu et al. 2015):

$$\begin{eqnarray}&&{\rm{\Gamma }}=\sqrt{-{g}_{{tt}}}{u}^{t},\end{eqnarray} \tag{ 5 }$$

where

$$\begin{eqnarray*}\begin{array}{rcl}{g}_{{tt}} & \equiv & -({\rm{\Delta }}-{a}^{2}{\sin }^{2}\theta )/{\rm{\Sigma }},\\ {\rm{\Delta }} & \equiv & {r}^{2}-2r+{a}^{2},{\rm{\Sigma }}\equiv {r}^{2}+{a}^{2}{\cos }^{2}\theta .\end{array}\end{eqnarray*}$$

We report another visible dependence on the BH spin a during the time evolution of the system in Figure 12, which exhibits Γ at the final stage t/t_g = 12,000, corresponding to each panel of Figure 8. It is notable that a high value of Γ (≲1.5 in a = 0.7, ≲2.3 in a = 0.9, and ≲2.8 in a = 0.99) is distributed between two outermost streamlines ($z\propto {R}^{2}$ and $z\propto {R}^{1.6}$) of the semianalytical FFE jet, which are anchored to the horizon (Figure 1).

Zoom In Zoom Out Reset image size

Qualitatively similar results (Γ becomes large at the outer layer in the funnel, while Γ ≃ 1 is sustained at the inner layer) are obtained in SRMHD simulations with a fixed jet boundary and a solid-body rotation at the jet inlet, which provides a good approximation of the behavior of field lines that thread the horizon (Komissarov et al. 2007, 2009). Furthermore, the trend can be seen in FFE (a = 1; Tchekhovskoy et al. 2008), GRFFE (a = 0.9–0.99; Tchekhovskoy et al. 2010), and GRMHD (a = 0.7–0.98; Penna et al. 2013) simulations with modified initial field geometries in the torus (Narayan et al. 2012). Note that the inner "cylinder-like" layer with a lower Γ < 1.5 for all spins arises from a suppression of the magnetic nozzle effect; an effective bunching poloidal flux takes place near the polar axis owing to the hoop stresses (see Figure 9), so that enough separation between neighboring poloidal field lines is suppressed (see also Komissarov et al. 2007, 2009).

Another notable feature is the inhomogeneous distribution of Γ (blobs) at the outer layer in the funnel; throughout the time evolution, knotty structures are formed with Γ ≳ 2 when a ≥ 0.9 (two right panels in Figure 12). Blobs appear around the funnel edge R/r_g ≳ 10 (z/r_g ≳ 20) with Γ ≳ 1.5, which could be observed as a superluminal motion ≳c with a viewing angle ≲25°. Note that no proper motion has been detected within a scale of 100r_g in deprojection in M87, where we interpret VLBI cores as an innermost jet emission. We discuss the formation of blobs and their observational counterparts in Section 5.3.2.

An evolution of the Lorentz factor in a highly magnetized MHD/FFE outflow can be expressed with the approximated formula (in the so-called "linear acceleration" regime²⁶ ):

$$\begin{eqnarray}&&{\rm{\Gamma }}\approx \sqrt{1+{(R{\rm{\Omega }})}^{2}}\approx R{\rm{\Omega }}\propto {z}^{1/\epsilon },\end{eqnarray} \tag{ 6 }$$

where

$$\begin{eqnarray}&&{\rm{\Omega }}\approx 0.5{{\rm{\Omega }}}_{{\rm{H}}}(a).\end{eqnarray} \tag{ 7 }$$

Ω_H denotes the angular frequency of the BH event horizon as

$$\begin{eqnarray}&&{{\rm{\Omega }}}_{{\rm{H}}}(a)=\displaystyle \frac{{ac}}{2{r}_{{\rm{H}}}}.\end{eqnarray} \tag{ 8 }$$

This formula is valid at a moderately relativistic regime Γ ≳ 2 (Komissarov et al. 2007, 2009; McKinney & Narayan 2007b, NMF07, TMN08) and/or all the way out to large distances (the jet radius is large enough: RΩ ≫ 1, where the curvature of the magnetic surface is unimportant; Beskin & Nokhrina 2006, 2009). This can be expected at $z/{r}_{{\rm{g}}}\gtrsim 100$ as is exhibited in the steady axisymmetric GRMHD solution (Pu et al. 2015). Numerical simulations of the FFE jets in both genuine parabolic and the BP82-type parabolic shapes provide a quantitative agreement with the analytical solution (TMN08).

Based on our result shown in Figure 12, at least up to r/r_g = 100, we could not find clear evidence of linear acceleration in the funnel jet. This is presumably due to a lack of a differential bunching/collimation of the poloidal magnetic flux as is expected in a highly magnetized MHD/FFE regime (e.g., Tchekhovskoy et al. 2009) during the lateral extension of the funnel (in the regime of RΩ ≫ 1). However, there are visible increases of distributed Γ in the funnel jet as a increases. We can identify a physical reason in Figure 13: a concentration of the poloidal magnetic field becomes strong as r increases, but it is also enhanced as a function of a ($a=0.5\to 0.99$). As a consequence of the differential operation of the magnetic nozzle effect, a larger Γ value in the funnel jet is obtained in the higher spin case. Note that blobs (knotty structures) also appear near the funnel boundary at a ≥ 0.9, and we discuss this feature in Section 5.3.2.

Zoom In Zoom Out Reset image size

Contours of the equipartition β_p = 1 are also displayed in Figure 12. We can see that one of the contours is elongated near the outermost BP82-type parabolic streamline. Thus, our boundary condition of the funnel edge (Equation (4)) is moderately sustained up to r_out/r_g = 100 (there is a departure of β_p = 1 from the funnel edge in a = 0.9, but this seems to be a temporal and/or boundary issue). Figure 12 also provides a clue of the velocity in the corona/wind region. At the funnel exterior, no coherent poloidal magnetic flux exists in the quasi-steady SANE state (Figure 9). The weakly magnetized (${b}^{2}/\rho \ll 1,\ {\beta }_{{\rm{p}}}\simeq 1$) coronal wind carries a substantial mass flux compared with the funnel jet (Figure 11; see also Sadowski et al. 2013; Yuan et al. 2015). Therefore, it is unlikely that the coronal wind could obtain a relativistic velocity; Γ = 1 is sustained in all cases with a = 0.5–0.99 (see also Yuan et al. 2015, with a = 0). Note that similar results are confirmed even in 3D simulations of the MAD state (e.g., Penna et al. 2013), in which the coherent poloidal magnetic flux is presumably arrested on the equatorial plane at R > r_H (Tchekhovskoy 2015).

With standard parameters in HARM, our GRMHD simulations provide several interesting results, which does not depend on the BH spin when a = 0.5–0.99. We, however, need further investigations to find features such as a linear acceleration of the underlying flow with an extended computational domain (r_out/r_g > 100). Also, one of the important issues is to investigate whether an unbound wind could surely exist on large scales. If so, it indicates that $\mathrm{Be}\simeq 0\ (-{u}_{t}\simeq 1)$ will not hold at the jet/wind boundary. How could the BP82-type parabolic funnel jet be maintained by an unbound wind and/or other external medium? How are the equipartition conditions ${b}^{2}/\rho \simeq 0.5$–1 and β_p ≃ 1 maintained (or modified)? These questions will be addressed in a forthcoming paper.

3.2.3. Jet Stagnation Surface

We present Figure 14 to examine the jet stagnation surface and the local value of b²/ρ, with respect to the BH spin (a). As is also shown in Figure 11, the jet stagnation surface (${u}^{r}=0$ inside the funnel) shifts toward the BH if a increases (see, e.g., Takahashi et al. 1990, for an analytical examination in the Kerr spacetime) owing to an increase of Ω (the outflow can be initiated at the inner side), but qualitatively similar structures of the surface are obtained (a = 0.5–0.99) as is clearly seen in Figure 14.

Zoom In Zoom Out Reset image size

Coherent poloidal magnetic field lines are regularly distributed inside the funnel edge along the outermost parabolic streamline (BP82: $z\propto {r}^{1.6}$), which is anchored to the event horizon. As the BH spin increases, the density of contours becomes high, indicating that the poloidal field strength goes up. Especially, as is examined in Figure 3, it is prominent around the polar axis (z) in the very vicinity of the event horizon ($r/{r}_{{\rm{g}}}\lesssim$ a few) as $a\to 1$ (e.g., Tchekhovskoy et al. 2010). The closest part (to the BH) of the jet stagnation is located around the funnel edge with a = 0.5–0.99. During the quasi-steady state, the jet stagnation surface is almost stationary in our GRMHD simulations. Outer edges of the stagnation surface in funnel jets (r/r_g ≃ 5–10 in the case of a = 0.5–0.99) give an (initial) jet half opening angle of ∼40°–50° in deprojection (see Figures 13 and 14).

Contours of b²/ρ with selected values are also displayed in each panel of Figure 14 for our reference. In the vicinity of the BH at r/r_g ≲ 20, the value of b²/ρ decreases monotonically (approximately independent of the colatitude angle θ inside the funnel) as r increases. Again, this can be interpreted as a consequence of no visible (but very weak) bunching of poloidal magnetic flux at r/r_g ≃ 5–10 (around the stagnation surface edges in the case of a = 0.5–0.99), as is shown in Figures 4, 13, and 14. Thus, we may not expect a significant concentration of the mass density toward the polar axis in the range of a few ≲ r/r_g ≲ 20 as examined with our fiducial run in Section 3.1 (see also Figures 2 and 3). Depending on the BH spin (a = 0.5, 0.7, 0.9, and 0.99), b²/ρ ≃ 2, 5, 10, and 20 are identified around the closest part (near the funnel edge) on the jet stagnation surface. This is located between two outermost streamlines (z ∝ r² and $z\propto {r}^{1.6}$) of the semianalytical solution of the FFE jet.

As is seen in Figure 12, a high value of Γ is distributed throughout an outer layer of the funnel between two outermost streamlines ($z\propto {R}^{2}$ and $z\propto {R}^{1.6}$), which are anchored to the event horizon. Having a high value of b²/ρ at the jet-launching point is suggestive that the flow will undergo bulk acceleration to relativistic velocities, as seen in Equation (1). This is a necessary but not sufficient condition, as the magnetic nozzle effect is also needed, which can be triggered by a differential bunching of poloidal flux toward the polar axis. As is suggested in Takahashi et al. (1990) and Pu et al. (2015), the location of the jet stagnation surface, where the magnetocentrifugal force is balanced by the gravity of the BH, is independent of the flow property, such as the rate of mass loading, because it is solely determined by a and Ω (Equation (7)). We point out that a departure of the jet stagnation surface from the BH at a higher colatitude ($\theta \to 0$) gives a prospective reason for the lateral stratification of Γ at large distances (where the sufficient condition for the bulk acceleration may be applied).

The above issue could be associated with the so-called limb-brightened feature in the M87 jet. Note that we identified the value of ${b}^{2}/\rho \simeq 0.5$–1 as the physical boundary at the funnel edge along the outermost BP82-type parabolic streamline, which is anchored to the event horizon. Thus, if this boundary condition holds even further downstream (r/r_g ≫ 100), the funnel edge is unlikely to be a site exhibiting a relativistic flow, as is examined in Figure 12. A highly Doppler-boosted emission may not be expected there, but an alternative process, such as the in situ particle acceleration, may be considered at the edge of the jet sheath (as a boundary shear layer) under the energy equipartition between the relativistic particles and the magnetic field (e.g., Stawarz & Ostrowski 2002). On the other hand, limbs in the M87 jet have a finite width δR inside their edges, and δR seems to increase in the downstream direction (e.g., Asada et al. 2016), suggesting a differential bunching of streamlines (e.g., Komissarov 2011).

In this paper, we identify the outer jet structure (limbs) as the jet sheath, while the inner jet structure (inside limbs) is identified as the jet spine. In the next section, our results are compared with VLBI observations, followed by related discussions in Section 5. In particular, we assign Section 5.2 for discussion of a limb-brightened feature in the context of MHD jets and Section 5.3.2 to describe the origin of knotty structures.

Figure 15 overviews the geometry of the M87 jet by compiling the data in the literature (see the caption for references). Multiwavelength observations by Asada & Nakamura (2012, hereafter AN12) revealed that the global structure of the jet sheath is characterized by the parabolic stream $z\propto {R}^{1.73}$ at z/r_g ∼ (400–4) × 10⁵ (see also Hada et al. 2013), while it changes into the conical stream $z\propto {R}^{0.96}$ beyond the Bondi radius of r_B/r_g ≃ 6.9 × 10⁵ (∼205 pc; Rafferty et al. 2006). Hada et al. (2013) and Nakamura & Asada (2013) examined the innermost jet region (z/r_g ≳ 10) by utilizing the VLBI core shift (Hada et al. 2011). Hada et al. (2013) suggest a possible structural change toward upstream at z/r_g ∼ 300, where the Very Long Baseline Array (VLBA) core at 5 GHz is located. The innermost jet sheath (z/r_g ≳ 200) is recently revealed with HSA 86 GHz (Hada et al. 2016). Based on our theoretical examinations presented in previous sections, we overlay the outermost streamlines of the semianalytical FFE jet model (NMF07, TMN08) with varying Kerr parameters (a = 0.5–0.99) on data points in Figure 15 for comparison.

Zoom In Zoom Out Reset image size

There are notable findings: (i) the inner jet radius (at z/r_g ≲ 100), which is represented by VLBI cores at 15–230 GHz, is traced by either outer parabolic or inner genuine parabolic streamlines of FFE jets, which are anchored to r = r_H with θ_fp = π/2. Within the uncertainties, we cannot distinguish the shape of the streamline, but there is a tendency that the mean values shift toward the genuine parabolic streamlines inside the funnel. Therefore, we consider that the millimeter VLBI core at 230 GHz with EHT observations (Doeleman et al. 2012; Akiyama et al. 2015; hereafter the EHT core) presumably shows the innermost jet to be in a highly magnetized (PFD) regime with b²/ρ ≫ 1 and Γ ≲ 1.5 (see Figures 8 and 12) inside the funnel. Note that the dominant magnetic energy of the VLBI core at 230 GHz is originally proposed by Kino et al. (2015). This may reflect the jet spine, however, rather than the jet sheath at the funnel edge (see Section 5.3.1 for discussions).

(ii) At the scale of 100 ≲ z/r_g ≲ 10⁴, we identify a clear coincidence between the radius of the jet sheath and the outermost BP82-type streamline of the FFE jet solution. We also confirm a reasonable overlap between the VLBI cores at 5 and 8 GHz and an extended emission of the jet sheath at VLBA 43 and HSA 86 GHz (see also Hada et al. 2013). Therefore, the hypothesis of the VLBI core as the innermost jet emission (Blandford & Königl 1979) is presumably correct at this scale, although a highly magnetized state of VLBI cores, suggested by Kino et al. (2014, 2015), and the frequency (ν) dependent VLBI core shift ${\rm{\Delta }}z(\nu )\propto 1/\nu$ taking place in the nonconical jet geometry in M87 (Hada et al. 2011) may conflict with original ideas in Blandford & Königl (1979), where an equipartition between the magnetic and synchrotron particle energy densities and a constant opening angle and constant velocity jet is considered. We also remind readers about our recent result on the jet geometry of blazars that examined VLBI cores (Algaba et al. 2017), suggesting that nonconical structures may exist inside the SOI r_SOI ∼ (10⁵–10⁶)r_g.

(iii) At z/r_g ≃ 10⁴–10⁵, it is visible that data points (the radius of the jet sheath) start to deviate slightly from $z\propto {R}^{1.6}$, but a parabolic shape is sustained. This may indicate a new establishment of the lateral force equilibrium between the funnel edge and the outer medium (wind/corona above the RIAF), or the jet sheath starts to be Doppler deboosted (see Figures 16 and 18, as well as our discussion in Section 5.2). Previous GRMHD simulations exhibit a conical shape of the funnel edge at z/r_g ≳ few × 10² (e.g., McKinney 2006), implying that the jet is overpressured against the outer medium. This, however, could not be the case in M87. The intrinsic half opening angle (θ_j) of the jet sheath attains the level of θ_j ≃ 05 at $z/{r}_{{\rm{g}}}\simeq 4\times {10}^{5}$.

Zoom In Zoom Out Reset image size

(iv) Data points are clearly deviated from the outermost BP82-type parabolic streamlines of the FFE jet solution beyond r_B. As is originally suggested by Cheung et al. (2007), a structured complex known as "HST-1" (Biretta et al. 1999) is located just downstream of r_B at ~10⁶r_g. AN12 suggest a geometrical transition of the M87 jet as a consequence of the overcollimation of the highly magnetized jet at the HST-1 complex, which can be initiated by forming the HST-1 complex, at ≃r_B. The jet exhibits the conical geometry with θ_j ≃ 05 (const) at the kiloparsec scale (z/r_g ≳ 3 × 10⁶), while θ_j ≃ 01 is obtained at HST-1. We can refer to another sample of the AGN jet structural transition at  ≃r_SOI in NGC 6251 (Tseng et al. 2016; see also Section 5.1).

Note that the magnetic pressures at HST-1 and several knots in the downstream (≃10⁻⁹–10⁻⁸ dyn cm⁻²; Owen et al. 1989; Perlman et al. 2003; Harris et al. 2009) are highly overpressured against the interstellar medium (ISM) pressure of ≃10⁻¹¹–10⁻¹⁰ dyn cm⁻² (Matsushita et al. 2002). Thus, the lateral pressure equilibrium between the conical jet sheath and the ambient medium does not seem to be sustained beyond r_B. The inner part of highly magnetized jets can be heavily overpressured with respect to the outer part owing to the hoop stress as is examined in numerical simulations and self-similar steady solutions (Nakamura et al. 2006; Zakamska et al. 2008). A conical expansion of the highly magnetized (with a dominant toroidal field component), overpressured jet sheath against the uniform ISM environment is reproduced in numerical simulations (e.g., Clarke et al. 1986).

As a summary of this section, we conclude that the edge of the jet sheath in M87 upstream of r_B can be approximately described as the outermost BP82-type streamline of the FFE jet solution with the Kerr parameter a > 0, which is anchored to the event horizon. Thus, we suggest that the parabolic jet sheath in M87 is likely powered by the spinning BH. Recent theoretical arguments clarified that the outward Poynting flux is generally nonzero (i.e., the BZ77 process generally works) along open magnetic field lines threading the ergosphere (Komissarov 2004; Toma & Takahara 2014). Thus, our findings support the existence of the ergosphere. We note, however, that there is an alternative suggestion that the jet sheath is launched in the inner part of the Keplerian disk at R ∼ 10r_g (Mertens et al. 2016).

Figure 16 overviews the jet kinematics by compiling the data in the literature (see the caption for references). Multiwavelength VLBI and optical observations reveal both subluminal and superluminal features in proper motion, providing a global distribution of the jet velocity field V in M87. We display the value of Γβ in Figure 16 by using simple algebraic formulae with the bulk Lorentz factor ${\rm{\Gamma }}\equiv {(1-{\beta }^{2})}^{-1/2}$ and $\beta ={\beta }_{\mathrm{app}}/({\beta }_{\mathrm{app}}\cos {\theta }_{{\rm{v}}}+\sin {\theta }_{{\rm{v}}})$, where β = V/c and β_app is the apparent speed of the moving component in units of c. The value of Γβ approaches β in the nonrelativistic regime (${\rm{\Gamma }}\to 1$) and represents Γ in the relativistic regime ($\beta \to 1$), thereby representing simultaneously the full dynamic range in velocity over both regimes.

Superluminal motions (β_app > 1) have been frequently observed at relatively large distances beyond r_B. Furthermore, these components seem to originate at the location of HST-1 (Biretta et al. 1999; Cheung et al. 2007; Giroletti et al. 2012). On the other hand, no prominent superluminal features inside r_B have been confirmed in VLBI observations over recent decades (Reid et al. 1989; Kellermann et al. 2004; Ly et al. 2007). Instead, subluminal features are considered as nonbulk motions, such as growing instability patterns and/or standing shocks (e.g., Kovalev et al. 2007). Thus, this discrepancy (a gap between subluminal and superluminal motions along the jet axial distance) has been commonly recognized. Asada et al. (2014) discovered a series of superluminal components upstream of HST-1 (z/r_g ∼ 10⁵–10⁶), providing the missing link in the jet kinematics of M87.

Very recently, superluminal motions on the scale of z/r_g ≃ 10³–10⁴ were finally discovered by Mertens et al. (2016) and Hada et al. (2017). These observations give a diversity to the velocity picture and suggest the hypothesis that the systematic bulk acceleration is taking place if the observed proper motions indeed represent the underlying bulk flow. A smooth acceleration from subliminal to superluminal motions upstream of HST-1 is argued in the context of the MHD jet with an expanding parabolic nozzle (Nakamura & Asada 2013; Asada et al. 2014; Mertens et al. 2016; Hada et al. 2017), while observed proper motions exhibit a systematic deceleration in the region downstream of HST-1 (Biretta et al. 1995, 1999; Meyer et al. 2013), where the jet forms a conical stream.

Paired sub-/superluminal motions in optical/radio observations at HST-1 (Biretta et al. 1999; Cheung et al. 2007; see Figure 16 at ∼10⁶r_g) are modeled by the quad relativistic MHD shock system with a coherent helical magnetic field (Nakamura et al. 2010; Nakamura & Meier 2014). Taking the complex 3D kinematic features of trailing knots downstream of HST-1 (Meyer et al. 2013) into account, a growing current-driven helical kink instability associated with forward/reverse MHD shocks in the highly magnetized relativistic jet (Nakamura & Meier 2004) may be responsible for organizing the conical jet in M87 at the kiloparsec scale.

We examine here the jet kinematics with observations far upstream of HST-1 at z/r_g ≃ 10³–10⁴ (Kellermann et al. 2004; Kovalev et al. 2007; Hada et al. 2016, 2017; Mertens et al. 2016). The distribution of Γβ reaches a maximum level of ≃3 and extends to a lower value by more than two orders of magnitude as is shown in Figure 16. Mertens et al. (2016) interpret that the flow consists of a slow, mildly relativistic (Γβ ∼ 0.6; subluminal) layer (the exterior of the jet sheath), associated with either instability patterns or winds, and a fast, relativistic (Γβ ∼ 2.3; superluminal) layer (the jet sheath), which is accelerating a cold MHD jet from the Keplerian disk (i.e., the BP mechanism). Note that β_app ≃ 1 corresponds to Γβ ≃ 1.46 with θ_v = 14° in M87.

In our numerical results, maximum values of Γβ ≃ 0.8–2.6 (the solid cyan vertical line in Figure 16) are obtained at r_out/r_g = 100 (θ ≲ 10°), depending on the BH spin. This range covers most of the higher part of observed proper motions. For sufficiently high spins (a ≥ 0.9), bulk speeds of Γβ ≃ 1.7–2.6 could be associated with knotty structures (see Figure 12). Thus, we give an additional interpretation that superluminal motions could be interpreted as moving blobs in the underlying flow of the jet sheath. Regarding highly subluminal motions, as is shown in Figures 11 and 12, we confirm that a nonrelativistic coronal wind universally exists for a = 0.5–0.99 with Γβ ≳ 0.1 (see also Yuan et al. 2015, for a = 0), which may be responsible for slow motions immediately exterior the jet sheath. An acceleration of winds is also of interest to us, but it is unclear in our numerical results (see Yuan et al. 2015, for their behaviors at r/r_g ≳ 100).

Under the assumption that an observed moving component (β_app) represents an underlying bulk flow (e.g., Lister et al. 2009), we compare observations with steady axisymmetric FFE jet solutions in Figure 16. Γβ with Equation (6) is displayed with different BH spins (a = 0.5–0.99). Our numerical simulations reveal that b²/ρ ≃ 0.5–1 is sustained at the funnel edge along the outermost BP82-type parabolic streamlines of $z\propto {R}^{1.6}$. Therefore, significant acceleration through the FFE mechanism or magnetic field conversion is not expected here. Instead, the inner part of the funnel (the jet sheath/limb), where a high ratio of magnetic to rest-mass energy density would be expected, is an appropriate region to apply the FFE jet solution. Parabolic streamlines of z ∝ R^1.8 (a = 0.5–0.99) with θ_fp = π/3 are chosen as our reference solutions, taking into account that a peak in Γ lies asymptotically between two outermost streamlines ($z\propto {R}^{2}$ and $z\propto {R}^{1.6}$ in Figure 12).

A linear acceleration of highly magnetized MHD/FFE outflows can be expected in the moderately relativistic regime (Γ ≳ 2) with ${\rm{\Gamma }}\beta \propto R\propto {z}^{0.56}\ (\beta \to 1)$ as is shown in Figure 16. Similar results are obtained by Beskin & Nokhrina (2006), McKinney & Narayan (2007b), and Pu et al. (2015). Maximum values of Γβ (≃0.8–2.6) in our GRMHD simulations (a = 0.5–0.99) are qualitatively consistent with those of the FFE jet at z/r_g = 100. We, however, consider that this would be by coincidence, as we cannot find any smooth increase in Γ well beyond ∼2 within r/r_g = 100 in Figure 12. FFE jet solutions with the parabolic shape of $z\propto {R}^{1.8}$ indicate Γβ ∼ 4–30 around the scale of z/r_g ∼ 10³–10⁴ (a ≥ 0.5). This velocity range corresponds to β_app ≈ 4–8 with θ_v = 14° in M87. There is a clear discrepancy between observed proper motions and theoretical expectations.

In reality, AGN jets at VLBI scale may not be exactly described by the FFE system. An agreement between the GRMHD results and FFE models is found to be good as far as r/r_g ∼ 10³; beyond this scale the matter inertia becomes non-negligible (Γ ≳ 10) in GRMHD simulations (McKinney & Narayan 2007a). As a consequence, a slower evolution than ${\rm{\Gamma }}\beta \propto R$ may presumably take place. Nonetheless, a departure from Γ ∼ 2 at z/r_g ≳ 100 could be expected in a GRMHD simulation if b²/ρ is sufficiently large (≫1) at the jet stagnation surface. To be fair enough, Γβ  ≲ 7 is achieved at z/r_g ≃ 700 in McKinney (2006),²⁷ which is quantitatively consistent with the FFE jet with a = 0.9–0.99 (see Figure 16). On the other hand, for a moderate case of b²/ρ ≲ 100, maximum values of Γβ ≃ 3–4 are reported at z/r_g ≃ 1000 with a = 0.7–0.98 (Penna et al. 2013), indicating that a slower evolution is taking place than in highly magnetized GRMHD/FFE outflows (see also Figure 16).

As is mentioned above, the detection of faster proper motions β_app ≳ 4 (∼15 mas yr⁻¹) and a signature of their accelerations at z/r_g ∼ 10³–10⁴, where the jet sheath maintains a parabolic shape, will be key to confirming our GRMHD parabolic jet hypothesis. A VLBI program with 15/22/43 GHz toward M87 with a high-cadence monitoring of less than a week (conducting each observation every few days over a few weeks) may be feasible to find motions faster than ≳0.3 mas week⁻¹.

Tseng et al. (2016, hereafter T16) analyzed multifrequency data (VLBA, EVN, and VLA) to investigate the jet structure in NGC 6251 and detect a structural transition of the jet radius from a parabolic to a conical shape at (2–4) × 10⁵r_g, which is close to r_SOI ≃ 10⁶r_g in this source. This is a remarkably similar result to M87 (AN12); one may consider the virial equilibrium at the center of the cooling core in the giant elliptical galaxies as a thermodynamically stable state, which gives r_B ≈ r_SOI. Furthermore, the jet radii (in units of 2r_g) before/after the transition are quantitatively overlapped with M87 as is shown in Figure 3 of T16. Obviously, this implies a tight correlation between the jet sheath and the outermost BP82-type parabolic streamline of the FFE jet solution as seen in M87 (Figure 15). Figure 17 confirms this at a quantitative level in NGC 6251.

Zoom In Zoom Out Reset image size

T16 performed the broken power-law fitting and obtained z  ∝ R^2.0 at ≲4.2 × 10⁵r_g and z ∝ R^0.94 far beyond. We hereby suggest that the inner jet could be the BP82-type parabolic geometry, which is similar to M87 (Figure 17), if a position offset of VLBI cores from the SMBH ≃ 8 × 10³r_g in T16 is taken into account. Note that the radius of the EVN core at 1.6 GHz is almost identical to the radii of the VLBA jet at 15 GHz, as is shown in Figure 17, so we also confirm that the VLBI core can be identified as the innermost jet emission given at these frequencies, which is also similar to M87 (Hada et al. 2013; Nakamura & Asada 2013). By comparing Figures 15 and 17, we realize that the data points of M87 (inside r_B) are distributed across more orders of magnitude than NGC 6251 (inside r_SOI). Thus, VLBI observations at higher frequencies are needed to confirm a precise power-law index in the parabolic stream inside the SOI.

The difference of the SMBH mass between M87 and NGC 6251 is about one order of magnitude. If the jet radial and axial sizes are normalized in units of r_g, they are remarkably identical (see Figures 15 and 17 in this paper, and Figure 3 in T16). This suggests that the structural transition may be a characteristic of AGN jets, at least in nearby radio galaxies. It is straightforward to seek a counterpart by studying nearby blazars with a relatively large SMBH mass of M ∼ 10⁹–10¹⁰ M_⊙. Historically, the upstream region of conical jets in blazars (i.e., inside the VLBI cores at ≲0.1 mas²⁸ at millimeter/centimeter wavelengths) has been unresolved. They are sometimes called the "pipeline" from the central energy generator to the jet, which is unknown, and it is even said that it may not exist (see the schematic view; Marscher & Gear 1985). Therefore, we expect that ultrahigh angular resolution VLBI at millimeter (HSA, GMVA, EHT) wavelengths with ALMA will explore the nonconical pipeline inside r_SOI for bright nearby blazars.

A limb-brightened feature, one of the unanswered issues in the M87 jet, is discussed by Hada et al. (2016); according to the jet spine-sheath scenario, by the presence of a velocity gradient transverse to the jet (e.g., Ghisellini et al. 2005; Clausen-Brown et al. 2013),²⁹ viewing the M87 jet from an angle θ_v = 14° would not cause a limb-brightened feature unless the jet spine was unrealistically faster than the jet sheath, indicating that alternative processes are involved (see Hada et al. 2016, and references therein). Let us revisit this issue based on Figure 16; a relativistic beaming effect in the M87 jet is diagnosed with the Doppler factor δ = {Γ(1 − β cos θ_v)}⁻¹ on the scale of z/r_g ∼ 10³–10⁴, where β and Γ are adopted from observed proper motions and the MHD jet theory.

The observed synchrotron flux density S_ν for a relativistically moving component is enhanced by a beaming factor δ^3−α (Jorstad et al. 2005; Savolainen et al. 2010), where α is the spectral index defined as ${S}_{\nu }\propto {\nu }^{+\alpha }$. We here adopt α = −0.5 for reference, and the corresponding beaming factor is displayed as a function of the 4-velocity Γβ in Figure 18. The beaming factor becomes less than unity in the highly relativistic regime Γβ ≳ 30. By taking into account the observed proper motions in HSA 86 GHz, VLBA 15/43 GHz, and KaVA 22 GHz on the scale of z/r_g ∼ 10³–10⁴, which corresponds to Γβ ≃ 10⁻²–3 in Figure 16, the beaming factor is expected to be ∼1–100 at the jet sheath. Note that a similar range of beaming factors can be expected at Γβ ≃ 4–30, which is expected in the FFE jet solutions for the parabolic geometry (z ∝ R^1.8) with a = 0.5–0.99 at z/r_g ∼ 10³–10⁴ in Figure 16.

Zoom In Zoom Out Reset image size

As Equation (6) suggests, the linear acceleration of a highly magnetized MHD/FFE outflow decreases toward the polar region of the funnel when the power-law index of the parabola () becomes large. A quasi-homogeneous distribution of the magnetization σ(≈b²/ρ) along the colatitude angle (θ) at r/r_g ≲ 20 and the structure of the jet stagnation surface, as is revealed in our GRMHD simulations (see Figure 14), may provide a feasible reasoning for an efficient bulk acceleration at the outer jet sheath, where the magnetic nozzle effect would be expected under the progress of a differential bunching of the poloidal magnetic flux toward the polar axis.

Therefore, the jet spine can be intrinsically less beamed than the jet sheath as the distribution of Γ exhibits (Figure 12; see also Tchekhovskoy et al. 2008, 2010; Penna et al. 2013). A lateral stratification of Γ is naturally expected, so that a limb-brightened feature may be fundamental in AGN jets if they consist of BH-driven GRMHD outflows. Interestingly, limb-brightened features are also observed in the best-known TeV BL Lac objects Mrk 421 (Piner et al. 2010) and Mrk 501 (Giroletti et al. 2008; Piner et al. 2009) on z/r_g ≃ 10⁴–10⁵ even with smaller angles of θ_v = 4° (Giroletti et al. 2004; Lico et al. 2012).

Based on our numerical results up to r_out/r_g = 100, we find no relevant evidence of a concentration of the poloidal magnetic flux at the funnel edge (Gracia et al. 2009) and/or a pileup of the material along the funnel edge (Zakamska et al. 2008), which might be related to a funnel-wall jet (De Villiers et al. 2003, 2005), as a possible mechanism of a limb-brightened feature. Note that these may conflict with the physical conditions necessary to accelerate MHD jets, e.g., a high ratio of magnetic to rest-mass energy density and the magnetic nozzle effect. It would, however, be necessary to conduct a further investigation of this issue at the corresponding scale (z/r_g ≳ 10³).

We comment on the power-law acceleration (see footnote 24) in the jet sheath (possibly a less-collimated parabolic stream than the genuine parabolic one). As is shown in Figure 16, steady axisymmetric FFE jet solutions for streamlines of z ∝ R^1.8 (a = 0.5–0.99) with θ_fp = π/3 (as the jet sheath) do not exhibit a transition from the linear to power-law acceleration at $z/{r}_{{\rm{g}}}\lt {10}^{5}$ (it takes place at z/r_g > 10¹⁰ for a = 0.99). Thus, the outer jet sheath is always faster than the inner jet spine, even if the jet spine is launched with a sufficiently high value of b²/ρ at the jet stagnation surface and Γ follows the linear acceleration due to an efficient magnetic nozzle effect. Both of these factors, however, are not supported by our GRMHD simulations. Therefore, we suggest that the limb-brightened feature in M87 may be associated with the intrinsic property of an MHD parabolic jet powered by the spinning BH, rather than the result of a special viewing angle as is previously discussed in Hada et al. (2016).

Finally, as is mentioned in Section 4.1, the radius of the jet sheath starts to deviate slightly (becoming narrower) from the outermost BP82-type streamline z ∝ R^1.6 at $z/{r}_{{\rm{g}}}\gtrsim {10}^{4}$ (see also Figure 15). If the jet sheath follows the linear acceleration up to this scale, as is examined in Figure 16, the underlying flow would reach Γβ ≃ 30 (a = 0.9–0.99) and result in a weaker Doppler deboosting (Figure 18). Furthermore, it is interesting to note that the emission of the parabolic jet sheath farther downstream disappears at z/r_g ≳ 4 × 10⁵ (Asada et al. 2014), where θ_j ≃ 05 is obtained (Figure 15). If the empirical relation Γθ_j ∼ 0.1–0.2 (Clausen-Brown et al. 2013) is applied, Γ ≃ 11–22 would be expected. This is close to the velocity range at which a Doppler deboosting may arise.

We now discuss the (sub)millimeter VLBI cores in M87, which are considered the innermost jet emission—at the given frequencies—in the vicinity of the SMBH (see also Figure 15). Figure 19 shows the radius and location of VLBI cores at millimeter bands (43, 86, and 230 GHz) and their expectation at submillimeter bands (345 and 690 GHz), by an extrapolation of the VLBI core at frequencies higher than 43 GHz (Hada et al. 2013; Nakamura & Asada 2013) and utilizing the frequency-dependent VLBI core shift (Blandford & Königl 1979) in M87 (Hada et al. 2011). Our GRMHD simulation result for a = 0.9 is overlaid for reference. What we currently know about the (sub)millimeter VLBI cores of M87 from observations are the size, the flux density, the brightness temperature (Doeleman et al. 2012; Akiyama et al. 2015), and the energetics (Kino et al. 2014, 2015).

Zoom In Zoom Out Reset image size

5.3.1. (Sub)millimeter VLBI Core as a Neighborhood of the Jet Origin

5.3.2. Origin of Superluminal Blobs and Shock-in-jet Hypothesis

5.3.3. Comparisons with Other Models and Future Studies