High-energy and Very High Energy Emission from Stellar-mass Black Holes Moving in Gaseous Clouds

Molecular clouds are generally gravitationally bound and occasionally contain several sites of star formation. Particularly, massive stars formed in a GMC can ionize the surrounding interstellar medium with their strong UV radiation. A combined action of such ionization, stellar winds, and supernova explosions blows off the gases in a GMC, leaving an OB association adjacent to dense molecular clouds. In this section, we focus on the physical parameters in such dense molecular clouds, which may be traversed by a BH formed in a neighboring OB association.

The physical parameters (e.g., temperature and density) of a molecular cloud can be examined by observing the strength, width, and profile of radio emission lines of probe molecules. A typical GMC has gas kinematic temperature between 30 and 50 K. Using the dense gas tracers, H₂CO and CS, we can infer the hydrogen molecule number density, ${n}_{{{\rm{H}}}_{2}}$, in the core of a GMC. A typical GMC core has density ${10}^{4}\,{\mathrm{cm}}^{-3}\lt {n}_{{{\rm{H}}}_{2}}\,\lt {10}^{6}\,{\mathrm{cm}}^{-3}$ and mass between 10 and 10³ M_⊙. Individual molecular clouds have lower densities, ${10}^{2}\,{\mathrm{cm}}^{-3}\lt {n}_{{{\rm{H}}}_{2}}\,\lt {10}^{5}\,{\mathrm{cm}}^{-3}$, and masses between 10³ and 10⁶ M_⊙. We use these values of ${n}_{{{\rm{H}}}_{2}}$ to estimate the mass accretion rate onto a BH in Section 3.3.

To consider the passage of a stellar-mass BH in a dense molecular cloud, let us briefly comment on the massive-star formation in a GMC. In a GMC, massive stars are formed in OB associations. A typical OB association contains 10¹–10² high-mass stars of type O and B and 10²–10³ stars of lower masses. The stellar line-of-sight velocity dispersion in an OB association is typically around 9 km s⁻¹ and usually less than 20 km s⁻¹ (Sitnik 2003). The strong winds and the supernovae resulting from these associations blow off the interstellar medium; thus, OB associations are found adjacent to molecular clouds (Blitz 1980). Depending on the mass and metallicity, such high-mass OB stars evolve into neutrons stars or BHs after core-collapse events. For example, if the progenitor has mass in the range of 25 M_⊙ < M < 40 M_⊙ with low or solar metallicity, it will evolve into a BH through a supernova explosion after the fallback of material onto an initial neutron star. In this case, the BH will acquire a certain kick velocity with respect to the star-forming region in a similar way to neutron stars. On the other hand, if it has M > 40 M_⊙ with low metallicity, it will evolve into a BH directly without a supernova explosion. These massive BHs will have smaller relative velocities with respect to the star-forming region compared to the lighter BHs formed through core-collapse supernovae. Nevertheless, we may expect that such BHs, whichever formed with or without supernovae, move into nearby dense molecular clouds. Thus, in the next subsection, we estimate the plasma accretion rate when such a stellar-mass BH enters a dense gaseous cloud.

To estimate the accretion rate onto a BH, we examine the Bondi accretion rate when a BH moves in a gaseous cloud. To this end, we begin with comparing the sound speed in a typical molecular cloud and the relative velocity between the cloud and the BH.

The sound speed of a cloud with kinetic temperature T (K) can be estimated by

$$\begin{eqnarray}&&{C}_{{\rm{s}}}\sim \sqrt{\displaystyle \frac{{k}_{{\rm{B}}}T}{{m}_{{\rm{p}}}}}=90\,{T}^{1/2}\,{\rm{m}}\,{{\rm{s}}}^{-1},\end{eqnarray} \tag{ 1 }$$

where k_B refers to the Boltzmann constant and m_p is the proton mass. For a typical GMC, we have T < 50 K, and for a typical dark cloud, we have T < 20 K. Thus, we obtain C_s < 630 m s⁻¹ as the upper limit of the sound speed in a molecular cloud.

A typical velocity V of a BH relative to the cloud may be estimated by the kick velocity in a supernova explosion. For a neutron star, it is typically a few hundred kilometers per second. For a BH, a greater fraction of mass is expected to be turned into a compact object; hence, we may expect the typical velocity to be around 10² km s⁻¹. Turbulent velocities measured from molecular line width are usually less than 10 km s⁻¹, and relative velocities among smaller-scale clouds are also within this small range. Thus, we can neglect such random or bulk motions and adopt the supernova kick velocity, V ∼ 10² km s⁻¹, as the typical velocity of a BH with respect to the molecular clouds. If a heavier BH (M > 40 M_⊙) is formed without a supernova explosion, the relative velocity will be less than this value.

On these grounds, we can safely put $V\gg {C}_{{\rm{s}}}$. In this case, the flow becomes supersonic with respect to the BH, and a shock wave is formed behind the hole. Accordingly, the gas particles within the Bondi radius (i.e., within the impact parameter) r_B ∼ GM/V² from the BH are captured, falling onto the BH with the Bondi accretion rate (Bondi & Hoyle 1944)

$$\begin{eqnarray}&&{\dot{M}}_{{\rm{B}}}=4\pi \lambda {({GM})}^{2}{V}^{-3}\rho ,\end{eqnarray} \tag{ 2 }$$

