Towards a Full MHD Jet Model of Spinning Black Holes--I: Framework and a split monopole example

TL;DR

This paper develops a comprehensive framework for modeling magnetohydrodynamical jets from spinning black holes, solving coupled equations numerically for specific magnetic configurations, and explores energy extraction mechanisms like the magnetic Penrose process.

Contribution

It introduces a self-consistent numerical approach to solve MHD jet equations for spinning black holes, including a split monopole example and analysis of energy extraction processes.

Findings

Magnetic field lines are dragged by fluid inertia.

Fluid alone does not extract energy from the black hole.

Magnetic Penrose process can operate for charged components at low plasma loading.

Abstract

In this paper, we construct a framework for investigating magnetohydrodynamical jet structure of spinning black holes (BHs), where electromagnetic fields and fluid motion are governed by the Grad-Shafranov equation and the Bernoulli equation, respectively. Assuming steady and axisymmetric jet structure, we can self-consistently obtain electromagnetic fields, fluid energy density and velocity within the jet, given proper plasma loading and boundary conditions. Specifically, we structure the two coupled governing equations as two eigenvalue problems, and develop full numerical techniques for solving them. As an example, we explicitly solve the governing equations for the split monopole magnetic field configuration and simplified plasma loading on the stagnation surface where the poloidal fluid velocity vanishes. As expected, we find the rotation of magnetic field lines is dragged down by fluid inertia, and the fluid as a whole does not contribute to energy extraction from the central BH, i.e., the magnetic Penrose process is not working. However, if we decompose the charged fluid as two oppositely charged components, we find the magnetic Penrose process does work for one of the two components when the plasma loading is low enough.

Equations

$$\mathrm{d}{s}^2$$

$$F_{r\varphi }$$

$$F_{tr}$$

$$F_{r\theta }$$

$$\frac{u^r}{{\Psi }_{,\theta }}=-\frac{u^{\theta }}{{\Psi }_{,r}}=\frac{(u^{\varphi }-\Omega u^t)}{F_{r\theta }},$$

$$\eta$$

$$=\frac{\sqrt{-g}n(u^{\varphi }-\Omega u^t)}{F_{r\theta }}.$$

$$T_{\mathrm{EM}}^{\mu \nu }$$

$$T_{\mathrm{MT}}^{\mu \nu }$$

$$E$$

$$L$$

$$D_{\Psi }^{\parallel }\equiv \frac{1}{\sqrt{-g}}({\Psi }_{,\theta }{\partial }_r-{\Psi }_{,r}{\partial }_{\theta }),$$

$$D_{\Psi }^{\parallel }\eta =(n{u}^{\mu }{)}_{;\mu },$$

$$T_{;\nu }^{\mu \nu }=S^{\mu },$$

$$D_{\Psi }^{\parallel }(\eta E)$$

$$D_{\Psi }^{\parallel }(\eta L)$$

$${\left(n{u}^{\mu }\right)}_{;\mu }$$

$$\mathcal{F}(u)=u_p^2+1-{\left(\frac{E}{m}\right)}^2{U}_g(r,\theta )=0,$$

$$U_g(r,\theta )$$

$$k_0$$

$$k_2$$

$$k_4$$

$${\mathcal{M}}^2=\frac{4\pi m{\eta }^2}{n}=4\pi mn\frac{u_p^2}{B_p^2}=4\pi m\eta \frac{u_p}{B_p},$$

$$(\sqrt{-g}{B}_p{)}^2=g_{rr}({\Psi }_{,\theta }{)}^2+g_{\theta \theta }({\Psi }_{,r}{)}^2.$$

$$k_{0}\big{|}_{r=r_{\rm LS}}=0\ .$$

$$-k_{0}+\mathcal{M}^{2}\big{|}_{r=r_{A}}=0\ .$$

$$\frac{E}{L}=-\frac{g_{tt}+g_{t\phi}\Omega}{g_{t\phi}+g_{\phi\phi}\Omega}\Bigg{|}_{r=r_{A}}\ ,$$

$$D^{\parallel}_{\Psi}k_{0}\big{|}_{r=r_{*}}=0\ .$$

$$\left(\frac{E}{m}\right)^{2}=\frac{k_{0}}{k_{2}}\Bigg{|}_{r=r_{*}}\ ,$$

$${\sigma }_{\ast }^{\mathrm{in};\mathrm{out}}=\frac{B_{p,\ast }}{4\pi m|\eta {|}_{\mathrm{in};\mathrm{out}}},$$

$$\sum _{i=0}^4{A}_i{u}_p^i$$

$$A_4$$

$$A_3$$

$$A_2$$

$$A_1$$

$$A_0$$

Full text

magnetic fields –magnetohydrodynamics (MHD) –black hole physics –galaxies: jets

††journal: ApJ

1 Introduction

Relativistic jets launched by accreting black holes (BHs) play an essential role in several energetic astrophysical phenomena, including stellar-mass BH X-ray binaries, active galactic nuclei and possibly gamma-ray bursts. After decades of debates among astrophysical communities, many open questions concerning the nature of the BH jets still remain to be answered (see e.g., Meier et al., 2001; Blandford et al., 2018, for reviews). To name a few of the most fundamental ones: what are the central engines of the jets; how the fluid within the jets is accelerated to the relativistic speed; how the jets are collimated.

