Analytic Properties of Force-free Jets in the Kerr Spacetime. III. Uniform Field Solution

The Blandford–Znajek (BZ) mechanism (Blandford & Znajek 1977) is believed to be one of most efficient ways to extract rotation energy from spinning black holes (BHs), which operates in BH systems on all mass scales, from the stellar-mass BHs of gamma-ray bursts to the supermassive BHs of active galactic nuclei. In the past decade, we have gained a better understanding of the BZ mechanism from general relativistic magnetohydrodynamic (GRMHD) simulations (e.g., Komissarov 2001, 2004a, 2004b, 2005; McKinney & Gammie 2004; Semenov et al. 2004; McKinney 2005; Komissarov & McKinney 2007; McKinney & Narayan 2007a, 2007b; Tchekhovskoy et al. 2008, 2010, 2011; Nagataki 2009, 2011; Palenzuela et al. 2011; Alic et al. 2012; Tchekhovskoy & McKinney 2012; McKinney et al. 2013; Penna et al. 2013), numerical solutions (e.g., Fendt 1997; Uzdensky 2004, 2005; Palenzuela et al. 2010; Contopoulos et al. 2013; Nathanail & Contopoulos 2014), and analytic perturbation solutions (e.g., Tanabe & Nagataki 2008; Beskin & Zheltoukhov 2013; Gralla & Jacobson 2014; Pan & Yu 2014, 2015a, 2015b, 2016; Gralla et al. 2015, 2016b; Penna 2015; Yang et al. 2015) to the steady axisymmetric force-free electrodynamics in the Kerr spacetime.⁵ Various studies converge to a common picture of how the BZ mechanism operates: the spinning BH distorts the poloidal magnetic field ${{\boldsymbol{B}}}_{{\rm{P}}}$ and induces the poloidal electric field ${{\boldsymbol{E}}}_{{\rm{P}}}$ and toroidal magnetic field ${{\boldsymbol{B}}}_{{\rm{T}}}$, which generate an outward Poynting flux ${{\boldsymbol{E}}}_{{\rm{P}}}\times {{\boldsymbol{B}}}_{{\rm{T}}}$ along the magnetic field lines threading the spinning BH. The rotation energy of the spinning BHs is extracted in the form of Poynting flux (Komissarov 2009; Beskin 2010).

Going a step further, it is natural to ask whether these different approaches give qualitatively and quantitatively consistent descriptions of the BH magnetosphere structure, e.g., the topology of magnetic fields, the electric current distributions, the angular velocities of the magnetic field lines, and the energy extraction rates. The answer is yes and no. The axisymmetric, steady-state, force-free magnetosphere around Kerr BHs is governed by the general relativistic Grad–Shafranov (GS) equation. For the simplest magnetic field configuration, split monopole field, both analytic (Pan & Yu 2015a) and numerical solutions (Nathanail & Contopoulos 2014) reproduce the simulated angular velocity of field lines Ω, poloidal electric current I, and energy extraction rate $\dot{E}$ to high precision (Tchekhovskoy et al. 2010). But for the asymptotically uniform field, different approaches do not even reach a consensus on the solution uniqueness. Time-dependent simulations (e.g., Komissarov 2005; Komissarov & McKinney 2007; Yang et al. 2015) seem to converge to a unique solution. Previous analytic studies (Beskin & Zheltoukhov 2013; Pan & Yu 2014; Gralla et al. 2016b) seem to find a unique perturbation solution that roughly agrees with GRMHD simulations. But in this paper, we will show that there are actually many of them, due to the superposition of the monopole component (and other possible components). According to the argument of Nathanail & Contopoulos (2014), solving the GS equation is actually an eigenvalue problem, with two eigenvalues ${\rm{\Omega }}({A}_{\phi })$ and $I({A}_{\phi })$ to be determined by requiring field lines to smoothly cross the light surfaces (LSs). For common field configurations, there usually exist two LSs, sufficing to determine two eigenvalues. With only one LS for the uniform field configuration and one more boundary condition, Nathanail & Contopoulos (2014) numerically found a unique solution, which, however, shows distinctive features from previous GRMHD simulations (Komissarov 2005).