where ρ denotes the mass density of the gas and λ is a constant of order unity. We have λ = 1.12 and 0.25 for an isothermal and adiabatic gas, respectively. Assuming molecular hydrogen gas, and normalizing with the Eddington accretion rate, we obtain the following dimensionless Bondi accretion rate:

$$\begin{eqnarray}&&{\dot{m}}_{{\rm{B}}}=5.39\times {10}^{-9}\lambda {n}_{{{\rm{H}}}_{2}}{M}_{1}{\left(\displaystyle \frac{\eta }{0.1}\right)}^{-1}{\left(\displaystyle \frac{V}{{10}^{2}\mathrm{km}{{\rm{s}}}^{-1}}\right)}^{-3},\end{eqnarray} \tag{ 3 }$$

where ${n}_{{{\rm{H}}}_{2}}$ is measured in units of ${\mathrm{cm}}^{-3}$, η ∼ 0.1 denotes the radiation efficiency of the accretion flow, and M₁ ≡ M/(10 M_⊙).

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 r_B. Thus, we neglect the mass loss as a disk wind between r_B and the innermost region, and evaluate the accretion rate near the BH, $\dot{m}$, with ${\dot{m}}_{{\rm{B}}}$. As will be shown in Section 5.3, the gap of a stellar-mass BH becomes most luminous when $6\times {10}^{-5}\lt \dot{m}\lt 2\times {10}^{-4}$. For $\dot{m}$ to reside in this range, a dense, isothermal molecular cloud core should have a density ${n}_{{{\rm{H}}}_{2}}\gt {10}^{4}\,{\mathrm{cm}}^{-3}$, if V = 100 km s⁻¹. If V = 50 km s⁻¹, however, a lower density, ${n}_{{{\rm{H}}}_{2}}\gt 1.2\times {10}^{3}\,{\mathrm{cm}}^{-3}$, is enough to activate the BH gap.

The Large Area Telescope (LAT) aboard the Fermi space gamma-ray observatory has detected pulsed signals in high-energy (0.1–10 GeV) gamma rays from more than 200 rotation-powered pulsars.⁹ Among them, 20 pulsars exhibit pulsed signals above 10 GeV, including 10 pulsars up to 25 GeV and another two pulsars above 50 GeV. Moreover, more than 99% of the LAT-detected young and millisecond pulsars exhibit phase-averaged spectra that are consistent with a pure-exponential or a subexponential cutoff above the cutoff energies at a few GeV. What is more, 30% of these young pulsars show subexponential cutoff, a slower decay than the pure-exponential functional form. These facts preclude the possibility of emissions from the inner magnetosphere as in the polar cap scenario (Harding et al. 1978; Daugherty & Harding 1982; Dermer & Sturner 1994; Timokhin & Arons 2013; Timokhin & Harding 2015), which predicts superexponential cutoff due to magnetic attenuation. That is, we can conclude that the pulsed emissions are mainly emitted from the outer magnetosphere, which is close to or outside the light cylinder.

One of the main scenarios of such outer-magnetospheric emissions is the outer gap model (Cheng et al. 1986a, 1986b, 2000; Chiang & Romani 1992; Romani 1996; Hirotani 2006). In the present paper, we apply this successful scenario to BH magnetospheres. Although the electrodynamics is mostly common between the pulsar outer gap model and the present BH gap model, there is a striking difference between them. In a pulsar magnetosphere, an outer gap arises because of the convex geometry of the dipolar-like magnetic field in the outer magnetosphere. However, in a BH magnetosphere, a gap arises because of the frame dragging in the vicinity of the event horizon. We describe this BH gap model below.

4.2.1. Background Spacetime Geometry

In a rotating BH magnetosphere, an electron–positron accelerator is formed in the direct vicinity of the event horizon. Thus, we start with describing the background spacetime in a fully general relativistic way. We adopt the geometrized unit, putting c = G = 1, where c and G denote the speed of light and the gravitational constant, respectively. Around a rotating BH, the spacetime geometry is described by the Kerr metric (Kerr 1963). In the Boyer–Lindquist coordinates, it becomes (Boyer & Lindquist 1967)

$$\begin{eqnarray}&&{{ds}}^{2}={g}_{{tt}}{{dt}}^{2}+2{g}_{t\varphi }{dtd}\varphi +{g}_{\varphi \varphi }d{\varphi }^{2}+{g}_{{rr}}{{dr}}^{2}+{g}_{\theta \theta }d{\theta }^{2},\end{eqnarray} \tag{ 4 }$$

where

$$\begin{eqnarray}&&{g}_{{tt}}\equiv -\displaystyle \frac{{\rm{\Delta }}-{a}^{2}{\sin }^{2}\theta }{{\rm{\Sigma }}},\quad {g}_{t\varphi }\equiv -\displaystyle \frac{2{Mar}{\sin }^{2}\theta }{{\rm{\Sigma }}},\end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray}&&{g}_{\varphi \varphi }\equiv \displaystyle \frac{A{\sin }^{2}\theta }{{\rm{\Sigma }}},\quad {g}_{{rr}}\equiv \displaystyle \frac{{\rm{\Sigma }}}{{\rm{\Delta }}},\quad {g}_{\theta \theta }\equiv {\rm{\Sigma }},\end{eqnarray} \tag{ 6 }$$