The Blandford-Znajek (BZ) mechanism (Blandford & Znajek, 1977; Znajek, 1977), which describes an electromagnetic (EM) process of extracting the rotation energy of the central BH in the form of Poynting flux, is believed to be the most promising candidate for the central engines of the BH jets. For understanding the jets powered by the BZ mechanism, one needs to study magnetohydrodynamic (MHD) process in the Kerr spacetime, where the EM fields and the fluid motion are coupled in a complicate way. Therefore people treat some components of the MHD process as dynamical variables with other components prescribed in many previous studies on this subject. For studying the EM fields of the jet, force-free electromagnetodynamics (FFE) is a convenient assumption, where the fluid energy is ignored and the EM fields are self-contained (e.g., Tanabe & Nagataki, 2008; Komissarov, 2001, 2002, 2004a, 2004b, 2005; Komissarov & McKinney, 2007; Beskin & Zheltoukhov, 2013; Contopoulos et al., 2013; Nathanail & Contopoulos, 2014; Gralla & Jacobson, 2014; Gralla et al., 2015, 2016; Pan & Yu, 2014, 2015a, 2015b, 2016; Pan et al., 2017; Pan, 2018; Yang & Zhang, 2014; Yang et al., 2015; East & Yang, 2018; Mahlmann et al., 2018). For studying the fluid motion within the jet, people usually treat the fluid as test fluid in prescribed EM fields (e.g., Takahashi et al., 1990; Beskin & Nokhrina, 2006; Globus & Levinson, 2013, 2014; Pu et al., 2012, 2015). There have also been some full MHD attempts where only outflow pattern and jet structure in weak gravity regime are addressed (see e.g., Polko et al., 2013, 2014; Ceccobello et al., 2018, for self-similar outflow solutions in pseudo-potential) and (see e.g., Beskin et al., 1998; Beskin & Malyshkin, 2000; Lyubarsky, 2009; Tchekhovskoy et al., 2009; Beskin & Nokhrina, 2010; Beskin et al., 2017, for outflow solutions in Minkowski spacetime). For understanding the physics of BH accretion systems, general relativistic MHD (GRMHD) simulation has been another powerful tool in the past two decades, in which full MHD equations in curved spacetime are solved. Nevertheless, GRMHD codes tend to become unstable to a vacuum and therefore a matter density floor is usually introduced (e.g., Gammie et al., 2003; Shibata & Sekiguchi, 2005; Porth et al., 2017), which may obscure our understanding of plasma loading and the flow within the jet.

Besides all the theoretical explorations summarized above, substantial progress in spatial resolution has been made on the observation side. Especially the Event Horizon Telescope (EHT) (e.g., Doeleman et al., 2008, 2012; Ricarte & Dexter, 2015; EHT Collaboration, 2019a, b) is expected to resolve the structure of supermassive BHs nearby (Sgr A∗ and M87) down to horizon scales. It is promising to unveil the physical nature of the jets in these systems if the coming EHT observations can be correctly deciphered. This motivates us to construct a full GRMHD jet model, considering that all the previous studies are confined by different limitations.

In this paper, we aim to construct a GRMHD framework for investigating the structure of steady and axisymmetric jets of spinning BHs, in which the EM fields and fluid motion are self-consistently taken care of. A cartoon picture in Fig. 1 is shown to illustrate the major elements of our jet model: a central BH, EM fields, a plasma loading zone, inflow and outflow. The magnetic field lines penetrate the event horizon of a spinning BH and extract the rotation energy from the BH in the form of Poynting flux. Quantifying plasma loading within the BH jet is also a complicate problem considering the rich plasma sources, including the accretion flow centered on the equatorial plane, pair production inside the jet (Levinson & Rieger, 2011; Broderick & Tchekhovskoy, 2015; Hirotani & Pu, 2016; Chen et al., 2018) and neutrino pair annihilation from an extremely hot accretion flow (see e.g. Popham et al., 1999; Narayan et al., 2001). In our jet model, we do not deal with these detailed processes. For convenience, we introduce a plasma loading zone where plasma is injected and prescribe the loading function, i.e., particle number flux per magnetic flux $\eta(r,\theta)$. Under the mutual influence of the central BH and the EM fields, the inflow and outflow pattern naturally forms. In summary: we aim to construct a framework for investigating MHD jet structure of spinning BHs, in which the EM fields and the fluid motion are self-consistently obtained given proper boundary conditions and a proper plasma loading function $\eta(r,\theta)$.

This paper is organized as follows. In Section 2, we summarize some basic equations and assumptions to be used in this paper. We derive the two governing equations: the Bernoulli equation and the MHD Grad-Shafranov (GS) equation in Section 3 and Section 4, respectively. We detail the numerical techniques for solving the governing equations in Section 5. The numerical solutions of MHD jet structure with split monopole magnetic field configuration are presented in Section 6. Summary and discussion are given in Section 7. For reference, we place some details for deriving the governing equations in Appendix A and B. Throughout this paper, we use the geometrical units $c=G=M=1$, where $M$ is the mass of the central BH.

2 Basic Setting Up

The background Kerr metric is written in the Boyer-Lindquist coordinates as follows,

[TABLE]

where $a$ and $M$ are the BH spin and mass, respectively, $\Sigma=r^{2}+a^{2}\cos^{2}\theta$, $\Delta=r^{2}-2Mr+a^{2}$, $\beta=(r^{2}+a^{2})^{2}-a^{2}\Delta\sin^{2}\theta$ and the square root of the determinant $\sqrt{-g}=\Sigma\sin\theta$.

We investigate the structure of a steady and axisymmetric BH jet and we assume the plasma within the jet is perfectly conducting, i.e., $\partial_{t}=\partial_{\phi}=0$ and $\mathbf{E}\cdot\mathbf{B}=0$, where $\mathbf{E}$ and $\mathbf{B}$ are the electric and the magnetic fields, respectively. Then all the non-vanishing components of Maxwell tensor are expressed as follows (see e.g., Pan & Yu, 2014)

