Particle-in-cell simulations of pair discharges in a starved magnetosphere of a Kerr black hole

We investigate the dynamics and emission of a starved magnetospheric region (gap) formed in the vicinity of a Kerr black hole horizon, using a new, fully general relativistic particle-in-cell code that implements Monte Carlo methods to compute gamma-ray emission and pair production through the interaction of pairs and gamma rays with soft photons emitted by the accretion flow. It is found that when the Thomson length for collision with disk photons exceeds the gap width, screening of the gap occurs through low-amplitude, rapid plasma oscillations that produce self-sustained pair cascades, with quasi-stationary pair and gamma-ray spectra, and with a pair multiplicity that increases in proportion to the pair production opacity. The gamma-ray spectrum emitted from the gap peaks in the TeV band, with a total luminosity that constitutes a fraction of about $10^{-5}$ of the corresponding Blandford-Znajek power. This stage is preceded by a prompt discharge phase of duration $\sim r_g/c$, during which the potential energy initially stored in the gap is released as a flare of curvature TeV photons. We speculate that the TeV emission observed in M87 may be produced by pair discharges in a spark gap.


Introduction
The activation of Blandford-Znajek (BZ) outflows requires continuous injection of plasma in the magnetospheric region enclosed between the inner and outer light cylinders. The origin of this plasma source is still an open issue. To fully screen out the magnetosphere, the plasma injection rate must be sufficiently high to maintain the density everywhere in the magnetosphere above the Goldreich-Julian (GJ) value. If the plasma source cannot accommodate this requirement, charge starved regions will be created, potentially leading to self-sustained pair discharges. A plausible plasma production mechanism in black hole (BH) engines is the annihilation of gamma rays that emanate from the accretion flow (e.g., Blandford & Znajek 1977;Levinson 2000;Neronov & Aharonian 2007;Levinson & Rieger 2011;Levinson & Segev 2017) or neutrino annihilation in the case of GRBs (Globus & Levinson 2014). Whether the pair density thereby produced is sufficiently high depends primarily on the luminosity of MeV photons emitted by the radiative inefficient accretion flow (Levinson & Rieger 2011;Katsoulakos & Rieger 2018) or by a putative corona in the case of sources accreting at a rate in excess of the critical ADAF rate (see, e.g., Levinson & Segev 2017).
It has been shown (Levinson & Rieger 2011;) that under conditions anticipated in many stellar and supermassive BH systems, the annihilation rate of disk photons is insufficient to maintain the charge density in the magnetosphere at the GJ value, giving rise to the formation of spark gaps. It has been further pointed out (Neronov & Aharonian 2007;Rieger 2011;Levinson & Rieger 2011;Lin et al. 2017;Katsoulakos & Rieger 2018) that the gap activity may be imprinted in the high-energy emission observed in these sources whereby the variable TeV emission detected in M87 (Aharonian et al. 2003;Acciari et al. 2009) and IC310 (Aleksić et al. 2014) was speculated to constitute examples of the signature of magnetospheric plasma production on horizon scales (Levinson 2000;Neronov & Aharonian 2007;Levinson & Rieger 2011;. In an attempt to compute the emission properties of an active gap, a fully general relativistic (GR) steady gap model has recently been developed Levinson & Segev 2017) that incorporates curvature emission, inverse Compton scattering, and pair creation via the interaction of gamma rays produced in the gap with the external photons emanating from the accretion disk. However, as argued in Levinson & Segev (2017), steady-state solutions are restricted to a narrow range of conditions that may not apply to most systems. Furthermore, it is unclear whether the steady gap solutions obtained in the works cited above are stable in the first place.
In this paper we explore the dynamics of a local, 1D gap using a new particle-in-cell (PIC) code, developed particularly for this purpose, that implements Monte Carlo methods to compute gamma-ray emission and pair production through the interaction of pairs and gamma rays with an external radiation field. The details are described in §2 and §3 below. We find ( §4) that after a prompt discharge phase that screens out the gap and produces a strong gamma-ray flare, the system relaxes to a state of self-sustained, rapid plasma oscillations that is independent of the initial conditions. The amplitude of the oscillating electric field in the relaxed state is regulated by pair production, such that the average pair creation rate inside the simulation box equals the rate at which pairs escape through the boundaries. The pair and photon spectra are quasi-stationary during this state. We find that the saturated pair multiplicity increases roughly linearly Article number, page 1 of 14 arXiv:1803.04427v2 [astro-ph.HE] 6 Jun 2018 A&A proofs: manuscript no. gap_pic_pap-AA with the opacity contributed by the ambient radiation field, while the average pair and photon energies decrease as the opacity increases. The gamma-ray luminosity radiated by the accelerating pairs following the prompt phase (after relaxation of the system) constitutes a fraction of about 10 −5 of the corresponding BZ power, weakly dependent on the pair production opacity and other model parameters. As pointed out in §5, this is in rough agreement with the TeV observations of M87.
Our approach is similar to that used by Timokhin (2010) and Timokhin & Arons (2013), with the following differences: it is fully GR, it has no external plasma source (the neutron star), and it includes inverse Compton scattering and pair production via interactions with an external radiation field in addition to curvature emission.