How can the relationship between the unique solution and the infinitely many possible candidates, and the discrepancy between previous numerical solutions and GRMHD simulations, be explained?⁶ Does the plasma inertia make a difference? The force-free condition is assumed in both analytic and numerical solutions, but the inertia cannot be completely ignored in simulations. Taking account of the plasma inertia, Takahashi et al. (1990) proposed the so-called MHD Penrose process, where the plasma particles within the ergosphere are projected onto negative-energy orbits by magnetic field and eventually are captured by the central BH. As a result, Alfvèn waves are generated along the magnetic field lines, and BH rotation energy is carried away by these Alfvèn waves. Koide et al. (2002) and Koide (2003) found that the MHD Penrose process was operating in GRMHD simulations (see e.g., Koide & Baba 2014; Lasota et al. 2014; Toma & Takahara 2014, 2016; Kojima 2015, for recent discussions on this issue). If the MHD Penrose process is the dominant energy extraction process, the unique solution found in simulations actually describes the MHD Penrose process instead of the BZ mechanism. However, later simulations showed that the MHD Penrose process is only a transient state, after which the Alfvèn waves decay, the system settles down into a steady state, and the BZ mechanism takes over (e.g., Komissarov 2005). Therefore, the plasma inertia seems to make little difference after the system settles down into the steady state.

Another possible explanation is that, among all these mathematically possible solutions, only the one found in simulations is stable. Yang & Zhang (2014) and Yang et al. (2015) analyzed the stability of these solutions, and no unstable mode was found at order O(a), where a is the dimensionless BH spin. Therefore, modes can be unstable with a growth rate at $\sim O({a}^{2})$ at most.⁷ But the relevant timescale is much longer than the transient timescales observed in simulations. Therefore, they concluded that the selection rule unlikely comes from instability.

In this paper, we show that the uniform field solution is unique, as strongly implied by previous GRMHD simulations and pointed out by Nathanail & Contopoulos (2014). Following the algorithm proposed by Contopoulos et al. (2013) and Nathanail & Contopoulos (2014), we numerically find a unique combination of Ω and I, ensuring both smooth field lines across the LS and a uniform field at infinity. Contrary to Nathanail & Contopoulos (2014), our numerical solution yields a discontinuity in the angular velocity of field lines and a current sheet along the last field line crossing the event horizon, which are features found in previous simulations.

We also investigate the applicability of the analytic perturbation approach to the GS equation, which relies on a fixed unperturbed solution and the priorly known asymptotic behavior of the magnetic field. The analytic approach breaks down if either of the two factors is violated. Both of them are satisfied for the monopole field in Kerr spacetime; therefore, we see the perfect match between high-order perturbation solutions and results from simulations and numerical solutions. But for the uniform field, the unperturbed background field is not fixed, due to the superposition of the monopole component; therefore, the perturbation approach cannot predict a unique solution.

The paper is organized as follows. In Section 2, we summarize the basic equations governing the steady axisymmetric force-free magnetospheres. In Section 3, we clarify the relation between constraint conditions, radiation conditions, and boundary conditions, as well as our numerical method to solve the GS equation. We apply the perturbation approach on the uniform field problem and clarify the applicability of the analytic perturbation approach in Section 4. Summary and discussions are given in Section 5. In the Appendix, we present a robust solver for the horizon regularity condition and its implication for the existence of electric current.

In the force-free approximation, electromagnetic energy greatly exceeds that of matter. Consequently, the force-free magnetosphere is governed by the energy conservation equation of the electromagnetic field, conventionally called the GS equation. In the Kerr spacetime, the axisymmetric and steady GS equation is written as (Pan & Yu 2014)

$$\begin{eqnarray} & & -{\rm{\Omega }}[{(\sqrt{-g}{F}^{{tr}})}_{,r}+{(\sqrt{-g}{F}^{t\theta })}_{,\theta }]+{F}_{r\theta }I^{\prime} ({A}_{\phi })\\ & & +[{(\sqrt{-g}{F}^{\phi r})}_{,r}+{(\sqrt{-g}{F}^{\phi \theta })}_{,\theta }]=0,\end{eqnarray} \tag{ 1 }$$

which expands as (see also, e.g., Contopoulos et al. 2013; Nathanail & Contopoulos 2014; Pan & Yu 2016, in slightly different forms)