${\rm{\Delta }}\equiv {r}^{2}-2{Mr}+{a}^{2}$, ${\rm{\Sigma }}\equiv {r}^{2}+{a}^{2}{\cos }^{2}\theta$, and $A\equiv {({r}^{2}+{a}^{2})}^{2}\,-{\rm{\Delta }}{a}^{2}{\sin }^{2}\theta$. At the horizon, we obtain Δ = 0, which gives the horizon radius, ${r}_{{\rm{H}}}\equiv M+\sqrt{{M}^{2}-{a}^{2}}$, where M corresponds to the gravitational radius, ${r}_{{\rm{g}}}\equiv {{GMc}}^{-2}=M$. The spin parameter a becomes a = M for a maximally rotating BH and a = 0 for a nonrotating BH. The spacetime dragging frequency is given by $\omega (r,\theta )=-{g}_{t\varphi }/{g}_{\varphi \varphi }$, which decreases outward as $\omega \propto {r}^{-3}$ at $r\gg {r}_{{\rm{g}}}=M$.

4.2.2. Poisson Equation for the Non-corotational Potential

We assume that the non-corotational potential Φ depends on t and φ only through the form $\varphi -{{\rm{\Omega }}}_{{\rm{F}}}t$, and we put

$$\begin{eqnarray}&&{F}_{\mu t}+{{\rm{\Omega }}}_{{\rm{F}}}{F}_{\mu \varphi }=-{\partial }_{\mu }{\rm{\Phi }}(r,\theta ,\varphi -{{\rm{\Omega }}}_{{\rm{F}}}t),\end{eqnarray} \tag{ 7 }$$

where Ω_F denotes the magnetic field line rotational angular frequency. We refer to such a solution as a "stationary" solution in the present paper.

Gauss's law gives the Poisson equation that describes Φ in a 3D magnetosphere (Hirotani 2006),

$$\begin{eqnarray}&&-\displaystyle \frac{1}{\sqrt{-g}}{\partial }_{\mu }\left(\displaystyle \frac{\sqrt{-g}}{{\rho }_{{\rm{w}}}^{2}}{g}^{\mu \nu }{g}_{\varphi \varphi }{\partial }_{\nu }{\rm{\Phi }}\right)=4\pi (\rho -{\rho }_{\mathrm{GJ}}),\end{eqnarray} \tag{ 8 }$$

where the general relativistic Goldreich–Julian (GJ) charge density is defined as (Hirotani 2006)

$$\begin{eqnarray}&&{\rho }_{\mathrm{GJ}}\equiv \displaystyle \frac{1}{4\pi \sqrt{-g}}{\partial }_{\mu }\left[\displaystyle \frac{\sqrt{-g}}{{\rho }_{{\rm{w}}}^{2}}{g}^{\mu \nu }{g}_{\varphi \varphi }({{\rm{\Omega }}}_{{\rm{F}}}-\omega ){F}_{\varphi \nu }\right].\end{eqnarray} \tag{ 9 }$$

Far away from the horizon, $r\gg M$, Equation (9) reduces to the ordinary, special relativistic expression of the GJ charge density (Goldreich & Julian 1969; Mestel 1971),

$$\begin{eqnarray}&&{\rho }_{\mathrm{GJ}}\equiv -\displaystyle \frac{{\boldsymbol{\Omega }}\cdot {\boldsymbol{B}}}{2\pi c}+\displaystyle \frac{({\boldsymbol{\Omega }}\times {\boldsymbol{r}})\cdot ({\rm{\nabla }}\times {\boldsymbol{B}})}{4\pi c}.\end{eqnarray} \tag{ 10 }$$

Therefore, the corrections due to magnetospheric currents, which are expressed by the second term of Equation (10), are included in Equation (9).

If the real charge density ρ deviates from the rotationally induced GJ charge density, ρ_GJ, in some region, Equation (8) shows that Φ changes as a function of position. Thus, an acceleration electric field, ${E}_{\parallel }=-\partial {\rm{\Phi }}/\partial s$, arises along the magnetic field line, where s denotes the distance along the magnetic field line. A gap is defined as the spatial region in which E_∥ is nonvanishing. At the null-charge surface, ρ_GJ changes sign by definition. Thus, a vacuum gap, in which $| \rho | \ll | {\rho }_{\mathrm{GJ}}|$, appears around the null-charge surface because $\partial {E}_{\parallel }/\partial s$ should have opposite signs at the inner and outer boundaries (Chiang & Romani 1992; Romani 1996; Cheng et al. 2000). As an extension of the vacuum gap, a nonvacuum gap, in which $| \rho |$ becomes a good fraction of $| {\rho }_{\mathrm{GJ}}|$, also appears around the null-charge surface (Section 2.3.2 of HP16), unless the injected current across either the inner or the outer boundary becomes a substantial fraction of the GJ value. If the injected current becomes non-negligible compared to the created current in the gap, the gap centroid position shifts from the null surface; however, the essential gap electrodynamics does not change.

In previous series of our papers (Hirotani & Pu 2016; Hirotani et al. 2016; Song et al. 2017), we have assumed ${\rm{\Delta }}\ll {M}^{2}$ in Equation (8), expanding the left-hand side in the series of Δ/M², and pick up only the leading orders. However, in the present paper, we discard this approximation and consider all the terms that arise at Δ ∼ M² or ${\rm{\Delta }}\gg {M}^{2}$.

It should be noted that ρ_GJ vanishes, and hence the null surface appears near the place where Ω_F coincides with the spacetime dragging angular frequency, ω (Beskin et al. 1992). The deviation of the null surface from this ω(r, θ) = Ω_F surface is, indeed, small, as Figure 1 of Hirotani & Okamoto (1998) indicates. Since ω can match Ω_F only near the horizon, the null surface, and hence the gap, generally appears within one or two gravitational radii above the horizon, irrespective of the BH mass.

