Effects of Inner Alfv\'en Surface Location on Black Hole Energy Extraction in the Limit of a Force-Free Magnetosphere

TL;DR

This study investigates how the position of the inner Alfvén surface influences the structure and energy extraction mechanisms of black hole magnetospheres in the force-free limit, revealing two main classes with different energy transmission characteristics.

Contribution

It provides a comprehensive numerical analysis of magnetospheres with varying Alfvén surface locations, highlighting their impact on energy extraction and jet formation near black holes.

Findings

Jet-like structures form near the horizon for Alfvén surfaces close to the ergosphere boundary.

Energy is directed towards the equatorial plane when the Alfvén surface is near the horizon.

Different Alfvén surface locations may lead to distinct observational signatures in high-energy phenomena.

Abstract

An energy extracting black hole magnetosphere can be defined by the location of its inner Alfv\'{e}n surface, which determines the rate of black hole energy extraction along a given magnetic field line. We study how the location of the inner Alfv\'{e}n surface can modify the structure of energy extracting black hole magnetospheres in the force-free limit. Hundreds of magnetospheres are numerically computed via a general relativistic extension of the Newtonian magnetofrictional method for a full range of black hole spins and flow parameters. We find that jet-like structures naturally form very close to the horizon for Alfv\'{e}n surfaces near the boundary of the ergosphere and that energy is extracted towards the equatorial plane for Alfv\'{e}n surfaces close to the horizon. This suggests that two broad classes of energy extracting black hole magnetospheres might exist; those that transmit extracted energy directly to distant observers, and those that transmit extracted energy to nearby accreting matter. Applied to transient high energy phenomena, we find that they might also differ in timescale by a factor of 20 or more.

Equations

$$d{s}^2$$

$$-\Sigma d{\theta }^2-\frac{A{\sin }^2\theta }{\Sigma }d{\varphi }^2,$$

$$\Sigma$$

$$\Delta$$

$$A$$

$${\left(\sqrt{-g}{F}^{\alpha \beta }\right)}_{,\beta }=-4\pi \sqrt{-g}{J}^{\alpha },$$

$$F_{\alpha \beta ,\gamma }+F_{\beta \gamma ,\alpha }+F_{\gamma \alpha ,\beta }=0.$$

$$F_{\alpha \beta }=A_{\beta ,\alpha }-A_{\alpha ,\beta }.$$

$$T^{\alpha \beta }$$

$$=\left(\rho +p\right)u^{\alpha }{u}^{\beta }-p{g}^{\alpha \beta }$$

$$+\frac{1}{4\pi }{g}^{\alpha \mu }{F}_{\mu \lambda }{F}^{\lambda \beta }+\frac{1}{16\pi }{g}^{\alpha \beta }{F}_{\mu \lambda }{F}^{\mu \lambda },$$

$$T^{\alpha \beta }{}_{;\beta }=0.$$

$$B^{\alpha }{\Psi }_{,\alpha }=0.$$

$$B^{\alpha }=\frac{1}{2\sqrt{-g}}{ϵ}^{\alpha \beta \mu \nu }{k}_{\beta }{F}_{\mu \nu },$$

$$F_{r\varphi }{\Omega }_{\text{F}}=F_{tr},$$

$$F_{\theta \varphi }{\Omega }_{\text{F}}=F_{t\theta }.$$

$$E$$

$$L$$

$$0<{\Omega }_{\text{F}}<{\omega }_{\text{H}}.$$

$$\sqrt{-g}{F}^{\theta r}=\pm \frac{\left(r_{+}^2+a^2\right)\left({\Omega }_{\text{F}}-{\omega }_{\text{H}}\right)\sin \theta }{\Sigma }{F}_{\theta \varphi }.$$

$$0$$

$$+\frac{1}{\sin \theta }{\left[\alpha F_{r\varphi }\right]}_{,r}+\frac{1}{\Delta }{\left[\frac{\alpha }{\sin \theta }{F}_{\theta \varphi }\right]}_{,\theta }$$

$$-\frac{G_{\varphi }}{\sin \theta }\left(F_{r\varphi }{\Omega }_{\text{F},r}+\frac{1}{\Delta }{F}_{\theta \varphi }{\Omega }_{\text{F},\theta }\right),$$

$$\alpha$$

$$G_{\varphi }$$

$$A_{\varphi }$$

$${\Omega }_{\text{F}}$$

$$\sqrt{-g}{F}^{\theta r}$$

$$h^{\prime }=-g\sin \theta .$$

$$A_{\varphi }$$

$${\Omega }_{\text{F}}$$

$$\sqrt{-g}{F}^{\theta r}$$

$$A_{\varphi }$$

$${\Omega }_{\text{F}}$$

$$\sqrt{-g}{F}^{\theta r}$$

$$-\{\frac{B_0{a}^3{\sin }^2\theta }{32m^4}[1+\frac{1}{2}\left(1-4R_2\right){\sin }^2\theta$$

$$+8R\left(r\right){\cos }^2\theta ]\},$$

$$R\left(r\right)$$

$$+\frac{m^3}{3r}+\frac{mr}{2}-\frac{r^2}{2}]$$

$$+\left(\frac{2r^3-3m{r}^2}{8m^3}\right)[{\text{Li}}_2\left(\frac{2m}{r}\right)$$

$$-\ln \left(1-\frac{2m}{r}\right)\ln \left(\frac{r}{2m}\right)],$$

Full text

I Introduction

When a rotating black hole is immersed in a magnetic field supported by a plasma, it is possible for that magnetic field to carry an electromagnetic Poynting flux of energy and angular momentum away from the black hole to a distant observer. Conservation of energy and angular momentum requires that the spin of the black hole be reduced; it is therefore said that such a magnetosphere extracts rotational energy from the black hole.

This type of electromagnetic black hole energy extraction is often referred to as the Blandford-Znajek mechanism after Blandford and Znajek (1977) who described it in terms of a stationary, axisymmetric, and force-free black hole magnetosphere. They showed that such a magnetosphere can efficiently extract an enormous amount of energy, and ever since black hole energy extraction has been considered as a possible mechanism powering a range of astrophysical phenomena, from stellar mass black holes in gamma-ray bursts to supermassive black holes in active galactic nuclei.

The Blandford-Znajek mechanism is not the only way to extract black hole rotational energy (Lasota et al. (2014) reviews others), nor is it strictly necessary for magnetic field lines to cross the horizon in order to extract black hole rotational energy (Komissarov (2004) and references therein). However the simple idea of a single magnetic field line connecting the horizon to distant observers can be seductive, as it is reminiscent of the more familiar and well understood case of a magnetic field line threading and torquing a rotating body. As one might naïvely expect from such an analogy, Blandford and Znajek (1977) showed that for black hole energy to be extracted along a given magnetic field line that field line must rotate more slowly than the black hole. In Takahashi et al. (1990) it was demonstrated that the condition of a relatively slowly rotating magnetic field line is actually the force-free limit of the requirement that a plasma inflow become super-Alfvénic inside the ergosphere.

Despite the defining role that the location of the inner Alfvén surface plays in determining whether or not a given magnetosphere extracts black hole rotational energy, it is currently unknown how the shape and location of the Alfvén surface might correspond to changes in the structure of a magnetosphere. The differences (if any) between a magnetosphere with an Alfvén surface near the horizon and a magnetosphere with an Alfvén surface near the boundary of the ergosphere are largely uncertain. This is due to the extreme intractability of the equations governing plasma flows in rotating black hole spacetimes, most especially the transfield equation that describes the force balance transverse to magnetic field lines (Nitta et al. (1991)). Currently no exact solutions of magnetospheres with plasma flows satisfying the transfield equation are known, so there is no robust analytic method of studying any potential correspondences between Alfvén surface location and magnetosphere structure.

In the force-free limit plasma inertial effects are ignored, the inner Alfvén surface coincides with the inner light surface described by the field line angular velocity of the magnetosphere, and the transfield equation simplifies to a purely electromagnetic condition. Due to those simplifications it is natural to consider the limit of force-free magnetospheres as a first step in attempting a study of the effects of Alfvén surface location. Unfortunately despite its relative simplicity the force-free limit of the transfield equation has very few analytic solutions, primarily limited to a handful of energy extracting magnetospheres found by perturbing in black hole spin (e.g. Blandford and Znajek (1977), McKinney and Gammie (2004), Pan and Yu (2015)) and a single class of fully exact but non energy extracting (and mostly unphysical) solutions (Menon and Dermer (2007) and Gralla and Jacobson (2014)). Therefore even in the force-free limit there is also no robust analytic method of studying the effects of Alfvén surface location.

Current numerical solutions and simulations also lack significant utility, as the nature of the Alfvén surface is ultimately an unknown function of the specific choice of boundary conditions and related initial assumptions that are made. This means that the Alfvén surface at best a result to be reported rather than an input to be explored, making it difficult to study how differing Alfvén surfaces might correspond to changes in magnetosphere structure.

The goal of this work was to study how Alfvén surface location might correspond to the structure of energy extracting black hole magnetospheres, thereby beginning to fill a gap in our understanding of black hole energy extraction. Due to the lack of useful analytic solutions it was necessary to calculate magnetospheres numerically. For simplicity we focused on the limit of force-free magnetospheres. We also focused on magnetospheres with a boundary condition in the equatorial plane compatible with a monopolar geometry; this avoids the confounding affects of the arbitrarily large forcings alternative boundary conditions can exert on the magnetosphere. This boundary condition is also compatible with those used to calculate magnetospheres in Contopoulos et al. (2013) and Nathanail and Contopoulos (2014) as well as the perturbed monopole solutions originally calculated by Blandford and Znajek (1977), allowing us to more easily compare our results with those works.

To form an initial survey it was necessary to calculate a large number of magnetospheres for a wide range of potential black hole spins and field line angular velocities (ultimately $\sim 500$ magnetospheres were found), so we developed numerical techniques that are efficient at calculating magnetospheres with similar structures.

Those numerical techniques are described in Section III. In Section II we provide some background on both ideal and force-free magnetospheres as well as the assumptions and equations we use. In Section IV we state some of our results, which we discuss in more depth in Section V before concluding.

II Background and Assumptions