$$\begin{eqnarray}&&\left[\displaystyle \frac{\beta }{{\rm{\Sigma }}}{{\rm{\Omega }}}^{2}{\sin }^{2}\theta -\displaystyle \frac{4{ra}}{{\rm{\Sigma }}}{\rm{\Omega }}{\sin }^{2}\theta -\left(1-\displaystyle \frac{2r}{{\rm{\Sigma }}}\right)\right]{A}_{\phi ,{rr}}\\ &&+\left[\displaystyle \frac{\beta }{{\rm{\Sigma }}}{{\rm{\Omega }}}^{2}{\sin }^{2}\theta -\displaystyle \frac{4{ra}}{{\rm{\Sigma }}}{\rm{\Omega }}{\sin }^{2}\theta -\left(1-\displaystyle \frac{2r}{{\rm{\Sigma }}}\right)\right]\displaystyle \frac{{\sin }^{2}\theta }{{\rm{\Delta }}}{A}_{\phi ,\mu \mu }\\ &&+\left[{{\rm{\Omega }}}^{2}{\sin }^{2}\theta {\left(\displaystyle \frac{\beta }{{\rm{\Sigma }}}\right)}_{,r}-{\rm{\Omega }}{\sin }^{2}\theta {\left(\displaystyle \frac{4{ra}}{{\rm{\Sigma }}}\right)}_{,r}+{\left(\displaystyle \frac{2r}{{\rm{\Sigma }}}\right)}_{,r}\right]{A}_{\phi ,r}\\ &&+\left[{{\rm{\Omega }}}^{2}{\left(\displaystyle \frac{\beta {\sin }^{2}\theta }{{\rm{\Sigma }}}\right)}_{,\mu }-{\rm{\Omega }}{\left(\displaystyle \frac{4{ra}{\sin }^{2}\theta }{{\rm{\Sigma }}}\right)}_{,\mu }+{\left(\displaystyle \frac{2r}{{\rm{\Sigma }}}\right)}_{,\mu }\right]\displaystyle \frac{{\sin }^{2}\theta }{{\rm{\Delta }}}{A}_{\phi ,\mu }\\ &&+\left[\left(\displaystyle \frac{\beta }{{\rm{\Sigma }}}\right){\rm{\Omega }}-\displaystyle \frac{2{ra}}{{\rm{\Sigma }}}\right]{\sin }^{2}\theta {\rm{\Omega }}^{\prime} \left({A}_{\phi ,r}^{2}+\displaystyle \frac{{\sin }^{2}\theta }{{\rm{\Delta }}}{A}_{\phi ,\mu }^{2}\right)\\ &&-\displaystyle \frac{{\rm{\Sigma }}}{{\rm{\Delta }}}{II}^{\prime} =0,\end{eqnarray} \tag{ 2 }$$

where ${\rm{\Sigma }}={r}^{2}+{a}^{2}\,{\mu }^{2}$, ${\rm{\Delta }}={r}^{2}-2r+{a}^{2}$, $\beta ={\rm{\Delta }}{\rm{\Sigma }}\ +2r({r}^{2}+{a}^{2})$, $\mu \equiv \cos \theta$, and the primes designate derivatives with respect to A_ϕ. For clarity, we may write the GS equation in a more illustrating form.

$$\begin{eqnarray}&&\left[{A}_{\phi ,{rr}}+\displaystyle \frac{{\sin }^{2}\theta }{{\rm{\Delta }}}{A}_{\phi ,\mu \mu }\right]{ \mathcal K }(r,\theta ;{\rm{\Omega }})\\ &&+\left[{A}_{\phi ,r}{\partial }_{r}^{{\rm{\Omega }}}+\displaystyle \frac{{\sin }^{2}\theta }{{\rm{\Delta }}}{A}_{\phi ,\mu }{\partial }_{\mu }^{{\rm{\Omega }}}\right]{ \mathcal K }(r,\theta ;{\rm{\Omega }})\\ &&+\displaystyle \frac{1}{2}\left[{A}_{\phi ,r}^{2}+\displaystyle \frac{{\sin }^{2}\theta }{{\rm{\Delta }}}{A}_{\phi ,\mu }^{2}\right]{\rm{\Omega }}^{\prime} {\partial }_{{\rm{\Omega }}}{ \mathcal K }(r,\theta ;{\rm{\Omega }})\\ &&-\displaystyle \frac{{\rm{\Sigma }}}{{\rm{\Delta }}}{II}^{\prime} =0,\end{eqnarray} \tag{ 3 }$$

where ${ \mathcal K }(r,\theta ;{\rm{\Omega }})$ is the prefactor of ${A}_{\phi ,{rr}}$ in Equation (2), ${\partial }_{i}^{{\rm{\Omega }}}(i=r,\mu )$ denotes the partial derivative with respect to coordinate i with Ω fixed, and ${\partial }_{{\rm{\Omega }}}$ is the derivative with respect to Ω. The GS equation written in this compact form manifests clear symmetry and therefore is beneficial in various aspects.⁸

