Very High-Energy Emission from the Direct Vicinity of Rapidly Rotating Black Holes

When a black hole accretes plasmas at very low accretion rate, an advection-dominated accretion flow (ADAF) is formed. In an ADAF, relativistic electrons emit soft gamma-rays via Bremsstrahlung. Some MeV photons collide with each other to materialize as electron-positron pairs in the magnetosphere. Such pairs efficiently screen the electric field along the magnetic field lines, when the accretion rate is typically greater than 0.03-0.3% of the Eddington rate. However, when the accretion rate becomes smaller than this value, the number density of the created pairs becomes less than the rotationally induced Goldreich-Julian density. In such a charge-starved magnetosphere, an electric field arises along the magnetic field lines to accelerate charged leptons into ultra-relativistic energies, leading to an efficient TeV emission via an inverse-Compton (C) process, spending a portion of the extracted hole's rotational energy. In this review, we summarize the stationary lepton accelerator models in black hole magnetospheres. We apply the model to super-massive black holes and demonstrate that nearby low-luminosity active galactic nuclei are capable of emitting detectable gamma-rays between 0.1 and 30 TeV with the Cherenkov Telescope Array.


Introduction
It is commonly accepted that every active galaxy harbors a supermassive black hole (BH), whose mass typically ranges between 10 6 Mʘ and 10 9.5 Mʘ, in its center (e.g., [1][2][3][4]; see also [5,6] for reviews). Compelling evidence of a supermassive BH was found by observing the line emission from water masers around the central region of galaxy NGC 4258 [7]. Moreover, evidence of supermassive BHs at individual galactic centers are shown by the tight correlation between the black hole mass and the velocity dispersion or the bulge mass (e.g., [8][9][10]), which has been confirmed for the measurements with reverberation mapping (e.g., [11]).
A likely mechanism for powering such an active galactic nucleus (AGN) is the release of the gravitational energy of accreting plasmas [12] or the electromagnetic extraction of the rotational energy of a rotating supermassive BH [13]. The latter mechanism, which is called the Blandford-Znajek (BZ) mechanism, works only when there is a plasma accretion, because a BH cannot have its own magnetic moment (e.g., [14]). As long as the magnetic field energy is in a rough equipartition with the gravitational binding energy of the accreting plasmas, both mechanisms contribute comparably in terms of luminosity. The former mechanism is supposed to power the mildly relativistic winds that are launched from the accretion disks [15][16][17]. There is, however, growing evidence that relativistic jets are energized by the latter BZ mechanism through numerical simulations [18][19][20] (see also [21] for an ergospheric disk jet model). Indeed, general relativistic (GR) magnetohydrodynamic (MHD) models show the existence of collimated and magnetically dominated jets in the polar regions [22][23][24], whose structures are similar to those in the force-free models [25][26][27]. Since the centrifugal-force barrier prevents On these grounds, to achieve bright gap emissions, we began to consider spatially extended gaps, which can be possible if the soft photon field is weak enough so that the photon-photon pair-production mean-free path becomes comparable to or greater than the gravitational radius, rg = M. For example, when the mass accretion rate is typically less than 1% of the Eddington rate, the accreting plasmas form an ADAF, emitting radio to infrared photons via synchrotron process and MeV photons via free-free and IC processes ( § 4). Under such a low accretion environment, gap-emitted TeV photons do not efficiently collide with the soft photons, leading to an un-screened, spatially extended gap, which has a much greater electric potential drop, and hence gap L than the denser soft-photon-field case.
In this context, Neronov and Aharonian [54] examined a BH gap emission and applied it to the central engine of M87, a nearby low luminosity AGN, which hosts a BH with mass M ≈ (3.2-6.6) × 10 9 Mʘ [55][56][57][58]. They assumed that the gap is extended, w ≈ 2M (= Schwarzschild radius), that the lepton number density, N  , is comparable to the typical GJ number density, which is obtained when the magnetic buoyancy balances the disk gravity; Edd / m M M  refers to the dimensionless accretion rate near the horizon; and M denotes the mass accretion rate. We have 2 Edd Edd L M c   , where the radiation efficiency can be estimated as 0.1   . They also assumed that the magnetic-field-aligned electric field is comparable to the perpendicular component, where this denotes the distance from the rotation axis, with the spacetime dragging angular frequency due to BH's rotation, BH does the magnetic field strength evaluated at the horizon. In the right-most near equality, we evaluate quantities near the horizon, where a gap is formed. Because of these four assumptions (i.e., Boltzmann equations in the two-dimensional (2D) poloidal plane, assuming a mono-energetic approximation for the lepton distribution functions [62]. Subsequently, Hirotani et al. (2017) considered an inhomogeneous soft photon field to examine the γ-ray emission properties of supermassive BHs, solving the distribution functions of the accelerated leptons explicitly from their Boltzmann equations, discarding the mono-energetic approximation [63]. They showed that the accelerated leptons emit copious photons via IC processes between 0.1 and 30 TeV for a distant observer, and that these IC fluxes will be detectable with IACTs such as the Cherenkov Telescope Array (CTA), provided that a low-luminosity active galactic nucleus is located within 1 Mpc for a million-solar-mass central BH or within 30 Mpc for a billion-solar-mass central BH. Lin et al. then applied this method to stellar-mass BHs and compare the prediction with the high-energy (HE) observations, re-analyzing the archival Fermi/LAT data [64]. In addition, Song et al. demonstrated that the gap emission of an aligned rotator is enhanced along the rotation axis if the BH is nearly maximally rotating (i.e., a→M) [65], because the magnetic fluxes concentrate polewards as a→M due to the magnetic pinch effect [24,66]. More recently, Hirotani et al. (2018a,b) investigated if stellar-mass BHs can emit detectable HE and VHE emissions when they encounter dense molecular clouds (6) [67,68]. In these stationary analysis [61][62][63][64][65]67,68], They found 14 gap BZ (10 10 ) LL   , depending on the accretion rate. The BZ power becomes [13,23]:  (4) is obtained in the slow rotating limit, a«M to the second order, (a/M) 2 . For rapidly rotating BHs, the BZ power is obtained by Tanabe and Nagataki [69] to the fourth order, and by Tchekhovskoy et al. [24] to the sixth order. The higher-order corrections suppress the BZ power and a flattening occurs as a function of (a/M) 2 for extremely rotating BHs: for example, equation (4) over-predicts the BZ power when a > 0.95M for geometrically thin disks.

