One-dimensional General Relativistic Particle-in-cell Simulations of Stellar-mass Black Hole Magnetospheres: A Semianalytic Model of Gamma-Rays from Gaps

We first demonstrate the results for the low-∣J₀∣ cases, i.e., j₀ = −1/2π. Figure 1 shows the evolution of the system for model LA1. An electric field triggered by the inward and outward drifts of the plasma develops around the null charge surface (r_null ≈ 2.0 for a_* = 0.9), and a quasiperiodic oscillation is observed, similarly as found in K20 and K22. Electrons linearly accelerated by the electric field become mono-energetic around u∼γ_e ∼ 10⁷ during the peak time of the electric field (t = 58.3r_g/c).

Zoom In Zoom Out Reset image size

Figure 2 shows the luminosities of curvature emission, IC emission, and particle kinetic energies for model LA1. Here, we normalize them by the BZ luminosity, written as ${L}_{{\rm{BZ}}}=({\kappa }_{{\rm{B}}}/4\pi c){\omega }_{{\rm{H}}}^{2}{{\rm{\Phi }}}_{{\rm{H}}}=({\kappa }_{{\rm{B}}}\pi c/4){{a}_{\ast }}^{2}{B}_{{\rm{H}}}^{2}{r}_{{\rm{g}}}^{2}$, where κ_B is set to be 0.053 for the split-monopole magnetic field configuration (Tchekhovskoy et al. 2010). We observe luminosity oscillations synchronized with that of the electric field, similarly as the SMBH cases. The peak values are L_cur,pk ∼ 10⁻² L_BZ, L_ic,pk ∼ 10⁻⁴ L_BZ, and L_kin,pk ∼ 10⁻⁵ L_BZ.

Zoom In Zoom Out Reset image size

We then investigate the dependence of the gap dynamics on the parameters related to the seed photon properties, τ₀, ${\epsilon }_{\min }$, and p. In the left panel of Figure 3 we show the ratio of rgw peak luminosity of curvature radiation to the BZ luminosity L_cur,pk/L_BZ for models LA–C for various values of τ₀. ${L}_{\mathrm{cur},\mathrm{pk}}/{L}_{\mathrm{BZ}}\propto {\tau }_{0}^{\alpha }(\alpha \sim -8/3)$ is found for models LA1–5, which indicates the sensitive dependence of the gap dynamics on the γ γ pair production efficiency, is similar to that shown in K22. We also find that the effect of changing the spectral hardness p on the gap behavior is minor. This is because most of the IC scatterings occur in the Klein–Nishina (KN) regime, for which the γ γ pair production optical depth is independent on p (see K20). A change in the peak energy ${\epsilon }_{\min }$, however, strongly affects the dynamics similarly to that found in K20 and K22. For ${\epsilon }_{\min }={10}^{-5}$, the width and period of the gaps are ∼10 times larger than those for ${\epsilon }_{\min }={10}^{-6}$ in the cases of τ₀ = 30 and 100 (compare models LC1 and LC2 with models LA1 and LA3, respectively). ⁸ In contrast, a very narrow and highly intermittent gap appears for ${\epsilon }_{\min }={10}^{-5}$ in the case of τ₀ = 300 (model LC3).

Zoom In Zoom Out Reset image size

We next perform the simulations for j₀ = −1/2 and −1. In these cases, the gap does not appear to be clearly periodic and the electric field strength is weaker than the case of j₀ = −1/2π for the same τ₀, p, and ${\epsilon }_{\min }$. This is because the higher electric current ∣J₀∣ leads to a higher density of the plasma flowing in the gap, which makes it easier to screen the electric field.

In Figure 3 we compare L_cur,pk and L_ic,pk for different values of j₀. L_cur,pk for j₀ = −1/2 is 10⁻¹–10⁻² times smaller than that for j₀ = −1/2π for all of the τ₀ values. This results from the weaker electric field and lower efficiency of particle acceleration. L_ic,pk, on the contrary, is slightly larger for lower j₀, due to an increase of the plasma number density in the gap. Its scaling is roughly L_ic,pk ∝ ∣j₀∣.

The dependence of the particle number density on j₀ at the peak time of electric field is ${\left.({\rm{\Delta }}/{{\rm{\Delta }}}_{\min })n/{n}_{\mathrm{GJ},\min }\right|}_{r={r}_{\mathrm{null}}}\sim 0.08$ for j₀ = −1/2π, ∼0.3 for j₀ = −1/2, and ∼0.5 for j₀ = −1, where ${n}_{\mathrm{GJ},\min }={\rho }_{\mathrm{GJ}}({r}_{\min })/e$ and ${{\rm{\Delta }}}_{\min }={\rm{\Delta }}({r}_{\min })$. Note that the density does not appear to be dependent on τ₀, p, and ${\epsilon }_{\min }$. This j₀ dependence can be understood analytically, by using Ampere's law (Equation (2)). When the development of E_r stops at the peak time, we have

$$\begin{eqnarray}&&{\rm{\Sigma }}{j}^{r}({r}_{\mathrm{null}})\approx {J}_{0}.\end{eqnarray} \tag{ 8 }$$

Since the positrons (electrons) carrying j^r are accelerated to relativistic velocities v_i ≈ −1(+1) around the null surface, j^r (r_null) can be written by using the relation between the number densities of electrons and positrons, n₊(r_null) ≈ n₋(r_null) = n(r_null)