In this section, we first clarify all the constraint conditions the GS equation satisfies, as well as their relation with boundary conditions at the horizon and at infinity. We find that the constraint conditions and boundary conditions are not independent. For a given ${\rm{\Omega }}({A}_{\phi })$, we can numerically find an $I({A}_{\phi })$ ensuring field lines smoothly crossing the LS, but the combination of ${\rm{\Omega }}({A}_{\phi })$ and $I({A}_{\phi })$ obtained this way usually is in conflict with the uniform field boundary condition at infinity. To be consistent with this boundary condition, ${\rm{\Omega }}({A}_{\phi })$ and $I({A}_{\phi })$ must satisfy one more constraint. Then, we numerically find the unique combination of ${\rm{\Omega }}({A}_{\phi })$ and $I({A}_{\phi })$ that ensures that field lines smoothly cross the LS and is consistent with the boundary condition at infinity. Finally, we compare our numerical solution with previous studies.

3.1. Constraint Conditions and Boundary Conditions

We want physically allowed solutions to be finite and smooth everywhere. At LS where ${ \mathcal K }=0$, the second-order GS equation degrades to a first-order equation

$$\begin{eqnarray}&&\left[{A}_{\phi ,r}{\partial }_{r}^{{\rm{\Omega }}}+\displaystyle \frac{{\sin }^{2}\theta }{{\rm{\Delta }}}{A}_{\phi ,\mu }{\partial }_{\mu }^{{\rm{\Omega }}}\right]{ \mathcal K }(r,\theta ;{\rm{\Omega }})\\ &&+\displaystyle \frac{1}{2}\left[{A}_{\phi ,r}^{2}+\displaystyle \frac{{\sin }^{2}\theta }{{\rm{\Delta }}}{A}_{\phi ,\mu }^{2}\right]{\rm{\Omega }}^{\prime} {\partial }_{{\rm{\Omega }}}{ \mathcal K }(r,\theta ;{\rm{\Omega }})=\displaystyle \frac{{\rm{\Sigma }}}{{\rm{\Delta }}}{II}^{\prime} .\end{eqnarray} \tag{ 4 }$$

Field lines smoothly crossing the LS must satisfy the above constraint, which we call the LS crossing constraint condition. At the horizon and infinity, the requirement of solution finiteness leads to the radiation conditions (e.g., Pan & Yu 2016), which read as

$$\begin{eqnarray}&&{I=\displaystyle \frac{2r({\rm{\Omega }}-{{\rm{\Omega }}}_{{\rm{H}}}){\sin }^{2}\theta }{{\rm{\Sigma }}}{A}_{\phi ,\mu }\right|}_{r={r}_{+}}\end{eqnarray} \tag{ 5 }$$

and

$$\begin{eqnarray}&&{I=-{\rm{\Omega }}{\sin }^{2}\theta {A}_{\phi ,\mu }| }_{r\to \infty },\end{eqnarray} \tag{ 6 }$$

where ${{\rm{\Omega }}}_{{\rm{H}}}$ is the angular velocity of the central BH.

But the radiation conditions and boundary conditions are not independent. For example, the radiation condition (5) uniquely determines the boundary values at the horizon if Ω and I are specified, and we will use it as the inner boundary condition in our numerical calculation. In the same way, the radiation condition (6) uniquely determines the boundary values at infinity, if Ω and I are specified; or the radiation condition (6) enforces a constraint on Ω and I, if the boundary condition at infinity is given. In our working example, the boundary condition at infinity

$$\begin{eqnarray}&&{A}_{\phi }(r\to \infty )={r}^{2}{\sin }^{2}\theta \end{eqnarray} \tag{ 7 }$$

is given. Plugging it into the radiation condition (6), we find that Ω and I must satisfy a new constraint (Nathanail & Contopoulos 2014; Pan & Yu 2014, 2016),

$$\begin{eqnarray}&&I=2{\rm{\Omega }}{A}_{\phi }.\end{eqnarray} \tag{ 8 }$$

Note that conditions (6), (7), and (8) are not independent, and we will use two of them, conditions (7) and (8), to close the GS equation.

Now we get two constraint conditions (4) and (8) and two boundary conditions in the r direction (5) and (7) ready, (where the inner boundary condition (5) is nontrivial; see the Appendix for details). The next step is to specify proper boundary conditions in the μ direction. According to the claim proved in Paper II, "In the steady axisymmetric force-free magnetosphere around a Kerr BH, all magnetic field lines that cross the infinite-redshift surface must intersect the event horizon."⁹ The possible field configuration in the steady state is shown in Figure 1 of Paper II. Consequently, we write boundary conditions in the μ direction as follows:

$$\begin{eqnarray}{A}_{\phi }(\mu =1) & = & 0,\\ {A}_{\phi }(\mu =0,{r}_{+}\leqslant r\leqslant 2) & = & {A}_{\phi }^{{\rm{H}}},\\ {A}_{\phi ,\mu }(\mu =0,r\geqslant 2) & = & 0,\end{eqnarray} \tag{ 9 }$$

where the horizon enclosed magnetic flux ${A}_{\phi }^{{\rm{H}}}$ is to be determined self-consistently.

Zoom In Zoom Out Reset image size

3.2. Numerical Method and Results

The algorithm for numerically solving the GS equation was proposed by Contopoulos et al. (2013) and was optimized by Nathanail & Contopoulos (2014). We slightly tailor their algorithm to accommodate the problem we are working on. We define a new radial coordinate $R=r/(1+r)$, confine our computation domain $R\times \mu$ in the region $[R({r}_{+}),1]\times [0,1]$, and implement a uniform 512 × 64 grid. The detailed numerical steps are as follows.

1. We choose some initial guess A_ϕ and trial functions Ω and I as follows:¹⁰

$$\begin{eqnarray}{A}_{\phi } & = & {r}^{2}{\sin }^{2}\theta ,\\ {\rm{\Omega }} & = & \displaystyle \frac{{{\rm{\Omega }}}_{{\rm{H}}}}{2}\cos \left(\displaystyle \frac{\pi }{2}\displaystyle \frac{{A}_{\phi }}{{A}_{\phi }^{{\rm{H}}}}\right),\\ I & = & {{\rm{\Omega }}}_{{\rm{H}}}{A}_{\phi }\cos \left(\displaystyle \frac{\pi }{2}\displaystyle \frac{{A}_{\phi }}{{A}_{\phi }^{{\rm{H}}}}\right).\end{eqnarray} \tag{ 10 }$$

2. We evolve A_ϕ using the relaxation method (Press et al. 1987), although this method does not work properly at the LS owing to the vanishing second-order derivatives. Fortunately, the directional derivative of A_ϕ is known as a function of ${II}^{\prime}$ there (see Equation (4)). We instead update A_ϕ at the LS using neighborhood grid points and the directional derivative. From the directional derivative and the grid points on the left/right side, we obtain ${A}_{\phi }({r}_{\mathrm{ILS}}^{-})$/${A}_{\phi }({r}_{\mathrm{ILS}}^{+})$. Usually the two are not equal and field lines are broken here. To smooth the field lines, we adjust $I({A}_{\phi })$ and update ${A}_{\phi }({r}_{\mathrm{ILS}})$ as follows:

$$\begin{eqnarray}{{II}}_{\mathrm{new}}^{\prime }({A}_{\phi ,\mathrm{new}}) & = & {{II}}_{\mathrm{old}}^{\prime }({A}_{\phi ,\mathrm{old}})\\ & & -0.02({A}_{\phi }({r}_{\mathrm{ILS}}^{+})-{A}_{\phi }({r}_{\mathrm{ILS}}^{-})),\end{eqnarray} \tag{ 11 }$$

where

$$\begin{eqnarray}&&{A}_{\phi ,\mathrm{new}}=0.5({A}_{\phi }({r}_{\mathrm{ILS}}^{+})+{A}_{\phi }({r}_{\mathrm{ILS}}^{-})).\end{eqnarray} \tag{ 12 }$$

Usually, ${II}^{\prime}$ obtained is not very smooth. To avoid possible numerical instabilities, we fit ${II}^{\prime} ({A}_{\phi })$ with eighth-order polynomials. In addition, ${II}^{\prime}$ consists of two pieces: a regular piece determined as described above, and a singular piece (the current sheet part)

$$\begin{eqnarray}&&-{\int }_{0}^{{A}_{\phi }^{{\rm{H}}}}{II}^{\prime} {{dA}}_{\phi }\ \delta ({A}_{\phi }-{A}_{\phi }^{{\rm{H}}}).\end{eqnarray} \tag{ 13 }$$

In our computation, we model the delta function as a parabola confined in a finite interval $[{A}_{\phi }^{{\rm{H}}},{A}_{\phi }^{{\rm{H}}}(1+\delta )]$ with $\delta =0.1$ (see e.g., Gruzinov 2005).

3. Repeat step 2 10 times, and then update ${\rm{\Omega }}({A}_{\phi })$ according to the constraint (8).