[TABLE]

where $\Psi=\Psi(r,\theta)$ is the magnetic flux and $\Omega=\Omega(\Psi)$ is the angular velocity of magnetic field lines. For convenience, we have defined poloidal electric current $I(r,\theta)\equiv\sqrt{-g}F^{\theta r}$. Therefore, the EM fields are completely determined by three quantities: $\{\Psi(r,\theta),\Omega(\Psi),I(r,\theta)\}$.

Before proceeding on, it is useful to define a few conserved quantities. From the perfectly conducting condition $F_{\mu\nu}u^{\nu}=0$, we find different components of fluid velocity are related by

[TABLE]

from which we can define the particle number flux per magnetic flux

[TABLE]

From the energy-momentum tensor $T^{\mu\nu}=T^{\mu\nu}_{\rm EM}+T^{\mu\nu}_{\rm MT}$, where the EM part and the matter (MT) part are

[TABLE]

we define total energy per particle $E$ and total angular moment $L$ per particle as follows,

[TABLE]

where $\rho$, $n$ and $m$ are the proper energy density, the proper number density and the particle rest mass, respectively; and we have assumed cold plasma.

Now let us examine the conservation property of these quantities along magnetic field lines. For this purpose, we define derivative along field lines

[TABLE]

and it is straightforward to obtain (see Appendix A)

[TABLE]

i.e., $D^{\parallel}_{\Psi}\eta$ quantifies the plasma loading rate. In general, we can write the energy-momentum conservation as

[TABLE]

where the source term $S^{\mu}$ comes from plasma loading. As a simple example, we assume $S^{\mu}=(D^{\parallel}_{\Psi}\eta)mu^{\mu}$ in this paper, i.e., the source term is contributed by the kinetic energy of newly loaded plasma. With a few steps of calculation as detailed in Appendix A, we obtain

[TABLE]

where $\eta E$ and $\eta L$ are the energy flux per magnetic flux and angular momentum flux per magnetic flux, respectively. Outside the plasma loading zone where there is no particle injection, the particle number conservation reads as

[TABLE]

therefore $\eta,E,L$ are conserved along field lines, i.e., $\eta=\eta(\Psi),E=E(\Psi),L=L(\Psi)$.

In summary: with assumptions of steady and axisymmetric jet structure and perfectly conducting plasma within the jet, we have obtained one conserved quantity $\Omega(\Psi)$ and three “quasi-conserved” quantities $\{\eta(\Psi),E(\Psi),L(\Psi)\}$ which are only conserved outside the plasma loading zone.

3 Bernoulli equation

From the normalization condition $u^{\mu}u_{\mu}=-1$ and Eqs. (4,6), we obtain the relativistic Bernoulli equation

[TABLE]

where $u_{p}^{2}\equiv u^{A}u_{A}$ with the dummy index $A=\{r,\theta\}$. In the Kerr space-time, the characteristic function $U_{g}$ is writen as (Camenzind, 1986a, b, 1987; Takahashi et al., 1990; Fendt & Greiner, 2001; Fendt & Ouyed, 2004; Levinson, 2006; Pu et al., 2015)

[TABLE]

where

[TABLE]

and the Alfvén  Mach number $\mathcal{M}$ is given by

[TABLE]

with the poloidal magnetic field $B_{p}$ defined by

[TABLE]

Several characteristic surfaces can be defined according to the critical points of the flow velocity (see e.g., Michel, 1969, 1982; Camenzind, 1986a, b; Takahashi et al., 1990; Beskin, 2009, for details). The light surface (LS) is defined by where the rotation velocity of field lines approaches light speed and where particles are forbidden to corotate with the field lines,

[TABLE]

The Alfvén  surface is defined by where the denominator of characteristic function $U_{g}(r,\theta)$ vanishes, i.e.,

[TABLE]

On the Alfvén  surface, we find

[TABLE]

where we have used Eqs.(12-3). The stagnation surface where $u_{p}=0$ is determined by

[TABLE]

The fast magnetosonic (FM) surface and slow magnetosonic (SM) surface are defined by where the denominator of $D_{\Psi}^{\parallel}u_{p}$ vanishes. In the cold plasma limit, the SM surface coincides with the stagnation surface. On the stagnation surface, where both $u_{p}$ and $\mathcal{M}$ vanish, we find

[TABLE]

where we have used Eqs.(12, 13).

Plugging Eq.(4) into Eq.(12), we find that the Bernoulli equation is a polynomial equation of fourth order in $u_{p}$ with to-be-determined eigenvalue $E/L$ (or equivalently the location of the Alfvén  surface $r_{A}$), given prescribed angular velocity $\Omega$ and particle number flux per field line $\eta(r,\theta)$ (see e.g., Camenzind, 1986a, b; Takahashi et al., 1990; Fendt & Greiner, 2001; Pu et al., 2012, 2015).

3.1 Single Loading Surface

As a first step towards a full MHD jet solution, we mathematically idealize the plasma loading zone as a single surface, and we choose the stagnation surface (Eq. (20)) as the plasma loading surface (see e.g., Broderick & Tchekhovskoy, 2015, for a detailed gap model) in this paper. To define the plasma loading for both inflow and outflow, we introduce a pair of dimensionless magnetization parameters on the loading surface,

[TABLE]

where $B_{p,*}$ is the poloidal field on the loading surface. In this way, the particle number flux per magnetic flux $\eta$ is completely determined by $\sigma_{*}^{\rm in;out}$, recalling that $\eta$ is a conserved quantity along field lines outside the loading zone. Note that $\eta_{\rm in}<0$ and $\eta_{\rm out}>0$, therefore there is a jump in $\eta$ at the loading surface, i.e., $D_{\Psi}^{\parallel}\eta\propto\delta(r-r_{*})$.

