# Synthetic Images of Magnetospheric Reconnection-Powered Radiation around Supermassive Black Holes

Benjamin Crinquand, Benoît Cerutti, Guillaume Dubus, Kyle Parfrey, Alexander Philippov

► To cite this version:

Benjamin Crinquand, Benoît Cerutti, Guillaume Dubus, Kyle Parfrey, Alexander Philippov. Synthetic Images of Magnetospheric Reconnection-Powered Radiation around Supermassive Black Holes. Physical Review Letters, 2022, 129, 10.1103/PhysRevLett.129.205101. insu-03860254

HAL Id: insu-03860254

https://insu.hal.science/insu-03860254

Submitted on 18 Nov 2022

**HAL** is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire **HAL**, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Synthetic Images of Magnetospheric Reconnection-Powered Radiation around
Supermassive Black Holes

Benjamin Crinquand●,¹,²,* Benoît Cerutti,¹ Guillaume Dubus,¹ Kyle Parfrey,¹ and Alexander Philippov⁴,⁵

¹Université Grenoble Alpes, CNRS, IPAG, 38000 Grenoble, France

²Department of Astrophysical Sciences, Peyton Hall, Princeton University, Princeton, New Jersey 08544, USA

³School of Mathematics, Trinity College Dublin, Dublin 2, Ireland

⁴Center for Computational Astrophysics, Flatiron Institute, 162 Fifth Avenue, New York, New York 10010, USA

⁵Department of Physics, University of Maryland, College Park, Maryland 20742, USA

<figure><img src="image_2.png" /></figure>

(Received 20 February 2022; revised 15 August 2022; accepted 20 September 2022; published 10 November 2022)

Accreting supermassive black holes can now be observed at the event-horizon scale at millimeter wavelengths. Current predictions for the image rely on hypotheses (fluid modeling, thermal electrons) which might not always hold in the vicinity of the black hole, so that a full kinetic treatment is in order. In this Letter, we describe the first 3D global general-relativistic particle-in-cell simulation of a black-hole magnetosphere. The system displays a persistent equatorial current sheet. Synthetic radio images are computed by ray-tracing synchrotron emission from nonthermal particles accelerated in this current sheet by magnetic reconnection. We identify several time-dependent features of the image at moderate viewing angles: a variable radius of the ring, and hot spots moving along it. In this regime, our model predicts that most of the flux of the image lies inside the critical curve. These results could help promote understanding of future observations of black-hole magnetospheres at improved temporal and spatial resolution.

DOI: 10.1103/PhysRevLett.129.205101

*Introduction*.—At least part of the formidable energy released by active galactic nuclei (AGN) may originate from the rotational energy of the central black hole, extracted by the means of strong magnetic fields threading its event horizon [1]. These magnetic fields are sustained by currents flowing in an accretion flow, which in the case of a low-luminosity AGN is hot, geometrically thick, and collisionless. Recently, for the first time, the Event Horizon Telescope (EHT) collaboration was able to spatially resolve the immediate vicinity of the low-luminosity AGN M87*. The image shows a circular asymmetric ring encompassing a depression in brightness [2]. From the high polarization fraction and brightness temperature of the radiation, we know that the light captured by the EHT is an optically thin synchrotron emission emitted by mildly relativistic leptons. These measurements are consistent with the scenario of dynamically important poloidal magnetic fields at event horizon scales, thus favoring the “magnetically arrested disk” (MAD) scenario for the accretion flow [3,4].

However, the location of the main emission sites, as well
as the mechanisms causing plasma heating and particle
acceleration, are still poorly constrained. A widely adopted

approach has been to interpret the emission as synchrotron radiation by thermal [e.g. [5,6]] or nonthermal [7–9] electrons, energized for example by subgrid kinetic turbulent dissipation or magnetic reconnection. General-relativistic (GR) magnetohydrodynamic (MHD) simulations can be educated with prescriptions for electron heating used to model the local synchrotron emissivity, which allows one to synthesize images when coupled with GR ray-tracing algorithms. In this framework, the image should prominently display a black-hole shadow determined by the properties of the spacetime [10–12]. It is still unclear, though, to what extent the features of the image are related to the astrophysical details of the emission [13–16].