In the stationary and axisymmetric limit, force-free black hole magnetospheres are naturally described by three scalar fields: the toroidal vector potential $A_{\phi}(r,\theta)$, the field line angular velocity $\Omega_{\text{F}}(A_{\phi})$, and the toroidal magnetic field $B_{\phi}(A_{\phi})$. In order to form a valid solution those scalar fields must satisfy a single differential equation, the force-free limit of the general relativistic Grad-Shafranov (or transfield) equation. For black hole rotational energy to be extracted along a given field line that field line must rotate slower than the black hole; $0<\Omega_{\text{F}}<\omega_{\text{H}}$, where $\omega_{\text{H}}$ is the angular velocity of the horizon. In this section we provide a brief review of those results and the assumptions they rely on. With the primary exception of $B_{\phi}$, which we define with opposite sign, the notation used is largely compatible with that of Blandford and Znajek (1977) and Takahashi et al. (1990) to facilitate comparison with those works; significant discrepancies will be noted. Throughout we use covariant 4-vector notation; Thorne and MacDonald (1982) and Komissarov (2004) provide reviews of and translations to the 3+1 formulations often used in other works.

II.1 Core Assumptions

We begin by assuming a stationary, axisymmetric, uncharged, and isolated black hole such that the spacetime may be described by the Kerr metric in Boyer-Lindquist coordinates. Using a $(+,-,-,-)$ metric signature and units where $c=G=1$, this corresponds to the line element:

[TABLE]

where

[TABLE]

The black hole is then immersed in a perfectly conducting (ideal) plasma and electromagnetic fields that together lack sufficient energy density to modify the background spacetime. The plasma and electromagnetic field configuration is also assumed to be stationary and axisymmetric, with an axis of symmetry corresponding to the rotational axis of the black hole. The electromagnetic fields must be consistent with Maxwell’s equations in a coordinate basis:

[TABLE]

Here $F^{\alpha\beta}$ is the field strength tensor, $J^{\alpha}$ is the current four vector, $g$ is the metric determinant, and a comma denotes a partial derivative. We use Greek letters to denote spacetime indices $(\alpha=t,r,\theta,\phi$), lowercase Latin letters to denote spatial indices ($a=r,\theta,\phi$), and uppercase Latin letters to denote poloidal indices ($A=r,\theta$). The field strength tensor may be expressed in terms of an electromagnetic vector potential $A_{\alpha}$ as:

[TABLE]

In this form it is clear that the toroidal electric field vanishes ($F_{t\phi}=0$) due to the assumptions of stationarity and axisymmetry, so “toroidal field” always references a toroidal magnetic field. The stress energy tensor outside the horizon is taken to be a linear combination of plasma and electromagnetic effects:

[TABLE]

where $\rho$ and $p$ reference the proper energy density and pressure of the plasma, and $u^{\alpha}$ its four velocity.

We finally assume that there is no external forcing (energy and momentum is conserved) over any region of interest such that the divergence of the stress energy tensor vanishes everywhere:

[TABLE]

Everything that follows rests upon the above assumptions and conditions.

II.2 Field Aligned Conserved Quantities

The three scalar fields that describe a force-free magnetosphere (vector potential $A_{\phi}$, field line angular velocity $\Omega_{\text{F}}$, and toroidal magnetic field $B_{\phi}$) are field-aligned conserved quantities in the sense that for a scalar field $\Psi$ we have:

[TABLE]

Here $B^{\alpha}$ is the magnetic field measured by a distant observer, which we define as:

[TABLE]

where $k_{\beta}$ is the Killing vector associated with the stationarity of the spacetime.

The toroidal vector potential $A_{\phi}$ directly describes the poloidal magnetic field, making it a useful flux function to trace poloidal magnetic field lines. The conservation of $A_{\phi}$ rests upon stationarity and axisymmetry, so $A_{\phi}$ is also a valid flux function in non force-free magnetospheres.

The conservation of field line angular velocity $\Omega_{\text{F}}$ is the general relativistic extension of Ferraro’s law of isorotation. The field line angular velocity is defined in terms of the field strength tensor as:

[TABLE]

This states that all electric fields are the result of the rotation of a purely magnetic configuration with angular velocity $\Omega_{\text{F}}$ (referenced to zero angular momentum frames). The rigid rotation of individual magnetic field lines through the conservation of $\Omega_{\text{F}}$ is a consequence of stationarity, axisymmetry, and a perfectly conducting plasma, so it is also valid in non force-free magnetospheres.

The toroidal field $B_{\phi}$ is defined through the definition of the magnetic field $B^{\alpha}$ in Equation 8 as $B_{\phi}=\sqrt{-g}F^{\theta r}$. This has opposite sign to the definition in Blandford and Znajek (1977), so for the remainder of this work the toroidal field will be referred to as $\sqrt{-g}F^{\theta r}$ to avoid confusion. The conservation of the toroidal field may be understood by considering Poynting fluxes through the magnetosphere. The only electric field is the result of rotating a purely magnetic configuration, and is therefore both entirely poloidal and perpendicular to the poloidal magnetic field. Therefore all poloidal Poynting fluxes are weighted only by the toroidal magnetic field and aligned with the poloidal magnetic field. As such the toroidal field is naturally interpreted as a measure of a conserved flux of energy $E$ and angular momentum $L$ per unit field line:

[TABLE]

Note that these definitions differ from those in Takahashi et al. (1990); in that work they are scaled by a conserved particle flux $\eta$ as $E\rightarrow\eta E$ and $L\rightarrow\eta L$. The conservation of the toroidal field along magnetic field lines is a consequence of stationarity, axisymmetry, and a perfectly conducting force-free plasma; in non force-free magnetospheres $E$ and $L$ include plasma parameters and $\sqrt{-g}F^{\theta r}$ is not conserved.

Varied derivations and discussions of these and other field-aligned conserved quantities may be found in Blandford and Znajek (1977), Bekenstein and Oron (1978), Camenzind (1986), and Takahashi et al. (1990).

II.3 Critical Points and Energy Extraction

In this section we explore key points along magnetic field lines and examine the conditions required for a magnetosphere to extract black hole rotational energy. We begin by making two assumptions that are unnecessary in general but greatly simplify the discussion. First, we assume a single magnetic field line that extends directly from the horizon to a distant observer. Second, we insist that this field line be “simple” in the sense that the distance from the horizon decreases monotonically as one travels along the field line towards the black hole; a useful mental image is a straight line connecting the horizon with distant regions.

Along any such field line there are seven key points, the first being the separation point. Inside the separation point gravitational forces dominate and draw plasma inwards to the horizon; outside the separation point centrifugal forces in the rotating frame of the magnetic field line dominate and accelerate plasma away from the black hole. The remaining six points are the ingoing and outgoing slow magnetosonic, Alfvén, and fast magnetosonic points through which the plasma accelerates as it travels from the separation point to the horizon or distant regions. In the force-free limit the ingoing slow magnetosonic point coincides with the separation point, the ingoing Alfvén point coincides with the inner light surface, and the ingoing fast magnetosonic point coincides with the horizon.

The inner and outer light surfaces are defined by the inner and outer surfaces exterior to the horizon where $\alpha=0$ ($\alpha$ is defined in Equation II.4 in the next section). The outer light surface is analagous to a pulsar light cylinder, where rigid rotation would cause a particle stuck to a magnetic field line to rotate at the speed of light. The particle would also rotate at the speed of light at the inner light surface as a consequence of gravitational time dilation, and it would be subluminal between the light surfaces.

A field line extracts black hole rotational energy when its ingoing Alfvén point is located inside the ergoregion. In the force-free limit this reduces to the requirement that the magnetic field line rotate more slowly than the black hole:

[TABLE]

Here $\omega_{H}=a/2mr_{+}$ is the angular velocity of a zero angular momentum observer on the horizon. In general there exist regularity conditions at each of the slow magnetosonic, Alfvén, and fast magnetosonic points. In the force-free limit the most useful of these is that of the ingoing fast magnetosonic point, which reduces to a condition on the fields at the horizon:

[TABLE]

The equations in use are insensitive to the sign of the toroidal field, but if an ingoing observer on the horizon is to measure finite electromagnetic fields Znajek (1977) showed that the negative branch must be chosen. This “Znajek condition” does not place any additional restrictions on the fields. It was already present in our core assumptions, most importantly the conservation of energy and momentum in the vanishing divergence of the stress energy tensor. However it can be a simple check to ensure that a solution is valid on the horizon.

Further discussion of the key points along a magnetic field line may be found in Takahashi et al. (1990) and Takahashi (2002).

II.4 The Force-Free Transfield Equation

In our core assumptions we demanded that there be no external forcing on the magnetosphere, as stated in the requirement that $T^{\alpha\beta}{}{}_{;\beta}=0$. If the poloidal components of this requirement are satisfied then the temporal and toroidal components are automatically satisfied. We therefore focus on the two poloidal components, which may be reduced to a single differential equation; this equation is often referred to as the “transfield equation” as it encapsulates the force balance transverse to poloidal magnetic field lines. In the force-free limit the transfield equation is given by (cf. Equation 3.14 of Blandford and Znajek (1977)):

[TABLE]

where

[TABLE]

The function $\alpha^{-1/2}$ may be interpreted as the “gravitational Lorentz factor” of a particle stuck to a magnetic field line, describing the effects of both rotation and gravitational time dilation, while $G_{\phi}$ may be interpreted as a measure of that particle’s rotational velocity with respect to the spacetime.

Solving the transfield equation is the core difficulty of finding black hole magnetospheres. In the stationary and axisymmetric limit of a rotating black hole, only one exact class of solutions is currently known (first published by Menon and Dermer (2007) and discussed in a more general framework by Gralla and Jacobson (2014)):

[TABLE]

Here $B_{0}$ is a scalar measure of the strength of the poloidal magnetic field; $h$ and $g$ are arbitrary functions of $\theta$ that must satisfy (prime denoting a derivative with respect to $\theta$):

[TABLE]

Solutions of this type generically contain an unphysical divergence in $\Omega_{\text{F}}$ in both the $\theta\rightarrow 0$ and low black hole spin $a\rightarrow 0$ limits, and are therefore mostly mathematical curiosities. They also do not extract rotational energy from the black hole along any field lines, further limiting their interest. However we still find them to be useful tools in determining the potential precision of a given numerical grid; see Sections III.4 and IV.4.

In the limit of a non-rotating ($a=0$) Schwarzschild black hole, there does exist a more reasonable variety of exact rotating “monopolar” solutions, given by:

[TABLE]

Such solutions do not extract rotational energy from the black hole (it has none to extract) but nontheless serve as a convenient starting point for exploring black hole energy extraction. These solutions are sometimes stated as “split-monopole” solutions by changing the sign of $B_{0}$ below the equatorial plane and appealing to an equatorial current sheet. For simplicity we will always refer to them as “monopolar” solutions, as regardless of interpretation there is still reflection symmetry across the equatorial plane. However it should be noted that “monopolar” references only the poloidal magnetic field; there is often a significant toroidal field, as can be seen in the $\Omega_{0}\neq 0$ solutions. Setting $\Omega_{0}=0$ for a static monopole and applying that solution to a slowly rotating black hole by perturbing in black hole spin (as in Blandford and Znajek (1977) and more recently Pan and Yu (2015)) yields:

