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 unscreened electric field can exist in the black hole magnetosphere, accelerating particles and producing high energy gamma-rays that can create $e^{\pm}$ 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^{\pm}$ 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 $r_g/c$, which coincides with the time scale 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
It is widely believed that powerful jets from Active Galactic Nuclei (AGN) are launched by rapidly spinning black holes through the Blandford-Znajek process (Blandford & Znajek 1977). When the black hole is threaded by a strong, coherent magnetic field, the rotating spacetime near the event horizon drags the field lines into rotation, resulting in a significant amount of Poynting flux flowing out along the magnetic field. The energy flow comes from the rotation of the black hole, and eventually gets dissipated further downstream along the jet, giving rise to the multi-wavelength emission we observe from jets. The launching and acceleration of Poynting flux dominated jets have been demonstrated in general relativistic magnetohydrodynamic (GRMHD) simulations (e.g., McKinney et al. 2012). However, it turns out that in the jet funnel, plasma tends to be evacuated either inward to the black hole or outward to infinity, and some source of continuous plasma supply is needed to carry the underlying electric current, which cannot be self-consistently captured in the MHD framework.
In reality, when the plasma supply runs low in the jet funnel, electrostatic gaps -regions with unscreened electric field -could develop and accelerate particles to high enough energies to initiate an e ± pair cascade, which would replenish the plasma needed in the jet circuit, similar to what happens around pulsars (e.g., Blandford & Znajek 1977;Beskin et al. 1992;Hirotani & Okamoto 1998). But unlike pulsars where the pair producing photons are mainly generated by curvature radiation, near the black hole the dominant process is inverse Compton (IC) scattering of soft photons emitted by the accretion disk. The electrostatic gaps are probably most important in low luminosity AGN. In high luminosity AGN, the accretion disk may already produce many MeV photons that could collide with each other and convert to pairs in the jet funnel. If these pairs are already enough to conduct the current, then a gap does not need to open up (e.g., Levinson & Rieger 2011).
The magnetospheric gap has been considered as possible candidate site for producing the rapid γ-ray variability seen from certain AGN, especially radio galaxies (e.g., Levinson & Rieger 2011;Aleksić et al. 2014;Aharonian et al. 2017;Katsoulakos & Rieger 2018). The most promising source seems to be M87 (Katsoulakos & Rieger 2018). It has a massive black hole with mass M ≈ 6.5 × 10 9 M (Event Horizon Telescope Collaboration et al. 2019), and an estimated jet power L jet ∼ 10 44 erg s −1 (e.g., Broderick & Tchekhovskoy 2015). Its very high energy (VHE, E > 100 GeV) γ-ray emission shows strong variability on timescales as short as a day (Abramowski et al. 2012). During the γ-ray flare, the VHE luminosity can reach 10 42 erg s −1 , and the flares have been observed to be accompanied by increased radio emission from the nucleus (Acciari et al. 2009), strongly suggesting particle acceleration in the vicinity of the black hole. In order to see whether the gap scenario can reproduce all the features of the γ-ray flares, a detailed study of the gap physics is needed.
So far in the literature there are quite a few works on the gap physics; most of them employ steady state/electrostatic models (Beskin et al. 1992;Hirotani & Okamoto 1998;Broderick & Tchekhovskoy 2015;Ptitsyna & Neronov 2016;Hirotani et al. , 2017Ford et al. 2018;Levinson & Segev 2017). However, the gap development and screening may be an intrinsically time dependent process (Broderick & Tchekhovskoy 2015;Levinson & Segev 2017), similar to the polar cap gap in pulsars (Timokhin & Arons 2013). Recently, Levinson & Cerutti (2018) used 1D general relativistic particle-in-cell (GRPIC) simulations to study the discharge of a gap that is depleted of particles at the beginning. They found that except for the initial transient, the long term evolution reaches a quasisteady state where no macroscopic gap is seen and pairs are generated everywhere randomly. Parfrey et al. (2019) carried out 2D GRPIC simulations of the jet launching process by immersing a rapidly spinning Kerr black hole in a uniform magnetic field. They used artificially simplified pair injection schemes, and found that the resulting plasma distribution depends sensitively on the pair creation physics.
In our previous paper, we also carried out 1D local simulations of the gap using a flat spacetime approximation (Chen et al. 2018, hereafter CY18). We found the gap development and screening to be a quasiperiodic process, a result qualitatively very different compared to Levinson & Cerutti (2018). We focused on the high opacity scenario where inverse Compton scattering mainly happens in the Thomson regime, and constructed a simple analytic model to extrapolate numerical simulations to M87 and Sgr A*. It was found that for M87, the gap power is not enough to power the observed TeV flares. In this paper, we properly implement General Relativity and a more complete numerical treatment for inverse Compton scattering. We extend our study into the deep Klein-Nishina (KN) regime and attempt to find an empirical scaling relation to extrap-olate our simulation results. We will first describe our method comprehensively in Section 2, then present the main results in Section 3. We will comment on the implication of these results for M87 in particular in Section 4, and finish with a discussion of potential caveats and future extensions in Section 5.

