Issue 
A&A
Volume 616, August 2018



Article Number  A184  
Number of page(s)  12  
Section  Astrophysical processes  
DOI  https://doi.org/10.1051/00046361/201832915  
Published online  11 September 2018 
Particleincell simulations of pair discharges in a starved magnetosphere of a Kerr black hole
^{1}
Univ. Grenoble Alpes, CNRS, IPAG, 38000 Grenoble, France
email: benoit.cerutti@univgrenoblealpes.fr
^{2}
Raymond and Beverly Sackler School of Physics & Astronomy, Tel Aviv University, Tel Aviv 69978, Israel
Received:
27
February
2018
Accepted:
21
May
2018
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 particleincell code that implements Monte Carlo methods to compute gammaray 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 lowamplitude, rapid plasma oscillations that produce selfsustained pair cascades, with quasistationary pair and gammaray spectra, and with a pair multiplicity that increases in proportion to the pair production opacity. The gammaray 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 ∼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.
Key words: black hole physics / acceleration of particles / radiation mechanisms: nonthermal / methods: numerical / galaxies: individual: M87 / gamma rays: galaxies
© ESO 2018
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
1. 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 selfsustained 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; Hirotani & Pu 2016; Hirotani et al. 2016; 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; Hirotani et al. 2016; 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; Hirotani et al. 2016) 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; Hirotani & Pu 2016; Hirotani et al. 2016; Lin et al. 2017; Katsoulakos & Rieger 2018) that the gap activity may be imprinted in the highenergy 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; Hirotani & Pu 2016).
In an attempt to compute the emission properties of an active gap, a fully general relativistic (GR) steady gap model has recently been developed (Hirotani & Pu 2016; Hirotani et al. 2016; 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), steadystate 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 particleincell (PIC) code, developed particularly for this purpose, that implements Monte Carlo methods to compute gammaray emission and pair production through the interaction of pairs and gamma rays with an external radiation field. The details are described in Sects. 2 and 3 below. We find (Sect. 4) that after a prompt discharge phase that screens out the gap and produces a strong gammaray flare, the system relaxes to a state of selfsustained, 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 quasistationary during this state. We find that the saturated pair multiplicity increases roughly linearly with the opacity contributed by the ambient radiation field, while the average pair and photon energies decrease as the opacity increases. The gammaray 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 Sect. 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.
2. 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 steadystate models (Hirotani & Pu 2016; 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 magnetohydrodynamic (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 forcefree 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.
2.1. Background geometry and choice of coordinates
The background spacetime is described by the Kerr metric, here given in Boyer−Lindquist coordinates with the notation (1) 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 . The angular velocity of the BH is defined as ω_{H} = ω(r = r_{H}) = ã/2r_{H}, where ã = a/r_{g} denotes the dimensionless spin parameter, and 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 ). It is related to r through (3) with . 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 (4) 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 (5) 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.
2.2. Basic equations
2.2.1. 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 . Its dynamics is governed by Eq. (A.6), here expressed in terms of the electric flux per steradian, , as (6) where is the global magnetospheric current, defined explicitly in Eq. (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 (7) where (8) is the GJ density, given explicitly in Eq. (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.
2.2.2. Particle motion
The plasma in the gap consists of electrons and positrons. Let denote the fourvelocity 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(9) subject to the normalization , where is a source term associated with curvature losses, satisfying , and dτ_{i} is the corresponding proper time interval. It is related to the time measured by a distant observer through . For the lowered index components, , we likewise have (10)
Poloidal motion is restricted to the radial direction in our 1D model, thus . Furthermore, since ∂_{φ} is a Killing vector it readily follows that , and since F_{φt} = F_{φr} = 0 we obtain . As the curvature drag term, s_{iφ}, is proportional to u_{φ}, it is evident from Eq. (10) that particles tend to a state of zero angular momentum, hence we set u_{φ} = s_{iφ} = 0, for which we obtain and . We find it convenient to use the fourvelocity components measured by a ZAMO, here denoted by , the particle Lorentz factor , and the threevelocity v_{i} = u_{i}/γ_{i}. The normalization condition yields . Using Eqs. (9) and (10) and the relations , , and , we arrive at (11) where −s^{t}/u^{t} = P_{cur}(γ_{i}) is the curvature power emitted by the particle, as measured by a ZAMO (Eq. (21) below), and (12) The particle trajectory is computed using (13) In terms of v_{i} the electric current density in Eq. (6) can be expressed as (14) Likewise, the charge density in Eq. (7) satisfies (15)
2.2.3. Radiation
Gamma rays are treated as a neutral species in our scheme. Let denote the momentum of the kth photon, as measured by a ZAMO. It is related to the coordinate momentum though , , , and . At 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 set . The normalization condition, , gives , as required by the fact that photons propagate at the speed of light in the ZAMO frame. The null geodesic equations reduce to (16) and the photon trajectory equation to (17)
We suppose that the gap is exposed to the emission of soft photons by the accretion flow from a putative source of size 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 (18) 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 (19) For 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 Eq. (B.1): (20) 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 Eq. (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 (21) 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 (22) 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 Eq. (11) we find (23) 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 (24) 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 gammaray 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 (25) 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 gammaray luminosity.
Fig. 1. Snapshots from a typical simulation of spark gap dynamics. Shown are the pair and photon densities (upper panels), electric flux (middle panels), 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. 
Fig. 2. Snapshots of electron (blue line), positron (red line), and photon (yellow line) spectra computed during the run presented in Fig. 1. 
2.3. Input parameters
The 1D gap model is characterized by the following input parameters: the BH 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 Eq. (A.5), and the fiducial optical depth (τ_{0}) defined in Eq. (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}.
3. Implementation of the 1D GRPIC code
3.1. 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 secondorder finitedifference scheme.
The electric field is evolved in time using Eq. (6), such that (26) 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 threevelocities 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 timecentered scheme gives (27) Assuming that , we obtain (28) The photon position at time t^{n+1} is then given by (29) For the charged particles, we have (see Eq. (11)) (30) We break up this equation into three components corresponding to each term on the righthand side (31) (32) (33) Assuming that and adding all three equations together yields (34) The next step is to estimate the particle Lorentz factor and threevelocity at time step t^{n}. To do this, we first compute using Eq. (32) where 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., , we compute (35) The last step is to inject these estimated values into Eq. (34), such that (36) The charged particle positions are updated in the same way as for photons, (37)
3.2. 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 nonuniform 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., (38) The simulation box is initially filled with a monoenergetic beam of gammaray 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 middleleft panel in Fig. 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 and en_{GJ} = ΩB/2πc, and adopting Ω = ω_{H}/2, we obtain (39) and (40) 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 gammaray photons per cell.
4. Results
4.1. 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 gammaray 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 selfsustained 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 semianalytic, twobeam 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 gammaray 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., subGJ 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 Figs. 1 and 2).
As seen in Fig. 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 quasistationary.
Figure 3 depicts the dependence of the pair multiplicity (upper panels), average Lorentz factor (middle panels) and average gammaray 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 gammaray luminosity emitted from the gap following the prompt phase depends weakly on τ_{0} (Fig. 4). For the BZ power used to normalize the luminosities exhibited in Figs. 4 and 5 we adopt (41) 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.
Fig. 3. Dependence of the pair multiplicity ⟨κ⟩, pair energy ⟨γ⟩, and gammaray energy ⟨ϵ⟩ averaged over the simulation domain as functions of the input parameters τ_{0} (right) and J_{0} (left). 
Fig. 4. Dependence of quiescent IC, curvature, and kinetic luminosities leaving the outer boundary of the box (r_{out} ≈ 4) on the fiducial opacity τ_{0}. 
Fig. 5. Gammaray light curve produced by the gap discharge. The red line corresponds to IC emission, the blue line to curvature emission, and the black dashed line to the sum of both components. 
4.2. 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 Fig. 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 Fig. 2 it is about 5 × 10^{7}. At this state curvature emission is completely negligible (see Fig. 5). As is evident from Fig. 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 . From Eq. (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 (42) 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 gammaray 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}.
5. Applications to M87
M87 exhibits TeV emission with a luminosity of L_{TeV} ∼ 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 (Aharonian et al. 2006; 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_{TeV}/L_{BZ} ∼ 10^{−4}. The SED exhibits a peak in the submm 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 from Eq. (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 Eq. (18) with p = 2, ϵ_{min} = 10^{−8}, the pair production optical depth is given approximately by . Thus, 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 gammaray 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 misaligned mini jets (Giannios et al. 2010), and star−jet interactions (Barkov et al. 2010, see also Aharonian et al. 2017).
6. 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 highenergy curvature photons. This episode is followed by a state of selfsustained, rapid plasma oscillations that lasts for the entire simulation time, during which the pair and gammaray spectra are quasistationary, with a gammaray luminosity that constitutes a fraction of about 10^{−5} of the BZ power depending weakly on input parameters. As pointed out in Sect. 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 steadystate 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 10^{−4} to 5×10^{−6} Eddington, assuming an equipartition magnetic field in the disk. Thus, if the gaps are unsteady, as our simulations seem to indicate, reevaluation of observational predictions 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 gammaray 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 (Lin et al. 2017; 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.
Acknowledgments
We thank Maxim Barkov and Alexander Philippov for their comments on the manuscript. AL acknowledges the kind hospitality of IPAG, where the essential part of this work was done, and the support from the visiting professor program of the Université Grenoble Alpes. BC acknowledges support from CNES and Labex OSUG@2020 (ANR10 LABX56). This work was granted access to the HPC resources of TGCC under the allocation t2016047669 made by GENCI.
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 forcefree 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 . Outside the gap, in the ideal MHD section, . Inside the gap, Gauss’s law, , reduces to (see, e.g., Levinson & Segev 2017) (A.1) where the GJ density is given by (A.2) and the last equality applies to the split monopole geometry invoked in our model, A_{φ} = C(1 − cosθ). We note that , which implies that , here B_{H} = B_{r}(r = r_{H}) and likewise A. Explicitly: (A.3) Since the gap dynamics is restricted to longitudinal oscillations we have in Eq. (A.1). Next, the radial component of Ampère’s law, , gives (A.4) and for our split monopole geometry F_{rφ} = 0. Outside the gap the flow is stationary (∂_{t} = 0), and charge conservation, , yields ∂_{r}(Σj^{r}) = 0 in the forcefree limit, so that the radial electric current is conserved along magnetic surfaces: Σj^{r} = J_{0} = const. From Eq. (A.4) we obtain (A.5) outside the gap. This conserved current must be flowing through the gap. Thus, the induction equation, Eq. (A.4), reduces to (A.6) Frame dragging couples the toroidal and poloidal components of the wave field: (A.7) (A.8) From the latter relations we find , where λ ≪ r is the characteristic wavelength of the oscillations, and likewise for , thus those fields can be neglected.
Appendix B: Monte Carlo scheme
B.1. Compton scattering
The Compton scattering opacity is given by (B.1) where γ is the Lorentz factor of the particle before scattering, I_{s}(ϵ_{s}) is the intensity of the ambient radiation field defined in Eq. (18), is the target photon energy in the rest frame of the particle, τ_{0} is the fiducial optical depth given in Eq. (19), and (B.2) is the Klein−Nishina cross section measured in units of σ_{T}. The opacity κ_{c}(γ) is used to draw scattering events.
Fig. B1. Cumulative probability distribution for selecting scattered photon energy, for different values of ϵ_{r} = γ(1+β)ϵ_{min}, with ϵ_{r} = 0.02 (solid lines), 0.2 (dashed lines), 2 (dotdashed lines), and 20 (longshort dashed lines). The black curves delineate the exact distribution, and the red curves the analytic approximation, Eqs. (B.9) and (B.10). 
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 in a single scattering can be expressed as (B.3) where A is a normalization coefficient, defined such that , , , with β holding for , and , with β holding for , and
(B.4) is the spectral density of target radiation field in the rest frame of the particle. The cross section for scattering of a photon of energy and direction to a final energy is given explicitly by (B.5) in terms of the differential crosssection (B.6) where is the angle between incident and scattered photons in the particle’s rest frame, and is obtained by substituting in Eq. (B.6).
To save computing time we use an approximate method to sample the probability density . We first note that to a good approximation . Then, substituting Eqs. (18), (B.4), and the latter relation into Eq. (B.3), we obtain (B.7)
The cumulative distribution is given by (B.8)
To shorten the notation we denote ϵ_{r} = γ(1 + β)ϵ_{min}. Approximate analytic expressions for F are obtained in terms of ϵ_{r} as follows:
For ϵ_{r} > 0.3, (B.10) where . The analytic distribution, Eqs. (B.9) and (B.10), is invertible and can be readily used to randomly select the energy . It is plotted in Fig. B.1 (red lines) for different values of ϵ_{r} = γ(1 + β)ϵ_{min}, and compared (black lines) with the exact expression computed by numerically integrating Eq. (B.3).
Drawing the scattered photon direction. Once the energy is selected it can be used to draw the direction of the scattered photon. The cross section for a scattering of target photon of energy into the final state is
(B.11) where (B.12) The conditional probability distribution reads (B.13) with , and given by Eq. (B.3). Due to the extremely strong beaming we can safely assume that all particles arrive from direction with a very small scatter. The latter relation then yields , and from Eq. (B.12) we have (B.14) Approximating the cross section by and recalling that , we find (B.15) The cumulative distribution, , explicitly given by
where if and F_{0} = 1 if , is used to randomly select .
B.2. Pair production
The full pair production cross section (measured in units of σ_{T}) is given by (B.16) where β_{cm} is the speed of the electron and positron in the center of momentum frame, given by (B.17) and ϵ_{γ}, ϵ_{s} are the dimensionless energies of the annihilating photons. The threshold energy for pair production through the interaction of a target photon with a gamma ray of energy ϵ_{γ} is given by the condition β_{cm} = 0 or ϵ_{th} = 2/(1 − μ)ϵ_{γ}. The pair production opacity reads (B.18)
References
 Abramowski, A., Acero, F., Aharonian, F., et al. 2012, ApJ, 746, 151 [NASA ADS] [CrossRef] [Google Scholar]
 Acciari, V. A., Aliu, E., Arlen, T., et al. 2009, Science, 325, 444 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Aharonian, F., Akhperjanian, A., Beilicke, M., et al. 2003, A&A, 403, L1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Aharonian, F., Akhperjanian, A. G., BazerBachi, A. R., et al. 2006, Science, 314, 1424 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Aharonian, F. A., Barkov, M. V., & Khangulyan, D. 2017, ApJ, 841, 61 [NASA ADS] [CrossRef] [Google Scholar]
 Albert, J., Aliu, E., Anderhub, H., et al. 2008, ApJ, 685, L23 [NASA ADS] [CrossRef] [Google Scholar]
 Aleksic, J., Ansoldi, S., Antonelli, L. A., et al. 2014, Science, 346, 1080 [NASA ADS] [CrossRef] [Google Scholar]
 Barkov, M. V., Aharonian, F. A., & BoschRamon, V. 2010, ApJ, 724, 1517 [NASA ADS] [CrossRef] [Google Scholar]
 Bicknell, G. V., & Begelman, M. C. 1996, ApJ, 467, 597 [NASA ADS] [CrossRef] [Google Scholar]
 Blandford, R. D., & Levinson, A. 1995, ApJ, 441, 79 [NASA ADS] [CrossRef] [Google Scholar]
 Blandford, R. D., & Znajek, R. L. 1977, MNRAS, 179, 433 [NASA ADS] [CrossRef] [Google Scholar]
 Bromberg, O., & Levinson, A. 2009, ApJ, 699, 1274 [NASA ADS] [CrossRef] [Google Scholar]
 Cerutti, B., Werner, G. R., Uzdensky, D. A., & Begelman, M. C. 2013, ApJ, 770, 147 [NASA ADS] [CrossRef] [Google Scholar]
 Giannios, D., Uzdensky, D. A., & Begelman, M. C. 2010, MNRAS, 402, 1649 [NASA ADS] [CrossRef] [Google Scholar]
 Globus, N., & Levinson, A. 2014, ApJ, 796, 26 [NASA ADS] [CrossRef] [Google Scholar]
 Hirotani, K., & Pu, H.Y. 2016, ApJ, 818, 50 [NASA ADS] [CrossRef] [Google Scholar]
 Hirotani, K., Pu, H.Y., Lin, L. C.C., et al. 2016, ApJ, 833, 142 [NASA ADS] [CrossRef] [Google Scholar]
 Hirotani, K., Pu, H.Y., Lin, L. C.C., et al. 2017, ApJ, 845, 77 [NASA ADS] [CrossRef] [Google Scholar]
 Katsoulakos, G., & Rieger, F. M. 2018, ApJ, 852, 112 [NASA ADS] [CrossRef] [Google Scholar]
 Levinson, A. 2000, Phys. Rev. Lett., 85, 912 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Levinson, A., & Rieger, F. 2011, ApJ, 730, 123 [NASA ADS] [CrossRef] [Google Scholar]
 Levinson, A., & Segev, N. 2017, Phys. Rev. D, 96, 123006 [NASA ADS] [CrossRef] [Google Scholar]
 Levinson, A., Melrose, D., Judge, A., & Luo, Q. 2005, ApJ, 631, 456 [NASA ADS] [CrossRef] [Google Scholar]
 Lin, L. C.C., Pu, H.Y., Hirotani, K., et al. 2017, ApJ, 845, 40 [NASA ADS] [CrossRef] [Google Scholar]
 Mahlmann, J. F., CerdáDurán, P., & Aloy, M. A. 2018, MNRAS, 477, 3927 [NASA ADS] [CrossRef] [Google Scholar]
 Nathanail, A., & Contopoulos, I. 2014, ApJ, 788, 186 [NASA ADS] [CrossRef] [Google Scholar]
 Neronov, A., & Aharonian, F. A. 2007, ApJ, 671, 85 [NASA ADS] [CrossRef] [Google Scholar]
 Owen, F. N., Eilek, J. A., & Kassim, N. E. 2000, ApJ, 543, 611 [NASA ADS] [CrossRef] [Google Scholar]
 Rieger, F. M. 2011, Int. J. Mod. Phys. D, 20, 1547 [NASA ADS] [CrossRef] [Google Scholar]
 Stawarz, L·., Aharonian, F., Kataoka, J., et al. 2006, MNRAS, 370, 981 [NASA ADS] [CrossRef] [Google Scholar]
 Timokhin, A. N. 2010, MNRAS, 408, 2092 [NASA ADS] [CrossRef] [Google Scholar]
 Timokhin, A. N., & Arons, J. 2013, MNRAS, 429, 20 [NASA ADS] [CrossRef] [Google Scholar]
All Figures
Fig. 1. Snapshots from a typical simulation of spark gap dynamics. Shown are the pair and photon densities (upper panels), electric flux (middle panels), 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. 

In the text 
Fig. 2. Snapshots of electron (blue line), positron (red line), and photon (yellow line) spectra computed during the run presented in Fig. 1. 

In the text 
Fig. 3. Dependence of the pair multiplicity ⟨κ⟩, pair energy ⟨γ⟩, and gammaray energy ⟨ϵ⟩ averaged over the simulation domain as functions of the input parameters τ_{0} (right) and J_{0} (left). 

In the text 
Fig. 4. Dependence of quiescent IC, curvature, and kinetic luminosities leaving the outer boundary of the box (r_{out} ≈ 4) on the fiducial opacity τ_{0}. 

In the text 
Fig. 5. Gammaray light curve produced by the gap discharge. The red line corresponds to IC emission, the blue line to curvature emission, and the black dashed line to the sum of both components. 

In the text 
Fig. B1. Cumulative probability distribution for selecting scattered photon energy, for different values of ϵ_{r} = γ(1+β)ϵ_{min}, with ϵ_{r} = 0.02 (solid lines), 0.2 (dashed lines), 2 (dotdashed lines), and 20 (longshort dashed lines). The black curves delineate the exact distribution, and the red curves the analytic approximation, Eqs. (B.9) and (B.10). 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.