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

In the absence of a sufficient amount of plasma injection into the black hole (BH) magnetosphere, the force-free state of the magnetosphere cannot be maintained, leading to the emergence of strong, time-dependent, longitudinal electric fields (i.e., spark gaps). Recent studies of supermassive BH magnetospheres using analytical methods and particle-in-cell (PIC) simulations propose the possibility of efficient particle acceleration and consequent gamma-ray emission in the spark gap. In this work, we perform 1D general relativistic PIC simulations to examine the gamma-ray emission from stellar-mass BH magnetospheres. We find that intermittent spark gaps emerge and particles are efficiently accelerated in a similar manner to the supermassive BH case. We build a semianalytic model of the plasma dynamics and radiative processes, which reproduces the maximum electron energies and peak gamma-ray luminosities of the simulation results. Based on this model, we show that the gamma-ray signals from stellar-mass BHs wandering through the interstellar medium could be detected by gamma-ray telescopes such as the Fermi Large Area Telescope or the Cherenkov Telescope Array.

Continuous injection of plasma into the magnetosphere is required to activate a steady force-free jet.A possible injection process is electron-positron (e + e − ) pair production by ambient MeV-energy photons from the accretion disk surrounding the magnetosphere, but it could be insufficient to maintain the force-free state in the case of a low accretion rate (Beskin et al. 1992;Hirotani & Okamoto 1998;Levinson & Rieger 2011;Mościbrodzka et al. 2011;Hirotani & Pu 2016;Wong et al. 2021;Yao et al. 2021; see Kimura & Toma 2020;Kuze et al. 2022 for a possible contribution of GeV photons).Magnetic reconnection around the equatorial plane close to the BH may also cause e + e − pair loading (Kimura et al. 2022;Ripperda et al. 2022;Chen et al. 2023;Hakobyan et al. 2023) or direct injection of accreting plasma.However, the former mechanism can inject a large amount of plasma only during a flaring state, and the latter mechanism may not screen the entire region of the magnetosphere (Niv et al. 2023).
Similar spark gaps may arise in charge-starved magnetospheres formed by accretion of interstellar medium (ISM) material onto stellar-mass BHs, which are the source candidates of unidentified gamma-ray objects (Hirotani et al. 2018).The size of the system, expected magnetic field strength, and characteristics of soft photons from the accretion disk differ significantly from SMBH magnetospheres, and thus the plasma dynamics for stellar-mass BHs should be examined separately.Hirotani et al. (2016), Song et al. (2017), andHirotani et al. (2018) present analytic steady gap models for stellar-mass BHs, but it is unlikely that the steady solutions are realized over broad ranges of parameters (Levinson & Segev 2017).Simulation-based modeling of time-dependent gaps is required to predict their emission properties.
In the present paper we examine the dynamics of a local spark gap of stellar-mass BHs and the subsequent radiation characteristics, by using 1D GRPIC simulations. 6We broadly explore the dependence of the gap activity on the ambient photon characteristics and the magnetospheric current (Section 3).Furthermore, we construct a semianalytic model of the particle acceleration and radiation processes in the gap, which reproduces the maximum particle energies and the peak gamma-ray luminosities in the simulations (Section 4).We also calculate the spectra of escaping gamma-rays in some cases with our semianalytic model (Section 5) and then briefly discuss detectability (Section 6).

The Numerical Simulation
We use the 1D GRPIC code described in Levinson & Cerutti (2018), Kisaka et al. (2020), andKisaka et al. (2022; the latter two are hereafter denoted as K20 and K22, respectively).In this code, we assume that the gap activity does not significantly affect the global magnetospheric structure, where the magnetic field is set to have a split-monopole configuration.Under this assumption, we solve Maxwell's equations for the longitudinal electric field E r measured in the frame of a zero angular momentum observer (ZAMO) in a background Kerr spacetime, given in Boyer-Lindquist coordinates (see Appendix A) as and where j t is the charge density and j r is the current density.The Goldreich-Julian (GJ) charge density ρ GJ is written as where B H denotes the magnetic field strength at the horizon, r g = GM/c 2 is the gravitational radius of the BH with mass M (G is the gravitational constant), r r a 1 1 is the horizon radius, ω H = a * c/2r H is the angular velocity of the BH (a * ≡ a/r g is the dimensionless spin parameter), and Ω F = ω H /2 is the angular velocity of the magnetic surface.J 0 is the magnetospheric current parameterized by where j 0 = −1 corresponds to the magnetospheric current in the BZ force-free solution (for the counterclockwise rotation). 7We note that the null charge surface, where ρ GJ = 0, r = r null (θ), is nearly spherical.The equation of motion for the electrons and positrons is where q i is the charge of the ith particle, u g u i i r rr

