Analytic properties of force-free jets in the Kerr spacetime -- III: uniform field solution

TL;DR

This paper investigates the uniqueness of solutions for force-free magnetospheres around Kerr black holes, revealing a unique solution with a discontinuity and current sheet when proper boundary and constraint conditions are applied.

Contribution

It clarifies the relationship between boundary conditions and solution uniqueness, providing a numerical method to find a unique force-free field configuration around Kerr black holes.

Findings

Multiple solutions approach uniform field at infinity

Numerical solution shows a unique configuration with a current sheet

Constraint and boundary conditions are interconnected via radiation conditions

Abstract

The structure of steady axisymmetric force-free magnetosphere of a Kerr black hole (BH) is governed by a second-order partial differential equation of $A_\phi$ depending on two "free" functions $\Omega(A_\phi)$ and $I(A_\phi)$, where $A_\phi$ is the $\phi$ component of the vector potential of the electromagnetic field, $\Omega$ is the angular velocity of the magnetic field lines and $I$ is the poloidal electric current. In this paper, we investigate the solution uniqueness. Taking asymptotically uniform field as an example, analytic studies imply that there are infinitely many solutions approaching uniform field at infinity, while only a unique one is found in general relativistic magnetohydrodynamic simulations. To settle down the disagreement, we reinvestigate the structure of the governing equation and numerically solve it with given constraint condition and boundary condition. We find that the constraint condition (field lines smoothly crossing the light surface (LS)) and boundary conditions at horizon and at infinity are connected via radiation conditions at horizon and at infinity, rather than being independent. With appropriate constraint condition and boundary condition, we numerically solve the governing equation and find a unique solution. Contrary to naive expectation, our numerical solution yields a discontinuity in the angular velocity of the field lines and a current sheet along the last field line crossing the event horizon. We also briefly discuss the applicability of the perturbation approach to solving the governing equation.

Equations

$$-$$

$$\left[\frac{\beta }{\Sigma }{\Omega }^2{\sin }^2\theta -\frac{4ra}{\Sigma }\Omega {\sin }^2\theta -\left(1-\frac{2r}{\Sigma }\right)\right]A_{\varphi ,rr}$$

$$+\left[\frac{\beta }{\Sigma }{\Omega }^2{\sin }^2\theta -\frac{4ra}{\Sigma }\Omega {\sin }^2\theta -\left(1-\frac{2r}{\Sigma }\right)\right]\frac{{\sin }^2\theta }{\Delta }{A}_{\varphi ,\mu \mu }$$

$$+\left[{\Omega }^2{\sin }^2\theta {\left(\frac{\beta }{\Sigma }\right)}_{,r}-\Omega {\sin }^2\theta {\left(\frac{4ra}{\Sigma }\right)}_{,r}+{\left(\frac{2r}{\Sigma }\right)}_{,r}\right]A_{\varphi ,r}$$

$$+\left[{\Omega }^2{\left(\frac{\beta {\sin }^2\theta }{\Sigma }\right)}_{,\mu }-\Omega {\left(\frac{4ra{\sin }^2\theta }{\Sigma }\right)}_{,\mu }+{\left(\frac{2r}{\Sigma }\right)}_{,\mu }\right]\frac{{\sin }^2\theta }{\Delta }{A}_{\varphi ,\mu }$$

$$+\left[\left(\frac{\beta }{\Sigma }\right)\Omega -\frac{2ra}{\Sigma }\right]{\sin }^2\theta {\Omega }^{\prime }\left(A_{\varphi ,r}^2+\frac{{\sin }^2\theta }{\Delta }{A}_{\varphi ,\mu }^2\right)$$

$$-\frac{\Sigma }{\Delta }I{I}^{\prime }=0,$$

$$\left[A_{\varphi ,rr}+\frac{{\sin }^2\theta }{\Delta }{A}_{\varphi ,\mu \mu }\right]\mathcal{K}(r,\theta ;\Omega )$$

$$+\left[A_{\varphi ,r}{\partial }_r^{\Omega }+\frac{{\sin }^2\theta }{\Delta }{A}_{\varphi ,\mu }{\partial }_{\mu }^{\Omega }\right]\mathcal{K}(r,\theta ;\Omega )$$