We iterate the initial guess solution following the above steps until field lines smoothly cross the LS and satisfy the boundary conditions. The numerical results are shown in Figure 1. In the left panel, we show the convergent field configuration, which, as expected, matches those of simulations (e.g., Komissarov & McKinney 2007). In the right panel, we show functions ${\rm{\Omega }}({A}_{\phi })$ and ${II}^{\prime} ({A}_{\phi })$. From this plot, we see that the angular velocity of the last field line crossing the event horizon ${\rm{\Omega }}({A}_{\phi }^{{\rm{H}}})$ is not vanishing, i.e., ${\rm{\Omega }}({A}_{\phi }^{{\rm{H}}})\simeq 0.28{{\rm{\Omega }}}_{{\rm{H}}}$, while we expect the angular velocity of field lines not crossing the BH to vanish, i.e., ${\rm{\Omega }}({A}_{\phi }\gt {A}_{\phi }^{{\rm{H}}})=0$.

3.3. Comparison with Previous Studies

In this section, we first recap the analytic perturbation approach to the GS equation and then apply it to the uniform field problem (Pan & Yu 2014, 2015a) and explain why this approach actually yields many solutions. Finally, we discuss the applicability of the perturbation method.

We start with an unperturbed solution A₀ in the Schwarzschild spacetime,

$$\begin{eqnarray}&&{ \mathcal L }{A}_{0}=0,\end{eqnarray} \tag{ 14 }$$

where the operator

$$\begin{eqnarray}&&{ \mathcal L }\equiv \displaystyle \frac{\partial }{\partial r}\left(1-\displaystyle \frac{2}{r}\right)\displaystyle \frac{\partial }{\partial r}+\displaystyle \frac{{\sin }^{2}\theta }{{r}^{2}}\displaystyle \frac{{\partial }^{2}}{\partial {\mu }^{2}}.\end{eqnarray} \tag{ 15 }$$

For the corresponding Kerr metric solution $\{{A}_{\phi },I({A}_{\phi }),{\rm{\Omega }}({A}_{\phi })\}$, we assume ${A}_{\phi }{| }_{r\to \infty }={A}_{0}$, define $i\ =I{| }_{r\to \infty },\omega ={\rm{\Omega }}{| }_{r\to \infty }$, and expand the solution to the leading order

$$\begin{eqnarray}&&{A}_{\phi }={A}_{0}+{a}^{2}{A}_{2},\quad \omega =a{\omega }_{1},\quad i={{ai}}_{1}.\end{eqnarray} \tag{ 16 }$$

Then we linearize the GS Equation (2) as

$$\begin{eqnarray}&&{ \mathcal L }{A}_{2}(r,\theta )={S}_{2}(r,\theta ;{i}_{1},{\omega }_{1}),\end{eqnarray} \tag{ 17 }$$

by only keeping leading-order perturbation terms, where the source function S₂ depends on i₁ and ${\omega }_{1}$, which can be figured out from the radiation conditions at horizon and at infinity (5)–(6). The solution to the linearized GS equation is written as

$$\begin{eqnarray}&&{A}_{2}(r,\theta )=\int \int {{dr}}_{0}d{\theta }_{0}\ {S}_{2}({r}_{0},{\theta }_{0})G(r,\theta ;{r}_{0},{\theta }_{0}),\end{eqnarray} \tag{ 18 }$$

where $G(r,\theta ;{r}_{0},{\theta }_{0})$ is the Green's function of operator ${ \mathcal L }$ (Petterson 1974; Blandford & Znajek 1977),

$$\begin{eqnarray}&&{ \mathcal L }G(r,\theta ;{r}_{0},{\theta }_{0})=\delta (r-{r}_{0})\delta (\theta -{\theta }_{0}).\end{eqnarray} \tag{ 19 }$$

In this way, for a given Schwarzschild metric solution A₀, the corresponding Kerr metric solution $\{{A}_{\phi }(a;,\,r,\theta ),I({A}_{\phi }),{\rm{\Omega }}({A}_{\phi })\}$ is uniquely determined order by order. Applying the method on the uniform field problem ${A}_{0}={r}^{2}{\sin }^{2}\theta$, we find (Beskin & Zheltoukhov 2013; Pan & Yu 2014; Gralla et al. 2015)