Equations of Motion
We consider a nearly force-free magnetosphere around the black hole, with open field lines threading the event horizon, in a more or less axisymmetric, quasi-steady state. Somewhere along the field line, gaps may open up when particle density runs low. If the spatial and temporal scale of the gap is much smaller than the global characteristic scale (∼ a few r g ≡ GM/c 2 , the gravitational radius), the gap dynamics can be well described using 1D approximation along the field line.
In our numerical simulations, we take the deviation from the background force-free configuration to be our dynamic variable, and solve the equations along a 1D flux tube labeled by the flux function ψ(r, θ) (see Appendix A for a description of the force-free configuration). The basic equations have a relatively simple form in the 3 + 1 formalism (e.g., Komissarov 2004). Here, the Kerr metric can be written as (we use Greek indices for 4-vectors, and Latin indices for spatial 3-vectors) where α is the lapse function and β i is the shift vector. It is convenient to use Boyer-Lindquist coordinates as γ ij is diagnal: γ ij = diag(g rr , g θθ , g φφ ) = diag(Σ/∆, Σ, A sin 2 θ/Σ), where Σ = r 2 + a 2 cos 2 θ, ∆ = r 2 − 2M r + a 2 , A = (r 2 + a 2 ) 2 − ∆a 2 sin 2 θ. Meanwhile, and β i only has φ component: where ω = −g 0φ /g φφ = 2M ar/A is the angular velocity of the zero angular momentum observer (ZAMO). To alleviate the divergence near the event horizon, we use the tortoise coordinates which is a transformation from the Boyer-Lindquist coordinates by letting dξ = dr/∆, such that where r ± = M ± √ M 2 − a 2 are the radius of the event horizons given by ∆ = 0. We can see that ξ → −∞ as r → r + . The field equations that we solve in (ξ, t) are: where the factor K 1 ≡ ∆ √ γ/(∂ θ ψ) accounts for field line divergence along the flux tube. We use the third equation to solve for the current due to particle movement, therefore Gauss's law (equation (6)) is automatically satisfied as long as it holds initially. We always set ρ = ρ ff initially and D ξ = 0, which trivially satisfies the Gauss's law. We observe excellent charge conservation as even after millions of time steps charge density is always driven back to ρ ff . For particles, we assume they are well magnetized and can only move freely along magnetic field lines. In other words, their motion is constrained to 1D, and u ξ = u r /∆ completely determines the full 4-velocity: where B ξ = B r /∆, and Ω is the angular velocity of the field line. We can then write down the Lagrangian in terms of coordinate ξ and velocity v ξ = dξ/dt: where u 0 is written in terms of v ξ as here The canonical momentum is and the Euler-Lagrange equations are where E ξ = −∂ ξ Φ. The equations can also be written purely in terms of p ξ , using In our field update equations D ξ is the dynamic variable, while Equation (14) needs E ξ . We can obtain the latter through This set of equations naturally reproduces the inner and outer light surfaces, outside which particles flow either to the horizon or to infinity. We place the boundaries of our simulation box to be slightly outside the two light surfaces, so that all particles near the boundary can only outflow.
Photons are not confined to field lines so they obey the general equations of motion. In our example, the background field lines on the poloidal plane is more or less radial. It is a good approximation to neglect θ component of the photon 4 velocity. Therefore the equations are where γ ξξ = γ rr /∆ 2 .

