Physics of Pair Producing Gaps in Black Hole Magnetospheres

In some low-luminosity accreting supermassive black hole systems, the supply of plasma in the funnel region can be a problem. It is believed that a local region with an unscreened electric field can exist in the black hole magnetosphere, accelerating particles and producing high-energy gamma-rays that can create e± pairs. We carry out time-dependent self-consistent 1D PIC simulations of this process, including inverse-Compton scattering and photon tracking. We find a highly time-dependent solution where a macroscopic gap opens quasi-periodically to create e± pairs and high-energy radiation. If this gap is operating at the base of the jet in M87, we expect an intermittency on the order of a few rg/c, which coincides with the timescale of the observed TeV flares from the same object. For Sagittarius A* the gap electric field can potentially grow to change the global magnetospheric structure, which may explain the lack of a radio jet at the center of our galaxy.


Introduction
The rotational energy of a spinning Kerr black hole can be electromagnetically extracted to launch powerful jets (Blandford & Znajek 1977). This process relies on the magnetic field being frozen into the plasma. However, the jet funnel is continuously evacuated dynamically, which can be understood by looking at a nearly force-free, steady monopolar magnetosphere ( Figure 1). 5 There are two light surfaces; particles moving in the strong magnetic field with a Larmor radius º  r r GM c g L 2 slide along the field line like beads on a wire and are flung outward toward infinity or inward toward the event horizon through the two light surfaces (Y. Yuan et al. 2018, in preparation). In between there is a stagnation surface, located at the maximum of an effective potential (compare with Takahashi et al. 1990;Broderick & Tchekhovskoy 2015). If disk material cannot cross field lines to enter the jet, and if there are not enough MeV photons to produce pairs through γ-γ collision, plasma density in the jet funnel will eventually become insufficient to conduct the current required by a force-free magnetosphere. As a result, E P will be induced, accelerating particles to high energies, producing gamma-ray photons and initiating a pair cascade (e.g., Blandford & Znajek 1977). This process is mostly possible in low accretion rate systems such as M87 and Sgr A * as they lack a reliable source of copious MeV photons.
The charge and current densities required to maintain a force-free magnetosphere have a few important properties for the monopole solution: (1) the poloidal current is constant along the flux tube; (2) the four-current is spacelike everywhere; (3) there exists a null surface where the zero angular momentum observer (ZAMO) measured charge density = J 0 0 ( Figure 1). The null surface has been regarded as a potential location for gap formation, analogous to pulsar outer gaps (e.g., Cheng et al. 1986;Beskin et al. 1992;Hirotani & Okamoto 1998). Meanwhile, other candidate places exist, e.g., the stagnation surface (Broderick & Tchekhovskoy 2015). Whether there is a preferred gap location has yet to be tested using kinetic simulations.
The magnetospheric gap has also been invoked to explain the fast γ-ray variability observed from some active galactic nuclei (AGNs; e.g., Levinson & Rieger 2011;Aleksić et al. 2014;Aharonian et al. 2017;Katsoulakos & Rieger 2018), but the gap physics has been treated under oversimplified assumptions, e.g., vacuum scaling, or, at best, steady-state models (Broderick & Tchekhovskoy 2015;Hirotani et al. , 2017Ford et al. 2017;Levinson & Segev 2017). In this work, we perform a detailed study of the microphysics and dynamics of the gap using radiative PIC simulations. Our approach is similar to Timokhin & Arons (2013) in spirit; however, we employ a novel setup appropriate for the black hole magnetosphere. Furthermore, the inverse-Compton (IC) and photon-photon pair production processes critical for the cascade are included self-consistently, as we explain in the following sections.