Recently, several GRMHD studies of MAD disks
(Ripperda *et al.* [17] and Chashkina *et al.* [18] in 2D and
Ripperda *et al.* [19] in 3D) have reported dynamical
transitions from an accreting state, during which magnetic
flux is accumulated onto the horizon, to a highly magnetized
quasi-force-free state, which is characterized by a stark drop
in the accretion rate and the existence of a thin equatorial
current sheet. This configuration favors fast dissipation of
the magnetic energy accumulated in the magnetosphere
through magnetic reconnection in the current sheet, trigger-
ing efficient particle acceleration and high-energy non-
thermal radiation [20,21]. Thus, it has been suggested
that gamma-ray flares from radio galaxies (such as M87*)
could naturally occur when the central supermassive black
hole is in this quasi-force-free state, hereafter labeled “high-
energy flaring state.”

<footnote>Published by the American Physical Society under the terms of
the Creative Commons Attribution 4.0 International license.
Further distribution of this work must maintain attribution to
the author(s) and the published article's title, journal citation,
and DOI.</footnote>
In this regime, with the innermost zone almost depleted of electron-ion plasma, predictions for the image based on GRMHD do not apply, whereas kinetic simulations can accurately capture the physics and dynamics of the highly magnetized and collisionless plasma. Global kinetic particle-in-cell (PIC) simulations of black-hole magnetospheres have only recently begun to be carried out [22–26]. Crinquand et al. [24] (hereafter, C21) previously studied the dynamics and high-energy gamma emission of the equatorial current sheet that develops generically within the ergosphere from a large class of ordered magnetic configurations [27]. Conversely, in this Letter, we characterize the low-energy emission at radio wavelengths from a black hole in a high-energy flaring state, embedded in a quasi-force-free magnetosphere, and discard any contribution from the accretion flow.

**Simulations.**—We use the GRPIC code GRZeltron to simulate magnetospheric pair plasma. The background spacetime is described by the Kerr metric with a dimensionless spin parameter $a = 0.99$. We use spherical horizon-penetrating Kerr-Schild coordinates $(r, \theta, \varphi)$, and model self-consistent plasma injection by including inverse Compton scattering and photon-photon annihilation [23]. We have run both a 2D and a 3D simulation, using the setup described in C21 with initially paraboloidal magnetic field lines switching polarity at the equator. The 3D simulation has a resolution of $1024(r) \times 256(\theta) \times 512(\varphi)$, whereas the 2D one has $1536(r) \times 1024(\theta)$. The simulation domain is $r \in [0.9r_h, 10r_g]$, $\theta \in [0.1, \pi - 0.1]$, $\varphi \in [0, 2\pi]$, with $r_g$ the gravitational radius of the black hole and $r_h$ the radius of the event horizon. Our fiducial 2D and 3D runs do not have synchrotron cooling. Because of a higher numerical cost, the values of the dimensionless input parameters are slightly less realistic in 3D than in 2D: the fiducial Larmor radius $r_L$ of particles is $r_L = 10^{-4}r_g$ in 3D instead of $2 \times 10^{-6}r_g$ in 2D, whereas the background soft photon energy is $\epsilon_0 = 10^{-2} m_e c^2$ instead of $10^{-3}$. In both cases, the scale separation is necessarily smaller than that of a realistic AGN. M87*, for example, is characterized by a ratio $r_L/r_g \approx 10^{-14}$, which is currently unfeasible numerically. We have also confirmed that 2D simulations are not sensitive to the precise initial distribution of plasma and photons.

The general behavior of the 3D simulation is very similar to that of the 2D one. The magnetic field develops a strong toroidal component, with opposite signs above and below the equator. Thus, a current sheet is needed to support the associated discontinuity. Pairs are accelerated by nonideal electric fields in the equatorial plane, triggering high-energy photon upscattering and even stronger pair creation. The densities are significantly larger in the current sheet than in the polar regions. The initial current sheet quickly fragments into flux ropes (see Fig. 1). The 3D evolution of the reconnecting current sheet is governed by the tearing instability, which produces flux ropes elongated in the

<figure><img src="image_3.png" /><figcaption>FIG. 1. Snapshot of the total plasma density, normalized by the fiducial Goldreich-Julian density, at $t = 40r_g/c$. Magnetic field lines are represented as blue solid lines. The black surface marks the event horizon.</figcaption></figure>

