The Force-Free Magnetosphere of a Rotating Black Hole

We revisit the Blandford&Znajek (1977) process and solve the fundamental equation that governs the structure of the steady-state force-free magnetosphere around a Kerr black hole. The solution depends on the distributions of the magnetic field angular velocity omega and the poloidal electric current I. These are not arbitrary. They are determined self-consistently by requiring that magnetic field lines cross smoothly the two singular surfaces of the problem, the inner `light surface' located inside the ergosphere, and the outer `light surface' which is the generalization of the pulsar light cylinder. We find the solution for the simplest possible magnetic field configuration, the split monopole, through a numerical iterative relaxation method analogous to the one that yields the structure of the steady-state axisymmetric force-free pulsar magnetosphere (Contopoulos, Kazanas&Fendt 1999). We obtain the rate of electromagnetic extraction of energy and confirm the results of Blandford and Znajek and of previous time dependent simulations. Furthermore, we discuss the physical applicability of magnetic field configurations that do not cross both `light surfaces'.


INTRODUCTION
The electromagnetic extraction of energy from a spinning black hole can be one of the most powerful and efficient engines in the Universe. It was Blandford & Znajek (1977) (hereafter BZ77) who first argued that a spinning black hole threaded by a sufficiently strong magnetic field and permeated by an electron-positron plasma generated from pair cascades, establishes a force-free magnetosphere, in direct analogy to the theory of the axisymmetric pulsar magnetosphere developed a few years earlier by Goldreich & Julian (1969). They derived the fundamental equation that governs the structure of this magnetosphere, the general relativistic force-free Grad-Shafranov equation, and obtained a second order perturbative solution for small values of the Kerr parameter a and a scaling for the electromagnetic power extracted from the rotating black hole.
The structure, singular surfaces and general properties of this equation for magnetic fields threading Schwarzschild and Kerr black holes were discussed recently by Uzdensky (2004Uzdensky ( , 2005. He emphasized the need for a self-consistent determination of the distributions of the magnetic field angular velocity and the poloidal electric current, something that early relaxation methods (Macdonald 1984) and later higher order perturbative solutions (Tanabe & Nagataki 2008) did not take into account. He also noted that this may be a very difficult task, and indeed, to the best of our knowledge a solution of this equation in the general case of an open magnetic field configuration threading a Kerr black hole is still lacking. Nevertheless, the broader astrophysical interest in the Blandford-Znajek process, has prompted in the past decade several researchers to seek and obtain the steady-state structure of the spinning black hole magnetospheres through time-dependent general relativistic force-free numerical simulations (e.g. Komissarov & McKinney 2007, Palenzuela et al. 2011, Lyutikov & McKinney 2011. These works established the universal applicability of the electromagnetic extraction of the rotation energy of a Kerr black hole in powering active galactic nuclei, jets in X-ray binaries, and even gamma-ray bursts. The presence of this extra source of energy in addition to that of the accretion disk may account for some of the phenomenology. One should bear in mind, however, that the situation is more complicated. For instance: (1) Almost 90% of extragalactic sources are radio quiet, i.e. they do not produce electromagnetic outflows in the form of radio jets, and yet they are believed to contain both constituents of the Blandford-Znajek process, namely a spinning black hole and a strong magnetic field threading it (e.g. Kukula 2003); (2) The formation and disruption of radio jets in black hole X-ray binaries over short time scales is at odds with this general notion (e.g. Belloni 2010); (3) There seems to exist strong observational evidence that the power in the radio jet is not directly related to black hole spin, contrary to the main result of BZ77 (Fender, Gallo & Russell 2010; however, see also Narayan & McClintock 2012). Notice, though, that this may be due to an erroneous determination of the black hole spin.
The above justify an independent re-evaluation of the Blanford-Znajek process through the solution of the general relativistic force-free Grad-Shafranov equation. The special relativistic form of that equation, the so called pulsar equation, was only solved in 1999 by Contopoulos, Kazanas and Fendt who showed that the magnetic field structure can only be obtained together with the determination of the unique electric current distribution that allows for a smooth and continuous structure to exist. In that respect, steady-state solutions offer a deeper insight in the physics of the problem than numerical timedependent ones. It should be noted that, while BZ77 provided an estimate of the energy extraction efficiency of this process, they did not obtain an exact solution of the structure of the black hole magnetosphere. Also, at that time, the critical role of the singular surfaces and the distributions of the magnetospheric electric current and field line angular velocity had not been appreciated. Today, we revisit this problem with all the knowledge we carry from our 13-year long investigation of the force-free pulsar magnetosphere.

THE GENERAL RELATIVISTIC PULSAR EQUATION
In order to derive the fundamental equation that governs the steady-state structure of the force-free magnetosphere around a Kerr black hole we follow closely the 3+1 formulation of Thorne & Macdonald (1982) used by most researchers in the astrophysical community (e.g. Uzdensky 2005). We restrict our analysis to steady-state and axisymmetric space times where (. . .) ,t = (. . .) ,φ = 0. In that case, the general 4-dimensional space-time geometry may be written in Boyer-Lindquist spherical coordinates x µ ≡ (t, r, θ, φ) as Here, α ≡ (∆Σ/A) 1/2 and Ω ≡ 2aM r/A are the lapse function and angular velocity of 'zero-angular momentum' fiducial observers (ZAMOs) respectively, ∆ ≡ r 2 − 2M r + a 2 , Σ ≡ r 2 + a 2 cos 2 θ , M is the black hole mass, and a its angular momentum (0 ≤ a ≤ M ). Throughout this paper we adopt geometric units where G = c = 1. Semicolon stands for covariant derivative, comma for partial derivative. Latin indices denote spatial components (1 − 3), Greek indices denote space-time components (0 − 3), and '∼' denotes the spatial part of vectors. ZAMOs move with 4-velocity U µ = (1/α , 0 , 0 , Ω/α) orthogonal to hypersurfaces of constant t. The force-free magnetosphere of a spinning black hole is characterized by the electromagnetic energymomentum tensor and the condition T µν ;ν = 0. Here, the rest mass and pressure contribution have been neglected. The electromagnetic field tensor F µν is related to the electric and magnetic fields E µ , B µ measured by ZAMOs through F µν = U µ E ν − U ν E µ + ǫ µνλρ B λ U ρ (ǫ µνλρ ≡ [µνλρ]|det(g µν )| −1/2 is the 4-dimensional Levi-Civita tensor). Under these conditions, the fundamental equation that governs the steady-state structure of the forcefree magnetosphere around a Kerr black hole becomes ρ e andJ are the electric charge and current densities respectively. Eq. (3) is supplemented by Maxwell's equations of electrodynamics∇ For several applications in astrophysics, perfect (infinite) conductivity is a valid approximation. In this case, and the electric and magnetic vector fields can be expressed in terms of three scalar functions, Ψ(r, θ), ω(Ψ), and I(Ψ) as ω is the angular velocity of the magnetic field lines, I is the poloidal electric current flowing through the circular loop r =const., θ =const., and Ψ is equal to (2π) −1 times the total magnetic flux enclosed in that loop. Notice that the electric field changes sign close to the horizon with respect to its sign at large distances. As explained in BZ77, a rotating observer (ZAMO) will in general see a Poynting flux of energy entering the horizon, but he will also see a sufficiently strong flux of angular momentum leaving the horizon. That ensures that energy is extracted from the black hole. The poloidal component of Eq.
(3) then yields the general relativistic force-free Grad-Shafranov equation (Eq. (3.14) of BZ77 re-written in our notation). Henceforth, primes will denote differentiation with respect to Ψ. One sees directly that if we set α = 0 and M = 0 in eq.
(2) we obtain which is the well known pulsar equation (Scharlemann & Wagoner 1973). The zeroing of the term multiplying the second order derivatives in eq. (2), yields the singular surfaces of the problem. We will henceforth call the singular surfaces 'light surfaces' (LS) 1 . When we set M = α = 0, it yields the standard pulsar light cylinder r sin θ = c/ω. In general, the shape of the outer LS is only asymptotically cylindrical as θ → 0 (see figure 2 below), and an inner LS appears inside the ergosphere. We will henceforth use the notation 'light cylinder' (LC) and 'outer LS' interchangeably. It is interesting to note that the outer boundary of the ergosphere corresponds to the solution of the singularity condition for ω = 0, whereas the inner boundary (the event horizon at r = r BH ≡ M + √ M 2 − a 2 ) corresponds to the solution of the singularity condition for ω = Ω BH ≡ a/(r 2 BH + a 2 ), where Ω BH is the angular velocity of the black hole. It is also interesting to note that the natural 'radiation condition' at infinity (energy must flow outwards along all field lines) requires that (BZ77), and therefore indeed the inner LS lies inside the ergosphere. Both Eqs.
(2) and (9) contain the two functions, ω(Ψ) and I(Ψ), which must be determined by the physics of the problem. In the case of an axisymmetric spinning neutron star, ω is usually taken to be equal to the neutron star angular velocity Ω NS 2 . In pulsars, I(Ψ) is self-consistently determined through an iterative numerical technique that implements a smooth crossing of the relativistic Alfvèn surface, which is the usual light cylinder where r sin θ = c/ω (Contopoulos, Kazanas & Fendt 1999, Timokhin 2006). In the case of a spinning black hole, the situation is qualitatively similar but more complicated. Contrary to a neutron star, the black hole does not have a solid surface, and therefore it cannot impose any restriction on ω other than the 'radiation condition' Eq. (11). The only natural restriction that is imposed by the physical problem is that magnetic field lines must be smooth and continuous everywhere. ω(Ψ) must therefore be determined together with I(Ψ) through the condition of smooth crossing of both LS of the problem, the inner one inside the ergosphere, and the outer one at large distances. In the next section, we will see that this can be achieved through an iterative numerical technique analogous to the one employed in the pulsar magnetosphere. As remarked previously, iterating with respect to two functions simultaneously is indeed a very difficult task. This is why Uzdensky (2005) opted to consider only field lines that connect the black hole horizon to the surrounding accretion disk where ω(Ψ) is determined by the Keplerian angular velocity of the disk. There exists one more interesting complication. Eq. (11) applies only to open field lines that cross the event horizon. In the case of a star, field lines that cross the stellar surface rotate with the angular velocity of the star, but they do not cross an event horizon (either because no event horizon forms, or because they simply avoid the horizon), and the above restriction does not apply. In the present case, we are interested in studying the fundamentals of the Blandford-Znajek process, and therefore, we will consider a black hole 'as clean as possible'. Obviously, a magnetized astrophysical black hole is not isolated. In that respect, we would like here to study an astrophysical black hole with the minimum number of extra elements, namely a surrounding thin disk of matter to hold the magnetic field through the necessary electric currents and charges. Charges and currents will be produced through pair production in the black hole magnetosphere. We will ignore magnetic field lines that do not cross both the inner and outer LS because we will have no physical way to determine their ω and the poloidal electric current I contained inside them. We will also ignore magnetic field lines that cross the equator outside the black hole event horizon because such field lines need a source of poloidal electric current on the equator (if they rotate, they will need to cross at least one LS where the condition of smooth crossing will in general require a nonzero electric current to flow along that line).
Therefore, in the present work we will only study black hole magnetospheres where all magnetic field lines reach the event horizon, and all magnetic field lines are open. As argued in Lyutikov (2012), such configuration is the most natural end stage of stellar collapse. Notice once again that it is sustained by electric currents in an equatorial current sheet. The effect of a surrounding disk of a certain height possibly threaded by the return magnetic field (as in the scenario predicted by the Cosmic Battery; Contopoulos & Kazanas 1998) will be considered in a future work.

THE NUMERICAL ITERATIVE SCHEME
We will here describe the technical numerical details that led to the solution of Eq. (2). Firstly, we change the radial variable from r to R(r) ≡ r/(r + M ). Our numerical integration extends from R min ≡ R(r BH ) (the event horizon) to R max = 1 (radial infinity). The θ coordinate extends from θ min = 0 (the axis of symmetry) to θ max = π/2 (the equatorial plane). We implemented a 256 × 64 numerical grid uniform in R and θ. Notice that this grid has a very high resolution in r around the black hole event horizon where the inner LS lies, but not as high around the outer LS, and in particular at low θ and a values.
Secondly, we specify boundary conditions on the axis of symmetry, the horizon, the equatorial plane, and infinity. Ψ(θ = 0) = 0. Ψ(r BH , θ) is determined by the integration over θ of the so called Znajek regularization condition (Znajek 1977). Notice that this condition is updated as the distribution I(Ψ) is updated along with the numerical solution for Ψ(r, θ). The equatorial boundary condition is interesting. As we discussed in the previous section, the equatorial region may contain an accretion disk, with or without its own magnetic field. Before investigating such a complex astrophysical system, we would like to understand first the simplest case, that of a magnetized black hole where the electric currents supporting its magnetic field are distributed on a thin equatorial current sheet. In that case, the equatorial boundary condition becomes Ψ(r, θ = π/2) = Ψ max . The outer boundary condition is also not important since the outer boundary has practically been moved to infinity. In practice, we chose Ψ ,r (R = 1) = 0, but any other boundary condition yields the same results. We initialize our numerical grid with a split monopole field configuration with ω(Ψ) = 0.5Ω BH and The reader can directly check that Eqs. (13)  Thirdly, we update the distributions of ω(Ψ) and I(Ψ) as follows: At each latitude θ, we check where the singularity condition (Eq.10) is satisfied in r. At each such radial position, we extrapolate Ψ inwards from larger r (Ψ(r + , θ)) and outwards from smaller r (Ψ(r − , θ)). In general, Ψ(r + , θ) and Ψ(r − , θ) differ. Then, at the inner LS we implement whereas at the outer LS we implement where at each LS. The reasoning here is that we impose weighted corrections on ω(Ψ) and I(Ψ) based on the non-smoothness of the Ψ(r, θ) distribution along all grid points inside and outside the two LS, and not through the regularization conditions as we did in Contopoulos, Kazanas & Fendt (1999). This is a very general procedure that may be applied to any similar singular equation. Our numerical scheme does not converge easily, and its parameters need to be empirically adjusted by following the progress of the iteration. The signs of the weighting coefficients are determined by trial and error: when we choose the wrong signs, the mismatches |Ψ(r + , θ)−Ψ(r − , θ)| increase and the Ψ(r, θ) distribution is very quickly destroyed. Their magnitudes are also obtained empirically: too large values (around unity) lead to over-corrections and instability, too small values (below 0.01) do not correct fast enough for the iteration to converge. Furthermore, in order to facilitate convergence, at the inner LS we update only ω, whereas at the outer LS we update both I and ω every other relaxation iteration for Ψ. Also, in order to avoid numerical instabilities, we smooth out the distributions of I and ω every 50 relaxation iterations for Ψ. We found that more frequent smoothing inhibits convergence. Our goal is to minimize the residuals in the discretized form of eq. (2) together with the non-smoothness in Ψ(r, θ), but our method fails for poor grid resolutions at the inner and/or the outer LS (see Discussion section).
The results of our numerical integration are shown in figure 1. We show here the distributions of ω(Ψ) (top panel), 2|I(Ψ)| (center panel), and Poynting flux at large distances 2|I(Ψ)|ω(Ψ) (bottom panel) after 6000 relaxation steps for Ψ(R, θ), for six values of the black hole spin parameter a. As is the case in pulsars, the magnetospheric electric current is non-zero at Ψ = Ψ max , and the global electric circuit closes through an equatorial current sheet. Notice that we have no way to update ω along the axis, and therefore we have chosen ω(Ψ = 0) = 0.5Ω BH . 3 On the equator we implemented ω ′ (Ψ = Ψ max ) = I ′ (Ψ = Ψ max ) = 0. The reader should check the qualitative agreement with the results of the time-dependent force-free general relativistic simulations performed by Tchekhovskoy, Narayan & McKinney (2010) (their figure 4).
In figure 2 we show the characteristic magnetic field configuration for a/M = 0.9999 near the black hole horizon (left panel) and at larger distances (right panel). We also plot the inner LS inside the ergosphere, and the outer LS which becomes asymptotically cylindrical as θ → 0. Notice that magnetic field lines are monopolar at large distances (beyond a few times the radius of the event horizon), and they bend upwards towards the axis as they approach the horizon. This effect is also seen in the numerical simulations of Lyutikov & McKinney (2011) for a = 0.99. It is not discernable in the field configurations of figure 1 of Tchekhovskoy, Narayan & McKinney (2010) because the scale is not adequate, but is implied in their figure 7. Notice that ω(Ψ) is equal to 0.5Ω BH to within 10%, and that I(Ψ) is very close to the split monopole expression I SM (Ψ) also to within 10% (Eqs. 13 & 14). Therefore, the Znajek condition (Eq. 12) can be rewritten using Eq. (6) as Finally, we calculate the total rate of electromagnetic energy extraction (from both hemispheres) as measured at infinity through Eq. (4.11) of Blanford & Znajek (1977) as We have here normalized our results to the split monopole value with ω SM = 0.5Ω BH . k ≈ 0.8 for a/M = 0.9999 and is practically indistinguishable from unity for a/M < ∼ 0.9 (see bottom panel of figure 1), thus confirming the main result of BZ77, namely that, in an open magnetic field configuration with an infinitely thin surrounding equatorial disk, electromagnetic energy is extracted at a rate proportional to both Ω 2 BH and Ψ 2 max .

DISCUSSION
We believe that the study of the steady-state open force-free black hole magnetosphere offers a deeper insight in the physics of the problem than time-dependent numerical simulations. In analogy to the axisymmetric force-free pulsar magnetosphere, we argued that steadystate solutions can only be obtained through a selfconsistent determination of the distributions of both the magnetic field angular velocity and the magnetospheric electric current ω(Ψ) and I(Ψ) respectively. These can only be determined when open field lines cross both singular LS, the inner LS inside the ergosphere outside the black hole horizon, and the outer LS which is a deformed light cylinder. We developed an iterative numerical procedure analogous to the one we implemented in the study of the pulsar magnetosphere. Iterating with respect to two functions simultaneously is a very difficult task and our method may certainly be improved in the future. Yet, we are confident enough to suggest that the distributions are unique for the particular boundary conditions of the problem (all magnetic field lines thread the black hole horizon, all are open, and the surrounding disk is infinitely thin) since different initial conditions for Ψ, ω and I evolve towards the same final solution (for each value of a that is). Our solutions may serve as test cases for time-dependent numerical codes. We plan to apply our method in the study of various astrophysically interesting situations where the surrounding disk is not infinitely thin but expands vertically with distance.
All field lines that cross the black hole horizon are bound to cross the inner LS. This is not the case, though, for the outer LS. In fact, several numerical simulations exist in the literature where the black hole is threaded by a vertical uniform magnetic field, and yet they converge to certain distributions of ω(Ψ) and I(Ψ) (e.g. Komissarov & McKinney 2007, Palenzuela et al. 2011. We expect that in that case, the solution is mathematically indeterminate, and the solutions shown depend on the particular choice of outer numerical boundary conditions. For example, the lateral boundary conditions in Palenzuela et al. (2011) are not perfectly absorbing but reflect outgoing waves as is manifested in the oscillations seen in their figure 2. We plan to study such magnetic field configurations in a future work and check whether we can indeed freely specify one of the two free functions, ω(Ψ) or I(Ψ). We will thus check whether the results of previous time-dependent numerical simulations are unique or not.
We would like to caution the reader about certain numerical complications that arise in all numerical simulations of the spinning black hole magnetosphere, time dependent and steady state. Numerical simulations that do not adequately resolve both the inner and outer LS should not be trusted. We faced a similar problem in our study of the pulsar magnetosphere. In that case, there exists no inner LS, and in order to study the transition through the outer LS (the light cylinder) we were forced to implement our inner boundary condition at a radial distance more than a hundred times larger than the actual neutron star radius (Spitkovsky 2006, McKinney 2006, Contopoulos & Kalapotharakos 2010. This, however, posed no problem since the angular velocity of the field lines is set by the rotation of the pulsar. The black hole problem is more complicated because the black hole surface is not infinitely conducting, and the numerical code must adequately resolve both the inner boundary and the light cylinder. Notice that our numerical method works best as a/M approaches unity since in that case the inner LS is well separated from the event horizon, and the outer LS is at its closest approach to the origin of the numerical grid where its radial resolution is the highest. In that case, a numerical resolution of 256 radial grid points corresponds to 34 grid cells between the horizon and the inner LS, and 94 grid cells between the two LS along the equator. When a/M > ∼ 0.9, our results are rather resolution independent. On the other hand, as a/M approaches zero, the ergosphere shrinks in width and the inner LS approaches so close to the event horizon that our numerical grid does not have enough radial resolution for the iterative procedure that yields ω(Ψ) to properly work 4 . At the same time, as a/M decreases so does Ω BH , pushing the outer LS so far out that our radially expanding numerical grid has very low radial resolution at the position of the outer LS for the iterative procedure to properly converge. Similarly, the resolution of our radially expanding grid on the light cylinder becomes worse as θ → 0, and this too complicates the convergence of our iterative method. Other numerical approaches certainly face similar complications at short and/or at large distances. Our method may be improved by a reformulation in radially expanding cylindrical coordinates, in Kerr-Schild coordinates, or in standard spherical coordinates with a numerical integration outer radial boundary.
In summary, the force-free spinning black hole magnetosphere is not very different from the force-free axisymmetric pulsar magnetosphere. In fact, beyond a distance on the order of a few times the black hole radius, the two problems are almost identical. An important corollary of this similarity is that a spinning black hole surrounded by a thin equatorial plasma disk does not generate a collimated outflow, only an uncollimated wind (as in Contopoulos, Kazanas & Fendt 1999). Thus, the problem of the black hole jet formation remains open. A possible resolution must involve a surrounding thick disk, or a surrounding magnetic disk wind that gradually collimate the initially monopolar black hole wind into a collimated jet outflow, as seen in the numerical simulations of Tchekhovskoy, Narayan & McKinney (2010). It is interesting to note here, again similarly to pulsars, that beyond the light cylinder the Lorentz factor Γ in the monopolar wind increases linearly with cylindrical dis-tance as Γ ∝ r sin θ c/ω ≈ rΩ BH sin θ 2c (20) (Contopoulos & Kazanas 2002), possibly up to the fast magnetosonic point of the outflow. Therefore, according to Eq. (20), collimated field/flow lines that do not extend as far from the axis as uncollimated ones are expected to reach lower Lorentz factors. Of course, the final Lorentz factor of such outflows will also depend on the mass loading of a given magnetic field lines, but such considerations are beyond the scope of the present work. We plan to continue our investigation of the interrelation between the black hole and disk outflows in the future.