4.2.3. Particle Boltzmann Equations

We outline the Boltzmann equations of e^± values, following Hirotani et al. (2017). Imposing a stationary condition, $\partial /\partial t+{{\rm{\Omega }}}_{{\rm{F}}}\partial /\partial \phi =0$, we obtain the following Boltzmann equations:

$$\begin{eqnarray}&&c\cos \chi \displaystyle \frac{\partial {n}_{\pm }}{\partial s}+\dot{p}\displaystyle \frac{\partial {n}_{\pm }}{\partial p}=\alpha ({S}_{\mathrm{IC},\pm }+{S}_{{\rm{p}},\pm }),\end{eqnarray} \tag{ 11 }$$

along each radial magnetic field line on the poloidal plane, where the upper and lower signs correspond to the positrons (with charge q = +e) and electrons (q = −e), respectively, and $p\equiv | {\boldsymbol{p}}| ={m}_{{\rm{e}}}c\sqrt{{\gamma }^{2}-1}$. The left-hand side is in dt basis, where t denotes the proper time of a distant static observer. Thus, the lapse α is multiplied on the right-hand side because both S_IC and S_p are evaluated in the zero-angular-momentum observer (ZAMO). Dimensionless lepton distribution functions per magnetic flux tube are defined by

$$\begin{eqnarray}&&{n}_{\pm }\equiv \displaystyle \frac{2\pi {ce}}{{{\rm{\Omega }}}_{{\rm{F}}}B}{N}_{\pm }(r,\theta ,\gamma ),\end{eqnarray} \tag{ 12 }$$

where N₊ and N₋ designate the distribution functions of positrons and electrons, respectively; γ refers to these leptons' Lorentz factor, and $B\equiv | {\boldsymbol{B}}|$. It is convenient to include the curvature emission as a friction term on the left-hand side; in this case, we obtain

$$\begin{eqnarray}&&\dot{p}\equiv {{qE}}_{\parallel }\cos \chi -\displaystyle \frac{{P}_{\mathrm{SC}}}{c},\end{eqnarray} \tag{ 13 }$$

where the pitch angle is assumed to be $\chi =0$ for outwardly moving positrons and χ = π for inwardly moving electrons. The curvature radiation force is given by (e.g., Harding 1981) ${P}_{\mathrm{SC}}/c=2{e}^{2}{\gamma }^{4}/(3{R}_{{\rm{c}}}{}^{2})$.

The inverse-Compton (IC) collision terms are expressed as

$$\begin{eqnarray}\begin{array}{rcl}{S}_{\mathrm{IC}} & \equiv & -{\displaystyle \int }_{{\epsilon }_{\gamma }\lt \gamma }d{\epsilon }_{\gamma }{\eta }_{\mathrm{IC}}^{\gamma }({\epsilon }_{\gamma },\gamma ,{\mu }_{\pm }){n}_{\pm }\\ & & +\ {\displaystyle \int }_{{\gamma }_{i}\gt \gamma }d{\gamma }_{i}{\eta }_{\mathrm{IC}}^{{\rm{e}}}({\gamma }_{i},\gamma ,{\mu }_{\pm }){n}_{\pm },\end{array}\end{eqnarray} \tag{ 14 }$$

where the IC redistribution function is defined by

$$\begin{eqnarray}&&{\eta }_{\mathrm{IC}}^{\gamma }\equiv (1-\beta {\mu }_{\pm }){\int }_{{E}_{\min }}^{{E}_{\max }}{{dE}}_{{\rm{s}}}\displaystyle \frac{{{dF}}_{{\rm{s}}}}{{{dE}}_{{\rm{s}}}}\displaystyle \frac{d{\epsilon }_{\gamma }^{* }}{d{\epsilon }_{\gamma }}{\int }_{-1}^{1}d{{\rm{\Omega }}}_{\gamma }^{* }\displaystyle \frac{d{\sigma }_{\mathrm{KN}}^{* }}{d{\epsilon }_{\gamma }^{* }d{{\rm{\Omega }}}_{\gamma }^{* }}.\end{eqnarray} \tag{ 15 }$$

${m}_{{\rm{e}}}{c}^{2}{\epsilon }_{\gamma }$ denotes the upscattered γ-ray energy. The asterisk denotes the quantity evaluated in the electron (or positron) rest frame, and $d{\sigma }_{\mathrm{KN}}^{* }/d{\epsilon }_{\gamma }^{* }d{{\rm{\Omega }}}_{\gamma }^{* }$ denotes the Klein–Nishina differential cross section. Energy conservation gives

$$\begin{eqnarray}&&{\eta }_{\mathrm{IC}}^{{\rm{e}}}({\gamma }_{i},\gamma ,{\mu }_{\pm })={\eta }_{\mathrm{IC}}^{\gamma }({\gamma }_{i}-\gamma ,{\gamma }_{i},{\mu }_{\pm }),\end{eqnarray} \tag{ 16 }$$