=
is the normalized four-velocity component in the r direction measured by a ZAMO, γ i is the Lorentz factor, v i = u i /γ i is the three-velocity component, and P e v r 2 3 is the curvature power with the curvature radius fixed as r g .For photons we solve where p j m ˜(μ = t, r) is the momentum components of the jth photon measured by the ZAMO.
We consider two main processes of photon interactions by using a Monte Carlo approach.One is the inverse Compton (IC) scattering of soft "seed" photons from the accretion disk and the other is γγ pair production by scattered photons and seed photons.We assume the intensity of the seed photon field to be steady, isotropic, and homogeneous with a simple powerlaw spectrum I I p s s 0 s min min s max ), where ò s is the seed photon energy normalized by the electron rest mass energy m e c 2 .The fiducial optical depth τ 0 is related to the intensity normalization I 0 as n r rI h c 4 0 s min T g T g 0 , where σ T is the Thomson cross section and h is the Planck constant.We compute the momenta of IC-scattered photons by Equation (6), while not those of curvature photons.The curvature luminosity at any given time is calculated by summing all the contribution of all particles in the simulation box, L r P , where r i is the position of the ith particle.
We apply this model to a BH with M = 10 M e and a * = 0.9.As for B H , we implicitly assume that the accretion disk surrounding the magnetosphere is in a highly magnetized state (MAD; Narayan et al. 2003;McKinney et al. 2012), likewise in our previous simulations for SMBHs (K20, K22).MADs are expected to be formed around stellar-mass BHs when the mass accretion rates  M are much lower than the Eddington rate Cao 2011;Ioka et al. 2017;Kimura et al. 2021bKimura et al. , 2021c)), and then we can write (Tchekhovskoy et al. 2011)  = -as a fiducial value of the peak energy, which corresponds to the optical energy band, consistent with the calculation of photon spectra from stellar-mass BH MADs (see Kimura et al. 2021b).
The computational domain is set from  r r 1.5 min g to  r r 4.3 max g , so that r min is fairly close to the horizon (r H ; 1.43r g for a * = 0.9).We impose open boundary conditions for both boundaries.The inclination of the simulating region with respect to the BH rotation axis is set as θ = 30°.In every run the number of grid cells exceeds 32,000, which was sufficient to resolve the skin depth, and the initial particles per cell number is 45 so that the simulation guarantees convergence (see K20).The time step is set to satisfy the Courant-Friedrichs-Lewy condition at the inner boundary.

Cases of Low
We first demonstrate the results for the low-|J 0 | cases, i.e., j 0 = −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 monoenergetic around u∼γ e ∼ 10 7 during the peak time of the electric field (t = 58.3rg /c).
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 , 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 −2 L BZ , L ic,pk ∼ 10 −4 L BZ , and L kin,pk ∼ 10 −5 L BZ .
We then investigate the dependence of the gap dynamics on the parameters related to the seed photon properties, τ 0 , 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 τ 0 .
) 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 min  , however, strongly affects the dynamics similarly to that found in K20 and K22.For 10 min 5  = -, the width and period of the gaps are ∼10 times larger than those for 10 min 6  = -in the cases of τ 0 = 30 and 100 (compare models LC1 and LC2 with models LA1 and LA3, respectively). 8In contrast, a very narrow and highly intermittent gap appears for 10 min 5  = -in the case of τ 0 = 300 (model LC3).

Cases of Intermediate and High |J 0 |
We next perform the simulations for j 0 = −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 0 = −1/2π for the same τ 0 , p, and min  .This is because the higher electric current |J 0 | 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 0 .L cur,pk for j 0 = −1/2 is 10 −1 -10 −2 times smaller than that for j 0 = −1/2π for all of the τ 0 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 0 , due to an increase of the plasma number density in the gap.Its scaling is roughly L ic,pk ∝ |j 0 |.
The dependence of the particle number density on j 0 at the peak time of electric field is n n 0.08 Note that the density does not appear to be dependent on τ 0 , p, and min  .This j 0 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 j r J .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 Then, we obtain 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.