direction orthogonal to the reconnecting field, and the drift-kink instability, which can corrugate the sheet along the reconnecting field [20,28]. Besides, the flux ropes are orbiting around the black hole and are constantly being sheared by differential rotation, and plasma can be pushed away from the reconnecting sheet and ejected along the separatrix (the last open field line connected to the black hole). These competing effects drive the current sheet toward a state of twisted and disordered flux ropes. It should be pointed out that our simulation ran for around $40r_g/c$, which is too short to observe the full cyclic behavior observed in C21, during which a giant magnetic island, formed between the current sheet and the separatrix, replenishes the black hole with magnetic flux. Although these flux ropes could lose coherence for longer integration times, we note that the simulation has run in steady state for about $20r_g/c$, which is longer than the characteristic growth time of the tearing and drift-kink instabilities which could disrupt it [29] and longer than the orbital period at the event horizon ($4\pi/\omega_h \sim 14r_g/c$, with $\omega_h = ac/2r_h$).

**Ray tracing.**—We use the public ray-tracing code GEOKERR [30] to perform geodesic integration of synchrotron photons from the emission point to a distant observer's screen, in postprocessing (see the Supplemental Material [31], which includes Ref. [32]). When analyzing the 2D simulation, we focus on the time dependence on the image. By virtue of axisymmetry, deposition of photons on the observer's screen is independent of the observer's azimuthal viewing angle. By contrast, when processing the 3D simulation, we effectively take a snapshot of the simulation but allow for different azimuthal viewing angles $\varphi_{\text{obs}}$. This way, we can separately assess the effects of time variability (inaccessible to the short 3D simulation) and nonaxisymmetry of the current sheet (inaccessible to the 2D axisymmetric simulation) on the image.
<figure><img src="image_4.png" /><figcaption>FIG. 2. Time-averaged synchrotron images from the 2D simulation. The flux is in arbitrary units. The three columns show three different viewing angles $\alpha_{\text{obs}}$, measured with respect to the spin axis of the black hole. The white cross marks the position of the black hole. The blue dashed line denotes the critical curve for $a = 0.99$ at each viewing angle.</figcaption></figure>

We only show the contribution of the current sheet. Indeed, synchrotron emission by particles accelerated in the polar gaps is unlikely to be observable. Only a small fraction of the Blandford-Znajek electromagnetic power is dissipated in these gaps, which probably achieve plasma densities too low to yield significant emission. Besides, unlike particles accelerated in the current sheet, particles emitted in the polar regions have low pitch angles, implying a comparatively low synchrotron emission. We have also run a 2D simulation with synchrotron cooling turned on, where the synchrotron radiation-reaction force is rescaled, so that the characteristic cooling time of an electron should be a few simulation time steps [33]. In that case, the morphology of the image is practically unchanged whereas the polar cap emission is quenched. This is consistent with previous PIC studies of magnetic reconnection with synchrotron cooling [34–36]: particles are accelerated near X points where the magnetic field vanishes, allowing particles to reach high energies before escaping the re-connecting layer and cooling.

**Averaged images.**—Figure 2 shows images from our 2D simulation averaged over $100r_g/c$, which are very similar to the images from our 3D simulation averaged over the azimuthal viewing angle. The leftmost panel represents a view close to face on, which is the most relevant to the M87* system, assuming that the spin of the black hole is aligned with the large-scale jet [37]. In general, the image is made up of two distinct components: an inner and an outer ring. The rings’ brightness is higher on the left side of the image, as a result of relativistic beaming: this is mainly a geometrical effect due to the Kerr metric, rather than Doppler boosting. The inner ring is the direct image of the current sheet, whereas the outer ring is the lensed image. On each image, we have plotted the critical curve corresponding to each inclination, i.e., the set of points on the screen hit by photons which have orbited an arbitrarily large number of times around the black hole close to bound, marginally stable spherical orbits [38,39]. This curve matches well the outer ring emission, confirming its

interpretation as the lensed image of the current sheet. However, the lensed image accounts for only a small portion (around 10%) of the flux of the image. Our model hence predicts that in a high-energy flaring state, most of the flux of the image lies inside the critical curve. In this configuration, the black-hole shadow prediction breaks down because most of the emission comes from within the ergosphere and is not spherically distributed. We find that a major fraction of the current sheet radiation originates from a very compact zone, within 2–3 $r_g$. Both rings become thinner as the observing frequency increases, since the emission originates from the most energetic particles which have been accelerated deeper within the current sheet.