where γ_i denotes the Lorentz factor before collision and μ₊ (or μ₋) denotes the cosine of the collision angle with the soft photon for outwardly moving electrons (or inwardly moving positrons). For more details, see Section 3.2.2 of Hirotani et al. (2003). The soft photons are emitted by the hot electrons within a radiatively inefficient accretion flow (RIAF). The effect of this inhomogeneous and anisotropic soft photon field is included in the differential soft photon flux, dF_s/dE_s, through the correction factor f_riaf (Section 3.4 of Hirotani et al. 2017). That is, we put ${{dF}}_{{\rm{s}}}/{{dE}}_{{\rm{s}}}={f}_{\mathrm{riaf}}\cdot {({{dF}}_{{\rm{s}}}/{{dE}}_{{\rm{s}}})}_{0}$.

The photon–photon pair creation term becomes

$$\begin{eqnarray}&&{S}_{{\rm{p}}}\equiv \int d{\nu }_{\gamma }{\alpha }_{\gamma \gamma }\displaystyle \frac{2\pi e}{{{\rm{\Omega }}}_{{\rm{F}}}B}\int \displaystyle \frac{{I}_{\omega }}{{\hslash }\omega }d{{\rm{\Omega }}}_{\gamma },\end{eqnarray} \tag{ 17 }$$

where

$$\begin{eqnarray}&&{\alpha }_{\gamma \gamma }=(1-{\mu }_{\pm }){\int }_{{E}_{\mathrm{th}}}^{\infty }{{dE}}_{{\rm{s}}}\displaystyle \frac{{{dF}}_{{\rm{s}}}}{{{dE}}_{{\rm{s}}}}\displaystyle \frac{d{\sigma }_{\gamma \gamma }}{d\gamma }.\end{eqnarray} \tag{ 18 }$$

The γ-ray specific intensity I_ω is integrated over the γ-ray propagation solid angle Ω_γ. For details, see Section 3.2.2 of Hirotani et al. (2003). Note that dF_s/dE_s in both ${\eta }_{\mathrm{IC}}^{\gamma }$ and α_γγ is evaluated in ZAMO (Section 3.4 of Hirotani et al. 2017).

It is noteworthy that the charge conservation ensures that the dimensionless total current density (per magnetic flux tube), ${j}_{\mathrm{tot}}\equiv \int (-{n}_{+}-{n}_{-})d\gamma$, conserves along the flow line. If we denote the created current density as J_cr, the injected current density across the inner and outer boundaries as J_in and J_out, respectively, and the typical GJ value as ${J}_{\mathrm{GJ}}\equiv {{\rm{\Omega }}}_{{\rm{F}}}{B}_{{\rm{H}}}/(2\pi )$, we obtain j_tot = j_cr + j_in + j_out, where ${j}_{\mathrm{cr}}\equiv {J}_{\mathrm{cr}}/{J}_{\mathrm{GJ}}$, ${j}_{\mathrm{in}}\equiv {J}_{\mathrm{in}}/{J}_{\mathrm{GJ}}$, and ${j}_{\mathrm{out}}\equiv {J}_{\mathrm{out}}/{J}_{\mathrm{GJ}}$.

4.2.4. Radiative Transfer Equation

In the same manner as Hirotani et al. (2016), 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 (Hirotani 2013),

$$\begin{eqnarray}&&\displaystyle \frac{{{dI}}_{\omega }}{{dl}}=-{\alpha }_{\omega }{I}_{\omega }+{j}_{\omega },\end{eqnarray} \tag{ 19 }$$

where ${dl}=\sqrt{{g}_{{rr}}}{dr}$ refers to the distance interval along the ray in ZAMO, and α_ω and j_ω are the absorption and emission coefficients evaluated in ZAMO, respectively. We consider both photon–photon and magnetic 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 Sections 4.2 and 4.3 of HP16 and Section 5.1.5 of H16.

4.2.5. Boundary Conditions

4.2.6. Gap Closure Condition

4.2.7. Advection-dominated Accretion Flow (ADAF)

Let us begin with the examination of the E_∥(r, θ) distribution along the individual magnetic field lines that are radial on the meridional plane. As demonstrated in Figure 3 of Hirotani et al. (2018), E_∥ peaks along the rotation axis, because the magnetic fluxes concentrate toward the rotation axis as the BH spin approaches its maximum value (i.e., as $a\to {r}_{{\rm{g}}}$; Komissarov & McKinney 2007; Tchekhovskoy et al. 2010). Therefore, to consider the greatest gamma-ray flux, we focus on the emission along the rotation axis, θ = 0°. The acceleration electric field, E_∥, decreases slowly outside the null surface in the same way as pulsar outer gaps (Hirotani & Shibata 1999). This is because the 2D screening effect of E_∥ works when the gap longitudinal (i.e., radial) width becomes non-negligible compared to its trans-field (i.e., meridional) thickness. In addition, in the present work, we pick up all the terms that contribute not only near the horizon (i.e., ${\rm{\Delta }}\ll {M}^{2}$) but also away from it (i.e., Δ ∼ M² or ${\rm{\Delta }}\gg {M}^{2}$). As a result, the exerted E_∥ is reduced from the case of ${\rm{\Delta }}\ll {M}^{2}$, which particularly reduces the curvature luminosity compared to our previous works (Lin et al. 2017).

In Figure 1, we present E_∥(s, θ = 0°) solved at four dimensionless accretion rates, $\dot{m}={10}^{-3.50}$, 10^−3.75, 10^−4.00, and 10^−4.25, where $s\equiv r-{r}_{0}(\theta )$ denotes the distance from the null surface, r = r₀(θ), and θ = 0° is adopted. As pointed out in previous BH gap models, the potential drop increases with decreasing $\dot{m}$. However, if the accretion further decreases as $\dot{m}\lt {10}^{-4.25}$, there are no stationary gap solutions. Below this lower bound accretion rate, the gap solution becomes inevitably nonstationary.