Numerical Setup
We would like to model the physical system described in Section 1 using the simplest physics possible while capturing all the salient features, namely, the existence of a null surface, a stagnation point, and two light surfaces. Spacetime correction to the equations of motion for particles and fields are responsible for these effects, while the magnitudes of these corrections are actually small compared to the electromagnetic forces. Therefore, we model the flux tube using a flat space model while trying to capture these features using different means.
In the flat spacetime model, we simply have 1D Maxwell equations (e.g., Chen & Beloborodov 2013;Timokhin & Arons 2013): where we take the background j B to be constant and use a spatial profile of r co similar to the GR background charge density (see the third row of Figure 2). We include the effect of the light surfaces in 1D using a model inspired by the light cylinder effect of a rotating neutron star. Consider the classic Michel monopole solution (Michel 1973): The field line rotates at an angular velocity Ω and becomes mostly toroidal outside the light cylinder where q W = r c sin , or β f =−1. For any particle outside the light cylinder, the corotation velocity is superluminal. The particle has to slide along the field line outward to keep the total velocity less than c.
We can easily derive the equation for the 1D constrained motion of a particle along a monopolar field line: where p r is the canonical momentum and is the Lorentz factor of the particle. One can immediately see that regardless of the value of p r , g < | | p mc 1 r , so when  b f 1 2 , v r >0 and the particle is only allowed to move in one direction. In our numerical simulations, we model the light surfaces in exactly the same way. We assume a profile for β f that is linear and antisymmetric across the center of the simulation box, reaching±1 at the light surfaces that are located at 0.05 and 0.95 of the simulation box. This way, the inner light surface is simply a mirror of the outer light surface; the midpoint of the box will be our "stagnation surface." We use a simplified version of the code Aperture developed by Alex Chen as a part of his PhD thesis (Chen 2017). The code evolves the 1D equations listed above, but keeps the charge-conserving current deposit scheme proposed by Esirkepov (2001). This ensures that Gauss's law is satisfied at all times if it is satisfied initially, so we only need to evolve the first part of Equation (1). Comparison between background and numerical values of ρ and j in the third row of Figure 2 confirms excellent charge conservation over time.

Choice of Units and Scales
In our flat spacetime model, we take the background current density j B to be constant, which naturally defines a plasma frequency and skin depth Thus, we measure time using wp 1 and length using λ p . Also choosing j B as the unit of current, the electric field will be measured by p w corresponds to a voltage drop of mc 2 over a single λ p , where the tilde denotes a dimensionless quantity. In this set of units, we naturally have the unit of energy being mc 2 and momentum being mc. We define pair multiplicity as  = + + -( ) n n ec j B . The profile of r co adds another parameter since it varies on the length scale of r g . The computational domain needs to accommodate several r g since we would like to include both light surfaces. For the physical systems we are interested in, e.g., M87, λ p /r g is on the order of 10 −8 ( Table 1). It is extremely difficult to have such scale separation in a PIC simulation. Therefore, we rescale this ratio, keeping l  r g p , and develop a semi-analytical model to infer what would happen at physical parameters.

Mechanism for Pair Production
The dominant mechanism for pair production in the black hole magnetosphere is the collision of high-energy photons with the low-energy photons from the disk. The high-energy photons come mostly from IC scattering of the background soft photons by energetic leptons. We carry out the full radiative transfer including IC scattering and the subsequent photonphoton collision assuming both happen on the same background photon field. We assume a soft photon spectral that cuts off at  min and extends up to ∼0.1 MeV. We use the Monte Carlo method to sample the photon energy from a single IC event, then we compute its free path by drawing from an exponential distribution with a mean gg ℓ . We track this photon until it is converted to an e ± pair at the end of its free path. When the photon is not energetic enough to convert within the box, we do not track it, but still cool the particle as if it emitted the photon.
The mean free path gg ℓ is energy dependent.
The modeling of the IC process introduces several new numerical parameters: the spectral index α, the peak soft photon energy  min , and a characteristic free path for IC scattering ℓ IC .  min sets the energy scale of the discharge, while ℓ IC puts a new length scale into the problem. The inferred characteristic values of these parameters are listed in Table 1. In our rescaling of the problem, we focus on the optically thick regime and ensure parameter ordering of l   ℓ r p g IC .

