© ESO 2018

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 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; 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 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; 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), 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 Sects. 2 and 3 below. We find (Sect. 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 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⁻⁵ 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 steady-state 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 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.

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²+a²−2r_gr, Σ = r²+a² cos²θ, A = (r²+a²)²−a²Δ sin²θ, and r_g = GM/c² = 1.5×10¹⁴ M₉ cm denotes the gravitational radius, and M = 10⁹ 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 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(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 four-velocity components measured by a ZAMO, here denoted by , the particle Lorentz factor , and the three-velocity 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_sL_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_ec² 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₀/hc, and can be used to define the fiducial optical depth (19) For we have l_s ≃ 10⁻² ϵ_minτ₀.

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_ec² units) is (22) where λ_c = ħ/m_ec 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₀ = ηB_H = η10³B₃ 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 τ₀ > 1. In a transparent gap, τ₀ < 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⁷ during the initial discharge (where γ ≃ γ_max) to about 10⁵ 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⁹ the characteristic energy of curvature photons becomes ϵ_c < 10³, 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 (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 gamma-ray 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, Σjr, normalized by the global magnetospheric current J0. 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⁹M₉ 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³ B₃ 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₀) defined in Eq. (A.5), and the fiducial optical depth (τ₀) defined in Eq. (19). The results presented below were computed using the following canonical choice of parameters (M₉ = 1, B₃ = 2π, θ = 30°, ã = 0.9, ϵ_s,min = 10⁻⁸, 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₀ and τ₀.

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 second-order finite-difference scheme.

The electric field is evolved in time using Eq. (6), such that (26) where δt = tⁿ⁺¹ − tⁿ = 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 (27) Assuming that , we obtain (28) The photon position at time tⁿ⁺¹ 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 right-hand 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 three-velocity at time step tⁿ. 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 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., (38) 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 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⁻⁹−10⁻⁷ for a range of opacities τ₀ = 1 to τ₀ = 10. For the size of our simulation box, h/r_g ≈ 3, at least 10⁴ cells are needed to resolve the skin depth for τ₀ = 10. Simulations with τ₀ < 10 showed good convergence for a total N_ξ = 16384 cells. For τ₀ = 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 (τ₀ = 1) to 10 (τ₀ = 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.

4. Results

4.1. Overall evolution

We ran simulations for a grid of models characterized by different values of the parameters τ₀ and J₀. Quite generally, we find that the evolution of the system depends on the value of τ₀, but is practically independent of J₀ (Fig. 3). The prime role of J₀ 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 long-term 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 mono-energetic 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 quasi-steady 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 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 τ₀ (left) and global current J₀ (right). As can be seen, the multiplicity increases roughly linearly with τ₀, while the average pair and photon energies decrease with increasing τ₀. On the other hand, those quantities are essentially independent of J₀. The overall gamma-ray luminosity emitted from the gap following the prompt phase depends weakly on τ₀ (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 gamma-ray energy ⟨ϵ⟩ averaged over the simulation domain as functions of the input parameters τ0 (right) and J0 (left).

Fig. 4. Dependence of quiescent IC, curvature, and kinetic luminosities leaving the outer boundary of the box (rout ≈ 4) on the fiducial opacity τ0.

Fig. 5. Gamma-ray 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⁷. 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_GJr_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 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)².

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⁻⁵ L_BZ.

5. Applications to M87

M87 exhibits TeV emission with a luminosity of L_TeV ∼ 10⁴⁰ erg s⁻¹ 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⁴³ to a few times 10⁴⁴ erg s⁻¹. Assuming that the jet is powered by the BZ mechanism implies L_TeV/L_BZ ∼ 10⁻⁴. The SED exhibits a peak in the sub-mm band at ν_peak ≃ 10¹² Hz, corresponding to ϵ_max ≃ 10⁻⁸ in our parametrization, with a bolometric luminosity of L_b≲10⁴¹ erg s⁻¹. 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⁹ M_⊙ this implies l_s < 10⁻⁷ 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⁻⁸, 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⁻¹ and τ₀ < 10².

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⁻⁵ 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 τ₀ = 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⁻⁸m_ec² (corresponding to a minimum frequency of about 10¹² 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⁻⁵ 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 steady-state models with extremely low accretion rates. For example, Hirotani et al. (2017) obtained values in the range 10⁻⁴−10⁻¹ for accretion rates in the range from 10⁻⁴ to 5×10⁻⁶ 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 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 (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⁻³ Eddington, can conceivably produce a TeV luminosity of about 10³¹ erg s⁻¹ 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 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 . 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 force-free limit, so that the radial electric current is conserved along magnetic surfaces: Σj^r = J₀ = 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, τ₀ 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 (dot-dashed lines), and 20 (long-short 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 cross-section (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.9) where

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₀ = 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

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, Σjr, normalized by the global magnetospheric current J0. 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 gamma-ray energy ⟨ϵ⟩ averaged over the simulation domain as functions of the input parameters τ0 (right) and J0 (left).
In the text

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

	Fig. 5. Gamma-ray 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 (dot-dashed lines), and 20 (long-short 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