$$\begin{eqnarray}{\rm{\Omega }} & = & I=0\quad ({A}_{\phi }\gt {A}_{\phi }^{{\rm{H}}}),\\ {\rm{\Omega }} & = & \displaystyle \frac{{{\rm{\Omega }}}_{{\rm{H}}}}{2}\displaystyle \frac{\sqrt{1-{A}_{\phi }/{A}_{\phi }^{{\rm{H}}}}}{1+\sqrt{1-{A}_{\phi }/{A}_{\phi }^{{\rm{H}}}}},\ I=2{\rm{\Omega }}{A}_{\phi }\quad ({A}_{\phi }\lt {A}_{\phi }^{{\rm{H}}}),\end{eqnarray} \tag{ 20 }$$

where ${A}_{\phi }^{{\rm{H}}}$ is the magnetic flux trapped by the event horizon (${A}_{\phi }^{{\rm{H}}}=4$ for the lowest-order perturbation solution and generally depends on BH spins and boundary conditions).

It seems that we have found the unique solution (at leading order) approaching a uniform field at infinity, but this is not the case. The Schwarzschild spacetime GS Equation (14) is linear. Both the uniform field and monopole field are its solutions, so do their linear combinations are also its solutions,

$$\begin{eqnarray}&&{A}_{0}(\epsilon )={r}^{2}{\sin }^{2}\theta +\epsilon (1-\cos \theta ),\end{eqnarray} \tag{ 21 }$$

where is some constant coefficient. The mixture of the monopole component generates a family of Schwarzschild metric solutions, ${A}_{0}(\epsilon )$, and all these solutions approach a uniform field at infinity. For each solution, the corresponding Kerr metric solution $\{{A}_{\phi }(\epsilon ),I(\epsilon ),{\rm{\Omega }}(\epsilon )\}$ can be obtained using the above perturbation method.

To summarize, the perturbation method depends on two main ingredients: the known asymptotic field behavior at infinity ${A}_{\phi }{| }_{r\to \infty }$ and the fixed underlying unperturbed field configuration A₀. But for the uniform field problem, there are many mathematically allowed unperturbed solutions due to the additional monopole component. That's why the perturbation approach cannot predict the unique solution.

The GS equation is a second-order differential equation with two to-be-determined functions Ω and I. Generally speaking, we need two constraint conditions to determine Ω and I, two boundary conditions in the r directions, and two boundary conditions in the μ direction to fix ${A}_{\phi }(r,\mu )$. For an asymptotically uniform field, we use constraint conditions (4) and (8), boundary conditions in the r direction (5) and (7), and boundary conditions in the μ direction (9) to close the GS equation. Our numerical solution of the uniform field yields a discontinuity in the ${\rm{\Omega }}({A}_{\phi })$ at ${A}_{\phi }^{{\rm{H}}}$; therefore, the ILS lies between the event horizon and the IRS, and there exists a current sheet along the last field line crossing the event horizon (Figure 1), which is as expected from the horizon regularity condition.

Following the same logic, let us reexamine two other well-studied field configurations: the monopole field in the Kerr spacetime and the dipole field in the flat spacetime (classical pulsars). For both field configurations, the number of LSs equals the number of to-be-determined functions I and Ω. For the monopole field in the Kerr spacetime, there are two LS crossing conditions and two radiation conditions. The former two determine Ω and I, and the latter two determine the inner and the outer boundaries. Hence, there is no more freedom for specifying a boundary condition at infinity, i.e., we actually do not know the solution at infinity before we really solve the GS equation. Previous simulations and numerical solutions indeed confirmed the asymptotic monopole field configuration ${A}_{\phi }\propto 1-\cos \theta$. For the pulsar dipole field, Ω is equal to the angular velocity of the central star and I is the only function to determine. The only LS uniquely determines I (Contopoulos et al. 1999), and two radiation conditions automatically determine boundary conditions. In the same way, there is no more freedom to impose a boundary condition at infinity. And previous numerical studies found that the field at infinity deviates from ${A}_{\phi }\propto 1-\cos \theta$ (e.g., Gralla et al. 2016a).

We also discuss the perturbation approach for solving the GS equation, whose applicability depends on two main ingredients: the known asymptotic field behavior at infinity ${A}_{\phi }{| }_{r\to \infty }$ and the fixed underlying unperturbed field configuration A₀. For the monopole field, both of them are satisfied; therefore, the perturbation approach is applicable and the high-order perturbation solutions show a good match with results from simulations (Pan & Yu 2015a). For both the uniform field in the Schwarzschild spacetime and the dipole field surrounding static stars, the superposition of the monopole component (and possible other components) generates many unperturbed solutions; as a result, the perturbation approach cannot predict a unique uniform field solution in the Kerr spacetime or a unique dipole field solution surrounding spinning stars.