Semianalytic Model of the Gap for Low Current Cases
K22 built a semianalytic model to estimate the peak curvature luminosity from the gap for SMBH cases with j 0 = −1/2π.In Section 3, we have seen the overall similarities of the gap dynamics for stellar-mass BH cases to SMBH cases, which motivate us to construct a semianalytic model applicable to stellar-mass BH cases.From the model of K22, we improve the model by taking into account general relativistic effects and the exact expression of the KN cross section.This model helps clarify the dynamics of gaps in BH magnetospheres and predict the gamma-ray luminosities and spectra from these gaps over broad ranges of mass M and accretion rate  M. We assume the electric field in the gap to be steady, since the oscillation periods in many cases are considerably longer than the light-crossing time of the gap.The nearly mono-energetic distribution of outgoing electrons around the peak time of electric field (see Figure 1) allows us to consider that all of the electrons would experience similar accelerations and radiative interactions.We solve the equation of motion and the scattered photon propagation for one representative electron, and determine the solutions for gap inner and outer boundaries r in and r out for which sufficient electron-positron pairs are created inside the gap.The corresponding maximum Lorentz factor e,max g and the gamma-ray peak luminosity L cur,pk are also calculated.Differences between our model and the previous analytic steady gap models (e.g., Hirotani et al. 2018) will be briefly discussed in Section 6.

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 f = 0 for all the particles, and the electron is rapidly accelerated to a highly relativistic velocity.Then we have dt g dr rr a » for the electron's motion, so that Equation (5) 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)) where we have imposed the boundary condition E r (r in ) = 0.The IC emitting power is written as Here σ ic is the KN cross section where The Lorentz factor of an electron at the kth step r k of the integration of Equation ( 11) is determined as  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 min  are most abundant when where l ic is the mean free path of IC scattering ic min e 0 rr s min ic We consider that all the scattered photons annihilate with the seed photons after propagating about the mean free path for the γγ pair annihilation l γγ where ò 2 is the energy of the annihilating seed photon, defined as max , } , and we use ) 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 where the radius of the "pair production front" r ppf of the scattered photons is given by The total number of newly created pairs from r in to position r is calculated as

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 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 which is satisfied in the simulation results, as shown in Figure 4. For the solution of r in and r out obtained above, we calculate the peak curvature luminosity by 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)).

Examples of Model Solutions
Figure 5 shows the solutions of our semianalytic model for τ 0 = 30 and 300, with the other parameters 10 min 6  = -, p = 2, and B H = 2π × 10 7 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 .
As can be seen in Figure 5, the maximum electron Lorentz factors 10 e, max 7 g ~for τ 0 = 30 and 2 10 e,max 6 g ~´for τ 0 = 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 τ 0 dependence of the IC-scattering mean free path; for τ 0 = 30, l ic becomes larger compared with the case of τ 0 = 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 7 .In this regime curvature cooling is more efficient than IC cooling because P cur e 4 g µ and the IC scattering is deeply in the KN regime.On the other hand, for τ 0 = 300, e,max g 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 τ 0 and L cur,pk  L ic,pk ∼ 10 −4 L BZ for τ 0  100 (see Figure 3).

Model-Simulation Comparison
We compare the simulation results and semianalytic model solutions of e,max g and L cur,pk /L BZ for j 0 = −1/2π in Figure 6.For models LA and LB, the semianalytic solutions show that e,max g and L cur,pk /L BZ are decreasing functions of τ 0 with a break around τ 0 ∼ 100-200, and well fit the simulation results particularly for τ 0 < 200.For models LC2 and LC3 we obtain ∼10 2 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 e,max g and L cur,pk /L BZ on τ 0 are similar to what we observe in the simulations.The breaks seen at τ 0  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 τ 0 dependences of the two cooling rates affect the dynamics and energetics of the gaps.

Model for j 0 = −1/2 Cases
In the simulations for j 0 = −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 e,max g and L cur,pk /L BZ in the simulation results by adding some modifications to the model for j 0 = −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 Note that the condition of Equation (24) for j 0 = −1/2π is consistent with Equation (27).We present the model solutions for j 0 = −1/2 (models I1-3) in Figure 6.As seen, the solutions well fit the simulation results.