Oscillating gap model
We study the dynamics of a local, 1D gap using fully general relativistic PIC simulations. We assume that radiation emanating from the inner regions of the accretion flow provides the dominant source of opacity for pair production and inverse Compton scattering. The intensity of this radiation is given as an input for the computations of the gap dynamics. Gamma rays generated via inverse Compton scattering of the ambient radiation by the pairs accelerated in the gap are treated as a neutral species in our PIC scheme (in addition to electrons and positrons). Pair production occurs through the interaction of the gamma rays thereby produced with the external photons emitted by the accretion flow. The various radiation processes are computed using a novel Monte Carlo method developed for this purpose (see Appendix B for details). In addition, we include curvature losses in the equations of motions and provide a rough estimate for the luminosity of curvature emission, as explained below. However, for numerical reasons we do not add the curvature photons to the pull of gamma rays, and therefore do not account for gamma ray production by the curvature photons. As shown below, for cases of interest curvature emission is only important during the initial discharge of the gap.
We implicitly assume that the gap constitutes a small disturbance in the magnetosphere, in the sense that its activity does not significantly affect the global structure, and in particular the magnetic field geometry and the angular velocity of magnetic surfaces Ω. The coupling between the gap and the global magnetosphere enters through the global electric current flowing in the magnetosphere, treated as a free input parameter, as in the steady-state models Levinson & Segev 2017). For simplicity we adopt a split monopole geometry, defined by A ϕ = C(1 − cos θ). For this choice F rϕ = 0, and from the ideal MHD condition we obtain F rt = −ΩF rϕ = 0 for the radial electric field outside the gap, in the ideal MHD sections of the magnetosphere. This effectively ignores any MHD waves that might be generated by the gap cycle and propagate throughout the force-free magnetosphere. The gap extends along a poloidal magnetic surface, characterized by an inclination angle θ. Inside the gap F rt 0 by virtue of the pair creating oscillations. This is the only wave field that appears in the dynamical equations derived below, hence the gap activity is restricted to longitudinal plasma oscillations in this model. Despite this restriction, this model captures the main features of plasma production and consequent emission. We consider it as a preliminary stage in our quest for the development of a 2D code that computes the global structure and dynamics of the magnetosphere.

Background geometry and choice of coordinates
The background spacetime is described by the Kerr metric, here given in Boyer-Lindquist coordinates with the notation ds 2 = −α 2 dt 2 + g ϕϕ (dϕ − ωdt) 2 + g rr dr 2 + g θθ dθ 2 , where (2) with ∆ = r 2 + a 2 − 2r g r, Σ = r 2 + a 2 cos 2 θ, A = (r 2 + a 2 ) 2 − a 2 ∆ sin 2 θ, and r g = GM/c 2 = 1.5×10 14 M 9 cm denotes the gravitational radius, and M = 10 9 M the BH mass. The parameter a = J/M represents the specific angular momentum. The determinant of the matrix g µν is given by √ −g = Σ sin θ. The angular velocity of the black hole is defined as ω H = ω(r = r H ) =ã/2r H , whereã = a/r g denotes the dimensionless spin parameter, and r H = r g + r 2 g − a 2 is the radius of the horizon. Hereafter, all lengths are measured in units of r g and time in units of t g = r g /c, so we set c = r g = 1, unless explicitly stated otherwise.
To avoid the singularity on the horizon, we find it convenient to transform to the tortoise coordinate ξ, defined by dξ = dr/∆ (in full units dξ = r 2 g dr/∆). It is related to r through with r ± = 1 ± √ 1 −ã 2 . We note that ξ → −∞ as r → r H = r + , and ξ → 0 as r → ∞. We prefer to use this version of the tortoise coordinate for the following reasons: (i) it pushes the horizon to −∞, thereby allowing us to conveniently choose the inner gap boundary as close to the horizon as needed; (ii) when using quantities measured in the frame of a zero angular momentum observer (ZAMO), the form of the equations in these coordinates is very similar to that in flat spacetime. This renders interpretation of the results more intuitive; and, most importantly, (iii) it does not mix the radial coordinate with the time and azimuthal coordinates, as it does in the case of Kerr-Schild coordinates, for example. This greatly simplifies the Monte Carlo computations of Compton scattering and pair production. When using Kerr-Schild coordinates additional transformations are needed in every computation step, because the emissivity and absorption coefficient are naturally defined in the ZAMO frame. This is an unnecessary complication that is avoided by our choice of the tortoise coordinates.
For a grid cell δr on a magnetic surface having an inclination angle θ, a volume element can be defined as Here, δV ξ defines the volume element with respect to the coordinate ξ, and is finite on the horizon. The average number density (of either pairs or gamma rays) within a grid cell, as measured by a distant observer, is given by where δN denotes the occupation number of the grid cell. The quantity n ξ = ∆n is finite on the horizon, and will be used hereafter to describe the density of the different plasma constituents.

