{'fraction_non_alphanumeric': 0.03370577281191806, 'fraction_numerical': 0.020074487895716947, 'mean_word_length': 4.373424054432659, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 1, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
The response of black hole spark gaps to external changes: A production mechanism of rapid TeV flares?

We study the response of a starved Kerr black hole magnetosphere to abrupt changes in the intensity of disk emission and in the global magnetospheric current, by means of 1D general relativistic particle-in-cell simulations. Such changes likely arise from the intermittency of the accretion process. We find that in cases where the pair production opacity contributed by the soft disk photons is modest, as in, e.g., M87, such changes can give rise to delayed, strong TeV flares, dominated by curvature emission of particles accelerated in the gap. The flare rise time, and the delay between the external variation and the onset of the flare emitted from the outer gap boundary, are of the order of the light crossing time of the gap. The rapid, large amplitude TeV flares observed in M87 and, perhaps, other AGNs may be produced by such a mechanism.


INTRODUCTION
Many active galactic nuclei (AGNs) occasionally exhibit rapidly variable TeV emission, the origin of which is yet unknown. A notable example is the large amplitude, rapid TeV flares observed in M87 (Aharonian et al. 2006;Abramowski et al. 2012). The short duration of these extreme flares, roughly the dynamical time of the central black hole (BH), as well as other indicators (e.g., coincidence with core radio emission) suggest that the gamma-ray emission originates from the very inner regions of the magnetosphere.
It has been proposed that charge starved magnetospheric regions (spark gaps) might provide sites for variable TeV emission (Levinson 2000;Levinson & Rieger 2011;Hirotani et al. 2016Hirotani et al. , 2017Levinson & Segev 2017). Such spark gaps are expected to form when direct plasma supply to the inner region of an active magnetosphere, where jets are formed, is too low to screen out the magnetosphere. Particle acceleration by the gap electric field then gives rise to copious pair creation through the interaction of particles accelerated in the gap with soft photons emitted by the accretion flow. Recent one-dimensional (1D) general relativistic Particle-in-Cell (GRPIC) simulations of local spark gaps (Levinson & Cerutti 2018;Chen & Yuan 2020;Kisaka et al. 2020) and global 2D GRPIC simulations of active Kerr BH magnetospheres (Parfrey et al. 2019;Crinquand et al. 2020Crinquand et al. , 2021 indicate that the spark process is highly intermittent, conceivably inherently cyclic under certain conditions. However, the amplitude of the variations and the overall luminosity were found to depend rather sensitively on the pair creation opacity contributed by the soft radiation emitted from the surrounding accretion flow. These simulations invoke steady external conditions, particularly the soft photon intensity that is given as input. In reality, it is likely that the accretion flow will undergo kisaka@hiroshima-u.ac.jp 1 Department of Physics, Hiroshima University, Higashi-Hiroshima, 739-8526, Japan 2 School of Physics & Astronomy, Tel Aviv University, Tel Aviv 69978, Israel 3 Frontier Research Institute for Interdisciplinary Sciences, Tohoku University, Sendai, 980-8578, Japan 4 Astronomical Institute, Graduate School of Science, Tohoku University, Sendai, 980-8578, Japan temporal variations that induce intermittencies in the pair production opacity and the strength (and geometry) of the magnetic field advected by the accretion flow. It is anticipated that such external changes might lead to a nonlinear response of the gap and its TeV emission, e.g., by the development of the electric field due to the change of the optical depth and the electric current density, and the screening of that electric field by the created pairs thus created. The variations in the gap emission should be delayed with respect to the variations of the target photon intensity. In principle, it might be possible to map such reverberations and extract information about the structure of the magnetosphere and the spark process, if the delays can be quantified and general relativistic time travel effects are taken into account.
In this paper we study the response of a local spark gap to external changes using 1D GRPIC simulations. We explore the gap dynamics induced by the change of different parameters, particularly target photon spectrum and luminosity, and the global magnetospheric current. We find that the gap response is highly non-linear by virtue of Klein-Nishina (KN) effects and the sensitive dependence of the curvature radiation power on the energy of emitting particles, and that such changes can lead to production of strong flares of curvature TeV emission by the gap accelerating particles.

THE GAP MODEL
Our analysis is restricted to the response of a 1D, local spark gap. The tacit assumption underlying the local gap model is that its activity does not affect significantly the global magnetospheric structure. The global electric current and the characteristics of soft-photon intensity are treated as input parameters in this model. The calculations are performed using the 1D GRPIC code described in Levinson & Cerutti (2018) and Kisaka et al. (2020). The dynamical equations are solved in a background Kerr spacetime, given in Boyer-Lindquist coordinates with the radial coordinate replaced by the tortoise coordinate ξ to avoid the singularity on the horizon (see Levinson & Cerutti 2018 for details).
The computational domain extends from a radius of ∼ 1.5r g to ∼ 4.3r g , where r g = GMc −2 = 1.5 × 10 14 M 9 cm is the gravitational radius, G is the Gravitational constant, c is light speed, and M ≡ 10 9 M ⊙ M 9 is a BH mass. For the results pre-sented below we used a BH mass of M = 10 9 M ⊙ and a dimensionless spin parameter a * = 0.9. The corresponding angular velocity of the BH is ω H = a * c/2r H ≈ 10 −4 M −1 9 Hz, where r H = r g (1 + 1 − a 2 * ) is the horizon radius. For the angular velocity of the magnetic surface along which the gap is located we adopt the value Ω = ω H /2, such that the null charge surface, on which the Goldreich-Julian (GJ) charge density vanishes, is located at r null = 2r g . For all cases considered below we use a split monopole magnetic field configuration with magnetic field strength B H /2π = 10 3 G on the horizon, and adopt an inclination angle θ = 30 • (with respect to the BH rotation axis) for the flux surface that contains the gap. An open boundary condition is imposed on both sides (boundaries) of the simulation box. The number of grid cells used in a typical simulation is 32768, sufficient to resolve the skin depth, and the initial particle per cell number (PPC) is 45, which guarantees convergence (see Kisaka et al. 2020 for details). Finally, the time step is set by the Courant-Friedrichs-Lewy condition at the inner boundary.

Emission and pair creation
Two processes dominate the emission by particles accelerated in the gap, inverse Compton (IC) scattering and curvature emission. In essentially all cases considered here, the mean energy of curvature photons is considerably smaller than that of IC scattered photons, rendering pair creation by curvature photons negligible, even in cases where the gamma ray luminosity is dominated by curvature emission (see Kisaka et al. 2020 for detailed estimates). We therefore treat these two processes in a distinct manner. The IC scattered photons are included as charge-neutral particles in the simulation and their momenta are computed by solving the appropriate equations of motions. The curvature emission, on the other hand, is computed at each time step using analytic formula for single particle emission, and summing the contributions from all particles in the simulation box, taking into account time travel effects (Levinson & Cerutti 2018;Kisaka et al. 2020). Pair creation is contributed solely by the interaction of IC photons with the background soft photons, envisioned to emerge from the surrounding accretion flow.
IC scattering and pair creation are computed using Monte-Carlo methods that are incorporated in the PIC code, as described in Levinson & Cerutti (2018). The full KN cross section is used for both processes. The background photon intensity is taken, for simplicity, to be isotropic and homogeneous, with a power-law spectrum: where ε is the photon energy normalized by an electron rest mass energy, m e c 2 . It is characterized by three key parameters: the minimum cutoff energy ε min , the slope p, and the normalization I 0 that we, alternatively, measure by the fiducial optical depth here σ T is the Thomson cross section and h is the Planck constant. As in Kisaka et al. (2020), to reduce the computational cost, we ignore the IC scattering for the particles with γ < 10 7 , since at these energies the created photons are below the pair production threshold.

Dichotomic gap behaviour
In Kisaka et al. (2020) we identified a dichotomy in the gap dynamics that depends on the sign of the magnetospheric current j 0 , henceforth normalized by [ω H B(1 + a 2 * ) cos θ /2π]. From this definition, the current carried by the outflowing plasma particles with GJ charge density at large distance is j 0 = −1.0. For negative values, j 0 < 0, the gap exhibits a cyclic discharge, with superposed plasma oscillations, confined to the region around the null surface (Kisaka et al. 2020). The duty cycle and amplitude of these oscillations depend quite sensitively on the fiducial optical depth τ 0 and minimum energy ε min . For positive currents we observed mainly rapid oscillations, with occasional opening of the gap at the boundaries in certain cases.
Despite the inherent intermittency of the discharge dynamics, the long-term behaviour has been found to be quasiregular. This is not surprising given that all key parameters were fixed to constant values. However, the high nonlinearity of the gap dynamics suggests that any temporal changes in any of these parameters should lead to a more chaotic, vigorous response. This is the main motivation for the analysis described in the following section. Since negative current is anticipated in the polar jet section we restrict our calculations to this case.

Simplistic modelling of external changes
Recent GR magnetohydrodynamic (MHD) simulations (e.g., Porth et al. 2019;Ripperda et al. 2020;Dexter et al. 2020;Nathanail et al. 2020;Chashkina et al. 2021) indicate a complex accretion dynamics that depends on the initial setup, particularly the initial magnetic field configuration. At sufficiently high resolution, drastic changes in the accretion rate and the structure of the inner magnetosphere have been observed during states of magnetically arrested disk, even for a simple initial magnetic field configuration. The flaring behaviour appears to be regulated by a reconnecting equatorial current sheet (Ripperda et al. 2020;Chashkina et al. 2021), and if small scale field is initially present also by sporadic reconnection at the boundaries of interacting flux tubes with opposite polarities (Nathanail et al. 2020;Chashkina et al. 2021). Formation of an equatorial current sheet with a feedback on the global magnetosphere has been observed also in recent global GRPIC simulations of active black holes (Crinquand et al. 2021). It is quite likely, therefore, that this complex behaviour will give rise to changes in disk emission and magnetospheric currents over dynamical timescales.
The radio observations indicate that the spectral peak flux of M87 would come from the inner region, r < 7r g (Fig 16 of EHT MWL Science Working Group et al. (2021)). The peak frequency is consistent with ∼ 10 −9 m e c 2 used as ε min in our model. Then, the soft photon field may be sensitive to the accretion rate and have a short variability, 10r g /c.
Within the framework of our simplistic model this can be captured by allowing temporal changes in the model parameters, j 0 , ε min , p, τ 0 . Our strategy is to study the gap response to abrupt changes in these parameters. For our reference model we invoke the parameters τ 0 = 100, ε min = 10 −9 , p = 2, and j 0 = −1.0. In each numerical experiment we change, instantaneously, one of those parameters relative to the reference model, at a transition time that we define as t = 0 (except for Fig. 1). The various models are listed in Table 1.   We consider first the gap response to changes in the disk emission for a constant magnetospheric current j 0 = −1. All luminosities are normalized by the Blandford-Znajek (BZ) power (Blandford & Znajek 1977), for the set of parameters (a * , r g , B H ) adopted in these simulations. We generally find that relatively modest changes in the parameters of the target photon intensity lead to a highly nonlinear response of the gap, that appears as large amplitude flares of curvature emission, in the regime of mild opacities (τ 0 300). The IC emission is far less sensitive. This is a consequence of the sensitivity of the curvature power to the Lorentz factor of emitting pairs (see appendix A for a detailed explanation), and the fact that IC scattering is in the deep KN regime.
To illustrate a typical time sequence of events, we show in Fig. 1 an example where the fiducial optical depth has changed from τ 0 = 100 to τ 0 = 30 (at t = 0) for a period of ∆t = 10 t g (left) and ∆t = 1t g (right), where t g ≡ r g /c, and then back to its initial value, keeping all other parameters fixed. As seen, a flare dominated by curvature radiation is emitted from the outer boundary of the simulation box, with a delay of about 7 t g in both cases, roughly the propagation time of photons from the gap to the boundary (including redshift effects). The detailed shape of the flare is clearly different in each case. In both cases the curvature power rises sharply initially by a factor of ∼ 10, with a rise time of about 1t g . In the case on the left there is a further, slower increase by a similar factor over time of about 5t g , whereas the sharp increase continues in the other case (since the opacity by this time has already switched back to τ 0 = 100, the number of emitting particle increases due to pair creation). The flare rise seen in both cases results from excessive emission by pairs accelerating in the growing gap electric field following the sudden decline in opacity at t = 0. The opening of the gap near the null surface commences almost immediately, followed by a growth of the electric field over several t g due to the decline in the density of photons and pairs, owing to the reduced opacity. This timescale determines the rise time of the flare, whereas the delay between the sudden change in the opacity and the flare onset is determined by the propagation time of the curvature photons from the gap to the outer boundary. In the case on left there is a subsequent plateau phase which we think indicates the saturation of the pair cascade before the opacity is switched back. This is not seen in the other case, presumably because the cycle was too short for saturation to occur before the opacity is switched back. Ultimately, there is a sharp drop by a large factor over a time of about 1 − 2t g , causing a large undershoot followed by a fast recovery to the pre-flare luminosity.
We find similar behaviour in other cases. The conclusion is that for abrupt changes the flare rise time is comparable to the gap size, although the exact shape of the light curve depends on details. We stress that our calculations do not take into account multi-dimensional as well as general relativistic effects on the trajectories of the emitted photons. In reality, if the gap is located close enough to the horizon, these effects might dominate. However, they can be computed accurately using ray tracing codes (e.g., Crinquand et al. 2021), so that the intrinsic timescales can be extracted. It is also worth noting that gradual changes of the external parameters (over times longer the the gap response time) are likely to smear out the intrinsic response time, but a delay is, nonetheless, expected. The extreme TeV flares observed in M87 (Aharonian et al. 2006;Abramowski et al. 2012) and IC 310 (Aleksić et al. 2014), if produced by such a process, imply gap size of ∼ r g , and changes in external parameters over time scales not much longer than the light crossing time of the gap (see right panel in Fig. 1). Figure 2 shows mosaic of lightcurves obtained for models A−F in table 1. The black, red, and blue lines correspond, respectively, to fluxes of curvature radiation, IC emission and particles kinetic energy at the outer boundary of the computation box. In these runs, we start with the fiducial model and then, at time t = 0 in the plots, change one of the parameters abruptly and let the system relax to the new state. As in Fig. 1, we observe a delay of about 7t g in the gap response in all cases. After the delay, the activity quickly (∼ t g ) reaches another quasi-periodic state for the new parameter set, which is consistent, in all cases, with the results reported in Kisaka et al. (2020). As seen, strong flares are produced in episodes in which the opacity effectively declines. This applies to the cases displayed in the three upper panels. Steepening of the spectrum (increased p) results in a reduced optical depth due to reduction in the number of photons above the peak (i.e., above ε min ) that interact with lower energy pairs. As for the change in ε min , the response arises since the particle kinetic energy flux increases following the transition. In the fiducial case with ε min = 10 −9 , particles lose their energy near the outer boundary by virtue of IC scattering in the Thomson regime. The typical energy of particles near the outer boundary is γ ∼ 10 8 . Then, after switching to ε min = 10 −8 , the same particles scatter in the KN regime, thereby cool much slower. As a result, their mean energy increase leading to significant enhancement in curvature emission.
In general, we find that the highly nonlinear response we observed in many runs is a combination of KN effects and the sensitive dependence of the curvature emission power on the energy of emitting particles. The dependence of the peak luminosity of curvature emission on the optical depth τ 0 is delineated by the black line in Figure 3. The red line corresponds to the analytic result derived in appendix A. This sensitivity of L cur to changes in opacity is the key factor that leads to the flaring behaviour seen in Figs. 1 and 2.

Change of magnetospheric current
In the second set of numerical experiments, models G−L in table 1, we keep the soft-photon intensity fixed and change the global current j 0 at time t = 0. Since the net electric current in the dynamical equation that governs the evolution of the gap electric field is ( j e − j 0 ), where j e is the current carried by pairs produce in the gap, changes in j 0 lead to a change in the net current and, consequently, an inductive response. For example, since the number density of particles in the gap is comparable to the GJ value times | j 0 |, when | j 0 | increases (as in case G, where j 0 = −1.0 → −2.0), the charge density in the gap is insufficient for compensating the current difference j e − j 0 . As a result, the flux of curvature photons temporarily increases, as seen in the left panels in Fig. 4. Figure 4 shows the resultant light curves. The upper three panels correspond to models G−I in table 1 (τ 0 = 100) and the lower panels to Models J−L (τ 0 = 30). In all cases shown the global current changes from j 0 = −1 to the value indicated in the top label. As seen, the curvature emission flux appears to be quite sensitive to j 0 for τ 0 = 100, less so for the lower optical depth run. Nonetheless, an increase by a factor of about 3 in the flux following the transition is observed also at τ 0 = 30 (lower panels in Fig. 4). The reason for that trend is not entirely clear to us. In the run with τ 0 = 100 we observed significantly larger amplitudes of the electric field oscillations at lower values of j 0 , whereas for τ 0 = 30 these oscillations appears much less sensitive to j 0 . We suspect that KN effects might play a role in dictating the dependence of gap opening on the global magnetospheric current.

DISCUSSION
We studied the response of a black hole spark gap to external changes by means of 1D GRPIC simulations. We focused on changes in the global magnetospheric current, and in the intensity of soft disk photons that dominate the pair-production and inverse-Compton opacities. Such changes likely arise from the intermittency of the accretion process, as observed in recent GRMHD simulations (e.g., Porth et al. 2019;Dexter et al. 2020;Ripperda et al. 2020;Nathanail et al. 2020;Porth et al. 2021;Chashkina et al. 2021), and conceivably the presence of small scale magnetic fields in the accretion flow (Nathanail et al. 2020;Chashkina et al. 2021). Our analysis is restricted to local 1D gaps and neglects the feedback on the global magnetosphere.
We generally find that in the regime of modest pairproduction opacity, as anticipated in, e.g., M87 and possibly other AGNs, sudden changes in the pair-production optical depth inside the gap give rise to strong flares of curvature TeV photons, with an intrinsic rise time of the order of the gap light crossing time and a similar delay between the variation of the external radiation intensity and the onset of the flare. Changes in the global magnetospheric current can also lead to flaring episodes under appropriate conditions. The highly nonlinear response of the gap stems from KN effects and the sensitive dependence of the curvature emission power on the energy of emitting particles.
Our analysis ignores time travel and general relativistic effects. In reality, if the gap is located close to the horizon, one expects strong lensing to significantly affect the lightcurve. The extreme TeV flares observed in M87 (Aharonian et al. 2006;Abramowski et al. 2012) and IC 310 (Aleksić et al. 2014), if produced by such a process, imply gap size of ∼ r g , and similar changes in disk emission. Alternatively, those flares can be produced by reconnection episodes in an equatorial current sheet near the horizon, hadronic process in the accretion disk (Kimura & Toma 2020), or cloud-jet interaction (Barkov et al. 2012).
We are grateful to the anonymous referee for constructive comments. Numerical computations were performed on Cray XC50 at cfca of National Astronomical Observatory of Japan, and on Cray XC40 and Yukawa-21 at Yukawa Institute for Theoretical Physics, Kyoto University. This work was supported by JSPS Grants-in-Aid for Scientific Research 18H01245 (K.T.), 18H01246, 19K14712 and 21H01078 (S.K.).

APPENDIX ANALYTICAL ESTIMATE OF THE CURVATURE LUMINOSITY
To elucidate the sensitive dependence of the luminosity of curvature radiation on the model parameters, we calculate the luminosity using a simple analytic model of the gap. The analytic results are compared with the curvature luminosity obtained from the simulations (Fig. 3). The analytic model assumes that the curvature luminosity emitted from the outer gap boundary is dominated by emission of particles accelerated in the open gap region around the null charge surface, where the electric field E || (r) is strongest and oscillations can be neglected. We ignore any velocity spread, as it is negligibly small, and suppose that at any given radius r the energy of an emitting particle, γ(r), is limited by curvature losses (as seen in Figs. 1,2 and 4, IC losses are negligible in most cases of interest). It is obtained by equating the curvature power of a single particle, Here e is a charge of an electron and R c is the curvature radius of the magnetic field line confining the accelerating particles, assumed to be constant throughout the gap. The total curvature luminosity is obtained by integrating P cur over all the particles inside the gap: where r out and r in are the outer and inner gap boundaries, respectively, α(r) is the lapse function, and n is the number density of emitting particles, taken to be equal to | j 0 |ω H B(1 + a 2 * ) cos θ /2πec at r = r null . The electric field is obtained from Gauss' law, noting that inside the gap the net charge density of the counter-streaming electrons and positrons is well below the GJ density, |ρ e | ≪ |ρ GJ |. With ρ GJ (r) ∝ [2(r/r H ) −2.5 − 1] used as an approximate fit to the GJ charge density around the null charge surface, we thus have The radius of the outer gap boundary, r out , is treated as a free parameter. For a given r out , we calculate the density of newly created pairs, as explained below, and use it to find the inner boundary of the gap, r in , which is taken to be at the radius at which the number density of created pairs is equal to the number density of particles injected into the gap. We then iteratively seek the radii r out and r in that satisfy (r out + r in )/2 = r null .
To estimate the density of created pairs, we only consider the annihilation of IC gamma-ray photons with soft disk photons. The number density of scattered IC photons in an interval ∆r is where l IC is the mean free path for IC scattering. Since the scattering is not in the deep KN regime, we simply use the following approximation for the mean free path: l IC ∼ r g /τ 0 (γε min ≤ 1) r g ε min γ/τ 0 (γε min > 1). (A6) The energy of a scattered photon satisfies ε IC ∼ 0.1γ 2 ε min (γε min ≤ 1) 0.1γ (γε min > 1).
We assume that a scattered photon is annihilated after propagating a distance comparable to the characteristic mean free path for the pair creation, where ε tar is the minimum energy of target photons, defined as ε tar = max{ε min , ε −1 IC }.
Using these estimates we integrate the transfer equation to obtain the number density of newly created pairs, and use it to determine the inner boundary of the gap, r in , assumed to be located where the integrated number density reaches the GJ value. As a rough estimate, using equations (A1) and (A3), the curvature luminosity is L cur ∝ γ 4 r gap , where r gap = r out − r in is the gap width. Using equations (A2) and (A4), the Lorentz factor of the particles is γ ∝ E 1/4 ∝ r 1/4 gap . The gap width is roughly the minimum of the sum of two mean free path, l IC and l γγ . In the KN regime, both mean free path scales as γ/τ 0 from equations (A6) and (A8). Then, the curvature luminosity scales as L cur ∝ r 2 gap ∝ τ −8/3 0 . Figure 3 shows the calculated curvature luminosity as a function of the optical depth τ 0 (red dashed curve), for j 0 = −1.0, ε min = 10 −9 , and p = 2. The analytic result is compared with the luminosity obtained from the numerical simulations. As seen, the agreement is extremely good. The curvature luminosity in the analytic result scales as L cur ∝ τ −3 0 , roughly consistent with the above rough estimate.