Radiative Transfer
We assume that the background soft photon field is isotropic in the ZAMO frame, and model the scattering event (IC and γγ pair production) in this frame. For simplicity, we assume that the number density of the soft photons in the ZAMO frame is a power law independent of the radius: We have chosen the normalization such that n 0 has the dimension of the total number of photons per volume. In addition, we consider the regime where the number density of the background soft photon field is much larger than those emitted by the particles in the gap, so that only scatterings on the background photons are important while nonlinear effects can be neglected. The rate of scattering between an electron and the background photons in the ZAMO frame is (e.g., Blumenthal & Gould 1970) Here σ IC, is the full KN cross section (e.g., Rybicki & Lightman 1986): where x is the incident photon energy in the electron rest frame measured in units of the electron rest mass: and γ = αu 0 is the energy of the electron in the ZAMO frame. At each time step, we determine the scattering probability P IC using the following relation: where ∆t is the global time step, and ∆t its corresponding value in ZAMO frame. If a scattering event occurs, we draw the scattered photon energy based on the IC spectrum (e.g., Jones 1968) where E 1 is the energy of the scattered photon in units of the initial electron energy: Since the photon is not tied to the field line, we compute its u ξ and u φ when it is created, and follow its trajectory using the photon equations of motion in section 2.1. For γγ pair production, the total cross section for photons of energy E, and incident angle θ is (Gould & Schréder 1967) where βc is the electron (positron) velocity in the center of mass frame, The threshold requires s > 1. A high energy photon of energy E traversing an isotropic background soft photon field with a number spectrum n( ) in the ZAMO frame will have an absorption probability per unit timê When a photon converts to a pair, we set the direction of motion of the particles to be along local B field, and the photon energy is evenly split to the pair in the local ZAMO frame. We tabulate all the cross sections and spectra and use a binary search table look-up to find the correct rates. This is an efficient procedure and we are able to process more than 10 9 scatterings per second on a single GPU.

Units and Scales
We adopt a numerical unit system where r g is the unit of length, and r g /c is the unit of time. Magnetic field in numerical units is measured in terms of cyclotron frequency. We use tilde to denote dimensionless quantities. For example,B = eBr g /(m e c 2 ).
There are three main parameters of the problem. B 0 is the background magnetic field strength, which determines the background current j ff and charge density ρ ff . It also sets a characteristic plasma skin depth which leads toB 0 = (r g /λ p ) 2 . In other words,B 0 sets the ratio between macroscopic and microscopic length scales. The peak of the soft photon spectrum min is the second parameter, which sets the energy scale of the radiative transfer. Particles with γ 0.1/˜ min interact with the soft photon field in the Thomson regime, while those with γ > 0.1/˜ min fall into the KN regime and suffer from reduced scattering cross section. Gammaray photons with energy 0.1/˜ min ˜ 1/˜ min have the smallest free path to γ-γ collision.
The third parameter τ 0 is the characteristic optical depth to inverse Compton scattering in the Thomson regime. It scales with the number density n 0 of the soft photon field: where IC = 1/(n 0 σ T ). τ 0 determines the efficiency of IC scattering as well as γ-γ collision. If any gap develops in the magnetosphere, it is difficult to screen it before it grows to h ∼ r g /τ 0 , simply because that is the minimum mean free path the high energy photons need to travel before creating pairs. If τ 0 is close to unity, then we expect the gap to at least grow to scale of r g . We use the code Aperture for simulations in this paper. This version of the code is highly optimized on GPUs, and provides a unified framework for simulating the magnetospheres of neutron stars and black holes. The resolution for the runs in this paper range from 16384 to 196608 grid cells, with up to 200 million particles and photons, evolved over several million time steps. All of these runs were carried out on a single GPU, either a Tesla P100 or a Quadro P6000.