We thank our referee Prof. I. Contopoulos for his insightful comments on the validity of the radiation condition for vertical field configurations and on the possible difficulties arising from the discontinuity in the angular velocity of field lines. We also thank A. Nathanail for clearly explaining his numerical algorithm. C.Y. is grateful for the support by the National Natural Science Foundation of China (grants 11173057, 11373064, 11521303), the Yunnan Natural Science Foundation (grants 2012FB187 and 2014HB048), and the Youth Innovation Promotion Association, CAS. Part of the computation was performed at the HPC Center, Yunnan Observatories, CAS, China. L.H. acknowledges the support by the National Natural Science Foundation of China (grant 11203055). This paper is supported in part by the Strategic Priority Research Program "The Emergence of Cosmological Structures" of the Chinese Academy of Sciences, grant no. XDB09000000, and the Key Laboratory for Radio Astronomy, CAS. This work made extensive use of the NASA Astrophysics Data System and of the astro-ph preprint archive at arXiv.org.

For a given combination of $I({A}_{\phi })$ and ${\rm{\Omega }}({A}_{\phi })$, the horizon regularity condition (Equation (5)) uniquely determines the boundary condition at horizon ${A}_{\phi }(r={r}_{+},\mu )$. We find that its numerical solution is not trivial, due to the nonlinearity. To construct a robust solver, we first rewrite Equation (5) as

$$\begin{eqnarray}&&{{ \mathcal I }=\displaystyle \frac{2r({\rm{\Omega }}-{{\rm{\Omega }}}_{{\rm{H}}}){\sin }^{2}\theta }{{\rm{\Sigma }}}{{ \mathcal A }}_{,\mu }\right|}_{r={r}_{+}},\end{eqnarray} \tag{ 22 }$$

where we have defined two normalized variables, ${ \mathcal I }=I/{A}_{\phi }^{{\rm{H}}}$ and ${ \mathcal A }={A}_{\phi }/{A}_{\phi }^{{\rm{H}}}$. Here ${ \mathcal A }$ runs from 0 to 1, and its values on boundaries are ${ \mathcal A }(\mu =0)=1$ and ${ \mathcal A }(\mu =1)=0$. Furthermore, we define $f({ \mathcal A })\equiv { \mathcal I }/2({{\rm{\Omega }}}_{{\rm{H}}}-{\rm{\Omega }})$, and the above equation is written in a variable separated form

$$\begin{eqnarray}&&\displaystyle \frac{{{ \mathcal A }}_{,\mu }}{f({ \mathcal A })}=-\displaystyle \frac{{r}_{+}{\sin }^{2}\theta }{{r}_{+}^{2}+{a}^{2}\,{\mu }^{2}},\end{eqnarray} \tag{ 23 }$$

which has a formal solution

$$\begin{eqnarray}&&{e}^{{\int }_{1}^{{ \mathcal A }(\mu )}}\tfrac{d{ \mathcal A }}{f({ \mathcal A })}=\displaystyle \frac{1-\mu }{1+\mu }\times {e}^{\tfrac{{a}^{2}}{{r}_{+}}\mu }.\end{eqnarray} \tag{ 24 }$$

In this form, numerically solving ${ \mathcal A }(\mu )$ is stable.

Here we show a general property of force-free magnetospheres read out of the formal solution (24): the horizon condition requires the existence of a current sheet at the equator. Enabling a nonsingular solution, the integral ${\int }_{1}^{{ \mathcal A }(\mu )}d{ \mathcal A }/f({ \mathcal A })$ must be finite, except at ${ \mathcal A }=0(\mu =1)$, where $f({ \mathcal A })=0$ owing to vanishing ${ \mathcal I }(\mu =1)$. At ${ \mathcal A }=1\,(\mu =0)$, the finite integral requires nonzero $f({ \mathcal A }=1)$, or quickly decreased $f({ \mathcal A })$ (e.g., $\sim \sqrt{1-{ \mathcal A }}$). Usually ${\rm{\Omega }}\leqslant {{\rm{\Omega }}}_{{\rm{H}}}/2$; therefore, $I({A}_{\phi }^{{\rm{H}}})$ must be nonzero, or quickly decrease to zero (e.g., $\sim \sqrt{{A}_{\phi }^{{\rm{H}}}-{A}_{\phi }}$), where the former is a current sheet at the equator and the latter is a divergent current density (but weaker than the former).