**Variability.**—The first row of Fig. 3 illustrates the time variability in the 2D axisymmetric simulation, for a viewing angle $\alpha_{\text{obs}} = 22.5^{\circ}$. Images are labeled using the observer’s Boyer-Lindquist coordinate time. At a certain time during the simulation, a bright circular feature appears between the inner and outer rings (top left panel). As time elapses, this feature shrinks in size until it overlaps the inner ring. Its typical lifetime is about $10r_g/c$. By computing the synchrotron emissivity of the plasma, we can relate this emission to a giant magnetic island accretion event (see C21). The second row of Fig. 3 illustrates the dependence on $\varphi_{\text{obs}}$ of images from the 3D simulation, for the same viewing angle. Both the inner and outer rings show non-axisymmetric features in the form of brighter spots. They are more pronounced on the inner ring, which is less impacted by gravitational lensing. The presence of hot spots remains visible even after blurring the image at a resolution $\approx 5r_g$. These hot spots can be traced back to the different equatorial orbiting flux ropes. They radiate relatively isotropically, so that every flux rope is visible regardless of $\varphi_{\text{obs}}$. This is why the positions of the hot spots along the ring shift continuously clockwise with increasing $\varphi_{\text{obs}}$ in Fig. 3. The hot spots appear to move counterclockwise along the ring, with a period $\sim 4\pi/\omega_{\text{h}} \sim 14r_g/c$, to a fixed observer at a given $\varphi_{\text{obs}}$.
<figure><img src="image_5.png" /><figcaption>FIG. 3. Top row: Snapshot of the image from the 2D simulation (leftmost panel) and differential maps illustrating the time evolution. The three maps show the difference between successive snapshots at later times and the first one. Bottom row: Image from the 3D simulation at a given azimuthal viewing angle $\varphi_{\text{obs}}$ (leftmost panel) and differential maps illustrating the dependence on $\varphi_{\text{obs}}$. The three maps show the difference between images at higher $\varphi_{\text{obs}}$ and the first one.</figcaption></figure>

Although the ratio between $r_g$ and the electron plasma skin depth is unrealistically small in these simulations, the hierarchical merger of magnetic islands should proceed until plasmons with macroscopic sizes are formed. The final number of synchrotron-emitting flux ropes, and hence of hot spots on the image, is determined by their characteristic escape time from the merging region [40]. The flux ropes are constantly reformed within the current sheet to compensate for their escape, along the separatrix or into the black hole. Cerutti and Giacinti [40] predict the final number of macroscopic plasmons to be ~$\pi/\beta_{\text{rec}}$, independently of the initial number of plasmons in the current sheet and the separation of scales. This is consistent with the number of hot spots we observe in the inner ring of the image.

**Observability in the radio band.**—The typical flux from the inner $10r_g$ in $M87^*$ at the frequency $\nu_0 = 230$ GHz is 1 Jy [2], translating into a radiated power $\nu F_\nu \approx 10^{40}$ erg s$^{-1}$. The total jet power $L_{\text{jet}}$ is likely to lie within the range $10^{42}-10^{45}$ erg s$^{-1}$ [41,42]. Assuming reconnection in the equatorial plane powers the emission, an upper bound on the total dissipated energy is $L_{\text{dis}} \sim \beta_{\text{rec}}L_{\text{jet}}$, with $\beta_{\text{rec}} \approx 0.1$ the collisionless reconnection rate [43]. The synchrotron cooling timescale $\tau_s$ is related to the inverse Compton cooling time $\tau_{\text{IC}}$ in the Thomson regime by $\tau_s/\tau_{\text{IC}} = U_{\text{IC}}/(B_0^2/8\pi) \approx 10^{-2}$, with $B_0 \approx 10-100$ G the typical magnetic field strength at horizon scale [4] and $U_{\text{IC}} \approx 10^{-2}$ erg cm$^{-3}$ the energy density of low-energy

photons in the inner zone (~$10r_g$) [44,45]. Since $\tau_s \ll \tau_{\text{IC}} \ll r_g/c$, most of the magnetic energy dissipated by magnetic reconnection, initially carried by nonthermal particles, is converted into synchrotron emission.