Zoom In Zoom Out Reset image size

We next consider the electrons' distribution function. Because of the the negative E_∥, electrons and positrons are accelerated outward and inward, respectively. As Figure 2 shows, the electrons' Lorentz factors concentrate around 3 × 10⁶ owing to the curvature radiation drag. At the same time, electrons distribute at lower Lorentz factors with a broad plateau typically between 6 × 10⁴ and 2 × 10⁶. Electrons stay at such relatively lower Lorentz factors because of the inverse Compton drag. Since the Klein–Nishina cross section increases with decreasing Lorentz factors, such lower-energy electrons with 6 × 10⁴ < γ < 2 × 10⁶ efficiently contribute to the VHE emission via the inverse Compton scatterings.

Zoom In Zoom Out Reset image size

Let us examine the gamma-ray spectra. In Figure 3, 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, $\theta \lt 15^\circ$, and that the gap luminosity rapidly decreases at θ < 30° if the gap equatorial boundary is located at θ = 60°. In what follows, to estimate the maximally possible VHE flux, we consider the emission along the rotation axis, θ = 0°.

Zoom In Zoom Out Reset image size

We also consider how the the SED depends on the BH spin. In Figure 4, we show the SEDs for a = 0.99M, 0.90M, and 0.50M; in each panel, SEDs for the four discrete accretion rates, $\dot{m}$ = 10^−3.50, 10^−3.75, 10^−4.00, and 10^−4.25, are plotted. It is clear that the gap luminosity increases with decreasing $\dot{m}$. It also follows that the gap emission could be detectable with CTA if a > 0.90M, provided that the distance is within 1 kpc and we observe nearly face-on. However, if the BH is moderately rotating as a = 0.50M, it is very difficult to detect its emission, unless it is located within 0.3 kpc.

Zoom In Zoom Out Reset image size

In Figure 5, we depict the individual emission components, selecting the case of $\dot{m}={10}^{-4.00}$. We find that the primary curvature component (magenta dashed line) dominates between 5 MeV and 0.5 GeV, while the primary IC component (magenta dotted-dashed line) dominates above 5 GeV. The secondary IC component (blue triple-dotted-dashed line) appears between 0.5 and 5 GeV. The primary IC component suffers absorption above 0.1 TeV.

Zoom In Zoom Out Reset image size

Let us demonstrate that the gap solution exists in a wide range of the created electric current within the gap, ${j}_{\mathrm{cr}}\equiv {J}_{\mathrm{cr}}/{J}_{\mathrm{GJ}}$, and that the resultant gamma-ray spectrum little depends on j_cr. In Figure 6, we show the solved E_∥(s) (left panels) and SEDs (right panels) for three discrete j_cr values: from the top, they correspond to j_cr = 0.3, 0.5, and 0.9. The case of j_cr = 0.7 is presented as Figure 1 and the top panel of Figure 4. 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. Note that there exist no stationary gap solutions if the created current density is set to be super-GJ, j_cr > 1.

Zoom In Zoom Out Reset image size

The gap solution exists in a wide range of the injected electric currents across the inner and outer boundaries. In this section, we consider only the current injected across the outer boundary, because the positrons created below the separation surface (Figure 2 of Hirotani & Pu 2016, hereafter HP16) may enter the gap across the outer boundary, whereas the electrons created below the inner boundary will fall onto the horizon. In the left panels of Figure 7, E_∥ is plotted as a function of the Boyer–Lindquist radial coordinate from the null-charge surface. The red dotted, blue dashed, black solid, and green dotted-dashed curves correspond to the cases of $\dot{m}={10}^{-3.50}$, 10^−3.75, 10^−4.00, and 10^−4.25, respectively. From the top, each panel show the results for j_out = 0.2, 0.4, 0.6, and 0.8. In the right panels, SEDs are presented for the same set of j_out values. The case of j_out = 0 is presented as Figure 1 and the top panel of Figure 4.

Zoom In Zoom Out Reset image size

It follows from the left panels that E_∥ shifts inward with increasing J_out and that $| {E}_{\parallel }|$ increases as the gap position shifts inward, as suggested by the outer gap solutions for rotation-powered pulsars (Hirotani & Shibata 2001a, 2001b, 2002). Note that the photon–photon collision mean free path does not decrease as the gap approaches the horizon, which forms a striking contrast from pulsars. As a result, the potential drop also increases as the gap shifts inward in the case of BHs. (In the case of pulsars, the soft photons are emitted from the cooling neutron star surface. Thus, as the gap approaches the stellar surface, the head-on collisions of inward γ-rays and the outward surface X-rays become more efficient, decreasing the gap longitudinal size and hence the potential drop.) In the present BH cases, the increased $| {E}_{\parallel }|$ and the potential drop lead to an increased gap luminosity and photon energies in the local reference frame. However, due to the gravitational redshift, the photon energy reduces for a distant static observer. Accordingly, as the right panels show, the final SEDs little change if j_out changes from 0 (top panel of Figure 4) to 0.8 (bottom right panel of Figure 7), although the E_∥(s) distribution changes significantly as Figure 1 and the left panels of Figure 7 show.