[TABLE]

where

[TABLE]

and $\text{Li}_{2}$ is the dilogarithm, defined as:

[TABLE]

The higher order corrections by Pan and Yu (2015) (and McKinney and Gammie (2004)) are shown in curly brackets. It is possible to extend the above to include $\mathcal{O}(a^{5})$ corrections, but the solutions achieve a level of complication far surpassing the point of reason. There are two things to note from the above expressions; first, and most obviously, is that attempting to find perturbed solutions is not a trivial undertaking. Perturbing one of the simplest possible magnetic field configurations is effectively impossible past $\mathcal{O}(a^{2})$ corrections to $A_{\phi}$. The second thing to note that the leading order correction to $\Omega_{0}=0$ is $\Omega_{\text{F}}=0.5\omega_{\text{H}}$ (to first order in $a$). This has partly contributed to the “rule of thumb” that energy extracting black hole magnetospheres should rotate at one half the rate of the black hole ($\Omega_{\text{F}}\approx 0.5\omega_{\text{H}}$ also extracts the maximum amount of black hole energy for low spin).

Exploring the $\Omega_{0}\neq 0$ solutions was one of the initial questions that prompted our work, as they naturally encompass a full range of inner Alfvén surfaces and corresponding magnetospheres. Unfortunately analytically perturbing around such solutions is even more intractable than the already non-trivial $\Omega_{0}=0$ case, so we developed numerical techniques to solve for the structure of black hole magnetospheres. Discussion of those techniques is the topic of Section III.

II.5 Additional Assumptions

In addition to the core assumptions of Section II.1, we make two additional simplifying assumptions guided by our goal of studying correspondences between inner Alfvén surface location and magnetosphere structure. The first is the assumption of a monopolar geometry, in the sense that the poloidal magnetic field line threading the horizon on the equator must remain on the equator. This assumption is enforced via a boundary condition on the flux function $A_{\phi}$; we require that $A_{\phi}$ be constant in the equatorial plane. In making this assumption we were guided by a desire to extend the Schwarzschild monopole solutions of Equation II.4 to rotating spacetimes and to avoid assuming a specific model for nearby accreting matter. Without the presence of such matter to restrict magnetic field lines, a configuration with varying vector potential in the equatorial plane could violate the condition that there can be no closed horizon loops in stationary force-free magnetospheres (see MacDonald and Thorne (1982) and additional discussion in Gralla and Jacobson (2014)). Additonally, boundary conditions encoding an assumed model of matter in the equatorial plane or elsewhere have the potential to exert arbitrarily large forcings on the magnetosphere, confounding attempts to explore basic correlations between Alfvén surface location and magnetosphere structure. The interested reader may see Uzdensky (2005) for some complications induced by horizon-disk connections and Fendt (1997) for disk-jet connections.

The second assumption we made is that of uniform field line angular velocity, as that is the simplest possible method of classifying inner light surfaces (inner Alfvén surfaces in the force-free limit) and allowed us to know exactly where the inner light surface would be prior to calculating the structure of the field lines that pass through it. While it is convenient, this assumption restricts us to a region outside the horizon that is bounded by the outer light surface. This is due to the fact that our dissipative numerical techniques cannot generically find magnetospheres that pass smoothly through both inner and outer light surfaces under the assumption of uniform field line angular velocity, even if such solutions exist (we discuss smoothness at light surfaces in Section III.2). Ultimately we do not believe this to be nearly as limiting an assumption as it might initially appear, but it is still the most awkward assumption that we have made. As such we return to it in Section V.1 with more detailed discussion.

III Numerical Techniques

In order to numerically solve for a large number of black hole magnetospheres, we have extended the Newtonian magnetofrictional method developed by Yang et al. (1986) to the general relativistic regime. The magnetofrictional method makes an initial guess as to the structure of a force-free magnetosphere and then gradually relaxes the fields towards a valid force-free state. This makes it a useful tool for exploring our parameter space; once a solution for a given black hole spin and field line angular velocity is known adjacent solutions for slightly different values are readily found. In this section we discuss the underlying theory of the magnetofrictional method as well as the computational specifics of our implementation.

III.1 The Magnetofrictional Method

Once initial guesses have been made for a vector potential $A_{\phi}$, field line angular velocity $\Omega_{\text{F}}$, and toroidal field $\sqrt{-g}F^{\theta r}$, they may be inserted into the transfield equation to see if those guesses form a valid solution. This is equivalent to determining if the poloidal components of the divergence of the stress energy tensor vanish:

[TABLE]

If the initial guesses formed a valid and self consistent force-free solution then $X^{A}$ would vanish, but in general $X^{A}\neq 0$ and there is an excess momentum flux that may be interpreted as an external forcing. The core idea of the magnetofrictional method is to take that external forcing and convert it into the coordinate velocity of a fictitious plasma via friction. In this way a physically invalid (or at least inconsistent) initial guess for a force-free magnetosphere may be converted into a valid non force-free magnetosphere, and the fields may be evolved towards a force-free state by applying the equations governing inertial plasma flows. The first step in this process is to evaluate the excess momentum flux $X^{A}$, which is found as a modification to the left hand side of the transfield equation (Equation II.4):

[TABLE]

Once the excess momentum flux has been determined, it may be converted into the coordinate velocity of a fictitious plasma $v^{A}$ ($v^{r}=\partial r/\partial t$, $v^{\theta}=\partial\theta/\partial t$) via friction:

[TABLE]

Here the coefficient $\nu$ measures the strength of the friction. We then assume that the fictitious plasma is a perfect conductor such that there is no electric field in its rest frame:

[TABLE]

Here $u^{\alpha}$ is the four velocity of the plasma, which we define in terms of its coordinate velocity as $v^{c}\equiv u^{c}/u^{t}$. The above may then be used to find:

[TABLE]

If we insert this relationship into Maxwell’s homogeneous equations (Equation 3) we find the relativistic analog of the ideal induction equation $\partial_{t}\mathbf{B}=\mathbf{\nabla}\mathbf{\times}(\mathbf{v}\mathbf{\times}\mathbf{B})$:

[TABLE]

Inserting either $F_{\theta\phi}$ or $F_{r\phi}$ into the above then shows that the time rate of change of the vector potential $A_{\phi}$ is given by:

[TABLE]

In Appendix A we prove that application of this equation inevitably leads to a force-free configuration and in Appendix B we expand it in Boyer-Lindquist coordinates. Although the coordinate velocity $v^{A}$ is a complicated function of the metric and electromagnetic fields, at a high level the magnetofrictional method reduces to the straightforward evolution of an advection equation. A variety of numerical techniques for evolving such equations exist. For simplicity we have implemented an upwind differencing scheme similar to that described in Hawley et al. (1984). More sophisticated techniques might converge on a valid solution more rapidly, but we found upwind differencing to be both very numerically stable and fast enough for our purposes.

We also note that the relaxation scheme implemented by Uzdensky (2005) is similar to the method described above, in that it uses $A_{\phi,t}\sim v^{A}$; the primary difference is in the weighting of $v^{A}$.

III.2 Fixing Kinks

The magnetofrictional method described above is limited by the development of kinks in the magnetic field, which we now discuss. The transfield equation (Equation II.4) changes character at the inner and outer light surfaces where the gravitational Lorentz factor $\alpha$ changes sign. This change in character is transmitted to the coordinate velocity $v^{A}$ of the fictitious plasma, which goes as:

[TABLE]

On the inner and outer light surfaces $\alpha$ vanishes; it is positive between them and negative outside. To avoid having our numerical scheme become anti-diffusive and unstable, outside the light surfaces we introduce an overall multiplicative factor of $-1$ to maintain stability. This factor may be thought of as a correction for the fact that the fictitious plasma introduced by the magnetofrictional method is superluminal outside the light surfaces, and does not affect convergence.

However the change in sign of $\alpha$ does split the computational domain into three regions, and in general those regions do not join smoothly. In other words for an initial guess of field line angular velocity $\Omega_{\text{F}}$ and toroidal field $\sqrt{-g}F^{\theta r}$, the solutions that the magnetofrictional method finds for $A_{\phi}$ will in general not match across light surfaces. Treating $A_{\phi}$ as a flux function, this leads to discontinuous “kinks” in the magnetic field at the inner and outer light surfaces. As there are two light surfaces, both the field line angular velocity and the toroidal field could be evolved simultaneously in order to diminish both kinks, as in Contopoulos et al. (2013) and Nathanail and Contopoulos (2014).

However, as stated in Section II.5, for simplicity we have chosen to use a fixed and uniform field line angular velocity so we do not have the freedom required to generically diminish both kinks using the inherently dissipative magnetofrictional method. We therefore divide the computational domain into two regions, one inside the inner light surface and one between the light surfaces, and only diminish the kink at the inner light surface by evolving the toroidal field. The limitations and advantages of this approach are discussed in Section V.1.

In order to measure the kink at the inner light surface, the magnetofrictional method is applied separately to the regions on either side until an empirical “zero level” in the velocity $v^{A}$ is reached. A vanishing plasma coordinate velocity corresponds to $X^{A}=0$ and a valid force-free configuration, so this procedure simply finds a “close enough” solution in both regions. The two solutions are then compared across the inner light surface and the magnitude of any kink is measured in terms of $\Delta A_{\phi}$. This measurement is then used to modify the toroidal field in a method similar to that developed by Contopoulos et al. (1999) for light cylinders in pulsar magnetospheres:

[TABLE]

Here $\gamma$ is an empirical constant, usually much less than one, and $A_{\phi-}$ denotes the value of $A_{\phi}$ just inside the inner light surface. After the above equation has been used to make a small correction to the toroidal field, the magnetofrictional method is again applied to the regions on either side of the inner light surface for a fixed number of time steps, and then (if the set “zero level” in the plasma coordinate velocity is still achieved) the kink is re-evaluated and toroidal field re-adjusted. The cycle repeats until matching solutions are found.

Finding the correct “zero level”, value for $\gamma$, and number of time steps to apply the magnetofrictional method between modifications of the toroidal field is akin to critically damping a simple harmonic oscillator. Poor choices can lead to either an underdamped, wildly oscillating kink or an overdamped kink that requires an excessive amount of computational time to diminish. We have found that a “zero level” of $10^{-4}$, a factor of $\gamma=0.0005$, and $3000$ time steps between toroidal field changes is typically a safe starting point, but optimal values can vary significantly depending upon everything from grid resolution to the current magnitude and nature of the kink and how good the initial guess for the toroidal field is.