Using Eq.(21), the Bernoulli equation (12) can be rewritten into a fourth-order polynomial equation as

[TABLE]

where the coefficients $A_{i}$ are functions expressed by

[TABLE]

As explored in several previous studies (e.g., Takahashi et al., 1990; Pu et al., 2015), solving the Bernoulli equation above is in fact an eigenvalue problem, where $(E/L)_{\rm in}$ is the to-be-determined eigenvalue ensuring inflow smoothly cross the FM surface, while $(E/L)_{\rm out}$ is given by the match condition on the loading surface. Eq.(10) provides conditions connecting the inflow and the outflow for the single surface loading,

[TABLE]

which give the match condition

[TABLE]

where $\delta\eta\equiv\eta_{\rm out}-\eta_{\rm in}$, and we have used the fact that $D_{\Psi}^{\parallel}\eta$ is a $\delta$-function centered on the loading surface in deriving Eq. (3.1). It is straightforward to see Eq. (3.1) guarantees a same jump in the total energy flux and its matter component, therefore the Poynting flux (and all the EM field components) is continuous across the loading surface.

As along as the Bernoulli equation is solved, i.e., both the eigenvalues $(E/L)_{\rm in,out}$ and the poloidal velocity field $u_{p}$ are obtained, $u^{r}$ and $u^{\theta}$ is obtained via Eq.(4) and $u_{p}^{2}=u^{A}u_{A}$, while $u_{t}$ and $u_{\phi}$ are obtained via relation $m(u_{t}+\Omega u_{\phi})=-(E-\Omega L)$ and the normalization condition $u\cdot u=-1$.

Before delving into the details of numerically solving the Bernoulli equation, we can now give an estimate of the eigenvalues. Combining the definitions of $E$ and $L$ Eqs.(6) with Eq.(12), we find

[TABLE]

plugging which back into Eqs.(6), we obtain

[TABLE]

which imply $E/L=\Omega[1+O(\sigma_{*}^{-1})]$ and we have used the fact $\eta_{\rm in}<0$ and $\eta_{\rm out}>0$.

4 MHD Grad-Shafranov equation

With the aid of the Maxwell’s equation

[TABLE]

the trans-field component of the energy conservation equation (9) is written as 111Eq. (30) only holds for the specific choice of source function $S^{\mu}=(nu^{\nu})_{;\nu}mu^{\mu}$.

[TABLE]

where we have used Eq. (8) and the source function $S^{\mu}$. The repeated Latin letters $A$ and $C$ run over the poloidal coordinates $r$ and $\theta$ only (Nitta et al., 1991; Beskin & Pariev, 1993; Beskin, 1997). This is known as the MHD GS equation, with the $1^{\rm st}$ and $2^{\rm nd}$ terms in the bracket being the fluid acceleration and the electromagnetic force, respectively.

After some tedious derivation (see Appendix B), we write the full MHD GS equation in a compact form

[TABLE]

Here $\mathcal{L}$ is a differential operator defined by

[TABLE]

where $\mu=\cos\theta$, $\mathcal{A}(r,\theta;\Omega)=-k_{0}(r,\theta;\Omega)+\mathcal{M}^{2}(r,\theta)$, and we have defined $\partial^{\Omega}_{A}(A=r,\mu)$ as the partial derivative with respect to coordinate $A$ with $\Omega$ fixed, $\partial_{\Omega}$ as the derivative with respect to $\Omega$, $D_{\Psi}^{\perp}$ as the derivative perpendicular to field lines

[TABLE]

which is equivalent to the ordinary derivative $d/d\Psi$ when acting on functions of $\Psi$. The two source terms are

[TABLE]

where $I=4\pi(\eta L-\eta mu_{\phi})$ [see Eq.(6)].

In the FFE limit, $\mathcal{M}^{2}=0$, $\mathcal{S}_{\rm MT}=0$, and the GS equation reduces to (Pan et al., 2017)

[TABLE]

The FFE solutions $\{\Psi|_{\rm FFE},\Omega|_{\rm FFE},(\eta L)|_{\rm FFE}\}$ have been well explored both analytically and numerically in many previous studies (see e.g., Blandford & Znajek, 1977; Tanabe & Nagataki, 2008; Contopoulos et al., 2013; Pan & Yu, 2015a, b). Similar to the FFE case, solving the MHD GS equation (31) is also eigenvalue problem, where $\Omega$ and $4\pi\eta L$ are the to-be-determined eigenvalues ensuring field lines smoothly cross the two Alfven surfaces (Contopoulos et al., 2013; Nathanail & Contopoulos, 2014; Pan et al., 2017; Mahlmann et al., 2018).

5 A split monopole example

As previewed in the Introduction, we aim to construct a framework for investigating MHD jet structure of spinning BHs, in which the EM fields $(F_{\mu\nu})$ and the fluid motion $(n,u^{\mu})$ are self-consistently obtained given a proper plasma loading function $\eta(r,\theta)$ and proper boundary conditions.

In this section, we detail the procedure of consistently solving the two governing equations for an example of the split monopole magnetic field configuration around a rapidly spinning central BH with a dimensionless spin $a=0.95$. For simplicity, we explore two different scenarios with magnetization parameters $\sigma_{*}^{\rm out}=2\sigma_{*}^{\rm in}$ and $\sigma_{*}^{\rm out}=\sigma_{*}^{\rm in}$, respectively. Remember that the loading function $\eta(r,\theta)$ is completely determined by the magnetization parameters via the definition (22).