The spectrum of reconnection-accelerated particles in pair plasma with zero guide field can largely be described as a power law $dN/d\gamma \propto \gamma^{-p}$ [20,21,46]. The index $p$ decreases with increasing upstream magnetization $\sigma$, going from $p \approx 2$ for $\sigma \gtrsim 10$ to 1 at $\sigma \gg 10$, in the absence of synchrotron cooling [46]. As $p$ gets closer to 1, more and more energy is carried by high-energy particles. The synchrotron spectrum emitted by such a population of leptons also follows a power law, with a spectral energy density depending on frequency as $\nu F_\nu \propto \nu^{(3-p)/2}$. The normalization is such that the total radiated power matches $L_{\text{dis}}$. The peak of the total synchrotron spectrum cannot extend significantly beyond the burn-off limit at the photon energy $h\nu_{\text{rad}} \approx 20(\beta_{\text{rec}}/0.1)$ MeV, with $h$ the Planck constant ($\nu_{\text{rad}} \sim 10^{10}\nu_0$). The characteristic synchrotron frequency of the photons emitted by the highest-energy particles $\nu_s \sim \sigma^2 eB_0/m_e c$ lies above $\nu_{\text{rad}}$ for $\sigma \gtrsim 10^6$. The value of $\sigma$ in $M87^*$ is likely to be much larger than 1, but its exact value is uncertain.

If $\sigma \gtrsim 10^6$, reconnection in the bare current sheet occurs in the strong cooling regime. Assuming $p \lesssim 2$, the radiated power $\nu F_\nu$ at the radio frequency $\nu_0 = 230$ GHz approximately reads as $L_{\text{dis}}(\nu_0/\nu_{\text{rad}})^{(3-p)/2}$, yielding
$$ \nu F_{\nu} \approx 10^{33} \left( \frac{L_{\text{jet}}}{10^{44} \text{ erg s}^{-1}} \right) \left( \frac{h\nu_{\text{rad}}}{20 \text{ MeV}} \frac{230 \text{ GHz}}{\nu_{0}} \right)^{\frac{p-1}{2}} \text{ erg s}^{-1}, \quad (1) $$

where we have neglected gravitational redshift. Here, the most energetic photons have energies in the MeV range, so that they can produce pairs in the upstream [36,47]. We show in the Supplemental Material [31] that these secondary pairs also cool efficiently through synchrotron radiation, and emit in the radio band. Their total luminosity is $10^{39} \text{ erg s}^{-1}$, which is much larger than the estimate from Eq. (1), but is still outshone by the accretion disk, though marginally. Hence, the radio image would be dimmer in the flaring state.

On the other hand, $\sigma$ could be smaller than $10^6$ but still large enough for reconnection to occur in the relativistic regime. Although the multiplicities resulting from electromagnetic cascades in the polar spark gaps are not very high [23,48], additional mass loading from the accretion flow could reduce $\sigma$ with respect to the values defined by the plasma densities produced in these gaps. In that case, the particle spectrum retains its nonradiative shape approximately up to $\gamma \sim \text{few} \times \sigma$. Beyond, the spectrum steepens, the exact shape depending on the relative strengths of cooling and secondary acceleration mechanisms. To give a simple estimate of the power at radio frequencies, we assume that most of the power is carried by particles below $\sigma$. The energetic constraint on the synchrotron radio flux is then relaxed with respect to the $\sigma > 10^6$ case: effectively, $\nu_{\text{rad}}$ must be replaced by $\nu_s \sim \sigma^2 eB_0/m_e c$ in Eq. (1). Fixing the power-law index at $p_0 = 1$, a typical value in the case of $e^{\pm}$ nonradiative reconnection up to $\sigma$, we obtain $\nu F_{\nu} \approx 10^{40} (L_{\text{jet}}/10^{44} \text{ erg s}^{-1})(\sigma/10^3)^{p_0-3} \text{ erg s}^{-1}$. The radiation emitted by the current sheet in the high-energy flaring state could then reach levels similar to the quiescent state and be observable by the EHT. To obtain this estimate, we have neglected the contribution of additional acceleration mechanisms, beyond the impulsive acceleration of particles by the reconnected electric field, to the energy budget [49–51].

**Conclusion.**—In this letter, we have provided a physically motivated, first-principles model for the image of a lowluminosity AGN. The validity of our results is limited to describing a high-energy flaring state, in which the accretion flow has receded drastically, rather than the quiescent state observed by the EHT. However, even if the current sheet that we obtain from our GRPIC simulations results from idealized initial conditions, very similar configurations are reached by GRMHD simulations, which also include accretion physics. As long as the upstream magnetization is much larger than 1, the morphology of the image should only depend on the black-hole spin.