$$+\frac{1}{2}\left[A_{\varphi ,r}^2+\frac{{\sin }^2\theta }{\Delta }{A}_{\varphi ,\mu }^2\right]{\Omega }^{\prime }{\partial }_{\Omega }\mathcal{K}(r,\theta ;\Omega )$$

$$-\frac{\Sigma }{\Delta }I{I}^{\prime }=0,$$

$$\left[A_{\varphi ,r}{\partial }_r^{\Omega }+\frac{{\sin }^2\theta }{\Delta }{A}_{\varphi ,\mu }{\partial }_{\mu }^{\Omega }\right]\mathcal{K}(r,\theta ;\Omega )$$

$$+\frac{1}{2}\left[A_{\varphi ,r}^2+\frac{{\sin }^2\theta }{\Delta }{A}_{\varphi ,\mu }^2\right]{\Omega }^{\prime }{\partial }_{\Omega }\mathcal{K}(r,\theta ;\Omega )=\frac{\Sigma }{\Delta }I{I}^{\prime }.$$

$$I=\frac{2r(\Omega-\Omega_{\rm H})\sin^{2}\theta}{\Sigma}A_{\phi,\mu}\Big{|}_{r=r_{+}},$$

$$I=-\Omega\sin^{2}\theta A_{\phi,\mu}\Big{|}_{r\rightarrow\infty},$$

$$A_{\varphi }(r\to \infty )=r^2{\sin }^2\theta$$

$$I=2\Omega A_{\varphi }.$$

$$A_{\varphi }(\mu =1)$$

$$A_{\varphi }\left(\mu =0,r_{+}\le r\le 2\right)$$

$$A_{\varphi ,\mu }\left(\mu =0,r\ge 2\right)$$

$$A_{\varphi }$$

$$\Omega$$

$$I$$

$$I{I}_{\mathrm{new}}^{\prime }(A_{\varphi ,\mathrm{new}})$$

$$-0.02(A_{\varphi }(r_{\mathrm{ILS}}^{+})-A_{\varphi }(r_{\mathrm{ILS}}^{-})),$$

$$A_{\varphi ,\mathrm{new}}=0.5(A_{\varphi }(r_{\mathrm{ILS}}^{+})+A_{\varphi }(r_{\mathrm{ILS}}^{-})).$$

$$-{\int }_0^{A_{\varphi }^{\mathrm{H}}}I{I}^{\prime }d{A}_{\varphi }\delta (A_{\varphi }-A_{\varphi }^{\mathrm{H}}).$$

$$\mathcal{L}{A}_0=0,$$

$$\mathcal{L}\equiv \frac{\partial }{\partial r}\left(1-\frac{2}{r}\right)\frac{\partial }{\partial r}+\frac{{\sin }^2\theta }{r^2}\frac{{\partial }^2}{\partial {\mu }^2}.$$

$$A_{\varphi }=A_0+a^2{A}_2,\omega =a{\omega }_1,i=a{i}_1.$$

$$\mathcal{L}{A}_2(r,\theta )=S_2(r,\theta ;i_1,{\omega }_1),$$

$$A_2(r,\theta )=\int \int d{r}_{0}d{\theta }_0{S}_2(r_0,{\theta }_0)G(r,\theta ;r_0,{\theta }_0),$$

$$\mathcal{L}G(r,\theta ;r_0,{\theta }_0)=\delta (r-r_0)\delta (\theta -{\theta }_0).$$

$$\Omega =I=0$$

$$\Omega =\frac{{\Omega }_{\mathrm{H}}}{2}\frac{\sqrt{1-A_{\varphi }/A_{\varphi }^{\mathrm{H}}}}{1+\sqrt{1-A_{\varphi }/A_{\varphi }^{\mathrm{H}}}},$$

$$A_0(ϵ)=r^2{\sin }^2\theta +ϵ(1-\cos \theta ),$$

$$\mathcal{I}=\frac{2r(\Omega-\Omega_{\rm H})\sin^{2}\theta}{\Sigma}\mathcal{A}_{,\mu}\Big{|}_{r=r_{+}}$$

$$\frac{{\mathcal{A}}_{,\mu }}{f(\mathcal{A})}=-\frac{r_{+}{\sin }^2\theta }{r_{+}^2+a^2{\mu }^2},$$