Simulation Results
We start from a plasma-filled initial condition where E=0 and r r = co . The initial pair multiplicity is 2, and all particles start at rest. Initially, a small plasma-scale electric field develops to help the current to flow, but since the box is leaky on both ends, plasma multiplicity drops over time. This happens fastest where dρ co /dx is largest. When   1, an  Notes. a L s , u s , and n s are soft photon luminosity, energy density, and number density at a few (∼3) r g from the black hole: electric gap opens locally to accelerate particles to high Lorentz factors and, subsequently, to initiate pair creation, screening this gap and launching macroscopic bunches of pair plasma to both directions. Screening of the electric field creates oscillations similar to those described by Levinson et al. (2005). Eventually, when the pairs are advected out of the light surfaces and multiplicity drops again, the same cycle is initiated. In the full length of one simulation, we are able to see several cycles of the gap formation. We also tried starting with a vacuum electric field, but obtained the same solution. Figure 2 shows one such gap cycle, where we used α=1.2, l = ℓ 10 p IC , lr 10 g p 4 , and  = -10 min 6 . It is the third time the gap develops in the simulation. As multiplicity drops from the previous cycle of pair creation, the system tries to maintain a macroscopic region as large as possible with   1 by drawing plasma from the side. When plasma flow cannot sustain this state, a gap opens quickly over the whole region where ~1. As a result, in all of our simulations, the gap size h∼r g , and it depends weakly on all parameters. The gap shown develops around the null surface, but it is not a guaranteed feature. It tends to develop wherever local multiplicity drops below unity, which can be anywhere due to plasma flow and delayed conversion of photons.
The photon spectrum shown in Figure 2 is not to be confused with the observable one. Due to limits of computational power, we only track photons that convert to pairs within the box, so the shown spectrum should be interpreted more as an absorption spectrum. In fact, most of the dissipated power in the gap goes into radiation that leaves the box, only a small fraction of it converting into e ± pairs. The peak multiplicity from the gap is usually ~10.
There are two well-defined spectral peaks for the particle energy shown in Figure 2. When the gap is screened, the lowenergy peak is a spectral break where the IC cooling becomes ineffective,l L cool,IC . When the gap is open, another spectral peak arises at higher energy. These are the primary particles accelerated in the gap, and the peak energy γ p is controlled by the gap electric field and IC cooling.

Physics of the Gap
Consider a region in the magnetosphere where plasma multiplicity  = 1 and j=j B at t=0. Electrons and positrons have to be counter-streaming at speed of light to provide the current, so the number density of each species evolves as Assuming r co varies linearly across this region with a slope r =k d dx j r c B g co , we see that the current decreases over time if > kj 0 , In particular, the timescale for the decrease of current depends on the spatial scale over which r co varies. As a result, the electric field at the center of the gap increases as t 2 : The gap starts to be screened when enough photons emitted by the primary particles convert to pairs within the gap. During the characteristic time t, a primary particle goes through a number of scatterings with target photons of energy ò (   min ), generating γ-rays at energy . Among these γ-rays, a number of κ convert to pairs within time t: , 1 3 p ph 2 2 IC 2 2 min 2 which turns out to be independent of ò. For the gap to be screened, κ needs to be on the order of 1-10 (we find that κ∼10 reproduces well the results of our production runs). For most parameter regimes of interest, when the screening happens the primary particles have short enough cooling lengths such that the electrostatic acceleration is balanced by the IC loss, so γ p is determined by Using Equations (11), (13), and (14), we can then obtain the peak electric field and γ p can be calculated from Equation (14). From E P and gap size h∼r g , we can estimate the maximum gap power We expect most of this power to be radiated away in gamma-rays. Figure 3 shows that for all runs below the Klein-Nishina regime, there is good agreement between the analytical model and the measured scaling from the simulations. However, the above calculation no longer holds well if primary particle energy approaches the Klein-Nishina regime:   g0.1 p min , which happens at a relatively small optical depth  t = r ℓ g 0 I C a few hundred. In that case, IC cooling becomes less efficient, and our argument for radiation-balanced acceleration breaks down. We expect γ p and the gap power to be much higher than our model would predict, which is what we see in the simulations. The gaps in this regime tend to be larger, but are still screened quasi-periodically as long as τ 0 >a few.
In the limit where  ℓ r g IC or  t 1 0 , we found that it is increasingly difficult to screen the gap, which develops to encompass the whole domain (this agrees with the discussion in, e.g., Broderick & Tchekhovskoy 2015). Particles are accelerated into the deep Klein-Nishina regime, where  g 1 p min and   E B. In this limit, we expect significant changes of the magnetospheric structure due to the gap, possibly killing the jet structure, and the 1D approximation we employed in this Letter is no longer appropriate.