III.3 Measuring Force-Freeness

In principle the magnetofrictional method could be applied to the limits of numerical precision. However even if this were practical the approximations and assumptions we have made do not justify that level of precision. Instead we evolve the fields until a “good enough” level is reached, which we now describe.

A force-free solution is a solution with a vanishing Lorentz force, $F_{\alpha\beta}J^{\beta}=0$. Therefore it is reasonable to take a ratio of the Lorentz force to the magnitude of electric current and field strength tensor in order to measure how force-free a solution is, an approach taken in McKinney and Gammie (2004):

[TABLE]

If $\xi\ll 1$ the configuration can be thought of as being force-free. Unfortunately such a technique is of limited utility in our case; many of our magnetospheres contain regions of vanishing or very small current $J_{\alpha}J^{\alpha}\ll 1$, leading to divergences in $\xi$.

In force-free magnetospheres with regions of vanishing current, one alternative to $\xi$ is to measure the ratio of the Lorentz force to the forces of magnetic pressure and tension, as in Malanushenko et al. (2014):

[TABLE]

As with $\xi$, if $\zeta\ll 1$ then the configuration may be taken to be force free. Unfortunately some of our magnetospheres contain very strong monopolar components for which the forces of magnetic pressure and tension vanish separately, again leading to divergences. We have therefore developed a more mathematical and less physically motivated measure of force-freeness. We begin by noting that the transfield equation may be broken up into seven terms:

[TABLE]

The separation of the transfield equation into $D_{i}$ components is done somewhat arbitrarily by grouping factors of $A_{\phi}$, field line angular velocity, and the toroidal field; $D_{1}\sim\sqrt{-g}F^{\theta r}$, $D_{3}\sim A_{\phi,rr}$, etc. The exact separations we use are expanded in Appendix C. Mathematically the above puts us in the position of having an expression of the form $\sum D_{i}=\delta$ and requiring that $\delta$ be “small”. As a measure of force-freeness, we therefore take the absolute maximum of the seven $D_{i}$ terms and require that $\delta$ be a factor of $\epsilon$ smaller than it:

[TABLE]

While such a method is disadvantaged by not having an obvious physical interpretation, it does not diverge when applied to magnetospheres with regions of vanishing current or with magnetic fields containing strong monopolar components. For this work we have selected $\epsilon=1\%$ as a level of sufficient force-freeness. In practice our numerical techniques result in only a few segments of the inner light surface approaching $\epsilon=1\%$; over the rest of the domain $\epsilon$ is significantly smaller.

III.4 Analytic Comparisons

In Section II.4 we discussed some analytic solutions to the force-free transfield equation. Here we select four specific solutions that are useful for making comparisons with our numerical solutions; in Section IV.4 we compare them to a numerical solution with $\Omega_{\text{F}}=0.5\omega_{\text{H}}$ and $a=0.3m$ in order to examine the consequences of the $\epsilon=1\%$ error level discussed in the previous section.

The first analytic solution that provides a useful comparison is the exact Schwarzschild monopole described by Equation II.4. The field line angular velocity $\Omega_{0}$ of this solution may be set to correspond to the field line angular velocity of a given numerical solution. The monopole weight $B_{0}$ is set to $B_{0}=1$ in order to match the numerical solution’s change in $A_{\phi}$ between the axis and the equator (discussed in the next section). This solution provides a comparison with all field line angular velocities but is most useful for very small spins.

The second useful solution is the first order perturbed monopole solution (Equation II.4 without the terms in curly brackets). This solution is very similar to the Schwarzschild monopole with $\Omega_{0}=0.5\omega_{\text{H}}$, but differs slightly due to the fact that its field line angular velocity is only equivalent to $0.5\omega_{\text{H}}$ to first order in spin parameter $a$. As with the Schwarzschild monopole solution we set $B_{0}=1$. The field line angular velocity of this solution is fixed, so it is most useful for comparison with numerical solutions that have $\Omega_{\text{F}}=0.5\omega_{\text{H}}$ and small black hole spins.

The third solution is the third order perturbed monopole solution described by Equation II.4. As with the first order solution we set $B_{0}=1$. This solution is most useful for comparison with numerical solutions of $\Omega_{\text{F}}=0.5\omega_{\text{H}}$ for slightly higher spins than the first order solutions, and for determining how a higher order perturbation in spin translates to changes in our measure of error $\epsilon$. It should be noted that the field line angular velocity of this solution is not uniform, as in the first order solution, but is a function of both $r$ and $\theta$.

The fourth solution is a fully exact solution, described by Equation II.4, with $h=\sin^{3}\theta$ and $B_{0}=1$. There are infinitely many possible choices for $h$. Our selection is motivated by a desire for physical plausibility, which is most clearly seen by taking the Newtonian limit of the fields described. In that limit, our choice of $h$ corresponds to magnetic and electric fields given by:

[TABLE]

Here the vector components correspond to a standard orthonormal spherical basis. We emphasize that the factor of $a$ above is from the definition of the fields; if the $a\rightarrow 0$ limit appropriate to a Newtonian transition were also taken in the electromagnetic fields they would diverge. In the above form it is clear that our choice of $h=\sin^{3}\theta$ was made to correspond to physically plausible fields, in the limited sense that they don’t diverge anywhere in our domain; less careful choices generally diverge on the azimuthal axis. This solution is most useful for estimating the potential precision of a given numerical grid for arbitrary black hole spin; we will refer to it as the HS3 solution when we compare it to our numerical results.

III.5 Computational Specifics

In this section we detail the computational specifics of our numerical methods, covering the grid resolutions used, domain boundaries and boundary conditions, initial conditions, and performance.

The poloidal plane is divided into an $(r,\theta)$ grid for the vector potential $A_{\phi}$. In the radial direction the grid extends from just inside the horizon to a radius $r_{\text{max}}$ outside the ergosphere. The value of $r_{\text{max}}$ is dependent upon the location of the outer light surface, which is in turn dependent upon the choice of field line angular velocity $\Omega_{\text{F}}$ and black hole spin. For small values of $\Omega_{\text{F}}$ (near 0) or black hole spin the outer light surface can be a very large distance from the black hole, so we artificially place a cap of $20m$ on the maximum radius. For large values of $\Omega_{\text{F}}$ (near $\omega_{\text{H}}$) and black hole spin the outer light surface approaches the ergosphere and $r_{\text{max}}\approx 2m$ is used.

It is critical to resolve the inner light surface, so for large values of $\Omega_{\text{F}}$ where the inner light surface approaches the horizon we use very small grid spacings in the radial direction near the horizon, then gradually relax that spacing as $r$ increases. For small values of $\Omega_{\text{F}}$ the horizon and inner light surface are well separated and we use larger radial grid spacings near the horizon, then increase that spacing as $r$ extends past the ergosphere to relatively larger $r_{\text{max}}$ values. In general around 400 grid squares of varied spacing in the radial direction are used to resolve the inner light surface and ergoregion and an additional $450$ of varied spacing are used to extend to a given $r_{\text{max}}$. This yields a rough average of 850 grid squares in $r$ per full magnetosphere; exactly how many are used varies from magnetosphere to magnetosphere.

In order to treat the value of $A_{\phi}$ on the horizon as an evolving entity, we extend three grid squares past the horizon towards the black hole. We cannot extend any further as our implementation becomes anti-diffusive inside the horizon and more grid squares allow numerical instabilites to develop. We are using Boyer-Lindquist coordinates, which are singular on the horizon. This is not a difficulty, as we scale the magnetofrictional coefficient $\nu$ in Equation 23 by the magnitude of the poloidal magnetic field (see Appendix B). The coordinate singularity in that scaling cancels with the coordinate singularity in the transfield equation such that the horizon is well behaved. As a solution is found the horizon naturally settles into a configuration consistent with the Znajek horizon condition (the force-free limit of the fast magnetosonic regularity condition, Equation 12), justifying the usage of Boyer-Lindquist coordinates.

In the $\theta$ direction we extend from $\theta=0$ on the azimuthal axis to $\theta=\pi/2$ on the equatorial plane using $200$ linearly spaced grid squares. The equations in use do diverge on the azimuthal axis, a consequence of axisymmetry requiring that the magnetic field be purely radial there. By using the azimuthal axis as a fixed boundary condition we enforce axisymmetry and avoid that divergence, and find that double precision arithmetic is sufficient for the grid squares immediately adjacent to the axis. The other fixed boundary condition that we use is the equatorial plane, in keeping with the assumption of perturbing around a monopolar geometry (Section II.5). On the azimuthal axis we fix $A_{\phi}=4$ and on the equatorial plane we fix $A_{\phi}=3$, for a $\Delta A_{\phi}$ between them consistent with a monopole of unit weight ($A_{\phi}=\cos\theta$ in Equation II.4). In principle any arbitrary range between the axis and equatorial plane could be used; we made our selection for convenience.

We do not use fixed boundary conditions for the vector potential $A_{\phi}$ at either $r_{\text{min}}$ or $r_{\text{max}}$. Instead, after every time step we shoot outwards to find a boundary that keeps $A_{\phi}$ smooth and use that boundary for the next iteration. The reason for this is practical; although the axis and equatorial plane are fixed by symmetry, there is no restriction on $A_{\phi}$ for a given $r_{\text{max}}$. The inner boundary could be fixed by applying the horizon regularity condition, but we chose to have the horizon evolve as a simple check that the numerical algorithm is converging properly.

In order to calculate the derivatives in $A_{\phi}$ required to calculate the coordinate velocity $v^{A}$ of the fictitious plasma (Equations B and B) we use centered finite difference approximations appropriate to the local grid spacing. One-sided finite difference derivatives in $A_{\phi}$ appropriate to an upwind differencing algorithm are then used to evolve the magnetofrictional advection equation (Equation 27).

Our initial conditions for each run were fairly simple. Any smooth $A_{\phi}$ that decreased monotonically from the azimuthal axis to the equatorial plane could be used; we found no dependence upon initial conditions. Lack of monotonicity caused “spikes” to develop in $A_{\phi}$ that required magnetic reconnection via numerical diffusion to dissipate, greatly extending computation time when the code remained stable. To speed convergence we found it desirable to begin with the vector potential of the closest already calculated magnetosphere, “close” being defined by either slightly different black hole spin or field line angular velocity. We do not directly use the toroidal field $\sqrt{-g}F^{\theta r}$, instead we propose a function of $A_{\phi}$ corresponding to the derivative with respect to $A_{\phi}$ of the square of the toroidal field (the left hand side of Equation 29). From the transfield equation it can be shown that this function must vanish on the axis and equatorial plane (at $A_{\phi}=4$ or $A_{\phi}=3$ using the boundary conditions described above) but is otherwise largely unrestricted. Beyond increasing computation time we found no dependence upon initial choice of this function. In order to decrease computation time we often found it desirable to use three neighboring magnetospheres with the same spin but different field line angular velocities to “shoot” an initial guess for this function for a given value of $A_{\phi}$.