Gamma-Ray Signals from Isolated Stellar-mass BHs
In this section, we estimate the gamma-ray luminosities from the magnetospheres of stellar-mass BHs in the Galaxy based on our semianalytic model.About 10 9 stellar-mass BHs are thought to exist in the Galaxy (e.Barkov et al. 2012;Ioka et al. 2017;Kimura et al. 2021b).The magnetic field at the BH vicinity in the case of rapidly spinning IBHs could be strong enough for their accretion disks to become MAD states (e.g., Kaaz et al. 2023).Then, the gamma-ray emission from the gaps in their magnetospheres could be detectable (Hirotani et al. 2016(Hirotani et al. , 2018;;Song et al. 2017).
When evaluating the input parameters of the semianalytic model, min  , τ 0 , and B H , we use the one-zone MAD model of Kimura et al. (2021b), in which  M is assumed to be the Bondi-Hoyle-Lyttleton accretion rate.We consider thermal synchrotron emission from the MAD act as the seed photons for IC scattering and γγ pair production in the gap.The magnetic field strength in the MAD B is derived as a function of M and the density of surrounding ISM n ISM for the fixed viscosity parameter α and plasma beta β.The electron temperature T e is also determined as a function of M and n ISM , by equating the electron heating rate in the MAD and the cooling rate via synchrotron radiation.Then, we can derive min where θ e = k B T e /m e c 2 denotes the normalized electron temperature, R r g  = is the size of MAD, and x M is a numerical factor (x M ∼ 25 for the synchrotron self-absorption thin limit).ν syn,pk and L syn,pk n are the synchrotron peak frequency and specific luminosity, respectively. 9We fix the spectral index of seed photons to p = 2, since the dependence of the gap dynamics on p is minor (see Section 3.1) and the spectral slope of the thermal component around the peak does not deviate much from We do not consider reduction of the mass accretion rate due to the disk wind for simplicity.The magnetic field strength at the BH horizon B H is estimated by Equation (7).
In Figure 7 we present the resultant color maps of L cur,pk /L BZ calculated for 3.2 M e M 10 3 M e and 1.0 cm −3 n ISM 10 4 cm −3 in the cases of j 0 = −1/2π and −1/2.As can be seen, L cur,pk /L BZ is higher in the parameter regions of smaller M and smaller n ISM for both cases of j 0 .(Note that L BZ has the opposite trend.)Combining the results of the two cases of j 0 , we find L cur,pk /L BZ  10 −4 for a limited range of n ISM : for M ∼ 10 M e , the range is 50 cm −3  n ISM  3.0 × 10 2 cm −3 , which corresponds to dense molecular clouds, while for M ∼ 50 M e , the range is 5 cm −3  n ISM  30 cm −3 , which corresponds to a cold neutral medium.The parameters for the white regions in Figure 7 yield τ 0  27 for the j 0 = −1/2π case and τ 0  14 for the j 0 = −1/2 case, for which we could not find the solution of the gap in r > r in .The gap would achieve its maximum efficiency of particle acceleration when τ 0 is just above this critical value, and correspondingly, the curvature luminosity is also the highest.These critical values of τ 0 are consistent with our simulation Figure 7.The color maps of L cur,pk /L BZ as functions of M and n ISM for j 0 = −1/2π (left) and j 0 = −1/2 (right).The white region indicates the parameter ranges that have no semianalytic solutions.The stars indicate the parameter sets which we use in the spectral calculations for Figures 8 and 9.We also plot contours of L BZ with dotted lines and B H = 10 7 G with a dotted-dashed line. 9Here and hereafter we use the convention Q Q 10 x x = ( ) in cgs units, except for the BH mass, M = 10 x M x M e .results: the simulation box slowly approaches the vacuum state for models with τ 0 = 10 for all three j 0 values (though not shown in Table 1), similar to what we saw in the SMBH cases. 10e present examples of broadband emission spectra from IBH MADs with magnetospheric gaps based on our semianalytic calculations in Figure 8.The parameters used for these examples are marked by the black stars in Figure 7.The spectral components from the magnetospheric gaps and the MADs, taking into account the internal attenuation by γγ pair annihilation of gamma-rays and MAD photons, are calculated by the methods in Appendix B and Kimura et al. (2021b), respectively.One thing to note is that we consider timeaveraged luminosities when calculating the gap emission: we observe that L cur and L ic oscillate in the simulations (Figure 2).The oscillation period T gap  10r g /c ; 5.0 × 10 −4 M 1 s, however, cannot be resolved by ongoing or forthcoming gammaray detectors for stellar-mass BHs.Hence, we calculate the observed spectral luminosities by reducing the intrinsic peak values by the typical duty cycle seen in the simulations, 60%.As seen in Figure 8, the gap emission peaks in the GeV-TeV energy band for both j 0 = −1/2π and −1/2.
The spectral cutoff at the sub-TeV energy band is the result of pair annihilation with synchrotron photons from the MAD thermal electrons (see Appendix B).For a given set of M and n ISM , the spectral luminosity in the GeV-TeV energy band is higher for j 0 = −1/2π, owing to the efficient curvature radiation.For j 0 = −1/2, the curvature peak luminosity and the characteristic frequency are 10 2 -10 3 times lower than those for j 0 = −1/2π, but IC emission is dominant in the GeV-TeV energy band due to the scaling L ic,pk ∝ |j 0 | (see Section 3.2).As a result, the gamma-ray signals from IBH magnetospheric gaps could be bright as L  10 −4 L BZ over wide ranges of M and n ISM .
For an illustrative purpose, we show the gamma-ray fluxes from the gaps of a nearby (d L = 0.5 kpc) IBH of M = 50 M e for various values of n ISM in Figure 9.The peak fluxes exceeds E F 10 erg s cm , which is high enough for detection with current gamma-ray detectors, including Fermi's Large Area Telescope (LAT; Atwood et al. 2009).

Summary and Discussion
We have conducted a comprehensive investigation of the dynamics of spark gaps in the charge-starved magnetospheres of stellar-mass Kerr BHs, using the 1D GRPIC code described in Levinson & Cerutti (2018), K20, and K22.The differences of dynamics for changing spectral properties of seed photons and the magnetospheric current have been examined.We have also presented a semianalytic model that reproduces the maximum electron Lorentz factors and peak gamma-ray The black and gray solid lines represent the total emission for j 0 = −1/2 and j 0 = −1/2π, respectively.The red dashed line represents radiation from the MAD, which consists of synchrotron radiation from MAD thermal electrons (which form the humps in the optical band) and that from MAD nonthermal electrons.The spectra of the curvature and IC radiation from the gap are shown by the blue and light blue lines, respectively.For each of them, the dashed lines correspond to those for the j 0 = −1/2 case and the dotted ones for the j 0 = −1/2π case.luminosities of the simulation results for various parameter values.Using this model, we have estimated the characteristics of gamma-ray signals from IBHs over broad ranges of BH masses and infalling ISM densities.The gamma-ray signals from the gaps of IBHs embedded in dense gas clouds located at distances less than a kiloparsec are suggested to be bright enough for detections in the GeV energy band.

Additional Pair Production Processes in the Magnetosphere
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 -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 r B is the magnetic field strength around the null charge surface.Furthermore, n ± /n GJ and τ 0 have the same dependence, M n 2 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) 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 5 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 3 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)  .B(r null ) ; 5 × 10 6 B H,7 G for M ∼ 10 M e and n ISM ∼ 10 2 cm −3 , 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 7 , and thus the photons produced in the gap via IC scattering will be ò γ ∼ 3.0 × 10 6 (see Equation ( 15)). 11The corresponding optical depth of the magnetic pair production is τ B ∼ α B r g ; 4 × 10 3 for sin 1.0 y = , and 2 × 10 −2 for sin 0.4 y = .In the latter estimate we have approximated ψ by the field pitch angle |B f |/B r ∼ 0.4 measured in the ZAMO frame (see Appendix A in Kimura et al. 2022).Such a 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  n 30 0 ISM 3 2 t µ , ò γ  10 6 and the gap width is much smaller than r g .