Boundary conditions used here are similar to those of force-free solutions. Explicitly, we choose $\Psi|_{\mu=0}=\Psi_{\rm max}$ on the equatorial plane, $\Psi|_{\mu=1}=0$ in the polar direction, $\Psi_{,r}|_{r=r_{\rm H}}=0$ and $\Psi_{,r}|_{r=\infty}=0$ for the inner and outer boundaries, respectively. Here $r_{\rm H}$ is the radius of the event horizon.

5.1 Numerical Techniques

We define a new radial coordinate $R=r/(r+1)$, confine our computation domain $R\times\mu$ in the region $[R(r_{\rm H}),R_{\rm max}]\times[0,1]$, and implement a uniform $256\times 64$ grid. In practice, we choose $R_{\rm max}=0.995$, i.e., $r_{\rm max}\approx 200M$.

The Bernoulli equation (12) and the MHD GS equation (31), governing the flow along the field lines and field line configuration, respectively, are coupled. So we solve them one by one in an iterative way:

[TABLE]

with $l=1,2,3,\cdots$ . In a given loop $l$, we solve the GS equation updating $\Psi$ and $\{\Omega,(\eta L)\}$ (with $\{n,u^{\mu}\}$ inherited from the previous loop $l-1$), ensuring field lines smoothly cross the two Alfvén  surfaces; in a similar way, we solve the Bernoulli equation updating $u^{\mu}$ and $(E/L)$ (with freshly updated $\Omega$ and $\Psi$ from solving the GS equation), ensuring a super-sonic inflow solution and an outflow solution satisfying the match condition (26). Combing solutions to both equations and definitions of $\{\eta,E,L\}$, we finally obtain all the desired quantities $\{F_{\mu\nu},n,u^{\mu}\}$ as functions of coordinates $r$ and $\theta$.

We activate the iteration with an initial guess

[TABLE]

where $\Omega_{\rm H}\equiv a/(r_{\rm H}^{2}+a^{2})$ is the BH angular velocity.

The numerical techniques for tackling the two eigenvalue problems are detailed as follows:

• Step 1

The MHD GS equation is a second-order differential equation which degrades to first order on the Alfvén  surfaces where $\mathcal{A}(r,\theta)=0$. Numerical techniques for dealing this problem have been well developed in previous force-free studies (Contopoulos et al., 2013; Nathanail & Contopoulos, 2014; Huang et al., 2016, 2018; Pan et al., 2017; Mahlmann et al., 2018), and we briefly recap them here.

In each loop $l$, we solve the GS equation (31) with the approximate solution obtained from the previous loop $\left\{\Omega^{(l-1)},(\eta L)^{(l-1)},\Psi^{(l-1)}\right\}$ as the initial guess. We evolve the flux function $\Psi^{(l)}$ using the overrelaxation technique with Chebyshev acceleration (Press, 1986), and $\Psi^{(l)}(r,\theta)$ is updated on grid points except those in the vicinity of the two Alfvén  surfaces. The flux function $\Psi^{(l)}(r,\theta)$ on the Alfvén  surfaces are obtained via interpolation from neighborhood grid points and the directional derivatives on the Alfvén  surfaces (Pan et al., 2017). Usually we obtain two different flux function $\Psi(r_{\rm A}^{-})$ versus $\Psi(r_{\rm A}^{+})$ on the Alfvén  surface via interpolations from grid points inside and outside, respectively. To decrease this discontinuity, we adjust $\Omega^{(l)}(\Psi)$ at the outer Alfvén (OA) surface:

[TABLE]

with $\Psi_{\rm new}=0.5[\Psi(r_{\rm OA}^{+})+\Psi(r_{\rm OA}^{-})]$, where the subscript old/new represents quantities before/after the above adjustment; and adjust both $\Omega^{(l)}(\Psi)$ and $(\eta L)^{(l)}(\Psi)$ at the inner Alfvén surface (IA):

[TABLE]

with $\Psi_{\rm new}=0.5[\Psi(r_{\rm IA}^{+})+\Psi(r_{\rm IA}^{-})]$.

After sufficient evolution, we obtain a converged solution $\{\Omega^{(l)},(\eta L)^{(l)},\Psi^{(l)}\}$ which ensures field lines smoothly cross the two Alfvén  surfaces.

• Step 2

The Bernoulli equation in the form of Eq.(23) is a fourth-order polynomial equation in $u_{p}$ (Camenzind, 1986a, b, 1987; Takahashi et al., 1990; Fendt & Greiner, 2001; Fendt & Ouyed, 2004; Levinson, 2006; Pu et al., 2015), where the FM point is a standard ‘X’-type singularity, while the Alfvén point turns out to be a higher-order singularity (Weber & Davis, 1967). Mathematically, a FM point is the location of a multiple root to the Bernoulli equation. The existence of FM point is very sensitive to the value of $(E/L)_{\rm in}^{(l)}$. For a slightly small value, there exists only sub-sonic solutions in the region $r<r_{A}^{\rm in}$, i.e., no multiple root. For a slightly larger value, on the other hand, there exists no global solution extending from $r_{\rm H}$ to $r_{A}^{\rm in}$. Only for some specific choice, there exist a global super-sonic solution that crossing the FM point (see Fig. 3 for numerical examples). We adjust $(E/L)_{\rm in}^{(l)}$ on each field line until an inflow solution that smoothly crosses the FM point is found.

With the inflow solution in hand, $(E/L)_{\rm out}^{(l)}$ on each field line is uniquely determined by the match condition Eq.(26), then it is straightforward to compute the outflow velocity.

• Step 3