$$\begin{eqnarray}&&{j}^{r}({r}_{\mathrm{null}})=-\displaystyle \frac{{\rm{\Delta }}}{\sqrt{A}}2{en}({r}_{\mathrm{null}}).\end{eqnarray} \tag{ 9 }$$

Then, we obtain

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{\Delta }}({r}_{\mathrm{null}})}{{{\rm{\Delta }}}_{\min }}\displaystyle \frac{n({r}_{\mathrm{null}})}{{n}_{\mathrm{GJ},\min }}\approx -0.58{j}_{0},\end{eqnarray} \tag{ 10 }$$

where we use Equations (3) and (4). Equation (10) is consistent with our findings in the simulations. We utilize this relation in our semianalytic model described in the next section.

4.1.1. Equation of Motion in the Gap

We first solve the equation of motion of an outgoing electron. We define r_in as the initial radius. Our 1D simulations assume u_ϕ = 0 for all the particles, and the electron is rapidly accelerated to a highly relativistic velocity. Then we have $\alpha {dt}\approx \sqrt{{g}_{\mathrm{rr}}}{dr}$ for the electron's motion, so that Equation (5) can be rewritten as

$$\begin{eqnarray}&&\displaystyle \frac{d{\gamma }_{{\rm{e}}}}{{dr}}=-\displaystyle \frac{\sqrt{{g}_{{\rm{rr}}}}}{{m}_{{\rm{e}}}{c}^{3}}\left[{{eE}}_{r}\left(r,{r}_{{\rm{in}}},{B}_{{\rm{H}}})c+{P}_{{\rm{cur}}}({\gamma }_{{\rm{e}}})+{P}_{{\rm{ic}}}({\gamma }_{{\rm{e}}}\right)\right],\end{eqnarray} \tag{ 11 }$$

where we have added the IC cooling term P_ic. (In the simulations IC cooling is treated via the Monte Carlo approach.) We confirmed that the gravitational force (first term of the right-hand side of Equation (5)) is negligible for the energy range achieved in our simulations, γ _e ≲ 10⁷.

The simulations show the number densities of electrons and positrons ${n}_{+}\approx {n}_{-}\ll {n}_{\mathrm{GJ},\min }$ around the peak time of the electric field in the gap. These enable us to approximate the electric field by the vacuum (j^t = 0) solution of Gauss' equation (Equation (1))

$$\begin{eqnarray}&&{E}_{{\rm{r}}}(r,{r}_{\mathrm{in}},{B}_{{\rm{H}}})=-\displaystyle \frac{4\pi }{\sqrt{A}}{\int }_{{r}_{\mathrm{in}}}^{r}{dr}^{\prime} \,{\rm{\Sigma }}{\rho }_{\mathrm{GJ}}(r^{\prime} ,{B}_{{\rm{H}}}),\end{eqnarray} \tag{ 12 }$$

where we have imposed the boundary condition E_r(r_in) = 0.

The IC emitting power is written as

$$\begin{eqnarray}&&{P}_{{\rm{ic}}}({\gamma }_{{\rm{e}}},{\tau }_{0})={\sigma }_{{\rm{ic}}}{n}_{{\rm{s}}}({\epsilon }_{min}){\epsilon }_{{\rm{ic}}}{m}_{{\rm{e}}}{c}^{3}=\displaystyle \frac{{\tau }_{0}}{{r}_{{\rm{g}}}}\displaystyle \frac{{\sigma }_{{\rm{ic}}}}{{\sigma }_{{\rm{T}}}}{\epsilon }_{{\rm{ic}}}{m}_{{\rm{e}}}{c}^{3}.\end{eqnarray} \tag{ 13 }$$

Here σ_ic is the KN cross section

$$\begin{eqnarray}&&\begin{array}{l}{\sigma }_{\mathrm{ic}}\\ =\displaystyle \frac{3}{8}{\sigma }_{{\rm{T}}}\displaystyle \frac{1}{{x}^{2}}\left[4+\displaystyle \frac{2{x}^{2}\left(1+x\right)}{{\left(1+2x\right)}^{2}}+\displaystyle \frac{{x}^{2}-2x-2}{x}\mathrm{log}(1+2x)\right],\end{array}\end{eqnarray} \tag{ 14 }$$

where $x={\gamma }_{{\rm{e}}}{\epsilon }_{{\rm{\min }}}$. _ic is the energy of the scattered photon, approximately given by