Possible Strategy for Candidate Search
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 E µ g g -G (Γ ∼ 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 whitedwarf-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 L cur,pk 0 ic,pk BZ µ µ , hence occasional luminosity variations by a factor of 10 2 (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, ) , where v ; 40 km s −1 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(Abdollahi et al. , 2020) ) if the source is bright.et al. (2016, 2017), Song et al. (2017), andHirotani 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).

Hirotani
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 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 2 10 ) , thus they will have a large mean free path for pair annihilation l R ) .Since sin 1 a < , 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 0 (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.
where 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 .

B.2. Inverse Compton Emission
For τ 0  10, for which the magnetosphere does not approach the vacuum, the primary IC photons produced by electrons with e e , m a x g g » 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 t t cool ic dyn cool cur outside the gap, where t cool ic , t cool 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 constant ic,pk ~outside the gap.The simulation results indicate L ic,pk ∼ 10 −3 |j 0 |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 is a kernel in the relativistic limit (Blumenthal & Gould 1970;Rybicki & Lightman 1986) and F x = −1/2π, ∼0.3 for j 0 = −1/2, and ∼0.5 for j 0 = −1,

Figure 1 .
Figure 1.Snapshots of the gap at three times in one period of the oscillation in model LA1.The three rows of panels from top to bottom for each of the times are (top) the number density of each particle species normalized by the GJ number density at the inner boundary n e GJ,min GJ,min r = ( r min min D = D( )), (middle) the electric field distribution normalized by B H , and (bottom) the energy distribution of outgoing particles of each species.
Figure 5), where γ e,acc,k is the value for the acceleration-dominant regime cur,k and γ e,ic,k are the values determined by the balances of the acceleration with the curvature and IC coolings, respectively,