Scaling to Real Systems
The most relevant systems where the spark gap might exist are low-luminosity AGNs. Of particular interest are M87 and Sgr A * (discussed extensively by, e.g., Broderick & Tchekhovskoy 2015), which are also primary targets for the Event Horizon Telescope. We list the physical parameters inferred from observation in Table 1, as well as predictions from our model. The observational parameters are based on Broderick & Tchekhovskoy (2015). For M87 we expect it to be well described by our model, and indeed, the predicted g~3 10 p 7 is well below the KN regime. The typical gamma-ray photons that are produced by these primary particles will be in the range 0.1 to a few TeV; depending on the actual location of the gap, some of them will escape the outer light surface. This coincides with the observed energy range of the TeV flares from M87 (Abramowski et al. 2012). Our time-dependent gap model also predicts time variability of several r g /c, which for M87 would be about ∼1 day, again coinciding with the observed timescale of the flares. However, the total gap power predicted by our model is at best only consistent with the quiescent state and is too low for the flares. Whether this mechanism can explain the origin of M87 flares will be investigated in a future paper. The typical luminosity we found here is much lower than, say, Broderick & Tchekhovskoy (2015), and that is mainly due to our dynamic model predicting a much lower  E . For Sgr A * , however,  g >0.1 p min , and we are in the Klein-Nishina regime. In this case, we expect the actual primary particle energy to be higher, and L gap might become comparable to L jet . As a result, the simplistic 1D approximation we adopted in this Letter is no longer applicable, as this gap should be able to significantly affect the global magnetosphere structure. This potentially can explain the lack of an apparent jet structure from the center of our galaxy. To properly treat this regime a global magnetospheric simulation will be needed.

Discussion
We have presented self-consistent 1D simulations of pair cascade in a magnetized plasma within the black hole magnetosphere. Informed by the numerical results, we developed a semi-analytical model for the electric gap, providing an estimate for gap power in systems that are optically thick to IC scattering.
Traditionally, the studies of the discharge problem in the BH magnetosphere were often based on a steady gap model around the null surface, drawing analogy to the outer gap model in pulsar magnetospheres (e.g., Ptitsyna & Neronov 2016;Ford et al. 2017;Hirotani et al. 2017). We have shown through numerical simulations that the physical conditions for such models are never realized: there is no steady gap. The electric field develops in a macroscopic region as a response to insufficient current and starts to be screened when pair cascade kicks in. The maximum electric field is controlled by the efficiency of pair creation and can be much lower than the corotation electric field (as assumed by Broderick & Tchekhovskoy 2015) if IC optical depth is high. Moreover, the gap can develop anywhere depending on plasma flow, not necessarily at the null surface or stagnation surface.
The gap opening and screening cycle we see are similar to those reported by Timokhin & Arons (2013). However, they did not investigate the situation when r co changes sign inside the simulation box. Their gaps all developed from the boundary of the pair-creating region because charge starvation occurs at the boundary when no pairs can flow from outside the region. They also only considered single-photon absorption, which is a different pair-creation mechanism than what is considered here.
We did not include the GR correction to the particle equations of motion. Instead, the GR effect is entirely captured by the varying background charge density r co and the presence of an inner light surface. Without GR, both features will be absent. We think this is an appropriate simplification that allows us to focus on the electrodynamics and microphysics. A proper general relativistic set of equations could, in principle, be implemented, as was recently done by Levinson & Cerutti (2018). However, they report an overall quasi-steady state that is different from what we observe. We believe the main difference is the treatment of light surfaces, and they focus on a low optical depth regime  r ℓ 10 g IC . The logical extension of the results in this Letter is to look at how the global structure of the magnetosphere will interact with the gap, especially when the gap power becomes comparable to the jet power.
We thank Matthew Kunz and Alexander Tchekhovskoy for very helpful discussions, and we thank Andrei Beloborodov, Anatoly Spitkovsky, and the anonymous referee for their constructive comments on the manuscript. Y.Y. acknowledges support from the Lyman Spitzer, Jr. Postdoctoral Fellowship awarded by the Department of Astrophysical Sciences at Princeton University. A.C. acknowledges the support of NASA grant NNX15AM30G. H.Y. is supported by NSERC and in part by the Perimeter Institute for Theoretical Physics. Research at Perimeter Institute is supported by the Government of Canada through the Department of Innovation, Science and Economic Development Canada and by the Province of Ontario through the Ministry of Research, Innovation and Science. The code used in this paper is available athttp://doi. org/10.5281/zenodo.1340677 under GPLv2.0.