Combining the inflow and the outflow solutions from Step 2, the global fluid velocity $u^{(l)}$ and therefore the number density $n^{(l)},$ and the Mach number $\mathcal{M}^{(l)}$ are obtained along each field line. We feed these quantities into the GS equation (31) for the next loop. We iterate Step 1 to Step 3 until all quantities converge to a given precision.

We should point out that there is an unphysical divergence arising from the idealized plasma loading on a single surface adopted in this paper. The particle number flux function $\eta$ is negative/positive for inflow/outflow and is not continous on the loading surface. According to the definition of $\eta$ in Eq.(4), the particle number density $n$ is proportional to $\eta/u_{p}$ and therefore diverges on the stagnation surface where $u_{p}=0$. The particle number density $n$ show up in the source terms of the MHD GS equation (31). To overcome the unphysical infinity, we smooth them in the vicinity of the stagnation surface before feeding into Eq.(31).222 In addition, we usually obtain two different angular momentum flux per magnetic flux $(\eta L)_{\rm Bern}$ from the Bernoulli equation and $(\eta L)_{\rm GS}$ from the GS equation, respectively. The former is obtained from $E/L$ and Eq.(21), while the latter is one of the eigenvalues of the GS equation (31). In general, the two do not match exactly, where $(\eta L)_{\rm GS}$ does approach $(\eta L)|_{\rm FFE}$ in the limit of $\sigma_{*}\rightarrow\infty$, while $(\eta L)_{\rm Bern}$ does not. For the cases investigated in Section 6, we find $(\eta L)_{\rm GS}$ is different from $(\eta L)_{\rm Bern}$ by $\lesssim 25\%$. This mismatch indicates that the particle number flux per magnetic flux $\eta(r,\theta)$ cannot be arbitrarily given, or perfectly conducting fluid is not a sufficient description here.

6 Numerical Results

6.1 Case Study

In this subsection, we explore two nearly force-free cases with $\sigma_{*}\gg 1$: Case 1 with magnetization parameters $\sigma_{*}^{\rm out}=\sigma_{*}^{\rm in}=75$ and Case 2 with $\sigma_{*}^{\rm out}=2\sigma_{*}^{\rm in}=60$.

In Fig. 2, we compare the magnetic field line configurations of MHD solutions with that of FFE version. For both cases, which are nearly force-free, the deviation from the FFE version is as expected small, though the deviation of Case 1 from the FFE solution is more obvious than that of Case 2, due to more efficient outflow acceleration in the former case (see Fig. 5).

In Fig. 3, we show the poloidal velocity $u_{p}$. To clarify the different singularity types of FM points and Alfvén  points, we explicitly show all the solutions (both physical one and unphysical ones) to the Bernoulli equation (23) for the field line with foot point $\mu(r_{\rm H})=0.508$, where the solutions with correct eigenvalue $E/L$ are shown in solid lines and the solutions with $E/L$ slightly off the correct value are shown in dotted and dashed lines. At the FM point which is a ‘X’-type singularity, there exists a multiple root of the Bernoulli equation for correct eigenvalue $E/L$, while no global solution (dotted lines) or no multiple root exists (long-dashed lines) for $E/L$ of slightly off value. At the Alfvén  point which is a higher-order singularity, there exists multiple root no matter $E/L$ takes the correct eigenvalue or not. For the physical solution, the inflow passes along $r_{*}\to r_{\rm A}\to r_{\rm FM}\to r_{\rm H}$ and the outflow passes along $r_{*}\to r_{\rm A}\to\infty$.

In Fig. 4, we show the fluid angular velocity $\Omega_{\rm MT}$, the angular velocity of magnetic field lines $\Omega$ (solid lines), and the eigvenvalues $(E/L)_{\rm in;out}$ (dashed grey lines). Consistent with our intuition, we find the rotation of magnetic field lines is dragged down by the fluid inertia compared with the FFE case, i.e., $\Omega|_{\rm MHD}<\Omega|_{\rm FFE}$. Due to nonzero poloidal velocity of both inflow and outflow, the fluid does not corotate with the field lines, i.e., $\Omega_{\rm MT}<\Omega$ for inflow and $\Omega_{\rm MT}>\Omega$ for outflow. Specifically, the fluid angular velocity on the event horizon $\Omega_{\rm MT}(r=r_{\rm H})$ slightly exceeds the BH angular velocity $\Omega_{\rm H}$, which guarantees the fluid energy to be positive on the horizon (see Fig. 5).

In Fig. 5, we show the specific particle energy $-u_{t}$ for the two cases. Both of them are positive everywhere, while the outflow of Case 1 gains more efficient acceleration.

6.2 Energy Extraction Rates

In this subsection, we investigate the energy extraction rate from the central BH via the MHD jet, which is defined as

[TABLE]

where we have used Eqs.(4-6) in the second line. In the third line, $E/L$ and $\eta L$ are the eigenvalues of the Bernoulli equation (12) and of the GS equation (31), respectively. In the similar way, we can define its matter/electromagnetic component as

[TABLE]

We measure these energy extraction rates at $r=r_{\rm H}$/ $r=\infty$, and quantify their dependence on the magnetization parameter $\sigma_{*}$. In practice, we find that these energy extraction rates are not sensitive to the value of $\sigma_{*}^{\rm out}$, except the matter component of energy rate at infinity $\dot{E}_{\rm MT}^{\infty}$. Without loss of generality, we only show the rates in relation to $\sigma_{*}^{\rm in}$ for the $\sigma_{*}^{\rm out}=2\sigma_{*}^{\rm in}$ scenario in Fig. 6, where all the rates are displayed in unit of the energy extraction rate in the force-free limit $\dot{E}_{\rm FFE}\approx 0.4(\Psi_{\rm max}^{2}/4\pi)$.