The overall number of time steps required to find a solution is sensitive to the accuracy of the initial guesses for both the vector potential and toroidal field, with the majority of computation time typically being spent reducing the kink at the inner light surface by finding a compatible toroidal field. In general with good initial guesses for the fields and well-tuned parameters a desktop computer can find a magnetosphere at the $\epsilon=1\%$ level (Section III.3) in a matter of hours; a $\epsilon\sim 10\%$ level can often be achieved in less than an hour. We have found that it is possible to significantly reduce computation time by bracketing the functional form of the derivative of the toroidal field using a kind of two dimensional root-finding algorithm with initial guesses taken from adjacent magnetospheres, but that technique requires careful setup and tuning to avoid instabilities. The naïve method of kink reduction described in Section III.2 is more robust, easier to implement, and while significantly slower is not prohibitively so for most purposes.

IV Results

We divide our results into four sections. First we explore the general structure of the magnetospheres as a function of black hole spin parameter and field line angular velocity. We then examine the rates of energy and angular momentum extraction. Lastly we explore the numerical error of our solutions.

IV.1 Field Line Structure

In this section we explore the behavior of poloidal magnetic field lines. We are limited to regions near to the black hole by our model assumptions and numerical techniques (Sections II.5 and III.2). However even over that limited domain interesting behaviors emerged, as shown in Figures 1 and 2.

Figure 1 plots poloidal magnetic field lines using spacing on the horizon corresponding to the strength of the radial magnetic field there (using $|B^{r}|\sim\csc\theta A_{\phi,\theta}$ appropriate to ingoing coordinates); denser field lines imply greater magnetic field strength. The first and last field line do not lie on the axis or equatorial plane, as those are fixed, but instead lie one grid square inwards on the horizon. The shading measures conserved momentum flux per unit field line ($L$ in Equation II.2), which is equivalent to the toroidal field $\sqrt{-g}F^{\theta r}$ multiplied by $4\pi$.

Figure 2 plots poloidal magnetic field lines using linear spacing in $\theta$ on the horizon. The first and last field line again lie one grid square inwards on the horizon. The shading measures conserved energy flux per unit field line ($E$ in Equation II.2).

The first thing to note is that as expected for very small black hole spin ($a\approx 0.1m$) all values of field line angular velocity yield magnetospheres that are almost indistinguishable from a Schwarzschild monopole (Equation II.4, noting that “monopole” refers only to the poloidal plane). For slightly larger values of spin ($a\approx 0.3m$), however, noticeable deviations from a monopolar structure rapidly emerge for both low and high field line angular velocities. As spin increases beyond $a\approx 0.3m$, only the $\Omega_{\text{F}}\approx 0.5\omega_{\text{H}}$ solutions remain roughly monopolar.

Low field line angular velocity solutions bend towards the azimuthal axis, with the strength of that bending increasing as black hole spin increases. For large spin ($a\approx 0.9m$) this bending can be severe, with the majority of poloidal magnetic field lines becoming nearly parallel to the azimuthal axis at a distance $r=20m$. The source of this bending is most easily seen in Figure 1; on the horizon the radial magnetic field is strong on the azimuthal axis and weak near the equatorial plane, while the toroidal magnetic field is strong on the equatorial plane and weak on the azimuthal axis. This discrepancy acts to naturally wind up the magnetic field and the magnetosphere self-collimates into a jet-like structure. The discrepancy becomes very large at high spin; for $a=m$ (not shown) the horizon radial magnetic field strength on the axis can be $\sim 10$ times that on the equator. For larger values of $\Omega_{\text{F}}$ (implying larger electric fields as viewed by distant observers, but not necessarily others such as ZAMOs) the toroidal field becomes much smaller, radial magnetic field strength becomes almost uniform on the horizon, and the field lines begin to bend toward the equatorial plane.

There are two separation surfaces shown in the figures (the separation surface may be defined as the point where $\alpha^{\prime}=0$, with the prime denoting differentiation along magnetic field lines); one for a monopolar geometry where $A_{\phi}$ is purely a function of $\theta$ and one corresponding to the numerically calculated $A_{\phi}$. For low field line angular velocities and high spins the separation surface moves away from the horizon at higher latitudes, while for high field line angular velocities it moves towards the horizon.

Figure 2 makes it clear that significantly more energy is extracted along field lines near the equator for all values of field line angular velocity, but low $\Omega_{\text{F}}$ and high spin magnetospheres redirect most of that extracted energy towards the azimuthal axis. We measure this behavior in the next section.

IV.2 Energy Extraction

In this section we explore the rate of black hole energy extraction. Energy flux is conserved along magnetic field lines, so the rate of black hole energy extraction is most easily calculated on the horizon. If we define $P$ as the power (energy per unit time as measured by a distant observer) leaving the horizon, we find (cf. Lasota et al. (2014)):

[TABLE]

As this is evaluated on the horizon, we may use the Znajek horizon condition (Equation 12) to find:

[TABLE]

Here $Q$ is a unitless scaling of the field line angular velocity to the horizon’s angular velocity; $\Omega_{\text{F}}=Q\omega_{\text{H}}$ (for energy extraction to take place we must have $0<Q<1$). In CGS units, the power $P$ becomes:

[TABLE]

where

[TABLE]

Here $a_{*}$ and $r_{+*}$ are dimensionless measures of black hole spin and horizon radius; $a=a_{*}m$ and $r_{+}=r_{+*}m$. The quantity $\chi$ is a dimensionless measure of the rate of black hole energy extraction that varies from magnetosphere to magnetosphere. The quantities $B_{x}$ and $r_{x*}$ are measures of monopolar magnetic field strength, in a sense that we now explore by considering the Newtonian limit.

In the Newtonian limit, a vector potential $A_{\phi}=B_{0}\cos\theta$ corresponds to a monopolar magnetic field that is given by:

[TABLE]

Here the vectors correspond to a standard orthonormal spherical basis. The weighting of $B_{0}$ on $A_{\phi}$ may therefore be interpreted as a measure of the monopolar magnetic field strength at a specified radius. For convenience our numerical solutions for $A_{\phi}$ assume unit spacing between the azimuthal axis and the plane ($B_{0}=1$), but any arbitrary weighting would yield identical results. The quantities $B_{x}$ and $r_{x*}$ are measures of that arbitrary weighting; $B_{x}$ is the field strength in Gauss of a Newtonian monopole at dimensionless radius $r_{x*}$. For example, $B_{x}=100\text{G}$ and $r_{x*}=20$ would correspond to a monopolar magnetic field strength of 100 Gauss at a distance $r=20m$ from a black hole of mass $m$. Some caution should be taken in applying this interpretation, however. None of our numerical solutions are exactly monopolar, and the strength of the magnetic field is an observer dependent quantity that does not have a simple translation from a flat space value to the spacetime of a rotating black hole. The above interpretation of $B_{x}$ and $r_{x*}$ is made purely for simplicity and convenience, and in general should be taken to be nothing more than a rough estimate.

In Figure 3 we plot the value of $\chi$ as a measure of the rate of black hole energy extraction for select values of black hole spin and field line angular velocity. We suppress spins greater than $a=0.99m$ as showing them would compress the curves of smaller spin that are of greater interest.

For a given value of $\Omega_{\text{F}}$, we note that increasing black hole spin always increases the rate of energy extraction but that changes in field line angular velocity can have a much larger effect than changing spin. The peak energy extracted for a given spin occurs at $\Omega_{\text{F}}\approx 0.5\omega_{\text{H}}$ for low spin and approaches $0.6\omega_{\text{H}}$ at high spin. The maximum energy extracted spans two orders of magnitude; $\chi_{\text{Max}}=0.03$ for $a\approx 0.99m$, $\chi_{\text{Max}}=0.003$ for $a\approx 0.5m$, and $\chi_{\text{Max}}=0.0003$ for $a\approx 0.2m$.

To a very good approximation, the functional form of $\chi(a)$ for fixed field line angular velocities $\Omega_{\text{F}}=0.1\omega_{\text{H}}$, $\Omega_{\text{F}}=0.95\omega_{\text{H}}$, and value for maximum energy extraction $\Omega_{\text{F Max}}$ are given by:

[TABLE]

In this form it is clear that for a given spin the maximum rate of energy extraction corresponding to $\Omega_{\text{F Max}}$ is roughly three times that of the lower values $\Omega_{\text{F}}\approx 0.1\omega_{\text{H}}$ or $\Omega_{\text{F}}\approx 0.95\omega_{\text{H}}$, an effect that can be much larger than changes in spin.

For low field line angular velocities the magnetic field bends towards the azimuthal axis, as shown in Figures 1 and 2. So while most energy is extracted on the horizon near the equatorial plane, a short distance from the horizon a large fraction of it ends up flowing outward along the azimuthal axis. In order to explore this behavior, we calculated the energy escaping through a spherical cone (both upper and lower hemispheres) at a radius $r=20m$ for $\Omega_{\text{F}}=0.1\omega_{\text{H}}$. The results are plotted in Figure 4. For high spins ($a_{*}\approx 0.85$ and greater), over 95% of the extracted energy escapes through a $45^{\circ}$ cone, and 80% escapes through a cone of less than $30^{\circ}$. This indicates that for high spin and low field line angular velocities a black hole magnetosphere is naturally inclined to extract energy via jet-like structures aligned with the rotational axis of the black hole.

IV.3 Angular Momentum Extraction

In this section we explore the rate of black hole angular momentum extraction. As with energy flux, angular momentum flux is conserved along magnetic field lines, and is most easily calculated on the horizon by exploiting the Znajek horizon condition. If we define $K$ as the rate of angular momentum extraction, we find:

[TABLE]

Note that this only differs from the rate of energy extraction (Equation IV.2) by a factor of field line angular velocity $\Omega_{\text{F}}$, consistent with the conserved energy $E$ and angular momentum $L$ per unit flux tube differing by the same factor (Equation II.2). In CGS units $K$ reduces to:

[TABLE]

where

[TABLE]