Full text

††slugcomment: Submitted to ApJ

1. Introduction

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 better understanding of the BZ mechanism from general relativistic magnetohydrodynamic (GRMHD) simulations (e.g. Komissarov, 2001, 2004a, 2004b, 2005; Semenov et al., 2004; McKinney & Gammie, 2004; McKinney, 2005; McKinney & Narayan, 2007a, b; Komissarov & McKinney, 2007; Tchekhovskoy et al., 2008, 2010, 2011; Palenzuela et al., 2011; Alic et al., 2012; Tchekhovskoy & McKinney, 2012; Penna et al., 2013; McKinney 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; Pan & Yu, 2014, 2015a, 2015b, 2016; Gralla & Jacobson, 2014; Gralla et al., 2015, 2016b; Yang et al., 2015; Penna, 2015) to the steady axisymmetric force-free electrodynamics in the Kerr spacetime.111A few families of exact solutions (e.g. Menon & Dermer, 2005, 2007, 2011; Brennan et al., 2013; Menon, 2015; Compère et al., 2016) to the equations of force-free electrodynamics in the Kerr spacetime have been found in the past decade. But these solutions have various limitations, e.g., not allowing energy extraction from the BH, being electrically dominated instead of magnetically dominated, lacking clear physical interpretation, or being time-dependent and not axisymmetric, which make it difficult to compare these exact solutions with simulations and numerical solutions. Various studies converge to a common picture of how the BZ mechanism operates: The spinning BH distorts the poloidal magnetic field $\bm{B_{\rm P}}$, and induces the poloidal electric field $\bm{E_{\rm P}}$ and toroidal magnetic field $\bm{B_{\rm T}}$, which generate an outward Poynting flux $\bm{E_{\rm P}\times 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).