We compare the gamma-ray emission scenarios from molecular clouds (Table 1). In the protostellar jet scenario (Bosch-Ramon et al. 2010), jets are ejected from massive protostars to interact with the surrounding dense molecular clouds, leading to an acceleration of electrons and protons at the termination shocks. Accordingly, the size of the emission region becomes comparable to the jet transverse thickness at the shock. In the hadronic cosmic-ray scenario (Ginzburg & Syrovatskii 1964; Blandford & Eichler 1987), protons and helium nuclei are accelerated in the supernova shock fronts, a portion of which propagate into dense molecular clouds. As a result of the proton–proton (or nuclear) collisions, neutral pions are produced and decay into gamma rays, whose spectrum becomes a single power law between 0.001 and 100 TeV. Since this interaction takes place most efficiently in a dense gaseous region, the size of the gamma-ray image will become comparable to the core of a dense molecular cloud. In the leptonic cosmic-ray scenario (Aharonian et al. 1997; Hillas et al. 1998; van der Swaluw et al. 2001), electrons are accelerated at pulsar wind nebulae or shell-type SNRs and radiate radio/X-rays and gamma rays via synchrotron and IC processes, respectively. Since the cosmic microwave background radiation provides the main soft photon field in the interstellar medium, the size may be comparable to the plerions, whose size increases with the pulsar age. In the BH gap scenario (see Section 1 for references), emission size does not exceed 10r_g. Since the angular resolution of the CTA is about five times better than the current IACTs, we propose that we can discriminate the present BH gap scenario from the three above-mentioned scenarios by comparing the gamma-ray image and spectral properties. Namely, if a VHE source has a point-like morphology like HESS J1800–2400C in a gaseous cloud (Section 2), and if the spectrum has two peaks around 0.01–1 GeV and 0.01–1 TeV but shows a (synchrotron) power-law component in neither radio nor X-ray wavelengths, we consider that the present scenario accounts for its emission mechanism.

Table 1.  Gamma-ray Emission Models from Molecular Clouds^a

Model	Emission Processes (Spectral Shape; Energy Range)	Size (cm)
Protostellar jets	e− synchrotron (power law; 10−6–102 eV);	1016–1017
	e− bremsstrahlung (power law; 0.1 MeV–TeV);
	pp collisions, π0 decaysb (power law; GeV–TeV)
Cosmic-ray hadrons	pp collisions, ${\pi }^{0}$ decaysb (power law; GeV–100 TeV)	1018–1019
Cosmic-ray leptonsc	e− synchrotron (power law; 10−6–102 eV);	1018–1020
	e− IC scatterings (broad peak; GeV–10 TeV)
BH gap	e− curvature process (broad peak; 0.01–1 GeV);	107
	e− IC scatterings (sharp peak; around 0.1 TeV)

Notes.

^aSee Section 6.1 for references. ^bNeutral pion (π⁰) decays follow proton–proton collisions. ^cGamma rays emitted by cosmic-ray leptons may not be associated with molecular clouds. Nevertheless, it is one of the main scenarios of the VHE emissions from massive-star-forming regions.

Download table as:  ASCIITypeset image

Although the magnetic (i.e., one-photon) pair production is also taken into account, most electron–positron pairs are found to be produced via photon–photon (i.e., two-photon) collisions, which take place via two paths. One path is through the collisions of two MeV photons, both of which were emitted from the equatorial ADAF. Another path is through the collisions of TeV and eV photons; the former photons were emitted by the gap-accelerated leptons via inverse Compton process, while the latter were emitted from the ADAF via the synchrotron process. There is, indeed, a third path, in which the gap-emitted GeV curvature photons collide with the ADAF-emitted keV inverse Compton photons; however, this path is negligible, particularly when $\dot{m}\ll 1$. If the pairs are produced via TeV–eV collisions (i.e., via the second path) outside the gap outer boundary, they have outward ultrarelativistic momenta to easily "climb up the hill" of the potential k₀ (Takahashi et al. 1990; see also Figure 2 of HP16) and propagate to large distances without turning back. However, if the pairs are produced via MeV–MeV collisions (i.e., via the first path), they are produced with subrelativistic outward momenta; thus, they eventually return to fall onto the horizon owing to the strong gravitational pull inside the separation surface (Figure 2 of HP16). When the returned pairs arrive at the gap outer boundary, only positrons can penetrate into the gap because of E_∥ < 0. Accordingly, electrons accumulate at the boundary, whose surface charge leads to the jump of the normal derivative of E_∥. Thus, although the stationary gap solutions show that the γ-ray spectrum little depends on the injected current density (Section 5.5), the gap solution inevitably becomes time dependent owing to the increasing discontinuity of $| {{dE}}_{\parallel }/{dr}|$ with an accumulated surface charge (in this case, electrons) at the outer boundary. If the injected current is much small compared to the GJ current, the time dependence will be mild. However, if the injected current becomes a good fraction of the GJ current, the assumption of the stationarity becomes invalid, as pointed out by Levinson & Segev (2017). In this sense, a caution should be made in the applicability of the stationary solutions presented in this paper, when the injected current is non-negligible compared to the current created within the gap.

In this subsection, we consider the stability of our stationary gap solutions. In the case of pulsar polar caps, it has been revealed that the pair-production cascade takes place in a highly time-dependent way by particle-in-cell (PIC) simulations (Timokhin 2010; Timokhin & Arons 2013; Timokhin & Harding 2015). Moreover, in the case of BH magnetospheres, it is recently demonstrated that a gap exhibits rapid spatial and temporal oscillations of the magnetic-field-aligned electric field and current by 1D PIC simulations (Levinson & Segev 2017; Chen et al. 2018). It is, however, out of the scope of the present paper to perform a PIC simulation or a linear perturbation analysis to examine the stability of a stationary gap solution. Instead, we will qualitatively discuss why we seek stationary solutions, comparing with pulsar outer-magnetospheric and polar cap gaps.