Force-Free Magnetosphere and the Necessity of a Gap
For the BZ process to work, there should exist a global electric current in the magnetosphere. However, plasmas flow inwards inside the inner light surface and flow outwards outside the outer light surface (e.g., [70]). Thus, plasmas should be continuously replenished between the two light surfaces (see §5.2 for details) so 2.2 10 mM   (7) the magnetosphere is kept force-free. In this case, there appears no gap. On the other hand, if the accretion rate stays below this value, the charge-starved magnetosphere inevitably has a gap, which may be either stationary or non-stationary.
Indeed, there is a growing consensus that the radio polarimetry in millimeter-submillimeter wavelengths provides useful diagnostics to infer the accretion rate of ADAF near the central engine. Through the Faraday rotation measure of the linear polarization, it is suggested that the dimensionless accretion rate becomes m =(1−10)×10 −6 within the radius r < 200 M (=100 Schwarzschild radii) for Sgr A* [71][72][73] and m < 10 −4 within r < 42 M for M87 [74]. In addition, the jet luminosity (~2 × 10 42 ergs s −1 ) of radio galaxy IC 310 [75] shows that the time-averaged accretion rate is m~10 −4 [76]. For low luminosity AGNs like M87 and IC 310 to sustain largescale jets under such small accretion rates, sufficient pair creation should be taking place by the other methods than the collisions of ADAF-emitted MeV photons. As demonstrated in previous BH gap models, BH gaps are indeed capable of supplying such plasmas in a charge-starved magnetosphere via the cascade of the gapemitted TeV photons into electron-positron pairs.

Detectability of Gap Emissions
As we have just seen, when the accretion rate becomes as small as:  (8)) as the thick solid line in Figure 1. Above this upper limit, the magnetosphere becomes force-free (blue shaded region). Thus, only under this thick solid line, BH gaps can be formed (white region). The dotted lines denote constant-LBZ lines as labeled.  (4)) for an extremely rotating case, a = M. From [62].
Second, to examine the upper bound of the gap luminosity, we substitute eq BB  into Equation (4) to obtain the BZ power, Putting a = M, we obtain the dotted lines in Figure 1.
Assuming that 100% of this power is converted into radiation, we obtain the upper limit of the gap flux at It follows that the IC component, which appears in VHE, may be detectable with ground-based IACTs. To quantify the actual flux as a function of the accretion rate, we must solve the gap electrodynamics. We will describe this issue in Sections 5-7.

Criticism of the Stationary Black Hole (BH) Gap Model
Levinson and Segev revisited the one-dimensional stationary BH gap model and found that the gap solutions are allowed only for small electric currents created within the gap [77]. They fully incorporated the GR effects in their basic equations, and solved the set of 1D Poisson equation, lepton Botzmann equations under mono-energetic approximation, and the radiative transfer equation, approximating all the IC photons gain the Galaxies 2018, 6, 122; doi:10.3390/galaxies6040122 7 of 37 initial lepton kinetic energy, assuming a homogeneous and isotropic soft photon field (emitted by an ADAF).
In their one-dimensional treatment, they found that E is proportional to 2 w , and the potential drop, gap V , is proportional to 3 w (because the Poisson equation is a second-order differential equation), as Figure 2a indicates. They imposed a boundary condition such that the created charge density matches the local GJ charge density at the outer boundary; accordingly, cr On the other hand, in the present review, we introduce the electric current per magnetic flux (5.4),    [78]). In this section, we present the result of a pulsar case, whose electrodynamics is essentially the same as BHs. In Figure 3, we show the pulsar's E solved for three discrete gap currents, cr 0.01 j  (dashed curve), 0.03 (solid curve), and 0.0592 (dotted curve). As the gap width approaches D  = 3 × 10 8 cm, E is substantially suppressed by the 2D screening effect. As a result, pulsar OG solutions exist in a wide range of cr j , namely 0 < cr j < 1. In the same way, for BHs, a stationary gap solution exists in a wide, 0 < cr j < 1, by virtue of the 2D screening effect, as will be demonstrated in §6.5. In another word, the difficulty of stationary BH gap model pointed out by [77], can be solved if we consider the 2D electrodynamic structure.

Background Space-Time Geometry
To examine the gap electrodynamics further quantitatively, let us describe the space time around a rotating BH. Since the self-gravity of the accreting plasmas and the electromagnetic field is negligible compared to the BH's gravitational field, the spacetime around a rotating BH is well described by the Kerr metric [79]. In the Boyer-Lindquist coordinate [80], it becomes: where, which does not depend on  , as expected. To lower an index of a tensor, we can use the following relations, vanishes at the horizon. We must choose appropriate vector and tensor components that are well-behaved at the horizon, as briefly discussed below in Equation (30).