The primary structural difference between these expressions and those for the rate of energy extraction (Equations 37 and 38) is the fact that $\varphi$ does not vanish for $Q=0$. This is a statement that it is possible to spin down a black hole without extracting energy. In Figure 5 we plot the value of $\varphi$ as a measure of the rate of black hole angular momentum extraction for the same values of black hole spin and field line angular velocity that were used to plot the rate of energy extraction.

To a very good approximation, the functional form of $\varphi(a)$ for fixed field line angular velocities $\Omega_{\text{F}}=0.1\omega_{\text{H}}$, $\Omega_{\text{F}}=0.95\omega_{\text{H}}$, and value for maximum energy extraction $\Omega_{\text{F Max}}$ are given by:

[TABLE]

In this form it is clear that the $\Omega_{\text{F}}=0.1\omega_{\text{H}}$ solutions extract nearly twice as much angular momentum as the solutions that maximize the rate of energy extraction. It is also clear that the $\Omega_{\text{F}}=0.95\omega_{\text{H}}$ solutions extract around 10% of the angular momentum of the $\Omega_{\text{F}}=0.1\omega_{\text{H}}$ solutions, while from Equation IV.2 we note that these solutions extract roughly the same amount of energy. This indicates that high $\Omega_{\text{F}}$ solutions can extract energy from a black hole for a much longer period of time than lower $\Omega_{\text{F}}$ solutions can, as it will take longer for high $\Omega_{\text{F}}$ magnetospheres to spin down the black hole.

IV.4 Numerical Error Estimation

In this section we explore the numerical error of our solutions. In Section III.4 we discussed four analytic solutions that can provide useful comparisons with our numerical results, most especially in determining the reliability of the $\epsilon=1\%$ error level we have set (discussed in Section III.3 and Appendix C).

The three perturbed monopole and exact HS3 solutions have varying regions of applicability. We examine the numerical solution for $a=0.3m$ and $\Omega_{\text{F}}=0.5\omega_{\text{H}}$ so that we can reasonably compare all four solutions simultaneously. For convenience we begin by comparing them along a slice of constant $\theta=45^{\circ}$ from just inside the horizon to $r_{\text{max}}=20m$ (Figure 6). Along that slice it is apparent that the Schwarzschild monopole and first order perturbed monopole are very similar, with both typically having errors of less than 5%. The third order perturbed solution does better in general, with an error of around 1% near the inner light surface and inside the ergosphere. Aside from the two grid squares immediately adjacent to the inner light surface the numerical solution has an error of less than 0.05% across the entire slice, and is sometimes “better” than the HS3 solution’s error of around 0.0005%. We place quotes around “better” because the error of the HS3 solution is a rough measure of the precision of our numerical grid and derivatives; significantly exceeding $\sim 0.0005\%$ likely involves unjustified precision. Slices along different $\theta$ values as well as different spin parameters and field line angular velocities yield qualitatively similar results.

In order to examine the increase in error along the inner light surface we plot the percent error of the numerical solution in the two grid squares immediately adjacent to the inner light surface in Figure 7. This is done for all values of $\theta$, again for $a=0.3m$ and $\Omega_{\text{F}}=0.5\omega_{\text{H}}$. Even though substantial portions of the inner light surface are well below the 1% level, there are spikes up to 1%. Reducing those spikes is what takes the largest amount of computational time. The perturbed solutions are either comparable to or greatly exceed the error in the numerical solutions along the entirety of the inner light surface. The HS3 solution is typically 10-100 times better than any other solution. The exception is near the azimuthal axis where the HS3 solution’s field line angular velocity diverges. Note that the error of the HS3 solution does not also diverge, as the combined equations are well behaved there; the error increase comes from amplification of the error in taking numerical derivatives. The implications of the varying error levels are discussed in more depth in Section V.3.

In order to reduce computation time we make a guess as to the ultimate structure of the fields to use as an initial condition. It is therefore reasonable to ask if our guesses bias our results towards solutions that are close to that initial guess and ignore other potentially valid magnetospheres. There are two unknown quantities that we make guesses for; the toroidal component of the vector potential, $A_{\phi}(r,\theta)$, and the derivative of the square of the toroidal field with respect to the vector potential as an unknown function of the vector potential, $d/dA_{\phi}(\sqrt{-g}F^{\theta r})^{2}=F(A_{\phi})$. In order to examine how sensitive our solutions might be to changes in the initial form of these functions, we examine the magnetospheres obtained for black hole spin parameter $a=0.8m$; the numerically obtained functions of the derivative of the toroidal field for those magnetospheres are shown in Figure 8.

To asses how initial conditions might modify our results, we first coupled the numerical solution for the derivative of the toroidal field for the $\Omega_{\text{F}}=0.05\omega_{\text{H}}$ magnetosphere with the vector potential $A_{\phi}$ obtained for every other value of field line angular velocity (extrapolating outward to $r=20m$ for the higher field line angular velocity solutions). In every case we found that the numerical code rapidly converged to a solution essentially indistinguishable from the original $\Omega_{\text{F}}=0.05\omega_{\text{H}}$ solution. There were some minor deviations at large radii, as the $1\%$ error level was achieved there last (in the original solution it was achieved last on the inner light surface), but significantly less than the deviation between adjacent $\Omega_{\text{F}}=0.1\omega_{\text{H}}$ and $\Omega_{\text{F}}=0.0\omega_{\text{H}}$ magnetospheres.

As modifying the initial vector potential seemed to have no effect, we next examined modifying our initial guess for the derivative of the toroidal field. Poor guesses for this function can result in differences (“kinks”) in $A_{\phi}$ across the inner light surface that can easily exceed the difference in $A_{\phi}$ between the azimuthal axis and the equatorial plane. Diminishing inner light surface kinks takes the most computation time, so good guesses for the derivative of the toroidal field are the most critical initial condition and most likely source of any sensitivity our numerical solutions might have to initial conditions.

For clarity we focus on two alternative initial guesses for the functional form of the derivative of the toroidal field, shown in Figure 8 as initial conditions A and B. They are symmetric quadratics in $A_{\phi}$ that straddle the numerically obtained function, with the addition of periodic oscillations in Case A in an attempt to make a very poor initial guess without being overly ridiculous. Both cases led to very large initial kinks across the inner light surface and significantly increased the computation time required to find a solution. However the solutions obtained in both cases were again essentially indistinguishable from the original $\Omega_{\text{F}}=0.05\omega_{\text{H}}$ solution; the minor deviations in the derivative of the toroidal field (shown in in the bottom panel of Figure 8) are dwarfed by the deviation between adjacent $\Omega_{\text{F}}=0.1\omega_{\text{H}}$ and $\Omega_{\text{F}}=0.0\omega_{\text{H}}$ magnetospheres.

V Discussion

The most limiting assumption underlying our solutions is that of uniform field line angular velocity. In this section we discuss that assumption in more detail, consider two extremes of the magnetospheres that we found, and briefly explore the numerical error of our solutions.

V.1 Assumption of Uniform $\Omega_{\text{F}}$

As stated in Section II.5, we have assumed uniform field line angular velocities $\Omega_{\text{F}}$ in order to simplify the task of studying how the location of the inner Alfvén surface might correspond to changes in the structure of energy extracting black hole magnetospheres. This restricted us to solving for the structure of magnetospheres inside the outer light surface. With solutions in hand we return and examine that restriction in more detail. We begin by examining how limiting the assumption of uniform $\Omega_{\text{F}}$ might be in the two extreme classes of magnetospheres ($\Omega_{\text{F}}/\omega_{\text{H}}\approx 0$ and $\Omega_{\text{F}}\approx\omega_{\text{H}}$) that we found. We then discuss the existence and uniqueness of solutions that pass smoothly through both light surfaces. We close with a discussion of how limiting uniform $\Omega_{\text{F}}$ might be in interpreting our results.

For low field line angular velocities ($\Omega_{\text{F}}/\omega_{\text{H}}\approx 0$) where the magnetic field lines bend towards the azimuthal axis the limitations imposed by solutions restricted to the interior of the outer light surface should not be a significant concern. In that case the outer light surface can be many thousands of gravitational radii away from the horizon near the azimuthal axis and formally infinitely far away exactly on the axis. Diminishing the kink on such a distant surface could slightly modify the structure of the magnetosphere near the horizon, but in reasonable application we would generally expect deviations from our core assumptions of stationarity, axisymmetry, a perfectly conducting force-free plasma over such an extended region to be far more significant. Should that not be the case, we also see a smooth transition from our $\Omega_{\text{F}}=0.05\omega_{\text{H}}$ solutions to our $\Omega_{\text{F}}=0$ solutions for which an outer light surface does not exist (formally located infinitely far away from the black hole). We see no reason to expect pathologies in the limit $\Omega_{\text{F}}\rightarrow 0$, implying that our low $\Omega_{\text{F}}$ solutions are close to ones those that pass smoothly through both light surfaces.

For high field line angular velocities ($\Omega_{\text{F}}\approx\omega_{\text{H}}$) where the magnetic field lines bend towards the equatorial plane we would expect to find a connection to nearby accreting matter (not necessarily a thin disk limited to the equatorial plane). Our boundary condition of a single magnetic field line in the equatorial plane might be reasonable close to the horizon, but moving away from the horizon it would become increasingly likely to find significant deviations depending upon the specific model of nearby accreting matter that was chosen. Finding solutions that passed smoothly through both light surfaces would therefore be somewhat ill-advised, as it would involve modifying the near horizon magnetosphere using restrictions found in unreasonably distant regions. A more realistic magnetosphere would include a description of nearby accreting matter that included a model of plasma injection into an inflow interior to or near the separation surface. Our solutions could be representative of such inflows, especially a thick disk or torus near the ergosphere. A solution that passed smoothly through an outer light surface, on the other hand, would also need to include some model of plasma injection into an outflow consistent with the field line geometry obtained, which might be difficult.

In arguing that it is possible to usefully interpret the above limits on $\Omega_{\text{F}}$ despite the limitations of our domain, we have ignored the question of whether or not similar solutions that pass smoothly though both light surfaces even exist. Using essentially identical boundary conditions on the vector potential $A_{\phi}$, Contopoulos et al. (2013) and Nathanail and Contopoulos (2014) (hereafter CKP13 and NC14) modified both the field line angular velocity and toroidal field to find solutions that pass smoothly through both light surfaces. In doing so they found a single nearly monopolar solution with varying $\Omega_{\text{F}}\approx 0.5\omega_{\text{H}}$ that is very similar in structure to our uniform $\Omega_{\text{F}}=0.5\omega_{\text{H}}$ solutions. It is therefore natural to ask if the solutions found by CKP13 and NC14 are unique. The exact solutions of Equation II.4, while non energy extracting and largely unphysical, indicate that they are not; there are in fact infinitely many solutions fulfilling our boundary conditions on $A_{\phi}$ that pass smoothly through both an inner and outer light surface. The question then becomes why we found the solutions that we did regardless of initial conditions, and why CKP13 and NC14 similarly found a single solution.