To 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 $\Omega$, 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 which roughly agrees with GRMHD simulations. But in this paper, we will show there are actually many of them due to the superposition of 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 $\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 exists usually 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 to explain the relationship between the unique solution and the infinitely many possible candidates, and the discrepancy between previous numerical solution and GRMHD simulations? 222It is definitely worthwhile examining the solution uniqueness problem of the BZ mechanism, considering its significance in the modern relativistic astrophysics. For example, rotating BHs are believed to be described by Kerr solution, but there is no solid evidence for it. Some efforts have been done to detect possible deviation from Kerr solution via BZ mechanism powered jet emission (e.g. Bambi, 2012a, b, 2015; Pei et al., 2016). A prerequisite for these BZ applications is its solution uniqueness. 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 the MHD Penrose process was operating in GRMHD simulations (see Lasota et al., 2014; Koide & Baba, 2014; Toma & Takahara, 2014; Kojima, 2015; Toma & Takahara, 2016, e.g. 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 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.333The stability problem of general force-free jets is a complicated story. Early studies implied that the force-free jets are vulnerable to various instabilities (e.g. Begelman, 1998; Lyubarskii, 1999; Li, 2000; Wang et al., 2004) while later studies showed that these instabilities are strongly suppressed by field rotation, poloidal field curvature, etc (e.g. Tomimatsu et al., 2001; McKinney & Blandford, 2009; Narayan et al., 2009, and references therein). But the relevant timescale is much longer than the transient time scales 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 $\Omega$ and $I$, ensuring both smooth field lines across the LS and 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 analytic perturbation approach to the GS equation, which relies on a fixed unperturbed solution and priorly known asymptotic behavior of magnetic field. Analytic approach breaks down if any of the two factors is violated. Both of the two are satisfied for 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; and our numerical method to solve the GS equation. We apply the perturbation approach on the uniform field problem and clarify the applicability of analytic perturbation approach in Section 4. Summary and Discussions are given in Section 5. In Appendix, we present a robust solver for the horizon regularity condition and its implication for the existence of electric current.

2. basic equations

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

[TABLE]

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

[TABLE]

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

[TABLE]

where $\mathcal{K}(r,\theta;\Omega)$ is the prefactor of $A_{\phi,rr}$ in Equation (2), $\partial_{i}^{\Omega}(i=r,\mu)$ denotes the partial derivative with respect to coordinate $i$ with $\Omega$ fixed, and $\partial_{\Omega}$ is the derivative with respect to $\Omega$. The GS equation written in this compact form manifests clear symmetry, therefore is beneficial in various aspects.444 For example, in flat spacetime, the classical pulsar equation (Scharlemann & Wagoner, 1973) is recovered by plugging $\mathcal{K}(r,\theta;\Omega)=\Omega^{2}r^{2}\sin^{2}\theta-1$ into Equation (2).

3. the solution uniqueness problem

In this section, we first clarify all the constraint conditions the GS equation satisfies, and their relation with boundary conditions at horizon and at infinity. We find the constraint conditions and boundary conditions are not independent. For a given $\Omega(A_{\phi})$, we can numerically find a $I(A_{\phi})$ ensuring field lines smoothly crossing the LS, but the combination of $\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 this boundary condition, $\Omega(A_{\phi})$ and $I(A_{\phi})$ must satisfy one more constraint. Then, we numerically find the unique combination of $\Omega(A_{\phi})$ and $I(A_{\phi})$ ensuring field lines smoothly cross the LS and being 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

[TABLE]

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

[TABLE]

and

[TABLE]

where $\Omega_{\rm H}$ is 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 horizon if $\Omega$ 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 $\Omega$ and $I$ are specified; or the radiation condition (6) enforces a constraint on $\Omega$ and $I$, if the boundary condition at infinity is given. In our working example, the boundary condition at infinity

[TABLE]

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

[TABLE]

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

Now we get two constraint conditions (4, 8), and two boundary conditions in the $r$ direction (5, 7) ready, (where the inner boundary condition (5) is nontrivial, see Appendix for details). The next step is to specify proper boundary conditions in the $\mu$ 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”, 555The claim depends on the relation $I\propto\Omega$, which can be derived from the radiation condition at infinity [Equation (6)]. But their are some debates about whether the radiation condition holds for vertical field configurations. As for the uniform field, there is no disagreement on the relation $I\propto\Omega$ [Equation (8)], though Pan & Yu (2016) interpreted it as the result of radiation condition, while Nathanail & Contopoulos (2014) interpreted differently. the possible field configuration in the steady state is shown in the Figure 1 of paper II. Consequently, we write boundary conditions in the $\mu$ direction as follows,

[TABLE]

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

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\times 64$ grid. The detailed numerical steps are as follows:

1. We choose some initial guess $A_{\phi}$, trial functions $\Omega$ and $I$ as follows, 666The convergent solution is independent of the trial field configuration or the grid resolution. For example, we tested different initial field configurations $A_{\phi}=r^{2}\sin^{2}\theta+\epsilon(1-\cos\theta)$, different initial trial functions $\Omega$ and $I$, and different grid resolutions.

[TABLE]

2. We evolve $A_{\phi}$ using relaxation method (Press et al., 1987), while this method does not work properly at the LS due to the vanishing second order derivatives. Fortunately, the directional derivative of $A_{\phi}$ is known as a function of $II^{\prime}$ there (see Equation (4)). We instead update $A_{\phi}$ 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_{\rm ILS}^{-})$/$A_{\phi}(r_{\rm 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_{\rm ILS})$ as follows:

[TABLE]

where

[TABLE]

Usually, $II^{\prime}$ obtained is not very smooth. To refrain possible numerical instabilities, we fit $II^{\prime}(A_{\phi})$ with a 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)

[TABLE]

In our compuation, 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 for $10$ times, then update $\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 $\Omega(A_{\phi})$ and $II^{\prime}(A_{\phi})$. From the plot, we see that the angular velocity of the last field line crossing the event horizon $\Omega(A_{\phi}^{\rm H})$ is not vanishing, i.e., $\Omega(A_{\phi}^{\rm H})\simeq 0.28\Omega_{\rm H}$, while we expect the angular velocity of field lines not crossing the BH vanishes, i.e., $\Omega(A_{\phi}>A_{\phi}^{\rm H})=0$.

3.3. Comparison with Previous Studies