Electrodynamics
As explained above, the only wave field in our model is the radial component of the electric field, F rt (see Appendix A for details). Measured in the ZAMO frame it reads E r = √ AF rt /Σ. Its dynamics is governed by Equation (A.6), here expressed in terms of the electric flux per steradian, where J 0 = Σ j r 0 is the global magnetospheric current, defined explicitly in Equation (A.5). This current is conserved along magnetic surfaces, and serves as an input parameter to the dynamic gap model. Gauss's law further yields where is the GJ density, given explicitly in Equation (A.3), B H is the strength of the magnetic field on the horizon, and A H ≡ A(r H ). We note that charge conservation readily implies that Eq. (7) is conserved in time. That is, if the initial state satisfies this equation, and the system then evolves in time according to Eq. (6), it is guaranteed that Eq. (7) will be satisfied at any given time. Thus, this equation is only used once at the beginning of each run to determine the initial state of the electric field, E r (r, t = 0), for a given choice of initial conditions.

Particle motion
The plasma in the gap consists of electrons and positrons. Let u µ i denote the four-velocity of the ith particle and q i its electric charge, where q i = −e for electrons and q i = +e for positrons. The equation of motion for the ith particle can then be expressed as subject to the normalization u µ i u i µ = −1, where s µ i is a source term associated with curvature losses, satisfying u iµ s µ i = 0 , and dτ i is the corresponding proper time interval. It is related to the time measured by a distant observer through dt = u t i dτ i . For the lowered index components, Poloidal motion is restricted to the radial direction in our 1D model, thus u θ i = 0. Furthermore, since ∂ ϕ is a Killing vector it readily follows that Γ α ϕ β u α i u β i = 0, and since F ϕt = F ϕr = 0 we obtain F ϕα u α i = F ϕθ u θ i = 0. As the curvature drag term, s iϕ , is proportional to u ϕ , it is evident from Equation (10) that particles tend to a state of zero angular momentum, hence we set u ϕ = s iϕ = 0, for which we obtain s ϕ i = ωs t i and s it = −α 2 s t i . We find it convenient to use the four-velocity components measured by a ZAMO, here denoted by u i = √ g rr u r i , the particle Lorentz factor γ i = αu t i , and the three-velocity v i = u i /γ i . The normalization condition yields γ 2 i = 1 + u 2 i . Using Equations (9) and (10) and the relations 2u i du i /dτ i = u ir du r i /dτ i + u r i du ri /dτ i , dτ i = dt/u t i , and u r i s ir + u t i s it = 0, we arrive at where −s t /u t = P cur (γ i ) is the curvature power emitted by the particle, as measured by a ZAMO (Equation (21) below), and The particle trajectory is computed using In terms of v i the electric current density in Equation (6) can be expressed as Likewise, the charge density in Equation (7) satisfies