Time-dependent Gap
In almost all cases, our simulations show quasiperiodic gap opening and screening, very similar to the results of CY18 . Typically the system will undergo an initial transient phase where a large gap opens, filling the box with plasma, then settle down to the quasiperiodic state. One such cycle is shown in Figure 1. Here a good criterion for forming a gap is based on the multiplicity, defined as the ratio between the actual number of charge carriers in the gap and the minimum amount needed to provide the current 1 : Whenever the local pair multiplicity drops below unity, parallel electric field starts to grow, accelerating particles and producing photons. It is remarkable that, except for the initial transient which depends on initial condition, the gap always grows from the null surface where ρ ff = 0. The quasi-periodic behavior of the gap can be seen from the integrated dissipation rate in the flux tube, L = αg rr D r j r ( √ γ/∂ θ ψ)dr. The time-evolution of this gap power can be interpreted as the light curve, and is plotted in Figure 2, where the power is normalized to the Poynting flux density in the same flux tube measured in the force-free simulation, L 0 ≈ 0.0124B 2 0 , which we can also identify as the jet power. Three runs of differentB 0 and˜ min but sameB 0˜ min and τ 0 are plotted, and they have near identical time evolution and gap dissipation. As will be discussed in section 3.2, we find that in general gap power and screening timescale only scale with the productB 0˜ min , but insensitive to the values ofB 0 and˜ min themselves.
In addition to highly time-dependent gap power, we observe the same variability reflected in the photon flux near the outer light surface. We keep track of about 5% of all the particles and photons in the simulation, and look at all the photons within 0.5r g of the outer light surface. The total energy as well as the spectrum of these outgoing photons vary according to gap activity. Figure 3 compares the spectra of the high state (when photon flux is highest) and low state (when photon flux is lowest). Both spectra form a power law, with the high state having a harder spectral index. The photon spectrum tends to extend a bit above 0.1/˜ min , at which point it is strongly absorbed. The lower end of the spectrum is somewhat artificial, since we need to remove some low energy photons in the simulation in order to make it computationally feasible. In this particular run, all photons with energy below 10 3 m e c 2 are removed at creation. However, photons created above Outgoing photon spectrum at high (blue curve) and low (green curve) states defined in section 3.1. k is the slop of the power law index in EdN/dE. The low energy cutoff is mostly due to our removal of photons with energies below ∼ 10 3 mec 2 at creation. This particular run hasB0 = 10 8 and˜ min = 10 −5 .
this energy threshold are allowed to propagate and can be redshifted to lower energies if they travel from near the gap to the outer boundary. Even with lower energy photons removed, the energy escaping from the box is still dominated by photons, suggesting most of the gap power goes into radiation. In our simulations the multiplicity of pairs escaping the box is usually between 10-50. However, most of the photons emitted in the gap do not convert in the box, and will continue to produce more particles if they travel outwards for a much larger distance, undergoing a post gap cascade which has been discussed by e.g. Broderick & Tchekhovskoy (2015). We expect the terminal multiplicity to be much higher than what we get in the local simulation.