Nathanail & Contopoulos (2014) also studied BH magnetosphere structure of uniform field and concluded that both $\Omega(A_{\phi})$ and $I(A_{\phi})$ must approach zero along the last field line crossing the event horizon, and therefore the ILS coincides with the IRS along the equator, and the electric current sheet does not appear. But as shown in the Appendix, the horizon regularity condition requires the existence of current sheet along the equator. And our numerical solution shows that there is a discontinuity in $\Omega(A_{\phi})$ at $A_{\phi}^{\rm H}$, and therefore the ILS lies between the event horizon and the IRS, and there exists a current sheet. The discrepancies here can be settled down by previous GRMHD simulations done by Komissarov (2005), where they observed a sharp transition in $\Omega(A_{\phi})$ at $A_{\phi}^{\rm H}$ and interpreted it as a discontinuity smeared by numerical viscosity. It is worth noting that the discontinuity in $\Omega(A_{\phi})$ does not lead to any physical difficulties, e.g., the continuity of $\bm{B}^{2}-\bm{E}^{2}$ across the LS is not affected.

In addition, we do not explicitly show the BZ power of the uniform field configuration here, because it is sensitive to the magnetic flux trapped by the event horizon $A_{\phi}^{\rm H}$ which is boundary condition dependent. In the real astrophysical environment, it is mainly determined by accretion process of the central BH (e.g. Garofalo, 2009).

4. Applicability of Analytic Perturbation Approach

In this section, we first recap the analytic perturbation approach to the GS equation, 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 a unperturbed solution $A_{0}$ in the Schwarzschild spacetime,

[TABLE]

where the operator

[TABLE]

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

[TABLE]

Then we linearize the GS equation (2) as

[TABLE]

by only keeping leading order perturbation terms, where the source function $S_{2}$ depends on $i_{1}$ 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

[TABLE]

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

[TABLE]

In this way, for a given Schwarzschild metric solution $A_{0}$, the corresponding Kerr metric solution $\{A_{\phi}(a;,r,\theta),I(A_{\phi}),\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),

[TABLE]

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 uniform field at infinity, but this is not the case. The Schwarzschild spacetime GS equation (14) is linear. Both uniform field and monopole field are its solutions, so do their linear combinations

[TABLE]

where $\epsilon$ is some constant coefficient. The mixture of monopole component generates a family of Schwarzschild metric solutions, $A_{0}(\epsilon)$, and all these solutions approach uniform field at infinity. For each solution, the corresponding Kerr metric solution $\{A_{\phi}(\epsilon),I(\epsilon),\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\rightarrow\infty}$ and the fixed underlying unperturbed field configuration $A_{0}$. 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.

5. Summary and Discussions

The GS equation is a second-order differential equation with two to-be-determined functions $\Omega$ and $I$. Generally speaking, we need two constraint conditions to determine $\Omega$ and $I$, two boundary conditions in the $r$ directions and two boundary conditions in the $\mu$ direction to fix $A_{\phi}(r,\mu)$. For asymptotically uniform field, we use constraint conditions (4,8), boundary conditions in the $r$ direction (5,7) and boundary conditions in the $\mu$ direction (9) to close the GS equation. Our numerical solution of the uniform field yields a discontinuity in the $\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 other two well studied field configurations: monopole field in the Kerr spacetime and 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 $\Omega$. For monopole field in the Kerr spacetime, there are two LS crossing conditions and two radiation conditions. The former two determine $\Omega$ and $I$, and the latter two determine the inner and the outer boundary. 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 pulsar dipole field, $\Omega$ is equal to the angular velocity of the central star and $I$ is the only function to determine. The only one 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 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\rightarrow\infty}$ and the fixed underlying unperturbed field configuration $A_{0}$. For the monopole field, both of them are satisfied, therefore the perturbation approach is applicable and the high-order perturbation solutions show good match with results from simulations (Pan & Yu, 2015a). For both uniform field in the Schwarzschild spacetime and dipole field surrounding static stars, the superposition of 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 radiation condition for vertical field configurations and on the possible difficulties arising form 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 Founda- tion of China (grants 11173057, 11373064, 11521303), Yunnan Natural Science Foundation (grants 2012FB187, 2014HB048), and the Youth Innovation Promotion Association, CAS. Part of the computation was performed at the HPC Center, Yunnan Observatories, CAS, China. L.H. thanks the support by the National Natural Science Foundation of China (grants 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 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.