Radiation
Gamma rays are treated as a neutral species in our scheme. Let p µ k denote the momentum of the kth photon, as measured by a ZAMO. It is related to the coordinate momentum p µ k though the energies of interest the gamma rays are highly beamed along the direction of motion of the emitting particles by virtue of momentum conservation, thus to a very good approximation we can setp θ k =p ϕ k = 0. The normalization condition, p kµ p µ k = 0, gives p t k = |p r k |, as required by the fact that photons propagate at the speed of light in the ZAMO frame. The null geodesic equations reduce to dp r and the photon trajectory equation to We suppose that the gap is exposed to the emission of soft photons by the accretion flow from a putative source of size R s = R s r g and luminosity L s = l s L Edd . For simplicity, we assume that the intensity of the seed radiation inside the simulation box is isotropic, constant in time, and homogeneous, with a power law spectrum Article number, page 3 of 14 A&A proofs: manuscript no. gap_pic_pap-AA with p > 1, where s is the photon energy in m e c 2 units, as measured by a ZAMO. The assumption that I s is isotropic is reasonable, except perhaps very near the horizon, since the size R s of the radiation source is typically larger than the gap dimensions. The total number density of photons is given approximately by n s = 4π (I s /hc s )d s 4πI 0 /hc, and can be used to define the fiducial optical depth ForR s = 10 we have l s 10 −2 min τ 0 . The interaction of pairs and gamma rays with the external radiation field is computed using a Monte Carlo approach. The details are given in Appendix B. In short, for every particle we compute the optical depth traversed by the particle between two consecutive time steps, t and t + δt, along its trajectory r(t), using the full Klein-Nishina opacity κ c , given in Equation (B.1): We then randomly draw a probability for a scattering event: 0 ≤ p sc ≤ 1. A scattering event occurs if p sc < 1 − e −δτ sc . The process is repeated at every time step for all particles in the simulation box. If a scattering event occurs, we transform to the rest frame of the particle and draw the energy and direction of the scattered gamma ray, as explained in Appendix B, and then transform back to the ZAMO frame and compute the new energy of the particle. The newly created gamma ray is added to the pull. Pair production via the interaction of gamma rays with the external radiation field is computed in a similar manner. We use Equation (20) with κ c replaced by the pair production opacity κ pp (Eq. (B.18)) to calculate the optical depth δτ ± traversed by a gamma ray between two consecutive time steps. Pair production occurs if p ± < 1 − e −δτ ± , upon drawing a probability 0 ≤ p ± ≤ 1. The process is repeated at every time step for all gamma rays in the simulation box. The pull of gamma rays and pairs is updated correspondingly. To simplify the analysis we suppose that in each pair production event the newly created electron and positron have identical energies. This is a good approximation near the threshold, where the cross section is at maximum (e.g., Blandford & Levinson 1995).
The curvature power emitted by an electron (positron) of energy γ in the ZAMO frame is given by where R c is the curvature radius of the particle trajectory. For the computations presented below we adopt R c = r g . The characteristic energy of a curvature photon (measured in m e c 2 units) is where λ c = /m e c is the Compton wavelength of the electron. The maximum Lorentz factor of the pairs is limited by back reaction. From Equation (11) we find where E 0 = ηB H = η10 3 B 3 G is the maximum strength of the electric field in the gap (see middle left panel in Fig. 1). This maximum value is delineated by the vertical dashed line in Fig.  2. As will be shown, in the initial discharge pairs indeed accelerate to this Lorentz factor; however, after the relaxation of the system the Lorentz factor is essentially limited by the oscillations of the electric field rather than curvature losses, provided τ 0 > 1. In a transparent gap, τ 0 < 1, we expect domination of curvature radiation during the entire evolution of the system. The number of curvature photons emitted by an electron (positron) of Lorentz factor γ over a dynamical time t g = r g /c is It varies between 10 7 during the initial discharge (where γ γ max ) to about 10 5 during the relaxed state. Thus, without proper resampling it is practically infeasible to track these photons in our PIC scheme. We note that once γ drops below 10 9 the characteristic energy of curvature photons becomes c < 10 3 , and their contribution to pair creation and to the gamma-ray power emitted from the gap can be neglected. However, curvature emission dominates the released power (and likely pair creation) during the initial discharge. In order to account for this power, yet avoiding tremendous numerical complications, the net curvature luminosity measured by a distant observer is computed, at any given time step, by summing up the contributions of all particles in the simulation box, accounting properly for redshift effects with α i ≡ α(r i ) denoting the lapse function of the ith particle, currently located at radius r i . This prescription lacks proper treatment of time travel effects, which are particularly important near the horizon. It merely provides a rough estimate of the contribution of curvature emission to the emitted gamma-ray luminosity.