Dependence on Parameters
In the high τ 0 or small IC case, e.g. τ 0 100, IC scatterings mainly happen in the Thomson regime. The behavior of the gap is identical to what was reported in CY18. Despite the high optical depth, the gap is macroscopic with a size h ∼ 0.1-1r g , sometimes even larger during the initial transient. Particles are predominantly radiation reaction limited, and gap screening happens when energies of the upscattered photons are increased such that γγ ∼ h. The fiducial optical depth, together with˜ min , sets the energy scale at which gap screening occurs. The scaling relation derived in CY18, namely their equation (15), still holds in the high τ 0 regime. Using the unit convention of this paper, the same relation can be written as (neglecting the order 1 constants) Notice that this ratio depends only on the product of B 0 and˜ min , not individually. With fixed τ 0 the gap size is also almost constant with respect to the product B 0˜ min . Although this equation should only be valid in the Thomson regime, we find that the dependence on the productB˜ min is more universal. This is because particle energy only has one characteristic scale, which is γ 0 = m e c 2 / min . Writing γ =γγ 0 , electric field E =ẼB 0 , particle number density n =ñm e c 2 /(e 2 r 2 g ), the main equations become (neglecting the geometric factors) ∂rẼ ∼ñ + −ñ + , and ∂tγ ∼ (B 0˜ min )Ẽ, which only depend on the productB 0˜ min .
When τ 0 100, particles are accelerated into the KN regime. Since the photon mean free path for γ-γ collision is smallest when the photon energy is comparable to m e c 2 /˜ min , the gap will be at least as large as the minimal γγ , which is γγ ∼ 5 IC = 5r g /τ 0 . We observe that the gap size scales very weakly with parameters, only becoming larger whenB 0˜ min is so small that the gap voltage ehDr max is comparable to m c c 2 /˜ min . Surprisingly, when τ 0 3, the particles are accelerated too far into the KN regime that the gap is no longer screened in our simulations.
The gap screening mechanism in the KN regime is different from what was described in CY18 for the Thomson regime. In the latter case IC loss is efficient so that most of the highest energy particles in the gap are radiation reaction limited, and the highest energy photons have the smallest free path. It is the highest energy photons that eventually produce the gap-screening pairs, therefore we focused on the kinematics of the highest energy particles to derive the gap scaling relation. However, in the KN regime this is no longer the case. Cooling becomes inefficient when the particle Lorentz factor reaches 1/˜ min . Particles are accelerated to much higher energies than the radiation reaction limit, and the first secondary pairs are created by photons with γ ∼ m e c 2 /˜ min , since they have the shortest mean free paths. However, γγ is large even for these photons, so there is a significant delay between the emission of these photons and the screening of the gap, during which time the primary particles are accelerated to even higher energies, e.g. γ 10/˜ min . The screening process in the KN regime depends strongly on the acceleration history of the primary particles in the gap. When γγ r g , the gap is global and the electric field grows so quickly that not enough photons with energy ∼ m e c 2 /˜ min are emitted. As a result, the gap can no longer be screened. Figure 4 shows the dependence of the gap power oñ B˜ min and τ 0 individually. In all the runs shown on the plot, we observe the same time-dependent gap described in 3.1, and measure the peak gap power averaged over several cycles. It can be seen that in general the power scales close to L ∼ (B˜ min ) −1 with fixed τ 0 . However, we have to note here that at highB˜ min the measured gap power could be smaller than it should be, due to numerical heating which will be discussed in section 3.3, since numerical heating tends to reduce the need for gap acceleration. Also in Figure 4 one can see that when τ 0 is around 10 the gap power scales weakly with τ −1/2 0 . However, above τ 0 ∼ 50 which is around the transition between Thomson and KN regime the power scales closer to τ −0.9 0 . This being steeper than the scaling ofDr max suggests that the gap size decreases with increasing τ 0 as well.

Dependence on Resolution
We discovered that with increasingB˜ min , numerical simulations become more and more challenging. Figure 5 shows light curve comparisons between simulations with identical parameters, but different numbers of particles per cell, or different resolutions, atB˜ min = 10 4 . With a smaller number of particles per cell or lower resolution, the simulation eventually settles to a quasisteady state not unlike what was reported by Levinson & Cerutti (2018). However, with more particles per cell or higher resolution, we observe the same kind of quasiperiodic opening and screening of the gap as shown in Section 3.1. Generally we find that when the system settles to a quasi-steady state, not only is the dissipation much less than the peak of the cyclic solution, the average power is also about a a factor of ∼ 1.5 less than the average power of the cyclic solution.
This resolution dependent behavior mainly happens with highB andB˜ min , where large electrostatic oscillations can accelerate particles close to γ ∼ 1/˜ min , thus sustain pair creation without forming a time-dependent gap. To minimize this numerical effect, one needs both high spatial resolution and large number of particles per cell. The spatial grid should ideally resolve the plasma skin depth λ p , otherwise the plasma will be heated to Lorentz factors γ T ∼ (∆x/λ p ) 2 , at which point the hot plasma skin depth can be resolved. If this energy is already capable of producing pairs, we will simply have a hot thermal plasma in pair creation equilibrium. Small number of particles per cell will lead to high discrete particle noises in the simulation, which also translates to larger plasma oscillations, broadening the spread of distribution function in momentum space. If this allows more particles to reach pair creation energies then again a quasi-steady state can be achieved. We believe this is the main reason CY18 and Levinson & Cerutti (2018) found qualitatively different solutions to the same problem. Both were using a similar numerical resolution, but Levinson & Cerutti (2018) were attempting to use parameters corresponding toB˜ min ∼ 5 × 10 6 while we found thatB˜ min = 10 4 is already very challenging for the typical numerical resolution used in this kind of simulations. For this reason, we also expect that the gap power we see in our highestB˜ min simulations may not be the most reliable (e.g. the last point in the left panel of Figure 4, whereB˜ min = 3×10 4 ), as numerical effects reduce the amount of E generated by the system.