As we see in Fig. 4, the rotation of magnetic field lines is dragged down by the loaded plasma, i.e., $\Omega|_{\rm MHD}<\Omega|_{\rm FFE}$, while the fluid that does not corotate with the field lines with angular velocity $\Omega_{\rm MT}|_{\rm inflow}>\Omega>\Omega_{\rm MT}|_{\rm outflow}$, tends to bend the field lines and induce a stronger $\phi$-component of magnetic field, i.e., $I|_{\rm MHD}>I|_{\rm FFE}$. The net result is that the Poynting energy extraction rate on the event horizon $\dot{E}_{\rm Poynting}^{\rm H}$ has little dependence on the magnetization. Going outward along the field lines, part of the Poynting flux is converted into the fluid kinetic energy. For the case with magnetization parameter $\sigma_{*}^{\rm in}=30$, the matter component makes up about $13\%$ of the total energy flux at infinity.

6.3 Penrose Process

An implicit assumption in our MHD jet model is the two-fluid description, since the electric current density $j^{\mu}$ is not proportional to the fluid velocity $u^{\mu}$. Therefore we can decompose the charged fluid as two oppositely charged components, positron ($e^{+}$) and electron ($e^{-}$). 333Though there is a degree of freedom in doing this decomposition, e.g., we can also decompose the fluid as electrons and ions, it does not change our conclusion qualitatively. We denote the number densities and the velocity fields as $n_{\pm}$ and $u^{\mu}_{\pm}$, respectively, which are related to $j^{\mu}$ and $nu^{\mu}$ via relations

[TABLE]

Consequently, we obtain

[TABLE]

and we can decompose the matter energy flux into two components $\dot{E}_{e^{\pm}}$. Here we are only interested in the energy extraction rates on the event horizon

[TABLE]

As an example, we choose the horizon enclosed magnetic flux $\Psi_{\rm max}=1000(m/e)$, and show $\dot{E}^{\rm H}_{e^{+}}/\dot{E}_{\rm FFE}$ and $\dot{E}^{\rm H}_{e^{-}}/\dot{E}_{\rm FFE}$ in relation to $\sigma_{*}^{\rm in}$ in Fig. 7. The energy extraction rate from positrons is always negative, while the energy extraction rate from electrons, become positive when the plasma loading is low enough. In this regime, denoted in shades in Fig. 7, the magnetic Penrose process is working, though, only for one of the two charged component. 444We should not do any quantitative interpretation for the results of this subsection, because the two-fluid decomposition done here is not accurate, e.g., there is no guarantee for the velocity of each component $u^{\mu}_{\pm}$ to be timelike and normalized. We will leave a more accurate two-fluid description of MHD jet structure (Koide, 2009; Liu et al., 2018) to future work . This finding is in good agreement with recent particle-in-cell simulations (Parfrey et al., 2018).

7 Summary and Discussion

7.1 Summary

To describe the MHD structure of BH jets, we need a minimum set of quantities as functions of spacetime: Maxwell tensor $F_{\mu\nu}$, fluid rest mass density $\rho$ (or equivalently particle number density $n$), and fluid four-velocity $u^{\mu}$. For determining all these quantities self-consistently, we constructed a full MHD framework, in which EM fields and fluid motion are governed by the MHD GS equation (31) and the Bernoulli equation (12), respectively. From these two governing equations, we can completely determine $\{F_{\mu\nu},\rho,u^{\mu}\}$ given proper boundary conditions and a proper plasma loading function $\eta(r,\theta)$ (see Eq.(4)). As an example, we consider a split monopole field configuration and idealized plasma loading on the stagnation surface.

Assuming steady and axisymmetric jet structure, and perfectly conductive plasma within the jet, the EM fields are completely determined by three functions: the magnetic flux $\Psi(r,\theta)$, the angular velocity of magnetic field lines $\Omega(\Psi)$ and the poloidal electric current $I(r,\theta)$ (see Eq.(2)). Given fluid energy density $\rho$ and velocity $u^{\mu}$, the MHD GS equation (31) turns out to be a second-order differential equation with respect to $\Psi(r,\theta)$ which degrades to be first-order on the two Alfvén  surfaces. Solving the GS equation is an eigenvalue problem, with eigenvalues $\Omega(\Psi)$ and $I(r,\theta)$ (or more precisely, the conserved quantity $4\pi\eta L(\Psi)$ defined in Eq.(6)) to be determined ensuring field lines smoothly cross the Alfvén  surfaces.

Given EM fields $F_{\mu\nu}$, the Bernoulli equation turns out to be a fourth-order polynomial equations in the poloidal fluid velocity $u_{p}$. Solving the Bernoulli equation is also an eigenvalue problem, with the eigenvalue $(E/L)_{\rm in}$ to be determined ensuring the inflow smoothly cross the FM surface, and $(E/L)_{\rm out}$ to be determined by the match condition (26) on the loading surface. With both $E/L$ and $u_{p}$ obtained, it is straightforward to obtain $n$ and $u^{\mu}$ via Eqs.(4,6) and the normalization condition $u\cdot u=-1$.