Since Sagittarius A* is expected to vary on timescales of several minutes, future observations will be able to observe changes in the morphology of the image within the EHT observation window. From our simulations, we expect the

radius of the bright ring to change with time, and possible hotspots to move along the ring, thereby providing a test of our model. In M87*, these hot spots should fully rotate along the ring in five days, whereas in Sagittarius A* it would take about five minutes.

The authors would like to thank the anonymous referees for insightful comments, as well as Bart Ripperda and Fabio Bacchini for helpful discussions. This project has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation program (Grant Agreement No. 863412). Computing resources were provided by TGCC and CINES under the allocation A0090407669 made by GENCI, and by Princeton Research Computing, a consortium of groups including the Princeton Institute for Computational Science and Engineering (PICSciE) and the Office of Information Technology’s High Performance Computing Center and Visualization Laboratory at Princeton University.

<footnote>*To whom all correspondence should be addressed. bcrinquand@princeton.edu</footnote>

<footnote>[1] R. D. Blandford and R. L. Znajek, Mon. Not. R. Astron. Soc. 179, 433 (1977).</footnote>

<footnote>[2] Event Horizon Telescope Collaboration et al., Astrophys. J. Lett. 875, L1 (2019).</footnote>

<footnote>[3] Event Horizon Telescope Collaboration et al., Astrophys. J. Lett. 910, L12 (2021).</footnote>

<footnote>[4] Event Horizon Telescope Collaboration et al., Astrophys. J. Lett. 910, L13 (2021).</footnote>

<footnote>[5] M. Mościbrodzka, H. Falcke, and H. Shiokawa, Astron. Astrophys. 586, A38 (2016).</footnote>

<footnote>[6] Event Horizon Telescope Collaboration et al., Astrophys. J. Lett. 875, L5 (2019).</footnote>

<footnote>[7] J. Davelaar, H. Olivares, O. Porth, T. Bronzwaer, M. Janssen, F. Roelofs, Y. Mizuno, C. M. Fromm, H. Falcke, and L. Rezzolla, Astron. Astrophys. 632, A2 (2019).</footnote>

<footnote>[8] N. Scepi, J. Dexter, and M. C. Begelman, Mon. Not. R. Astron. Soc. 511, 3536 (2022).</footnote>

<footnote>[9] C. M. Fromm, A. Cruz-Osorio, Y. Mizuno, A. Nathanail, Z. Younsi, O. Porth, H. Olivares, J. Davelaar, H. Falcke, M. Kramer, and L. Rezzolla, Astron. Astrophys. 660, A107 (2022).</footnote>

<footnote>[10] H. Falcke, F. Melia, and E. Agol, Astrophys. J. Lett. 528, L13 (2000).</footnote>

<footnote>[11] R. Narayan, M. D. Johnson, and C. F. Gammie, Astrophys. J. Lett. 885, L33 (2019).</footnote>

<footnote>[12] T. Bronzwaer, J. Davelaar, Z. Younsi, M. Mościbrodzka, H. Olivares, Y. Mizuno, J. Vos, and H. Falcke, Mon. Not. R. Astron. Soc. 501, 4722 (2021).</footnote>

<footnote>[13] S. E. Gralla, D. E. Holz, and R. M. Wald, Phys. Rev. D 100, 024018 (2019).</footnote>

<footnote>[14] K. Nalewajko, M. Sikora, and A. Różańska, Astron. Astrophys. 634, A38 (2020).</footnote>

<footnote>[15] F. H. Vincent, M. Wielgus, M. A. Abramowicz, E. Gourgoulhon, J. P. Lasota, T. Paumard, and G. Perrin, Astron. Astrophys. 646, A37 (2021).</footnote>
[16] W. Lockhart and S. E. Gralla, Mon. Not. R. Astron. Soc. **509**, 3643 (2022).

[17] B. Ripperda, F. Bacchini, and A. A. Philippov, Astrophys. J. **900**, 100 (2020).

[18] A. Chashkina, O. Bromberg, and A. Levinson, Mon. Not. R. Astron. Soc. **508**, 1241 (2021).

[19] B. Ripperda, M. Liska, K. Chatterjee, G. Musoke, A. A. Philippov, S. B. Markoff, A. Tchekhovskoy, and Z. Younsi, Astrophys. J. Lett. **924**, L32 (2022).