IMPLICATIONS FOR M87
M87, as a low luminosity AGN, is one of the most promising candidates where this process can operate. One can obtain a reasonable estimate for the magnetic field from the estimated jet power, B 0 ∼ 200 G (e.g., Levinson & Rieger 2011;Broderick & Tchekhovskoy 2015), which translates toB 0 ∼ 1.2 × 10 14 . However, the parameters for the soft photon field in the vincinity of the central black hole is a lot more uncertain. The SED for the radio flux seems to peak around 1.2 meV, which translates to˜ min ∼ 2 × 10 −9 , but the soft photon spectrum near the horizon may well be very different from what is observed. Assuming all photons are coming from near the horizon, the upper bound for photon energy density is u s ∼ 0.1 erg/cm 3 which implies that τ 0 ∼ 5 × 10 3 . But since most photons are produced in the disk, u s near the horizon could be orders of magnitude lower.
In CY18 we found that if we take τ 0 ∼ 5 × 10 3 then the gap power in the Thomson regime would be L ∼ 3 × 10 39 erg s −1 . If we assume much lower fiducial optical depth, e.g. τ 0 ∼ 10 as in our simulations, and simply scale L/L 0 according to Figure 4 to realistic value ofB˜ min ∼ 2.4 × 10 5 , we can expect the power dissipated in the gap to be L ∼ 10 −3 L 0 . The jet power of M87 is usually estimated to be L 0 ∼ 10 43 -10 44 erg s −1 , so our estimated gap power should be around L gap ∼ 10 40 -10 41 erg s −1 , somewhat smaller than, but not entirely inconsistent with the reported power of the VHE flares from the same object (e.g. Abramowski et al. 2012).
We expect the spectrum from the gap to be similar to what is shown in Figure 3, which is a power law extending up to ∼ 0.1/˜ min = 5 × 10 7 , followed by an exponential cutoff. This translates to energy up to around 25 TeV, also consistent with the observed flare spectrum. The exact slope of the spectrum will likely be subject to further modification due to scattering with the ambient soft photon field and radiation from the secondary pairs.

DISCUSSION
We presented results from first-principle 1D GR PIC simulations of the pair discharge process in the vincinity of supermassive black holes. In general we observe a highly time-dependent solution where an electric gap opens quasi-periodically near the null surface. The variability time scale of this activity is usually on the order of several r g /c, and the outgoing photon spectrum forms a power law with index that depends on the gap activity.
In this work we ignored synchrotron and curvature radiation. In reality, when τ is order unity, photon free path is significant compared to the macroscopic length scale. As they convert to a pair, the particles should have a relatively large pitch angle with respect to the magnetic field. It will be damped quickly through synchrotron radiation, which should dissipate a sizable fraction of the particle energy and produce detectable radia-tion. The highest energy particles that have γ ∼ 10/˜ min can approach curvature radiation reaction limit, and curvature loss can be a significant energy sink as well. We plan to carry out a more complete spectral analysis in the future, including the effect of synchrotron and curvature radiation in this framework, which could produce a unified spectrum from X-ray to VHE gamma rays out of the pair producing gaps.
In this paper we considered γ-γ collision to be the only pair creation mechanism. As pointed out by Petropoulou et al. (2019), when particles are in the deep KN regime, triplet pair production may play a role in the screening of the gap. Due to the highly uncertain nature of the soft photon field in the vincinity of the supermassive black hole of M87, there is a parameter regime where triplet pair production could be the main mechanism that regulates the gap physics. This can potentially further increase the gap luminosity and could be powering the strongest flares that we see up to L ∼ 10 42 erg s −1 . A radiation module that handles triplet pair creation is already implemented in Aperture and we expect to carry out a separate study for the importance of this process in M87, and whether this produces any observable signatures on the spectrum. This paper studied the gap microphysics in 1D, whereas in many runs the gap evolves to have a size comparable to r g , at which point the 1D approximation is not completely valid. Working 2D axisymmetric GR PIC codes now exist (Parfrey et al. 2019). It is an important and interesting question whether the results presented in this paper are robust when the 1D approximation is relaxed. However, we caution any future attempt at solving this problem in 2D full GR that resolution is important. Significant rescaling and some abstraction of the radiative physics must be employed, and extreme care must be taken in order to draw any meaningful comparisons with observations.