The answer to that question lies in our numerical techniques. We show in Appendix A that the magnetofrictional method works by minimizing the energy in the electromagnetic fields. In minimizing the kink across the inner light surface we are finding matched minimum energy solutions. The apparent uniqueness of our results is a suggestion of the uniqueness of a minimum energy solution, not the uniqueness of our solution in full generality. The majority of numerical codes are common in the behavior of being dissipative and seeking minimum energy states; any code that allows energy to be added at will is likely to be numerically unstable. For example there is a full range of valid rotating Schwarzschild monopole solutions (Equation II.4), but it would be a very unique numerical code that would be capable of blindly finding all of them. Most codes would always converge on the minimum energy $\Omega_{0}=0$ solution.

For any value of black hole spin, the minimum energy magnetosphere consistent with our boundary conditions will be as close to monopolar as possible ($\Omega_{\text{F}}\approx 0.5\omega_{\text{H}}$). The bunching of magnetic field lines towards the axis or equator seen in our low and high $\Omega_{\text{F}}$ solutions requires the addition of energy to spin the magnetosphere away from that monopole. It is therefore not surprising that CKP13 and NC14 found a single roughly monopolar solution that passes smoothly through both inner and outer light surfaces; without the explicit and very careful addition of energy to move to and maintain a rotating magnetosphere any stable code would likely converge on that solution. Our utilization of uniform $\Omega_{\text{F}}$ may be interpreted as a way of exploring the structure of arbitrarily rotating magnetospheres using a dissipative numerical code, and we believe that a large fraction of our solutions could be taken to be good approximations of solutions that pass smoothly through both light surfaces. In realistic application, however, they should largely be taken to be representative of inflow solutions. Outflow solutions would require the additional description of plasma injection mechanisms, and even if our solutions passed smoothly through an outer light surface there is no guarantee that plasma parameters would be either continuous or conserved across an extended plasma injection region (as implied by single solutions that pass smoothly through both light surfaces).

We have argued that our limited domain might not prevent useful interpretations and that our solutions might approximate solutions that pass smoothly through both light surfaces (even if such solutions might be of dubious value). It could still be asked how representative our solutions might be of energy extracting black hole magnetospheres, as completely uniform $\Omega_{\text{F}}$ magnetospheres are unlikely to exist. We cannot fully answer that question, but we emphasize that exploring uniform $\Omega_{\text{F}}$ magnetospheres is not the goal of this work. The goal of this work is to explore the effects of inner Alfvén surface location. Uniform $\Omega_{\text{F}}$ in the force-free limit is merely a useful tool to do so, both in its obvious simplicity as well as in allowing us to know exactly where the inner light surface will be prior to solving for the structure of the magnetosphere that passes through it. Slightly deforming the location of the Alfvén surfaces used here will not have significant effects, and in general we expect the Alfvén surface to be continuous, so the solutions obtained here could be merged to explore more complex scenarios. For example, if one desired an ingoing solution with a connection to accreting matter near the equator and the launching of a jet-like structure from the horizon along the rotational axis, an appeal could be made to an Alfvén surface that lies closer to the boundary of the ergosphere at higher latitudes and closer to the horizon near the equator. On the horizon this would correspond to lower $\Omega_{\text{F}}$ for small values of $\theta$ and higher $\Omega_{\text{F}}$ for large values of $\theta$, compatible with the simulations conducted by McKinney and Gammie (2004).

V.2 Two Types of Magnetospheres

Our solutions indicate two extreme methods for black hole energy to be extracted and transmitted to distant observers. For low field line angular velocities ($\Omega_{\text{F}}/\omega_{\text{H}}\approx 0$) jet-like structures naturally form to transmit extracted energy directly from the horizon to distant observers along the azimuthal axis. For high field line angular velocities ($\Omega_{\text{F}}\approx\omega_{\text{H}}$) we find magnetospheres that are consistent with the transmission of energy from the horizon to nearby accreting matter, which in turn might transmit the extracted energy to distant observers. In this section we discuss potential consequences of those two types of magnetospheres, noting at that outset that any potential distinction between the two in realistic astrophysical scenarios is likely to be fuzzy.

If we examine the rate of energy and angular momentum extraction through Equations IV.2 and IV.3, we see that the two types might have significant time-dependent distinctions due to their differing rates of angular momentum extraction and concurrent black hole spindown. Therefore transient high energy phenomena powered by black hole energy extraction might also have two distinct signatures. For specificity, consider a gamma-ray burst (GRB). Black hole energy extraction is a candidate for the central engine of a GRB (e.g. Gehrels and Mészáros (2012)), and the exponential decay associated with black hole energy extraction might be compatible with some GRB observations (e.g. Nathanail et al. (2016)). We do not demand a specific type or model of GRB, however. We only require a transient high energy event that might be powered by black hole energy extraction, and “GRB” is a convenient specific stand-in for“transient high energy event” even if our treatment is too crude to completely describe or differentiate between realistic GRBs of any specific type.

If we take $B_{x}=10^{16}\text{G}$ and $r_{x}=3$ as constants in Equations 37 and 42 (i.e. a magnetic field strength of roughly $10^{16}\text{G}$ just outside the ergosphere), assume that $a=0.95m$ and $m=4M_{\odot}$ at time $t=0$, and for simplicity insist upon a state of suspended accretion (such as a magnetically arrested disk in a low accretion state), we find the rates of black hole energy extraction plotted in Figure 9.

We immediately note that all three cases exhibit roughly exponential decay, but the $\Omega_{\text{F}}=0.95\omega_{\text{H}}$ case decays relatively slowly. After roughly 20-30 seconds both the $\Omega_{\text{F}}=0.1\omega_{\text{H}}$ and $\Omega_{\text{F Max}}$ luminosities have decreased by factors of $\sim 10^{3}$ due to the black hole spin dropping below $a\approx 0.05m$. The rate of angular momentum extraction for the $0.95\omega_{\text{H}}$ case is far lower, however, and after 30 seconds its luminosity has not significantly decreased. Note that in this model larger black hole masses would take less time to spin down due to the increase in magnetic field strength implied by $r_{x}$ corresponding to a larger dimensional radius. For example, a black hole with twice the mass but otherwise identical parameters would be spun down in around half the time. Also note that the scaling on the axes is somewhat irrelevant; they were chosen for crude correspondence to a GRB, but a wide range of systems with different parameters, luminosities, and timescales exhibit qualitatively identical curves.

It takes over 350 seconds for the $0.95\omega_{\text{H}}$ luminosity to drop to the same level that the $0.1\omega_{\text{H}}$ luminosity fell to in around 20 seconds. We therefore expect that GRBs with magnetic field lines that directly connect the horizon to distant observers might be at least 10-20 times shorter than GRBs with magnetic field lines connecting to a disk or other matter before connecting to distant observers. The addition of more realistic effects such as accretion actively spinning up the black hole would obviously complicate matters, and as noted at the outset the boundary between both types of magnetospheres is unlikely to be very distinct. However, if GRBs (again emphasizing the limited sense in which we are using the term) are powered by black hole energy extraction it would not be unreasonable to expect two broad classes with differing characteristic durations as well as different radiation signatures; the longer lived $0.95\omega_{\text{H}}$ type might also exhibit a greater degree of thermalization due to the deposition of the extracted energy into inflowing matter before that energy is transmitted to distant observers.

The analysis presented here is too crude to reasonably claim that distinctions between GRBs can arise due to differences in the location of their inner Alfvén surface, encapsulated here in differences in field line angular velocity. A more reasonable model is beyond the scope of our present work, but the possibility that such a model could allow finer classification of some GRBs or other transient high energy phenomena via correlation of timescale and degree of thermalization is interesting.

V.3 Numerical Error Analysis

In this section we discuss how reliable the numerical calculations underlying our solutions might be.

With the possible exception of the kink at the inner light surface, we find that our numerical solutions are consistent with the precision of their numerical grid. For the case of $a=0.3m$ and $\Omega_{\text{F}}=0.5\omega_{\text{H}}$, the error of our solution is generally around 100 times smaller than the error in any of the perturbed monopole solutions that are considered to be good approximations at low spin. While the numerical error involved in implementing the exact HS3 solution (arising mostly from taking numerical derivatives) can be another 100 times smaller still, there are large regions in which the numerical solution is actually “better” than the HS3 solution and likely exceeds the precision justified by our numerical grid.

Along the inner light surface the error in our numerical solution 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 solution generally exceeds the numerical solution by another factor of 10-100. However significantly better solutions would be difficult to obtain. At the $1\%$ level the “kink” (change in $A_{\phi}$) in our solution across the inner light surface in the radial direction is around $10^{-6}$, much smaller than the required linear change in $A_{\phi}$ of $5\times 10^{-3}$ in the $\theta$ direction implied by our usage of 200 grid squares between the azimuthal axis (where $A_{\phi}=4$) and the equator (where $A_{\phi}=3$). This means that it becomes very difficult to accurately quantify the kink, much less diminish it further.

Some experimentation reveals that a much higher error level of around 10% would not have changed our results in any significant way. This is consistent with the error of the perturbed solutions, which should be accurate enough for most purposes and have similarly sized errors. Reducing the error of our solutions to the 1% level increased the required computation time significantly, and achieves a precision that almost certainly far exceeds what our assumptions allow for in any reasonable application. However reducing the error to the 1% level allows us to confidently state that any deficiencies in our results rest within our assumptions and not in our numerical calculations.

Due to the importance of selecting appropriate initial conditions in diminishing computation time, it might be asked if our initial guesses in any way affect or drive our results. As shown in Figure 8, even with somewhat ridiculous choices of initial condition we ultimately arrive at effectively identical solutions that are appreciably different from neighboring solutions. We cannot (and do not) claim that this implies that our solutions are unique in general; it is possible that there are many solutions that satisfy our boundary conditions and assumptions. The only potentially persuasive uniqueness in our solutions lies in unique minimum energy solutions that are matched across the inner light surface, as a consequence of the magnetofrictional method being an energy minimizing algorithm (shown in Appendix A). In developing our numerical code we conducted many more trials and experiments than are shown here, and we never observed any indication that it might be possible to converge on two different solutions. Therefore while we cannot prove that our solutions are unique matched minimum energy solutions consistent with our boundary conditions and assumptions, we do believe that it is a reasonable assumption to make.

VI Conclusions

