Radiative reconnection-powered TeV flares from the black hole magnetosphere in M87

Active Galactic Nuclei in general, and the supermassive black hole in M87 in particular, show bright and rapid gamma-ray flares up to energies of 100 GeV and above. For M87, the flares show multiwavelength components, and the variability timescale is comparable to the dynamical time of the event horizon, suggesting that the emission may come from a compact region nearby the nucleus. However, the emission mechanism for these flares is not well understood. Recent high-resolution general-relativistic magnetohydrodynamic simulations show the occurrence of episodic magnetic reconnection events that can power flares nearby the black hole event horizon. In this work, we analyze the radiative properties of the reconnecting current layer under the extreme plasma conditions applicable to the black hole in M87 from the first principles. We show that abundant pair production is expected in the vicinity of the reconnection layer, to the extent that the produced secondary pair plasma dominates the reconnection dynamics. Using analytic estimates backed by two-dimensional particle-in-cell simulations we demonstrate that in the presence of strong synchrotron cooling, reconnection can produce a hard power-law distribution of pair plasma imprinted in the outgoing synchrotron (up to a few tens of MeV) and the inverse-Compton signal (up to TeV). We produce synthetic radiation spectra from our simulations, which can be directly compared with the results of future multiwavelength observations of M87* flares.


2012)
, is non-negligible compared to the total jet power of 10 42 − 10 44 erg/s (e.g., Prieto et al. 2016). The rapid variability timescales of the order of a few light-crossing times of the event horizon associated with the TeV emission from M87 imply that the emitting region is of the order of the Schwarzschild radius, within the uncertainty of an unknown beaming factor (Harris et al. 2011). Additionally, the existence of an X-ray emission counterpart originating from nearby the nucleus for some of the events (Harris et al. 2011;Abramowski et al. 2012;Sun et al. 2018) suggests that the origin of the flares is in the vicinity of the central black hole. The conversion of magnetic energy into kinetic energy and radiation through magnetic reconnection is an obvious and often proposed candidate to power flares from magnetized plasma nearby the event horizon. The ex-act emission mechanism that powers these TeV flares is however not well understood. To study whether magnetic reconnection can power the observed radiation, it is essential to know the geometry and strength of the magnetic field nearby the event horizon. Event Horizon Telescope (EHT) observations of polarized radio emission (Event Horizon Telescope Collaboration et al. 2021) show that the accretion disk of M87 is likely to be in a particular scenario uncovered by general relativistic magnetohydrodynamics (GRMHD) simulations: the magnetically arrested state (Narayan et al. 2003;Tchekhovskoy et al. 2011). In this state, the magnetic flux on the black hole accumulates onto the event horizon until its pressure becomes dynamically important. When the magnetic pressure due to flux accumulation becomes comparable to the pressure of accreting gas, the flux is expelled through magnetic reconnection events in an equatorial current sheet that separates the northern and southern jet (Ripperda et al. 2020(Ripperda et al. , 2022. Based on their unambiguous occurrence in global GRMHD simulations, these reconnection events have been conjectured to power high-energy flares from nearby the event horizon of the black hole (Ripperda et al. 2020;Chashkina et al. 2021;Bransgrove et al. 2021;Ripperda et al. 2022).
GRMHD simulations, however, cannot predict the properties of the reconnection-powered emission. To understand the mechanisms that power the multiwavelength flaring emission, it is essential to rely on kinetic physics that describes particle acceleration from first principles. Moreover, the interaction between accelerated leptons and the emitted photons, which involves both classical (such as the radiation drag), and quantum effects (such as pair-production and Compton scattering), is unexplored in the regimes applicable for radiatively inefficient accretion flows in AGN.
Possible radiative mechanisms discussed in the literature to explain this emission often invoke curvature and inverse Compton (IC) radiation by leptons (electrons and positrons), accelerated to high energies in vacuum "gap" regions (Levinson & Rieger 2011;Broderick & Tchekhovskoy 2015;Chen & Yuan 2020;Crinquand et al. 2020). While being a promising source of TeV flares, gap models require an external trigger mechanism, which can substantially modify the properties of the background radiation field (Kisaka et al. 2022). Considerable changes can lead to significant particle acceleration and flaring emission due to the development of large regions with an unscreened electric field, that appear because of the change in the optical depth to γγ pair-production. Magnetic reconnection, on the other hand, self-consistently occurs as a consequence of the magnetic flux accumulation near the black-hole horizon during accretion. However, since the reconnection is a microscopic process dissipating the magnetic field energy and accelerating particles on scales much smaller than the horizon scale, it is unclear whether it can power the observed TeV flares both in terms of the total luminosity and in terms of the maximum photon energy.
In this Letter, we build a semi-analytic model for the flares in M87* and other AGN powered by the intermittent magnetic reconnection in the equatorial plane. We propose that short-duration high-luminosity TeV flares can indeed be produced via IC radiation powered by the upscattering of soft photons from the accretion flow by the electrons and positrons accelerated during magnetic reconnection. Our model further predicts a broadband spectrum ranging from radio to TeV. To quantify our findings, we conduct first-principles kinetic simulations of radiative reconnection in collisionless pair plasma, including dynamically important synchrotron radiation, and taking into account IC scattering in post-processing, for parameters applicable nearby the event horizon of the supermassive black hole in the center of M87.
2. M87* TEV FLARES: THEORY Any theoretical model aiming to explain the underlying physics behind the TeV emission is bound to incorporate three empirical characteristics of these flares: energetics, luminosity, and fast rise and fall timescales. First of all, in order to produce TeV emission it is necessary to at least have an efficient mechanism for particle energization able to accelerate leptons to energies E e ± 10 6 m e c 2 . In the context of reconnection -the process our discussion is centered upon -the characteristic particle energy depends on the available magnetic energy per particle, which in turn depends on how much plasma is supplied to the reconnection region. On top of that, the produced TeV flux, regardless of the underlying emission mechanism, has to be able to escape the magnetosphere. This fact puts additional upper limits for the optical depth to the Breit-Wheeler pair production process which can potentially inhibit the TeV photon flux, converting its energy into secondary e ± pairs (see, e.g., Levinson & Rieger 2011). Subsequently, given the relatively high luminosity of the flares (for M87 they reach 10 39 -10 40 erg/s) and the rapid rise/fall timescales (few days for M87), the model needs to invoke a very fast and efficient energy dissipation mechanism. As we discuss further in this section, radiative fast magnetic reconnection is an attractive solution in that regard, since it is able to dissipate the supplied magnetic energy rapidly and subsequently convert most of the liberated energy into the high-energy emission.
Recent GRMHD simulations of magnetically arrested black hole accretion flow (Ripperda et al. 2020;Chashkina et al. 2021;Ripperda et al. 2022) show clear indications of episodic magnetic reconnection events occurring close to the black hole event horizon. During such a flux eruption, a large part of the accretion flow can be ejected, and a low-density highly magnetized "magnetosphere" forms near the black hole. In GRMHD, this magnetospheric region upstream of the reconnection layer, i.e., at the jet base, consists of low-density plasma with a magnetization set by the density floor, indicating that pair plasma from the broadened base of the jet supplies the matter into the current sheet (Ripperda et al. 2020(Ripperda et al. , 2022. From first-principles kinetic simulations we know that even in the presence of strong radiation drag, magnetic reconnection renders itself as an effective process for particle acceleration Guo et al. 2014;Werner et al. 2016;Hakobyan et al. 2019a;Werner et al. 2019;Sironi & Beloborodov 2020). As we demonstrate below, reconnection near the event horizon of M87 * occurs in the radiative regime, when the synchrotron cooling timescale is comparable to the particle acceleration time. In this case, most of the dissipated magnetic field energy is radiated as synchrotron photons with a broad spectrum peaking around ε c = (3β rec /α F )m e c 2 ≈ 20 MeV, where α F is the finestructure constant, and β rec ≈ 0.1 is a dimensionless number characterizing the rate at which reconnection occurs in the collisionless regime. If the luminosity of the synchrotron emission is large, the plasma content of the current sheet will self-consistently be regulated by pair production fueled by these synchrotron photons spanning from keV to GeV energies. The important di-mensionless parameter that determines the efficiency of pair acceleration in reconnection is (twice) the available magnetic energy per particle, otherwise referred to as the magnetization parameter : σ = B 2 /(4πρ e ± c 2 ), where B is the characteristic strength of the unreconnected magnetic field (i.e., the field upstream of the current sheet, in the jet), and ρ e ± is the density of pair-plasma.
The characteristic optical depth for the photons interacting with the other synchrotron photons of average energy ε 1 and luminosity L is: τ γγ ≈ f γγ σ Tṅ (w/c)w ≈ 3f γγ σ T L/(4πcwε), where we assumed that (4πw 3 /3)ṅ ≈ L/ε, w ∼ 10 r g is the typical length scale of the reconnecting region, and f γγ ≈ 0.25 (set by the maximum γγ cross section; see, e.g., Akhiezer & Berestetskij 1985). To get a rough estimate on the fraction of the radiated photons being converted back into pairs, we evaluate the fiducial optical depth, assuming L ≈ f rec L jet (fraction of the jet power), and ε ≈ m e c 2 : τ 0 γγ ≈ 3f rec f γγ σ T L/(4πm e c 3 w) ∼ 10 −4 1. As it becomes clear in sec. 2.1, most of the photons that participate in the pair-production process also happen to be the ones that carry most of the dissipated power. This suggests that the fiducial value for the optical depth is a good proxy for the actual population-averaged optical depth, τ 0 γγ ≈ τ γγ . The amount of produced pairs (over the whole magnetosphere) can then be estimated asṄ e ± ≈ (1/2)σ T L 2 /(wcε 2 c ) (for a detailed derivation see the appendix A), where we implicitly assumed that most of the synchrotron energy is carried by the photons with energies around a few-to-ten MeV, i.e., ε c . 1 To obtain the multiplicity of the produced plasma with respect to the Goldreich-Julian number density, we can compareṄ e ± with the Goldreich-Julian number flux,Ṅ GJ ≈ (cL jet ) 1/2 /|e|: 2 N e ± N GJ ∼ 10 6 f rec 0.1 2 L jet 10 43 erg/s 1 The justification of this ansatz will become clear later in the section 2.1, where we show that the peak of the emission is slightly above MeV. Even for a relatively broad distribution this formula is a reasonable approximation for the total pair-production efficiency, as long as the peak is around mec 2 and the optical depth is small (see, e.g., Hakobyan et al. 2019b). 2 To estimate the Goldreich-Julian number flux, we takeṄ GJ ≈ I GJ /|e|, where I GJ is the magnetospheric current induced by the rotation of the black hole in the jet region. From Maxwell's equation for the toroidal component of the magnetic field, (4π/c)I GJ ≈ 2πwB φ , where B φ also enters into the expression for the Poynting flux: S P ≈ (c/4π)B 2 φ . Connecting the Poynting flux with the jet power, L jet ≈ πw 2 S P , we then arrive at where we assume that a fraction of the Poynting flux carried by jet, L rec ≈ f rec L jet , is dissipated during the reconnection event and is radiated away, and M BH and r g = GM BH /c 2 are the mass and the gravitational radius of the black hole, normalized to fiducial values for M87. This estimate is well justified in the regime where the effective optical depth to pair production is small, τ γγ ∼ 10 −4 , and both high and low energy photons participating in pair production are produced by the same dissipation mechanism. We demonstrate a more rigorous estimation for the f rec parameter in sec. 2.3.1. Based on this estimate, one can clearly see that the dynamics of the reconnecting current sheet is controlled by the density of the pair plasma n e ± ≈ 10 6 n GJ produced in-situ, which is much larger compared to the density produced by gaps or pair "drizzles" in the jet (Moscibrodzka et al. 2011;Crinquand et al. 2020;Chen & Yuan 2020;Wong et al. 2021;Yao et al. 2021). This yields a value for the plasma magnetization: where ω B = |e|B/m e c, and Ω BH ≈ (2r g /c) −1 . In the situation when synchrotron cooling is dynamically important, pairs are able to accelerate to energies E e ± ≈ σm e c 2 , forming a hard power-law energy distribution f e ± (γ) ∝ γ −1 -γ −1.5 , with γ being the Lorentz factor of the pairs (we demonstrate this directly from simulations in sec. 3; also see, e.g., Hakobyan et al. 2019a, and Chernoglazov et al. in prep.). A complete schematic illustration of the current sheet with all the sources for photons is shown in Figure 1. In the next sections we discuss the main radiation mechanisms and their energetics during the transient reconnection event.

The role of the radiation drag
To quantify the effects of cooling and its dynamical importance with respect to the acceleration by reconnection, it is helpful to define a dimensionless parameter, γ rad , as E rec |e| ≡ 2σ T U γ 2 rad , where the left-hand side is the accelerating force experienced by particles in the reconnection sites, while the right-hand side corresponds to a radiation drag force either due to synchrotron or IC (where we assume the Thomson regime and σ T is the Thomson cross section) radiation. Here, E rec is the electric field strength in the so-called x-point regions of the reconnecting current sheet, where the bulk of particle energization and energy dissipation takes place. Its value is controlled by the strength of the reconnecting magnetic field, and the rate at which the magnetic energy is dissipated (reconnection rate); it can be written as E rec ≈ β rec B. We take β rec = 0.1, which is characteristic for collisionless relativistic reconnection. The value of U in the definition of γ rad can either be interpreted as the energy density of the magnetic field, U = U B ≡ B 2 /8π, or as the energy density of the background (seed) photon field, depending on whether we are considering synchrotron or IC cooling. When the forces are balanced, particles with energies γ γ rad will cool down faster (either via synchrotron or IC) than they are able to accelerate, whereas cooling will be negligible on acceleration timescales for particles with γ γ rad .
Thus, comparing γ rad with the magnetization that powers the acceleration, σ, offers an important insight for understanding the cooling efficiency in the context of acceleration by reconnection. For synchrotron radiation, which for M87 is by far the strongest cooling mechanism, the parameter can be estimated as γ syn rad = (4πβ rec |e|/(σ T B)) 1/2 ∼ 4 · 10 6 (B/100 G) −1/2 . One can already draw two important conclusions from these estimates. First, since σ γ syn rad , the reconnection proceeds in the radiative regime where most of the dissipated magnetic energy is quickly converted into synchrotron radiation on timescales much shorter than the dynamical time of the system. This suggests that the energy density of the synchrotron photons in the steady state is comparable to the dissipated magnetic energy density, U syn rad ≈ β rec U B . Note that the synchrotron cooling is negligible in x-points, where the magnetic field vanishes, meaning that the fast acceleration process is uninhibited by synchrotron radiation; the highest energy particles will still accelerate to E e ± ≈ σm e c 2 . The peak of the synchrotron emission will then be set by the burnoff limit: ε c ≈ ω B (γ syn rad ) 2 ∼ 20 MeV (Uzdensky & Spitkovsky 2014). 3

Radiation from secondary pairs
The optical depth for synchrotron photons produced in reconnection is small τ γγ 1 for M87 (Ripperda et al. 2022), yielding a very non-local pair production in the whole region of size w in the black hole magnetosphere. These pairs, produced by photon-photon collisions, have arbitrary pitch angles to the local magnetic field and thus radiate rapidly giving rise to a distinct synchrotron spectral component different from the synchrotron continuum from reconnection. As the main fuel for these secondary pairs are the synchrotron photons from reconnection, their characteristic Lorentz factors at the time of birth can be estimated as γ sec ≈ (few · ε c )/(m e c 2 ) ∼ 200. The synchrotron cooling timescale for the secondary pairs, t sec , compared to the light crossing time, is t sec c/w ≈ (γ sec l B ) −1 ∼ 10 −3 , where l B = σ T U B w/(m e c 2 ) ∼ 1 is the magnetic compactness parameter (Beloborodov 2017), which means that a significant fraction of the energy contained in these pairs is radiated away inside the magnetosphere. The synchrotron peak energy associated with these pairs is ε sec ≈ ω B γ 2 sec ∼ 0.01-0.1 eV, with the energy den-   sity of resulting photons being controlled by the optical depth U sec rad ≈ U syn rad τ γγ ∼ 0.006 erg/cm 3 , and the luminosity being close to L sec rad ≈ L rec τ γγ ∼ 10 38 erg/s.

Generation of the ultra-high energy signal
Pairs accelerated by reconnection to energies E e ± 10 6 -10 7 m e c 2 can Compton-upscatter lower energy photons to TeV energies (the upper limit will be determined by the maximum pair energy). There will be several variations of this TeV signal. First of all, the reconnection region is filled with a soft background photon field from the accretion disk with average energies around ε soft ≈ 0.001 eV (300 GHz) and a characteristic energy density U soft rad ≈ 0.01 erg/cm 3 (Broderick & Tchekhovskoy 2015; EHT MWL Science Working Group et al. 2021). On the other hand, the region is also filled with synchrotron photons both from the pairs accelerated in reconnection and from the secondary pairs produced in the magnetosphere (as discussed in the previous section). Primary synchrotron photons from reconnection possess a wide range of energies up to a few ε c ≈ 10-100 MeV. The energy density of these photons is comparable to magnetic field energy density U syn rad ≈ β rec U B ∼ 40 erg/cm 3 , whereas the energy density of the secondary synchrotron photons is the aforementioned U sec rad ∼ 0.006 erg/cm 3 . A TeV signal can be produced by upscattering of either of the photon components populating the reconnection region by the highest energy pairs with γ 10 6 : the soft background (we will refer to this channel as IC), the primary synchrotron photons (synchrotron self-Compton, SSC), and the secondary synchrotron photons (secondary synchrotron self-Compton, SSC2). In all of these cases most of the power will be focused in the narrow frequency range corresponding to the highest energy pairs, i.e., above a few hundred GeV (Blumenthal & Gould 1970).
We will assume that the reconnection produces a power-law distribution of pairs with a cutoff, f e ± (γ) ∝ γ −p e −γ/γc , where the cutoff is γ c ≈ σ, and 1 < p < 1.5 (Chernoglazov et al. in prep.), an assumption that we later verify using kinetic simulations. Further in this section, we will be comparing the Compton luminosi-ties peaking in TeV (IC, SSC, and SSC2) with the synchrotron luminosity which peaks at few tens of MeV. The power density of this emission is equal to P syn ≈ 2σ T cn e ± γ 2B2 ⊥ /(8π) , where the averaging is done over the particle distribution x ≡ dγf e ± (γ)x, andB ⊥ is the effective magnetic field strength perpendicular to the motion of the particle (see eq. (4) for the full definition). As demonstrated in our simulations in sec. 3, we can accurately approximate γ 2B2 ⊥ /(8π) ≈ χ 2 (γ syn rad ) 2 U B (as long as the cooling is strong), with U B being the magnetic energy density in the unreconnected upstream, and χ ≈ 1 being a dimensionless parameter of order unity. Using that, one finds P syn ≈ 2σ T cn e ± χ 2 (γ syn rad ) 2 U B .

Total synchrotron luminosity
In the previous section, we assumed that the radiated synchrotron power during the reconnection is a non-negligible fraction of the jet power (L syn ≈ L rec ≈ f rec L jet , where f rec ∼ 0.1). To obtain this estimate more rigorously let us multiply the emitted synchrotron power density by the volume to get L syn = P syn απw 2 δ cs /(2π), where α ∼ 1 is the azimuthal extent of the region undergoing reconnection, w and δ cs are the radial extent and the thickness of the current sheet. As we demonstrate from the simulations in sec. 3, the characteristic temperature (or the mean energy) of particles in the current sheet is dictated by the synchrotron burnoff limit: 3T e ± = (κγ syn rad ) m e c 2 , where κ is a dimensionless parameter defined as follows: γ = κγ syn rad . Given that, we can eliminate the number density, n e ± , using the pressure balance condition: n e ± ≈ U B /T e ± . The width of the current sheet in collisionless reconnection is set by the characteristic Larmor radius of particles (e.g., Uzdensky & Spitkovsky 2014): δ cs ≈ γ m e c 2 /(|e|B). The total radiated synchrotron power then reads: where L BZ ≡ B 2 0 r 2 g c/24 is the Blandford & Znajek (1977) luminosity for a maximally spinning black hole, and we further assumed that the magnetic field falls with the distance from the horizon: B = B 0 (r/r g ) −1 . Notice that the dependence on w has disappeared, due to the fact that the magnetic field strength decreases linearly with distance, and so does the effective value of γ syn rad . Provided that the jet power is close to the BZ value L jet ≈ L BZ , we can state that radiatively efficient reconnection has a potential of emitting about 10% of the jet power, i.e., indeed, f rec ∼ 0.1.
As discussed earlier in this section, the bulk of the radiated power, L syn , will be carried by the photons with energies close to the burnoff limit, ε c ∼ 20 MeV. The very-high-energy (TeV) signal in our model is produced via Compton upscattering of lower-energy photons by the accelerated pairs, and will thus contain only a small fraction of the total synchrotron power. To estimate the power in the TeV component of the observed emission, in the following sections 2.3.2, 2.3.3, and 2.3.4 we investigate three different channels, which essentially consist of three populations of photons that can be upscattered to produce the TeV emission.

Upscattering of soft photons from the disk
The IC scattering of soft disk photons occurs in the Thomson regime, as the typical energies of photons in the lepton rest-frame, are small, i.e., Γ ε = 4γε/(m e c 2 ) 1, with ε being the energies of soft background photons in the lab frame before the scattering, ε/(m e c 2 ) ≈ 10 −9 , and γ ≈ σ ∼ 10 7 being the characteristic Lorentz factors of the leptons. The power per unit volume for the TeV signal can, thus, be found to be (see, e.g., Rybicki & Lightman 1979) P IC ≈ (4/3)σ T cn e ± γ 2 U soft rad . For γ c 1 and 1 < p < 1.5 we can approximate the average lepton energy to be the cutoff energy, i.e., γ 2 ≈ γ 2 c . We may then compare this power density to the total synchrotron power density to obtain P IC /P syn ∼ 10 −3 U (soft) 0.01 B −2 100 F −2 0.1 , where we assumed that most of the IC power is generated by particles close to the cutoff energy, γ c ∼ σ. 4

Upscattering of synchrotron photons from secondary pairs
Photons from the fast cooling of the secondary pairs have energies ε sec ≈ 0.01-0.1 eV. Collisions with the most energetic particles will occur in the regime intermediate between the Thompson and Klein-Nishina scattering approximations, where Γ ε ≈ 0.1-1. As an upper bound for the power we can use the same Thomson scattering relation as before, P SSC2 ≈ (4/3)σ T cn e ± γ 2 U sec rad . More realistically, one may use an approximation for the Γ ε ≈ 1 regime: P SSC2 ≈ (1/48)σ T cn e ± (m e c 2 ) 2 U sec rad /ε 2 sec ; see eq. (B3) in the appendix (B). From these two estimates we find that P SSC2 /P sync ∼ 10 −6 -10 −3 .

Upscattering of reconnection-produced synchrotron photons
The synchrotron photon field contains a much larger energy density than the soft background photons, However, most of this energy is contained in photons around a few to a few tens of MeV, ε c . In this case most of the SSC power is generated by the upscattering of low energy photons with energies ε s 1 eV by the highest energy particles, so that Γ ε ≈ 1. 5 In this regime neither the Thomson (Γ ε 1) nor the Klein-Nishina (Γ ε 1) approximation is valid. For a power-law distribution of particles one can estimate the power to be close to where ε s ≈ ε c is the average energy of synchrotron photons. The value of this power is very sensitive to the shape of the power-law, p, as the latter determines the relative amount of photons with energies ε s 1 eV. From this approximation we deduce that the SSC power compared to the synchrotron power for values of the power law 1 < p < 1.5 is P SSC /P sync ∼ 10 −9 -10 −7 , which is clearly subdominant compared to the IC signal.

Escape of the high-energy signal
Pair production may greatly inhibit the observed luminosity of both the synchrotron MeV and the inverse-Compton TeV signal if the optical depth for these photons is large. As we estimated above, the optical depth for MeV photons interacting with each other is small, τ γγ ∼ 10 −4 . Outgoing TeV photons, however, will interact with ∼ 1 eV photons from the soft background. We may estimate the optical depth to that process similarly: τ TeV ∼ 0.1σ T wU eV rad /(1 eV), where U eV rad is the energy density of background photons with energies around 1 eV. Given the total energy density of the soft radiation, U soft rad , as well as the slope of the distribution beyond ε soft ≈ 10 −3 eV, εdn/dε ∝ ε −1.2 (Broderick & Tchekhovskoy 2015), we can take U eV rad ≈ U soft rad (ε soft /1 eV) 0.2 . We then find an upper limit for the optical depth: τ TeV ≈ 1, meaning that a large fraction of the TeV signal escapes the inner few r g .
Our predictions are shown in Figure 2, where we plot the emerging spectra from all the emission components of the reconnecting sheet. In this plot we fix the main (synchrotron) luminosity to 10 43 erg/s, while other components scale appropriately. The synchrotron component from the secondary pairs (magenta line) is roughly estimated by taking a power-law distribution of pairs with a cutoff around γ sec (see sec. 2.2).
We vary a 5 Here we assume that the distribution of angles between the momenta of pairs and those of the photons is uniform. This assumption does not have a big impact on the amplitude of the outgoing very-high-energy flux, as we argue in the appendix D.
few parameters of the problem: the effective magnetic field strength, B, in the region where most of the synchrotron emission takes place, 6 the power-law index of the accelerated pairs, p, and their cutoff energy, γ c . For comparison, we overplot the spectrum of soft photons from the disk at energies 10 16 Hz from Broderick & Tchekhovskoy (2015).
The main takeaway of these estimates is that the outgoing TeV signal is comprised of the IC scattering of soft background photons from the disk, synchrotron photons produced in-situ by the reconnecting sheet, and from the cooling of secondary pairs outside the reconnecting sheet. Our estimates show that the Compton scattering of soft background photons (IC component) contributes most to the very high energy signal, with the total TeV luminosity up to ∼ 0.01-0.1% of the jet power, or about 10 39 -10 41 erg/s. Synchrotron radiation from secondary pairs appears as a distinct component at energies ∼ 0.01-0.1 eV, with a luminosity up to 10 38 -10 39 erg/s.

SIMULATIONS
Our reconnection-powered flaring model described in the previous section relied on multiple quantitative assumptions and relations, the validity of which can only be tested by first-principles simulations. In this section we describe the results from 2D particle-in-cell (PIC) simulations of relativistic radiative reconnection, performed to verify the validity of these assumptions. The simulations are performed using Tristan v2 -a multispecies radiative PIC code developed by Hakobyan & Spitkovsky (2020).

Setup
All of our simulations are initialized with a cold (T e m e c 2 ) uniform background pair-plasma of density n 0 , and a hot, thin, overdense layer of width ∆ along the xcoordinate in the middle of the box (at y = 0). The magnetic field is initialized as B = B 0x tanh (y/∆) (where B 0 is the magnetic field strength in the far upstream), and E = 0 everywhere. The width, the temperature, and the density of the layer, as well as the in-and outof-plane velocities of the particles in the layer are chosen in such a way that the setup is initially in the Harris equilibrium (Harris 1962). Near the x-boundaries of the box we impose absorbing boundary conditions that In our calculations, we fix the luminosity of the synchrotron emission, and vary the effective magnetic field, 10 < B [G] < 10 3 , the power-law index, 1 < p < 1.5, and the cutoff energy of the accelerated pairs, 10 7 < γc < 10 8 .
damp the fields and absorb outgoing particles . Along the x-axis at a fixed distance from the layer we inject leptons to replenish the upstream plasma. These boundary conditions allow us to run our simulations for many light-crossing times of the box, ensuring the results are insensitive to the initial conditions.
In all of our simulations, the value of the magnetization in the upstream is fixed: σ = B 2 0 / 4πn 0 m e c 2 = 10 3 . The size of the box is fixed at a value of 160 × 80 in units of σρ 0 , where ρ 0 = m e c 2 / (|e|B 0 ) is the Larmor radius of fiducial particles with a velocity γβ = 1. The value σρ 0 can thus be interpreted as the characteristic Larmor radius of particles with γβ ∼ σ in the field equal to the background value; the size of the box is chosen large enough to contain many gyro-orbits of the most energetic particles. Our simulations have a resolution of d 0 = c/ω p0 = 5∆x, with d 0 and ω p0 being the skin-depth and the plasma frequency of the background plasma; this means that the size of the box is 24000 × 13000 cells. The particle distribution function is sampled with 5 particles per cell (PPC) in the upstream, with significantly more in the reconnection layer; the results are virtually unchanged for higher PPC. Deposited currents are smoothed with 8 passes of a digital filter with a stencil (1/4, 1/2, 1/4) at each timestep.
To model synchrotron cooling, our simulations also include the radiative drag force in the Landau & Lifshitz (1975) formulation: assuming a particle moves with a four-velocity γβ in the background electromagnetic fields, E and B. To be able to properly model the effect of the drag force at realistic timescales, we upscale its magnitude introducing the dimensionless parameter, γ syn rad , discussed in sec. 2.1. The magnitude of the force is scaled in such a way, that |f rad | = 0.1B 0 |e| for γ = γ syn rad , |B| = B 0 , |E| = 0, and β ⊥ B (with B 0 being the strength of the upstream magnetic field). The ratio of γ syn rad /σ determines the regime of the synchrotron cooling, with γ syn rad /σ 1 corresponding to the dynamically strong cooling regime.

Results
To start the reconnection process we remove the pressure balance in the middle of the box, which triggers the formation of two magnetic islands propagating away from the center of the box. After about one light-crossing time, fast plasmoid-mediated magnetic reconnection is triggered across the whole current layer.
In Figure 3 we plot three different quantities for three simulations with varying cooling strength (γ syn rad /σ). Panels a1-a3 show the overall pair-plasma density in units normalized to the upstream density. In these panels, one can clearly see the main features of the reconnecting current sheet: the plasmoids -circular magnetic structures containing the accelerated plasma, and the x-points in between them, where the magnetic field vanishes and the main energy dissipation takes place. Plasmoids, which travel along the current sheet, adiabatically grow and merge with each other, are held intact via the balance between the magnetic and plasma pressure. From panels a1-a3 one can clearly see that when the synchrotron cooling "removes" some of that plasma pressure, densities (and sizes) of the plasmoids have to adjust to maintain the proper balance.
In panels b1-b3 we plot the mean energy of particles in each simulation cell, normalized to m e c 2 . From these three panels, it is evident that the dimensionless parameter γ syn rad controls the relativistic temperature inside the plasmoids. As the cooling strength increases, i.e., the ratio γ syn rad /σ drops, the temperature inside the plasmoids decreases, while the density increases, retaining the pressure (the product of the two) roughly constant. We see no evidence of the variation in the plasmoid magnetic field strength for different cooling regimes.
In panels c1-c3 of Figure 3 we show the total synchrotron emissivity (synchrotron power radiated from a unit volume) for the highest-energy photons. This quantity is defined as dγf ± e γ 2B2 ⊥ , where f ± e is the distribution function of high-energy particles (with γ > σ/4), andB ⊥ is the effective magnetic field component perpendicular to the direction of the particle velocity, defined in eq. (4). When there is no synchrotron cooling, or when the cooling is dynamically weak (i.e., the regime least applicable to the reconnection layer in the magnetosphere of M87*), plasmoids carry most of the highenergy particles. These particles lose their energies at timescales much longer than the characteristic injection timescale from the current sheet, and the plasmoids get rapidly replenished with fresh accelerated plasma. In this regime most of the high-energy photons are thus produced inside the largest plasmoids (c1-c2). In the strong cooling case, on the other hand, the highestenergy particles are rapidly cooled after they leave the acceleration sites (x-points). Because of that in panel c3 we can only observe the smallest (youngest) plasmoids, as well as the relativistic outflows in the sheet carrying the freshly accelerated pairs, as sources of energetic synchrotron photons.
In the radiatively efficient regime, γ syn rad /σ 1, a fraction which we will denote χ 2 , of the dissipated magnetic energy inside the current layer, p diss ≈ |e|E rec c, is radiated via the synchrotron mechanism: p syn ≈ 2σ T c γ 2B2 ⊥ /(8π) = χ 2 p diss (both powers are computed per particle). Here E rec ≈ β rec B 0 is the accelerating electric field in the layer. Using the definition of γ syn rad one can find: γ 2B2 ⊥ ≈ χ 2 (γ syn rad ) 2 B 2 0 , a formula that was used in sec. 2.3.1 when computing the total radiated synchrotron power. To find the value of χ we compare four different simulations with varying cooling strength: one with weak cooling (γ syn rad /σ = 2.5), one with moderate cooling (γ syn rad /σ = 1), and the other two with strong cooling (γ syn rad /σ = 0.2, 0.5). The values for the dimensionless parameter χ measured from simulations are shown in Table 1; the results indicate that for a wide range of values of γ syn rad our approximate analytic argument of χ being constant holds rather well (we inspect this more closely in the appendix C). 7 We also cite the values of κ ≡ γ /γ syn rad , which was used to estimate the effective temperature of the current sheet. And again, varying the cooling strength within an order of magnitude only marginally affects this parameter, which is close to 0.1-0.3.  Table 1. Characteristic values of physical parameters that have observational significance measured directly from our simulations with different cooling strengths. The first column γ syn rad /σ indicates the cooling strength, with the infinity corresponding to the run without radiation; βrec is the time-averaged reconnection rate (normalized to c); the ratio χ ≡ γ 2B2 ⊥ 1/2 / (γ syn rad B0) provides an estimate for the total radiated synchrotron power; κ ≡ γ /γ syn rad correlates the effective temperature with the cooling burnoff limit; p and γc are the best fit power-law index and the cutoff energy for particle distribution functions f ± e ∝ γ −p e −γ/γc . In panels where we cite two values of p we fit a broken power-law with an ankle near γ ≈ γ syn rad .
We also directly measure the reconnection rate by evaluating the component of the E × B/|B| 2 velocity in the upstream, pointing towards the current sheet (in  . Snapshots taken at the same exact timestep from three different 2D simulations of the magnetic reconnection with varying synchrotron cooling strength. Panels from top to bottom show: (a1-a3) the pair-plasma density overplotted with magnetic field lines, (b1-b3) the mean Lorentz factors of particles, (c1-c3) synchrotron emissivity, dγf e ± γ 2B2 ⊥ , for particles with γ > σ/4. In all three simulations, the upstream plasma magnetization is σ = 10 3 . The ratio γ syn rad /σ quantifies the strength of synchrotron cooling with smaller values corresponding to stronger cooling. As the cooling strength increases, plasma inside the plasmoids efficiently radiates away the excess momentum perpendicular to the magnetic field. As a result, plasmoids in the strongly cooled case become more compressed (see a1-a3), while the mean energy of the particles inside plasmoids corresponds to the radiation limit, γ ≈ γ syn rad /4. When the cooling is weak, high-energy particles inside plasmoids contribute the most to the synchrotron emission (c1-c2), while in the strong cooling case only the smallest plasmoids and the relativistic outflows in the current sheet are visible. The x-coordinate in these snapshots is measured in units of the Larmor radius, ρL of particles with energies ∼ σmec 2 , i.e., ρL ≈ σρ0, where ρ0 = mec 2 /(|e|B0), with B0 being the magnetic field strength in the far upstream. Figure 3 that corresponds to the ±y direction). We find that the strength of the synchrotron cooling is unimportant in determining the reconnection rate, with only marginal variations likely caused by the intermittency. The value of the rate itself is close to β rec ∼ 0.25-0.28, in agreement with what has been found by Guo et al. (2015), and Sironi et al. (2016) in 2D simulations.  Runs with varying cooling strength, γ syn rad /σ, are shown with different colors, while the magnetization, σ, is fixed for all the runs. Panel (b): time-averaged synchrotron and IC spectra based on the simulation particles. The synchrotron spectrum (peaked near ωB (γ syn rad ) 2 ) is computed on-the-fly, and includes the self-consistent values forB ⊥ . The IC spectrum (peaked near TeV) is computed in post-processing, assuming an isotropic soft photon background, with the particle Lorentz factors rescaled to realistic values. Dashed lines on both panels show the fiducial case of f e ± ∝ γ −1 e −γ/γc (the cutoff is fixed at γc = σ), which is expected to be applicable for σ ≈ 10 7 in the M87 case, and the analytic prediction for the corresponding synchrotron and IC spectrum.
Particle-in-cell simulations also enable us to directly measure both the energy distribution of pairs and the spectra of both synchrotron and IC photons from first principles (shown in Figure 4). Particle distribution functions for different values of the cooling strength, γ syn rad /σ, are shown in panel (a) with different colors. Deduced power-law indices and cutoff energies for the best fit power-law (or double-power-law) of the form f e ± ∝ γ −p e −γ/γc (or ∝ γ −p1 for γ < γ syn rad , and ∝ γ −p2 e −γ/γc for γ γ syn rad ) are shown in the Table 1. In the weak-to-marginal cooling regime, when γ syn rad /σ 1, the distribution of pairs is best described by a single power-law, with an index varying between 1 and 1.5 (which justifies the values used in our analytical model in the previous section). When the cooling becomes stronger, we see the second power law generated beyond γ γ syn rad with a slightly steeper index of ≈ 2.5. In all of the cases, however, except for the uncooled regime, the energy cutoff is strictly controlled by the value of σ, since the highest-energy pairs are produced in the x-points where there is no synchrotron cooling. In the uncooled case, which is unrealistic for our astrophysical scenario, the spectrum extends much further than σ due to various secondary acceleration channels (Petropoulou & Sironi 2018;Hakobyan et al. 2021;Zhang et al. 2021).
In Figure 4b we show both the synchrotron spectra computed on-the-fly, as well as the IC spectra with an assumed soft isotropic background radiation computed in post-processing. The peaks of the synchrotron emission in all of the cases correspond to ω B (γ syn rad ) 2 , which agrees with the earlier predictions by Uzdensky & Spitkovsky (2014). The peaks for the IC signal are selfconsistently dictated by the energies of leptons. There are several important differences to the analytic expectation (shown with a dashed line). First, at lower energies our simulation underpredicts the synchrotron flux. This is likely caused by a combination of two reasons: (a) limited scale separation leading to smaller number of low-energy particles in the sheet (we only consider radiation from particles in the current layer), and (b) in 2D, plasmoids, where most of the radiation takes place, are compressed, and magnetic field strengths are stronger than the upstream, which can yield an on average marginally higher radiation frequency. Also notice, that at higher frequencies, ε ω B (γ syn rad ) 2 , in the synchrotron spectrum we see slightly elevated luminosity, especially for the strongest cooling cases. This is due to the fact that in the strong cooling regime, particles beyond γ γ syn rad cool very rapidly. This makes the spectrum more intermittent, and, when averaged over time, results in more flux at high energies (cf. Hakobyan et al. 2019a). To see this, let us estimate the characteristic maximum photon energy, or the cutoff energy. We can use the peak synchrotron energy formula: E cut ≈ (3/2) γ 2 max (|e|B ⊥ )/(m e c 2 ), where γ max ≈ σ, andB ⊥ is the perpendicular component of the magnetic field for the highest-energy particles. Since σ γ syn rad (strong cooling regime), particles with Lorentz factors close to γ max rapidly lose their energy via synchrotron radiation, emitting the perpendicular component of their momenta. Because of this, the effectivẽ B ⊥ is smaller than the value of the upstream field, and may be approximated as γ maxB⊥ ≈ γ syn rad B 0 (see also appendix C). We may then simplify the expression for the spectral cutoff: E cut ≈ (3/2) ω B (γ syn rad ) 2 (σ/γ syn rad ) ≈ 16 MeV (σ/γ syn rad ). We can then conclude that the spectral cutoff for the synchrotron emission is at higher energies for the stronger cooling regimes.
In our simulations, we had to downscale some of the physical parameters to be able to tackle the problem computationally. In particular, the expected magnetization parameter of the upstream plasma, σ, during a flare from the magnetosphere of M87* is expected to be close to 10 7 (see eq. (2)), whereas in our simulations we employed σ = 10 3 . At the same time the synchrotron limit, γ syn rad , was downscaled from about 10 6 to 200. In light of this, there might be a concern that some of the arguments we have made in this section might not extrapolate when the problem is rescaled to realistic parameters. First of all, it is important to notice that we still retain the strong-cooling hierarchy, γ syn rad < σ, which conserves the balance between the competing mechanisms (acceleration and cooling) when the dimensionless values are downscaled. On top of that, employing σ ∼ 10 3 allows us to have at least a decade of scale separation between σ, γ syn rad , and u A /c ∼ √ σ (Alfvénic four-velocity). We also do not expect the overall dynamics of the layer to change drastically when u A /c 1, as for these values the Alfvénic three-velocity, which dictates the characteristic velocities of the flow in the sheet, approaches the speed of light.
Our simulations are two-dimensional, while the real current sheet in the vicinity of the black hole is inherently 3D. In particular, 3D reconnection has been demonstrated (see Zhang et al. 2021) to enable a fast acceleration channel for particles beyond γ ∼ σ, which in 2D operates at a much slower rate (Petropoulou & Sironi 2018;Hakobyan et al. 2021). Ironically, the peculiarity of radiative reconnection helps us in this scenario, since in this regime any secondary acceleration channel outside the x-points is essentially forbidden due to the strong cooling, that vanishes only in these microscopic regions comparable in size to the plasma skin-depth.

DISCUSSION
The black hole in the center of the M87 galaxy is a perfect target to study not only the dynamics of the large-scale accretion flows but also the extreme plasma-physical processes that intermittently occur in the near vicinity of its event horizon. In particular, long timescales associated with the variability of its core as well as relatively high luminosity makes it a prime target for the observations of TeV flares. At the same time, independent measurements of the magnetic field strength and constraints on the jet power make analytic estimations well-constrained and a lot less ambiguous. While for these reasons the focus of the current paper has been on the M87* black hole, it is important to note that the results can be extrapolated to any system which transiently hosts a magnetospheric, radiatively efficient, reconnecting current sheet (such as potentially Sgr A* or Cen A, and even young energetic pulsars, such as the Crab pulsar; see, e.g., Lyubarskii 1996).

Summary
In this paper, we have shown that large-scale magnetic reconnection, occurring in the episodically forming current sheets in magnetically arrested accretion flows around black holes, can power bright multiwavelength flares extending to TeV energies. We have shown that the reconnecting sheet can dissipate a significant fraction of the jet power on timescales controlled by the universal reconnection rate: β rec ≈ 0.1-0.3; see eq. (3). While the efficient dissipation of the Poynting flux during magnetic reconnection is not a novel concept, 8 plasma and magnetic field parameters that control the dissipation determine the observational footprint of the process in a wide range of wavelengths.
Our findings indicate that synchrotron radiation is the main cooling mechanism for the leptons, through which the dissipated magnetic energy is converted into radiation. The synchrotron photons from the accelerated plasma, which peak at around a few tens of MeV, have finite optical depth and are thus susceptible to abundant pair production. This process in turn loads the current sheet with electron-positron pairs; the expected resulting pair-multiplicity with respect to the Goldreich-Julian density is given by eq. (1), and its estimated value is rather large: ∼ 10 7 .
We have shown that the expected magnetization parameter in the upstream of the reconnecting current sheet, which determines the available magnetic energy per lepton, is large: σ ≈ 10 7 , see eq. (2). Despite strong synchrotron cooling, reconnection can acceler-ate pairs to energies close to few-to-tens of TeV,9 or E e ± ≈ σm e c 2 . The highest-energy pairs can then Compton-scatter both the soft photons from the disk and the synchrotron photons from the sheet, producing the very-high-energy TeV component of the emission. Our estimations suggest that the upscattering of the soft thermal photons produced in the bulk of the accretion flow is the most efficient way of producing a TeV signal, with a predicted peak luminosity of 0.01-0.1% of the jet power, or about 10 39 -10 41 erg/s.
We have also verified some of our assumptions using first-principles 2D particle-in-cell (PIC) simulations of reconnecting current sheets with synchrotron cooling included self-consistently. In particular, our simulations support the claim that the high-energy cutoff of the distribution of particles is not sensitive to the synchrotron cooling strength (see table 1), since the particle acceleration sites (x-points) have zero magnetic field strength, and are thus not susceptible to synchrotron cooling. On top of that, both the power-law index (which we found to be close to p ≈ 1.5) and the reconnection rate (β rec ≈ 0.28) are also insensitive to the synchrotron cooling, as long as the cooling is strong.

Observational implications
Using particle distributions from our simulations we have generated both synchrotron and IC spectra (shown in fig. 4b), with particle Lorentz factors rescaled to realistic values. For the weak cooling regime our generated spectra overall match analytic predictions; both the spectral peak and the cutoff of the synchrotron component are determined by the value of γ syn rad , and are thus close to E peak ≈ ω B (γ syn rad ) 2 ≈ 20 MeV. However, in the strong cooling regime while the peak energy is still set by the value of γ syn rad , the spectral cutoff is higher E cut ≈ E peak (σ/γ syn rad ). This intermittent synchrotron excess at higher energies leads to a much less eminent gap between the synchrotron component (between fewto-ten MeV and GeV) and the inverse-Compton component (above 100 GeV), with the overall spectrum being more continuous.
GRMHD simulations indicate that during the formation of the current sheet the accretion flow in the 9 If there is an additional (unknown) efficient mechanism for plasma loading, the actual magnetization in the current sheet near the M87 black hole may be smaller than our calculations suggest, i.e., σ γ syn rad . In this case, particles can additionally accelerate via 3D secondary mechanisms after they leave the xpoints (see, e.g., Zhang et al. 2021;Chernoglazov et al. in prep.). Their energies will eventually be limited by the burnoff limit, γmec 2 γ syn rad mec 2 , i.e., above TeV, which still allows production of TeV emission by the IC scattering. magnetospheric region where reconnection takes place is repelled by the magnetic stresses beyond roughly w ∼ 10r g . As a consequence, the radio to near-IR flux from the disk is expected to dim significantly during the flare (Jia & et al. in prep.). In sec. 2.2 we discussed the possible emergence of a millimeter-to-infrared synchrotron signal from the secondary pairs produced in the reconnection upstream (shown with the pink in Figure 2). The luminosity of this signal is unlikely to exceed 0.1%-0.01% of the primary synchrotron counterpart (or about 10 38 -10 39 erg/s), as the pair production optical depth, τ γγ ∼ 10 −4 , is rather small. Distinguishing it from the dimmed background emission of the accretion flow will thus be challenging, however, a possible feature one should be looking for is enhanced values of the polarization degree in near-IR during the flare. The synchrotron radiation from secondary pairs is produced in the reconnection upstream (i.e., at the jet base), where the magnetic field is coherent, as the turbulent accretion flow would have been repelled. 10 To summarize, we expect that during the TeV flare, the emission of the disk below eV should dim, as the disk is repelled, and the reconnection region is evacuated, while the polarization degree in the same frequency range can increase, owing to the contribution from the synchrotron emission of secondary pairs.
Let us also briefly discuss the temporal characteristics of the flare at different wavelengths. The dimensionless reconnection rate, β rec , controls both the timescale of the eruption during which the magnetic energy dissipation takes place, as well as its peak synchrotron energy; ε c ≈ (3β rec /α F )m e c 2 , see eq. (3). In the relativistic collisionless regime, studied using PIC simulations for almost a decade, the value of the reconnection rate varies between 0.1-0.3 depending on the magnetization parameter, the dimensionality of the simulation, and the boundary conditions (our 2D simulations find values close to 0.25-0.3, as shown in the Table 1; 3D simulations tend to have smaller values, Guo et al. 2015). This prediction is an order of magnitude faster than what has been found by relativistic MHD simulations (Ripperda et al. 2019). For the largest current sheets, 3D GRMHD simulations show a flux decay period, i.e., a reconnection event that can power a flare, of ∼ 100r g /c (Ripperda et al. 2022), which is approximately 30 days for M87*. Due to an order of magnitude faster reconnection rate we expect the flux decay period, and hence the characteristic flare time scale in collisionless plasma to be shorter (Bransgrove et al. 2021), which aligns well with observational flare durations of the order of few days for M87*. 11 The magnetic compactness of the inner region near the horizon is significant: l B = σ T U B w/(m e c 2 ) ≈ few (Beloborodov 2017), with w ≈ 10r g being the length of the current sheet. This corresponds to the cooling timescale of t cool ≈ (w/c)/(γl B ) for particles with a Lorentz factor γ. Since the cooling time of even the marginally relativistic particles is short, the variability on timescales comparable to ∼ r g /c can only be expected from global effects, such as the beaming due to the orientation of the current sheet. At lower energies, the flare can also have shorter intermittency caused by the dynamics of relativistic outflows, the growth and collisions of plasmoids, and particle acceleration and cooling. For example, in the radio range, which is expected to be produced by the particles with Lorentz factors of the order of γ ∼ 100, the variability timescale can be comparable to 0.1r g /c (few hours for M87*). At higher energies the same estimate yields t int ≈ t cool ( ω B γ 2 = ε s ) ≈ 10 min · (ε s /eV) −1/2 . While this short-timescale intermittency will be virtually undetectable at the peak energies of MeV (around a second), its detection is possible in the optical to X-ray band, where we expect the characteristic cooling timescale to be around a few tens of seconds to minutes. the process more carefully. We assume that the reconnection produces a distribution of synchrotron photons N γñγ (ε) ∝ ε −(p+1)/2 exp −(2/3)(ε/ε max ) 1/3 (Zirakashvili & Aharonian 2007; Lefa et al. 2012), where N γ is the normalization of the distribution function, i.e., the total number of photons, whereas ñ γ (ε)dε = 1 (the distribution of pairs producing this synchrotron spectrum is assumed to beñ γ ∝ γ −p exp [−γ/γ max ], where ε max = (3/2) ω B γ 2 max ). We can then estimate the total pair production yield with the following relation (Svens-son 1987):Ṅ where V is the volume of the pair-production region in the reconnection upstream. Here a γγ (ε) = η p (m e c 2 /ε)ñ γ (m e c 2 /ε) is a dimensionless function describing the angle-averaged cross section for the twophoton pair production, and η p is a dimensionless constant that depends on the slope of the photon distribution. Performing the integration for ε > m e c 2 yields:Ṅ e ± ≈ (σ T c/V )N 2 γ η p Λ(ε max ), where Λ(ε max ) ≈ log (0.6 ε max /m e c 2 ) ∼ 3 (for ε max ≈ 20 MeV ∼ 40 m e c 2 ). For a power-law index p = 1 we can take, η p = 7/12 (Svensson 1987) to then conclude:Ṅ e ± ≈ 2(σ T c/V )N 2 γ . We can further take N γ ≈ V U γ / ε γ , with U γ being the average energy density of the synchrotron photons, or expressing the same in terms of the outgoing synchrotron luminosity: N γ ≈ L γ (w/c)/ ε γ (w is the characteristic size of the pair-production region), where we can also take ε γ ≈ ε max .

B. SSC LUMINOSITY
To obtain the synchrotron self-Compton power for the parameters applicable for M87* we cannot use the classical Thomson or Klein-Nishina approximations, as the characteristic energies of photons that contribute the most to the outgoing Compton-luminosity in the restframe of the highest-energy pairs is close to m e c 2 . Instead, we may derive an estimate for the luminosity in this regime using the following formula for the total power per lepton radiated by upscattering photons of energies ε 0 by leptons having Lorentz factors γ (see, e.g., Blumenthal & Gould 1970, Lefa et al. 2012: Here, the distribution functions of photons, dñ 0 /dε 0 , and pairs, dñ e ± /dγ, are normalized to 1, and n 0 is the number density of incoming synchrotron photons in the steady state. In this equation Γ ε ≡ 4ε 0 γ/(m e c 2 ), and G 0 (Γ ε ) is a slowly varying function of Γ with the asymptotes G 0 (Γ ε ) ≈ 1/9 for Γ ε 1, and G 0 (Γ ε ) ∝ (1/2) ln Γ ε −(11/6). Integration of the kernel in eq. (B2) over all values of ε 0 and γ yields the total SSC power radiated per unit volume.
If the synchrotron peak energy, ε s , is close to Γ ε = 4ε s γ c /(m e c 2 ) ≈ 1, we may safely use the Γ ε ≈ 1 (and G 0 ≈ 1/9) approximation in eq. (B2). In this case the total power per lepton is equal to On the other hand, for extended power-law distributions of pairs, and the synchrotron peak deep in the Klein-Nishina regime, we may no longer directly substitute Γ ε = 1, as the main contribution to the total flux comes from the photons with much lower energies than the peak. To simplify integrating eq. (B2), it is useful to substitute ε 0 = Γ ε /(4γ), and assume a power-law distribution for pairs, dñ e ± /dγ ∝ γ −p , and synchrotron photons dñ 0 /dε 0 ∝ ε −(p+1)/2 0 . Integration over ε 0 may then be substituted by the integration over Γ ε , and the total SSC power then reads Integration over Γ ε can be well approximated by a delta function near Γ ε ≈ 1 (for values of 1 ≤ p 3), to arrive at the following estimate for the radiated power per lepton: In sect 2.3.1 we approximated the distributionaveraged synchrotron power γ 2B2 ⊥ as χ 2 γ 2 rad B 2 0 , where χ is a dimensionless parameter of order unity, and B 0 is the strength of the upstream magnetic field. Indeed, as shown in table 1, the χ parameter varies only marginally for a range of different cooling strengths (γ rad /σ).
In figure 5 we inspect this more closely by plotting 2D histograms for all the particles in our simulations, binning them according to their values of γ and B ⊥ . It is evident, that both the distribution of particles, as well as the compression of plasmoids (i.e., the strength of the magnetic field and, most notably, the density; also evident from figure 3a) adjust in such a way that the majority of particles are piled near the line corresponding to the cooling-acceleration balance, i.e., σ T γ 2B2 ⊥ /(4π) ≈ |e|E rec . In the non-cooling case, unsurprisingly, the distribution is a lot more uniform with particles accelerating freely regardless ofB ⊥ .
In 2D as soon as the energized particles enter the plasmoids, they remain confined inside them until the plasmoids escape the simulation box. In more realistic 3D simulations particles may leave the plasmoids, while traveling in the longitudinal direction (along the plasmoid axis), and the plasmoids themselves can get disrupted due to the kink instability . Thus, in 3D we do not expect the magnetic field strengths inside the plasmoids to reach values above B/B 0 ∼ few. The radiation/injection balance discussed above, γB ⊥ ≈ γ syn rad B 0 , has to still be maintained, and thus the general implications of our 2D simulations will hold even in 3D. However, the details on how particle distribution, as well as the structures of plasmoids, adjust to the cooling-imposed balance condition remains to be seen in future 3D simulations.

D. BEAMED SYNCHROTRON SELF-COMPTON
In the estimate of the luminosity of the synchrotron self-Compton component, we assumed a random distribution of pitch angles, θ, between the photons and the scattering leptons. Because of this assumption the effective scattering cross section was suppressed by a Klein-Nishina factor γε (for photons with energies ≈ ε; in this section all the energies are measured in units of m e c 2 ) with respect to the Thomson cross-section. Here, we discuss implications of this assumption.
First, let us briefly reformulate our previous result on the energy yield of the SSC emission in terms of the effective scattering cross-section. The number of scatterings is proportional to the effective IC cross section, σ eff , while the energy yield scales with the maximum energy of scattered photons: P SSC ∝ σ eff ε max ≈ σ eff γ c . Comparing this with the synchrotron luminosity, P syn ∝ σ T ε c , yields: P SSC /P syn ≈ σ eff γ/(σ T ε c ) 1, where the exact value of σ eff depends on the energy distribution of incoming photons. From eq. (B5) we can see that for the applicable regime it is close to σ eff ≈ σ T γ −(p+3)/2 c (this is due to the fact that most of the outgoing energy is ac-cumulated by upscattering the low-energy synchrotron photons far below the peak, ε c ). This finally yields: P SSC /P syn ≈ 1/(γ (p+1)/2 c ε c ) 1. In the more general case the effective cross-section will be determined by the distribution of pitch angles of the incoming synchrotron photons dn εe ± /dθ: σ eff = σ ε /(γε) ≈ π 0 [σ θ (1 − β cos θ)(dn εe ± /dθ)]dθ (primed quantities are measured in the lepton reference frame: ε = γε(1 − β cos θ); see, e.g., Landau & Lifshitz 1975). Here σ θ ≈ σ T when θ θ T ≡ 1/ √ γε (Thomson regime), and ≈ σ T /(γε) -otherwise (Klein-Nishina regime). If the pitch-angle distribution is uniform, dn εe ± /dθ = 1/π, the effective cross section is suppressed by a factor of 1/(γε c ) (in agreement with our original conclusion for γ ≈ γ c and p = 1). One might expect that since the emission of the synchrotron photons is beamed in the direction of particles, Compton scattering of the same population will be more efficient, as the photons will have smaller energy in the lepton rest frame. In the most optimistic scenario when all the photons are beamed in the direction of particles, i.e., dn εe ± /dθ is a delta function in θ, we find σ eff = σ T (1 − β) ≈ σ T /(2γ 2 ). In terms of the outgoing luminosity this results in P SSC /P syn ≈ 1/(γ c ε c ). Since the optical depth for the synchrotron photons to Compton scattering in the large-scale reconnecting sheet of the black-hole magnetosphere is finite, the scatterings are not entirely local, and only a fraction of the scatterings can be considered beamed. The actual effective cross section is then σ T /(2γ 2 ) σ eff σ T /(γε), which is still σ T , and thus P SSC P syn .