Input parameters
The 1D gap model is characterized by the following input parameters: the black hole mass (M = 10 9 M 9 M ), the BH spin parameter (ã), the angular velocity (Ω) of the magnetic surface along which the gap lies, the strength of the magnetic field on the horizon (B H = 10 3 B 3 G), the inclination angle of magnetic surface (θ), the minimum energy and slope of the target spectrum ( s,min and p) (Eq. 18), the global electric current (J 0 ) defined in Equation (A.5), and the fiducial optical depth (τ 0 ) defined in Equation (19). The results presented below were computed using the following canonical choice of parameters (M 9 = 1, B 3 = 2π, θ = 30 • ,ã = 0.9, s,min = 10 −8 , p = 2), which correspond to supermassive BHs accreting in the RIAF regime. For this choice of parameters we explore how the behavior of the solutions depends on J 0 and τ 0 .

Numerical methods
The overall structure of the code is based on the 1D version of the ZELTRON code (Cerutti et al. 2013), a highly parallelized (special) relativistic electromagnetic PIC code, in addition to which GR corrections and the Monte Carlo scheme introduced in the previous section were implemented for this study. The code solves the equations of motions and Maxwell's equation using an explicit second-order finite-difference scheme.
The electric field is evolved in time using Eq. (6), such that where δt = t n+1 − t n = t n+1/2 − t n−1/2 is the time step. The electric field is defined at full time steps n and n + 1, while the current is defined at half time steps n + 1/2. This offset in time ensures a second order accuracy of the scheme. Both E r and j r are defined at the nodes of the grid given by the index iξ. Currents from individual charged particles are deposited on the grid using Eq. (14) where the three-velocities are known at the half time step n+1/2. In the 1D limit, the equations of motions are also straightforward to integrate even though more steps are needed than for the field. Starting with the equation of motion of the photons (Eq. 16), a time-centered scheme gives Assuming thatp r,n k = (p r,n+1/2 k +p r,n−1/2 k )/2, we obtaiñ The photon position at time t n+1 is then given by For the charged particles, we have (see Eq. 11) We break up this equation into three components corresponding to each term on the right-hand side The next step is to estimate the particle Lorentz factor and threevelocity at time step t n . To do this, we first compute u n+1/2 i,L using Eq. (32) where E n r is already known at the nodes of the grid and linearly interpolated at the particle position. Using the same trick as for the photons, i.e., u n i,L = (u n+1/2 i,L The last step is to inject these estimated values into Eq. (34), such that The charged particle positions are updated in the same way as for photons,

Numerical setup
The spatial grid is fixed in time and uniform in ξ, ranging from ξ min = −3 (r min ≈ 1.5) to ξ max = −0.3 (r max ≈ 4.3). Therefore, the grid is highly non-uniform in radius, refined near the BH horizon and sparse in the outer regions of the box. The time step is set by the Courant-Friedrichs-Lewy condition defined at the inner boundary, i.e., δt ≤ δt CFL ≡ r 2 min +ã 2 − 2r min δξ.
The simulation box is initially filled with a monoenergetic beam of gamma-ray photons uniformly distributed along the ξ-grid for reasons explained below, but it is empty of charged particles. The initial (vacuum) electric field is obtained by numerically integrating Eq. (7) with j t = 0. The solution is shown in the middle-left panel in figure 1. Gauss's law is integrated only once at the beginning of the simulation. We have checked that it is well satisfied throughout the simulation by virtue of charge conservation. The choice of boundary conditions for the fields are trivial in this problem because Ampère's law given in Eq. (6) does not involve any spatial derivative in 1D. The electric field is free to evolve throughout the box. The particles, whether charged or neutral, are simply deleted from the memory as soon as they cross the boundaries to mimic an open boundary on both sides. Therefore, we assume that no plasma injection from outside is permitted.
Spatial and temporal resolutions are harder quantities to define in this setup because they essentially depend on the energy and the density of particles, which are the unknowns that we are trying to measure here. Thus, resolution and numerical convergence can be checked a posteriori only. More specifically, it is crucial to resolve the collisionless plasma skin depth, l p , otherwise the plasma will heat up artificially until it is resolved by the grid. It is instructive to give an estimate of the skin depth in terms of the model parameters, the pair multiplicity, κ = n ± /n GJ , where n ± denotes the pair density, and the mean Lorentz factor of pairs γ . Recalling that the plasma frequency is ω p = 4πe 2 n ± /m e γ and en GJ = ΩB/2πc, and adopting Ω = ω H /2, we obtain and r g l p 10 7 κM 9 B 3 2 γ .
The number of cells needed to resolve the skin depth in a simulation box of size h is approximately h/l p = (h/r g )(r g /l p ). In the cases studied below we find κ/ γ in the range 10 −9 − 10 −7 for a range of opacities τ 0 = 1 to τ 0 = 10. For the size of our simulation box, h/r g ≈ 3, at least 10 4 cells are needed to resolve the skin depth for τ 0 = 10. Simulations with τ 0 < 10 showed good convergence for a total N ξ = 16384 cells. For τ 0 = 10, we had to run with up to N ξ = 65536 cells to see convergence. For the results described below, the initial beam of gamma rays is modelled with five particles per cell. At the end of the simulation, the average number of particles per cell varies from 2 (τ 0 = 1) to 10 (τ 0 = 10) for the pairs and from 2 to 50 for the photons. We have also checked the good convergence of our results with respect to the initial number of injected gamma-ray photons per cell.

Overall evolution
We ran simulations for a grid of models characterized by different values of the parameters τ 0 and J 0 . Quite generally, we find that the evolution of the system depends on the value of τ 0 , but is practically independent of J 0 (Fig. 3). The prime role of J 0 is to fix the time average value of the oscillating current. In all cases explored we find an initial discharge of the gap that produces a gamma-ray flare, dominated by curvature losses, with a duration of about ∆t ∼ r g /c and a luminosity that approaches the maximum allowed power, L γ ∼ χL BZ , where χ < ∼ 1 depends on the initial gap width (see Eq. (42) and text below), followed by rapid, small amplitude oscillations that last for the entire simulation time, during which the pair plasma is continuously replenished through self-sustained pair creation bursts. The longterm behavior of the system (after the decay of the prompt spark) is essentially independent of the initial condition, provided that the initial pair density is well below the GJ density inside the gap. This behavior is reminiscent of the longitudinal oscillations found in the semi-analytic, two-beam model of Levinson et al. (2005).
A typical example is shown in Figs. 1 and 2. Figure 1 exhibits four snapshots of the pair and photon densities (upper panels), electric flux (middle panels), and radial electric current (lower panels), and Fig. 2 the corresponding spectral energy distributions of pairs and photons. The initial pair density was taken to be zero in this example. To ignite the discharge process a monoenergetic gamma-ray distribution (left panel of Fig. 2) was injected at time t = 0 inside the simulation box, with a uniform distribution in ξ space. The initial energy of injected photons was chosen to optimize pair creation in order to speed up the computations. We also performed several runs with the same setup but different initial conditions, e.g., sub-GJ pair density and no photons, and verified that the overall evolution of the system is independent of the initial condition. In all cases, we find that after several crossing times the system approaches a state of quasisteady oscillations with essentially the same properties (as in the rightmost panels in figures 1 and 2).
As seen in figure 1, after about one crossing time from the start of the simulation the pair density exceeds the GJ value, giving rise to a nearly complete screening of the electric field. The sporadic pair creation leads to rapid spatial and temporal oscillations of the electric field and radial current inside the simulations box. The pair density and the amplitude of the electric field oscillations ultimately approach their saturation levels (after several crossing times), at which time pair creation inside the simulation box balances pair escape through its boundaries. At this stage the pair and photon spectra are quasi-stationary. Figure 3 depicts the dependence of the pair multiplicity (upper panels), average Lorentz factor (middle panels) and average gamma-ray energy (lower panels) on the fiducial opacity τ 0 (left) and global current J 0 (right). As can be seen, the multiplicity increases roughly linearly with τ 0 , while the average pair and photon energies decrease with increasing τ 0 . On the other hand, those quantities are essentially independent of J 0 . The overall gamma-ray luminosity emitted from the gap following the prompt phase depends weakly on τ 0 (figure 4). For the BZ power used to normalize the luminosities exhibited in figures 4 and 5 we adopt This is close to values obtained from solutions to the Grad-Shafranov equation (Nathanail & Contopoulos 2014;Mahlmann et al. 2018), but may somewhat overestimate the values expected in realistic jets.

Flaring states
During the initial discharge (at t < r g /c), when the gap electric field is still sufficiently intense, the newly created pairs quickly accelerate to the terminal Lorentz factor γ max (Eq. (23)), at which time the energy gain is balanced by curvature losses (indicated by the vertical dashed line in figure 2). We note that at these energies Compton scattering is in the deep Klein-Nishina regime and is highly suppressed. As time passes the average pair energy declines, ultimately approaching a final value. For the case shown in figure 2 it is about 5 × 10 7 . At this state curvature emission is completely negligible (see Fig. 5). As is evident from figure 2, during the initial discharge most particles quickly accelerate to the terminal Lorentz factor, γ i = γ max , where P cur (γ i ) eE. The peak luminosity occurs roughly when the pair multiplicity approaches unity, while the electric field is not yet significantly screened out. Thus, the curvature luminosity can be approximated as L cur P cur n GJ 4πr 3 H /3 eEn GJ 4πr 3 H /3. From Equation (7) we estimate that E en GJ r H before screening ensues. With en GJ ΩB H /2π and Ω = ω H /2 =ã/4r H , we obtain where L BZ is the corresponding Blandford-Znajeck power given explicitly in Eq. (41). Figure 5 exhibits the light curve computed from the PIC simulations. As seen, the gamma-ray luminosity indeed approaches the BZ power during the initial spark, and then decays to the terminal value as the gap electric field is screened out. It can be readily shown that if the initial gap width is small, h r H , then L cur is reduced by a factor χ (h/r H ) 2 . The above considerations suggest that strong rapid flares should be produced every time a magnetospheric gap is restored, for example by an accretion episode. The flaring episode is then followed by quiescent emission with a luminosity L γ ∼ 10 −5 L BZ .

Applications to M87
M87 exhibits TeV emission with a luminosity of L T eV ∼ 10 40 erg s −1 in the quiescent state. Several strong flares with durations of ∆t t g have been recorded in the past decade Albert et al. 2008;Acciari et al. 2009;Abramowski et al. 2012). Various estimates of the average jet power in M87 (see, e.g., Bicknell & Begelman 1996;Owen et al. 2000;Stawarz et al. 2006;Bromberg & Levinson 2009) yield a range of a few times 10 43 to a few times 10 44 erg s −1 . Assuming that the jet is powered by the BZ mechanism implies L T eV /L BZ ∼ 10 −4 . The SED exhibits a peak in the sub-mm band at ν peak 10 12 Hz, corresponding to max 10 −8 in our parametrization, with a bolometric luminosity of L b < ∼ 10 41 erg s −1 . This sets an upper limit on the luminosity of the putative RIAF emission, L s = η d L b , as the source of the observed SED is yet unresolved. For a BH mass of M = 6 × 10 9 M this implies l s < 10 −7 and τ 0 = 10 3 η d (R s /10) −2 < 10 3 from Equation (19). If the observed TeV emission originates from a spark gap, then the pair production opacity must be sufficiently low to allow TeV photons to escape the system. For the target photon spectrum invoked in Equation (18) with p = 2, min = 10 −8 , the pair production optical depth is given approximately by τ γγ 0.1τ 0Rs ( max γ ) 2 = 10 −13 η d (R s /10) −1 2 γ . Thus, , and radial electric current, Σ j r , normalized by the global magnetospheric current J 0 . The leftmost panels delineate the initial state, at t = 0. The rightmost panels show the relaxed state, following the prompt discharge. We note the scale change on the vertical axis in the middle panels. ct/r g =5.07 10 5 10 6 10 7 10 8 10 9 10 10 10 11 u 10 3 10 4 10 5 10 6 10 7 ct/r g =20.06 it is transparent at energies below about 3 TeV. However, the observed spectrum appears to extend up to ∼ 10 TeV, implying η d < ∼ 10 −1 and τ 0 < 10 2 . Speculating that the observed TeV emission is produced in a spark gap, we note that the gamma-ray power released by the gap, as predicted by our model, L γ ∼ 10 −5 L BZ , is somewhat lower, but still consistent with the observed emission in the quiescent state given the various model uncertainties. Figure  2 confirms that for τ 0 = 10 the spectrum extends up to about 10 TeV. The strong flares recorded in the past two decades might be caused by some episodes during which the gap is suddenly restored. We note that a small percent of the full gap width is sufficient to account for the strongest flares observed by a sudden discharge. Alternative models for the M87 flares include mis-  aligned mini jets (Giannios et al. 2010), and star-jet interactions (Barkov et al. 2010, see also Aharonian et al. 2017).