$$\begin{eqnarray}&&\begin{array}{c}{\epsilon }_{{\rm{ic}}}=0.3\times \left\{\begin{array}{c}{\gamma }_{{\rm{e}}}^{2}{\epsilon }_{min}({\gamma }_{{\rm{e}}}{\epsilon }_{min}\lesssim 1),\\ {\gamma }_{{\rm{e}}}({\gamma }_{{\rm{e}}}{\epsilon }_{min}\gtrsim 1).\end{array}\right.\end{array}\end{eqnarray} \tag{ 15 }$$

The Lorentz factor of an electron at the kth step r_k of the integration of Equation (11) is determined as ${\gamma }_{{\rm{e}},k}\,\equiv {\gamma }_{{\rm{e}}}({r}_{k})=\min \{{\gamma }_{{\rm{e}},\mathrm{acc},k},{\gamma }_{{\rm{e}},\mathrm{cur},k},{\gamma }_{{\rm{e}},\mathrm{ic},k}\}$ (see Figure 5), where γ_e,acc,k is the value for the acceleration-dominant regime

$$\begin{eqnarray}&&{\gamma }_{{\rm{e}},\mathrm{acc},k}={\gamma }_{{\rm{e}},k-1}-\sqrt{{g}_{\mathrm{rr}}}\displaystyle \frac{{{eE}}_{{\rm{r}}}}{{m}_{{\rm{e}}}c}{dr}.\end{eqnarray} \tag{ 16 }$$

γ_e,cur,k and γ_e,ic,k are the values determined by the balances of the acceleration with the curvature and IC coolings, respectively,

$$\begin{eqnarray}&&-{{eE}}_{{\rm{r}}}(\,{r}_{k-1},{r}_{{\rm{in}}},{B}_{{\rm{H}}})c=\left\{\begin{array}{l}{P}_{{\rm{cur}}}({\gamma }_{{\rm{e}},{\rm{cur}},k}),\\ {P}_{{\rm{ic}}}({\gamma }_{{\rm{e}},{\rm{ic}},k},{\tau }_{0}).\end{array}\right.\end{eqnarray} \tag{ 17 }$$

4.1.2.  γ γ Pair Production

Next, we compute γ γ pair production in the gap. The number of photons scattered by the electron as it propagates from r_k−1 to r_k is calculated by considering that seed photons with ${\epsilon }_{\min }$ are most abundant when

$$\begin{eqnarray}&&\delta {N}_{\gamma ,k}=\displaystyle \frac{{dr}}{{l}_{\mathrm{ic}}({\epsilon }_{\min },{\gamma }_{{\rm{e}},k},{\tau }_{0})},\end{eqnarray} \tag{ 18 }$$

where l_ic is the mean free path of IC scattering

$$\begin{eqnarray}&&{l}_{{\rm{ic}}}({\epsilon }_{min},{\gamma }_{{\rm{e}}},{\tau }_{0})=\displaystyle \frac{1}{\sqrt{{g}_{{\rm{rr}}}}}\displaystyle \frac{1}{{n}_{{\rm{s}}}({\epsilon }_{min}){\sigma }_{{\rm{ic}}}}=\displaystyle \frac{{r}_{{\rm{g}}}}{\sqrt{{g}_{{\rm{rr}}}}{\tau }_{0}}\displaystyle \frac{{\sigma }_{{\rm{T}}}}{{\sigma }_{{\rm{ic}}}}.\end{eqnarray} \tag{ 19 }$$

We consider that all the scattered photons annihilate with the seed photons after propagating about the mean free path for the γ γ pair annihilation l_{γ γ}

$$\begin{eqnarray}\begin{array}{rcl}{l}_{\gamma \gamma }({\epsilon }_{\min },{\gamma }_{{\rm{e}}},{\tau }_{0}) & = & \displaystyle \frac{1}{\sqrt{{g}_{\mathrm{rr}}}}\displaystyle \frac{1}{{n}_{{\rm{s}}}({\epsilon }_{2}){\sigma }_{\gamma \gamma }}\\ & = & \displaystyle \frac{10\,{\epsilon }_{\mathrm{ic}}{\epsilon }_{2}{r}_{{\rm{g}}}}{\sqrt{{g}_{\mathrm{rr}}}{\tau }_{0}}{\left(\displaystyle \frac{{\epsilon }_{2}}{{\epsilon }_{\min }}\right)}^{p},\end{array}\end{eqnarray} \tag{ 20 }$$

where ₂ is the energy of the annihilating seed photon, defined as ${\epsilon }_{2}=\max \{{\epsilon }_{\min },{\epsilon }_{\mathrm{ic}}^{-1}\}$, and we use ${\sigma }_{\gamma \gamma }\sim 0.1{\sigma }_{{\rm{T}}}{({\epsilon }_{\mathrm{ic}}{\epsilon }_{2})}^{-1}$ for the γ γ pair annihilation cross section. Thus, the number of pairs created by photons scattered between the k − 1 step and the kth step is

$$\begin{eqnarray}&&\delta {N}_{\mathrm{pair},k}({r}_{\mathrm{ppf},k},{\gamma }_{{\rm{e}},k})=\delta {N}_{\gamma ,k},\end{eqnarray} \tag{ 21 }$$

where the radius of the "pair production front" r_ppf of the scattered photons is given by

$$\begin{eqnarray}&&{r}_{\mathrm{ppf},k}\equiv {r}_{\mathrm{ppf}}({r}_{k},{\gamma }_{{\rm{e}},k})={r}_{k}+{l}_{\gamma \gamma }({\epsilon }_{\min },{\gamma }_{{\rm{e}},k}).\end{eqnarray} \tag{ 22 }$$

The total number of newly created pairs from r_in to position r is calculated as

$$\begin{eqnarray}&&{N}_{\mathrm{pair}}(r)=\displaystyle \sum _{{r}_{\mathrm{ppf},k}\leqslant r}\delta {N}_{\mathrm{pair},k}({r}_{\mathrm{ppf},k}).\end{eqnarray} \tag{ 23 }$$

4.1.3. Gap Width and Luminosity

We iteratively calculate the total amount of created pairs for trial values of r_in with r_H < r_in < r_null, until we find the outer boundary of the gap r_out that satisfies the following two criteria: (1) the total number of newly created pairs between r_in and r_out is

$$\begin{eqnarray}&&{N}_{\mathrm{pair}}({r}_{\mathrm{out}})=2.0.\end{eqnarray} \tag{ 24 }$$

It is expected that this number should be of order unity to have the electric field stop growing. We found that with a value of 2.0 our semianalytic model well fits the simulation results. (2) The gap is symmetric in the ξ coordinate (see Appendix A) is, i.e.,

$$\begin{eqnarray}&&\displaystyle \frac{\xi ({r}_{\mathrm{out}})+\xi ({r}_{\mathrm{in}})}{2}=\xi ({r}_{\mathrm{null}}),\end{eqnarray} \tag{ 25 }$$

which is satisfied in the simulation results, as shown in Figure 4.

Zoom In Zoom Out Reset image size

For the solution of r_in and r_out obtained above, we calculate the peak curvature luminosity by

$$\begin{eqnarray}&&\begin{array}{c}{L}_{{\rm{cur}}}=2\pi {\int }_{{r}_{{\rm{in}}}}^{{r}_{{\rm{out}}}}{\alpha }^{2}(r){P}_{{\rm{cur}}}\left({\gamma }_{{\rm{e}}}\left(r\right)\right)n({r}_{{\rm{null}}}){\rm{\Sigma }}(r){dr}\\ \,\approx 0.39\displaystyle \frac{{e}^{2}c}{{r}_{{\rm{g}}}^{2}}\displaystyle \frac{{{\rm{\Delta }}}_{min}}{\left.{\rm{\Delta }}({r}_{{\rm{null}}}\right)}{n}_{{\rm{GJ}},min}{\int }_{{r}_{{\rm{in}}}}^{{r}_{{\rm{out}}}}{\alpha }^{2}{\gamma }_{{\rm{e}}}^{4}{\rm{\Sigma }}{dr},\end{array}\end{eqnarray} \tag{ 26 }$$

where we have assumed that the number density of emitting particles is constant throughout the gap, and evaluate it by n(r_null) (see Equation (10)).

Figure 5 shows the solutions of our semianalytic model for τ₀ = 30 and 300, with the other parameters ${\epsilon }_{\min }={10}^{-6}$, p = 2, and B_H = 2π × 10⁷ G (same as the LA set of models). The top panels show the energy evolution of the electron. The vertical dashed lines indicate r_in (magenta) and r_out (cyan) determined as the solutions, and the white region corresponds to the gaps. The dashed curves show the evolution of the Lorentz factor calculated by Equations (16) and (17) for the solutions. The determined γ_e at each r is shown by the red solid line. We also show r_ppf for each r and γ_e.

Zoom In Zoom Out Reset image size

As can be seen in Figure 5, the maximum electron Lorentz factors ${\gamma }_{{\rm{e}},\ \max }\sim {10}^{7}$ for τ₀ = 30 and ${\gamma }_{{\rm{e}},\max }\sim 2\times {10}^{6}$ for τ₀ = 300 are determined by curvature cooling and the IC cooling, respectively. These result from the different particle energy dependences of the two cooling mechanisms as well as the τ₀ dependence of the IC-scattering mean free path; for τ₀ = 30, l_ic becomes larger compared with the case of τ₀ = 300, which leads to lower δ N_pair. Thus, a solution with a larger gap width is obtained, in which the maximum Lorentz factor reaches γ _e ∼ 10⁷. In this regime curvature cooling is more efficient than IC cooling because ${P}_{\mathrm{cur}}\propto {\gamma }_{{\rm{e}}}^{4}$ and the IC scattering is deeply in the KN regime. On the other hand, for τ₀ = 300, ${\gamma }_{{\rm{e}},\max }$ is regulated by IC cooling in the Thomson regime or marginally in the KN regime. This transition of the dominant cooling process matches with the trend of the gap dynamics and energetics in the simulations, in which we see that the gaps become narrower for higher τ₀ and L_cur,pk ≲ L_ic,pk ∼ 10⁻⁴ L_BZ for τ₀ ≳ 100 (see Figure 3).

We compare the simulation results and semianalytic model solutions of ${\gamma }_{{\rm{e}},\max }$ and L_cur,pk/L_BZ for j₀ = −1/2π in Figure 6. For models LA and LB, the semianalytic solutions show that ${\gamma }_{{\rm{e}},\max }$ and L_cur,pk/L_BZ are decreasing functions of τ₀ with a break around τ₀ ∼ 100–200, and well fit the simulation results particularly for τ₀ < 200. For models LC2 and LC3 we obtain ∼10² times higher L_cur,pk/L_BZ than the simulation results, and for model LC1 we could not find a solution, but the overall dependence of ${\gamma }_{{\rm{e}},\max }$ and L_cur,pk/L_BZ on τ₀ are similar to what we observe in the simulations.

Zoom In Zoom Out Reset image size

The breaks seen at τ₀ ≳ 100 in all the solutions are due to the transition from the curvature-dominant state to the IC-dominant state, as explained in Section 4.2. The luminosity drop along with this transition is prominent for models LC1–3, which is consistent with the simulation results. We can conclude that the energy and τ₀ dependences of the two cooling rates affect the dynamics and energetics of the gaps.

In the simulations for j₀ = −1/2, the time evolution of the gap does not exhibit clear periodicity. Also, we observe n₊ ≠ n₋ around the peak time except at the middle point of the gap, hence the assumption of a vacuum electric field is not necessarily valid. Nevertheless, we can build an effective model consistent with the average values of ${\gamma }_{{\rm{e}},\max }$ and L_cur,pk/L_BZ in the simulation results by adding some modifications to the model for j₀ = −1/2π. We use 1/10 times weaker E_r than that calculated by Equation (12) and use a different condition on N_pair instead of Equation (24)

$$\begin{eqnarray}&&{N}_{\mathrm{pair}}({r}_{\mathrm{out}})=\displaystyle \frac{1}{\pi | {j}_{0}| }.\end{eqnarray} \tag{ 27 }$$

Note that the condition of Equation (24) for j₀ = −1/2π is consistent with Equation (27). We present the model solutions for j₀ = −1/2 (models I1–3) in Figure 6. As seen, the solutions well fit the simulation results.

In this work we have assumed the absence of sufficient plasma injection into the magnetosphere and considered only the pair production process by the annihilation of IC-scattered photons and soft photons from the MAD. Here we discuss other pair production processes in the magnetosphere that might screen the gap. One possibility we need to consider is a steady plasma injection produced by the annihilation of MeV photons from the MAD (Levinson & Rieger 2011; Kimura & Toma 2020). The expected pair number density around the null surface is ${n}_{\pm }\sim ({r}_{\mathrm{null}}{\sigma }_{\gamma \gamma }/16{\pi }^{2}{R}^{4}{c}^{2}){\left({L}_{\mathrm{MeV}}/{E}_{\mathrm{MeV}}\right)}^{2}\simeq {10}^{5}$ ${{ \mathcal R }}_{1}^{-4}{M}_{1}{n}_{\mathrm{ISM},2}^{2}\,{\mathrm{cm}}^{-3}$, where ${L}_{\mathrm{MeV}}\simeq {10}^{32}{M}_{1}^{2}{n}_{\mathrm{ISM},2}\,\mathrm{erg}\,{{\rm{s}}}^{-1}$ and E_MeV are the spectral luminosity (see Figure 8) and the characteristic energy in the MeV energy band, respectively. This is much lower than the GJ density ${n}_{\mathrm{GJ}}({r}_{\mathrm{null}})\sim {{\rm{\Omega }}}_{{\rm{F}}}B({r}_{\mathrm{null}})/2\pi {ec}\simeq 1.4\,\times \,{10}^{8}\,({\phi }_{\mathrm{BH}}/50){M}_{1}^{-1}{n}_{\mathrm{ISM},2}^{1/2}\,{\mathrm{cm}}^{-3}$, where $B{({r}_{\mathrm{null}})\sim {B}_{{\rm{H}}}({r}_{\mathrm{null}}/{r}_{{\rm{H}}})}^{-2}\sim {B}_{{\rm{H}}}/2$ is the magnetic field strength around the null charge surface. Furthermore, n_±/n_GJ and τ₀ have the same dependence, $\propto {M}^{2}{n}_{\mathrm{ISM}}^{3/2}$, hence MeV photon annihilation will not inject a sufficient amount of pairs to screen the gap for the ranges of M and n_ISM for which the gamma-ray radiation is bright.

Alternatively, as briefly mentioned in Section 1, sporadic flares by magnetic reconnection around the equatorial plane at r_rec ∼ 2r_g can produce MeV photons in the magnetosphere (Chashkina et al. 2021; Kimura et al. 2022; Ripperda et al. 2022; Chen et al. 2023; Hakobyan et al. 2023). During each flare, the magnetic energy is dissipated in a localized reconnection region of size l_rec ∼ r_g with reconnection velocity β_rec ∼ 0.1 (e.g., Guo et al. 2020), and the accelerated electrons are efficiently cooled by emitting MeV photons. Then the luminosity is ${L}_{\mathrm{MeV}}\sim 2(B{({r}_{\mathrm{rec}})}^{2}/8\pi ){l}_{\mathrm{rec}}^{2}{\beta }_{\mathrm{rec}}c\sim ({\phi }_{\mathrm{BH}}^{2}/64{\pi }^{3}){\beta }_{\mathrm{rec}}\dot{M}{c}^{2}$ ≃$\,8\times {10}^{33}{({\phi }_{\mathrm{BH}}/50)}^{2}{\beta }_{\mathrm{rec},-1}{M}_{1}^{2}{n}_{\mathrm{ISM},2}\,\mathrm{erg}\,{{\rm{s}}}^{-1}$. Thus, one obtains ${n}_{\pm }\simeq 6\times {10}^{12}{(d/{r}_{{\rm{g}}})}^{-4}{({\phi }_{\mathrm{BH}}/50)}^{4}{\beta }_{\mathrm{rec},-1}^{2}{M}_{1}^{5}$ ${n}_{\mathrm{ISM},2}^{2}{{cm}}^{-3}$ in a similar manner to the above discussion for MAD MeV photons, replacing R with the distance between the emitting region and the gap d ∼ r_g. This is ∼10⁵ times higher than n_GJ(r_null), so that the gap will be completely screened during the flare. Ripperda et al. (2022) suggests that each flare is followed by a quiescent phase of timescale ∼10³ r_g/c. The injected pairs will escape from the null surface in a timescale ∼r_g/c after the flare, and hence, the gap is expected to develop in the quiescent phase.

The possibility of single photon pair production, or magnetic pair production, in the magnetosphere also needs to be considered. This opacity can be written as a function of the interacting photon energy _γ , the local magnetic field B, and the angle between the photon momentum and the magnetic field ψ (Erber 1966) as

$$\begin{eqnarray}&&{\alpha }_{B}({\epsilon }_{\gamma },B,\psi )\approx 0.23\displaystyle \frac{{\alpha }_{f}}{-{\lambda }_{c}}b\sin \psi \exp \left(-\displaystyle \frac{4}{3\chi }\right),\end{eqnarray} \tag{ 30 }$$

where α_f = e²/ℏ c ≃ 1/137 is the fine structure constant, $-{\lambda }_{c}$ = ℏ/m_e c ≃ 3.9 × 10⁻¹¹ cm is the reduced Compton wavelength, b = B/B_q is the local magnetic field strength normalized by ${B}_{q}\equiv e/{\alpha }_{f}-{\lambda }_{c}^{2}\simeq 4.4\times {10}^{13}\,G$, and $\chi =(1/2){\epsilon }_{\gamma }b\sin \psi$. B(r_null) ≃ 5 × 10⁶ B_H,7G for M ∼ 10 M_⊙ and n_ISM ∼ 10² cm⁻³, for which the gamma-ray emission is bright and the magnetic field is strong (see Figure 7 and the left panel of Figure 8). The particle Lorentz factor in this case is γ_e ∼ 10⁷, and thus the photons produced in the gap via IC scattering will be _γ ∼ 3.0 × 10⁶ (see Equation (15)). ¹¹ The corresponding optical depth of the magnetic pair production is τ_B ∼ α_B r_g ≃ 4 × 10³ for $\sin \psi =1.0$, and 2 × 10⁻² for $\sin \psi =0.4$. In the latter estimate we have approximated ψ by the field pitch angle ∣B_ϕ ∣/B_r ∼ 0.4 measured in the ZAMO frame (see Appendix A in Kimura et al. 2022). Such a sensitive dependence of τ_B on ψ implies that the detailed magnetic field structure and the general relativistic lensing effect significantly alter the pair loading efficiency. To analyze them is beyond the scope of this paper. Nevertheless, we confirmed that magnetic pair production is unimportant for the models shown in the right panel of Figure 8 and all the models in Figure 9, regardless of ψ: for smaller n_ISM, the magnetic field is weaker. For ${\tau }_{0}\propto {n}_{\mathrm{ISM}}^{3/2}\gg 30$, _γ ≲ 10⁶ and the gap width is much smaller than r_g.

Our expected main targets are Fermi-LAT unidentified objects (unIDs). The catalog from 12 yr operations of Fermi-LAT (4FGL-DR3; Ajello et al. 2022) contains more than 2000 objects lacking associations with known pulsars/supernovae/AGNs. Nearly a half of them are lying at a low Galactic latitude (∣b∣ < 10), aligning with our expectation of molecular cloud associations. Our IBH candidates can be distinguished by a hard spectral index (${dN}/{{dE}}_{\gamma }\propto {E}_{\gamma }^{-{\rm{\Gamma }}}$ (Γ ∼ 2/3)) with a break around the GeV–TeV range and a point-source like morphology. Since objects with a hard photon index are rare in the entire population of unIDs, we may be able to narrow down the candidates solely by the Fermi-LAT data. Furthermore, we can identify IBHs by taking correlations with survey catalogs conducted at other wavelengths. Our model also predicts bright radiation from the MADs in the optical to MeV bands (see Figure 8; Kimura et al. 2021b). An optical signal from a MAD resembles that from a white dwarf, which can be found by Gaia (Gaia Collaboration et al. 2018), and MAD X-ray signals can be detected by X-ray observations including eROSITA (Predehl et al. 2021) and ROSAT (Truemper 1982). Thus, we could specify our IBH candidates by examining associations among Fermi-LAT unIDs, X-ray sources, and Gaia white-dwarf-like objects.

Time variations of the gamma-ray luminosity due to a fluctuation in the mass accretion rate can also be useful to distinguish IBHs from other objects. ${L}_{\mathrm{cur},\mathrm{pk}}\propto {\tau }_{0}^{-8/3}{L}_{\mathrm{BZ}}\,\propto {\dot{M}}^{-3}(\propto {n}_{\mathrm{ISM}}^{-3})$ and ${L}_{\mathrm{ic},\mathrm{pk}}\propto {L}_{\mathrm{BZ}}\propto \dot{M}$, hence occasional luminosity variations by a factor of ≲10² (as shown in Figure 9) are expected. The timescale of such variations would be no shorter than the interaction timescale of IBHs and ISM clouds, ${T}_{\mathrm{var}}\sim {GM}/{v}^{3}\simeq 2.1\times {10}^{7}{M}_{1}{(v/40\,\mathrm{km}\,{{\rm{s}}}^{-1})}^{-3}\,{\rm{s}}$, where v ≃ 40 km s⁻¹ is the IBH proper motion velocity. Multiepoch observations with the Cherenkov Telescope Array (Cherenkov Telescope Array Consortium et al. 2019) can be utilized to capture luminosity variations in such cases. Fermi-LAT can also detect such gamma-ray variability (Abdollahi et al. 2017, 2020) if the source is bright.

Hirotani et al. (2016, 2017), Song et al. (2017), and Hirotani et al. (2018) also examined gap GeV–TeV radiation for various masses of BHs with their analytical 2D gap models. They assume that the gap is steadily maintained. However, various numerical simulations have shown that the gap is time variable, and our model reproduces the characteristics of the gap at the peak time. These authors obtain a nonzero charge density in the gap, while we assume the vacuum solution of the electric field, consistent with the simulation results. Such differences in the model settings and assumptions result in various differences in the calculations results. For example, the IC component in their calculations peaks at a higher energy range with a higher luminosity than our model. Furthermore, the luminosities do not significantly depend on the magnetospheric current in their model (see Figure 6 of Hirotani et al. 2018).

It is uncertain whether local stellar-mass BHs can acquire a high spin. One possibility is spin-up via a binary BH merger (Zlochower & Lousto 2015), but it would be still challenging to achieve an extremely high spin rate as ${a}_{\ast }=0.9$. We thus need to take into account the fraction of rapidly rotating stellar-mass BHs when discussing the expected number of IBH detections. In addition, it would be important to understand the BH spin dependence of the gap dynamics and radiation features.

To calculate the spectrum of IC emission escaping from the soft photon field in the magnetosphere (r ≤ R = 10r_g), we have assumed that L_ic,pk ∼ constant outside the gap by neglecting synchrotron loss (Section 5 and Appendix B.2). However, secondary pairs have momenta in the direction of incident IC photons, which is not necessarily aligned with the local magnetic field vector, so that they could lose their momenta perpendicular to the local field via synchrotron emission. Indeed, for secondary pairs with a finite pitch angle α and Lorentz factor ${\gamma }_{\sec }\sim {\epsilon }_{\mathrm{ic}}/2\sim {10}^{6}$, we have ${t}_{\mathrm{cool}}^{\mathrm{syn}}\ll {t}_{\mathrm{cool}}^{\mathrm{ic}}$. ¹² The synchrotron photon's typical energy is ${\epsilon }_{\mathrm{syn}}=(3{heB}/4\pi {m}_{{\rm{e}}}^{2}{c}^{3}){\gamma }_{\sec }^{2}\sin \alpha \simeq 1.7\times {10}^{5}\sin \alpha (B/5.1\times {10}^{6}{\rm{G}}){\gamma }_{\sec ,6}^{2}$, thus they will have a large mean free path for pair annihilation $l/R\sim {\tau }_{0}^{-1}{({\epsilon }_{\mathrm{syn}}^{-1}/{\epsilon }_{\min })}^{2}$ ∼$\,1.2{(\sin \alpha )}^{-2}{\epsilon }_{\min ,-6}^{-2}$ ${(B/5.1\times {10}^{6}{\rm{G}})}^{-2}{\gamma }_{\sec ,6}^{-4}{({\tau }_{0}/30)}^{-1}$. Since $\sin \alpha \lt 1$, most of the synchrotron photons may escape from the region of r < R. Since _syn is also in the GeV–TeV energy band, the peak luminosity of the IC plus synchrotron radiation may not substantially change from our rough estimate of IC luminosity in Section 5. Detailed consideration of multidimensional magnetic field structure and α would be required to calculate the IC and synchrotron spectra more accurately.

Our 1D calculations do not take account of a general relativistic lensing effect on the photon transfer. This affects the light curve, which will be considered in future work. Differences between 1D and 2D GRPIC simulation results on the gap dynamics should also be further investigated. We have shown that J₀ (or the toroidal magnetic field) substantially affects the gap dynamics. The toroidal field is determined by external pressure on the magnetosphere. Including effects of variable external pressure in 2D simulations would be interesting.

We calculate the observed spectrum of curvature emission from the magnetospheric gaps at the peak time, for which we may assume that the emitting electrons have a mono-energetic energy distribution ${\gamma }_{{\rm{e}}}\approx {\gamma }_{{\rm{e}},max}$. The radiated power of one electron with Lorentz factor γ_e per unit energy is

$$\begin{eqnarray}&&{P}_{\epsilon }^{{\rm{cur}}}({\gamma }_{{\rm{e}}})=\displaystyle \frac{\sqrt{3}{e}^{2}{\gamma }_{{\rm{e}}}}{{{hr}}_{{\rm{g}}}}F\left(\displaystyle \frac{\epsilon }{{\epsilon }_{{\rm{cur}}}}\right),\end{eqnarray} \tag{ B1 }$$

where

$$\begin{eqnarray}&&F(x)=x{\int }_{x}^{\infty }d\xi {K}_{5/3}(\xi ),\end{eqnarray} \tag{ B2 }$$

K_5/3(ξ) is the modified Bessel function of 5/3 order, and

$$\begin{eqnarray}&&{\epsilon }_{{\rm{cur}}}=\displaystyle \frac{3}{4\pi }\displaystyle \frac{{hc}}{{r}_{{\rm{g}}}}{\gamma }_{{\rm{e}}}^{3},\end{eqnarray} \tag{ B3 }$$

is the characteristic energy (Rybicki & Lightman 1986). Let ${N}_{{\rm{e}}}^{{\rm{cur}}}$ denote the number of electrons with ${\gamma }_{{\rm{e}}}\approx {\gamma }_{{\rm{e}},max}$, then we can derive the specific luminosity as

$$\begin{eqnarray}&&{L}_{\epsilon }^{\mathrm{cur}}=\displaystyle \frac{\sqrt{3}{e}^{2}{\gamma }_{{\rm{e}},\max }}{{{hr}}_{{\rm{g}}}}{N}_{e}^{\mathrm{cur}}F\left(\displaystyle \frac{\epsilon }{{\epsilon }_{\mathrm{cur}}({\gamma }_{{\rm{e}},\max })}\right).\end{eqnarray} \tag{ B4 }$$

${N}_{{\rm{e}}}^{{\rm{cur}}}$ is determined so that ${L}_{\mathrm{cur},\mathrm{pk}}={\int }_{0}^{\infty }{L}_{\epsilon }^{\mathrm{cur}}d\epsilon$. The attenuation of gamma-rays via pair annihilation is evaluated by multiplying ${L}_{\epsilon }^{\mathrm{cur}}$ by the factor $\exp [-{\tau }_{\gamma \gamma }(\epsilon ,{\epsilon }_{{\rm{s}}})]$. Here τ_{γ γ} is the pair annihilation optical depth (see, e.g., Kimura et al. 2019, and references therein) at r = R, within which the seed photons are dense (R = 10r_g is the size of the MAD, see Section 5). We show only the attenuated spectrum in Figure 8. The overall spectrum is gravitationally redshifted by α(r_emit), where we set the emission radius of curvature photons to be r_emit = r_null.

For τ₀ ≳ 10, for which the magnetosphere does not approach the vacuum, the primary IC photons produced by electrons with ${\gamma }_{{\rm{e}}}\approx {\gamma }_{{\rm{e}},max}$ annihilate with the seed photons to create e⁺ e⁻ pairs in the gap. These pairs produce gamma-rays via IC scatterings even outside the gap. Here we do not solve the detailed radiation transfer, but simply estimate the spectrum of the IC emission escaping from r = R of the magnetosphere.

The IC gamma-rays induce pair cascade while propagating in the magnetosphere. (Note that we do not consider the cascade of IC gamma-rays injected into the MAD.) We find ${t}_{\mathrm{cool}}^{\mathrm{ic}}\ll {t}_{\mathrm{dyn}}\ll {t}_{\mathrm{cool}}^{\mathrm{cur}}$ outside the gap, where ${t}_{\mathrm{cool}}^{\mathrm{ic}}$, ${t}_{\mathrm{cool}}^{\mathrm{cur}}$, and t_dyn are the IC cooling timescale, the curvature cooling timescale, and the dynamical timescale, respectively. The particle kinetic luminosity thus keeps subdominant outside the gap as shown in Figure 2, and energy loss via curvature emission is negligible, so that one may assume ${L}_{{\rm{ic}},{\rm{pk}}}\sim {\rm{constant}}$ outside the gap. The simulation results indicate L_ic,pk ∼ 10⁻³∣j₀∣L_BZ. The characteristic energy of escaping IC photons _ic,esc can be estimated as the energy for which the pair annihilation optical depth τ_{γ γ} (_ic,esc, _s) = 1 (see Timokhin & Harding 2015; Kisaka & Tanaka 2017; K20). The Lorentz factor of electrons that produce IC photons of energy _{ic, esc} is evaluated as ${\gamma }_{\mathrm{ic},\mathrm{esc}}=\min \{{\gamma }_{{\rm{e}},\max },\sqrt{{\epsilon }_{\mathrm{ic},\mathrm{esc}}/0.3{\epsilon }_{\min }}\}$.

The spectral luminosity is then given by

$$\begin{eqnarray}&&{L}_{\epsilon }^{{\rm{ic}}}=4{\sigma }_{{\rm{T}}}{N}_{{\rm{e}}}^{{\rm{ic}}}{\int }_{0}^{1}{dx}\,g(x){F}_{{\epsilon }_{{\rm{s}}}}(x),x=\displaystyle \frac{\epsilon }{4{\gamma }_{{\rm{ic}},{\rm{esc}}}^{2}{\epsilon }_{{\rm{s}}}},\end{eqnarray} \tag{ B5 }$$

where $g(x)=1\,+\,x\,+\,2x\mathrm{ln}x-2{x}^{2}$ is a kernel in the relativistic limit (Blumenthal & Gould 1970; Rybicki & Lightman 1986) and ${F}_{{\epsilon }_{{\rm{s}}}}(x)$ denotes the seed photon flux. ${N}_{{\rm{e}}}^{{\rm{ic}}}$ is determined to satisfy ${L}_{\mathrm{ic},\mathrm{pk}}={\int }_{0}^{\infty }{L}_{\epsilon }^{\mathrm{ic}}d\epsilon$. To evaluate the attenuation at > _ic,esc, we multiply ${L}_{\epsilon }^{\mathrm{ic}}$ by $[1-\exp (-{\tau }_{\gamma \gamma }(\epsilon ,{\epsilon }_{{\rm{s}}}))]/{\tau }_{\gamma \gamma }$. We show only the attenuated spectrum in Figure 8. The spectrum is gravitationally redshifted by α(r_emit), where we set r_emit = R.