Effects of Inner Alfvén Surface Location on Black Hole Energy Extraction in the Limit of a Force-Free Magnetosphere
EEffects of Inner Alfv´en Surface Location on Black Hole Energy Extraction in theLimit of a Force-Free Magnetosphere
Kevin Thoelecke and Sachiko Tsuruta
Department of Physics, Montana State University, Bozeman, Montana 59717, USA
Masaaki Takahashi
Department of Physics and Astronomy, Aichi University of Education, Kariya, Aichi 448-8542, Japan
An energy extracting black hole magnetosphere can be defined by the location of its inner Alfv´ensurface, which determines the rate of black hole energy extraction along a given magnetic field line.We study how the location of the inner Alfv´en surface can modify the structure of energy extractingblack hole magnetospheres in the force-free limit. Hundreds of magnetospheres are numericallycomputed via a general relativistic extension of the Newtonian magnetofrictional method for a fullrange of black hole spins and flow parameters. We find that jet-like structures naturally form veryclose to the horizon for Alfv´en surfaces near the boundary of the ergosphere and that energy isextracted towards the equatorial plane for Alfv´en surfaces close to the horizon. This suggests thattwo broad classes of energy extracting black hole magnetospheres might exist; those that transmitextracted energy directly to distant observers, and those that transmit extracted energy to nearbyaccreting matter. Applied to transient high energy phenomena, we find that they might also differin timescale by a factor of 20 or more.
I. INTRODUCTION
When a rotating black hole is immersed in a magneticfield supported by a plasma, it is possible for that mag-netic field to carry an electromagnetic Poynting flux ofenergy and angular momentum away from the black holeto a distant observer. Conservation of energy and angu-lar momentum requires that the spin of the black hole bereduced; it is therefore said that such a magnetosphereextracts rotational energy from the black hole.This type of electromagnetic black hole energy extrac-tion is often referred to as the
Blandford-Znajek mecha-nism after Blandford and Znajek [1] who described it interms of a stationary, axisymmetric, and force-free blackhole magnetosphere. They showed that such a magneto-sphere can efficiently extract an enormous amount of en-ergy, and ever since black hole energy extraction has beenconsidered as a possible mechanism powering a range ofastrophysical phenomena, from stellar mass black holesin gamma-ray bursts to supermassive black holes in ac-tive galactic nuclei.The Blandford-Znajek mechanism is not the only wayto extract black hole rotational energy (Lasota et al. [2]reviews others), nor is it strictly necessary for magneticfield lines to cross the horizon in order to extract blackhole rotational energy (Komissarov [3] and referencestherein). However the simple idea of a single magneticfield line connecting the horizon to distant observers canbe seductive, as it is reminiscent of the more familiarand well understood case of a magnetic field line thread-ing and torquing a rotating body. As one might na¨ıvelyexpect from such an analogy, Blandford and Znajek [1]showed that for black hole energy to be extracted along agiven magnetic field line that field line must rotate moreslowly than the black hole. In Takahashi et al. [4] it wasdemonstrated that the condition of a relatively slowly ro- tating magnetic field line is actually the force-free limitof the requirement that a plasma inflow become super-Alfv´enic inside the ergosphere.Despite the defining role that the location of the in-ner Alfv´en surface plays in determining whether or nota given magnetosphere extracts black hole rotational en-ergy, it is currently unknown how the shape and locationof the Alfv´en surface might correspond to changes in thestructure of a magnetosphere. The differences (if any)between a magnetosphere with an Alfv´en surface nearthe horizon and a magnetosphere with an Alfv´en surfacenear the boundary of the ergosphere are largely uncer-tain. This is due to the extreme intractability of theequations governing plasma flows in rotating black holespacetimes, most especially the transfield equation thatdescribes the force balance transverse to magnetic fieldlines (Nitta et al. [5]). Currently no exact solutions ofmagnetospheres with plasma flows satisfying the trans-field equation are known, so there is no robust analyticmethod of studying any potential correspondences be-tween Alfv´en surface location and magnetosphere struc-ture.In the force-free limit plasma inertial effects are ig-nored, the inner Alfv´en surface coincides with the innerlight surface described by the field line angular velocityof the magnetosphere, and the transfield equation simpli-fies to a purely electromagnetic condition. Due to thosesimplifications it is natural to consider the limit of force-free magnetospheres as a first step in attempting a studyof the effects of Alfv´en surface location. Unfortunatelydespite its relative simplicity the force-free limit of thetransfield equation has very few analytic solutions, pri-marily limited to a handful of energy extracting mag-netospheres found by perturbing in black hole spin (e.g.Blandford and Znajek [1], McKinney and Gammie [6],Pan and Yu [7]) and a single class of fully exact but a r X i v : . [ a s t r o - ph . H E ] F e b non energy extracting (and mostly unphysical) solutions(Menon and Dermer [8] and Gralla and Jacobson [9]).Therefore even in the force-free limit there is also no ro-bust analytic method of studying the effects of Alfv´ensurface location.Current numerical solutions and simulations also lacksignificant utility, as the nature of the Alfv´en surface isultimately an unknown function of the specific choice ofboundary conditions and related initial assumptions thatare made. This means that the Alfv´en surface at best aresult to be reported rather than an input to be explored,making it difficult to study how differing Alfv´en surfacesmight correspond to changes in magnetosphere structure.The goal of this work was to study how Alfv´en sur-face location might correspond to the structure of energyextracting black hole magnetospheres, thereby beginningto fill a gap in our understanding of black hole energyextraction. Due to the lack of useful analytic solutions itwas necessary to calculate magnetospheres numerically.For simplicity we focused on the limit of force-free mag-netospheres. We also focused on magnetospheres witha boundary condition in the equatorial plane compat-ible with a monopolar geometry; this avoids the con-founding affects of the arbitrarily large forcings alter-native boundary conditions can exert on the magneto-sphere. This boundary condition is also compatible withthose used to calculate magnetospheres in Contopoulos et al. [10] and Nathanail and Contopoulos [11] as well asthe perturbed monopole solutions originally calculatedby Blandford and Znajek [1], allowing us to more easilycompare our results with those works.To form an initial survey it was necessary to calculatea large number of magnetospheres for a wide range ofpotential black hole spins and field line angular veloci-ties (ultimately ∼
500 magnetospheres were found), sowe developed numerical techniques that are efficient atcalculating magnetospheres with similar structures.Those numerical techniques are described in SectionIII. In Section II we provide some background on bothideal and force-free magnetospheres as well as the as-sumptions and equations we use. In Section IV we statesome of our results, which we discuss in more depth inSection V before concluding.
II. BACKGROUND AND ASSUMPTIONS
In the stationary and axisymmetric limit, force-freeblack hole magnetospheres are naturally described bythree scalar fields: the toroidal vector potential A φ ( r, θ ),the field line angular velocity Ω F ( A φ ), and the toroidalmagnetic field B φ ( A φ ). In order to form a valid solu-tion those scalar fields must satisfy a single differentialequation, the force-free limit of the general relativisticGrad-Shafranov (or transfield) equation. For black holerotational energy to be extracted along a given field linethat field line must rotate slower than the black hole;0 < Ω F < ω H , where ω H is the angular velocity of the horizon. In this section we provide a brief review ofthose results and the assumptions they rely on. Withthe primary exception of B φ , which we define with op-posite sign, the notation used is largely compatible withthat of Blandford and Znajek [1] and Takahashi et al. [4] to facilitate comparison with those works; significantdiscrepancies will be noted. Throughout we use covari-ant 4-vector notation; Thorne and MacDonald [12] andKomissarov [3] provide reviews of and translations to the3+1 formulations often used in other works. A. Core Assumptions
We begin by assuming a stationary, axisymmetric, un-charged, and isolated black hole such that the spacetimemay be described by the Kerr metric in Boyer-Lindquistcoordinates. Using a (+ , − , − , − ) metric signature andunits where c = G = 1, this corresponds to the line ele-ment: ds = (cid:18) − mr Σ (cid:19) dt + 4 mar sin θ Σ dtdφ − Σ∆ dr − Σ dθ − A sin θ Σ dφ , (1)where Σ = r + a cos θ, ∆ = r − mr + a ,A = (cid:0) r + a (cid:1) − ∆ a sin θ. (2)The black hole is then immersed in a perfectly con-ducting (ideal) plasma and electromagnetic fields that to-gether lack sufficient energy density to modify the back-ground spacetime. The plasma and electromagnetic fieldconfiguration is also assumed to be stationary and ax-isymmetric, with an axis of symmetry corresponding tothe rotational axis of the black hole. The electromagneticfields must be consistent with Maxwell’s equations in acoordinate basis: (cid:0) √− gF αβ (cid:1) ,β = − π √− gJ α ,F αβ,γ + F βγ,α + F γα,β = 0 . (3)Here F αβ is the field strength tensor, J α is the currentfour vector, g is the metric determinant, and a commadenotes a partial derivative. We use Greek letters to de-note spacetime indices ( α = t, r, θ, φ ), lowercase Latinletters to denote spatial indices ( a = r, θ, φ ), and upper-case Latin letters to denote poloidal indices ( A = r, θ ).The field strength tensor may be expressed in terms ofan electromagnetic vector potential A α as: F αβ = A β,α − A α,β . (4)In this form it is clear that the toroidal electric field van-ishes ( F tφ = 0) due to the assumptions of stationarityand axisymmetry, so “toroidal field” always references atoroidal magnetic field. The stress energy tensor outsidethe horizon is taken to be a linear combination of plasmaand electromagnetic effects: T αβ = T αβ Plasma + T αβ EM = ( ρ + p ) u α u β − pg αβ + 14 π g αµ F µλ F λβ + 116 π g αβ F µλ F µλ , (5)where ρ and p reference the proper energy density andpressure of the plasma, and u α its four velocity.We finally assume that there is no external forcing (en-ergy and momentum is conserved) over any region of in-terest such that the divergence of the stress energy tensorvanishes everywhere: T αβ ; β = 0 . (6)Everything that follows rests upon the above assump-tions and conditions. B. Field Aligned Conserved Quantities
The three scalar fields that describe a force-free magne-tosphere (vector potential A φ , field line angular velocityΩ F , and toroidal magnetic field B φ ) are field-aligned con-served quantities in the sense that for a scalar field Ψ wehave: B α Ψ ,α = 0 . (7)Here B α is the magnetic field measured by a distant ob-server, which we define as: B α = 12 √− g (cid:15) αβµν k β F µν , (8)where k β is the Killing vector associated with the sta-tionarity of the spacetime.The toroidal vector potential A φ directly describes thepoloidal magnetic field, making it a useful flux functionto trace poloidal magnetic field lines. The conservation of A φ rests upon stationarity and axisymmetry, so A φ is alsoa valid flux function in non force-free magnetospheres.The conservation of field line angular velocity Ω F is thegeneral relativistic extension of Ferraro’s law of isorota-tion. The field line angular velocity is defined in termsof the field strength tensor as: F rφ Ω F = F tr ,F θφ Ω F = F tθ . (9)This states that all electric fields are the result of therotation of a purely magnetic configuration with angu-lar velocity Ω F (referenced to zero angular momentumframes). The rigid rotation of individual magnetic fieldlines through the conservation of Ω F is a consequence of stationarity, axisymmetry, and a perfectly conduct-ing plasma, so it is also valid in non force-free magneto-spheres.The toroidal field B φ is defined through the definitionof the magnetic field B α in Equation 8 as B φ = √− gF θr .This has opposite sign to the definition in Blandfordand Znajek [1], so for the remainder of this work thetoroidal field will be referred to as √− gF θr to avoid con-fusion. The conservation of the toroidal field may beunderstood by considering Poynting fluxes through themagnetosphere. The only electric field is the result of ro-tating a purely magnetic configuration, and is thereforeboth entirely poloidal and perpendicular to the poloidalmagnetic field. Therefore all poloidal Poynting fluxes areweighted only by the toroidal magnetic field and alignedwith the poloidal magnetic field. As such the toroidalfield is naturally interpreted as a measure of a conservedflux of energy E and angular momentum L per unit fieldline: E = 14 π √− gF θr Ω F ,L = 14 π √− gF θr . (10)Note that these definitions differ from those in Takahashi et al. [4]; in that work they are scaled by a conserved par-ticle flux η as E → ηE and L → ηL . The conservationof the toroidal field along magnetic field lines is a con-sequence of stationarity, axisymmetry, and a perfectlyconducting force-free plasma; in non force-free magneto-spheres E and L include plasma parameters and √− gF θr is not conserved.Varied derivations and discussions of these and otherfield-aligned conserved quantities may be found in Bland-ford and Znajek [1], Bekenstein and Oron [13], Camen-zind [14], and Takahashi et al. [4]. C. Critical Points and Energy Extraction
In this section we explore key points along magneticfield lines and examine the conditions required for a mag-netosphere to extract black hole rotational energy. Webegin by making two assumptions that are unnecessaryin general but greatly simplify the discussion. First, weassume a single magnetic field line that extends directlyfrom the horizon to a distant observer. Second, we insistthat this field line be “simple” in the sense that the dis-tance from the horizon decreases monotonically as onetravels along the field line towards the black hole; a use-ful mental image is a straight line connecting the horizonwith distant regions.Along any such field line there are seven key points,the first being the separation point. Inside the separa-tion point gravitational forces dominate and draw plasmainwards to the horizon; outside the separation point cen-trifugal forces in the rotating frame of the magnetic fieldline dominate and accelerate plasma away from the blackhole. The remaining six points are the ingoing and out-going slow magnetosonic, Alfv´en, and fast magnetosonicpoints through which the plasma accelerates as it trav-els from the separation point to the horizon or distantregions. In the force-free limit the ingoing slow mag-netosonic point coincides with the separation point, theingoing Alfv´en point coincides with the inner light sur-face, and the ingoing fast magnetosonic point coincideswith the horizon.The inner and outer light surfaces are defined by theinner and outer surfaces exterior to the horizon where α = 0 ( α is defined in Equation 14 in the next sec-tion). The outer light surface is analagous to a pulsarlight cylinder, where rigid rotation would cause a parti-cle stuck to a magnetic field line to rotate at the speed oflight. The particle would also rotate at the speed of lightat the inner light surface as a consequence of gravita-tional time dilation, and it would be subluminal betweenthe light surfaces.A field line extracts black hole rotational energy whenits ingoing Alfv´en point is located inside the ergoregion.In the force-free limit this reduces to the requirementthat the magnetic field line rotate more slowly than theblack hole: 0 < Ω F < ω H . (11)Here ω H = a/ mr + is the angular velocity of a zero an-gular momentum observer on the horizon. In generalthere exist regularity conditions at each of the slow mag-netosonic, Alfv´en, and fast magnetosonic points. In theforce-free limit the most useful of these is that of the ingo-ing fast magnetosonic point, which reduces to a conditionon the fields at the horizon: √− gF θr = ± (cid:0) r + a (cid:1) (Ω F − ω H ) sin θ Σ F θφ . (12)The equations in use are insensitive to the sign of thetoroidal field, but if an ingoing observer on the horizonis to measure finite electromagnetic fields Znajek [15]showed that the negative branch must be chosen. This“Znajek condition” does not place any additional restric-tions on the fields. It was already present in our coreassumptions, most importantly the conservation of en-ergy and momentum in the vanishing divergence of thestress energy tensor. However it can be a simple checkto ensure that a solution is valid on the horizon.Further discussion of the key points along a magneticfield line may be found in Takahashi et al. [4] and Taka-hashi [16]. D. The Force-Free Transfield Equation
In our core assumptions we demanded that there beno external forcing on the magnetosphere, as stated inthe requirement that T αβ ; β = 0. If the poloidal compo-nents of this requirement are satisfied then the temporal and toroidal components are automatically satisfied. Wetherefore focus on the two poloidal components, whichmay be reduced to a single differential equation; thisequation is often referred to as the “transfield equation”as it encapsulates the force balance transverse to poloidalmagnetic field lines. In the force-free limit the transfieldequation is given by (cf. Equation 3.14 of Blandford andZnajek [1]):0 = 12 Σ∆ sin θ ddA φ (cid:0) √− gF θr (cid:1) + 1sin θ [ αF rφ ] ,r + 1∆ (cid:104) α sin θ F θφ (cid:105) ,θ − G φ sin θ (cid:18) F rφ Ω F ,r + 1∆ F θφ Ω F ,θ (cid:19) , (13)where α = g tt + 2 g tφ Ω F + g φφ Ω ,G φ = g tφ + g φφ Ω F . (14)The function α − / may be interpreted as the “gravita-tional Lorentz factor” of a particle stuck to a magneticfield line, describing the effects of both rotation and grav-itational time dilation, while G φ may be interpreted as ameasure of that particle’s rotational velocity with respectto the spacetime.Solving the transfield equation is the core difficulty offinding black hole magnetospheres. In the stationary andaxisymmetric limit of a rotating black hole, only one ex-act class of solutions is currently known (first publishedby Menon and Dermer [8] and discussed in a more generalframework by Gralla and Jacobson [9]): A φ = B h ( θ ) , Ω F = 1 a sin θ , √− gF θr = B a g ( θ ) . (15)Here B is a scalar measure of the strength of the poloidalmagnetic field; h and g are arbitrary functions of θ thatmust satisfy (prime denoting a derivative with respect to θ ): h (cid:48) = − g sin θ. (16)Solutions of this type generically contain an unphysicaldivergence in Ω F in both the θ → a → a = 0) Schwarzschildblack hole, there does exist a more reasonable variety ofexact rotating “monopolar” solutions, given by: A φ = B cos θ, Ω F = Ω , √− gF θr = B Ω sin θ. (17)Such solutions do not extract rotational energy from theblack hole (it has none to extract) but nontheless serveas a convenient starting point for exploring black holeenergy extraction. These solutions are sometimes statedas “split-monopole” solutions by changing the sign of B below the equatorial plane and appealing to an equatorialcurrent sheet. For simplicity we will always refer to themas “monopolar” solutions, as regardless of interpretationthere is still reflection symmetry across the equatorialplane. However it should be noted that “monopolar”references only the poloidal magnetic field; there is oftena significant toroidal field, as can be seen in the Ω (cid:54) = 0solutions. Setting Ω = 0 for a static monopole andapplying that solution to a slowly rotating black hole byperturbing in black hole spin (as in Blandford and Znajek[1] and more recently Pan and Yu [7]) yields: A φ = B cos θ − (cid:26) a m B R ( r ) cos θ sin θ (cid:27) , Ω F = a m + (cid:26) a m (cid:20) − R ) sin θ (cid:21)(cid:27) , √− gF θr = − B a m sin θ − (cid:26) B a sin θ m (cid:20) − R ) sin θ + 8 R ( r ) cos θ (cid:21)(cid:27) , (18)where R ( r ) = 1 m (cid:20) m + 3 mr − r
12 ln (cid:16) r m (cid:17) + 11 m m r + mr − r (cid:21) + (cid:18) r − mr m (cid:19) (cid:20) Li (cid:18) mr (cid:19) − ln (cid:18) − mr (cid:19) ln (cid:16) r m (cid:17)(cid:21) ,R = 172 (cid:0) π − (cid:1) , (19)and Li is the dilogarithm, defined as:Li ( x ) = (cid:90) x t ln (1 − t ) dt. (20)The higher order corrections by Pan and Yu [7] (andMcKinney and Gammie [6]) are shown in curly brackets.It is possible to extend the above to include O ( a ) cor-rections, but the solutions achieve a level of complication far surpassing the point of reason. There are two thingsto note from the above expressions; first, and most ob-viously, is that attempting to find perturbed solutions isnot a trivial undertaking. Perturbing one of the simplestpossible magnetic field configurations is effectively im-possible past O ( a ) corrections to A φ . The second thingto note that the leading order correction to Ω = 0 isΩ F = 0 . ω H (to first order in a ). This has partly con-tributed to the “rule of thumb” that energy extractingblack hole magnetospheres should rotate at one half therate of the black hole (Ω F ≈ . ω H also extracts the max-imum amount of black hole energy for low spin).Exploring the Ω (cid:54) = 0 solutions was one of the ini-tial questions that prompted our work, as they naturallyencompass a full range of inner Alfv´en surfaces and cor-responding magnetospheres. Unfortunately analyticallyperturbing around such solutions is even more intractablethan the already non-trivial Ω = 0 case, so we developednumerical techniques to solve for the structure of blackhole magnetospheres. Discussion of those techniques isthe topic of Section III. E. Additional Assumptions
In addition to the core assumptions of Section II A,we make two additional simplifying assumptions guidedby our goal of studying correspondences between innerAlfv´en surface location and magnetosphere structure.The first is the assumption of a monopolar geometry, inthe sense that the poloidal magnetic field line threadingthe horizon on the equator must remain on the equa-tor. This assumption is enforced via a boundary condi-tion on the flux function A φ ; we require that A φ be con-stant in the equatorial plane. In making this assumptionwe were guided by a desire to extend the Schwarzschildmonopole solutions of Equation 17 to rotating spacetimesand to avoid assuming a specific model for nearby ac-creting matter. Without the presence of such matter torestrict magnetic field lines, a configuration with vary-ing vector potential in the equatorial plane could violatethe condition that there can be no closed horizon loopsin stationary force-free magnetospheres (see MacDonaldand Thorne [17] and additional discussion in Gralla andJacobson [9]). Additonally, boundary conditions encod-ing an assumed model of matter in the equatorial planeor elsewhere have the potential to exert arbitrarily largeforcings on the magnetosphere, confounding attempts toexplore basic correlations between Alfv´en surface loca-tion and magnetosphere structure. The interested readermay see Uzdensky [18] for some complications inducedby horizon-disk connections and Fendt [19] for disk-jetconnections.The second assumption we made is that of uniformfield line angular velocity, as that is the simplest possiblemethod of classifying inner light surfaces (inner Alfv´ensurfaces in the force-free limit) and allowed us to know ex-actly where the inner light surface would be prior to cal-culating the structure of the field lines that pass throughit. While it is convenient, this assumption restricts usto a region outside the horizon that is bounded by theouter light surface. This is due to the fact that our dissi-pative numerical techniques cannot generically find mag-netospheres that pass smoothly through both inner andouter light surfaces under the assumption of uniform fieldline angular velocity, even if such solutions exist (we dis-cuss smoothness at light surfaces in Section III B). Ul-timately we do not believe this to be nearly as limitingan assumption as it might initially appear, but it is stillthe most awkward assumption that we have made. Assuch we return to it in Section V A with more detaileddiscussion. III. NUMERICAL TECHNIQUES
In order to numerically solve for a large number ofblack hole magnetospheres, we have extended the Newto-nian magnetofrictional method developed by Yang et al. [20] to the general relativistic regime. The magnetofric-tional method makes an initial guess as to the structureof a force-free magnetosphere and then gradually relaxesthe fields towards a valid force-free state. This makes ita useful tool for exploring our parameter space; once asolution for a given black hole spin and field line angularvelocity is known adjacent solutions for slightly differ-ent values are readily found. In this section we discussthe underlying theory of the magnetofrictional methodas well as the computational specifics of our implemen-tation.
A. The Magnetofrictional Method
Once initial guesses have been made for a vector poten-tial A φ , field line angular velocity Ω F , and toroidal field √− gF θr , they may be inserted into the transfield equa-tion to see if those guesses form a valid solution. This isequivalent to determining if the poloidal components ofthe divergence of the stress energy tensor vanish: T Aβ ; β = X A . (21)If the initial guesses formed a valid and self consistentforce-free solution then X A would vanish, but in general X A (cid:54) = 0 and there is an excess momentum flux that maybe interpreted as an external forcing. The core idea of themagnetofrictional method is to take that external forcingand convert it into the coordinate velocity of a fictitiousplasma via friction. In this way a physically invalid (orat least inconsistent) initial guess for a force-free mag-netosphere may be converted into a valid non force-freemagnetosphere, and the fields may be evolved towards aforce-free state by applying the equations governing in-ertial plasma flows. The first step in this process is toevaluate the excess momentum flux X A , which is found as a modification to the left hand side of the transfieldequation (Equation 13): − π Σ sin θ X A A φ,A = 12 Σ∆ sin θ ddA φ (cid:0) √− gF θr (cid:1) + . . . (22)Once the excess momentum flux has been determined, itmay be converted into the coordinate velocity of a ficti-tious plasma v A ( v r = ∂r/∂t , v θ = ∂θ/∂t ) via friction: X A = 1 ν v A . (23)Here the coefficient ν measures the strength of the fric-tion. We then assume that the fictitious plasma is aperfect conductor such that there is no electric field inits rest frame: F αβ u β = 0 . (24)Here u α is the four velocity of the plasma, which we definein terms of its coordinate velocity as v c ≡ u c /u t . Theabove may then be used to find: F tc = F cb v b . (25)If we insert this relationship into Maxwell’s homogeneousequations (Equation 3) we find the relativistic analog ofthe ideal induction equation ∂ t B = ∇× ( v × B ): F ab,t = ( F bc v c ) ,a − ( F ac v c ) ,b . (26)Inserting either F θφ or F rφ into the above then showsthat the time rate of change of the vector potential A φ is given by: A φ,t = − v A A φ,A . (27)In Appendix A we prove that application of this equationinevitably leads to a force-free configuration and in Ap-pendix B we expand it in Boyer-Lindquist coordinates.Although the coordinate velocity v A is a complicatedfunction of the metric and electromagnetic fields, at ahigh level the magnetofrictional method reduces to thestraightforward evolution of an advection equation. A va-riety of numerical techniques for evolving such equationsexist. For simplicity we have implemented an upwinddifferencing scheme similar to that described in Hawley et al. [21]. More sophisticated techniques might convergeon a valid solution more rapidly, but we found upwinddifferencing to be both very numerically stable and fastenough for our purposes.We also note that the relaxation scheme implementedby Uzdensky [18] is similar to the method describedabove, in that it uses A φ,t ∼ v A ; the primary differenceis in the weighting of v A . B. Fixing Kinks
The magnetofrictional method described above is lim-ited by the development of kinks in the magnetic field,which we now discuss. The transfield equation (Equa-tion 13) changes character at the inner and outer lightsurfaces where the gravitational Lorentz factor α changessign. This change in character is transmitted to the co-ordinate velocity v A of the fictitious plasma, which goesas: v A ∼ α ( A φ,rr + A φ,θθ ) . (28)On the inner and outer light surfaces α vanishes; it ispositive between them and negative outside. To avoidhaving our numerical scheme become anti-diffusive andunstable, outside the light surfaces we introduce an over-all multiplicative factor of − α does split the compu-tational domain into three regions, and in general thoseregions do not join smoothly. In other words for an ini-tial guess of field line angular velocity Ω F and toroidalfield √− gF θr , the solutions that the magnetofrictionalmethod finds for A φ will in general not match acrosslight surfaces. Treating A φ as a flux function, this leadsto discontinuous “kinks” in the magnetic field at the in-ner and outer light surfaces. As there are two light sur-faces, both the field line angular velocity and the toroidalfield could be evolved simultaneously in order to diminishboth kinks, as in Contopoulos et al. [10] and Nathanailand Contopoulos [11].However, as stated in Section II E, for simplicity wehave chosen to use a fixed and uniform field line angularvelocity so we do not have the freedom required to generi-cally diminish both kinks using the inherently dissipativemagnetofrictional method. We therefore divide the com-putational domain into two regions, one inside the innerlight surface and one between the light surfaces, and onlydiminish the kink at the inner light surface by evolvingthe toroidal field. The limitations and advantages of thisapproach are discussed in Section V A.In order to measure the kink at the inner light surface,the magnetofrictional method is applied separately to theregions on either side until an empirical “zero level” inthe velocity v A is reached. A vanishing plasma coordi-nate velocity corresponds to X A = 0 and a valid force-free configuration, so this procedure simply finds a “closeenough” solution in both regions. The two solutions arethen compared across the inner light surface and the mag-nitude of any kink is measured in terms of ∆ A φ . Thismeasurement is then used to modify the toroidal field ina method similar to that developed by Contopoulos et al. [22] for light cylinders in pulsar magnetospheres: ddA φ (cid:0) √− gF θr (cid:1) = ddA φ (cid:0) √− gF θr (cid:1) (cid:12)(cid:12)(cid:12)(cid:12) A φ − − γ ∆ A φ . (29)Here γ is an empirical constant, usually much less thanone, and A φ − denotes the value of A φ just inside the inner light surface. After the above equation has beenused to make a small correction to the toroidal field, themagnetofrictional method is again applied to the regionson either side of the inner light surface for a fixed numberof time steps, and then (if the set “zero level” in theplasma coordinate velocity is still achieved) the kink isre-evaluated and toroidal field re-adjusted. The cyclerepeats until matching solutions are found.Finding the correct “zero level”, value for γ , and num-ber of time steps to apply the magnetofrictional methodbetween modifications of the toroidal field is akin to criti-cally damping a simple harmonic oscillator. Poor choicescan lead to either an underdamped, wildly oscillatingkink or an overdamped kink that requires an excessiveamount of computational time to diminish. We havefound that a “zero level” of 10 − , a factor of γ = 0 . C. Measuring Force-Freeness
In principle the magnetofrictional method could be ap-plied to the limits of numerical precision. However even ifthis were practical the approximations and assumptionswe have made do not justify that level of precision. In-stead we evolve the fields until a “good enough” level isreached, which we now describe.A force-free solution is a solution with a vanishingLorentz force, F αβ J β = 0. Therefore it is reasonableto take a ratio of the Lorentz force to the magnitudeof electric current and field strength tensor in order tomeasure how force-free a solution is, an approach takenin McKinney and Gammie [6]: ξ = (cid:12)(cid:12)(cid:12)(cid:12) F µν J ν F µk J k J µ J µ F κλ F κλ (cid:12)(cid:12)(cid:12)(cid:12) . (30)If ξ (cid:28) J α J α (cid:28) ξ .In force-free magnetospheres with regions of vanishingcurrent, one alternative to ξ is to measure the ratio ofthe Lorentz force to the forces of magnetic pressure andtension, as in Malanushenko et al. [23]: ζ = | F L || F mp | + | F mt | . (31)As with ξ , if ζ (cid:28) T Aβ ; β = − F Aβ J β = A φ,A (cid:88) i =1 D i . (32)The separation of the transfield equation into D i compo-nents is done somewhat arbitrarily by grouping factorsof A φ , field line angular velocity, and the toroidal field; D ∼ √− gF θr , D ∼ A φ,rr , etc. The exact separationswe use are expanded in Appendix C. Mathematically theabove puts us in the position of having an expression ofthe form (cid:80) D i = δ and requiring that δ be “small”. As ameasure of force-freeness, we therefore take the absolutemaximum of the seven D i terms and require that δ be afactor of (cid:15) smaller than it: | δ | < (cid:15) · Max ( | D i | ) . (33)While such a method is disadvantaged by not having anobvious physical interpretation, it does not diverge whenapplied to magnetospheres with regions of vanishing cur-rent or with magnetic fields containing strong monopolarcomponents. For this work we have selected (cid:15) = 1% as alevel of sufficient force-freeness. In practice our numeri-cal techniques result in only a few segments of the innerlight surface approaching (cid:15) = 1%; over the rest of thedomain (cid:15) is significantly smaller. D. Analytic Comparisons
In Section II D we discussed some analytic solutionsto the force-free transfield equation. Here we select fourspecific solutions that are useful for making comparisonswith our numerical solutions; in Section IV D we comparethem to a numerical solution with Ω F = 0 . ω H and a =0 . m in order to examine the consequences of the (cid:15) = 1%error level discussed in the previous section.The first analytic solution that provides a useful com-parison is the exact Schwarzschild monopole describedby Equation 17. The field line angular velocity Ω of thissolution may be set to correspond to the field line angu-lar velocity of a given numerical solution. The monopoleweight B is set to B = 1 in order to match the nu-merical solution’s change in A φ between the axis and theequator (discussed in the next section). This solutionprovides a comparison with all field line angular veloci-ties but is most useful for very small spins.The second useful solution is the first order perturbedmonopole solution (Equation 18 without the terms incurly brackets). This solution is very similar to theSchwarzschild monopole with Ω = 0 . ω H , but differsslightly due to the fact that its field line angular velocityis only equivalent to 0 . ω H to first order in spin param-eter a . As with the Schwarzschild monopole solution we set B = 1. The field line angular velocity of this solutionis fixed, so it is most useful for comparison with numeri-cal solutions that have Ω F = 0 . ω H and small black holespins.The third solution is the third order perturbedmonopole solution described by Equation 18. As withthe first order solution we set B = 1. This solution ismost useful for comparison with numerical solutions ofΩ F = 0 . ω H for slightly higher spins than the first ordersolutions, and for determining how a higher order per-turbation in spin translates to changes in our measure oferror (cid:15) . It should be noted that the field line angular ve-locity of this solution is not uniform, as in the first ordersolution, but is a function of both r and θ .The fourth solution is a fully exact solution, describedby Equation 15, with h = sin θ and B = 1. There areinfinitely many possible choices for h . Our selection ismotivated by a desire for physical plausibility, which ismost clearly seen by taking the Newtonian limit of thefields described. In that limit, our choice of h correspondsto magnetic and electric fields given by: B = − B r sin θ cos θ ˆ r − B ar cos θ ˆ φ, E = 3 B ar cos θ ˆ θ. (34)Here the vector components correspond to a standardorthonormal spherical basis. We emphasize that the fac-tor of a above is from the definition of the fields; if the a → h = sin θ was made to correspond to physically plau-sible fields, in the limited sense that they don’t divergeanywhere in our domain; less careful choices generally di-verge on the azimuthal axis. This solution is most usefulfor estimating the potential precision of a given numeri-cal grid for arbitrary black hole spin; we will refer to it asthe HS3 solution when we compare it to our numericalresults.
E. Computational Specifics
In this section we detail the computational specificsof our numerical methods, covering the grid resolutionsused, domain boundaries and boundary conditions, ini-tial conditions, and performance.The poloidal plane is divided into an ( r, θ ) grid forthe vector potential A φ . In the radial direction the gridextends from just inside the horizon to a radius r max out-side the ergosphere. The value of r max is dependent uponthe location of the outer light surface, which is in turndependent upon the choice of field line angular velocityΩ F and black hole spin. For small values of Ω F (near 0)or black hole spin the outer light surface can be a verylarge distance from the black hole, so we artificially placea cap of 20 m on the maximum radius. For large values ofΩ F (near ω H ) and black hole spin the outer light surfaceapproaches the ergosphere and r max ≈ m is used.It is critical to resolve the inner light surface, so forlarge values of Ω F where the inner light surface ap-proaches the horizon we use very small grid spacings inthe radial direction near the horizon, then gradually re-lax that spacing as r increases. For small values of Ω F the horizon and inner light surface are well separated andwe use larger radial grid spacings near the horizon, thenincrease that spacing as r extends past the ergosphere torelatively larger r max values. In general around 400 gridsquares of varied spacing in the radial direction are usedto resolve the inner light surface and ergoregion and anadditional 450 of varied spacing are used to extend toa given r max . This yields a rough average of 850 gridsquares in r per full magnetosphere; exactly how manyare used varies from magnetosphere to magnetosphere.In order to treat the value of A φ on the horizon asan evolving entity, we extend three grid squares past thehorizon towards the black hole. We cannot extend anyfurther as our implementation becomes anti-diffusive in-side the horizon and more grid squares allow numericalinstabilites to develop. We are using Boyer-Lindquist co-ordinates, which are singular on the horizon. This is nota difficulty, as we scale the magnetofrictional coefficient ν in Equation 23 by the magnitude of the poloidal mag-netic field (see Appendix B). The coordinate singularityin that scaling cancels with the coordinate singularity inthe transfield equation such that the horizon is well be-haved. As a solution is found the horizon naturally settlesinto a configuration consistent with the Znajek horizoncondition (the force-free limit of the fast magnetosonicregularity condition, Equation 12), justifying the usageof Boyer-Lindquist coordinates.In the θ direction we extend from θ = 0 on the az-imuthal axis to θ = π/ A φ = 4 and on the equa-torial plane we fix A φ = 3, for a ∆ A φ between themconsistent with a monopole of unit weight ( A φ = cos θ inEquation 17). In principle any arbitrary range betweenthe axis and equatorial plane could be used; we made ourselection for convenience.We do not use fixed boundary conditions for the vectorpotential A φ at either r min or r max . Instead, after everytime step we shoot outwards to find a boundary thatkeeps A φ smooth and use that boundary for the next it-eration. The reason for this is practical; although the axisand equatorial plane are fixed by symmetry, there is no restriction on A φ for a given r max . The inner boundarycould be fixed by applying the horizon regularity condi-tion, but we chose to have the horizon evolve as a simplecheck that the numerical algorithm is converging prop-erly.In order to calculate the derivatives in A φ requiredto calculate the coordinate velocity v A of the fictitiousplasma (Equations B4 and B5) we use centered finitedifference approximations appropriate to the local gridspacing. One-sided finite difference derivatives in A φ ap-propriate to an upwind differencing algorithm are thenused to evolve the magnetofrictional advection equation(Equation 27).Our initial conditions for each run were fairly simple.Any smooth A φ that decreased monotonically from theazimuthal axis to the equatorial plane could be used; wefound no dependence upon initial conditions. Lack ofmonotonicity caused “spikes” to develop in A φ that re-quired magnetic reconnection via numerical diffusion todissipate, greatly extending computation time when thecode remained stable. To speed convergence we found itdesirable to begin with the vector potential of the closestalready calculated magnetosphere, “close” being definedby either slightly different black hole spin or field lineangular velocity. We do not directly use the toroidalfield √− gF θr , instead we propose a function of A φ cor-responding to the derivative with respect to A φ of thesquare of the toroidal field (the left hand side of Equa-tion 29). From the transfield equation it can be shownthat this function must vanish on the axis and equatorialplane (at A φ = 4 or A φ = 3 using the boundary con-ditions described above) but is otherwise largely unre-stricted. Beyond increasing computation time we foundno dependence upon initial choice of this function. Inorder to decrease computation time we often found itdesirable to use three neighboring magnetospheres withthe same spin but different field line angular velocities to“shoot” an initial guess for this function for a given valueof A φ .The overall number of time steps required to find asolution is sensitive to the accuracy of the initial guessesfor both the vector potential and toroidal field, with themajority of computation time typically being spent re-ducing the kink at the inner light surface by finding acompatible toroidal field. In general with good initialguesses for the fields and well-tuned parameters a desk-top computer can find a magnetosphere at the (cid:15) = 1%level (Section III C) in a matter of hours; a (cid:15) ∼
10% levelcan often be achieved in less than an hour. We havefound that it is possible to significantly reduce computa-tion time by bracketing the functional form of the deriva-tive of the toroidal field using a kind of two dimensionalroot-finding algorithm with initial guesses taken from ad-jacent magnetospheres, but that technique requires care-ful setup and tuning to avoid instabilities. The na¨ıvemethod of kink reduction described in Section III B ismore robust, easier to implement, and while significantlyslower is not prohibitively so for most purposes.0
IV. RESULTS
We divide our results into four sections. First we ex-plore the general structure of the magnetospheres as afunction of black hole spin parameter and field line an-gular velocity. We then examine the rates of energy andangular momentum extraction. Lastly we explore thenumerical error of our solutions.
A. Field Line Structure
In this section we explore the behavior of poloidal mag-netic field lines. We are limited to regions near to theblack hole by our model assumptions and numerical tech-niques (Sections II E and III B). However even over thatlimited domain interesting behaviors emerged, as shownin Figures 1 and 2.Figure 1 plots poloidal magnetic field lines using spac-ing on the horizon corresponding to the strength of theradial magnetic field there (using | B r | ∼ csc θA φ,θ ap-propriate to ingoing coordinates); denser field lines im-ply greater magnetic field strength. The first and lastfield line do not lie on the axis or equatorial plane, asthose are fixed, but instead lie one grid square inwardson the horizon. The shading measures conserved momen-tum flux per unit field line ( L in Equation 10), which isequivalent to the toroidal field √− gF θr multiplied by 4 π .Figure 2 plots poloidal magnetic field lines using linearspacing in θ on the horizon. The first and last field lineagain lie one grid square inwards on the horizon. Theshading measures conserved energy flux per unit field line( E in Equation 10).The first thing to note is that as expected for very smallblack hole spin ( a ≈ . m ) all values of field line angularvelocity yield magnetospheres that are almost indistin-guishable from a Schwarzschild monopole (Equation 17,noting that “monopole” refers only to the poloidal plane).For slightly larger values of spin ( a ≈ . m ), however,noticeable deviations from a monopolar structure rapidlyemerge for both low and high field line angular velocities.As spin increases beyond a ≈ . m , only the Ω F ≈ . ω H solutions remain roughly monopolar.Low field line angular velocity solutions bend towardsthe azimuthal axis, with the strength of that bendingincreasing as black hole spin increases. For large spin( a ≈ . m ) this bending can be severe, with the ma-jority of poloidal magnetic field lines becoming nearlyparallel to the azimuthal axis at a distance r = 20 m .The source of this bending is most easily seen in Figure1; on the horizon the radial magnetic field is strong onthe azimuthal axis and weak near the equatorial plane,while the toroidal magnetic field is strong on the equato-rial plane and weak on the azimuthal axis. This discrep-ancy acts to naturally wind up the magnetic field andthe magnetosphere self-collimates into a jet-like struc-ture. The discrepancy becomes very large at high spin;for a = m (not shown) the horizon radial magnetic field strength on the axis can be ∼
10 times that on the equa-tor. For larger values of Ω F (implying larger electric fieldsas viewed by distant observers, but not necessarily otherssuch as ZAMOs) the toroidal field becomes much smaller,radial magnetic field strength becomes almost uniform onthe horizon, and the field lines begin to bend toward theequatorial plane.There are two separation surfaces shown in the fig-ures (the separation surface may be defined as the pointwhere α (cid:48) = 0, with the prime denoting differentiationalong magnetic field lines); one for a monopolar geometrywhere A φ is purely a function of θ and one correspond-ing to the numerically calculated A φ . For low field lineangular velocities and high spins the separation surfacemoves away from the horizon at higher latitudes, whilefor high field line angular velocities it moves towards thehorizon.Figure 2 makes it clear that significantly more energy isextracted along field lines near the equator for all valuesof field line angular velocity, but low Ω F and high spinmagnetospheres redirect most of that extracted energytowards the azimuthal axis. We measure this behavior inthe next section. B. Energy Extraction
In this section we explore the rate of black hole energyextraction. Energy flux is conserved along magnetic fieldlines, so the rate of black hole energy extraction is mosteasily calculated on the horizon. If we define P as thepower (energy per unit time as measured by a distantobserver) leaving the horizon, we find (cf. Lasota et al. [2]): P = (cid:90) r + T rt √− gdθdφ = 12 (cid:90) r + Ω F A φ,θ √− gF θr dθ. (35)As this is evaluated on the horizon, we may use the Zna-jek horizon condition (Equation 12) to find: P = 12 Q (1 − Q ) a r + a (cid:90) π A φ,θ sin θdθ. (36)Here Q is a unitless scaling of the field line angular ve-locity to the horizon’s angular velocity; Ω F = Qω H (forenergy extraction to take place we must have 0 < Q < P becomes: P = 5 . × · χ · r x ∗ B x G m M (cid:12) ergs , (37)where χ = 12 Q (1 − Q ) a ∗ (cid:0) r ∗ + a ∗ (cid:1) (cid:90) π A φ,θ sin θr ∗ + a ∗ cos θ dθ. (38)1 FIG. 1. The structure of poloidal magnetic field lines for various values of black hole spin and field line angular velocity.The magnitude of the outward conserved angular momentum flux per unit field line ( L in Equation 10) is colored, scaled to∆ A φ = 1 between the axis and the equatorial plane ( B = 1 in Equation 17). The inner light surface (green), ergosphere (red),horizon (black), monopole separation surface (cyan), and calculated separation surface (dotted red) are overplotted. 80% ofthe total extracted energy flows outward between the azimuthal axis and the magnetic field line drawn in solid magenta; 95%flows between the dotted magenta field line and the axis. The eight black field lines are spaced according to the magnitude ofthe radial magnetic field on the horizon; if they fall outside the 80% line they are dotted. FIG. 2. The same structure of poloidal magnetic field lines as in Figure 1, but with the magnitude of the outward conservedenergy flux per unit field line ( E in Equation 10) colored. The inner light surface (green), ergosphere (red), horizon (black),monopole separation surface (cyan), and calculated separation surface (dotted red) are overplotted. 80% of the total extractedenergy flows outward between the azimuthal axis and the magnetic field line drawn in solid magenta; 95% flows between thedotted magenta field line and the axis. The eight black field lines are spaced evenly on the horizon; if they fall outside the 80%line they are dotted. a ∗ and r + ∗ are dimensionless measures of black holespin and horizon radius; a = a ∗ m and r + = r + ∗ m . Thequantity χ is a dimensionless measure of the rate of blackhole energy extraction that varies from magnetosphere tomagnetosphere. The quantities B x and r x ∗ are measuresof monopolar magnetic field strength, in a sense that wenow explore by considering the Newtonian limit.In the Newtonian limit, a vector potential A φ = B cos θ corresponds to a monopolar magnetic field thatis given by: B = B r ˆ r. (39)Here the vectors correspond to a standard orthonormalspherical basis. The weighting of B on A φ may thereforebe interpreted as a measure of the monopolar magneticfield strength at a specified radius. For convenience ournumerical solutions for A φ assume unit spacing betweenthe azimuthal axis and the plane ( B = 1), but anyarbitrary weighting would yield identical results. Thequantities B x and r x ∗ are measures of that arbitraryweighting; B x is the field strength in Gauss of a New-tonian monopole at dimensionless radius r x ∗ . For ex-ample, B x = 100G and r x ∗ = 20 would correspond toa monopolar magnetic field strength of 100 Gauss at adistance r = 20 m from a black hole of mass m . Somecaution should be taken in applying this interpretation,however. None of our numerical solutions are exactlymonopolar, and the strength of the magnetic field is anobserver dependent quantity that does not have a simpletranslation from a flat space value to the spacetime of arotating black hole. The above interpretation of B x and r x ∗ is made purely for simplicity and convenience, and ingeneral should be taken to be nothing more than a roughestimate.In Figure 3 we plot the value of χ as a measure ofthe rate of black hole energy extraction for select valuesof black hole spin and field line angular velocity. Wesuppress spins greater than a = 0 . m as showing themwould compress the curves of smaller spin that are ofgreater interest.For a given value of Ω F , we note that increasing blackhole spin always increases the rate of energy extractionbut that changes in field line angular velocity can have amuch larger effect than changing spin. The peak energyextracted for a given spin occurs at Ω F ≈ . ω H for lowspin and approaches 0 . ω H at high spin. The maximumenergy extracted spans two orders of magnitude; χ Max =0 .
03 for a ≈ . m , χ Max = 0 .
003 for a ≈ . m , and χ Max = 0 . a ≈ . m .To a very good approximation, the functional form of χ ( a ) for fixed field line angular velocities Ω F = 0 . ω H ,Ω F = 0 . ω H , and value for maximum energy extractionΩ F Max are given by: χ . ω H ≈ . × − (cid:0) . a . ∗ + 1 . a ∗ (cid:1) ,χ . ω H ≈ . × − (cid:0) . a . ∗ + 1 . a ∗ (cid:1) ,χ Max ≈ . × − (cid:0) . a . ∗ + 1 . a ∗ (cid:1) . (40) Ω F / ω H . . . . . . χ a * = 0.99 a * = 0.95 a * = 0.875 a * = 0.75 a * = 0.5 FIG. 3. The rate of black hole energy extraction for selectmagnetospheres; χ is defined in Equation 38. For each value ofblack hole spin 20 magnetospheres were calculated, indicatedfor a ∗ = 0 .
875 with blue dots. The approximate maximumvalues of χ for various spin values are indicated with red stars;they corresponding to the 23 spin parameters a ∗ χ occurs near Ω F = 0 . ω H ;for a ∗ = 0 .
99 the maximum occurs near Ω F = 0 . ω H . For ex-tremal spin ( a ∗ = 1, not shown) the rate of energy extractionpeaks slightly above χ = 0 . In this form it is clear that for a given spin the maxi-mum rate of energy extraction corresponding to Ω
F Max isroughly three times that of the lower values Ω F ≈ . ω H or Ω F ≈ . ω H , an effect that can be much larger thanchanges in spin.For low field line angular velocities the magnetic fieldbends towards the azimuthal axis, as shown in Figures 1and 2. So while most energy is extracted on the horizonnear the equatorial plane, a short distance from the hori-zon a large fraction of it ends up flowing outward alongthe azimuthal axis. In order to explore this behavior, wecalculated the energy escaping through a spherical cone(both upper and lower hemispheres) at a radius r = 20 m for Ω F = 0 . ω H . The results are plotted in Figure 4.For high spins ( a ∗ ≈ .
85 and greater), over 95% of theextracted energy escapes through a 45 ◦ cone, and 80%escapes through a cone of less than 30 ◦ . This indicatesthat for high spin and low field line angular velocities ablack hole magnetosphere is naturally inclined to extractenergy via jet-like structures aligned with the rotationalaxis of the black hole. C. Angular Momentum Extraction
In this section we explore the rate of black hole angu-lar momentum extraction. As with energy flux, angularmomentum flux is conserved along magnetic field lines,4
10 20 30 40 50 60 70 80 90
Degrees from Azimuthal Axis (at r = 20 m ) . . . . χ a * = 0.99 a * = 0.95 a * = 0.875 a * = 0.75 a * = 0.590%80%50% 95% FIG. 4. The rate of black hole energy extracted as a functionof the angle of a spherical cone at a radius of 20 m (in bothupper and lower hemispheres) for a field line angular velocityΩ F = 0 . ω H . The angles through which 50%, 80%, 90%, and95% of the total extracted energy escapes are overplotted. and is most easily calculated on the horizon by exploitingthe Znajek horizon condition. If we define K as the rateof angular momentum extraction, we find: K = − (cid:90) r + T rφ √− gdθdφ = 12 (cid:90) r + A φ,θ √− gF θr dθ. (41)Note that this only differs from the rate of energy ex-traction (Equation 35) by a factor of field line angularvelocity Ω F , consistent with the conserved energy E andangular momentum L per unit flux tube differing by thesame factor (Equation 10). In CGS units K reduces to: K = 2 . × · ϕ · r x ∗ B x G m M (cid:12) erg , (42)where ϕ = 12 (1 − Q ) a ∗ (cid:90) π A φ,θ sin θr ∗ + a ∗ cos θ dθ. (43)The primary structural difference between these expres-sions and those for the rate of energy extraction (Equa-tions 37 and 38) is the fact that ϕ does not vanish for Q = 0. This is a statement that it is possible to spindown a black hole without extracting energy. In Figure5 we plot the value of ϕ as a measure of the rate of blackhole angular momentum extraction for the same values ofblack hole spin and field line angular velocity that wereused to plot the rate of energy extraction.To a very good approximation, the functional form of ϕ ( a ) for fixed field line angular velocities Ω F = 0 . ω H ,Ω F = 0 . ω H , and value for maximum energy extraction Ω F / ω H . . . . ϕ a * = . a * = . a * = . a * = . a * = . FIG. 5. The rate of black hole angular momentum extrac-tion for the selection of magnetospheres used in Figure 3; ϕ is defined in Equation 43. For each spin the angular momen-tum extracted at the maximum rate of energy extraction isindicated by red stars. The blue dots indicate the 20 valuesof field line angular velocity used for the a ∗ = 0 .
875 case.The red stars indicate the 23 different values of spin listed inFigure 3. Ω F Max are given by: ϕ . ω H ≈ . × − (cid:0) . a . ∗ + 2 . a . ∗ (cid:1) ,ϕ . ω H ≈ . × − (cid:0) . a . ∗ + 9 . a . ∗ (cid:1) ,ϕ Max ≈ . × − (cid:0) . a . ∗ + 2 . a . ∗ (cid:1) . (44)In this form it is clear that the Ω F = 0 . ω H solutionsextract nearly twice as much angular momentum as thesolutions that maximize the rate of energy extraction.It is also clear that the Ω F = 0 . ω H solutions extractaround 10% of the angular momentum of the Ω F = 0 . ω H solutions, while from Equation 40 we note that these so-lutions extract roughly the same amount of energy. Thisindicates that high Ω F solutions can extract energy froma black hole for a much longer period of time than lowerΩ F solutions can, as it will take longer for high Ω F mag-netospheres to spin down the black hole. D. Numerical Error Estimation
In this section we explore the numerical error of oursolutions. In Section III D we discussed four analyticsolutions that can provide useful comparisons with ournumerical results, most especially in determining the re-liability of the (cid:15) = 1% error level we have set (discussedin Section III C and Appendix C).The three perturbed monopole and exact
HS3 solu-tions have varying regions of applicability. We examinethe numerical solution for a = 0 . m and Ω F = 0 . ω H so that we can reasonably compare all four solutions si-multaneously. For convenience we begin by comparing5them along a slice of constant θ = 45 ◦ from just insidethe horizon to r max = 20 m (Figure 6). Along that sliceit is apparent that the Schwarzschild monopole and firstorder perturbed monopole are very similar, with bothtypically having errors of less than 5%. The third orderperturbed solution does better in general, with an errorof around 1% near the inner light surface and inside theergosphere. Aside from the two grid squares immediatelyadjacent to the inner light surface the numerical solutionhas an error of less than 0.05% across the entire slice,and is sometimes “better” than the HS3 solution’s errorof around 0.0005%. We place quotes around “better” be-cause the error of the
HS3 solution is a rough measure ofthe precision of our numerical grid and derivatives; sig-nificantly exceeding ∼ . θ values as well as differ-ent spin parameters and field line angular velocities yieldqualitatively similar results. Radius (m) -4 -3 -2 -1 ε ( % ) st Order MonopoleSchwarzschild MonopoleNumericalHS33 rd Order Monopole
FIG. 6. The percent error along a θ = 45 ◦ slice for a = 0 . m and Ω F = 0 . ω H ; the radial spacing corresponds to uniformspacing in the numerical grid. The horizon (H), inner lightsurface (ILS), ergosphere (E), and separation point (SP) areindicated by vertical lines. With the exception of the two gridsquares immediately adjacent to the inner light surface thenumerical solution is significantly better than the perturbedsolutions. The various solutions used are listed in SectionIII D In order to examine the increase in error along the in-ner light surface we plot the percent error of the numeri-cal solution in the two grid squares immediately adjacentto the inner light surface in Figure 7. This is done for allvalues of θ , again for a = 0 . m and Ω F = 0 . ω H . Eventhough substantial portions of the inner light surface arewell below the 1% level, there are spikes up to 1%. Re-ducing those spikes is what takes the largest amount ofcomputational time. The perturbed solutions are eithercomparable to or greatly exceed the error in the numeri-cal solutions along the entirety of the inner light surface.The HS3 solution is typically 10-100 times better thanany other solution. The exception is near the azimuthal axis where the
HS3 solution’s field line angular velocitydiverges. Note that the error of the
HS3 solution doesnot also diverge, as the combined equations are well be-haved there; the error increase comes from amplificationof the error in taking numerical derivatives. The impli-cations of the varying error levels are discussed in moredepth in Section V C.
10 20 30 40 50 60 70 80 90
Degrees from Azimuthal Axis -3 -2 -1 ε ( % ) Schwarzschild Monopole1 st Order Monopole 3 rd Order MonopoleNumerical (ILS +/-)HS3
FIG. 7. The percent error along the inner light surface for a = 0 . m and Ω F = 0 . ω H . The grid squares on either side ofthe inner light surface (red and black) have almost identicalerror; reducing the spikes in that error to the 1% threshold wehave set is what takes the most computational time. Despitethose spikes the numerical solution still generally does as goodor better than the perturbed solutions. The upticks in thethird order monopole and exact HS3 solutions near θ = 0 arefrom the conversion of the toroidal field to a differentiatedfunction of A φ . Most of the terms in the transfield equationare very small there, so numerical errors in that conversionare amplified. In order to reduce computation time we make a guessas to the ultimate structure of the fields to use as aninitial condition. It is therefore reasonable to ask ifour guesses bias our results towards solutions that areclose to that initial guess and ignore other potentiallyvalid magnetospheres. There are two unknown quanti-ties that we make guesses for; the toroidal component ofthe vector potential, A φ ( r, θ ), and the derivative of thesquare of the toroidal field with respect to the vector po-tential as an unknown function of the vector potential, d/dA φ ( √− gF θr ) = F ( A φ ). In order to examine howsensitive our solutions might be to changes in the initialform of these functions, we examine the magnetospheresobtained for black hole spin parameter a = 0 . m ; thenumerically obtained functions of the derivative of thetoroidal field for those magnetospheres are shown in Fig-ure 8.To asses how initial conditions might modify our re-sults, we first coupled the numerical solution for thederivative of the toroidal field for the Ω F = 0 . ω H mag-netosphere with the vector potential A φ obtained for ev-6 - . - . - . - . D e r i v a ti v e o f T o r o i d a l F i e l d ω H ω H ω H ω H Initial Condition B Initial Condition A( θ = π /2) ( θ = 0) φ -8 -6 -4 -2 R e l . D i ff e r e n ce FIG. 8. The derivatives of the toroidal field d/dA φ ( √− gF θr ) obtained numerically for black hole spin parameter a = 0 . m and various values of field line angular velocity Ω F (in incre-ments of 0 . ω H ). The blue and red dashed lines indicate twodifferent initial guesses for the derivative of the toroidal fieldin the Ω F = 0 . ω H case. The blue and red dotted line is thefinal result in both cases; they completely overlap and obscurethe original result obtained using a much better initial guess.The relative difference between the results obtained using theA and B initial conditions and the original result are shownin the bottom panel. ery other value of field line angular velocity (extrapolat-ing outward to r = 20 m for the higher field line angularvelocity solutions). In every case we found that the nu-merical code rapidly converged to a solution essentiallyindistinguishable from the original Ω F = 0 . ω H solu-tion. There were some minor deviations at large radii, asthe 1% error level was achieved there last (in the originalsolution it was achieved last on the inner light surface),but significantly less than the deviation between adjacentΩ F = 0 . ω H and Ω F = 0 . ω H magnetospheres.As modifying the initial vector potential seemed tohave no effect, we next examined modifying our initialguess for the derivative of the toroidal field. Poor guessesfor this function can result in differences (“kinks”) in A φ across the inner light surface that can easily exceedthe difference in A φ between the azimuthal axis and theequatorial plane. Diminishing inner light surface kinkstakes the most computation time, so good guesses forthe derivative of the toroidal field are the most criticalinitial condition and most likely source of any sensitivityour numerical solutions might have to initial conditions.For clarity we focus on two alternative initial guessesfor the functional form of the derivative of the toroidalfield, shown in Figure 8 as initial conditions A and B.They are symmetric quadratics in A φ that straddle thenumerically obtained function, with the addition of peri-odic oscillations in Case A in an attempt to make a verypoor initial guess without being overly ridiculous. Bothcases led to very large initial kinks across the inner light surface and significantly increased the computation timerequired to find a solution. However the solutions ob-tained in both cases were again essentially indistinguish-able from the original Ω F = 0 . ω H solution; the minordeviations in the derivative of the toroidal field (shownin in the bottom panel of Figure 8) are dwarfed by thedeviation between adjacent Ω F = 0 . ω H and Ω F = 0 . ω H magnetospheres. V. DISCUSSION
The most limiting assumption underlying our solutionsis that of uniform field line angular velocity. In this sec-tion we discuss that assumption in more detail, considertwo extremes of the magnetospheres that we found, andbriefly explore the numerical error of our solutions.
A. Assumption of Uniform Ω F As stated in Section II E, we have assumed uniformfield line angular velocities Ω F in order to simplify thetask of studying how the location of the inner Alfv´ensurface might correspond to changes in the structure ofenergy extracting black hole magnetospheres. This re-stricted us to solving for the structure of magnetospheresinside the outer light surface. With solutions in handwe return and examine that restriction in more detail.We begin by examining how limiting the assumption ofuniform Ω F might be in the two extreme classes of mag-netospheres (Ω F /ω H ≈ F ≈ ω H ) that we found.We then discuss the existence and uniqueness of solutionsthat pass smoothly through both light surfaces. We closewith a discussion of how limiting uniform Ω F might bein interpreting our results.For low field line angular velocities (Ω F /ω H ≈
0) wherethe magnetic field lines bend towards the azimuthal axisthe limitations imposed by solutions restricted to the in-terior of the outer light surface should not be a signif-icant concern. In that case the outer light surface canbe many thousands of gravitational radii away from thehorizon near the azimuthal axis and formally infinitelyfar away exactly on the axis. Diminishing the kink onsuch a distant surface could slightly modify the struc-ture of the magnetosphere near the horizon, but in rea-sonable application we would generally expect deviationsfrom our core assumptions of stationarity, axisymmetry,a perfectly conducting force-free plasma over such an ex-tended region to be far more significant. Should that notbe the case, we also see a smooth transition from ourΩ F = 0 . ω H solutions to our Ω F = 0 solutions for whichan outer light surface does not exist (formally located in-finitely far away from the black hole). We see no reasonto expect pathologies in the limit Ω F →
0, implying thatour low Ω F solutions are close to ones those that passsmoothly through both light surfaces.7For high field line angular velocities (Ω F ≈ ω H ) wherethe magnetic field lines bend towards the equatorial planewe would expect to find a connection to nearby accret-ing matter (not necessarily a thin disk limited to theequatorial plane). Our boundary condition of a singlemagnetic field line in the equatorial plane might be rea-sonable close to the horizon, but moving away from thehorizon it would become increasingly likely to find sig-nificant deviations depending upon the specific model ofnearby accreting matter that was chosen. Finding solu-tions that passed smoothly through both light surfaceswould therefore be somewhat ill-advised, as it would in-volve modifying the near horizon magnetosphere usingrestrictions found in unreasonably distant regions. Amore realistic magnetosphere would include a descrip-tion of nearby accreting matter that included a modelof plasma injection into an inflow interior to or near theseparation surface. Our solutions could be representativeof such inflows, especially a thick disk or torus near theergosphere. A solution that passed smoothly through anouter light surface, on the other hand, would also need toinclude some model of plasma injection into an outflowconsistent with the field line geometry obtained, whichmight be difficult.In arguing that it is possible to usefully interpret theabove limits on Ω F despite the limitations of our do-main, we have ignored the question of whether or notsimilar solutions that pass smoothly though both lightsurfaces even exist. Using essentially identical bound-ary conditions on the vector potential A φ , Contopoulos et al. [10] and Nathanail and Contopoulos [11] (hereafterCKP13 and NC14) modified both the field line angu-lar velocity and toroidal field to find solutions that passsmoothly through both light surfaces. In doing so theyfound a single nearly monopolar solution with varyingΩ F ≈ . ω H that is very similar in structure to our uni-form Ω F = 0 . ω H solutions. It is therefore natural to askif the solutions found by CKP13 and NC14 are unique.The exact solutions of Equation 15, while non energy ex-tracting and largely unphysical, indicate that they arenot; there are in fact infinitely many solutions fulfill-ing our boundary conditions on A φ that pass smoothlythrough both an inner and outer light surface. The ques-tion then becomes why we found the solutions that wedid regardless of initial conditions, and why CKP13 andNC14 similarly found a single solution.The answer to that question lies in our numerical tech-niques. We show in Appendix A that the magnetofric-tional method works by minimizing the energy in theelectromagnetic fields. In minimizing the kink across theinner light surface we are finding matched minimum en-ergy solutions. The apparent uniqueness of our results isa suggestion of the uniqueness of a minimum energy solu-tion, not the uniqueness of our solution in full generality.The majority of numerical codes are common in the be-havior of being dissipative and seeking minimum energystates; any code that allows energy to be added at willis likely to be numerically unstable. For example there is a full range of valid rotating Schwarzschild monopolesolutions (Equation 17), but it would be a very uniquenumerical code that would be capable of blindly findingall of them. Most codes would always converge on theminimum energy Ω = 0 solution.For any value of black hole spin, the minimum energymagnetosphere consistent with our boundary conditionswill be as close to monopolar as possible (Ω F ≈ . ω H ).The bunching of magnetic field lines towards the axis orequator seen in our low and high Ω F solutions requiresthe addition of energy to spin the magnetosphere awayfrom that monopole. It is therefore not surprising thatCKP13 and NC14 found a single roughly monopolar so-lution that passes smoothly through both inner and outerlight surfaces; without the explicit and very careful addi-tion of energy to move to and maintain a rotating mag-netosphere any stable code would likely converge on thatsolution. Our utilization of uniform Ω F may be inter-preted as a way of exploring the structure of arbitrar-ily rotating magnetospheres using a dissipative numericalcode, and we believe that a large fraction of our solutionscould be taken to be good approximations of solutionsthat pass smoothly through both light surfaces. In real-istic application, however, they should largely be takento be representative of inflow solutions. Outflow solu-tions would require the additional description of plasmainjection mechanisms, and even if our solutions passedsmoothly through an outer light surface there is no guar-antee that plasma parameters would be either contin-uous or conserved across an extended plasma injectionregion (as implied by single solutions that pass smoothlythrough both light surfaces).We have argued that our limited domain might not pre-vent useful interpretations and that our solutions mightapproximate solutions that pass smoothly through bothlight surfaces (even if such solutions might be of dubi-ous value). It could still be asked how representative oursolutions might be of energy extracting black hole magne-tospheres, as completely uniform Ω F magnetospheres areunlikely to exist. We cannot fully answer that question,but we emphasize that exploring uniform Ω F magneto-spheres is not the goal of this work. The goal of this workis to explore the effects of inner Alfv´en surface location.Uniform Ω F in the force-free limit is merely a useful toolto do so, both in its obvious simplicity as well as in allow-ing us to know exactly where the inner light surface willbe prior to solving for the structure of the magnetospherethat passes through it. Slightly deforming the locationof the Alfv´en surfaces used here will not have significanteffects, and in general we expect the Alfv´en surface tobe continuous, so the solutions obtained here could bemerged to explore more complex scenarios. For example,if one desired an ingoing solution with a connection toaccreting matter near the equator and the launching ofa jet-like structure from the horizon along the rotationalaxis, an appeal could be made to an Alfv´en surface thatlies closer to the boundary of the ergosphere at higherlatitudes and closer to the horizon near the equator. On8the horizon this would correspond to lower Ω F for smallvalues of θ and higher Ω F for large values of θ , compat-ible with the simulations conducted by McKinney andGammie [6]. B. Two Types of Magnetospheres
Our solutions indicate two extreme methods for blackhole energy to be extracted and transmitted to distantobservers. For low field line angular velocities (Ω F /ω H ≈
0) jet-like structures naturally form to transmit extractedenergy directly from the horizon to distant observersalong the azimuthal axis. For high field line angular ve-locities (Ω F ≈ ω H ) we find magnetospheres that are con-sistent with the transmission of energy from the horizonto nearby accreting matter, which in turn might transmitthe extracted energy to distant observers. In this sectionwe discuss potential consequences of those two types ofmagnetospheres, noting at that outset that any poten-tial distinction between the two in realistic astrophysicalscenarios is likely to be fuzzy.If we examine the rate of energy and angular momen-tum extraction through Equations 40 and 44, we seethat the two types might have significant time-dependentdistinctions due to their differing rates of angular mo-mentum extraction and concurrent black hole spindown.Therefore transient high energy phenomena powered byblack hole energy extraction might also have two distinctsignatures. For specificity, consider a gamma-ray burst(GRB). Black hole energy extraction is a candidate forthe central engine of a GRB (e.g. Gehrels and M´esz´aros[24]), and the exponential decay associated with blackhole energy extraction might be compatible with someGRB observations (e.g. Nathanail et al. [25]). We do notdemand a specific type or model of GRB, however. Weonly require a transient high energy event that might bepowered by black hole energy extraction, and “GRB” isa convenient specific stand-in for“transient high energyevent” even if our treatment is too crude to completelydescribe or differentiate between realistic GRBs of anyspecific type.If we take B x = 10 G and r x = 3 as constants inEquations 37 and 42 (i.e. a magnetic field strength ofroughly 10 G just outside the ergosphere), assume that a = 0 . m and m = 4 M (cid:12) at time t = 0, and for simplic-ity insist upon a state of suspended accretion (such as amagnetically arrested disk in a low accretion state), wefind the rates of black hole energy extraction plotted inFigure 9.We immediately note that all three cases exhibitroughly exponential decay, but the Ω F = 0 . ω H case de-cays relatively slowly. After roughly 20-30 seconds boththe Ω F = 0 . ω H and Ω F Max luminosities have decreasedby factors of ∼ due to the black hole spin droppingbelow a ≈ . m . The rate of angular momentum ex-traction for the 0 . ω H case is far lower, however, andafter 30 seconds its luminosity has not significantly de- Time (s) L u m i no s it y ( e r g / s ) ω H ω H Ω F Max
FIG. 9. Rate of black hole energy extraction for constantfield line angular velocities Ω F = 0 . ω H , Ω F = 0 . ω H , andthe field line angular velocity corresponding to the maximumrate of energy extraction for a black hole with initial spin a = 0 . m . The curves for Ω F = 0 . ω H and Ω F = Ω F Max areterminated when the black hole’s spin falls below a = 0 . m ;it takes around 520 seconds for the Ω F = 0 . ω H case to reachthat level. creased. Note that in this model larger black hole masseswould take less time to spin down due to the increase inmagnetic field strength implied by r x corresponding to alarger dimensional radius. For example, a black hole withtwice the mass but otherwise identical parameters wouldbe spun down in around half the time. Also note thatthe scaling on the axes is somewhat irrelevant; they werechosen for crude correspondence to a GRB, but a widerange of systems with different parameters, luminosities,and timescales exhibit qualitatively identical curves.It takes over 350 seconds for the 0 . ω H luminosity todrop to the same level that the 0 . ω H luminosity fell toin around 20 seconds. We therefore expect that GRBswith magnetic field lines that directly connect the horizonto distant observers might be at least 10-20 times shorterthan GRBs with magnetic field lines connecting to a diskor other matter before connecting to distant observers.The addition of more realistic effects such as accretionactively spinning up the black hole would obviously com-plicate matters, and as noted at the outset the bound-ary between both types of magnetospheres is unlikely tobe very distinct. However, if GRBs (again emphasizingthe limited sense in which we are using the term) arepowered by black hole energy extraction it would not beunreasonable to expect two broad classes with differingcharacteristic durations as well as different radiation sig-natures; the longer lived 0 . ω H type might also exhibita greater degree of thermalization due to the depositionof the extracted energy into inflowing matter before thatenergy is transmitted to distant observers.The analysis presented here is too crude to reason-ably claim that distinctions between GRBs can arise due9to differences in the location of their inner Alfv´en sur-face, encapsulated here in differences in field line an-gular velocity. A more reasonable model is beyond thescope of our present work, but the possibility that sucha model could allow finer classification of some GRBs orother transient high energy phenomena via correlation oftimescale and degree of thermalization is interesting. C. Numerical Error Analysis
In this section we discuss how reliable the numericalcalculations underlying our solutions might be.With the possible exception of the kink at the innerlight surface, we find that our numerical solutions areconsistent with the precision of their numerical grid. Forthe case of a = 0 . m and Ω F = 0 . ω H , the error of oursolution is generally around 100 times smaller than theerror in any of the perturbed monopole solutions thatare considered to be good approximations at low spin.While the numerical error involved in implementing theexact HS3 solution (arising mostly from taking numericalderivatives) can be another 100 times smaller still, thereare large regions in which the numerical solution is ac-tually “better” than the
HS3 solution and likely exceedsthe precision justified by our numerical grid.Along the inner light surface the error in our numericalsolution is generally a few times smaller than the third-order perturbed monopole solution, and exceeds the first-order solution by a factor of 10-100. Again the
HS3 solu-tion generally exceeds the numerical solution by anotherfactor of 10-100. However significantly better solutionswould be difficult to obtain. At the 1% level the “kink”(change in A φ ) in our solution across the inner light sur-face in the radial direction is around 10 − , much smallerthan the required linear change in A φ of 5 × − in the θ direction implied by our usage of 200 grid squares be-tween the azimuthal axis (where A φ = 4) and the equator(where A φ = 3). This means that it becomes very diffi-cult to accurately quantify the kink, much less diminishit further.Some experimentation reveals that a much higher er-ror level of around 10% would not have changed ourresults in any significant way. This is consistent withthe error of the perturbed solutions, which should beaccurate enough for most purposes and have similarlysized errors. Reducing the error of our solutions to the1% level increased the required computation time signifi-cantly, and achieves a precision that almost certainly farexceeds what our assumptions allow for in any reason-able application. However reducing the error to the 1%level allows us to confidently state that any deficienciesin our results rest within our assumptions and not in ournumerical calculations.Due to the importance of selecting appropriate initialconditions in diminishing computation time, it might beasked if our initial guesses in any way affect or driveour results. As shown in Figure 8, even with somewhat ridiculous choices of initial condition we ultimately ar-rive at effectively identical solutions that are appreciablydifferent from neighboring solutions. We cannot (and donot) claim that this implies that our solutions are uniquein general; it is possible that there are many solutionsthat satisfy our boundary conditions and assumptions.The only potentially persuasive uniqueness in our solu-tions lies in unique minimum energy solutions that arematched across the inner light surface, as a consequenceof the magnetofrictional method being an energy mini-mizing algorithm (shown in Appendix A). In developingour numerical code we conducted many more trials andexperiments than are shown here, and we never observedany indication that it might be possible to converge ontwo different solutions. Therefore while we cannot provethat our solutions are unique matched minimum energysolutions consistent with our boundary conditions andassumptions, we do believe that it is a reasonable as-sumption to make. VI. CONCLUSIONS
The defining feature of energy extracting black holemagnetospheres is the location of their inner Alfv´en sur-face, which coincides with the inner light surface in theforce-free limit. Despite its importance, no comprehen-sive studies of the effects of modifying the location of theinner Alfv´en surface have been accomplished. We havebegun addressing that deficiency by studying how simpleforce-free magnetospheres can be modified as a functionof uniform field line angular velocity, which together withblack hole spin determines the location of the inner lightsurface. In order to do so we extended the Newtonianmagnetofrictional method for computing force-free mag-netospheres into the general relativistic regime, whichallowed us to efficiently calculate hundreds of energy ex-tracting black hole magnetospheres as functions of innerAlfv´en surface location. In so doing we found that innerAlfv´en surfaces near the horizon cause extracted energyto flow towards the equatorial plane, while inner Alfv´ensurfaces near the boundary of the ergosphere cause ex-tracted energy to flow outwards via jet-like structuresaligned with the azimuthal axis. Applied to transienthigh energy phenomena, those magnetospheres imply twotimescales that might differ by a factor of 20 or more.This suggests that two classes of transient high energyphenomena might exist; shorter ones that directly con-nect the horizon to distant observers, and longer onesthat have a disk or other significant matter separatingthe energy flow from the horizon to distant observers thatpotentially thermalizes the extracted energy, creating adifferent radiation signature.0
ACKNOWLEDGMENTS
K.T. gratefully acknowledges useful conversations withDana Longcope. M.T. was supported by JSPS KAK-ENHI Grant Number 24540268. [1] R. D. Blandford and R. L. Znajek, Mon. Not. R. Astron.Soc. , 433 (1977).[2] J.-P. Lasota, E. Gourgoulhon, M. Abramowicz,A. Tchekhovskoy, and R. Narayan, Phys. Rev. D ,024041 (2014), arXiv:1310.7499 [gr-qc].[3] S. S. Komissarov, Mon. Not. R. Astron. Soc. , 427(2004).[4] M. Takahashi, S. Nitta, Y. Tatematsu, and A. Tomi-matsu, Astrophys. J. , 206 (1990).[5] S.-Y. Nitta, M. Takahashi, and A. Tomimatsu, Phys.Rev. D , 2295 (1991).[6] J. C. McKinney and C. F. Gammie, Astrophys. J. ,977 (2004), astro-ph/0404512.[7] Z. Pan and C. Yu, Astrophys. J. , 57 (2015),arXiv:1504.04864 [astro-ph.HE].[8] G. Menon and C. D. Dermer, General Relativity andGravitation , 785 (2007), astro-ph/0511661.[9] S. E. Gralla and T. Jacobson, Mon. Not. R. Astron. Soc. , 2500 (2014), arXiv:1401.6159 [astro-ph.HE].[10] I. Contopoulos, D. Kazanas, and D. B. Papadopoulos,Astrophys. J. , 113 (2013), arXiv:1212.0320 [astro-ph.HE].[11] A. Nathanail and I. Contopoulos, Astrophys. J. , 186(2014), arXiv:1404.0549 [astro-ph.HE].[12] K. S. Thorne and D. MacDonald, Mon. Not. R. Astron.Soc. , 339 (1982).[13] J. D. Bekenstein and E. Oron, Phys. Rev. D , 1809(1978).[14] M. Camenzind, Astron. Astrophys. , 137 (1986).[15] R. L. Znajek, Mon. Not. R. Astron. Soc. , 457 (1977).[16] M. Takahashi, Astrophys. J. , 264 (2002), astro-ph/0201327.[17] D. MacDonald and K. S. Thorne, Mon. Not. R. Astron.Soc. , 345 (1982).[18] D. A. Uzdensky, Astrophys. J. , 889 (2005), astro-ph/0410715.[19] C. Fendt, Astron. Astrophys. , 1025 (1997).[20] W. H. Yang, P. A. Sturrock, and S. K. Antiochos, As-trophys. J. , 383 (1986).[21] J. F. Hawley, L. L. Smarr, and J. R. Wilson, Astrophys.J. Suppl. Ser. , 211 (1984).[22] I. Contopoulos, D. Kazanas, and C. Fendt, Astrophys.J. , 351 (1999), astro-ph/9903049.[23] A. Malanushenko, C. J. Schrijver, M. L. DeRosa,and M. S. Wheatland, Astrophys. J. , 102 (2014),arXiv:1312.5389 [astro-ph.SR].[24] N. Gehrels and P. M´esz´aros, Science , 932 (2012),arXiv:1208.6522 [astro-ph.HE].[25] A. Nathanail, A. Strantzalis, and I. Contopoulos, Mon.Not. R. Astron. Soc. , 4479 (2016), arXiv:1507.02143[astro-ph.HE]. Appendix A: Proof of Convergence
In this section we show that the magnetofrictionalmethod inevitably drives a given vector potential A φ to-wards a force-free state. We do this in two steps; first,we show that a force-free state is a minimum (extremum)in the electromagnetic energy of a given volume. Sec-ond, we show that the magnetofrictional method reducesthe electromagnetic energy contained in a given volume,thereby inevitably driving that volume towards a mini-mum energy and force-free state. For simplicity we as-sume Boyer-Lindquist coordinates, that a force-free stateexists, and that no pathologies develop (e.g. magneticreconnection, if necessary, is assumed to occur via nu-merical diffusion).We define the energy W contained within the fields ata given coordinate time to be: W ≡ (cid:90) T tt √ γd x. (A1)Here √ γd x is the appropriate proper volume elementand T αβ is the electromagnetic stress energy tensor. Dif-ferentiating with respect to coordinate time, we find: W ,t = (cid:90) T tt ; t √ γd x = (cid:90) ( T tα ; α − T ta ; a ) √ γd x = − (cid:90) F tβ J β √ γd x − (cid:73) T tA d Σ A . (A2)Here d Σ A is the appropriate directed surface elementbounding the region such that the second term is a mea-sure of the net Poynting flux through the region’s bound-ary. We have assumed no field pathologies, so any fieldline entering the volume also exits the volume. Energyflux per unit field line is conserved, so there can be nonet Poynting flux into the region along any magnetic fieldlines and the second term necessarily vanishes. In a force-free configuration F αβ J β = 0 and the first term also van-ishes; therefore a force-free configuration extremises thefield energy W in coordinate time.The core conceit of the magnetofrictional method isthat the excess momentum flux of a non force-free con-figuration may be converted into the coordinate velocityof a fictitious plasma; mathematically, this may be statedas: − F aβ J β = 1 ν v a . (A3)1Here the coordinate velocity v a of the plasma is defined interms of its four velocity u a as v a ≡ u a /u t . Both the fieldline angular velocity Ω F and the toroidal field √− gF θr are assumed to be functions of A φ (i.e. conserved alongmagnetic field lines), so v φ = 0. We now make the sim-plifying assumption that Ω F vanishes (i.e. F tr = F tθ = 0and there is no electric field from the perspective of adistant observer); we will address the Ω F (cid:54) = 0 case later.We then add in a fictitious electric field that allows v A to fulfill the condition of a perfectly conducting plasma, F αβ u β = 0: F tr = − F θr v θ ,F tθ = F θr v r ,F tφ = − F rφ v r − F θφ v θ . (A4)This is analagous to defining E = − v × B . FromMaxwell’s equations we note that ∇× E = − ∂ t B , sodefining the electric field in this way should be inter-preted as defining the time rate of change of the magneticfield. Now that the condition of a perfectly conductingplasma has been met, we simplify the first integrand forthe rate of change of the field energy W ,t (Equation A2)using F αβ u β = 0 to find: F tβ J β = ( F bc v c ) J b = v c (cid:0) F βc J β − F tc J t (cid:1) = 1 ν v c v c . (A5)Therefore the rate of change of the field energy for anon force-free configuration under application of the mag-netofrictional method is given by: W ,t = − (cid:90) ν v A v A √ γd x. (A6)For ν < A φ via the conditionof a perfectly conducting plasma inevitably drives thefield energy towards an extremum and a force-free state.We now address the case of Ω F (cid:54) = 0, where the electricfield according to a distant observer does not originallyvanish. Treating this case separately is not necessary,but makes the above somewhat more transparent andallows us to explore the meaning of the field line angularvelocity. First, consider a subvolume over which Ω F maybe considered to be a constant. In that subvolume thevector potential A α = ( A t , A r , A θ , A φ ) is given by: A α = ( − Ω F A φ , A r , A θ , A φ ) . (A7)Suppose we now make a boost to a new frame via φ (cid:48) = φ − Ω F t . Then A α (cid:48) is given by: A α (cid:48) = (0 , A r , A θ , A φ ) . (A8)Noting that F tr = − A t,r and F tθ = − A t,θ , we see thatin the new frame F tr , F tθ , and Ω F all vanish (sugges-tive of the fact that Ω F may be interpreted as a mea-sure of the rotation of magnetic field lines referenced to zero angular momentum frames). In this frame we arefree to define an electric field that will evolve the mag-netic field towards a force-free configuration using theprocedure outlined above. Noting that t = t (cid:48) , r = r (cid:48) , θ = θ (cid:48) , and v A (cid:48) = v A we find that A φ (cid:48) ( r (cid:48) , θ (cid:48) ) = A φ ( r, θ ),so A φ (cid:48) ,t (cid:48) = − v A (cid:48) A φ (cid:48) ,A (cid:48) = − v A A φ,A = A φ,t . Thereforeapplication of the magnetofrictional method in the origi-nal frame inevitably moves the vector potential towards aforce-free configuration for all field line angular velocities. Appendix B: Expanded Magnetofrictional Terms
In this section we expand the advection equation ofthe magnetofrictional method, A φ,t = − v A A φ,A , in orderto explicitly show the form of the equations being used.We begin by noting that the advection equation may berewritten as: A φ,t = − νf (cid:0) g rr A φ,r + g θθ A φ,θ (cid:1) . (B1)Here we have exploited the fact that the coordinate ve-locity may be rewritten as v A = − νF Aβ J β = νf A φ,A for a function f ( r, θ, A φ , Ω F , √− gF θr ). We then weightthe coefficient of friction ν by a measure of the poloidalmagnetic field strength: ν = ν | B p | . (B2)Here ν is a constant; it is negative so that the mag-netofrictional method converges and its magnitude is se-lected for numerical stability. We define the magnitudeof the poloidal field | B p | as: | B p | ≡ θ (cid:0) A φ,θ + ∆ A φ,r (cid:1) . (B3)This weighting of ν is chosen primarily for convenience,in that it prevents regions with relatively sharp gradientsin A φ (i.e. large poloidal field) from evolving too rapidlyand its factor of ∆ removes the coordinate singularity onthe horizon. Using this weighting, the poloidal velocities v r and v θ are given by: v r = − ν f Σ | B p | ∆ A φ,r ,v θ = − ν f Σ | B p | A φ,θ . (B4)The common prefactor may be expanded as: f Σ | B p | = 1 A φ,θ + ∆ A φ,r π Σ (cid:20) C B φ ddA φ (cid:0) √− gF θr (cid:1) + C r A φ,r + C rr A φ,rr + C Ω r Ω F ,r + C θ A φ,θ + C θθ A φ,θθ + C Ω θ Ω F ,θ (cid:21) . (B5)2The factor C B φ is given by: C B φ = −
12 Σ . (B6) C r is given by: C r = − Σ∆ α ,r = − m ∆Σ (cid:0) r − a cos θ (cid:1) + 4 am ∆Σ (cid:0) r − a cos θ (cid:1) sin θ · Ω F + 2∆Σ (cid:0) r + 2 a r cos θ − a mr sin θ + a r cos θ + a m cos θ sin θ (cid:1) sin θ · Ω − amr ∆ sin θ · Ω F ,r + ∆ (cid:2) r + a r (cid:0) − sin θ (cid:1) + 2 a mr sin θ + a cos θ (cid:3) sin θ · F Ω F ,r . (B7) C rr is given by: C rr = − Σ∆ α = − ∆ (cid:0) r − mr + a cos θ (cid:1) − amr ∆ sin θ · Ω F + ∆ (cid:2) r + a r (cid:0) − sin θ (cid:1) + 2 a mr sin θ + a cos θ (cid:3) sin θ · Ω . (B8) C θ is given by: C θ = − Σ sin θ (cid:104) α sin θ (cid:105) ,θ = 1Σ (cid:2) r − mr + 2 a r cos θ − a mr (cid:0) − θ (cid:1) + a cos θ (cid:3) cos θ sin θ − amr Σ (cid:2) r + a (cid:0) θ (cid:1)(cid:3) sin θ cos θ · Ω F + 1Σ (cid:2) r + a r (cid:0) − θ (cid:1) + 6 a mr sin θ + a r (cid:0) − sin θ (cid:1) cos θ + a mr (cid:0) − θ (cid:1) sin θ + a cos θ (cid:3) sin θ cos θ · Ω − amr sin θ · Ω F ,θ + (cid:2) r + a r (cid:0) − sin θ (cid:1) + 2 a mr sin θ + a cos θ (cid:3) sin θ · F Ω F ,θ . (B9) C θθ is given by: C θθ = − Σ α = − (cid:0) r − mr + a cos θ (cid:1) − amr sin θ · Ω F + (cid:2) r + a r (cid:0) − sin θ (cid:1) + 2 a mr sin θ + a cos θ (cid:3) sin θ · Ω . (B10) C Ω r is given by: C Ω r = Σ∆ G φ A φ,r = Σ∆ ( g tφ + g φφ Ω F ) A φ,r = 2 amr ∆ A φ,r sin θ − ∆ (cid:2) r + a r (cid:0) − sin θ (cid:1) + 2 a mr sin θ + a cos θ (cid:3) sin θ · Ω F A φ,r . (B11) C Ω θ is given by: C Ω θ = Σ G φ A φ,θ = Σ ( g tφ + g φφ Ω F ) A φ,θ = 2 amrA φ,θ sin θ − (cid:2) r + a r (cid:0) − sin θ (cid:1) + 2 a mr sin θ + a cos θ (cid:3) sin θ · Ω F A φ,θ . (B12)In this form it is obvious that there are no divergences onthe horizon; as mentioned above this is a consequence ofscaling ν by the magnitude of the poloidal field. Thereis however a sin θ divergence in C θ ; this is enforcing A φ,θ = 0 on the axis as a consequence axisymmetry. Notealso that inside the horizon C rr and C θθ differ in sign;this leads to the magnetofrictional method becoming in-trinsically anti-diffusive and numerically unstable there.As mentioned in the main text, outside the light surfacesan overall factor of − α inthe C rr and C θθ terms leads to anti-diffusive numericalinstabilities. We finally note that all derivatives of A φ used to evaluate v A may be taken using central finite dif-ferencing. A one-sided derivative appropriate to upwinddifferencing is only taken to evolve the A φ,A derivativein A φ,t = − v A A φ,A .The derivative of the toroidal field (and field line an-gular velocity, in cases where it is not uniform) is takento be an unknown function of A φ . Practically we im-plement it using a lookup table that is modified duringruntime to diminish any kinks that develop on the innerlight surface. Appendix C: Expanded Percent Terms
In this section we detail our measure of how force-free agiven configuration is. Note that the force-free equationsmay be written as: − F αβ J β = X α . (C1)For a configuration to be force-free we must have X α = 0.If the toroidal field √− gF θr is conserved along magneticfield lines, as it must be when implemented as a functionof A φ , then to within numerical error X t = X φ = 0 . The poloidal components of the momentum flux may bewritten as: X A = f · A φ,A . (C2)Here f is a function of r , θ , Ω F , √− gF θr , and A φ ; setting f = 0 yields the transfield equation. To measure how3force-free a given solution is, we must measure how closeto zero f is; explicitly, f is given by: f = − π Σ∆ sin θ [ D + D + D + D + D + D + D ] . (C3)The functions D i are given by: D = 12 Σ sin θ ddA φ (cid:0) √− gF θr (cid:1) = − sin θ Σ · C B φ ddA φ (cid:0) √− gF θr (cid:1) D = ∆ α ,r sin θA φ,r = − sin θ Σ · C r A φ,r D = ∆ α sin θA φ,rr = − sin θ Σ · C rr A φ,rr D = sin θ (cid:16) α sin θ (cid:17) ,θ A φ,θ = − sin θ Σ · C θ A φ,θ D = α sin θA φ,θθ = − sin θ Σ · C θθ A φ,θθ D = − ∆ G φ sin θA φ,r Ω F ,r = − sin θ Σ · C Ω r Ω F ,r D = − G φ sin θA φ,θ Ω F ,θ = − sin θ Σ · C Ω θ Ω F ,θ (C4) To establish a measure of force-freeness, we take the sum (cid:80) i D i and compare it to the (absolute) maximum D i term, and insist that the ratio fall below a given thresh-hold: (cid:15) = (cid:80) i D i Max ( | D i | ) . (C5)In this work we insisted that (cid:15) <
1% over the entiredomain. In practice our algorithm yields (cid:15) (cid:28)
1% overmost of the domain with (cid:15) ≈≈