Conclusions
We explored the dynamics of pair discharges in a starved magnetosphere of a Kerr BH, using 1D PIC simulations. Our analysis takes into account inverse Compton scattering and pair creation via the interaction of pairs and gamma rays with an ambient radiation field, assumed to be emitted by the putative accretion flow, as well as curvature losses. In the computations presented here the intensity of the external radiation field is taken to be a power law with a slope of −2 and a minimum energy of 10 −8 m e c 2 (corresponding to a minimum frequency of about 10 12 Hz). We find, quite generally, that the initial gap electric field is screened out by a prompt discharge of duration ∼ r g /c that produces a strong flare of very high-energy curvature photons. This episode is followed by a state of self-sustained, rapid plasma oscillations that lasts for the entire simulation time, during which the pair and gamma-ray spectra are quasi-stationary, with a gamma-ray luminosity that constitutes a fraction of about 10 −5 of the BZ power depending weakly on input parameters. As pointed out in §5, this value is in rough agreement with observations of M87. It is worth noting that the ratio L γ /L BZ obtained during the quiescent state in our model is typically smaller than the values obtained in steady-state models with extremely low accretion rates. For example, Hirotani et al. (2017) obtained values in the range 10 −4 − 10 −1 for accretion rates in the range from may be needed. On the other hand, for the same accretion rates the BZ power can be larger by up to an order of magnitude if the magnetic field approaches the saturation level predicted by MAD models, as seems to be suggested by observations of M87. This will give rise to correspondingly higher gamma-ray luminosities.
Our choice of parameters in the present analysis was motivated by the applications to M87 and conceivably other AGNs. However, the gap emission may also be relevant to galactic BH transients Levinson & Segev 2017), which span a different regime in parameter space. If similar L γ /L BZ ratios are produced by Galactic BHs, then a 10 solar mass BH accreting at a rate of 10 −3 Eddington, can conceivably produce a TeV luminosity of about 10 31 erg s −1 that can be detected by the Cherenkov Telescope Array out to a distance of about 1 kpc. However, it is plausible that in these sources the emission will be dominated by curvature losses due to the smaller curvature radius, and it remains to be seen how the shape of the emitted spectrum scales with source parameters. We leave this problem to a future work.
Ultimately, global PIC simulations are required to compute the full structure of the magnetosphere and its response to pair discharges in starved region. Moreover, our 1D model is restricted to longitudinal plasma oscillations, while in reality transverse modes might be excited by the gap activity, which can only be accounted for in 2D simulations. Nonetheless, our model captures the basic features of plasma production and gap emission.