We start by discussing the pulsar outer(-magnetospheric) gaps, because they have essentially the same electrodynamics as BH gaps. Since the neutron star's dipole magnetic field lines have convex geometry, the magnetic field becomes perpendicular with respect to the star's rotation axis at a good fraction of the so-called "light-cylinder radius." In this case, if the magnetosphere is highly vacuum in the sense $| \rho | \ll | {\rho }_{\mathrm{GJ}}|$ in Equation (8), ρ − ρ_GJ ≈ −ρ_GJ > 0 (or <0) holds in the lower (or upper) half of the gap. As a result, E_∥ has a positive (or negative) gradient in the lower (or upper) half; thus, the acceleration electric field naturally closes. Note that E_∥ > 0 is realized when the magnetic axis resides in the same hemisphere as the rotation axis. Without loss of any generality, we can adopt such a positive E_∥ in the outer gap model.

Because ${E}_{\parallel }\gt 0$, positrons (or electrons) are accelerated outward (or inward). Thus, ρ has a positive gradient along individual magnetic field lines. When the gap closure condition is satisfied, there exists a stationary solution whose ρ/B distribution can be illustrated as the left panel of Figure 8. If pair production increases perturbatively from this stationary solution, ρ − ρ_GJ decreases its absolute value, which leads to a decrease of E_∥ due to the reduced gradient of $| \rho -{\rho }_{\mathrm{GJ}}|$, and hence a decrease of the pair production (right panel of Figure 8). Note that this negative feedback effect works because of E_∥ > 0, which stems from the fact that there exists a null-charge surface in the gap. As a result, the outer gap solutions exist for a wide range of pulsar parameters such as the period, period derivative, neutron star surface temperature, inclination of the magnetic axis with respect to the star's rotation axis, and magnetospheric current, from young, middle-aged to millisecond pulsars. An analogous (but still qualitative) argument of gap stability is possible if we use the gap closure condition; see Section 5.2 of Hirotani (2001) for details.

Zoom In Zoom Out Reset image size

Because null-charge surfaces also exist around rotating BHs, the gap electrodynamics little changes between pulsar outer gaps and BH gaps. Around a rotating BH, a null surface is formed near the event horizon by the frame dragging. However, around a rotating neutron star, it is formed in the outer magnetosphere (i.e., far away from the neutron star) by the convex magnetic field geometry. Accordingly, a negative feedback effect also works in BH gaps. That is, the gap closure condition, which is required for a BH gap to be stationary, is accommodated for a wide range of magnetospheric current values, as explicitly demonstrated in Section 5.

Since there frequently appears to be confusion between pulsar polar cap and outer gap (and hence BH gap) electrodynamics, particularly on the stability argument, it is helpful to describe also the pulsar polar cap accelerator. In a pulsar polar cap, there exists no null surface. Thus, if a 3D polar cap region is charge starved in the sense $| \rho -{\rho }_{\mathrm{GJ}}| \ll {\rho }_{\mathrm{GJ}}$, a positive −ρ_GJ leads to a negative E_∥ when the magnetic axis resides in the same hemisphere as the rotation axis. Accordingly, electrons are drawn from the neutron star surface as a space-charge-limited flow. In the direct vicinity of the neutron star surface, the nonrelativistic electrons produces a large negative ρ (the green vertical lines along the ordinates in Figure 9) such that ρ − ρ_GJ ≈ ρ < 0. In the outer part of a polar cap accelerator, on the other hand, relativistic electrons produce a moderate negative charge density such that ρ − ρ_GJ > 0 (left panel of Figure 9). Thus, E_∥ has a negative gradient along the magnetic field in the direct vicinity of the neutron star, but it has a positive gradient near the upper boundary where ρ − ρ_GJ ≈ 0. Accordingly, we obtain a negative E_∥ in the pulsar polar caps.

Zoom In Zoom Out Reset image size

However, if a small-amplitude pair production takes place in the polar gap, inwardly migrating positrons (or outwardly migrating electrons) result in an increased (or decreased) ρ − ρ_GJ in the lower (or upper) half of the gap (right panel of Figure 9). Accordingly, $| {E}_{\parallel }|$ increases (or decreases) in the lower part (or the uppermost part). The increased $| {E}_{\parallel }|$ further enhances pair production in the uppermost part, and $| {E}_{\parallel }|$, and hence the pair production, increases with time. Because of this positive feedback effect, instability sets in, as demonstrated by PIC simulations (Timokhin 2010; Timokhin & Arons 2013). In short, pulsar polar cap accelerators are inherently unstable for pair production because of E_∥ < 0, which stems from the fact that there exists no null-charge surface in the pulsar polar cap region.

On these grounds, we cannot readily conclude that the stationary BH gap solutions presented in this paper are unstable because of the highly time-dependent nature of pulsar polar gaps. A careful examination with a PIC simulation is needed for BH gaps, as performed recently by Levinson & Segev (2017). Since the saturated solution given by Levinson & Segev (2017) is much less violently time dependent compared to pulsar polar cap accelerators (Timokhin & Arons 2013), we have sought stationary solutions of BH gaps as the first step in the present paper.