Non-Stationary BH Gap Model
The aforementioned BH gap scenarios, however, are based on the stationary assumption. Thus, Levinson and Cerutti [81] examined the lepton acceleration and gamma-ray emission within a gap around a rotating BH, performing one-dimensional particle-in-cell (PIC) simulations. From the homogeneous part of the Maxwell equations (i.e., the Faraday law), one obtains: where the electro-magnetic field strength tensor is related to the scalar and vector potentials by . They then investigated one-dimensional disturbances along the radial direction, setting 0      . Thus, they obtained: It follows from 0 t F   that the radial component of the magnetic field should be time-independent.
As an example, they assumed a split-monopole solution, fixing the magnetic flux function to be ( , ) ( where J  refers to the electric four current. Putting r   and, and noting that and, holds, we obtain: and,  (24), (29) and (30). The electric current density J  in Ampere's law is constructed from the actual motion of the charged particles in the PIC scheme.
By this one-dimensional, general-relativistic (GR) PIC simulation scheme, they demonstrated that E is efficiently screened out by the discharge of created pairs within the duration that is comparable to the dynamical time scale, 3 g // r c GM c M  , that the residual E leads to a strong flare of VHE curvature radiation, which is followed by a state of self-sustained, rapid plasma oscillations, and that the resultant quasistationary gamma-ray luminosity attains only about 10 −5 of the BZ power; that is, the time-averaged luminosity becomes less than the stationary analysis. See also Levinson et al., who semi-analytically found longitudinal oscillations in their two beam model [82]. Figure 4 exhibits the snapshots at four discrete simulation times. The upper panel shows the variation of the pair and photon densities. Pair density increases due to pair production from zero and saturates at about 10 times GJ value within t = 5 GM/c 3 = 5 M, that is, within a few dynamical time scales, where we may define the dynamical time scale to be 2 GM/c 3 = 2 M. Accordingly, E is nearly screened out. The sporadic pair creation leads to rapid spatial and temporal oscillations of E (middle panel) and r J (bottom panel). Above nearly 10 dynamical time scales, t > 20 M, pair and photon distribution function attain quasi-stationary state. Figure 5 exhibits the light curve evolution. It follows that the gamma-ray luminosity approaches the BZ power during the initial spike but decays to the terminal value as E is screened out. It is possible to argue that strong, rapid flares should be produced every time a magnetospheric gap is restored. The flaring episode is then followed by quiescent emission with a luminosity,  More recently, Chen et al. performed 1D PIC simulations on the BH gap model [83]. They assumed a Minkowski spacetime in the Maxwell equations, lepton equation of motion, and radiative transfer equation, except for the GJ charge density, whose GR expression is essential for the formation of a null surface near the event horizon, whereas full GR expressions were implemented in all the basic equations in [81]. Cheng et al. adopted a Newtorian split-monopole solution [84], and solved the time evolution of r EE  , E  , and B  . It is found that a BH gap opens quasi-periodically and the screening of E creates oscillations. Eventually, when the created pairs escape across the light surfaces, the charge-starved magnetosphere recur in the formation of a gap. Figure 6 shows the third cycle of such a gap development. The middle panel in the top indicates that E arises around the null surface, which is shown in the third panels as the crossing of the GR charge density (green curve in the orange band) and 0 (horizontal green line). However, the gap can also appear at different positions, as the right two panels in the top row show. The gap (shown in the top row) appears in the low multiplicity region, as indicated in the second row. Note that the green curve in the bottom panels show the spectra of the photons that materialize within the simulation box; thus, they do not show the gamma-ray spectra to be observed.  In short, time-dependent, one-dimensional GR PIC simulations revealed that gaps are unsteady but tend to a quasi-steady state or a quasi-periodical state. The former, quasi-steady solution shows that the timeaveraged gap gamma-ray luminosity is of an order of magnitude smaller than stationary analysis because of the efficient screening [81]. Nevertheless, in the present review we focus on the stationary BH gap scenario in order to elucidate the detailed gap electrodynamics, which will be also important to understand the timedependent phenomena.

Advection-Dominated Accretion Flow
Before moving onto the gap electrodynamics, let us briefly describe the radiatively inefficient accretion flow in the lower latitudes (i.e., near the equator), because it is essential to quantify the photon-photon pair production and IC scatterings taking place inside and outside the gap in the higher latitudes (i.e., in the polar funnel). Since we are concerned with the soft photon field near the BH, we do not consider the structure of accretion flows away from the BH, focusing on an ADAF, a representative model of radiatively inefficient accretion flow. We consider that plasma accretion takes place near the rotational equator with a certain vertical thickness (e.g., [20,85,86]). Such an accretion cannot penetrate into the higher latitudes, that is, in the polar funnel because the centrifugal-force barrier prevents plasma accretion toward the rotation axis, and because the timescale for magnetic Rayleigh-Taylor instability or turbulent diffusion is long compared to the dynamical timescale of accretion. In this evacuated funnel, the poloidal magnetic field lines resemble a split monopole in a time-averaged sense [22]. The lower-latitude accretion emits MeV photons into the higher latitudes, supplying electron-positron pairs in the funnel. If the pair density in the funnel becomes less than NGJ, E arises. Thus, a plasma accretion in the lower altitudes and a gap formation in the higher altitudes are compatible in a BH magnetosphere.
Since the gap luminosity increases with decreasing accretion rate [61][62][63]67,68], we consider a dimensionless accretion rate that is much small compared to unity. Such a low accretion status corresponds to the quiescence or low-hard state for X-ray binaries, and low luminosity active galactic nuclei for super-massive BHs. At such a low accretion rate, the equatorial accretion flow becomes optically thin for Bremsstrahlung absorption and radiatively inefficient because of the weak Coulomb interaction between the ions and electrons. Accordingly, the ions store the heat within the flow itself and accretes onto the BH without losing the thermal energy as radiation. This radiatively inefficient flow can be described by an ADAF [28,29,[87][88][89][90][91][92], and provides the target soft photons for the IC-scattering and the photon-absorption processes in the polar funnel. Thus, to tabulate the redistribution functions for these two processes, we compute the specific intensity of the ADAFemitted photons. For this purpose, we adopt the analytical self-similar ADAF spectrum presented in [90]. The spectrum includes the contribution of the synchrotron, IC, and Bremsstrahlung processes. These three cooling mechanisms balance with the heating due to the viscosity and the energy transport form ions, and determine the temperature of the electrons in an ADAF to be around Te~10 9 K. In radio wavelength, the ADAF radiation field is dominated by the synchrotron component whose peak frequency, c,syn  , varies with the accretion rate , the seed pair density becomes less than the GJ value [59], thereby leading to an occurrence of a vacuum gap in the funnel.
It is noteworthy, however, that at low accretion rate, namely 3 10 m   , less efficient thermalization may lead to a hybrid thermal-nonthermal energy distribution. Observationally, it is indeed known that non-thermal electrons are needed to reproduce the quiescent low-frequency radio emission in Sgr A* [93][94][95] and other lowluminosity AGNs [96]. Theoretically, the direct heating of electrons (i.e., not via Coulomb collisions with ions) is proposed via several mechanisms, such as magnetic reconnection [97][98][99][100][101], MHD turbulence [98,[102][103][104][105], or dissipation of pressure anisotropy in a collisionless plasma [106]. At the present stage, there is no consensus on the theory of electron heating yet. In the present review, we neglect such direct heating of electrons and assume that only less than 1% of the ions thermal energy is transferred to electrons. However, if more gravitational energy is converted into electrons' heat by direct heating (i.e., if the fraction of electron heating becomes » 0.01), the increased ADAF soft photon field would significantly reduce the gap luminosity at a fixed accretion rate. We should keep this point in mind when we construct a BH gap model, irrespective of whether the gap is stationary or non-stationary. We will present the input ADAF spectrum when we show the predicted gammaray spectrum of the gap emissions.

Stationary, Two-Dimensional Gap Electrodynamics
Next, we formulate the BH gap, which will be formed when up mm  .

Poisson Equation for the Non-Corotational Potential
In a stationary and axisymmetric spacetime, Gauss's law for the electrostatic potential gives where  denotes the covariant derivative, In the present formalism, we assume that the electromagnetic fields depend on t and  only through  We adopt the static observer whose four velocity is given by the Killing vector, The acceleration electric field is given by Substituting Equation (30) into (29), we obtain the Poisson equation for the non-corotational potential, where the general-relativistic Goldreich-Julian (GJ) charge density is defined by Not only F  and  but also  may depend on t and  through F t   . In this case, equation (34) gives "stationary" gap solution in the "co-rotational" frame, in the sense that  and  are a function of r,  , and

The Null-Charge Surface
Equation (34) shows that In a stationary gap, the derivative of E along B should change sign at the outer and the inner boundaries, where the outer boundary is located at larger distance from the hole, while the inner boundary is located near the hole. Thus, a vacuum gap (with 0   ) is located around the null-charge surface, on which GJ  changes sign. If the gap becomes non-vacuum (i.e.,  is slight smaller than GJ  ), the gap position is essentially unchanged (6.6). Thus, the null-charge surface becomes the natural place for a gap to arise, in the same way as the pulsar outer-gap model (for an analytic argument of the gap position, see 2 of [112]).
In the case of BHs, the null-charge surface appears near the surface on which  coincides with F  [52].
Since  matches F  only near the horizon, the gap inevitably appears near the horizon, irrespective of the BH mass. In Figure 7, we plot the distribution of the null surface as the thick red solid curve. In the left panel, we assume that the magnetic field is split-monopole, adopting (1 cos ) AC    as the magnetic flux function [13], where C is a constant. In the right panel, on the other hand, we assume a parabolic magnetic field line with    [20,113].
It follows that the null surface distributes nearly spherically, irrespective of the poloidal magnetic field configuration. This is because its position is essentially determined by the condition because ω has weak dependence on  , and because we assume F  is constant for A  for simplicity. If F  decreases (or increases) toward the polar region, the null surface shape becomes prolate (or oblate). The field angular frequency F  is, indeed, deeply related to the accretion conditions. For instance, in order to get adequate jet efficiency from the magnetic field in the funnel above the BH for a radio-loud active galactic nuclei (AGNs), one needs significant flux compression from the lateral boundary imposed by the accretion flow and corona [114]. Although it is not a self-consistent solution, we adopt a constant F  in the present review for simplicity. Let us quickly take a look at the case of slower and faster rotations. Figure 8 shows the GJ charge density and the null-charge surface distribution when the poloidal magnetic field line is radial (i.e., split-monopole); the BH mass and spin are the same as Figure 7.

The Stagnation Surface
In a stationary and axisymmetric BH magnetosphere, an MHD fluid flows from a greater k0 region to a smaller k0 one. It follows that both the inflows and outflows start from the two-dimensional surface on which k0 maximizes along the poloidal magnetic field line. Thus, putting k0' = 0, we obtain the position of the stagnation surface (thick green solid curve in Figure 9), where the prime denotes the derivative along the poloidal magnetic field line. The condition k0' = 0 is equivalent to imposing a balance among the gravitational, centrifugal, and Lorentz forces on the poloidal plane. In Figure 9, we plot the contours of k0 for split-monopole and parabolic poloidal magnetic field lines.
It is clear that the stagnation surface is located at r < 5 M in the lower latitudes but at r > 5 M in the higher latitudes. As the poloidal field geometry changes from radial to parabolic, the stagnation surface moderately moves away from the rotation axis in the higher latitudes. This analytical argument of the position of the stagnation surface [115], was confirmed by general relativistic (GR) MHD simulations [60]. In these numerical works, the stagnation surface is time-dependent but stably located at 5-10 M with a prolate shape, as depicted in Figure 9. The two-dimensional surfaces on which k0 vanishes, are called the outer and inner light surfaces (thick red solid curves in Figure 9). Outside the outer light surface (or inside the inner light surface), plasma particles must flow outwards (or inwards).

The Gap Position
To examine the Poisson equation (34), it is worth noting that the distribution of GJ  has relatively small dependence on the magnetic field geometry near the horizon, because its distribution is essentially governed by the space-time dragging effect, as indicated by the red solid, dashed, and dotted contours in Figure 7 for split-monopole and parabolic fields. As a result, the gap solution little depends on the magnetic field line configuration, particularly in the higher latitudes. In what follows, we thus assume a split-monopole field and adopt a constant F  for A  for simplicity, as described in the left panels of Figures 7 and 9.
In a pulsar magnetosphere, the gap position is located around the null-charge surface if there is no leptonic injection across either of the inner and outer boundaries, and that the gap position shifts outwards (or inwards) when leptons are injected into the gap across the inner (or the outer) boundary [116,117]. Since the gap electrodynamics is essentially the same between pulsar outer gaps and BH gaps, we can expect that the same conclusion holds in the case of BH magnetospheres.
In a BH magnetosphere, a gap indeed distributes around the null-charge surface, as will be shown in sections 6.3 and 7.2.

Advection-Dominated Accretion Flow (ADAF)-Emitted Soft Photon Field in the Magnetosphere
To quantify the gap electrodynamics, we need to compute the pair creation rate. To this end, we must tabulate the specific intensity of the soft photons at each position in the polar funnel. In this section, we therefore consider how the soft photons are emitted in an ADAF and propagate around a rotating BH. We assume that the soft-photon field is axisymmetric with respect to the BH rotation axis. When the mass accretion rate is much smaller than the Eddington rate, the accreting plasmas form an ADAF with a certain thickness in the equatorial region. For simplicity, in this paper we approximate that such plasmas rotate around the BH with the GR Keplerian angular velocity, and that their motion is dominated by this rotation. That is, we neglect the motion of the soft-photon-emitting plasmas on the poloidal plane, (r, ), for simplicity. Let us introduce the local rest frame (LRF) of such rotating plasmas. The orthonormal tetrad of the LRF is given by [63], where, and the redshift factor between the LRF and the distant static observer becomes: Note that the GR Keplerian angular frequency on the equator, K  , is different from the magnetic-field angular frequency, F  . Thus, D is not k0. Note also that equation (42) is the reciprocal of what appears in equation (38).
For further details on how to compute the specific intensity of the ADAF soft photon field, see [63]. Between the distant static observer (i.e., us) and the LRF, the photon energy changes by the redshift factor,   for a fixed ( t k  , k  ), which enables us to ray-trace the soft photons in the Kerr space time [63].
In what follows, we compute the collision frequencies of two-photon pair creation and IC scatterings in the frame of a zero-angular-momentum observer (ZAMO), which rotates with the same angular frequency as the space-time dragging frequency, / t gg     . Putting    , we obtain the lapse (section 3.3 of [63]): Between the distant static observer and ZAMO, the photon energy changes by the redshift factor, Combining equations (46) and (58), we find that the photon energy changes from the LRF (in the disk) to ZAMO (in the magnetosphere including the polar funnel) by the factor: By the ray-tracing technique described above (see also [118]), we can now compute the specific intensity of the soft photons emitted from the ADAF at each place in the magnetosphere in the ZAMO frame of reference. Integrating this ZAMO-measured specific intensity over the propagation solid angle in each directional bin at each point in ZAMO, we obtain the soft photon flux at each point in each direction. Finally, we use this flux to compute the IC optical depth and the photon-photon collision optical depth in ZAMO.
Let us define s s 0 ( / ) dF dE to be the differential soft photon flux that would be obtained if the specific intensity were homogeneous and isotropic in a Minkowski space time within the radius 6M, and if the total soft photon luminosity were given by radial flux at 6M multiplied by

The Poisson Equation in the Tortoise Coordinate
To solve the radial dependence of Φ in the Poisson equation (34), we introduce the following dimensionless tortoise coordinate, In this coordinate, the horizon corresponds to the "inward infinity," *    . In this paper, we set whereas the dimensionless GJ number density is defined by: N± denotes the distribution function of positrons and electrons in ordinary definition. We have introduced the dimensionless non-corotational potential, Without loss of any generality, we can assume 0 F   in the northern hemisphere. In this case, a negative E‖ arises in the gap, which is consistent with the direction of the global current flow pattern.

Boltzmann Equations for Electrons and Positrons
Next, let us describe the lepton Boltzmann equations. Imposing a stationary condition, Where  denotes the pitch angle (with 0   for out-going leptons and   for in-coming ones), s and p show the position and dimensionless momentum along the magnetic field, SIC denotes the rate of positrons or electrons transferred into Lorentz factor  from another Lorentz factor via IC scatterings, Sp refers to that which appeared into Lorentz factor  via photon-photon and magnetic pair production, and  refers to the lapse function, or equivalently, the gravitational redshift factor. In the Boyer-Lindquist coordinate, we obtain It follows that 0   at the horizon and 1   at infinity. We consider both photon-photon and magnetic pair production process; however, only the former contributes in all the cases we consider. The upper and lower sign of n, SIC, and Sp corresponds to the positrons (with charge q = +e) and electrons (q = −e), respectively. The dimensionless momentum is related to the Lorentz factor by where the curvature radiation force is given by (e.g., [119]), . The IC and pairproduction redistribution functions are given in [50,63]. It is noteworthy that the charge conservation ensures that the dimensionless total current density (per magnetic flux tube),

Radiative Transfer Equation
We assume that all photons are emitted with vanishing angular momenta and hence propagate on a constant- cone. Under this assumption of radial propagation, we obtain the radiative transfer equation: where rr dl g dr  refers to the distance interval along the ray in ZAMO,   and j  the absorption and emission coefficients evaluated in ZAMO, respectively. We consider only photon-photon collisions for absorption, pure curvature and IC processes for primary lepton emissions, and synchrotron and IC processes for the emissions by secondary and higher-generation pairs. For more details of the computation of absorption and emission, see Section 4.2 and 4.3 of [61] and Section 5.1.5 of [62].

Boundary Conditions
The elliptic type second-order partial differential Equation (49) (54), we need to specify the particle injection rate across the inner and outer boundaries. Since the magnetospheric current is to be constrained by a global condition including the distant dissipative region, we should treat cr j , in j and out j as free parameters, when we focus on the local gap electrodynamics. For simplicity, we assume that there is no electron injection across the inner boundary and put in 0 j  throughout this paper. In what follows, we examine stationary gap solutions for several representative values of cr j and out j . The radiative-transfer Equation (58), a first-order ordinary differential, contains no photon injection across neither the outer nor the inner boundaries.

Gap Closure Condition
We compute the multiplicity (Eq.

Lepton Accelerator around Stellar-Mass Black Holes
Before we come to the main subject, the case of super-massive BHs, it is desirable to investigate stellarmass BHs, because important gap electrodynamics have so far been elucidated in this case [67,68].

Black Hole Accretion in a Gaseous Cloud
When a BH moves in an ambient medium, accretion takes place onto the BH. Since the temperature is very low in a molecular cloud, the BH's velocity V becomes supersonic with respect to the gaseous cloud. For a uniform medium, particles within the impact parameter for a molecular hydrogen gas, where  Figure 11a. Since the accreting gases have little angular momentum as a whole with respect to the BH, they form an accretion disk only within a radius that is much less than the Bondi radius B r . Thus, we neglect the mass loss as a disk wind between B r and the inner-most region, and evaluate the accretion rate near the BH, m with B m . In section 6, we consider a 10-solar-mass BH, which is typical as a stellar-mass BH [122,123].

Lepton Accelerator around Stellar-Mass Black Holes
In a vacuum, rotating magnetosphere, an acceleration electric field, E , inevitably arises. Accordingly, leptons (red arrows in Figure 11b) are accelerated into ultra-relativistic energies to emit high-energy gammarays (wavy line with middle wavelength) via the curvature process and VHE gamma-rays (wavy line with shortest wavelength) via the inverse-Compton (IC) scatterings of the soft photons (wavy line with longest wavelength) that were emitted from the ADAF. A fraction of such VHE photons collide with the ADAF soft photons to materialize as e ± pairs, which partially screen the original E‖ when they polarize. To compute the actual strength of E , we solve the e ± pair production cascade in a stationary and axisymmetric magnetosphere on the meridional plane (r, ). The magnetic field lines are twisted into the counter-rotational direction due to the poloidal currents, and the curvature radius of such a toroidal field is assumed to be g rM  in the local reference frame.

Results: Acceleration Electric Field
To estimate the maximum luminosity of a gap, we consider an extremely rotating case of a = 0.99 ( Figure  4 of [68]). It follows from Figure 12a that E peaks along the rotation axis, because the magnetic fluxes concentrate towards the rotation axis as a→M. In what follows, we thus focus on the emission along the rotation axis, 0   . E decreases slowly outside the null surface in the same way as pulsar outer gaps (e.g., Figure 12 of [124]). This is due to the two-dimensional screening effect of E .  [67,68].
In Figure 12b, we present ( , 0 ) Es   at four dimensionless accretion rates (as indicated in the box).
We find that the potential drop increases with decreasing m . However, if the accretion further decreases as m < 10 −4.25 , there exists no stationary gap solutions. Below this lower bound accretion rate, the gap becomes inevitably non-stationary (section 8.1).

Results: Gamma-Ray Spectra
Let us examine the gamma-ray spectra. In Figure 13a, we present the spectral energy distribution (SED) of the gap emission along five discrete viewing angles. It follows that the gap luminosity maximizes if we observe the BH almost face-on,  < 15°, and that the gap luminosity rapidly decreases at  < 30°. However, if the BH is moderately rotating as a = 0.50 M, it is very difficult to detect its emission, unless it is located within 0.3 kpc [68]. It also follows that the spectrum has two peaks. The curvature photons are emitted by the gap-accelerated primary electrons and appear between 5 MeV and 0.5 GeV. The same electrons up-scatter the ambient soft photons, forming the second peak between 0.5 GeV and 5 GeV. The latter IC photons suffers absorption above 0.1 TeV ( Figure 5 of [68]).

Results: Dependence on the Current Created in the Gap
Let us demonstrate that the gap solution exists in a wide range of the created electric current within the gap, cr cr GJ / j J J  , and that the resultant gamma-ray spectrum little depends on cr j . In Figure 14, we show the solved () Es (left panels) and SEDs (right panels) for three discrete cr j 's: from the top, they correspond to cr j = 0.3, 0.5, and 0.9. The case of cr j = 0.7 is presented as Figure 13b. It is clear that the gap spectra modestly depend on the created current within the gap as long as the created current is sub-GJ. There exist no stationary gap solutions if cr 1 j  . Note that the maximum value of E saturates below 4.5 × 10 5 statvolt cm −1 because of the 2D screening effect, which becomes important when the gap width become comparable to or greater than the transverse thickness ~M (3.4). It can also be shown that the SEDs depend little on the injected current across the outer boundary ( Figure  7 of [68]). For further details, including the authors' comments to a criticism raised by [77] regarding the locations of stationary gaps and the stagnation surface, see section 6.2 of [68]. the GJ value. In all the three cases, the injected current density across the inner or outer boundaries is set to be zero. The curves corresponds to the accretion rate of m =3.16×10 −4 (red dotted), 1.77 × 10 −4 (blue dashed), 1.00 × 10 −4 (black solid), and 5.62 × 10 −5 (green dash-dotted). From [68].

Lepton Accelerator around Super-Massive Black Holes
In this section, we apply the BH gap model to supermassive BHs.

Very High-Energy Observations of Active Galaxies
The High-Energy Stereoscopic System (H.E.S.S.), the Major Atmospheric Gamma Imaging Cherenkov (MAGIC), and the Very Energetic Radiation Imaging Telescope Array System (VERITAS) have so far detected 219 TeV sources [125]. Among them, 70 sources have been identified as AGNs. Although most of them are blazer type, NGC 1275, 3C264, M87, and Cen A are classified as Fanaroff-Riley (FR) I radio galaxies. In addition, another radio galaxy, IC 310 appears to be of a transitional type between FR I and BL Lac. For these five nonblazer radio galaxies, their inner jet (on parsec or sub-parsec scales) are moderately misaligned with respect to the line of sight. Thus, their TeV emission from the jets will not be highly Doppler-boosted. Accordingly, the emission from the direct vicinity of the central supermassive BH, if it exists, may not be hidden by the strong jet emissioin. We are, therefore, motivated to study these non-blazer radio galaxies, comparing theoretical predictions with VHE observations. See Rieger and Levinson [126] for a comprehensive review on the VHE observations of these non-blazer radio galaxies, and associated theoretical arguments.
From NGC 1275, a TeV flare has been reported by MAGIC [127] with the shortest variation time scale of var t ~10 h. Since the light crossing time of the event horizon is about 3 × 10 3 s for this source, the dimension of the TeV emitting region should be less than 10 times Schwarzschild radii. M87 exhibited rapidly varying VHE flares that are observed in 2008 and 2010 [128][129][130][131], where the shortest variation time scale was observed by all the three IACTs, and was found to be var t ~1 day [128]. Its horizoncrossing time scale, 4 × 10 4 s, indicates that the TeV emitting region should be less than a few Schwarzschild radii. Also, VHE photons were detected during the other periods, including the low-emission state [132][133][134].
Results of multi-wavelength observations of M87 were summarized in [135].
IC 310 showed very rapid variabilities whose shortest time scale attained var t ~5 m = 300 s, which was detected by MAGIC [77]. Compared with its horizon-crossing time scale of 3 × 10 3 s, we find that the TeV photons should be emitted from a sub-horizon length scale, presumably from the direct vicinity of the event horizon. For a comparison of representative emission models for such rapidly varying extragalactic sources in the TeV sky, see [136].
In the rest of this review, we consider the BH gap model as a possible explanation of such horizon-scale or sub-horizon-scale TeV flares observed from M87 and IC 310.

Results: Acceleration Electric Field
To investigate the VHE emission from the direct vicinity of a supermassive BH, we apply the BH gap model, and comapre the results obtained for M = 10 9 Mʘ. and 10 6 Mʘ. Unless explicitly mentioned, we adopt is adopted. Near the lower-latitude boundary,  =60°, a small E extends into higher altitudes, because a smaller E requires a greater width for the gap closure condition to be satisfied. This behavior is common with pulsar lower-altitude slot-gap models [137,138]. However, near the rotation axis,  = 0°, stronger E arises, mainly because the polar region is far from the meridional boundary at  = 60°, and partly because the magentic fluxes more or less concentrate toward the rotation axis even for a non-extreme rotator with a = 0.

Results: Gap Width
Let us briefly examine how the gap width, w, depends on the ADAF soft-photon field as the accretion rate changes. In Figure 16

Results: Gamma-Ray Spectra
The predicted photon spectra are depicted in Figure 18 for six m values. The thin curves on the left denote the input ADAF spectra, while the thick lines on the right denote the output spectra from the gap. The left panel shows the result in the case of a billion solar-mass BH assuming the distance of d = 10 Mpc, while the right panel does that in the case of a million solar-mass BH with d = 1 Mpc. In both cases, the emitted γ-ray flux increases with decreasing m , because the decreased soft photon field increases the photon-photon collision mean free path, and hence the gap width.
In the case of M = 10 9 Mʘ (left panel), the spectra peak between 1 and 30 TeV, because the IC process dominates the curvature one. It is clear that the gap HE flux lies well below the detection limit of the Fermi/LAT (three thin solid curves labeled with "LAT 10 year" [139]). Nevertheless, its VHE flux appears above the CTA detection limits (dashed and dotted curves labeled with "CTA 50 h" [140]). Is is possible to argue that such a large VHE flux will be detected within one night with CTA from a nearby low-luminosity like M87, when the accretion rate coincidentally resides in the range, 6 × 10 −6 < m < 3 × 10 −5 .
In the case of M = 10 6 Mʘ (right panel), the spectra peak between 30 GeV and 10 TeV. The IC process dominates the curvature one also for million solar-mass BHs. When the accretion rate is in the narrow range 2×10 −5 < m < 3×10 −5 , the gap emission may be detectable with CTA, particularly if the source is located in the southern sky.

The Case of Very Small Accretion Rate
Let us discuss the case of very small accretion rates, m < low ( , ) m M a . The ADAF photon field peaks around eV for stellar-mass BHs and around meV for supermassive BHs. These photons are emitted from the innermost region, r∼Rmin∼6M, which is located well inside the gap outer boundary when low mm  ; thus, pair production is sustained only marginally in an extended gap. However, at low mm  , stationary pair production can be no longer sustained and a vacuum region develops in the entire polar funnel. In this vacuum region, migratory leptons are accelerated by the vacuum E and cascade into copious primary electrons and positrons that are accelerated in the opposite directions. The resultant emission will become inevitably non-stationary.

Gap Emission Versus Jet Emission
Next, let us discuss how to discriminate the gap and jet emissions. As we have seen in Sections 6 and 7, the gap gamma-ray luminosity increases with decreasing m . Thus, we can predict an anti-correlation between the ADAF-emitted infrared (IR)/optical and the gap-emitted HE/VHE fluxes for stellar-mass BHs, and between the ADAF-emitted radio and the gap-emitted VHE fluxes for supermassive BHs. This forms a contrast to the standard shock-in-jet scenario, whose natural prediction will be a correlation between the radio/IR/optical and the HE/VHE fluxes (for reviews or catalogues of gamma-ray observations of AGNs, see e.g., [126,[141][142][143][144][145]). This anti-correlation or correlation can be detectable both for stellar-mass and supermassive BHs. For stellar-mass BHs, BH transients may exhibit a correlation between IR/optical and HE/VHE fluxes during quiescence or lowhard state. For supermassive BHs, nearby low-luminosity AGNs may show an anti-correlation between submillimeter wavelength and VHE. In X-rays, the gap emission is very weak. Thus, if X-ray photons are detected, they are probably emitted from the jet or from the accretion flow.

Comparison with Pulsar Gap Models
Let us compare the present BH gap model with the pulsar outer (OG) gap model. In both gap models, the gap appears around the null-charge surface where the GJ charge density vanishes, as the stationary solution of the Maxwell-Boltzmann equations. Moreover, the gap width, and hence the γ-ray efficiency, which is defined by the ratio between the γ-ray and spin-down luminosities, increases with decreasing soft-photon flux. This is because the pair-production mean-free path increases with decreasing photon flux. For instance, the OG γ-ray efficiency increases with decreasing NS surface emission, or equivalently with increasing NS age [45,146]. In the same manner, the BH gap γ-ray efficiency increases with decreasing accretion rate of ADAF. In addition, electron-positron pairs are supplied by photon-photon pair production in both pulsar and BH magnetosphere, although ions can be drawn into the pulsar OG from the neutron-star surface as a space-charge limit flow [147], in the same manner as in PC models (2).
However, there are differences, as described below. First, in a pulsar magnetosphere, the null surface appears because of the convex geometry of a dipole magnetic field. On the other hand, in a BH magnetosphere, a null surface appears due to the frame-dragging effect around a rotating BH. Thus, for the same magnetic field polarity, the sign of E is opposite.
Second, in the OG model, a soft-photon field is provided by the cooling NS thermal emission and/or the heated polar-cap thermal emission. The NS surface emissions peak in the X-ray energies. Thus, the curvature 1-10 GeV photons efficiently materialize as pairs within the gap, thereby contributing to the gap closure. However, the thermal IR photon field is much weaker than the thermal X-ray field; thus, the IC photons do not contribute to the gap closure in pulsars. On the other hand, around a BH it is the ADAF that provides the softphoton field. For a supermassive BH, the ADAF spectrum peaks in submillimeter wavelengths. Thus, the photons having energies around 100 TeV or above efficiently materialize as pairs to close the gap. The 0.1-10 GeV curvature photons, on the other hand, do not materialize efficiently because the IR photon number flux is five to six orders of magnitude weaker than the submillimeter one.
Let us briefly compare the BH gap model with the pulsar polar-cap (PC) model. In a pulsar magnetosphere, electrons may be drawn outward as a space-charge-limited flow at the NS surface in the PC region. Thus, in a stationary gap [31,32] or in a non-stationary gap [34], γ-ray emission could be realized without pair creation within the polar cap, although pairs are indeed created via magnetic pair creation (e.g., at least at the outer boundary where E‖ should be screened). However, in BH magnetospheres, causality prevents any plasma emission across the horizon. Thus, a gap could not be sustained without pair creation.
If we compare the stationary BH gap model with pulsar OG and PC models, it is possible to argue for the stability of gaps in a qualitative manner ( §6.3 of [68]). In this review, finally we discuss the fact that the nonstationary nature of pulsar PCs [34] cannot readily be extended to the stability arguments of stationary BH gap solutions. In a pulsar PC, there exists no null surface. Thus, if a three-dimensional PC region is charge-starved in the sense simulations [34]. In short, pulsar PC accelerators are inherently unstable for pair production because of 0 E  , which stems from the fact that there exists no null-charge surface in the PC region due to the relatively weak frame-dragging effect (i.e., the Lense-Thirring effect [148]) around a rotating neutron star.
On these grounds, we cannot readily conclude that the stationary BH gap solutions are unstable because of the highly time-dependent nature of the pulsar PC accelerators. Careful examinations with PIC simulations are needed for BH gaps, as performed recently by [81,83]. Since the saturated solution in recent PIC simulation is much less violently time-dependent compared to pulsar PC cases [34,81], it may be reasonable to start with stationary BH gap models as the first step.
Funding: This research received no external funding.