Figure 2 .
Figure 2. Luminosities of curvature emission (black line), IC emission (red line), and kinetic energy (blue line) as ratios to the BZ lumionsity L BZ .The vertical dashed lines correspond to the times of the snapshots in Figure 1.

Figure 3 .
Figure 3. Plots of the peak curvature luminosity normalized by L BZ (left) in our simulations for various values of p, min  , and j 0 as a function of τ 0 , and the peak IC luminosity normalized by L BZ (right) in our simulations with 10 min 6  = -and p = 2.0 as a function of |j 0 |.

Figure 4 .
Figure 4. Middle points of the gaps in the ξ coordinate from the simulation results.The horizontal dashed line represents the position of the null charge surface.

Figure 5 .
Figure 5. Results of our semianalytic model for different τ 0 (left: τ 0 = 30, right: τ 0 = 300), with p 10 , 2.0, min 6  = = and B H = 2π × 10 7 G.The top panels show the evolution of γ e,acc (orange dashed line), γ e,cur (blue dashed line), and γ e,ic (light-green dashed line).The red solid curve represents the solution of γ e .The black solid line represents the radius at which the scattered photons annihilate with the seed photons, r ppf .The bottom panel show the evolution of the number of pairs N pair .The vertical dashed lines indicate r in (where N pair = 0) and r out (where N pair = 2.0), and the white region corresponds to the gap.

Figure 6 .
Figure 6.Comparison between the simulation results (points) and the semianalytic model solutions (lines).We plot e,max g (left) and L cur /L BZ (right) as functions of τ 0 .

Figure 8 .
Figure 8. Broadband emission spectra from IBH MADs with magnetospheric gaps in the cases of M = 10 M e and n ISM = 100 cm −3 (left panel), and M = 50 M e and n ISM = 10 cm −3 (right panel).The valued of τ 0 and    m M M Edd = are shown in each panel.The black and gray solid lines represent the total emission for j 0 = −1/2 and j 0 = −1/2π, respectively.The red dashed line represents radiation from the MAD, which consists of synchrotron radiation from MAD thermal electrons (which form the humps in the optical band) and that from MAD nonthermal electrons.The spectra of the curvature and IC radiation from the gap are shown by the blue and light blue lines, respectively.For each of them, the dashed lines correspond to those for the j 0 = −1/2 case and the dotted ones for the j 0 = −1/2π case.
α f = e 2 /ÿc ; 1/137 is the fine structure constant, c l -= ÿ/m e c ; 3.9 × 10 −11 cm is the reduced Compton wavelength, b = B/B q is the local magnetic field strength normalized by /

Table 1
Model Parameters of the Simulation x , and the accelerated electrons are efficiently cooled by emitting MeV photons.Then the luminosity is Timokhin & Harding 2015;Kisaka & Tanaka 2017; K20).The Lorentz factor of electrons that produce IC photons of energy ò ic, esc is evaluated