The two governing equations are coupled, therefore we numerically solved them in an iterative way (see Sec. 5.1). As a result, we find the rotation of magnetic field lines is dragged down by the plasma loaded, i.e., $\Omega|_{\rm MHD}<\Omega|_{\rm FFE}$; for the fluid angular velocity, we find $\Omega_{\rm MT}|_{\rm outflow}<\Omega<\Omega_{\rm MT}|_{\rm inflow}$ ; the non-corotating fluid tends to bend the field lines and induce a stronger $\phi$-component of magnetic field, therefore a stronger poloidal electric current i.e., $I|_{\rm MHD}>I|_{\rm FFE}$. The net result is that the Poynting energy extraction on the horizon is insensitive to the magnetization, i.e., $\dot{E}_{\rm Poynting}^{\rm H}|_{\rm MHD}\approx\dot{E}_{\rm Poynting}^{\rm H}|_{\rm FFE}$ (see Fig. 6). Going outward along the field lines, part of the Poynting flux is converted to the fluid kinetic energy. For the case we explored with $\sigma_{*}^{\rm in}=30$, the matter component makes up $\sim 13\%$ of the total energy flux at infinity.

Finally, we examined the MHD Penrose process for the cases we numerically solved. We found the specific fluid energy $-mu_{t}$ is always positive on the event horizon, i.e., the MHD Penrose process is not working and therefore the BZ mechanism defines fully the jet energetics. However, if we decompose the charged fluid as two oppositely charged components ($e^{\pm})$, we found the magnetic Penrose process does work for one of the two components when the plasma loading is low enough (see Fig. 7).

7.2 Discussion

As a first step towards a full MHD jet model, we have investigated the MHD jet structure of split monopole geometry assuming an idealized plasma loading on the stagnation surface. This simplified plasma loading gives rise to a few unphysical problems in the vicinity of the loading surface, including divergence of particle number density $n(r_{*})$, which shows up in the source terms of the MHD GS equation (34). To avoid the singularity arising from the unphysical divergence, we smoothed the function $n(r,\theta)$ in the vicinity of the loading surface. Another consequence of the simplified plasma loading is that we must impose the continuity equation (26) to ensure the EM fields to be continuous across the loading surface. As a result, $(E/L)_{\rm out}$ is specified by $(E/L)_{\rm in}$, i.e., $r_{A,\rm out}$ is specified by $r_{A,\rm in}$. Therefore, we lose the freedom to adjust $(E/L)_{\rm out}$ until a supersonic outflow solution is found as we did for the inflow solution. Consequently, all the outflow solutions obtained in this paper are subsonic (see Fig. 3).

In future work, we aim to investigate a full MHD jet model with a more realistic extending loading zone where the plasma injection is described by a continuous function $\eta(r,\theta)$. Then all the unphysical discontinuity and divergence described above would be avoided. For the extending plasma loading, the smooth EM fields would be naturally preserved, and the continuity requirement would not be a constraint. As a result, we can adjust $(E/L)_{\rm out}$ for finding a supersonic outflow solution, which is more consistent with recent observations (Hada et al., 2016; Mertens et al., 2016). In addition to the plasma loading, the BH surroundings also play an important role in shaping the jet structure (e.g. Tchekhovskoy et al., 2010; Beskin et al., 2017). The role of more realistic BH environment, including accretion flows and hot plasma with non-zero pressure will also be considered in future work.

Acknowledgements

We thank the referee for his/her careful reading of this manuscript and giving insightful suggestions. L.H. thanks the support by the National Natural Science Foundation of China (grants 11590784 and 11773054), and Key Research Program of Frontier Sciences, CAS (grant No. QYZDJ-SSW-SLH057). Z.P. thanks Hung-Yi Pu for his invaluable help throughout this research. Z.P. is supported by Perimeter Institute for Theoretical Physics. Research at Perimeter Institute is supported by the Government of Canada through the Department of Innovation, Science and Economic Development Canada and by the Province of Ontario through the Ministry of Research, Innovation and Science. C.Y. has been supported by the National Natural Science Foundation of China (grants 11521303, 11733010 and 11873103). This work made extensive use of the NASA Astrophysics Data System and of the astro-ph preprint archive at arXiv.org.

Appendix A Derivation of $\{D_{\Psi}^{\parallel}\eta,D_{\Psi}^{\parallel}(\eta E),D_{\Psi}^{\parallel}(\eta L)\}$

Using Eq.(4), it is straightforward to see

[TABLE]

from which we conclude that $D_{\Psi}^{\parallel}\eta$ is the source function of particle number density; i.e., $D_{\Psi}^{\parallel}(\eta)$ vanishes outside the plasma loading zone and is positive inside. Due to the existence of plasma loading, the energy conservation of the electromagnetic fields and plasma system is written as $T^{\mu\nu}_{\ \ ;\nu}=S^{\mu}$ [Eq. (9)], where the source term $S^{\mu}$ comes from plasma loading and in this paper we have assumed $S^{\mu}=(D_{\Psi}^{\parallel}\eta)mu^{\mu}$. As a result,

[TABLE]

where $\xi=\partial_{t}$ is the timelike Killing vector. On the other hand,

[TABLE]

Therefore we arrive at $D_{\Psi}^{\parallel}(\eta E)=(D_{\Psi}^{\parallel}\eta)(-mu_{t})$. In the similar way, we can derive $D_{\Psi}^{\parallel}(\eta L)=(D_{\Psi}^{\parallel}\eta)(mu_{\phi})$.

Appendix B Derivation of the MHD-GS equation (31)

We now expand Eq.(30) in terms of $\Psi(r,\theta),\Omega(\Psi),\eta(\Psi),\mathcal{M}(r,\theta),u_{t}(r,\theta),u_{\phi}(r,\theta)$, where the electromagnetic force have been derived in many previous FFE studies

[TABLE]

where $\mathcal{K}(r,\theta;\Omega)=-k_{0}$. With the aid of the normalization condition $u^{\mu}u_{\mu}=-1$ and $(u^{\mu}u_{\mu})_{;A}=0$, the matter acceleration is rewritten as

[TABLE]