[20] L. Sironi and A. Spitkovsky, Astrophys. J. Lett. **783**, L21 (2014).

[21] F. Guo, H. Li, W. Daughton, and Y.-H. Liu, Phys. Rev. Lett. **113**, 155005 (2014).

[22] K. Parfrey, A. Philippov, and B. Cerutti, Phys. Rev. Lett. **122**, 035101 (2019).

[23] B. Crinquand, B. Cerutti, A. Philippov, K. Parfrey, and G. Dubus, Phys. Rev. Lett. **124**, 145101 (2020).

[24] B. Crinquand, B. Cerutti, G. Dubus, K. Parfrey, and A. Philippov, Astron. Astrophys. **650**, A163 (2021).

[25] A. Bransgrove, B. Ripperda, and A. Philippov, Phys. Rev. Lett. **127**, 055101 (2021).

[26] I. El Mellah, B. Cerutti, B. Crinquand, and K. Parfrey, Astron. Astrophys. **663**, A169 (2022).

[27] S. S. Komissarov, Mon. Not. R. Astron. Soc. **350**, 427 (2004).

[28] B. Cerutti, G. R. Werner, D. A. Uzdensky, and M. C. Begelman, Astrophys. J. **782**, 104 (2014).

[29] S. Zenitani and M. Hoshino, Astrophys. J. **670**, 702 (2007).

[30] J. Dexter and E. Agol, Astrophys. J. **696**, 1616 (2009).

[31] See Supplemental Material at http://link.aps.org/supplemental/10.1103/PhysRevLett.129.205101 for a description of the numerical techniques and a theoretical estimate of the secondary pairs' luminosity.

[32] B. Carter, Phys. Rev. **174**, 1559 (1968).

[33] B. Cerutti, A. A. Philippov, and A. Spitkovsky, Mon. Not. R. Astron. Soc. **457**, 2401 (2016).

[34] B. Cerutti, G. R. Werner, D. A. Uzdensky, and M. C. Begelman, Astrophys. J. **770**, 147 (2013).

[35] D. Kagan, E. Nakar, and T. Piran, Astrophys. J. **833**, 155 (2016).

[36] H. Hakobyan, A. Philippov, and A. Spitkovsky, Astrophys. J. **877**, 53 (2019).

[37] R. C. Walker, P. E. Hardee, F. B. Davies, C. Ly, and W. Junor, Astrophys. J. **855**, 128 (2018).

[38] C. T. Cunningham and J. M. Bardeen, Astrophys. J. **183**, 237 (1973).

[39] E. Teo, Gen. Relativ. Gravit. **35**, 1909 (2003).

[40] B. Cerutti and G. Giacinti, Astron. Astrophys. **656**, A91 (2021).

[41] A. E. Broderick, R. Narayan, J. Kormendy, E. S. Perlman, M. J. Rieke, and S. S. Doeleman, Astrophys. J. **805**, 179 (2015).

[42] M. A. Prieto, J. A. Fernández-Ontiveros, S. Markoff, D. Espada, and O. González-Martín, Mon. Not. R. Astron. Soc. **457**, 3801 (2016).

[43] L. Sironi, D. Giannios, and M. Petropoulou, Mon. Not. R. Astron. Soc. **462**, 48 (2016).

[44] A. E. Broderick and A. Tchekhovskoy, Astrophys. J. **809**, 97 (2015).

[45] EHT MWL Science Working Group et al., Astrophys. J. Lett. **911**, L11 (2021).

[46] G. R. Werner, D. A. Uzdensky, B. Cerutti, K. Nalewajko, and M. C. Begelman, Astrophys. J. Lett. **816**, L8 (2016).

[47] J. M. Mehlhaff, G. R. Werner, D. A. Uzdensky, and M. C. Begelman, Mon. Not. R. Astron. Soc. **508**, 4532 (2021).

[48] A. Y. Chen and Y. Yuan, Astrophys. J. **895**, 121 (2020).

[49] M. Petropoulou and L. Sironi, Mon. Not. R. Astron. Soc. **481**, 5687 (2018).

[50] H. Hakobyan, M. Petropoulou, A. Spitkovsky, and L. Sironi, Astrophys. J. **912**, 48 (2021).

[51] H. Zhang, L. Sironi, and D. Giannios, Astrophys. J. **922**, 261 (2021).