Appendix A: Derivation of the gap electrodynamic equations
We restrict our analysis to axisymmetric systems. We suppose that outside the gap the flow is stationary, and that inside the gap particles can only move along magnetic surfaces; i.e., the gap oscillations are longitudinal (electrostatic). In the force-free section outside the gap the angular velocity of field lines, Ω, is fixed. If the gap forms a small perturbation in the magnetosphere, then variation in Ω across the gap can be ignored. It is then appropriate to define the electric field in the corotating frame as F µt = F µt + ΩF µϕ . Outside the gap, in the ideal MHD section, F µt = 0. Inside the gap, Gauss's law, ∂ µ ( √ −gF tµ ) = √ −g j t , reduces to (see, e.g., Levinson & Segev 2017) ∆ sin θ α 2 F rt ,r where the GJ density is given by and the last equality applies to the split monopole geometry invoked in our model, A ϕ = C(1 − cos θ). We note that √ AB r sin θ = F θϕ = C sin θ, which implies that √ AB r = √ A H B H = C, here B H = B r (r = r H ) and likewise A. Explicitly: Since the gap dynamics is restricted to longitudinal oscillations we have F θt = 0 in Eq. (A.1). Next, the radial component of and for our split monopole geometry F rϕ = 0. Outside the gap the flow is stationary (∂ t = 0), and charge conservation, ∂ µ ( √ −g j µ ) = 0, yields ∂ r (Σ j r ) = 0 in the force-free limit, so that the radial electric current is conserved along magnetic surfaces: Σ j r = J 0 = const. From Equation (A.4) we obtain outside the gap. This conserved current must be flowing through the gap. Thus, the induction equation, Eq. (A.4), reduces to Frame dragging couples the toroidal and poloidal components of the wave field: From the latter relations we find |E ϕ | ∼ (λ/r)|E r | |E ϕ |, where λ r is the characteristic wavelength of the oscillations, and likewise for B ϕ , thus those fields can be neglected.

Drawing the scattered photon energy
Once a scattering event occurs, the energy of the scattered photon is drawn upon transforming to the rest frame of the scatterer. The probability density that an electron (positron) will produce a gamma ray of rest frame energy γ = γ(1 − βµ γ ) γ in a single scattering can be expressed as with −β holding for s > min /γ(1 − β), and βµ max = min {β, ( max /γ s ) − 1}, with β holding for s < max /γ(1 + β), and where cos Θ = cos θ s cos θ γ −sin θ s sin θ γ cos(ϕ s −ϕ γ ) is the angle between incident and scattered photons in the particle's rest frame, and dσ( s , γ )/dΩ is obtained by substituting cos Θ = 1 + 1/ s − 1/ γ in Equation (B.6).