The defining feature of energy extracting black hole magnetospheres is the location of their inner Alfvén surface, which coincides with the inner light surface in the force-free limit. Despite its importance, no comprehensive studies of the effects of modifying the location of the inner Alfvén surface have been accomplished. We have begun addressing that deficiency by studying how simple force-free magnetospheres can be modified as a function of uniform field line angular velocity, which together with black hole spin determines the location of the inner light surface. In order to do so we extended the Newtonian magnetofrictional method for computing force-free magnetospheres into the general relativistic regime, which allowed us to efficiently calculate hundreds of energy extracting black hole magnetospheres as functions of inner Alfvén surface location. In so doing we found that inner Alfvén surfaces near the horizon cause extracted energy to flow towards the equatorial plane, while inner Alfvén surfaces near the boundary of the ergosphere cause extracted energy to flow outwards via jet-like structures aligned with the azimuthal axis. Applied to transient high energy phenomena, those magnetospheres imply two timescales that might differ by a factor of 20 or more. This suggests that two classes of transient high energy phenomena might exist; shorter ones that directly connect the horizon to distant observers, and longer ones that have a disk or other significant matter separating the energy flow from the horizon to distant observers that potentially thermalizes the extracted energy, creating a different radiation signature.

Acknowledgments

K.T. gratefully acknowledges useful conversations with Dana Longcope. M.T. was supported by JSPS KAKENHI Grant Number 24540268.

Appendix A Proof of Convergence

In this section we show that the magnetofrictional method inevitably drives a given vector potential $A_{\phi}$ towards 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. Second, we show that the magnetofrictional method reduces the electromagnetic energy contained in a given volume, thereby inevitably driving that volume towards a minimum energy and force-free state. For simplicity we assume Boyer-Lindquist coordinates, that a force-free state exists, and that no pathologies develop (e.g. magnetic reconnection, if necessary, is assumed to occur via numerical diffusion).

We define the energy $W$ contained within the fields at a given coordinate time to be:

[TABLE]

Here $\sqrt{\gamma}d^{3}x$ is the appropriate proper volume element and $T^{\alpha\beta}$ is the electromagnetic stress energy tensor. Differentiating with respect to coordinate time, we find:

[TABLE]

Here $d\Sigma_{A}$ is the appropriate directed surface element bounding the region such that the second term is a measure of the net Poynting flux through the region’s boundary. We have assumed no field pathologies, so any field line entering the volume also exits the volume. Energy flux per unit field line is conserved, so there can be no net Poynting flux into the region along any magnetic field lines and the second term necessarily vanishes. In a force-free configuration $F_{\alpha\beta}J^{\beta}=0$ and the first term also vanishes; therefore a force-free configuration extremises the field energy $W$ in coordinate time.

The core conceit of the magnetofrictional method is that the excess momentum flux of a non force-free configuration may be converted into the coordinate velocity of a fictitious plasma; mathematically, this may be stated as:

[TABLE]

Here the coordinate velocity $v^{a}$ of the plasma is defined in terms of its four velocity $u^{a}$ as $v^{a}\equiv u^{a}/{u^{t}}$. Both the field line angular velocity $\Omega_{\text{F}}$ and the toroidal field $\sqrt{-g}F^{\theta r}$ are assumed to be functions of $A_{\phi}$ (i.e. conserved along magnetic field lines), so $v^{\phi}=0$. We now make the simplifying assumption that $\Omega_{\text{F}}$ vanishes (i.e. $F_{tr}=F_{t\theta}=0$ and there is no electric field from the perspective of a distant observer); we will address the $\Omega_{\text{F}}\neq 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_{\alpha\beta}u^{\beta}=0$:

[TABLE]

This is analagous to defining $\mathbf{E}=-\mathbf{v}\mathbf{\times}\mathbf{B}$. From Maxwell’s equations we note that $\mathbf{\nabla}\mathbf{\times}\mathbf{E}=-\partial_{t}\mathbf{B}$, so defining the electric field in this way should be interpreted as defining the time rate of change of the magnetic field. Now that the condition of a perfectly conducting plasma has been met, we simplify the first integrand for the rate of change of the field energy $W_{,t}$ (Equation A) using $F_{\alpha\beta}u^{\beta}=0$ to find:

[TABLE]

Therefore the rate of change of the field energy for a non force-free configuration under application of the magnetofrictional method is given by:

[TABLE]

For $\nu<0$ the right hand side is always negative (recall the signature of our metric) and using a fictitious plasma and electric field to evolve $A_{\phi}$ via the condition of a perfectly conducting plasma inevitably drives the field energy towards an extremum and a force-free state.

We now address the case of $\Omega_{\text{F}}\neq 0$, where the electric field according to a distant observer does not originally vanish. Treating this case separately is not necessary, but makes the above somewhat more transparent and allows us to explore the meaning of the field line angular velocity. First, consider a subvolume over which $\Omega_{\text{F}}$ may be considered to be a constant. In that subvolume the vector potential $A_{\alpha}=(A_{t},A_{r},A_{\theta},A_{\phi})$ is given by:

[TABLE]

Suppose we now make a boost to a new frame via $\phi^{\prime}=\phi-\Omega_{\text{F}}t$. Then $A_{\alpha^{\prime}}$ is given by:

[TABLE]

Noting that $F_{tr}=-A_{t,r}$ and $F_{t\theta}=-A_{t,\theta}$, we see that in the new frame $F_{tr}$, $F_{t\theta}$, and $\Omega_{\text{F}}$ all vanish (suggestive of the fact that $\Omega_{\text{F}}$ may be interpreted as a measure of the rotation of magnetic field lines referenced to zero angular momentum frames). In this frame we are free to define an electric field that will evolve the magnetic field towards a force-free configuration using the procedure outlined above. Noting that $t=t^{\prime}$, $r=r^{\prime}$, $\theta=\theta^{\prime}$, and $v^{A^{\prime}}=v^{A}$ we find that $A_{\phi^{\prime}}(r^{\prime},\theta^{\prime})=A_{\phi}(r,\theta)$, so $A_{\phi^{\prime},t^{\prime}}=-v^{A^{\prime}}A_{\phi^{\prime},A^{\prime}}=-v^{A}A_{\phi,A}=A_{\phi,t}$. Therefore application of the magnetofrictional method in the original frame inevitably moves the vector potential towards a force-free configuration for all field line angular velocities.

Appendix B Expanded Magnetofrictional Terms

In this section we expand the advection equation of the magnetofrictional method, $A_{\phi,t}=-v^{A}A_{\phi,A}$, in order to explicitly show the form of the equations being used. We begin by noting that the advection equation may be rewritten as:

[TABLE]

Here we have exploited the fact that the coordinate velocity may be rewritten as $v_{A}=-\nu F_{A\beta}J^{\beta}=\nu fA_{\phi,A}$ for a function $f(r,\theta,A_{\phi},\Omega_{\text{F}},\sqrt{-g}F^{\theta r})$. We then weight the coefficient of friction $\nu$ by a measure of the poloidal magnetic field strength:

[TABLE]

Here $\nu_{0}$ is a constant; it is negative so that the magnetofrictional method converges and its magnitude is selected for numerical stability. We define the magnitude of the poloidal field $\left|B_{\text{p}}\right|^{2}$ as:

[TABLE]

This weighting of $\nu$ is chosen primarily for convenience, in that it prevents regions with relatively sharp gradients in $A_{\phi}$ (i.e. large poloidal field) from evolving too rapidly and its factor of $\Delta$ removes the coordinate singularity on the horizon. Using this weighting, the poloidal velocities $v^{r}$ and $v^{\theta}$ are given by:

[TABLE]

The common prefactor may be expanded as:

[TABLE]

The factor $C_{B_{\phi}}$ is given by:

[TABLE]

$C_{r}$ is given by:

[TABLE]

$C_{rr}$ is given by:

[TABLE]

$C_{\theta}$ is given by:

[TABLE]

$C_{\theta\theta}$ is given by:

[TABLE]

$C_{\Omega r}$ is given by:

[TABLE]

$C_{\Omega\theta}$ is given by:

[TABLE]

In this form it is obvious that there are no divergences on the horizon; as mentioned above this is a consequence of scaling $\nu$ by the magnitude of the poloidal field. There is however a $\sin\theta$ divergence in $C_{\theta}$; this is enforcing $A_{\phi,\theta}=0$ on the axis as a consequence axisymmetry. Note also that inside the horizon $C_{rr}$ and $C_{\theta\theta}$ differ in sign; this leads to the magnetofrictional method becoming intrinsically anti-diffusive and numerically unstable there. As mentioned in the main text, outside the light surfaces an overall factor of $-1$ is added to the above in order to maintain numerical stability, as otherwise the sign of $\alpha$ in the $C_{rr}$ and $C_{\theta\theta}$ terms leads to anti-diffusive numerical instabilities. We finally note that all derivatives of $A_{\phi}$ used to evaluate $v^{A}$ may be taken using central finite differencing. A one-sided derivative appropriate to upwind differencing is only taken to evolve the $A_{\phi,A}$ derivative in $A_{\phi,t}=-v^{A}A_{\phi,A}$.

The derivative of the toroidal field (and field line angular velocity, in cases where it is not uniform) is taken to be an unknown function of $A_{\phi}$. Practically we implement it using a lookup table that is modified during runtime to diminish any kinks that develop on the inner light surface.

Appendix C Expanded Percent Terms

In this section we detail our measure of how force-free a given configuration is. Note that the force-free equations may be written as:

[TABLE]

For a configuration to be force-free we must have $X_{\alpha}=0$. If the toroidal field $\sqrt{-g}F^{\theta r}$ is conserved along magnetic field lines, as it must be when implemented as a function of $A_{\phi}$, then to within numerical error $X_{t}=X_{\phi}=0.$ The poloidal components of the momentum flux may be written as:

[TABLE]

Here $f$ is a function of $r$, $\theta$, $\Omega_{\text{F}}$, $\sqrt{-g}F^{\theta r}$, and $A_{\phi}$; setting $f=0$ yields the transfield equation. To measure how force-free a given solution is, we must measure how close to zero $f$ is; explicitly, $f$ is given by:

[TABLE]

The functions $D_{i}$ are given by:

[TABLE]

To establish a measure of force-freeness, we take the sum $\sum_{i}D_{i}$ and compare it to the (absolute) maximum $D_{i}$ term, and insist that the ratio fall below a given threshhold:

[TABLE]

In this work we insisted that $\epsilon<1\%$ over the entire domain. In practice our algorithm yields $\epsilon\ll 1\%$ over most of the domain with $\epsilon\approx 1\%$ only on some segments of the inner light surface.