yet been accessible to TeV instruments. From this point of view, the results of our model remain, for the moment, predictions. The coming online of the Cherenkov Telescope Array (CTA) over the next few years will provide increased sensitivity in the TeV band (Cherenkov Telescope Array Consortium et al. 2019), enabling enhanced temporal resolution and, hence, a more thorough probing of the IC(HDR) scenario.

We note that BL Lacs are typically much brighter than FSRQs at TeV energies, and thus many more IACT observations of BL Lacs exist – and typically with much greater temporal resolution – than of FSRQs (Wakely & Horan 2008). However, without strong external photon fields supplied by the BLR or the HDR, the source of seed photons for reconnection in BL Lacs is far less constrained. If these are the synchrotron photons produced during reconnection (a single-zone synchrotron self-Compton, or SSC, setup, cf. Maraschi, Ghisellini & Celotti 1992; Bloom & Marscher 1996), a much more detailed seed photon model – taking into account multichromaticity, time-dependence, anisotropy, and spatial inhomogeneity – needs to be folded in to the simulations in order to make robust predictions. However, if the seed photons impinge on the putative reconnection region from elsewhere in the jet (a multizone SSC paradigm), then our static seed photon model is more appropriate, but it would still potentially need to be generalized to the non-monochromatic case. We note that in the famous minute-scale flare of PKS 2155–304 presented by Aharonian et al. (2007), no evidence of spectral variability was found. This is consistent with our reconnection model, however the photon indices, even after correcting for the background absorption (e.g. Aharonian et al. 2009) are steeper than one would expect on the decaying side of the flare, wherein Klein–Nishina and pair-production physics should harden the spectrum closer to the universal shape of Section 4.5. Ultimately, not enough is known about the seed photons to draw firm conclusions, however.

We close this discussion by examining the potential pair yield of Klein–Nishina reconnection in the IC(BLR) and IC(HDR) scenarios. Suppose that the jet is launched with an electron–proton composition (no positrons). Let us also assume that the jet evolves so as to be moderately magnetized in its rest frame, $B_0^2/4\pi n_0 m_p c^2 \sim 1-10$, at the parsec scale (cf. Giannios 2013; Giannios & Uzdensky 2019; Mehlhaff et al. 2021), which is near the transition point where the HDR overtakes the BLR as the dominant seed photon source (Nalewajko et al. 2012; Mehlhaff et al. 2021). The cold electron magnetization is then $\sigma_{c,0} \sim (m_p/m_e)(B_0^2/4\pi n_0 m_p c^2) \sim 2 \times 10^{3-4}$. In our earlier work (Mehlhaff et al. 2021), we estimate (in the jet’s rest frame) $\gamma_{KN, BLR} \sim 300$ and $\gamma_{KN, HDR} \sim 1 \times 10^4$. This means that the pair yield control parameter (Section 5) is $\gamma_{pp, max}/4\sigma_{c,0} \sim 30\gamma_{KN}/4\sigma_{c,0} \sim 0.1-1$ for reconnection illuminated by the BLR and 4–40 for HDR irradiation. This control parameter only becomes small (the regime of high pair yield) in the IC(BLR) scenario, which is, incidentally, also the scenario for which available flaring observations, as discussed in this section, best corroborate a Klein–Nishina reconnection model. If pair production is predominantly decided by leptonic physics as unveiled by our pair-plasma simulations, then the present estimates suggest that reconnection functions as an *in situ* source of antimatter in FSRQ jets strongly illuminated by the BLR, possibly creating more than 1 positron per electron and effectively transforming any initially electron–proton plasma into a plasma with a prominent positron component. This is an important result in blazar studies because the jet composition is notoriously difficult to ascertain observationally (Madejski & Sikora 2016); if it is true, it could mean that BLR-illuminated FSRQ jets generically carry a strong pair-plasma component downstream of the GeV emission zone.

### 6.1.3 Synopsis of FSRQ discussion

In summary, the bright Doppler-boosted emission from blazars highlights the physics of particle acceleration in their jets. In particular, for the FSRQ sources, GeV and TeV observations function as respective probes of scenarios where magnetic reconnection is coupled, through Klein–Nishina and pair-production physics, to soft seed photons produced by the BLR and the HDR. Observations of FSRQ flares in the GeV band are in broad agreement with the spectral-temporal signatures of Klein–Nishina reconnection expected from this study (Section 4.5), constituting potential evidence that this type of reconnection occurs within the BLR of FSRQs. The expected pair yield (Section 5) of IC(BLR)-coupled reconnection could also contribute significantly to the antimatter content of FSRQ jets far away from their central engines. The IC(HDR) scenario cannot be probed at the same level of detail as the IC(BLR) case with the current generation of IACTs, a limitation which the CTA will help to overcome. Finally, while Klein–Nishina reconnection could also operate in BL Lacs, there are a lot more uncertainties concerning the seed photons in those sources, pushing detailed observational comparisons into the domain of future work.

## 6.2 Black hole accretion disc coronae

### 6.2.1 Background information

In our previous work, Mehlhaff et al. (2021), we made a case for Klein–Nishina reconnection operating in the coronae of accreting BHXRBs in their high/soft states. Conducting basic estimates, we showed that, for the case where an underlying optically thick, geometrically thin accretion disc (Shakura & Sunyaev 1973) supplies

$$ \epsilon_{\text{disc}} \sim 1 \text{ keV} \qquad (58) $$

seed photons to reconnection in a highly magnetized coronal plasma, the radiative scale hierarchy is $\gamma_{\text{cool}} < \gamma_{\text{KN}} < \gamma_{\text{rad, T}} < \gamma_{\text{max}}$ [equation (50)], pushing reconnection into the Klein–Nishina regime. Due to the intense accretion disc radiation bath, absorption fiducially kicks in for gamma-rays above the threshold energy,

$$ \epsilon_{\text{th,disc}} = m_e^2 c^4 / \epsilon_{\text{disc}} = 260 \text{ MeV}. \qquad (59) $$

Meanwhile, the typical photon energy emitted by a $\gamma_{\text{KN}}$ particle is $\gamma_{\text{KN,disc}}m_e c^2/2 = \epsilon_{\text{th,disc}}/8 = 30 \text{ MeV}$, where $\gamma_{\text{KN,disc}} \sim 100$.

To our knowledge, such energies have only been detected in the high/soft state of an accreting BHXRB – for which their origin in an ejected jet is not expected – in one object: Cyg X-1. This was during an approximately 100-Ms exposure by the Fermi-LAT, presented by Zdziarski et al. (2017), who report the detection of gamma-rays up to a cut-off energy of about 20–40 MeV, somewhat below (59). As pointed out in Mehlhaff et al. (2021), this cut-off may be consistent with gamma-ray absorption because the coronal region is likely highly radiatively compact (which translates into a high fiducial pair-production optical depth, $\tau_{\gamma\gamma} \gg 1$), meaning that absorption is still prominent at energies below $\epsilon_{\text{th,disc}} \sim 300 \text{ MeV}$ by seed photons in the exponential tail of the disc spectrum.

### 6.2.2 Connection to this work

Given the long exposure time necessary for the gamma-ray detection (Zdziarski et al. 2017), the prospects for examining temporal signatures of Klein–Nishina reconnection in high/soft BHXRBs at the energy scales near $\gamma_{\text{KN,disc}}m_e c^2$, where Klein–Nishina effects are likely most prominent, are not promising. At lower, X-ray energies,
millisecond variability was detected from the high/soft state of Cyg X-1 by Gierliński & Zdziarski (2003). Unlike the case of the IC(BLR) and IC(HDR) scenarios discussed in the preceding section, these X-ray observations exhibit a strong harder-when-brighter correlation between the flaring flux and the spectral shape. This could hint at the importance of synchrotron cooling in this system. Alternatively, due again to the high radiative compactness, the radiative cooling time-scale may be so short that, even in the presence of Klein–Nishina and pair-production effects, the reconnection-energized particles cool down instantaneously on the observationally resolved time-scales, masking the radiatively stabilized (on presumably faster time-scales; Section 4.5) spectrum.

Using our previous estimate of the coronal magnetization (Mehlhaff et al. 2021), $\sigma_{c,0} \sim 10^{2-4}$, we can infer a pair yield control parameter of $\gamma_{pp, \max}/4\sigma_{c,0} \sim 30\gamma_{KN, disc}/4\sigma_{c,0} \sim 0.08-8$, indicating Klein–Nishina reconnection as a potentially important source of electron–positron pairs in the coronae of BHXRBs in their high/soft states. However, due to the high radiative compactness, another likely source of pair production is the collision of disc-Comptonized photons with each other (Beloborodov 2017). This interaction occurs predominantly among photons both with energies $< \epsilon_{th, disc} = m_e^2 c^4 / \epsilon_{disc}$, since those with energies beyond $\epsilon_{th, disc}$ are most likely to be absorbed by the much denser disc-supplied radiation field. The total pair yield is then the sum of the contributions from both channels.

### 6.2.3 Synopsis of BHXRB discussion

To summarize, the case for coronal Klein–Nishina reconnection in the high/soft states of BHXRBs is excellent on theoretical grounds (Mehlhaff et al. 2021). It is likely, for example, that reconnection in this regime contributes to the pair content in these objects (Section 5). However, given that even in an exceptionally bright source, Cyg X-1, the gamma-ray flux is far too low to probe corona-scale reconnection-powered flaring variability, firm connections with the characteristic temporal signatures of Klein–Nishina reconnection (Section 4.5) will likely remain out of reach for the foreseeable future.

## 6.3 The M87* magnetosphere

### 6.3.1 Background information

The M87 galaxy has been monitored in TeV gamma-rays for nearly 20 yr (e.g. Aharonian et al. 2006; MAGIC Collaboration 2020; EHT MWL Science Working Group et al. 2021). This includes three major flaring periods, one in 2006 (Aharonian et al. 2006), one in 2008 (MAGIC Collaboration 2008; Acciari et al. 2009), and one in 2010 (Abramowski et al. 2012; Aliu et al. 2012), that exhibited rapid, one-day variability time-scales – of order the light crossing time of the galaxy’s central supermassive black hole, M87*. For two of these TeV-loud periods, observations at longer, more spatially resolved wavelengths revealed (nearly) contemporaneous flux enhancements from the galactic core (Acciari et al. 2009; Abramowski et al. 2012). Thus, variability considerations combined with the multiwavelength context suggest M87* itself (more precisely, its immediate plasma environment) as a viable site of TeV emission.

Direct polarized imaging of M87* suggests that accretion proceeds in a so-called magnetically arrested (MAD) state (Event Horizon Telescope Collaboration 2021), characterized, as revealed in large part by magnetohydrodynamics simulations, by quasi-periodic cycles of gradual accumulation of magnetic flux onto the black hole, eventual flux saturation, and finally abrupt and violent flux expulsion

(e.g. Igumenshchev 2008; Tchekhovskoy, Narayan & McKinney 2011; Avara, McKinney & Reynolds 2016; Ripperda et al. 2020, 2022; Chashkina, Bromberg & Levinson 2021; Porth et al. 2021; Mizuno 2022; Scepi et al. 2022). These expulsion events are mediated in the black-hole magnetosphere by large-scale reconnecting current sheets (Ripperda et al. 2022). The current sheets are irradiated by the larger-scale accretion flow, which provides a target photon bath for reconnection-accelerated particles to Comptonize up to the observed TeV energies (Hakobyan et al. 2023b). If particle acceleration is efficient enough, the tail of the Comptonized radiation spectrum falls above pair-production threshold with the accretion flow seed photons, triggering potentially copious pair creation (Crinquand et al. 2021, 2022). Thus, reconnection in the M87* magnetosphere is: (1) a potential source of the observed TeV emission, and (2) likely coupled to the same radiative physics as treated in this study.

### 6.3.2 Connection to this work

We first demonstrate that, in line with point (2) above, the Klein–Nishina hierarchy (50) is likely realized in the M87* magnetosphere, pointing to the importance of Klein–Nishina effects on reconnection there. We adopt a black hole mass $M_{BH} = 6 \times 10^9 M_{\odot}$ (Gebhardt et al. 2011; Event Horizon Telescope Collaboration 2019) and horizon-scale reconnecting magnetic field strength $B_{BH} = 100$ G (Yao et al. 2021; Ripperda et al. 2022). We also use the seed photon radiation field calculated by Yao et al. (2021) that is produced in a MAD accretion state and illuminates the inner black hole magnetosphere with overall energy density $U_{BH} = 10 \text{ erg cm}^{-3}$ and characteristic photon energy $\epsilon_{BH} \sim 10^{-2} \text{ eV}$. Then, using equation (7), we arrive at

$$ \gamma_{KN,BH} = \frac{m_e c^2}{4 \epsilon_{BH}} \sim 1 \times 10^7, \qquad (60) $$

while, via equation (48), we estimate

$$ \gamma_{\text{rad,T}}^{\text{(BH)}} = \sqrt{\frac{0.3eB_{BH}}{4\sigma_{\text{T}}U_{BH}}} \sim 2 \times 10^7. \qquad (61) $$

In addition, adopting a fiducial reconnection current sheet length of order the gravitational radius, $r_g = GM_{BH}/c^2 \simeq 9 \times 10^{14}$ cm, gives, by equation (38),

$$ \gamma_{\text{max,BH}} = \frac{0.1eB_{BH}r_g}{m_e c^2} \sim 5 \times 10^{12} \qquad (62) $$

and, through $\gamma_{\text{rad,T}}^2 = \gamma_{\text{max}}\gamma_{\text{cool}}$ [equation (49)],

$$ \gamma_{\text{cool,BH}} \sim 100, \qquad (63) $$

which can be used, via equation (51), to evaluate

$$ \tau_{\gamma\gamma,\text{BH}} = \frac{3\gamma_{\text{KN,BH}}}{5\gamma_{\text{cool,BH}}} \sim 7 \times 10^4. \qquad (64) $$

Thus, the Klein–Nishina scale hierarchy (50), $\gamma_{\text{cool}} < \gamma_{\text{KN}} < \gamma_{\text{rad,T}} < \gamma_{\text{max}}$, likely holds (though $\gamma_{\text{KN,BH}}$ in these estimates somewhat encroaches on $\gamma_{\text{rad,T}}^{\text{(BH)}}$ from below).

Although the above estimates suggest a coupling of reconnection to Klein–Nishina physics, synchrotron losses, unlike in our simulations, are non-negligible in the M87* magnetosphere. To show this, we note that the synchrotron power radiated per particle is

$$ P_{\text{syn}}(\gamma) = 2\sigma_{\text{T}}\beta^2\gamma^2 U_{\text{B}} \sin^2 \alpha, \qquad (65) $$

where $U_{\text{B}} \equiv B^2/8\pi$ is the local magnetic field energy density and $\alpha$ is the pitch angle between the radiating particle’s velocity and
millisecond variability was detected from the high/soft state of Cyg X-1 by Gierliński & Zdziarski (2003). Unlike the case of the IC(BLR) and IC(HDR) scenarios discussed in the preceding section, these X-ray observations exhibit a strong harder-when-brighter correlation between the flaring flux and the spectral shape. This could hint at the importance of synchrotron cooling in this system. Alternatively, due again to the high radiative compactness, the radiative cooling time-scale may be so short that, even in the presence of Klein–Nishina and pair-production effects, the reconnection-energized particles cool down instantaneously on the observationally resolved time-scales, masking the radiatively stabilized (on presumably faster time-scales; Section 4.5) spectrum.

Using our previous estimate of the coronal magnetization (Mehlhaff et al. 2021), $\sigma_{c,0} \sim 10^{2-4}$, we can infer a pair yield control parameter of $\gamma_{pp, \max}/4\sigma_{c,0} \sim 30\gamma_{KN, disc}/4\sigma_{c,0} \sim 0.08-8$, indicating Klein–Nishina reconnection as a potentially important source of electron–positron pairs in the coronae of BHXRBs in their high/soft states. However, due to the high radiative compactness, another likely source of pair production is the collision of disc-Comptonized photons with each other (Beloborodov 2017). This interaction occurs predominantly among photons both with energies $< \epsilon_{th, disc} = m_e^2 c^4 / \epsilon_{disc}$, since those with energies beyond $\epsilon_{th, disc}$ are most likely to be absorbed by the much denser disc-supplied radiation field. The total pair yield is then the sum of the contributions from both channels.

### 6.2.3 Synopsis of BHXRB discussion

To summarize, the case for coronal Klein–Nishina reconnection in the high/soft states of BHXRBs is excellent on theoretical grounds (Mehlhaff et al. 2021). It is likely, for example, that reconnection in this regime contributes to the pair content in these objects (Section 5). However, given that even in an exceptionally bright source, Cyg X-1, the gamma-ray flux is far too low to probe corona-scale reconnection-powered flaring variability, firm connections with the characteristic temporal signatures of Klein–Nishina reconnection (Section 4.5) will likely remain out of reach for the foreseeable future.

## 6.3 The M87* magnetosphere

### 6.3.1 Background information

The M87 galaxy has been monitored in TeV gamma-rays for nearly 20 yr (e.g. Aharonian et al. 2006; MAGIC Collaboration 2020; EHT MWL Science Working Group et al. 2021). This includes three major flaring periods, one in 2006 (Aharonian et al. 2006), one in 2008 (MAGIC Collaboration 2008; Acciari et al. 2009), and one in 2010 (Abramowski et al. 2012; Aliu et al. 2012), that exhibited rapid, one-day variability time-scales – of order the light crossing time of the galaxy’s central supermassive black hole, M87*. For two of these TeV-loud periods, observations at longer, more spatially resolved wavelengths revealed (nearly) contemporaneous flux enhancements from the galactic core (Acciari et al. 2009; Abramowski et al. 2012). Thus, variability considerations combined with the multiwavelength context suggest M87* itself (more precisely, its immediate plasma environment) as a viable site of TeV emission.

Direct polarized imaging of M87* suggests that accretion proceeds in a so-called magnetically arrested (MAD) state (Event Horizon Telescope Collaboration 2021), characterized, as revealed in large part by magnetohydrodynamics simulations, by quasi-periodic cycles of gradual accumulation of magnetic flux onto the black hole, eventual flux saturation, and finally abrupt and violent flux expulsion

(e.g. Igumenshchev 2008; Tchekhovskoy, Narayan & McKinney 2011; Avara, McKinney & Reynolds 2016; Ripperda et al. 2020, 2022; Chashkina, Bromberg & Levinson 2021; Porth et al. 2021; Mizuno 2022; Scepi et al. 2022). These expulsion events are mediated in the black-hole magnetosphere by large-scale reconnecting current sheets (Ripperda et al. 2022). The current sheets are irradiated by the larger-scale accretion flow, which provides a target photon bath for reconnection-accelerated particles to Comptonize up to the observed TeV energies (Hakobyan et al. 2023b). If particle acceleration is efficient enough, the tail of the Comptonized radiation spectrum falls above pair-production threshold with the accretion flow seed photons, triggering potentially copious pair creation (Crinquand et al. 2021, 2022). Thus, reconnection in the M87* magnetosphere is: (1) a potential source of the observed TeV emission, and (2) likely coupled to the same radiative physics as treated in this study.

### 6.3.2 Connection to this work

We first demonstrate that, in line with point (2) above, the Klein–Nishina hierarchy (50) is likely realized in the M87* magnetosphere, pointing to the importance of Klein–Nishina effects on reconnection there. We adopt a black hole mass $M_{BH} = 6 \times 10^9 M_{\odot}$ (Gebhardt et al. 2011; Event Horizon Telescope Collaboration 2019) and horizon-scale reconnecting magnetic field strength $B_{BH} = 100$ G (Yao et al. 2021; Ripperda et al. 2022). We also use the seed photon radiation field calculated by Yao et al. (2021) that is produced in a MAD accretion state and illuminates the inner black hole magnetosphere with overall energy density $U_{BH} = 10 \text{ erg cm}^{-3}$ and characteristic photon energy $\epsilon_{BH} \sim 10^{-2} \text{ eV}$. Then, using equation (7), we arrive at

$$ \gamma_{KN,BH} = \frac{m_e c^2}{4\epsilon_{BH}} \sim 1 \times 10^7, \qquad (60) $$

while, via equation (48), we estimate

$$ \gamma_{\text{rad,T}}^{\text{(BH)}} = \sqrt{\frac{0.3eB_{BH}}{4\sigma_{\text{T}}U_{\text{BH}}}} \sim 2 \times 10^7. \qquad (61) $$

In addition, adopting a fiducial reconnection current sheet length of order the gravitational radius, $r_g = GM_{BH}/c^2 \simeq 9 \times 10^{14}$ cm, gives, by equation (38),

$$ \gamma_{\text{max,BH}} = \frac{0.1eB_{\text{BH}}r_g}{m_e c^2} \sim 5 \times 10^{12} \qquad (62) $$

and, through $\gamma_{\text{rad,T}}^2 = \gamma_{\text{max}}\gamma_{\text{cool}}$ [equation (49)],

$$ \gamma_{\text{cool,BH}} \sim 100, \qquad (63) $$

which can be used, via equation (51), to evaluate

$$ \tau_{\gamma\gamma,\text{BH}} = \frac{3\gamma_{\text{KN,BH}}}{5\gamma_{\text{cool,BH}}} \sim 7 \times 10^4. \qquad (64) $$

Thus, the Klein–Nishina scale hierarchy (50), $\gamma_{\text{cool}} < \gamma_{\text{KN}} < \gamma_{\text{rad,T}} < \gamma_{\text{max}}$, likely holds (though $\gamma_{\text{KN,BH}}$ in these estimates somewhat encroaches on $\gamma_{\text{rad,T}}^{\text{(BH)}}$ from below).

Although the above estimates suggest a coupling of reconnection to Klein–Nishina physics, synchrotron losses, unlike in our simulations, are non-negligible in the M87* magnetosphere. To show this, we note that the synchrotron power radiated per particle is

$$ P_{\text{syn}}(\gamma) = 2\sigma_{\text{T}}c^2\gamma^2U_{\text{B}}\sin^2\alpha, \qquad (65) $$

where $U_B \equiv B^2/8\pi$ is the local magnetic field energy density and $\alpha$ is the pitch angle between the radiating particle’s velocity and
the local magnetic field. An approximate lower bound on the ratio of synchrotron radiation to IC radiation is then, via equations (65), (41), and (42),

$$ \frac{P_{\text{syn}}(\gamma)}{P_{\text{IC}}(\gamma)} \geq \frac{P_{\text{syn}}(\gamma)}{P_{\text{T}}(\gamma)} = \frac{3U_{\text{B}} \sin^2 \alpha}{2U_{\text{BH}}}, \qquad (66) $$

where, in the first step, we ignore potential Klein–Nishina suppression of the IC cross-section. Thus, excepting special regions (such as reconnection X-points) where the magnetic field weakens ($B^2 \ll B_{BH}^2$) or the pitch angle shrinks ($\sin^2\alpha \ll 1$), synchrotron cooling typically exceeds IC losses by the factor $B_{BH}^2/8\pi U_{BH} \sim 40$.

Having recognized their importance, let us now reason about how strong synchrotron losses may change the observable aspects of reconnection as unveiled in this work. The fact that equation (65) is identical to the Thomson IC formula (41) but with $U_{ph}$ replaced by $3U_{B}\sin^2\alpha/2$ suggests that efficient synchrotron radiation and Thomson IC cooling may play dynamically similar roles. We might therefore expect, akin to our simulations of Thomson IC reconnection, radiation from the M87* magnetosphere to exhibit a much tighter coupling between spectral shape and total luminosity – that is, ‘harder-when-brighter’ – than when only Klein–Nishina IC cooling and pair production dominate the radiative physics (cf. Section 4.5). This would hold whether one observes the synchrotron or the IC radiation, as both probe the same underlying distribution of reconnection-energized particles. Such an argument for a prominent harder-when-brighter trend is, however, in mild tension with the few M87 TeV flares, which indicate this trend only mildly (Aliu et al. 2012) and not in every event (Acciari et al. 2010).

We next consider the potential ramifications of pair production between IC photons and the seed photons from the accretion flow. To provide a context for this discussion, we briefly summarize recent semi-analytic works, Kimura et al. (2022), Hakobyan et al. (2023b), and Chen et al. (2023), on the matter content of the putative magnetospheric M87* reconnection layer. Hakobyan et al. (2023b) conduct a detailed analysis of the possible radiation and pair-production channels in this system. They predict that the brightest radiation emerges through the synchrotron band and peaks between roughly 1 and 20 MeV. Though the synchrotron photons are not above pair-production threshold with the ambient accretion-flow photons, they are above threshold with themselves. There is, however, little chance for an individual synchrotron photon to be absorbed by another since the optical depth, $\tau_{\gamma\gamma}^{(\text{syn})} \sim 10^{-4}$, presented by the synchrotron radiation is small. Thus, whereas the high optical thickness furnished by the accretion-flow light causes nearly all of the above-threshold IC photons to be absorbed close to the reconnection current sheet, covering it with a thin pair coat, only a small fraction of the synchrotron radiation is absorbed, leading to diffuse pair production throughout the magnetosphere. The resulting pairs feed the reconnection layer with a highly magnetized, $\sigma_{c,0} \sim 5 \times 10^7$, plasma. The pictures presented by Kimura et al. (2022) and Chen et al. (2023) are similar, but those authors predict more copious synchrotron–synchrotron pair production, yielding the reduced magnetizations, $\sigma_{c,0} \sim 8.7 \times 10^4$ and $\sim 6 \times 10^4$, respectively.

In the context of our reconnection model, the diffuse synchrotron–synchrotron pair production predicted by Kimura et al. (2022), Hakobyan et al. (2023b), and Chen et al. (2023) supplies the background plasma on top of which pair production between IC photons and the radiation impinging from the accretion flow (hereafter, *IC pair production*) may add supplementary pairs localized

near the current sheet. If we naively apply our findings from Section 5 on the reconnection-powered IC pair yield (temporarily ignoring potential changes due to synchrotron cooling), we infer that a substantial amplification of the background (synchrotron–synchrotron) pair density requires $4\sigma_{c,0} \geq \gamma_{\text{pp,max}} \simeq 30\gamma_{\text{KN}} \sim 4 \times 10^8$. Comparing to $\sigma_{c,0} \sim 5 \times 10^7$ from Hakobyan et al. (2023b), we estimate $\gamma_{\text{pp,max}}/4\sigma_{c,0} \simeq 2$, which is slightly too high to attain an order-unity IC pair yield according to this study (Section 5). One should keep in mind, however, that the uncertainties on these order-of-magnitude estimates are high, while even a reduction in the ratio $\gamma_{\text{pp,max}}/4\sigma_{c,0}$ by a factor of 4 is enough to bring the IC pair yield up to unity. However, if synchrotron–synchrotron pair production is more efficient, as predicted by Kimura et al. (2022) and Chen et al. (2023), then the inflowing plasma is probably not sufficiently magnetized to lead to appreciable IC pair production. In the event that a substantial IC pair yield is realized, it induces a transverse density gradient on the magnetic field lines forming the jet funnel, with a higher plasma load on the field lines that participate in reconnection near the jet walls. This would then be potentially important for jet-boundary interactions (e.g. Ripperda et al. 2020; Sironi, Rowan & Narayan 2021; El Mellah et al. 2022, 2023; Chow et al. 2023a, b), such as may power observed limb-brightening (Ly, Walker & Junor 2007; Kim et al. 2018; Walker et al. 2018).

In the preceding discussion, we simply apply the forecasted pair yield obtained from this study (Section 5) to the M87* case. We now anticipate, as we did earlier for radiative signatures, how strong synchrotron losses may alter the picture of pair production as revealed by this work. First, rapid synchrotron cooling disposes the IC pairs born into the upstream region to cool down much more before entering the reconnection layer, which may change them from a hot population to a merely warm or even cold one (thus reducing their contribution to the upstream plasma energy density; cf. Section 4.6). Second, the IC pair yield may be restricted since high-energy particles are likely to emit fewer pair-producing IC photons before radiating away their energy as synchrotron light. Interestingly, however, the global magnetospheric simulations of Crinquand et al. (2022), which include the same processes of IC emission and pair creation studied here, find that IC pair production is efficient enough to fuel the reconnection layer with plasma even when synchrotron losses are made as strong as numerically possible. This is likely connected to the fact that TeV emitting particles are accelerated near reconnection X-points (Hakobyan et al. 2023b) where their synchrotron losses are suppressed. Despite these arguments, then, previous numerical work hints that synchrotron radiation does not quench IC pair production.

### 6.3.3 *Synopsis of M87*\* discussion

To summarize, the basic Klein–Nishina hierarchy (50) seems to be satisfied in the M87* magnetosphere. Nevertheless, synchrotron cooling is also an important radiative mechanism for reconnection in this system. Because we neglect this in our simulations, we cannot directly apply our results to the M87* case. We therefore, in this section, attempt to anticipate the modifications synchrotron cooling could induce on the radiative signatures and pair yield of Klein–Nishina reconnection as uncovered in this work (Sections 4.5 and 5). When formulating these expectations, we find that the resulting conclusions are not convincingly supported by recent studies. For example, we argue that synchrotron losses might induce a tighter correlation between spectral shape and total luminosity, making
the local magnetic field. An approximate lower bound on the ratio of synchrotron radiation to IC radiation is then, via equations (65), (41), and (42),

$$ \frac{P_{\text{syn}}(\gamma)}{P_{\text{IC}}(\gamma)} \ge \frac{P_{\text{syn}}(\gamma)}{P_{\text{T}}(\gamma)} = \frac{3U_{\text{B}} \sin^2 \alpha}{2U_{\text{BH}}}, \quad (66) $$

where, in the first step, we ignore potential Klein–Nishina suppression of the IC cross-section. Thus, excepting special regions (such as reconnection X-points) where the magnetic field weakens ($B^2 \ll B_{BH}^2$) or the pitch angle shrinks ($\sin^2\alpha \ll 1$), synchrotron cooling typically exceeds IC losses by the factor $B_{BH}^2/8\pi U_{BH} \sim 40$.

Having recognized their importance, let us now reason about how strong synchrotron losses may change the observable aspects of reconnection as unveiled in this work. The fact that equation (65) is identical to the Thomson IC formula (41) but with $U_{ph}$ replaced by $3U_{B}\sin^2\alpha/2$ suggests that efficient synchrotron radiation and Thomson IC cooling may play dynamically similar roles. We might therefore expect, akin to our simulations of Thomson IC reconnection, radiation from the M87* magnetosphere to exhibit a much tighter coupling between spectral shape and total luminosity – that is, ‘harder-when-brighter’ – than when only Klein–Nishina IC cooling and pair production dominate the radiative physics (cf. Section 4.5). This would hold whether one observes the synchrotron or the IC radiation, as both probe the same underlying distribution of reconnection-energized particles. Such an argument for a prominent harder-when-brighter trend is, however, in mild tension with the few M87 TeV flares, which indicate this trend only mildly (Aliu et al. 2012) and not in every event (Acciari et al. 2010).

We next consider the potential ramifications of pair production between IC photons and the seed photons from the accretion flow. To provide a context for this discussion, we briefly summarize recent semi-analytic works, Kimura et al. (2022), Hakobyan et al. (2023b), and Chen et al. (2023), on the matter content of the putative magnetospheric M87* reconnection layer. Hakobyan et al. (2023b) conduct a detailed analysis of the possible radiation and pair-production channels in this system. They predict that the brightest radiation emerges through the synchrotron band and peaks between roughly 1 and 20 MeV. Though the synchrotron photons are not above pair-production threshold with the ambient accretion-flow photons, they are above threshold with themselves. There is, however, little chance for an individual synchrotron photon to be absorbed by another since the optical depth, $\tau_{\gamma\gamma}^{(\text{syn})} \sim 10^{-4}$, presented by the synchrotron radiation is small. Thus, whereas the high optical thickness furnished by the accretion-flow light causes nearly all of the above-threshold IC photons to be absorbed close to the reconnection current sheet, covering it with a thin pair coat, only a small fraction of the synchrotron radiation is absorbed, leading to diffuse pair production throughout the magnetosphere. The resulting pairs feed the reconnection layer with a highly magnetized, $\sigma_{c,0} \sim 5 \times 10^7$, plasma. The pictures presented by Kimura et al. (2022) and Chen et al. (2023) are similar, but those authors predict more copious synchrotron–synchrotron pair production, yielding the reduced magnetizations, $\sigma_{c,0} \sim 8.7 \times 10^4$ and $\sim 6 \times 10^4$, respectively.

In the context of our reconnection model, the diffuse synchrotron–synchrotron pair production predicted by Kimura et al. (2022), Hakobyan et al. (2023b), and Chen et al. (2023) supplies the background plasma on top of which pair production between IC photons and the radiation impinging from the accretion flow (hereafter, IC pair production) may add supplementary pairs localized

near the current sheet. If we naively apply our findings from Section 5 on the reconnection-powered IC pair yield (temporarily ignoring potential changes due to synchrotron cooling), we infer that a substantial amplification of the background (synchrotron–synchrotron) pair density requires $4\sigma_{c,0} \ge \gamma_{\text{pp,max}} \simeq 30\gamma_{\text{KN}} \sim 4 \times 10^8$. Comparing to $\sigma_{c,0} \sim 5 \times 10^7$ from Hakobyan et al. (2023b), we estimate $\gamma_{\text{pp,max}}/4\sigma_{c,0} \simeq 2$, which is slightly too high to attain an order-unity IC pair yield according to this study (Section 5). One should keep in mind, however, that the uncertainties on these order-of-magnitude estimates are high, while even a reduction in the ratio $\gamma_{\text{pp,max}}/4\sigma_{c,0}$ by a factor of 4 is enough to bring the IC pair yield up to unity. However, if synchrotron–synchrotron pair production is more efficient, as predicted by Kimura et al. (2022) and Chen et al. (2023), then the inflowing plasma is probably not sufficiently magnetized to lead to appreciable IC pair production. In the event that a substantial IC pair yield is realized, it induces a transverse density gradient on the magnetic field lines forming the jet funnel, with a higher plasma load on the field lines that participate in reconnection near the jet walls. This would then be potentially important for jet-boundary interactions (e.g. Ripperda et al. 2020; Sironi, Rowan & Narayan 2021; El Mellah et al. 2022, 2023; Chow et al. 2023a, b), such as may power observed limb-brightening (Ly, Walker & Junor 2007; Kim et al. 2018; Walker et al. 2018).

In the preceding discussion, we simply apply the forecasted pair yield obtained from this study (Section 5) to the M87* case. We now anticipate, as we did earlier for radiative signatures, how strong synchrotron losses may alter the picture of pair production as revealed by this work. First, rapid synchrotron cooling disposes the IC pairs born into the upstream region to cool down much more before entering the reconnection layer, which may change them from a hot population to a merely warm or even cold one (thus reducing their contribution to the upstream plasma energy density; cf. Section 4.6). Second, the IC pair yield may be restricted since high-energy particles are likely to emit fewer pair-producing IC photons before radiating away their energy as synchrotron light. Interestingly, however, the global magnetospheric simulations of Crinquand et al. (2022), which include the same processes of IC emission and pair creation studied here, find that IC pair production is efficient enough to fuel the reconnection layer with plasma even when synchrotron losses are made as strong as numerically possible. This is likely connected to the fact that TeV emitting particles are accelerated near reconnection X-points (Hakobyan et al. 2023b) where their synchrotron losses are suppressed. Despite these arguments, then, previous numerical work hints that synchrotron radiation does not quench IC pair production.

### 6.3.3 *Synopsis of M87*^* discussion*

To summarize, the basic Klein–Nishina hierarchy (50) seems to be satisfied in the M87* magnetosphere. Nevertheless, synchrotron cooling is also an important radiative mechanism for reconnection in this system. Because we neglect this in our simulations, we cannot directly apply our results to the M87* case. We therefore, in this section, attempt to anticipate the modifications synchrotron cooling could induce on the radiative signatures and pair yield of Klein–Nishina reconnection as uncovered in this work (Sections 4.5 and 5). When formulating these expectations, we find that the resulting conclusions are not convincingly supported by recent studies. For example, we argue that synchrotron losses might induce a tighter correlation between spectral shape and total luminosity, making
the radiative signatures of reconnection more Thomson-like (with a more prominent harder-when-brighter trend). However, observational evidence for this is somewhat ambiguous, providing, at best, limited support (Aliu et al. 2012) and, at worst, mild tension (Acciari et al. 2010; Beilicke & VERITAS Collaboration 2012). We further argue that synchrotron cooling may shut down the IC pair yield, but this seems to be in conflict with first-principles global simulations (Crinquand et al. 2022). In short, comparing to observations (Acciari et al. 2010; Aliu et al. 2012; Beilicke & VERITAS Collaboration 2012) and simulations (Crinquand et al. 2022) suggests that our simple arguments do not satisfactorily generalize the results of this study to reconnection where efficient synchrotron cooling also comes into play. This creates fertile ground for future reconnection simulations to self-consistently incorporate synchrotron radiation on an equal footing with Klein–Nishina effects to illuminate how the two conspire together to shape the gamma-ray signatures and pair yield in the M87* magnetosphere.

## 6.4 Gamma-ray binaries

### 6.4.1 Background information

Gamma-ray binaries consist of a relativistic compact object (neutron star or black hole) and a massive stellar companion (generally of type O or Be) and are defined by a spectral energy density peaking, in the $\nu F(\nu)$ representation, above 1 MeV (e.g. Dubus 2013; Dubus et al. 2017). Of the handful (less than a dozen) of known gamma-ray binaries, only two are directly observed to host pulsars (Dubus 2013; Tam et al. 2018; Chernyakova et al. 2019; Chernyakova & Malyshev 2020). Nevertheless, several general observed features suggest that the compact object in these systems is generically a rotation-powered pulsar (Dubus 2006, 2013).

Adopting this view, two often-invoked gamma-ray emission sites are the pre- and post-shocked pulsar wind, where the shock in question interfaces between the winds of the pulsar and the massive companion (e.g. Kirk, Ball & Skjæraasen 1999; Ball & Kirk 2000; Ball & Dodd 2001; Sierpowska & Bednarek 2005; Cerutti et al. 2008; Sierpowska-Bartosik & Bednarek 2008; Khangulyan et al. 2012) and not, as would be the case in isolated pulsars, between the pulsar wind and the interstellar medium. As argued in the review by Dubus (2013), the fact that the high-energy (0.1–10 GeV) spectra of gamma-ray binaries are often similar in terms of slope and cut-off to isolated pulsars could hint at a similar emission mechanism between the two object classes, motivating an investigation of the unshocked pulsar wind, and perhaps even of the pulsar magnetosphere, as gamma-ray production sites in gamma-ray binaries. However, the GeV spectra also exhibit modulations on the binary orbital period, which would seem to disfavour the magnetosphere as the dominant emission zone, since it is insensitive to the orbit of the binary (Dubus 2013). Therefore, in this section, we consider the possibility that the unshocked pulsar wind significantly contributes to the observed high-energy gamma-rays in gamma-ray binaries (cf. Ball & Kirk 2000; Ball & Dodd 2001; Cerutti et al. 2008; Khangulyan et al. 2012).

The unshocked pulsar wind behaves exactly as that of an isolated pulsar except for one key difference: it is illuminated from beyond by the hot massive companion star. The temperature, $T_*$ ~ 40 000 K, of the companion's surface produces a characteristic blackbody photon energy,

$$ \epsilon_* \sim 3k_B T_* \sim 10 \text{ eV}, \quad (67) $$

and radiation energy density,

$$ U_* = \frac{\sigma_{SB} T_*^4}{c} \left(\frac{R_*}{d}\right)^2 \sim 1 \times 10^3 \text{ erg cm}^{-3}, \quad (68) $$

where $d = 0.1$ au is the typical separation at periastron and $R_* = 10 R_\odot$ is the stellar radius (Dubus 2013). Equation (67) implies a critical Klein–Nishina Lorentz factor in the pulsar wind of

$$ \gamma_{KN,w} \equiv \frac{m_e c^2}{4 \epsilon_*} \sim 1 \times 10^4 \quad (69) $$

and that pair production becomes possible above the threshold

$$ \epsilon_{\text{th,w}} = m_e^2 c^4 / \epsilon_* = 30 \text{ GeV}. \quad (70) $$

The typical photon energy emitted by particles with $\gamma = \gamma_{KN,w}$ is then $\gamma_{KN,w} m_e c^2 / 2 = \epsilon_{\text{th,w}} / 8 \sim 3$ GeV.

### 6.4.2 Connection to this work

Let us consider the effect that illumination by the companion may have on the pulsar wind. Here, we adopt the theoretical picture (Coroniti 1990; Michel 1994; Bogovalov 1999; Lyubarsky & Kirk 2001; Kirk & Skjæraasen 2003; Lyubarsky 2003; Kirk, Lyubarsky & Petri 2009), brought into sharper focus by recent first-principles kinetic simulations (Cerutti & Philippov 2017; Philippov & Spitkovsky 2018; Cerutti et al. 2020), that this wind is not purely cold, but is instead striped – laced with a large-scale reconnecting current sheet that expands radially outward while undulating about the pulsar’s rotational equator. The angular excursions about the equator approximately equal the obliquity angle between the magnetic and spin pulsar axes. Reconnection converts the outgoing wind Poynting flux into a combination of bulk acceleration and NTPA.

An isolated pulsar is not bathed in the intense light of a companion star, and so cooling of accelerated wind particles remains dominated by synchrotron losses. In the present case, however, cooling via IC scattering of the intense radiation bath (68) likely outpaces synchrotron cooling once the local magnetic field energy density falls below $U_*$. This occurs at a critical magnetic field strength,

$$ B_{\text{IC}} \equiv \sqrt{8\pi U_*} \sim 200 \text{ G}. \quad (71) $$

For a pulsar rotational period $P \sim 100$ ms (characteristic of the two confirmed pulsars in gamma-ray binaries: PSR J2032 + 4127, Abdo et al. 2009a; and PSR B1259–63, Johnston et al. 1992) and a surface magnetic field $B_{\text{surf}} = 10^{12}$ G, the magnetic field is diluted to $B_{\text{LC}} \simeq B_{\text{surf}}(R_{\text{psr}}/R_{\text{LC}})^3 \sim 9 \times 10^3$ G, at the light cylinder, $R_{\text{LC}} = cP/2\pi$, where $R_{\text{psr}} = 10$ km is the assumed pulsar radius. Beyond the light cylinder, the magnetic field falls off slower, as $R_{\text{LC}}/R$ where $R$ is the cylindrical radius measured from the pulsar’s spin axis, and, hence, even for the strong surface field $B_{\text{psr}} = 10^{12}$ G, the striped wind’s radiative losses become IC dominated at $R/R_{\text{LC}} \simeq (B_{\text{surf}}/B_{\text{IC}})(R_{\text{psr}}/R_{\text{LC}})^3 \sim 60$. This is far before the shock with the companion’s wind, which is expected to occur on scales $R \sim d \sim 10^3 R_{\text{LC}}$ (Dubus 2013). Moreover, $R \sim 60 R_{\text{LC}}$ is also before the point where the pulsar wind’s electromagnetic flux is expected to be fully dissipated, which kinetic simulations anticipate at $R/R_{\text{LC}}$ roughly between $10^2$ and $10^4$ (Cerutti et al. 2020). Taken together, these estimates suggest that most of the pulsar wind’s dissipation takes place before the shock with the companion’s wind, but after the critical radius where IC losses surpass synchrotron losses.

Not only does most of the pulsar wind dissipate through reconnection subject to strong IC cooling, but, as we show now, reconnection may occur in the Klein–Nishina regime of this study, attaining the
the radiative signatures of reconnection more Thomson-like (with a more prominent harder-when-brighter trend). However, observational evidence for this is somewhat ambiguous, providing, at best, limited support (Aliu et al. 2012) and, at worst, mild tension (Acciari et al. 2010; Beilicke & VERITAS Collaboration 2012). We further argue that synchrotron cooling may shut down the IC pair yield, but this seems to be in conflict with first-principles global simulations (Crinquand et al. 2022). In short, comparing to observations (Acciari et al. 2010; Aliu et al. 2012; Beilicke & VERITAS Collaboration 2012) and simulations (Crinquand et al. 2022) suggests that our simple arguments do not satisfactorily generalize the results of this study to reconnection where efficient synchrotron cooling also comes into play. This creates fertile ground for future reconnection simulations to self-consistently incorporate synchrotron radiation on an equal footing with Klein–Nishina effects to illuminate how the two conspire together to shape the gamma-ray signatures and pair yield in the M87* magnetosphere.

## 6.4 Gamma-ray binaries

### 6.4.1 Background information

Gamma-ray binaries consist of a relativistic compact object (neutron star or black hole) and a massive stellar companion (generally of type O or Be) and are defined by a spectral energy density peaking, in the $\nu F(\nu)$ representation, above 1 MeV (e.g. Dubus 2013; Dubus et al. 2017). Of the handful (less than a dozen) of known gamma-ray binaries, only two are directly observed to host pulsars (Dubus 2013; Tam et al. 2018; Chernyakova et al. 2019; Chernyakova & Malyshev 2020). Nevertheless, several general observed features suggest that the compact object in these systems is generically a rotation-powered pulsar (Dubus 2006, 2013).

Adopting this view, two often-invoked gamma-ray emission sites are the pre- and post-shocked pulsar wind, where the shock in question interfaces between the winds of the pulsar and the massive companion (e.g. Kirk, Ball & Skjæraasen 1999; Ball & Kirk 2000; Ball & Dodd 2001; Sierpowska & Bednarek 2005; Cerutti et al. 2008; Sierpowska-Bartosik & Bednarek 2008; Khangulyan et al. 2012) and not, as would be the case in isolated pulsars, between the pulsar wind and the interstellar medium. As argued in the review by Dubus (2013), the fact that the high-energy (0.1–10 GeV) spectra of gamma-ray binaries are often similar in terms of slope and cut-off to isolated pulsars could hint at a similar emission mechanism between the two object classes, motivating an investigation of the unshocked pulsar wind, and perhaps even of the pulsar magnetosphere, as gamma-ray production sites in gamma-ray binaries. However, the GeV spectra also exhibit modulations on the binary orbital period, which would seem to disfavour the magnetosphere as the dominant emission zone, since it is insensitive to the orbit of the binary (Dubus 2013). Therefore, in this section, we consider the possibility that the unshocked pulsar wind significantly contributes to the observed high-energy gamma-rays in gamma-ray binaries (cf. Ball & Kirk 2000; Ball & Dodd 2001; Cerutti et al. 2008; Khangulyan et al. 2012).

The unshocked pulsar wind behaves exactly as that of an isolated pulsar except for one key difference: it is illuminated from beyond by the hot massive companion star. The temperature, $T_*$ ~ 40 000 K, of the companion's surface produces a characteristic blackbody photon energy,

$$ \epsilon_* \sim 3k_B T_* \sim 10 \text{ eV}, \qquad (67) $$

and radiation energy density,

$$ U_* = \frac{\sigma_{SB} T_*^4}{c} \left(\frac{R_*}{d}\right)^2 \sim 1 \times 10^3 \text{ erg cm}^{-3}, \qquad (68) $$

where $d = 0.1$ au is the typical separation at periastron and $R_* = 10 R_\odot$ is the stellar radius (Dubus 2013). Equation (67) implies a critical Klein–Nishina Lorentz factor in the pulsar wind of

$$ \gamma_{KN,w} \equiv \frac{m_e c^2}{4 \epsilon_*} \sim 1 \times 10^4 \qquad (69) $$

and that pair production becomes possible above the threshold

$$ \epsilon_{\text{th,w}} = m_e^2 c^4 / \epsilon_* = 30 \text{ GeV}. \qquad (70) $$

The typical photon energy emitted by particles with $\gamma = \gamma_{KN,w}$ is then $\gamma_{KN,w} m_e c^2 / 2 = \epsilon_{\text{th,w}} / 8 \sim 3$ GeV.

### 6.4.2 Connection to this work

Let us consider the effect that illumination by the companion may have on the pulsar wind. Here, we adopt the theoretical picture (Coroniti 1990; Michel 1994; Bogovalov 1999; Lyubarsky & Kirk 2001; Kirk & Skjæraasen 2003; Lyubarsky 2003; Kirk, Lyubarsky & Petri 2009), brought into sharper focus by recent first-principles kinetic simulations (Cerutti & Philippov 2017; Philippov & Spitkovsky 2018; Cerutti et al. 2020), that this wind is not purely cold, but is instead striped – laced with a large-scale reconnecting current sheet that expands radially outward while undulating about the pulsar’s rotational equator. The angular excursions about the equator approximately equal the obliquity angle between the magnetic and spin pulsar axes. Reconnection converts the outgoing wind Poynting flux into a combination of bulk acceleration and NTPA.

An isolated pulsar is not bathed in the intense light of a companion star, and so cooling of accelerated wind particles remains dominated by synchrotron losses. In the present case, however, cooling via IC scattering of the intense radiation bath (68) likely outpaces synchrotron cooling once the local magnetic field energy density falls below $U_*$. This occurs at a critical magnetic field strength,

$$ B_{\text{IC}} \equiv \sqrt{8\pi U_*} \sim 200 \text{ G}. \qquad (71) $$

For a pulsar rotational period $P \sim 100$ ms (characteristic of the two confirmed pulsars in gamma-ray binaries: PSR J2032 + 4127, Abdo et al. 2009a; and PSR B1259–63, Johnston et al. 1992) and a surface magnetic field $B_{\text{surf}} = 10^{12}$ G, the magnetic field is diluted to $B_{\text{LC}} \simeq B_{\text{surf}}(R_{\text{psr}}/R_{\text{LC}})^3 \sim 9 \times 10^3$ G, at the light cylinder, $R_{\text{LC}} = cP / 2\pi$, where $R_{\text{psr}} = 10$ km is the assumed pulsar radius. Beyond the light cylinder, the magnetic field falls off slower, as $R_{\text{LC}} / R$ where $R$ is the cylindrical radius measured from the pulsar’s spin axis, and, hence, even for the strong surface field $B_{\text{psr}} = 10^{12}$ G, the striped wind’s radiative losses become IC dominated at $R / R_{\text{LC}} \simeq (B_{\text{surf}}/B_{\text{IC}})(R_{\text{psr}}/R_{\text{LC}})^3 \sim 60$. This is far before the shock with the companion’s wind, which is expected to occur on scales $R \sim d \sim 10^3 R_{\text{LC}}$ (Dubus 2013). Moreover, $R \sim 60 R_{\text{LC}}$ is also before the point where the pulsar wind’s electromagnetic flux is expected to be fully dissipated, which kinetic simulations anticipate at $R / R_{\text{LC}}$ roughly between $10^2$ and $10^4$ (Cerutti et al. 2020). Taken together, these estimates suggest that most of the pulsar wind’s dissipation takes place before the shock with the companion’s wind, but after the critical radius where IC losses surpass synchrotron losses.

Not only does most of the pulsar wind dissipate through reconnection subject to strong IC cooling, but, as we show now, reconnection may occur in the Klein–Nishina regime of this study, attaining the
critical scale hierarchy (50), $\gamma_{\text{cool}} < \gamma_{\text{KN}} < \gamma_{\text{rad, T}} < \gamma_{\text{max}}$. Assuming $B_{\text{LC}} = 9 \times 10^3$ G, $P = 100$ ms, and that the reconnecting magnetic field strength, $B_w$, is $B_w = B_{\text{LC}}R_{\text{LC}}/R$, as well as identifying the length, $L$, of the reconnection layer with the local radius, $R$, in the pulsar wind, we have, by equation (38),

$$ \gamma_{\text{max},w} \equiv \frac{0.1eB_w R}{m_e c^2} \sim 3 \times 10^8, \qquad (72) $$

which is independent of $R$. Similarly, we have, by equation (48),

$$ \gamma_{\text{rad, T}}^{(w)} \equiv \sqrt{\frac{0.3eB_w}{4\sigma_{\text{T}}U_{*}}} \sim 3 \times 10^6 \left(\frac{R}{60R_{\text{LC}}}\right)^{-1/2}. \qquad (73) $$

We then can use $\gamma_{\text{rad, T}}^2 = \gamma_{\text{cool}}\gamma_{\text{max}}$ [equation (49)], to estimate

$$ \gamma_{\text{cool},w} \sim 3 \times 10^4 \left( \frac{R}{60 R_{\text{LC}}} \right)^{-1}. \qquad (74) $$

Finally, using equation (51), we estimate the pair-production optical depth as

$$ \tau_{\gamma\gamma,w} \equiv \frac{3\gamma_{\text{KN},w}}{5\gamma_{\text{cool},w}} \sim 0.3 \left(\frac{R}{60R_{\text{LC}}}\right). \qquad (75) $$

Note that we have normalized $R$ to the critical radius where IC losses overtake synchrotron losses. The pulsar wind extends much farther than this, which pushes $\tau_{\gamma\gamma,w}$ above unity at the largest radii. The geometry of the striped wind is somewhat peculiar, however, in that the spacing between stripes is $\sim R_{\text{LC}} \ll R$ and, hence, even if the overall current sheet length is optically thick to pair production, the spacing between current sheets is thin, such that photons emitted in one stripe may be absorbed inside another.

Let us now examine what the implications of Klein–Nishina reconnection are on the pre-shocked pulsar wind of gamma-ray binaries. We discuss first the potential effect of pair production on the wind. The cold magnetization at the light cylinder is

$$ \sigma_{c,0}^{(w)} = \frac{B_{\text{LC}}^2}{4\pi\kappa n_{\text{GJ}} m_e c^2} = \frac{eP B_{\text{LC}}}{4\pi\kappa m_e c} = 1 \times 10^5 \left( \frac{B_{\text{LC}}}{9 \times 10^3 \text{ G}} \right) \left( \frac{\kappa}{10^4} \right)^{-1}, \qquad (76) $$

where $n_{\text{GJ}} \equiv B_{\text{LC}}/ecP$ is the Goldreich–Julian number density and $\kappa = 10^4$ is a fiducial pair multiplicity (e.g. Timokhin & Harding 2015, 2019). This magnetization is ‘frozen-in’ at the light cylinder in the sense that, because $B_w \propto 1/R$ and $n_0 \propto 1/R^2$ beyond the light cylinder, $\sigma_{c,0}^{(w)}$ remains constant (modulo pair production) in the unreconnected plasma of the wind throughout its expansion. The pair yield control parameter in the pre-shocked wind is then $\gamma_{\text{pp,max}}/4\sigma_{c,0}^{(w)} \simeq 30\gamma_{\text{KN,w}}/4\sigma_{c,0}^{(w)} \sim 0.7$, which is in the regime of order-unity pair yield. Furthermore, if the pulsar magnetosphere fails to launch the pulsar wind with $\kappa \sim 10^4$, but instead with a lower multiplicity, the magnetization $\sigma_{c,0}^{(w)}$ increases, leading to potentially copious in situ pair production in the expanding wind. Using our results from Section 5, we infer a critical self-regulated magnetization of $\sigma_c^{(\gamma\gamma)} \simeq \gamma_{\text{pp,max}} \simeq 30\gamma_{\text{KN,w}} \sim 4 \times 10^5$. If the plasma is injected with a higher magnetization than this (e.g. by virtue of an underdense multiplicity), Klein–Nishina pair production fills in the plasma deficit, pulling up the multiplicity toward $\kappa = eP B_{\text{LC}}/4\pi m_e c \sigma_c^{(\gamma\gamma)} \sim 3 \times 10^3$. This critical multiplicity depends solely on the pulsar properties and those of the ambient radiation field.

Let us close by examining prospects for observing signatures of Klein–Nishina reconnection in gamma-ray binaries. The temporal observables uncovered in Section 4.5 likely manifest themselves most clearly during transient flares. Of the known gamma-ray binaries, PSR B1259–63 is known to flare once per 3.4-yr orbit (e.g. Abdo et al. 2011; Tam et al. 2011, 2018; Caliandro et al. 2015; Chang et al. 2021). However, this binary is far less compact

than others, with $d \sim 0.9$ au $\gg 0.1$ au, even at periastron. Moreover, the flares occur significantly later in the orbit than the time of periastron. Such wide separations dilute the radiation field from the stellar companion, tending to move reconnection out of the Klein–Nishina regime. The binaries LS 5039, LS I+61°303, and 1FGL J1018.6–5856 are all much more compact ($d \sim 0.1$ au at periastron; Dubus 2013), but, instead of sudden flares, GeV gamma-ray observations reveal smooth modulations on each binary’s orbital period (for LS 5039, Abdo et al. 2009c; for LS I+61°303, Abdo et al. 2009b; for 1FGL J1018.6–5856, Fermi LAT Collaboration 2012). Such modulations likely probe quasistatic changes to the unshocked pulsar wind, including its orbit-dependent illumination by the companion, rather than variability intrinsic to reconnection (a situation anticipated in the concluding remarks of Section 4.5). In order to interpret these modulations, one needs to go beyond the flaring-centric treatment of observable signatures adopted in this work and conduct explicit global modelling of the steady-state wind.

### 6.4.3 Synopsis of gamma-ray binary discussion

To summarize, we have shown that Klein–Nishina reconnection likely takes place in short-period gamma-ray binaries if the compact object in the binary is a pulsar. In that case, reconnection occurs in the pre-shocked striped pulsar wind and is immersed in a bright seed photon bath supplied by the companion star. As a result, the reconnecting stripes produce a minimum pair-plasma density, setting a multiplicity/density floor in the wind (even if it is launched underdense from the pulsar). Observational support for this scenario is somewhat limited, as the few observed bright flares occur in configurations where the stellar companion is too far removed to supply a sufficiently dense radiation bath for Klein–Nishina reconnection. More compact systems, on the other hand, do not exhibit rapid flares, but rather smooth modulations to their GeV signal over the entire orbit. Detailed predictions of this signal require global modelling, which we leave to a future study.

## 6.5 Overall discussion summary

We supply a graphic recapitulation of the discussion in this section in Table 4. In all four types of systems – FSRQs, (high/soft-states of) black hole accretion disc coronae, the M87* magnetosphere, and gamma-ray binaries – a strong case can be made that the basic Klein–Nishina reconnection scale hierarchy (50), $\gamma_{\text{cool}} < \gamma_{\text{KN}} < \gamma_{\text{rad, T}} < \gamma_{\text{max}}$, is realized. However, only for the FSRQs [particularly the IC(BLR) scenario] do suitable observations exist for comparing with the expected temporal signatures of Klein–Nishina reconnection. For these objects, we find that Fermi-LAT data are, on the whole, consistent with the picture of a reconnection-powered flaring spectral energy density stabilized by pair production (Section 4.5). For the other object classes, more detailed observational connections are inhibited by: current instrumental sensitivity for the FSRQ IC(HDR) scenario (see Section 6.1; the coming online of the CTA may soon alleviate this); intrinsically dim sources for the accretion disc coronae of BHXRBs (see Section 6.2; this difficulty will probably not be overcome); and modelling details for the M87* magnetosphere and gamma-ray binaries (see Sections 6.3 and 6.4; the most pressing needs here are to account for synchrotron cooling in the case of M87* and for the global geometry of gamma-ray binaries).
critical scale hierarchy (50), $\gamma_{\text{cool}} < \gamma_{\text{KN}} < \gamma_{\text{rad, T}} < \gamma_{\text{max}}$. Assuming $B_{\text{LC}} = 9 \times 10^3$ G, $P = 100$ ms, and that the reconnecting magnetic field strength, $B_w$, is $B_w = B_{\text{LC}}R_{\text{LC}}/R$, as well as identifying the length, $L$, of the reconnection layer with the local radius, $R$, in the pulsar wind, we have, by equation (38),

$$ \gamma_{\text{max},w} \equiv \frac{0.1eB_w R}{m_e c^2} \sim 3 \times 10^8, \qquad (72) $$

which is independent of $R$. Similarly, we have, by equation (48),

$$ \gamma_{\text{rad, T}}^{(w)} \equiv \sqrt{\frac{0.3eB_w}{4\sigma_{\text{T}}U_{\star}}} \sim 3 \times 10^6 \left(\frac{R}{60R_{\text{LC}}}\right)^{-1/2}. \qquad (73) $$

We then can use $\gamma_{\text{rad, T}}^2 = \gamma_{\text{cool}}\gamma_{\text{max}}$ [equation (49)], to estimate

$$ \gamma_{\text{cool},w} \sim 3 \times 10^4 \left( \frac{R}{60 R_{\text{LC}}} \right)^{-1}. \qquad (74) $$

Finally, using equation (51), we estimate the pair-production optical depth as

$$ \tau_{\gamma\gamma,w} \equiv \frac{3\gamma_{\text{KN},w}}{5\gamma_{\text{cool},w}} \sim 0.3 \left(\frac{R}{60R_{\text{LC}}}\right). \qquad (75) $$

Note that we have normalized $R$ to the critical radius where IC losses overtake synchrotron losses. The pulsar wind extends much farther than this, which pushes $\tau_{\gamma\gamma,w}$ above unity at the largest radii. The geometry of the striped wind is somewhat peculiar, however, in that the spacing between stripes is $\sim R_{\text{LC}} \ll R$ and, hence, even if the overall current sheet length is optically thick to pair production, the spacing between current sheets is thin, such that photons emitted in one stripe may be absorbed inside another.

Let us now examine what the implications of Klein–Nishina reconnection are on the pre-shocked pulsar wind of gamma-ray binaries. We discuss first the potential effect of pair production on the wind. The cold magnetization at the light cylinder is

$$ \sigma_{c,0}^{(w)} = \frac{B_{\text{LC}}^2}{4\pi\kappa n_{\text{GJ}} m_e c^2} = \frac{eP B_{\text{LC}}}{4\pi\kappa m_e c} = 1 \times 10^5 \left( \frac{B_{\text{LC}}}{9 \times 10^3 \text{ G}} \right) \left( \frac{\kappa}{10^4} \right)^{-1}, \qquad (76) $$

where $n_{\text{GJ}} \equiv B_{\text{LC}}/ecP$ is the Goldreich–Julian number density and $\kappa = 10^4$ is a fiducial pair multiplicity (e.g. Timokhin & Harding 2015, 2019). This magnetization is ‘frozen-in’ at the light cylinder in the sense that, because $B_w \propto 1/R$ and $n_0 \propto 1/R^2$ beyond the light cylinder, $\sigma_{c,0}^{(w)}$ remains constant (modulo pair production) in the unreconnected plasma of the wind throughout its expansion. The pair yield control parameter in the pre-shocked wind is then $\gamma_{\text{pp,max}}/4\sigma_{c,0}^{(w)} \simeq 30\gamma_{\text{KN,w}}/4\sigma_{c,0}^{(w)} \sim 0.7$, which is in the regime of order-unity pair yield. Furthermore, if the pulsar magnetosphere fails to launch the pulsar wind with $\kappa \sim 10^4$, but instead with a lower multiplicity, the magnetization $\sigma_{c,0}^{(w)}$ increases, leading to potentially copious *in situ* pair production in the expanding wind. Using our results from Section 5, we infer a critical self-regulated magnetization of $\sigma_c^{(\gamma\gamma)} \simeq \gamma_{\text{pp,max}} \simeq 30\gamma_{\text{KN,w}} \sim 4 \times 10^5$. If the plasma is injected with a higher magnetization than this (e.g. by virtue of an underdense multiplicity), Klein–Nishina pair production fills in the plasma deficit, pulling up the multiplicity toward $\kappa = eP B_{\text{LC}}/4\pi m_e c \sigma_c^{(\gamma\gamma)} \sim 3 \times 10^3$. This critical multiplicity depends solely on the pulsar properties and those of the ambient radiation field.

Let us close by examining prospects for observing signatures of Klein–Nishina reconnection in gamma-ray binaries. The temporal observables uncovered in Section 4.5 likely manifest themselves most clearly during transient flares. Of the known gamma-ray binaries, PSR B1259–63 is known to flare once per 3.4-yr orbit (e.g. Abdo et al. 2011; Tam et al. 2011, 2018; Caliandro et al. 2015; Chang et al. 2021). However, this binary is far less compact

than others, with $d \sim 0.9$ au $\gg 0.1$ au, even at periastron. Moreover, the flares occur significantly later in the orbit than the time of periastron. Such wide separations dilute the radiation field from the stellar companion, tending to move reconnection out of the Klein–Nishina regime. The binaries LS 5039, LS I+61°303, and 1FGL J1018.6–5856 are all much more compact ($d \sim 0.1$ au at periastron; Dubus 2013), but, instead of sudden flares, GeV gamma-ray observations reveal smooth modulations on each binary’s orbital period (for LS 5039, Abdo et al. 2009c; for LS I+61°303, Abdo et al. 2009b; for 1FGL J1018.6–5856, Fermi LAT Collaboration 2012). Such modulations likely probe quasistatic changes to the unshocked pulsar wind, including its orbit-dependent illumination by the companion, rather than variability intrinsic to reconnection (a situation anticipated in the concluding remarks of Section 4.5). In order to interpret these modulations, one needs to go beyond the flaring-centric treatment of observable signatures adopted in this work and conduct explicit global modelling of the steady-state wind.

### 6.4.3 Synopsis of gamma-ray binary discussion

To summarize, we have shown that Klein–Nishina reconnection likely takes place in short-period gamma-ray binaries if the compact object in the binary is a pulsar. In that case, reconnection occurs in the pre-shocked striped pulsar wind and is immersed in a bright seed photon bath supplied by the companion star. As a result, the reconnecting stripes produce a minimum pair-plasma density, setting a multiplicity/density floor in the wind (even if it is launched underdense from the pulsar). Observational support for this scenario is somewhat limited, as the few observed bright flares occur in configurations where the stellar companion is too far removed to supply a sufficiently dense radiation bath for Klein–Nishina reconnection. More compact systems, on the other hand, do not exhibit rapid flares, but rather smooth modulations to their GeV signal over the entire orbit. Detailed predictions of this signal require global modelling, which we leave to a future study.

## 6.5 Overall discussion summary

We supply a graphic recapitulation of the discussion in this section in Table 4. In all four types of systems – FSRQs, (high/soft-states of) black hole accretion disc coronae, the M87* magnetosphere, and gamma-ray binaries – a strong case can be made that the basic Klein–Nishina reconnection scale hierarchy (50), $\gamma_{\text{cool}} < \gamma_{\text{KN}} < \gamma_{\text{rad, T}} < \gamma_{\text{max}}$, is realized. However, only for the FSRQs [particularly the IC(BLR) scenario] do suitable observations exist for comparing with the expected temporal signatures of Klein–Nishina reconnection. For these objects, we find that Fermi-LAT data are, on the whole, consistent with the picture of a reconnection-powered flaring spectral energy density stabilized by pair production (Section 4.5). For the other object classes, more detailed observational connections are inhibited by: current instrumental sensitivity for the FSRQ IC(HDR) scenario (see Section 6.1; the coming online of the CTA may soon alleviate this); intrinsically dim sources for the accretion disc coronae of BHXRBs (see Section 6.2; this difficulty will probably not be overcome); and modelling details for the M87* magnetosphere and gamma-ray binaries (see Sections 6.3 and 6.4; the most pressing needs here are to account for synchrotron cooling in the case of M87* and for the global geometry of gamma-ray binaries).
## 7 CONCLUSIONS

In this work, we construct (Section 3) a numerical model of Klein–Nishina reconnection. We consider a reconnection system immersed in a background radiation bath that is static, homogeneous, isotropic, and monochromatic, characterized solely by its total energy density $U_{ph}$ and individual photon energy $\epsilon_{ph}$ [equation (4)]. We restrict to a 2D reconnection setup, a consequence of the large box sizes needed to maintain an adequate separation between the many scales in the problem, several of which are introduced by the QED physics. In particular, we need to satisfy the critical hierarchy of Lorentz factor energy scales (50), $\gamma_{\text{cool}} \ll \gamma_{\text{KN}} \ll \gamma_{\text{rad, T}} \ll \gamma_{\text{max}}$. This hierarchy, which is apparently satisfied in a variety of astrophysical systems (Section 6; Mehlhaff et al. 2021), permits: (1) efficient IC radiative losses on time-scales much shorter than the light crossing time, $L/c$; (2) copious particle acceleration above the energy, $\sim \gamma_{\text{KN}}$, where particles emit photons above pair-production threshold; and (3) a high optical depth, $\tau_{\gamma\gamma} \sim \gamma_{\text{KN}}/\gamma_{\text{cool}} \gg 1$, such that nearly all of the above-threshold radiation is absorbed inside the system.

To model Klein–Nishina effects, we use a QED-enabled version of the ZELTRON PIC code, the development of which is presented in Section 2. We employ two main numerical optimizations to speed up the Monte Carlo IC scattering and pair production (Section 2.2.3): the null collision method and an approximation to the IC cross-section. Of these two, the former is widely used, while, prior to this work, the latter has, to our knowledge, only been implemented in GRZELTRON (see the supporting information of Crinquand et al. 2020), and we supply here the first exhaustive description.

We deploy the capabilities developed in Section 2 to run simulations of Klein–Nishina reconnection, which we present in Sections 4 and 5. These include various control runs – one non-radiative, one Thomson-radiative, and one Klein–Nishina-radiative with pair production artificially suppressed – to elicit the unique properties induced by the novel QED physics. We draw the following main conclusions for Klein–Nishina reconnection:

(i) From Sections 4.1 and 4.2: The added Klein–Nishina and pair-production physics does not substantially change the large-scale qualitative aspects of reconnection: the reconnection rate and hierarchical plasmoid chain remain essentially unchanged.

(ii) From Section 4.3: Reconnection-powered NTPA remains efficient but is somewhat inhibited by radiative cooling. On time-average, the reconnection-energized plasma is colder, and the particle energy distribution's extended non-thermal tail steeper, compared to the non-radiative regime, but the plasma is not as cold, nor the distribution as steep, as in the Thomson radiative regime.

(iii) From Section 4.4: The time-averaged IC emission spectrum is steeper than in the non-radiative regime but exhibits little change in slope from the Thomson-cooled case. This owes to Klein–Nishina effects suppressing the IC cross-section, and, hence, largely can-celling the hardening of the underlying particle energy distribution with respect to the Thomson case (cf. Moderski et al. 2005).

(iv) From Section 4.5: The intrinsic IC emission peaks above pair-production threshold, meaning that most of the initially radiated energy is reprocessed (by radiative cooling and pair production) to longer wavelengths before escaping the system. This leads, despite similarities in time-averaged spectra, to profound differences in the time-dependent signatures of reconnection between the Klein–Nishina and Thomson radiative regimes. In the latter case, the observed spectrum exhibits a tight correlation between total luminosity and prominence of the high-energy non-thermal tail (i.e. harder-than-brighter). In the Klein–Nishina case, however, the spectral shape becomes virtually independent of total luminosity, including in

the decaying phase after reconnection-powered particle energization has ceased.

(v) From Section 4.5: When particle energization is shut off, Klein–Nishina IC cooling and pair production conspire to produce a universal spectral shape, with the particle energy distribution scaling as $dN/d\gamma \propto \gamma^{-2}$ for $\gamma < \gamma_{\text{KN}}$ and IC emission spectrum as $\epsilon F_{\text{IC}}(\epsilon) \propto \epsilon^{1/2}$ for $\epsilon < \gamma_{\text{KN}}m_e c^2$.

(vi) From Section 4.6: The upstream plasma enthalpy – and, hence, the upstream hot magnetization – is not sufficiently loaded by newborn pairs to trigger the limit cycles sketched by Mehlhaff et al. (2021), according to the pair-loading efficiency requirement derived by those authors. This conclusion, however, applies only when the newborn upstream pairs are hot and tenuous, contributing significantly to the pressure of the reconnection inflow plasma but not to its number density. In regimes (suggested to exist by our findings concerning the pair yield below) where pair production loads the upstream plasma number density – and, thus, the cold magnetization in addition to the hot magnetization – the possibility of limit cycles remains an open issue.

(vii) From Section 5: The pair yield (per reconnection-processed lepton) follows an empirical exponential decay law, equation (53), in the parameter, $\gamma_{\text{pp, max}}/4\sigma_{c,0} \simeq 30\gamma_{\text{KN}}/4\sigma_{c,0}$. When $\gamma_{\text{pp, max}} \lesssim 4\sigma_{c,0}$, reconnection-powered NTPA cuts off at $\gamma \sim \gamma_X \equiv 4\sigma_{c,0}$, far beyond the energy, $\gamma_{\text{pp, max}}$, where typical Comptonized photons lie at peak pair-production cross-section with the seed photons, enabling an order-unity pair yield. We speculate that the exponential formula breaks down in the limit, $4\sigma_{c,0} \gg \gamma_{\text{pp, max}}$ (not probed by our simulations), where reconnection would instead produce copious pairs. In such a case, the population of newborn pairs would no longer be hot and tenuous as in Section 4.6, but rather hot and abundant.

(viii) From Section 6: The Klein–Nishina reconnection scale hierarchy (50) is likely satisfied in at least four classes of astrophysical systems: FSRQs, high/soft states of BHXRBs, the magnetosphere of M87*, and gamma-ray binaries. Of note are FSRQs, where the often-observed spectral stability during GeV flaring states could be due to Klein–Nishina effects as discussed (and summarized above) in Section 4.5. For the other object classes, further modelling or instrumental development will enable firmer observational connections. A more detailed summary of Section 6 can be found in Section 6.5 and Table 4.

A few broad implications of Section 4.5's results merit additional discussion. The spectral stability properties in that section stem from the similar shape of the Klein–Nishina reconnection IC spectrum to that produced when an initial population of particles is allowed to passively cool through IC emission and pair production. This introduces a fundamental degree of degeneracy in efforts to leverage astrophysical observations to learn about plasma physical particle acceleration processes. Namely, any Klein–Nishina-coupled particle accelerator with an IC spectrum resembling that to which the radiative physics inevitably relaxes (i.e. after the accelerator is turned off) may yield the same spectral stability as reconnection. However, this degeneracy is, thankfully, incomplete, because acceleration processes with softer intrinsic spectra would instead likely yield an asymmetry between the rising side of a gamma-ray flare (where the steep particle acceleration spectrum would be visible) and the decaying phase (where the shallower Klein–Nishina-decaying spectrum would probably dominate).

In the face of such degeneracy, observations of temporally resolved spectra provide precious information. Already in this study, such spectra provide crucial and obvious distinctions between Thomson and Klein–Nishina reconnection – two regimes whose time-averaged
## 7 CONCLUSIONS

In this work, we construct (Section 3) a numerical model of Klein–Nishina reconnection. We consider a reconnection system immersed in a background radiation bath that is static, homogeneous, isotropic, and monochromatic, characterized solely by its total energy density $U_{ph}$ and individual photon energy $\epsilon_{ph}$ [equation (4)]. We restrict to a 2D reconnection setup, a consequence of the large box sizes needed to maintain an adequate separation between the many scales in the problem, several of which are introduced by the QED physics. In particular, we need to satisfy the critical hierarchy of Lorentz factor energy scales (50), $\gamma_{\text{cool}} \ll \gamma_{\text{KN}} \ll \gamma_{\text{rad, T}} \ll \gamma_{\text{max}}$. This hierarchy, which is apparently satisfied in a variety of astrophysical systems (Section 6; Mehlhaff et al. 2021), permits: (1) efficient IC radiative losses on time-scales much shorter than the light crossing time, $L/c$; (2) copious particle acceleration above the energy, $\sim\gamma_{\text{KN}}$, where particles emit photons above pair-production threshold; and (3) a high optical depth, $\tau_{\gamma\gamma} \sim \gamma_{\text{KN}}/\gamma_{\text{cool}} \gg 1$, such that nearly all of the above-threshold radiation is absorbed inside the system.

To model Klein–Nishina effects, we use a QED-enabled version of the ZELTRON PIC code, the development of which is presented in Section 2. We employ two main numerical optimizations to speed up the Monte Carlo IC scattering and pair production (Section 2.2.3): the null collision method and an approximation to the IC cross-section. Of these two, the former is widely used, while, prior to this work, the latter has, to our knowledge, only been implemented in GRZELTRON (see the supporting information of Crinquand et al. 2020), and we supply here the first exhaustive description.

We deploy the capabilities developed in Section 2 to run simulations of Klein–Nishina reconnection, which we present in Sections 4 and 5. These include various control runs – one non-radiative, one Thomson-radiative, and one Klein–Nishina-radiative with pair production artificially suppressed – to elicit the unique properties induced by the novel QED physics. We draw the following main conclusions for Klein–Nishina reconnection:

(i) From Sections 4.1 and 4.2: The added Klein–Nishina and pair-production physics does not substantially change the large-scale qualitative aspects of reconnection: the reconnection rate and hierarchical plasmoid chain remain essentially unchanged.

(ii) From Section 4.3: Reconnection-powered NTPA remains efficient but is somewhat inhibited by radiative cooling. On time-average, the reconnection-energized plasma is colder, and the particle energy distribution's extended non-thermal tail steeper, compared to the non-radiative regime, but the plasma is not as cold, nor the distribution as steep, as in the Thomson radiative regime.

(iii) From Section 4.4: The time-averaged IC emission spectrum is steeper than in the non-radiative regime but exhibits little change in slope from the Thomson-cooled case. This owes to Klein–Nishina effects suppressing the IC cross-section, and, hence, largely can-celling the hardening of the underlying particle energy distribution with respect to the Thomson case (cf. Moderski et al. 2005).

(iv) From Section 4.5: The intrinsic IC emission peaks above pair-production threshold, meaning that most of the initially radiated energy is reprocessed (by radiative cooling and pair production) to longer wavelengths before escaping the system. This leads, despite similarities in time-averaged spectra, to profound differences in the time-dependent signatures of reconnection between the Klein–Nishina and Thomson radiative regimes. In the latter case, the observed spectrum exhibits a tight correlation between total luminosity and prominence of the high-energy non-thermal tail (i.e. harder-than-brighter). In the Klein–Nishina case, however, the spectral shape becomes virtually independent of total luminosity, including in

the decaying phase after reconnection-powered particle energization has ceased.

(v) From Section 4.5: When particle energization is shut off, Klein–Nishina IC cooling and pair production conspire to produce a universal spectral shape, with the particle energy distribution scaling as $dN/d\gamma \propto \gamma^{-2}$ for $\gamma < \gamma_{\text{KN}}$ and IC emission spectrum as $\epsilon F_{\text{IC}}(\epsilon) \propto \epsilon^{1/2}$ for $\epsilon < \gamma_{\text{KN}}m_e c^2$.

(vi) From Section 4.6: The upstream plasma enthalpy – and, hence, the upstream hot magnetization – is not sufficiently loaded by newborn pairs to trigger the limit cycles sketched by Mehlhaff et al. (2021), according to the pair-loading efficiency requirement derived by those authors. This conclusion, however, applies only when the newborn upstream pairs are hot and tenuous, contributing significantly to the pressure of the reconnection inflow plasma but not to its number density. In regimes (suggested to exist by our findings concerning the pair yield below) where pair production loads the upstream plasma number density – and, thus, the cold magnetization in addition to the hot magnetization – the possibility of limit cycles remains an open issue.

(vii) From Section 5: The pair yield (per reconnection-processed lepton) follows an empirical exponential decay law, equation (53), in the parameter, $\gamma_{\text{pp, max}}/4\sigma_{c,0} \simeq 30\gamma_{\text{KN}}/4\sigma_{c,0}$. When $\gamma_{\text{pp, max}} \lesssim 4\sigma_{c,0}$, reconnection-powered NTPA cuts off at $\gamma \sim \gamma_X \equiv 4\sigma_{c,0}$, far beyond the energy, $\gamma_{\text{pp, max}}$, where typical Comptonized photons lie at peak pair-production cross-section with the seed photons, enabling an order-unity pair yield. We speculate that the exponential formula breaks down in the limit, $4\sigma_{c,0} \gg \gamma_{\text{pp, max}}$ (not probed by our simulations), where reconnection would instead produce copious pairs. In such a case, the population of newborn pairs would no longer be hot and tenuous as in Section 4.6, but rather hot and abundant.

(viii) From Section 6: The Klein–Nishina reconnection scale hierarchy (50) is likely satisfied in at least four classes of astrophysical systems: FSRQs, high/soft states of BHXRBs, the magnetosphere of M87*, and gamma-ray binaries. Of note are FSRQs, where the often-observed spectral stability during GeV flaring states could be due to Klein–Nishina effects as discussed (and summarized above) in Section 4.5. For the other object classes, further modelling or instrumental development will enable firmer observational connections. A more detailed summary of Section 6 can be found in Section 6.5 and Table 4.

A few broad implications of Section 4.5's results merit additional discussion. The spectral stability properties in that section stem from the similar shape of the Klein–Nishina reconnection IC spectrum to that produced when an initial population of particles is allowed to passively cool through IC emission and pair production. This introduces a fundamental degree of degeneracy in efforts to leverage astrophysical observations to learn about plasma physical particle acceleration processes. Namely, any Klein–Nishina-coupled particle accelerator with an IC spectrum resembling that to which the radiative physics inevitably relaxes (i.e. after the accelerator is turned off) may yield the same spectral stability as reconnection. However, this degeneracy is, thankfully, incomplete, because acceleration processes with softer intrinsic spectra would instead likely yield an asymmetry between the rising side of a gamma-ray flare (where the steep particle acceleration spectrum would be visible) and the decaying phase (where the shallower Klein–Nishina-decaying spectrum would probably dominate).

In the face of such degeneracy, observations of temporally resolved spectra provide precious information. Already in this study, such spectra provide crucial and obvious distinctions between Thomson and Klein–Nishina reconnection – two regimes whose time-averaged
spectra have practically identical photon indices – and, in the case of FSRQs (Section 6), supply compelling evidence for the operation of the latter over the former. In the future, time-resolved spectra may be necessary to distinguish among different possible Klein–Nishina-coupled particle acceleration mechanisms. Hence, this work strongly motivates collection of time-resolved (versus simply time-averaged) spectra whenever possible, as well as new instruments for which such observations are more frequently feasible.

This study also lays the groundwork for a broad range of further theoretical exploration. As detailed in Section 6, future efforts can target the M87* magnetosphere by explicitly treating synchrotron radiative cooling; firmer connections to gamma-ray binaries can be made by considering global properties of irradiated striped pulsar winds. Additional modelling may also verify the potential existence, outlined here (Section 5), of a Klein–Nishina reconnection regime featuring copious electron–positron pair production, which would have important implications for systems with a high fiducial cold magnetization, $\sigma_{c,0} \gg \gamma_{pp, max}/4$. However, even with a more modest order-unity pair yield (realized, in our simulations, when $\gamma_{pp, max}/4\sigma_{c,0} \lesssim 0.25$), Klein–Nishina reconnection could still function as an important source of *in situ* antimatter, taking, for example, an initially electron–proton plasma and injecting a prominent positron count. Such potential modifications to the composition of the reconnection plasma environment can be rigorously probed by applying the Klein–Nishina radiative physics studied here to the case of an initial electron–ion plasma.

More broadly, the theoretical framework established by this work (and the preceding article, Mehlhaff et al. 2021), especially the systematic recasting of Klein–Nishina physics as dimensionless energy scales, provides a paradigm for coupling the same physics to other particle acceleration processes (e.g. turbulence and shocks). That is, apart from factors of the reconnection rate, $\beta_{\text{rec}} = 0.1$, entering into $\gamma_{\text{max}}$ and $\gamma_{\text{rad, T}}$, the Klein–Nishina scale hierarchy (50) also applies to other plasma physical processes, indicating when they, too, may couple to Klein–Nishina effects. The physical and observable consequences of this coupling, however, may be very different than in the context of reconnection as uncovered here and can only be elucidated via future dedicated simulations. This work thus represents only the first step toward a much larger suite of numerical investigations into Klein–Nishina-coupled particle acceleration, and, even more broadly, constitutes one of a growing number of theoretical and technological stepping stones (e.g. Hakobyan et al. 2019, 2023b; Schoeffler et al. 2019, 2023; Crinquand et al. 2020, 2021; Cruz et al. 2021; Crinquand et al. 2022; Chen et al. 2023; Groselj et al. 2023; Xie et al. 2023) toward a richer understanding of QED-coupled plasma physics in high-energy astrophysical environments.

## ACKNOWLEDGEMENTS

The authors gratefully acknowledge Guillaume Dubus, Gilles Henri, and Hayk Hakobyan for stimulating discussions. The authors also extend their gratitude to the anonymous referee for a thorough review that greatly improved the manuscript. This project has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement no. 863412). This work was also supported by NASA and the NSF, grant nos NASA ATP NNX17AK57G, NASA ATP 80NSSC20K0545, NASA ATP 80NSSC22K0828, NSF AST-1806084, and NSF AST-1903335. This research has made use of the TeVCat online source catalog (http://tevcat.uchicago.edu).

## DATA AVAILABILITY

The simulation data underlying this article were generated at the XSEDE/TACC STAMPEDE2 supercomputer and are archived at the TACC/Ranch storage facility. As long as the data remain in the archive, they will be shared upon reasonable request to the corresponding author. The code used to generate the data, ZELTRON, is available at https://ipag.osug.fr/~ceruttbe/Zeltron/. The public version of the code does not yet contain the QED module underlying this work; a QED-enabled version will be shared upon reasonable request to the corresponding author.

## REFERENCES

* Abdo A. A. et al., 2009a, *Science*, 325, 840

* Abdo A. A. et al., 2009b, *ApJ*, 701, L123

* Abdo A. A. et al., 2009c, *ApJ*, 706, L56

* Abdo A. A. et al., 2011, *ApJ*, 736, L11

* Abdollahi S. et al., 2020, *ApJS*, 247, 33

* Abdollahi S. et al., 2023, *ApJS*, 265, 31

* Abeysekara A. U. et al., 2015, *ApJ*, 815, L22

* Abramowski A. et al., 2012, *ApJ*, 746, 151

* Acciari V. A. et al., 2009, *Science*, 325, 444

* Acciari V. A. et al., 2010, *ApJ*, 716, 819

* Ackermann M. et al., 2015, *ApJ*, 810, 14

* Ackermann M. et al., 2016, *ApJ*, 824, L20

* Aharonian F. et al., 2006, *Science*, 314, 1424

* Aharonian F. et al., 2007, *ApJ*, 664, L71

* Aharonian F. et al., 2009, *A&A*, 502, 749

* Ahnen M. L. et al., 2016, *A&A*, 595, A98

* Ahnen M. L. et al., 2017, *A&A*, 603, A29

* Ajello M. et al., 2020, *ApJ*, 892, 105

* Aleksić J. et al., 2011, *ApJ*, 730, L8

* Aleksić J. et al., 2014, *A&A*, 569, A46

* Aliu E. et al., 2012, *ApJ*, 746, 141

* Atwood W. B. et al., 2009, *ApJ*, 697, 1071

* Avara M. J., McKinney J. C., Reynolds C. S., 2016, *MNRAS*, 462, 636

* Ball L., Dodd J., 2001, *Publ. Astron. Soc. Aust.*, 18, 98

* Ball L., Kirk J. G., 2000, *Astropart. Phys.*, 12, 335

* Ball D., Özel F., Psaltis D., Chan C.-k., 2016, *ApJ*, 826, 77

* Ball D., Sironi L., Özel F., 2018, *ApJ*, 862, 80

* Begelman M. C., Sikora M., 1987, *ApJ*, 322, 650

* Beilicke M., VERITAS Collaboration, 2012, in Aharonian F. A., Hofmann W., Rieger F. M., eds, AIP Conf. Ser. Vol. 1505, High Energy Gamma-Ray Astronomy: 5th International Meeting on High Energy Gamma-Ray Astronomy. Am. Inst. Phys., New York, p. 586

* Beloborodov A. M., 2017, *ApJ*, 850, 141

* Bird G. A., 1989, Prog. Astronaut. Aeronaut., 117, 211

* Birdsall C. K., 1991, *IEEE Trans. Plasma Sci.*, 19, 65

* Birdsall C. K., Langdon A. B., 1991, Plasma Physics via Computer Simulation. Adam Hilger, Bristol, England

* Blackman E. G., Field G. B., 1994, *Phys. Rev. Lett.*, 72, 494

* Blandford R., Meier D., Readhead A., 2019, *ARA&A*, 57, 467

* Bloom S. D., Marscher A. P., 1996, *ApJ*, 461, 657

* Blumenthal G. R., Gould R. J., 1970, *Rev. Mod. Phys.*, 42, 237

* Boeuf J. P., Marode E., 1982, *J. Phys. D Appl. Phys.*, 15, 2169

* Bogovalov S. V., 1999, *A&A*, 349, 1017

* Caliandro G. A., Cheung C. C., Li J., Scargle J. D., Torres D. F., Wood K. S., Chernyakova M., 2015, *ApJ*, 811, 68

* Cerutti B., Philippov A. A., 2017, *A&A*, 607, A134

* Cerutti B., Werner G., 2019, Astrophysics Source Code Library, record ascl:1911.012

* Cerutti B., Dubus G., Henri G., 2008, *A&A*, 488, 37

* Cerutti B., Uzdensky D. A., Begelman M. C., 2012a, *ApJ*, 746, 148

* Cerutti B., Werner G. R., Uzdensky D. A., Begelman M. C., 2012b, *ApJ*, 754, L33
spectra have practically identical photon indices – and, in the case of FSRQs (Section 6), supply compelling evidence for the operation of the latter over the former. In the future, time-resolved spectra may be necessary to distinguish among different possible Klein–Nishina-coupled particle acceleration mechanisms. Hence, this work strongly motivates collection of time-resolved (versus simply time-averaged) spectra whenever possible, as well as new instruments for which such observations are more frequently feasible.

This study also lays the groundwork for a broad range of further theoretical exploration. As detailed in Section 6, future efforts can target the M87* magnetosphere by explicitly treating synchrotron radiative cooling; firmer connections to gamma-ray binaries can be made by considering global properties of irradiated striped pulsar winds. Additional modelling may also verify the potential existence, outlined here (Section 5), of a Klein–Nishina reconnection regime featuring copious electron–positron pair production, which would have important implications for systems with a high fiducial cold magnetization, $\sigma_{c,0} \gg \gamma_{pp, \text{max}}/4$. However, even with a more modest order-unity pair yield (realized, in our simulations, when $\gamma_{pp, \text{max}}/4\sigma_{c,0} \lesssim 0.25$), Klein–Nishina reconnection could still function as an important source of *in situ* antimatter, taking, for example, an initially electron–proton plasma and injecting a prominent positron count. Such potential modifications to the composition of the reconnection plasma environment can be rigorously probed by applying the Klein–Nishina radiative physics studied here to the case of an initial electron–ion plasma.

More broadly, the theoretical framework established by this work (and the preceding article, Mehlhaff et al. 2021), especially the systematic recasting of Klein–Nishina physics as dimensionless energy scales, provides a paradigm for coupling the same physics to other particle acceleration processes (e.g. turbulence and shocks). That is, apart from factors of the reconnection rate, $\beta_{\text{rec}} = 0.1$, entering into $\gamma_{\text{max}}$ and $\gamma_{\text{rad, T}}$, the Klein–Nishina scale hierarchy (50) also applies to other plasma physical processes, indicating when they, too, may couple to Klein–Nishina effects. The physical and observable consequences of this coupling, however, may be very different than in the context of reconnection as uncovered here and can only be elucidated via future dedicated simulations. This work thus represents only the first step toward a much larger suite of numerical investigations into Klein–Nishina-coupled particle acceleration, and, even more broadly, constitutes one of a growing number of theoretical and technological stepping stones (e.g. Hakobyan et al. 2019, 2023b; Schoeffler et al. 2019, 2023; Crinquand et al. 2020, 2021; Cruz et al. 2021; Crinquand et al. 2022; Chen et al. 2023; Groselj et al. 2023; Xie et al. 2023) toward a richer understanding of QED-coupled plasma physics in high-energy astrophysical environments.

## ACKNOWLEDGEMENTS

The authors gratefully acknowledge Guillaume Dubus, Gilles Henri, and Hayk Hakobyan for stimulating discussions. The authors also extend their gratitude to the anonymous referee for a thorough review that greatly improved the manuscript. This project has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement no. 863412). This work was also supported by NASA and the NSF, grant nos NASA ATP NNX17AK57G, NASA ATP 80NSSC20K0545, NASA ATP 80NSSC22K0828, NSF AST-1806084, and NSF AST-1903335. This research has made use of the TeVCat online source catalog (http://tevcat.uchicago.edu).

## DATA AVAILABILITY

The simulation data underlying this article were generated at the XSEDE/TACC STAMPEDE2 supercomputer and are archived at the TACC/Ranch storage facility. As long as the data remain in the archive, they will be shared upon reasonable request to the corresponding author. The code used to generate the data, ZELTRON, is available at https://ipag.osug.fr/~ceruttbe/Zeltron/. The public version of the code does not yet contain the QED module underlying this work; a QED-enabled version will be shared upon reasonable request to the corresponding author.

## REFERENCES

Abdo A. A. et al., 2009a, *Science*, 325, 840
Abdo A. A. et al., 2009b, *ApJ*, 701, L123
Abdo A. A. et al., 2009c, *ApJ*, 706, L56
Abdo A. A. et al., 2011, *ApJ*, 736, L11
Abdollahi S. et al., 2020, *ApJS*, 247, 33
Abdollahi S. et al., 2023, *ApJS*, 265, 31
Abeysekara A. U. et al., 2015, *ApJ*, 815, L22
Abramowski A. et al., 2012, *ApJ*, 746, 151
Acciari V. A. et al., 2009, *Science*, 325, 444
Acciari V. A. et al., 2010, *ApJ*, 716, 819
Ackermann M. et al., 2015, *ApJ*, 810, 14
Ackermann M. et al., 2016, *ApJ*, 824, L20
Aharonian F. et al., 2006, *Science*, 314, 1424
Aharonian F. et al., 2007, *ApJ*, 664, L71
Aharonian F. et al., 2009, *A&A*, 502, 749
Ahnen M. L. et al., 2016, *A&A*, 595, A98
Ahnen M. L. et al., 2017, *A&A*, 603, A29
Ajello M. et al., 2020, *ApJ*, 892, 105
Aleksić J. et al., 2011, *ApJ*, 730, L8
Aleksić J. et al., 2014, *A&A*, 569, A46
Aliu E. et al., 2012, *ApJ*, 746, 141
Atwood W. B. et al., 2009, *ApJ*, 697, 1071
Avara M. J., McKinney J. C., Reynolds C. S., 2016, *MNRAS*, 462, 636
Ball L., Dodd J., 2001, *Publ. Astron. Soc. Aust.*, 18, 98
Ball L., Kirk J. G., 2000, *Astropart. Phys.*, 12, 335
Ball D., Özel F., Psaltis D., Chan C.-k., 2016, *ApJ*, 826, 77
Ball D., Sironi L., Özel F., 2018, *ApJ*, 862, 80
Begelman M. C., Sikora M., 1987, *ApJ*, 322, 650
Beilicke M., VERITAS Collaboration, 2012, in Aharonian F. A., Hofmann W., Rieger F. M., eds, AIP Conf. Ser. Vol. 1505, High Energy Gamma-Ray Astronomy: 5th International Meeting on High Energy Gamma-Ray Astronomy. Am. Inst. Phys., New York, p. 586
Beloborodov A. M., 2017, *ApJ*, 850, 141
Bird G. A., 1989, Prog. Astronaut. Aeronaut., 117, 211
Birdsall C. K., 1991, *IEEE Trans. Plasma Sci.*, 19, 65
Birdsall C. K., Langdon A. B., 1991, Plasma Physics via Computer Simulation. Adam Hilger, Bristol, England
Blackman E. G., Field G. B., 1994, *Phys. Rev. Lett.*, 72, 494
Blandford R., Meier D., Readhead A., 2019, *ARA&A*, 57, 467
Bloom S. D., Marscher A. P., 1996, *ApJ*, 461, 657
Blumenthal G. R., Gould R. J., 1970, *Rev. Mod. Phys.*, 42, 237
Boeuf J. P., Marode E., 1982, *J. Phys. D Appl. Phys.*, 15, 2169
Bogovalov S. V., 1999, *A&A*, 349, 1017
Caliandro G. A., Cheung C. C., Li J., Scargle J. D., Torres D. F., Wood K. S., Chernyakova M., 2015, *ApJ*, 811, 68
Cerutti B., Philippov A. A., 2017, *A&A*, 607, A134
Cerutti B., Werner G., 2019, Astrophysics Source Code Library, record ascl:1911.012
Cerutti B., Dubus G., Henri G., 2008, *A&A*, 488, 37
Cerutti B., Uzdensky D. A., Begelman M. C., 2012a, *ApJ*, 746, 148
Cerutti B., Werner G. R., Uzdensky D. A., Begelman M. C., 2012b, *ApJ*, 754, L33
Cerutti B., Werner G. R., Uzdensky D. A., Begelman M. C., 2013, ApJ, 770, 147

Cerutti B., Werner G. R., Uzdensky D. A., Begelman M. C., 2014a, Phys. Plasmas, 21, 056501

Cerutti B., Werner G. R., Uzdensky D. A., Begelman M. C., 2014b, ApJ, 782, 104

Cerutti B., Philippov A. A., Dubus G., 2020, A&A, 642, A204

Chang Z., Zhang S., Chen Y.-P., Ji L., Kong L.-D., Wang P.-J., 2021, Universe, 7, 472

Chashkina A., Bromberg O., Levinson A., 2021, MNRAS, 508, 1241

Chen A. Y., Uzdensky D., Dexter J., 2023, ApJ, 944, 173

Cherenkov Telescope Array Consortium et al., 2019, Science with the Cherenkov Telescope Array. World Scientific Publishing Co. Pte. Ltd., USA

Chernoglazov A., Hakobyan H., Philippov A. A., 2023, ApJ, 959, 122

Chernyakova M., Malyshev D., 2020, in Multifrequency Behaviour of High Energy Cosmic Sources – XIII. 3–8 June 2019. Palermo. p. 45, preprint (arXiv:2006.03615)

Chernyakova M. et al., 2019, A&A, 631, A177

Chow A., Rowan M. E., Sironi L., Davelaar J., Bodo G., Narayan R., 2023a, MNRAS, 524, 90

Chow A., Davelaar J., Rowan M. E., Sironi L., 2023b, ApJ, 951, L23

Christie I. M., Petropoulou M., Sironi L., Giannios D., 2019, MNRAS, 482, 65

Christie I. M., Petropoulou M., Sironi L., Giannios D., 2020, MNRAS, 492, 549

Comisso L., Sironi L., 2021, Phys. Rev. Lett., 127, 255102

Coroniti F. V., 1990, ApJ, 349, 538

Costamante L., Cutini S., Tosti G., Antolini E., Tramacere A., 2018, MNRAS, 477, 4749

Crinquand B., Cerutti B., Philippov A., Parfrey K., Dubus G., 2020, Phys. Rev. Lett., 124, 145101

Crinquand B., Cerutti B., Dubus G., Parfrey K., Philippov A., 2021, A&A, 650, A163

Crinquand B., Cerutti B., Dubus G., Parfrey K., Philippov A., 2022, Phys. Rev. Lett., 129, 205101

Cruz F., Grismayer T., Chen A. Y., Spitkovsky A., Silva L. O., 2021, ApJ, 919, L4

Del Gaudio F., Grismayer T., Fonseca R. A., Silva L. O., 2020, J. Plasma Phys., 86, 905860516

Derouillat J. et al., 2018, Comput. Phys. Commun., 222, 351

Di Matteo T., 1998, MNRAS, 299, L15

Dubus G., 2006, A&A, 456, 801

Dubus G., 2013, A&AR, 21, 64

Dubus G., Guillard N., Petrucci P.-O., Martin P., 2017, A&A, 608, A59

EHT MWL Science Working Group et al., 2021, ApJ, 911, L11

El Mellah I., Cerutti B., Crinquand B., Parfrey K., 2022, A&A, 663, A169

El Mellah I., Cerutti B., Crinquand B., 2023, A&A, 677, A67

Esirkepov T. Z., 2001, Comput. Phys. Commun., 135, 144

Event Horizon Telescope Collaboration, 2019, ApJ, 875, L6

Event Horizon Telescope Collaboration, 2021, ApJ, 910, L13

Fedeli L. et al., 2022, in SC22: International Conference for High Performance Computing, Networking, Storage and Analysis. IEEE Computer Society, Los Alamitos, CA, USA, p. 1

Fermi LAT Collaboration, 2012, Science, 335, 189

Fonseca R. A. et al., 2002, in Sloot P. M. A., Hoekstra A. G., Tan C. J. K., Dongarra J. J., eds, Computational Science – ICCS 2002. Springer, Berlin, Heidelberg, p. 342

Fossati G., Maraschi L., Celotti A., Comastri A., Ghisellini G., 1998, MNRAS, 299, 433

French O., Guo F., Zhang Q., Uzdensky D. A., 2023, ApJ, 948, 19

Galeev A. A., Rosner R., Vaiana G. S., 1979, ApJ, 229, 318

Galishnikova A., Philippov A., Quataert E., Bacchini F., Parfrey K., Ripperda B., 2023, Phys. Rev. Lett., 130, 115201

Gebhardt K., Adams J., Richstone D., Lauer T. R., Faber S. M., Gültekin K., Murphy J., Tremaine S., 2011, ApJ, 729, 119

Ghisellini G., 2011, in Aharonian F. A., Hofmann W., Rieger F. M., eds, AIP Conf. Proc. Vol. 1381, 25th Texas Symposium on Relativistic AstroPhysics (Texas 2010). Am. Inst. Phys., New York, p. 180

Ghisellini G., Righi C., Costamante L., Tavecchio F., 2017, MNRAS, 469, 255

Giannios D., 2013, MNRAS, 431, 355

Giannios D., Uzdensky D. A., 2019, MNRAS, 484, 1378

Giannios D., Uzdensky D. A., Begelman M. C., 2009, MNRAS, 395, L29

Giannios D., Uzdensky D. A., Begelman M. C., 2010, MNRAS, 402, 1649

Gierliński M., Zdziarski A. A., 2003, MNRAS, 343, L84

Goodbred M., Liu Y.-H., 2022, Phys. Rev. Lett., 129, 265101

Goodman J., Uzdensky D., 2008, ApJ, 688, 555

Gould R. J., Schréder G. P., 1967, Phys. Rev., 155, 1404

Greenwood A. D., Cartwright K. L., Luginsland J. W., Baca E. A., 2004, J. Comput. Phys., 201, 665

Grismayer T., Vranic M., Martins J. L., Fonseca R. A., Silva L. O., 2016, Phys. Plasmas, 23, 056706

Groselj D., Hakobyan H., Beloborodov A. M., Sironi L., Philippov A., 2023, preprint (arXiv:2301.11327)

Guo F., Li H., Daughton W., Liu Y.-H., 2014, Phys. Rev. Lett., 113, 155005

Guo F., Liu Y.-H., Daughton W., Li H., 2015, ApJ, 806, 167

Guo F., Li X., Daughton W., Kilian P., Li H., Liu Y.-H., Yan W., Ma D., 2019, ApJ, 879, L23

Guo F., Li X., Daughton W., Li H., Kilian P., Liu Y.-H., Zhang Q., Zhang H., 2021, ApJ, 919, 111

H. E. S. S. Collaboration, 2019, A&A, 627, A159

H. E. S. S. Collaboration, 2020, A&A, 633, A162

H. E. S. S. Collaboration, 2021, A&A, 648, A23

Hakobyan H., Philippov A., Spitkovsky A., 2019, ApJ, 877, 53

Hakobyan H., Petropoulou M., Spitkovsky A., Sironi L., 2021, ApJ, 912, 48

Hakobyan H., Spitkovsky A., Chernoglazov A., Philippov A., Groselj D., Mahlmann J., 2023a, PrincetonUniversity/tristan-mp-v2: v2.6, Zenodo. Available at:https://zenodo.org/records/7566725.

Hakobyan H., Ripperda B., Philippov A. A., 2023b, ApJ, 943, L29

Harrison M. G., Neukirch T., 2009, Phys. Rev. Lett., 102, 135003

Hayashida M. et al., 2012, ApJ, 754, 114

Hayashida M. et al., 2015, ApJ, 807, 79

Hockney R. W., Eastwood J. W., 1988, Computer Simulation using Particles. Hilger, Bristol

Hoshino M., Lyubarsky Y., 2012, Space Sci. Rev., 173, 521

Igumenshchev I. V., 2008, ApJ, 677, 317

Johnston S., Manchester R. N., Lyne A. G., Bailes M., Kaspi V. M., Qiao G., D'Amico N., 1992, ApJ, 387, L37

Jones F. C., 1968, Phys. Rev., 167, 1159

Khangulyan D., Aharonian F. A., Bogovalov S. V., Ribó M., 2012, ApJ, 752, L17

Kim J. Y. et al., 2018, A&A, 616, A188

Kimura S. S., Toma K., Noda H., Hada K., 2022, ApJ, 937, L34

Kirk J. G., Skjæræasen O., 2003, ApJ, 591, 366

Kirk J. G., Ball L., Skjæræasen O., 1999, Astropart. Phys., 10, 31

Kirk J. G., Lyubarsky Y., Petri J., 2009, in Becker W.ed., Astrophysics and Space Science Library, Vol. 357, Neutron Stars and Pulsars. Springer, Berlin Heidelberg, p. 421

Levinson A., Cerutti B., 2018, A&A, 616, A184

Li Y.-P., Yuan F., Wang Q. D., 2017, MNRAS, 468, 2552

Li X., Guo F., Li H., Li S., 2018, ApJ, 866, 4

Lin S. L., Bardsley J. N., 1977, J. Chem. Phys., 66, 435

Ly C., Walker R. C., Junor W., 2007, ApJ, 660, 200

Lyubarsky Y. E., 2003, MNRAS, 345, 153

Lyubarsky Y. E., 2005, MNRAS, 358, 113

Lyubarsky Y., Kirk J. G., 2001, ApJ, 547, 437

Lyutikov M., Uzdensky D., 2003, ApJ, 589, 893

MAGIC Collaboration, 2008, Science, 320, 1752

MAGIC Collaboration, 2018, A&A, 619, A159

MAGIC Collaboration, 2020, MNRAS, 492, 5354
Cerutti B., Werner G. R., Uzdensky D. A., Begelman M. C., 2013, ApJ, 770, 147

Cerutti B., Werner G. R., Uzdensky D. A., Begelman M. C., 2014a, Phys. Plasmas, 21, 056501

Cerutti B., Werner G. R., Uzdensky D. A., Begelman M. C., 2014b, ApJ, 782, 104

Cerutti B., Philippov A. A., Dubus G., 2020, A&A, 642, A204

Chang Z., Zhang S., Chen Y.-P., Ji L., Kong L.-D., Wang P.-J., 2021, Universe, 7, 472

Chashkina A., Bromberg O., Levinson A., 2021, MNRAS, 508, 1241

Chen A. Y., Uzdensky D., Dexter J., 2023, ApJ, 944, 173

Cherenkov Telescope Array Consortium et al., 2019, Science with the Cherenkov Telescope Array. World Scientific Publishing Co. Pte. Ltd., USA

Chernoglazov A., Hakobyan H., Philippov A. A., 2023, ApJ, 959, 122

Chernyakova M., Malyshev D., 2020, in Multifrequency Behaviour of High Energy Cosmic Sources – XIII. 3–8 June 2019. Palermo. p. 45, preprint (arXiv:2006.03615)

Chernyakova M. et al., 2019, A&A, 631, A177

Chow A., Rowan M. E., Sironi L., Davelaar J., Bodo G., Narayan R., 2023a, MNRAS, 524, 90

Chow A., Davelaar J., Rowan M. E., Sironi L., 2023b, ApJ, 951, L23

Christie I. M., Petropoulou M., Sironi L., Giannios D., 2019, MNRAS, 482, 65

Christie I. M., Petropoulou M., Sironi L., Giannios D., 2020, MNRAS, 492, 549

Comisso L., Sironi L., 2021, Phys. Rev. Lett., 127, 255102

Coroniti F. V., 1990, ApJ, 349, 538

Costamante L., Cutini S., Tosti G., Antolini E., Tramacere A., 2018, MNRAS, 477, 4749

Crinquand B., Cerutti B., Philippov A., Parfrey K., Dubus G., 2020, Phys. Rev. Lett., 124, 145101

Crinquand B., Cerutti B., Dubus G., Parfrey K., Philippov A., 2021, A&A, 650, A163

Crinquand B., Cerutti B., Dubus G., Parfrey K., Philippov A., 2022, Phys. Rev. Lett., 129, 205101

Cruz F., Grismayer T., Chen A. Y., Spitkovsky A., Silva L. O., 2021, ApJ, 919, L4

Del Gaudio F., Grismayer T., Fonseca R. A., Silva L. O., 2020, J. Plasma Phys., 86, 905860516

Derouillat J. et al., 2018, Comput. Phys. Commun., 222, 351

Di Matteo T., 1998, MNRAS, 299, L15

Dubus G., 2006, A&A, 456, 801

Dubus G., 2013, A&AR, 21, 64

Dubus G., Guillard N., Petrucci P.-O., Martin P., 2017, A&A, 608, A59

EHT MWL Science Working Group et al., 2021, ApJ, 911, L11

El Mellah I., Cerutti B., Crinquand B., Parfrey K., 2022, A&A, 663, A169

El Mellah I., Cerutti B., Crinquand B., 2023, A&A, 677, A67

Esirkepov T. Z., 2001, Comput. Phys. Commun., 135, 144

Event Horizon Telescope Collaboration, 2019, ApJ, 875, L6

Event Horizon Telescope Collaboration, 2021, ApJ, 910, L13

Fedeli L. et al., 2022, in SC22: International Conference for High Performance Computing, Networking, Storage and Analysis. IEEE Computer Society, Los Alamitos, CA, USA, p. 1

Fermi LAT Collaboration, 2012, Science, 335, 189

Fonseca R. A. et al., 2002, in Sloot P. M. A., Hoekstra A. G., Tan C. J. K., Dongarra J. J., eds, Computational Science – ICCS 2002. Springer, Berlin, Heidelberg, p. 342

Fossati G., Maraschi L., Celotti A., Comastri A., Ghisellini G., 1998, MNRAS, 299, 433

French O., Guo F., Zhang Q., Uzdensky D. A., 2023, ApJ, 948, 19

Galeev A. A., Rosner R., Vaiana G. S., 1979, ApJ, 229, 318

Galishnikova A., Philippov A., Quataert E., Bacchini F., Parfrey K., Ripperda B., 2023, Phys. Rev. Lett., 130, 115201

Gebhardt K., Adams J., Richstone D., Lauer T. R., Faber S. M., Gültekin K., Murphy J., Tremaine S., 2011, ApJ, 729, 119

Ghisellini G., 2011, in Aharonian F. A., Hofmann W., Rieger F. M., eds, AIP Conf. Proc. Vol. 1381, 25th Texas Symposium on Relativistic AstroPhysics (Texas 2010). Am. Inst. Phys., New York, p. 180

Ghisellini G., Righi C., Costamante L., Tavecchio F., 2017, MNRAS, 469, 255

Giannios D., 2013, MNRAS, 431, 355

Giannios D., Uzdensky D. A., 2019, MNRAS, 484, 1378

Giannios D., Uzdensky D. A., Begelman M. C., 2009, MNRAS, 395, L29

Giannios D., Uzdensky D. A., Begelman M. C., 2010, MNRAS, 402, 1649

Gierliński M., Zdziarski A. A., 2003, MNRAS, 343, L84

Goodbred M., Liu Y.-H., 2022, Phys. Rev. Lett., 129, 265101

Goodman J., Uzdensky D., 2008, ApJ, 688, 555

Gould R. J., Schréder G. P., 1967, Phys. Rev., 155, 1404

Greenwood A. D., Cartwright K. L., Luginsland J. W., Baca E. A., 2004, J. Comput. Phys., 201, 665

Grismayer T., Vranic M., Martins J. L., Fonseca R. A., Silva L. O., 2016, Phys. Plasmas, 23, 056706

Groselj D., Hakobyan H., Beloborodov A. M., Sironi L., Philippov A., 2023, preprint (arXiv:2301.11327)

Guo F., Li H., Daughton W., Liu Y.-H., 2014, Phys. Rev. Lett., 113, 155005

Guo F., Liu Y.-H., Daughton W., Li H., 2015, ApJ, 806, 167

Guo F., Li X., Daughton W., Kilian P., Li H., Liu Y.-H., Yan W., Ma D., 2019, ApJ, 879, L23

Guo F., Li X., Daughton W., Li H., Kilian P., Liu Y.-H., Zhang Q., Zhang H., 2021, ApJ, 919, 111

H. E. S. S. Collaboration, 2019, A&A, 627, A159

H. E. S. S. Collaboration, 2020, A&A, 633, A162

H. E. S. S. Collaboration, 2021, A&A, 648, A23

Hakobyan H., Philippov A., Spitkovsky A., 2019, ApJ, 877, 53

Hakobyan H., Petropoulou M., Spitkovsky A., Sironi L., 2021, ApJ, 912, 48

Hakobyan H., Spitkovsky A., Chernoglazov A., Philippov A., Groselj D., Mahlmann J., 2023a, PrincetonUniversity/tristan-mp-v2: v2.6, Zenodo. Available at:https://zenodo.org/records/7566725.

Hakobyan H., Ripperda B., Philippov A. A., 2023b, ApJ, 943, L29

Harrison M. G., Neukirch T., 2009, Phys. Rev. Lett., 102, 135003

Hayashida M. et al., 2012, ApJ, 754, 114

Hayashida M. et al., 2015, ApJ, 807, 79

Hockney R. W., Eastwood J. W., 1988, Computer Simulation using Particles. Hilger, Bristol

Hoshino M., Lyubarsky Y., 2012, Space Sci. Rev., 173, 521

Igumenshchev I. V., 2008, ApJ, 677, 317

Johnston S., Manchester R. N., Lyne A. G., Bailes M., Kaspi V. M., Qiao G., D'Amico N., 1992, ApJ, 387, L37

Jones F. C., 1968, Phys. Rev., 167, 1159

Khangulyan D., Aharonian F. A., Bogovalov S. V., Ribó M., 2012, ApJ, 752, L17

Kim J. Y. et al., 2018, A&A, 616, A188

Kimura S. S., Toma K., Noda H., Hada K., 2022, ApJ, 937, L34

Kirk J. G., Skjæræasen O., 2003, ApJ, 591, 366

Kirk J. G., Ball L., Skjæræasen O., 1999, Astropart. Phys., 10, 31

Kirk J. G., Lyubarsky Y., Petri J., 2009, in Becker W.ed., Astrophysics and Space Science Library, Vol. 357, Neutron Stars and Pulsars. Springer, Berlin Heidelberg, p. 421

Levinson A., Cerutti B., 2018, A&A, 616, A184

Li Y.-P., Yuan F., Wang Q. D., 2017, MNRAS, 468, 2552

Li X., Guo F., Li H., Li S., 2018, ApJ, 866, 4

Lin S. L., Bardsley J. N., 1977, J. Chem. Phys., 66, 435

Ly C., Walker R. C., Junor W., 2007, ApJ, 660, 200

Lyubarsky Y. E., 2003, MNRAS, 345, 153

Lyubarsky Y. E., 2005, MNRAS, 358, 113

Lyubarsky Y., Kirk J. G., 2001, ApJ, 547, 437

Lyutikov M., Uzdensky D., 2003, ApJ, 589, 893

MAGIC Collaboration, 2008, Science, 320, 1752

MAGIC Collaboration, 2018, A&A, 619, A159

MAGIC Collaboration, 2020, MNRAS, 492, 5354
MAGIC Collaboration, 2021, A&A, 647, A163

Madejski G. G., Sikora M., 2016, ARA&A, 54, 725

Maraschi L., Ghisellini G., Celotti A., 1992, ApJ, 397, L5

Mehlhaff J. M., Werner G. R., Uzdensky D. A., Begelman M. C., 2020, MNRAS, 498, 799

Mehlhaff J. M., Werner G. R., Uzdensky D. A., Begelman M. C., 2021, MNRAS, 508, 4532

Melia F, Königl A., 1989, ApJ, 340, 162

Melzani M., Walder R., Folini D., Winisdoerffer C., Favre J. M., 2014, A&A, 570, A112

Meyer M., Scargle J. D., Blandford R. D., 2019, ApJ, 877, 39

Michel F. C., 1994, ApJ, 431, 397

Mizuno Y., 2022, Universe, 8, 85

Moderski R., Sikora M., Coppi P. S., Aharonian F., 2005, MNRAS, 363, 954

Nalewajko K., 2016, Galaxies, 4, 28

Nalewajko K., Giannios D., Begelman M. C., Uzdensky D. A., Sikora M., 2011, MNRAS, 413, 333

Nalewajko K., Begelman M. C., Cerutti B., Uzdensky D. A., Sikora M., 2012, MNRAS, 425, 2519

Nalewajko K., Yuan Y., Chrušlińska M., 2018, J. Plasma Phys., 84, 755840301

Nätilil J., Beloborodov A. M., 2021, ApJ, 921, 87

Nenkova M., Sirocky M. M., Ivezić Ž., Elitzur M., 2008a, ApJ, 685, 147

Nenkova M., Sirocky M. M., Nikutta R., Ivezić Ž., Elitzur M., 2008b, ApJ, 685, 160

Neukirch T., Wilson F., Harrison M. G., 2009, Phys. Plasmas, 16, 122102

Olguín-Iglesias A. et al., 2016, MNRAS, 460, 3202

Ortuño-Macías J., Nalewajko K., 2020, MNRAS, 497, 1365

Petropoulou M., Sironi L., 2018, MNRAS, 481, 5687

Petropoulou M., Giannios D., Sironi L., 2016, MNRAS, 462, 3325

Philippov A. A., Spitkovsky A., 2018, ApJ, 855, 94

Porth O., Mizuno Y., Younsi Z., Fromm C. M., 2021, MNRAS, 502, 2023

Pozdnyakov L. A., Sobol I. M., Syunyaev R. A., 1983, Astrophys. Space
Phys. Res., 2, 189

Prandini E., Ghisellini G., 2022, Galaxies, 10, 35

Rees H. D., 1968, Phys. Lett. A, 26, 416

Ripperda B., Bacchini F., Philippov A. A., 2020, ApJ, 900, 100

Ripperda B., Liska M., Chatterjee K., Musoke G., Philippov A. A., Markoff
S. B., Tchekhovskoy A., Younsi Z., 2022, ApJ, 924, L32

Rybicki G. B., Lightman A. P., 1979, Radiative Processes in Astrophysics.
Wiley, New York

Scepi N., Dexter J., Begelman M. C., 2022, MNRAS, 511, 3536

Schoeffler K. M., Grismayer T., Uzdensky D., Fonseca R. A., Silva L. O., 2019, ApJ, 870, 49

Schoeffler K. M., Grismayer T., Uzdensky D., Silva L. O., 2023, MNRAS, 523, 3812

Shakura N. I., Sunyaev R. A., 1973, A&A, 24, 337

Sierpowska A., Bednarek W., 2005, MNRAS, 356, 711

Sierpowska-Bartosik A., Bednarek W., 2008, MNRAS, 385, 2279

Sikora M., Begelman M. C., Rees M. J., 1994, ApJ, 421, 153

Sikora M., Stawarz Ł., Moderski R., Nalewajko K., Madejski G. M., 2009, ApJ, 704, 38

Sironi L., 2022, Phys. Rev. Lett., 128, 145102

Sironi L., Beloborodov A. M., 2020, ApJ, 899, 52

Sironi L., Spitkovsky A., 2014, ApJ, 783, L21

Sironi L., Giannios D., Petropoulou M., 2016, MNRAS, 462, 48

Sironi L., Rowan M. E., Narayan R., 2021, ApJ, 907, L44

Sitarek J. et al., 2015, in 34th International Cosmic Ray Conference (ICRC2015). p. 825, preprint (arXiv:1508.04580)

Sobacchi E., Nätilil J., Sironi L., 2021, MNRAS, 503, 688

Sridhar N., Sironi L., Beloborodov A. M., 2021, MNRAS, 507, 5625

Sridhar N., Sironi L., Beloborodov A. M., 2023, MNRAS, 518, 1301

Tam P. H. T., Huang R. H. H., Takata J., Hui C. Y., Kong A. K. H., Cheng K. S., 2011, ApJ, 736, L10

Tam P. H. T., He X. B., Pal P. S., Cui Y., 2018, ApJ, 862, 165

Tamburini M., Pegoraro F., Di Piazza A., Keitel C. H., Macchi A., 2010, New J. Phys., 12, 123005

Tavecchio F., Ghisellini G., 2008, MNRAS, 386, 945

Tchekhovskoy A., Narayan R., McKinney J. C., 2011, MNRAS, 418, L79

Timokhin A. N., Harding A. K., 2015, ApJ, 810, 144

Timokhin A. N., Harding A. K., 2019, ApJ, 871, 12

Uzdensky D. A., 2011, Space Sci. Rev., 160, 45

Uzdensky D. A., 2016, Radiative Magnetic Reconnection in Astrophysics, Vol. 427. Springer International Publishing, Cham, p. 473

Uzdensky D. A., 2022, J. Plasma Phys., 88, 905880114

Uzdensky D. A., Goodman J., 2008, ApJ, 682, 608

Uzdensky D. A., Loureiro N. F., Schekochihin A. A., 2010, Phys. Rev. Lett., 105, 235002

Uzdensky D. et al., 2019, BAAS, 51, 362

Vranic M., Grismayer T., Martins J. L., Fonseca R. A., Silva L. O., 2015, Comput. Phys. Commun., 191, 65

Wakely S. P., Horan D., 2008, in Caballero R., D'Olivo J. C., Medina-Tanco G., Nellen L., Sánchez F. A., Valdés-Galicia J. F., eds, International Cosmic Ray Conference, Vol. 3. Universidad Nacional Autónoma de México, Mexico City, Mexico, p. 1341, Available at: http://tevcat.uchicago.edu

Walker R. C., Hardee P. E., Davies F. B., Ly C., Junor W., 2018, ApJ, 855, 128

Watanabe N., Yokoyama T., 2006, ApJ, 647, L123

Werner G. R., Uzdensky D. A., 2017, ApJ, 843, L27

Werner G. R., Uzdensky D. A., Cerutti B., Nalewajko K., Begelman M. C., 2016, ApJ, 816, L8

Werner G. R., Uzdensky D. A., Begelman M. C., Cerutti B., Nalewajko K., 2018, MNRAS, 473, 4840

Werner G. R., Philippov A. A., Uzdensky D. A., 2019, MNRAS, 482, L60

Xie Y. et al., 2023, Sci. Bull., 68, 1857

Yao P. Z., Dexter J., Chen A. Y., Ryan B. R., Wong G. N., 2021, MNRAS, 507, 4864

Zacharias M. et al., 2017, in AIP Conf. Proc. Vol. 1792, 6th International Symposium on High Energy Gamma-Ray Astronomy. Am. Inst. Phys., New York, p. 050023

Zdziarski A. A., Malyshev D., Chernyakova M., Pooley G. G., 2017, MNRAS, 471, 3657

Zenitani S., Hoshino M., 2001, ApJ, 562, L63

Zenitani S., Hoshino M., 2008, ApJ, 677, 530

Zhang H., Sironi L., Giannios D., 2021, ApJ, 922, 261

Zhang H., Li X., Giannios D., Guo F., Thiersen H., Böttcher M., Lewis T., Venters T., 2022, ApJ, 924, 90

Zhang H., Sironi L., Giannios D., Petropoulou M., 2023, ApJ, 956, L36

Zhdankin V., Uzdensky D. A., Werner G. R., Begelman M. C., 2020, MNRAS, 493, 603

Zhdankin V., Uzdensky D. A., Kunz M. W., 2021, ApJ, 908, 71

APPENDIX A: APPROXIMATE SCATTERED
PHOTON ENERGY ACCURACY

Here, we demonstrate that approximation (26) preserves the average energy, ⟨ε<sub>scat</sub⟩, to which seed photons are IC scattered by a particle of a given Lorentz factor, γ. This is done in Fig. A1, where we compare the exact mean scattered photon energy to the one obtained after many Monte Carlo trials using equation (26). The curves overlap well across the presented range of γ/γ<sub>KN</sub>, which covers the entire range of γ/γ<sub>KN</sub> encountered across all simulations presented in the main text. Because the total cross-section (23) is not approximated, the preservation of ⟨ε<sub>scat</sub>⟩ demonstrated here also implies the correct total IC power (42).
MAGIC Collaboration, 2021, A&A, 647, A163

Madejski G. G., Sikora M., 2016, ARA&A, 54, 725

Maraschi L., Ghisellini G., Celotti A., 1992, ApJ, 397, L5

Mehlhaff J. M., Werner G. R., Uzdensky D. A., Begelman M. C., 2020, MNRAS, 498, 799

Mehlhaff J. M., Werner G. R., Uzdensky D. A., Begelman M. C., 2021, MNRAS, 508, 4532

Melia F, Königl A., 1989, ApJ, 340, 162

Melzani M., Walder R., Folini D., Winisdoerffer C., Favre J. M., 2014, A&A, 570, A112

Meyer M., Scargle J. D., Blandford R. D., 2019, ApJ, 877, 39

Michel F. C., 1994, ApJ, 431, 397

Mizuno Y., 2022, Universe, 8, 85

Moderski R., Sikora M., Coppi P. S., Aharonian F., 2005, MNRAS, 363, 954

Nalewajko K., 2016, Galaxies, 4, 28

Nalewajko K., Giannios D., Begelman M. C., Uzdensky D. A., Sikora M., 2011, MNRAS, 413, 333

Nalewajko K., Begelman M. C., Cerutti B., Uzdensky D. A., Sikora M., 2012, MNRAS, 425, 2519

Nalewajko K., Yuan Y., Chrušlińska M., 2018, J. Plasma Phys., 84, 755840301

Nätilil J., Beloborodov A. M., 2021, ApJ, 921, 87

Nenkova M., Sirocky M. M., Ivezić Ž., Elitzur M., 2008a, ApJ, 685, 147

Nenkova M., Sirocky M. M., Nikutta R., Ivezić Ž., Elitzur M., 2008b, ApJ, 685, 160

Neukirch T., Wilson F., Harrison M. G., 2009, Phys. Plasmas, 16, 122102

Olguín-Iglesias A. et al., 2016, MNRAS, 460, 3202

Ortuño-Macías J., Nalewajko K., 2020, MNRAS, 497, 1365

Petropoulou M., Sironi L., 2018, MNRAS, 481, 5687

Petropoulou M., Giannios D., Sironi L., 2016, MNRAS, 462, 3325

Philippov A. A., Spitkovsky A., 2018, ApJ, 855, 94

Porth O., Mizuno Y., Younsi Z., Fromm C. M., 2021, MNRAS, 502, 2023

Pozdnyakov L. A., Sobol I. M., Syunyaev R. A., 1983, Astrophys. Space
Phys. Res., 2, 189

Prandini E., Ghisellini G., 2022, Galaxies, 10, 35

Rees H. D., 1968, Phys. Lett. A, 26, 416

Ripperda B., Bacchini F., Philippov A. A., 2020, ApJ, 900, 100

Ripperda B., Liska M., Chatterjee K., Musoke G., Philippov A. A., Markoff
S. B., Tchekhovskoy A., Younsi Z., 2022, ApJ, 924, L32

Rybicki G. B., Lightman A. P., 1979, Radiative Processes in Astrophysics.
Wiley, New York

Scepi N., Dexter J., Begelman M. C., 2022, MNRAS, 511, 3536

Schoeffler K. M., Grismayer T., Uzdensky D., Fonseca R. A., Silva L. O., 2019, ApJ, 870, 49

Schoeffler K. M., Grismayer T., Uzdensky D., Silva L. O., 2023, MNRAS, 523, 3812

Shakura N. I., Sunyaev R. A., 1973, A&A, 24, 337

Sierpowska A., Bednarek W., 2005, MNRAS, 356, 711

Sierpowska-Bartosik A., Bednarek W., 2008, MNRAS, 385, 2279

Sikora M., Begelman M. C., Rees M. J., 1994, ApJ, 421, 153

Sikora M., Stawarz Ł., Moderski R., Nalewajko K., Madejski G. M., 2009, ApJ, 704, 38

Sironi L., 2022, Phys. Rev. Lett., 128, 145102

Sironi L., Beloborodov A. M., 2020, ApJ, 899, 52

Sironi L., Spitkovsky A., 2014, ApJ, 783, L21

Sironi L., Giannios D., Petropoulou M., 2016, MNRAS, 462, 48

Sironi L., Rowan M. E., Narayan R., 2021, ApJ, 907, L44

Sitarek J. et al., 2015, in 34th International Cosmic Ray Conference (ICRC2015). p. 825, preprint (arXiv:1508.04580)

Sobacchi E., Nätilil J., Sironi L., 2021, MNRAS, 503, 688

Sridhar N., Sironi L., Beloborodov A. M., 2021, MNRAS, 507, 5625

Sridhar N., Sironi L., Beloborodov A. M., 2023, MNRAS, 518, 1301

Tam P. H. T., Huang R. H. H., Takata J., Hui C. Y., Kong A. K. H., Cheng K. S., 2011, ApJ, 736, L10

Tam P. H. T., He X. B., Pal P. S., Cui Y., 2018, ApJ, 862, 165

Tamburini M., Pegoraro F., Di Piazza A., Keitel C. H., Macchi A., 2010, New
J. Phys., 12, 123005

Tavecchio F., Ghisellini G., 2008, MNRAS, 386, 945

Tchekhovskoy A., Narayan R., McKinney J. C., 2011, MNRAS, 418, L79

Timokhin A. N., Harding A. K., 2015, ApJ, 810, 144

Timokhin A. N., Harding A. K., 2019, ApJ, 871, 12

Uzdensky D. A., 2011, Space Sci. Rev., 160, 45

Uzdensky D. A., 2016, Radiative Magnetic Reconnection in Astrophysics, Vol. 427. Springer International Publishing, Cham, p. 473

Uzdensky D. A., 2022, J. Plasma Phys., 88, 905880114

Uzdensky D. A., Goodman J., 2008, ApJ, 682, 608

Uzdensky D. A., Loureiro N. F., Schekochihin A. A., 2010, Phys. Rev. Lett., 105, 235002

Uzdensky D. et al., 2019, BAAS, 51, 362

Vranic M., Grismayer T., Martins J. L., Fonseca R. A., Silva L. O., 2015, Comput. Phys. Commun., 191, 65

Wakely S. P., Horan D., 2008, in Caballero R., D'Olivo J. C., Medina-Tanco G., Nellen L., Sánchez F. A., Valdés-Galicia J. F., eds, International Cosmic Ray Conference, Vol. 3. Universidad Nacional Autónoma de México, Mexico City, Mexico, p. 1341, Available at: http://tevcat.uchicago.edu

Walker R. C., Hardee P. E., Davies F. B., Ly C., Junor W., 2018, ApJ, 855, 128

Watanabe N., Yokoyama T., 2006, ApJ, 647, L123

Werner G. R., Uzdensky D. A., 2017, ApJ, 843, L27

Werner G. R., Uzdensky D. A., Cerutti B., Nalewajko K., Begelman M. C., 2016, ApJ, 816, L8

Werner G. R., Uzdensky D. A., Begelman M. C., Cerutti B., Nalewajko K., 2018, MNRAS, 473, 4840

Werner G. R., Philippov A. A., Uzdensky D. A., 2019, MNRAS, 482, L60

Xie Y. et al., 2023, Sci. Bull., 68, 1857

Yao P. Z., Dexter J., Chen A. Y., Ryan B. R., Wong G. N., 2021, MNRAS, 507, 4864

Zacharias M. et al., 2017, in AIP Conf. Proc. Vol. 1792, 6th International Symposium on High Energy Gamma-Ray Astronomy. Am. Inst. Phys., New York, p. 050023

Zdziarski A. A., Malyshev D., Chernyakova M., Pooley G. G., 2017, MNRAS, 471, 3657

Zenitani S., Hoshino M., 2001, ApJ, 562, L63

Zenitani S., Hoshino M., 2008, ApJ, 677, 530

Zhang H., Sironi L., Giannios D., 2021, ApJ, 922, 261

Zhang H., Li X., Giannios D., Guo F., Thiersen H., Böttcher M., Lewis T., Venters T., 2022, ApJ, 924, 90

Zhang H., Sironi L., Giannios D., Petropoulou M., 2023, ApJ, 956, L36

Zhdankin V., Uzdensky D. A., Werner G. R., Begelman M. C., 2020, MNRAS, 493, 603

Zhdankin V., Uzdensky D. A., Kunz M. W., 2021, ApJ, 908, 71

APPENDIX A: APPROXIMATE SCATTERED
PHOTON ENERGY ACCURACY

Here, we demonstrate that approximation (26) preserves the average energy, $\langle \epsilon_{scat} \rangle$, to which seed photons are IC scattered by a particle of a given Lorentz factor, $\gamma$. This is done in Fig. A1, where we compare the exact mean scattered photon energy to the one obtained after many Monte Carlo trials using equation (26). The curves overlap well across the presented range of $\gamma/\gamma_{KN}$, which covers the entire range of $\gamma/\gamma_{KN}$ encountered across all simulations presented in the main text. Because the total cross-section (23) is not approximated, the preservation of $\langle \epsilon_{scat} \rangle$ demonstrated here also implies the correct total IC power (42).
<figure><img src="image_19.png" /><figcaption>**Figure A1.** The exact (solid line) and approximate (points) mean final photon energy, $\langle \epsilon_{scat} \rangle$, after a seed photon (energy $\epsilon_{ph}$) undergoes IC scattering with a particle of Lorentz factor $\gamma$. The exact curve (cf. fig. 3 of Mehlhaff et al. 2021) is computed by dividing equation (42) by equation (5). The approximate data are obtained by averaging over $10^5$ Monte Carlo trials per $\gamma$ value, each one using the procedure described in Section 2.2.3 that employs approximation (26).</figcaption></figure>

## APPENDIX B: LIKELIHOOD OF LIMIT CYCLES

Here, we expand on Section 4.6, using measurements from our IC(KN)+PP run (Section 3) to fill in the main quantitative uncertainties from the work of Mehlhaff et al. (2021) on the pair-loading efficiency – the parameter $\xi$ in what follows – including its implications for the pair-loaded magnetization, $\sigma_h^*$, and on the possibility of $\sigma_h$-mediated limit cycles. We report first the expressions for the produced-particle energy density, $U_{prod.}$, and the pair-loaded hot magnetization, $\sigma_h$, flowing into the reconnection layer from the upstream region (thus, in the context of Fig. 14, both quantities should be evaluated along, or perhaps just upstream of, the reconnection separatrix). Mehlhaff et al. (2021) found that $U_{prod.}$ evaluated at this location can be written as

$$U_{prod.} = \xi \mathcal{F} \frac{B_0^2}{8\pi}. \qquad (B1)$$

Here, $\mathcal{F}$ is the fraction of the IC power, emitted by particles in the reconnection layer, that is radiated above pair-production threshold, and $\xi$ is the pair-loading efficiency, which we refer to from here onward according to a more precise name signifying its role in equation (B1), the energy recapture efficiency: how much of the above-threshold radiation is recaptured by the reconnection layer from the inflow region as hot pairs. Accounting for the newborn pair energy density $U_{prod.}$, the hot magnetization, $\sigma_h$, can be written in terms of $\mathcal{F}$ and $\xi$ as

$$\sigma_h = \frac{\sigma_{h,0}}{1 + 2\mathcal{F}\xi\sigma_{h,0}/3} \approx \frac{3}{2\mathcal{F}\xi}, \qquad (B2)$$

where the approximation holds when $\mathcal{F}\xi\sigma_{h,0} \gg 1$. Provided that the reconnection layer responds only to the present – as opposed to the past – magnetization of the plasma feeding it, equation (B2) encodes a universal ($\sigma_{h,0}$-independent) pair-loaded magnetization,

$$\sigma_h^* = \frac{3}{2\mathcal{F}\xi}, \qquad (B3)$$

which we first introduced qualitatively in Section 4.6. The pair-loaded magnetization $\sigma_h^*$ is determined self-consistently by the (potentially $\sigma_h^*$-dependent) values of $\mathcal{F}$ and $\xi$. As discussed in Section 4.6, any Klein–Nishina reconnection layer with initial magnetization

$\sigma_{h,0} > \sigma_h^*$ will try to self-regulate to $\sigma_h^*$, but, if the self-regulation is too efficient, the system will enter a limit cycle between a high and low magnetization on either side of $\sigma_h^*$.

One major uncertainty in the model developed by Mehlhaff et al. (2021) was the energy recapture efficiency, $\xi$. This parameter is particularly important for two reasons. First, it is one of the key quantities deciding the pair-loaded magnetization, $\sigma_h^*$. Second, it determines how strongly the reconnection layer is coupled, via the upstream plasma, to its own NTPA, and, in particular, whether this coupling is strong enough to drive the system into a limit cycle. For a wide range of assumptions, Mehlhaff et al. (2021) found that limit cycles require $\xi$ of order unity – at least 0.3 or so in the most lenient case, but often even higher. In order to make contact with these two issues – the value of $\sigma_h^*$ and the existence (or not) of limit cycles – we now proceed to measure $\xi$ and $\sigma_h^*$ from our simulation.

We estimate the above-threshold power fraction, $\mathcal{F}$, as the part of the average IC emission spectrum in Fig. 9 radiated above pair-production threshold, measuring $\mathcal{F} \simeq 0.5$. This agrees with the prediction one arrives at by taking, from Fig. 7, a particle distribution power law, $dN/d\gamma \propto \gamma^{-2}$, with a cut-off at $200\gamma_{KN}$ and consulting the corresponding $\mathcal{F}$-value from fig. 10 of Mehlhaff et al. (2021). Armed with $\mathcal{F}$, we can estimate $\xi$ by measuring $U_{prod.}$ along the reconnection separatrix in Fig. 14 and plugging the result into equation (B1). From the figure, we see that, once the pair coat is built up around the reconnection layer, it presents an inflowing produced-particle energy density of about $0.01 B_0^2/8\pi$ at the separatrix crossing. This implies that $\xi \mathcal{F} \sim 0.01$ and, hence, $\xi \sim 0.02$. Finally, inserting $\xi \mathcal{F} \sim 0.01$ into equation (B3) yields $\sigma_h^* \sim 150$. This is about a factor of 10 smaller than $\sigma_{h,0} = 1250$ in the IC(KN)+PP run, and, hence, the contour $\sigma_h = 0.1\sigma_{h,0}$ in Fig. 14 nearly overlaps the $U_{prod.} = 0.01 B_0^2/8\pi$ contour.<fnref n="8" />

According to Mehlhaff et al. (2021), our measured value of $\xi$ is much too low (by about an order of magnitude) for limit cycles to occur: the coupling between the inflow region and the layer is too mild. This is consistent with our simulations, from which we identify no evidence of cyclic behaviour, neither in the reconnection-powered NTPA, nor in the pair loading and resulting magnetization, $\sigma_h$, presented to the reconnection layer.

Let us now examine how our measured value of $\xi$, including its implications on the existence of limit cycles, may generalize under changes of the reconnection parameters (such as to those of real astrophysical systems). To inform this discussion, we summarize here the basic physics that determines $\xi$. Mehlhaff et al. (2021) explain that $\xi$ can be written as the product, $\xi = f_{inj}f_{nocoot}f_{noesc}$, where the three factors on the right-hand side correspond to the three main loss mechanisms that inhibit above-threshold radiation from being recaptured by the reconnection layer from the upstream region. First, $f_{inj}$ is the fraction of the layer-produced above-threshold radiation that successfully traverses the reconnection separatrix into the upstream region. The factor $f_{inj}$ falls below unity when some photons – for example, in the case of extreme beaming of reconnection-accelerated particles along the reconnection layer (e.g. Cerutti, Uzdensky & Begelman 2012a; Cerutti et al. 2012b, 2013, 2014a, b; Mehlhaff

<footnote n="8">⁸In fact, the $\sigma_h = 0.1\sigma_{h,0}$ contour lies slightly farther from the main reconnection X-point than the one for $U_{prod.} = 0.01 B_0^2/8\pi$. The small discrepancy originates from the in-plane magnetic field slightly weakening near reconnection X-points, an effect visible in the cold magnetization maps of Fig. 14. This effect owes to the inflowing magnetic field lines draping themselves across plasmons, hanging from them like the cables of a suspension bridge, and, hence, thinning out as they sink toward the X-points in between.</footnote>
<figure><img src="image_20.png" /><figcaption>**Figure A1.** The exact (solid line) and approximate (points) mean final photon energy, $\langle \epsilon_{scat} \rangle$, after a seed photon (energy $\epsilon_{ph}$) undergoes IC scattering with a particle of Lorentz factor $\gamma$. The exact curve (cf. fig. 3 of Mehlhaff et al. 2021) is computed by dividing equation (42) by equation (5). The approximate data are obtained by averaging over $10^5$ Monte Carlo trials per $\gamma$ value, each one using the procedure described in Section 2.2.3 that employs approximation (26).</figcaption></figure>

## APPENDIX B: LIKELIHOOD OF LIMIT CYCLES

Here, we expand on Section 4.6, using measurements from our IC(KN)+PP run (Section 3) to fill in the main quantitative uncertainties from the work of Mehlhaff et al. (2021) on the pair-loading efficiency – the parameter $\xi$ in what follows – including its implications for the pair-loaded magnetization, $\sigma_h^*$, and on the possibility of $\sigma_h$-mediated limit cycles. We report first the expressions for the produced-particle energy density, $U_{prod.}$, and the pair-loaded hot magnetization, $\sigma_h$, flowing into the reconnection layer from the upstream region (thus, in the context of Fig. 14, both quantities should be evaluated along, or perhaps just upstream of, the reconnection separatrix). Mehlhaff et al. (2021) found that $U_{prod.}$ evaluated at this location can be written as

$$U_{prod.} = \xi \mathcal{F} \frac{B_0^2}{8\pi}. \qquad (B1)$$

Here, $\mathcal{F}$ is the fraction of the IC power, emitted by particles in the reconnection layer, that is radiated above pair-production threshold, and $\xi$ is the pair-loading efficiency, which we refer to from here onward according to a more precise name signifying its role in equation (B1), the energy recapture efficiency: how much of the above-threshold radiation is recaptured by the reconnection layer from the inflow region as hot pairs. Accounting for the newborn pair energy density $U_{prod.}$, the hot magnetization, $\sigma_h$, can be written in terms of $\mathcal{F}$ and $\xi$ as

$$\sigma_h = \frac{\sigma_{h,0}}{1 + 2\mathcal{F}\xi\sigma_{h,0}/3} \approx \frac{3}{2\mathcal{F}\xi}, \qquad (B2)$$

where the approximation holds when $\mathcal{F}\xi\sigma_{h,0} \gg 1$. Provided that the reconnection layer responds only to the present – as opposed to the past – magnetization of the plasma feeding it, equation (B2) encodes a universal ($\sigma_{h,0}$-independent) pair-loaded magnetization,

$$\sigma_h^* = \frac{3}{2\mathcal{F}\xi}, \qquad (B3)$$

which we first introduced qualitatively in Section 4.6. The pair-loaded magnetization $\sigma_h^*$ is determined self-consistently by the (potentially $\sigma_h^*$-dependent) values of $\mathcal{F}$ and $\xi$. As discussed in Section 4.6, any Klein–Nishina reconnection layer with initial magnetization

$\sigma_{h,0} > \sigma_h^*$ will try to self-regulate to $\sigma_h^*$, but, if the self-regulation is too efficient, the system will enter a limit cycle between a high and low magnetization on either side of $\sigma_h^*$.

One major uncertainty in the model developed by Mehlhaff et al. (2021) was the energy recapture efficiency, $\xi$. This parameter is particularly important for two reasons. First, it is one of the key quantities deciding the pair-loaded magnetization, $\sigma_h^*$. Second, it determines how strongly the reconnection layer is coupled, via the upstream plasma, to its own NTPA, and, in particular, whether this coupling is strong enough to drive the system into a limit cycle. For a wide range of assumptions, Mehlhaff et al. (2021) found that limit cycles require $\xi$ of order unity – at least 0.3 or so in the most lenient case, but often even higher. In order to make contact with these two issues – the value of $\sigma_h^*$ and the existence (or not) of limit cycles – we now proceed to measure $\xi$ and $\sigma_h^*$ from our simulation.

We estimate the above-threshold power fraction, $\mathcal{F}$, as the part of the average IC emission spectrum in Fig. 9 radiated above pair-production threshold, measuring $\mathcal{F} \simeq 0.5$. This agrees with the prediction one arrives at by taking, from Fig. 7, a particle distribution power law, $dN/d\gamma \propto \gamma^{-2}$, with a cut-off at $200\gamma_{KN}$ and consulting the corresponding $\mathcal{F}$-value from fig. 10 of Mehlhaff et al. (2021). Armed with $\mathcal{F}$, we can estimate $\xi$ by measuring $U_{prod.}$ along the reconnection separatrix in Fig. 14 and plugging the result into equation (B1). From the figure, we see that, once the pair coat is built up around the reconnection layer, it presents an inflowing produced-particle energy density of about $0.01 B_0^2/8\pi$ at the separatrix crossing. This implies that $\xi \mathcal{F} \sim 0.01$ and, hence, $\xi \sim 0.02$. Finally, inserting $\xi \mathcal{F} \sim 0.01$ into equation (B3) yields $\sigma_h^* \sim 150$. This is about a factor of 10 smaller than $\sigma_{h,0} = 1250$ in the IC(KN)+PP run, and, hence, the contour $\sigma_h = 0.1\sigma_{h,0}$ in Fig. 14 nearly overlaps the $U_{prod.} = 0.01 B_0^2/8\pi$ contour.<fnref n="8" />

According to Mehlhaff et al. (2021), our measured value of $\xi$ is much too low (by about an order of magnitude) for limit cycles to occur: the coupling between the inflow region and the layer is too mild. This is consistent with our simulations, from which we identify no evidence of cyclic behaviour, neither in the reconnection-powered NTPA, nor in the pair loading and resulting magnetization, $\sigma_h$, presented to the reconnection layer.

Let us now examine how our measured value of $\xi$, including its implications on the existence of limit cycles, may generalize under changes of the reconnection parameters (such as to those of real astrophysical systems). To inform this discussion, we summarize here the basic physics that determines $\xi$. Mehlhaff et al. (2021) explain that $\xi$ can be written as the product, $\xi = f_{inj}f_{nocoo}f_{noesc}$, where the three factors on the right-hand side correspond to the three main loss mechanisms that inhibit above-threshold radiation from being recaptured by the reconnection layer from the upstream region. First, $f_{inj}$ is the fraction of the layer-produced above-threshold radiation that successfully traverses the reconnection separatrix into the upstream region. The factor $f_{inj}$ falls below unity when some photons – for example, in the case of extreme beaming of reconnection-accelerated particles along the reconnection layer (e.g. Cerutti, Uzdensky & Begelman 2012a; Cerutti et al. 2012b, 2013, 2014a, b; Mehlhaff

<footnote n="8">⁸In fact, the $\sigma_h = 0.1\sigma_{h,0}$ contour lies slightly farther from the main reconnection X-point than the one for $U_{prod.} = 0.01 B_0^2/8\pi$. The small discrepancy originates from the in-plane magnetic field slightly weakening near reconnection X-points, an effect visible in the cold magnetization maps of Fig. 14. This effect owes to the inflowing magnetic field lines draping themselves across plasmons, hanging from them like the cables of a suspension bridge, and, hence, thinning out as they sink toward the X-points in between.</footnote>
et al. 2020) – produce pairs in the downstream region.<fnref n="9" /> Second, $f_{\text{nocool}}$ is that part of the energy deposited into the upstream region as newborn pairs that is not radiated away while those pairs are re-advected toward the layer. Finally, $f_{\text{noesc}}$ is the fraction of newborn upstream particles that do not escape the system (e.g. by travelling along a field line in the $\pm$x-directions) before being swept (in the $\pm$y-directions) into the layer.

To determine how $\xi$ may change with the reconnection parameters, we sketch the dependence of the three governing factors $f_{\text{inj}}$, $f_{\text{nocool}}$, and $f_{\text{noesc}}$, on these parameters as revealed both by analytic theory (Mehlhaff et al. 2021) and by our simulations. Given the periodic boundaries of our setup, particle escape is impossible. Hence, $f_{\text{noesc}}$, perhaps artificially, equals 1 in the simulations, which are then restricted to probing $f_{\text{nocool}}$ and $f_{\text{inj}}$. Of these two, Mehlhaff et al. (2021) provide analytic estimates of $f_{\text{nocool}}$, showing that, similarly to $\mathcal{F}$, it depends only on the shape of the reconnection-energized particle distribution – for a power law, $dN/d\gamma \propto \gamma^{-p}$, on the index, $p$, and on the high-energy cut-off, $\gamma_c$, normalized by $\gamma_{\text{KN}}$. Those authors found that, for $p > 2$, $f_{\text{nocool}}$ is controlled by the low-energy particles, and its value, $f_{\text{nocool}} \sim 0.1$, is therefore independent of $\gamma_c$. For $p < 2$, $f_{\text{nocool}}$ acquires a weak dependence on $\gamma_c$ – signalling the enhanced importance of particles in the high-energy tail – but even then does not reach order unity unless $\gamma_c \gg 10^3 \gamma_{\text{KN}}$ and $p \simeq 1$. In the context of our IC(KN)+PP run, the Mehlhaff et al. (2021) estimates (e.g. their fig. B2) suggest that $f_{\text{nocool}} \sim 0.1$ for $p \simeq 2$ and $\gamma_c \sim 200 \gamma_{\text{KN}}$ (as in Fig. 7). If we adopt this value, we can derive the empirical measurement, $f_{\text{inj}} \sim \xi/f_{\text{nocool}} \sim 0.2$.

Now, while $f_{\text{nocool}}$ depends on the details of reconnection-powered NTPA, $f_{\text{noesc}}$ and $f_{\text{inj}}$ are dictated instead mostly by the large-scale geometry and kinematics of reconnection. For example, $f_{\text{noesc}}$ depends on the free-streaming time, $\le L/c$, for a newborn upstream particle to vacate the system by following an unreconnected magnetic field line (a more complete discussion of the factors influencing $f_{\text{noesc}}$ can be found in the appendix C of Mehlhaff et al. 2021). As another example, the factor, $f_{\text{inj}}$, depends on the effective width of the reconnection layer and on beaming, both that associated with bulk (e.g. plasmoid-chain) motion and that stemming from kinetic effects near reconnection X-points. These remarks suggest that $f_{\text{noesc}}$ and $f_{\text{inj}}$ are constants (at least for $\sigma_{\text{h,0}} \gg 1$), because the processes deciding them are either generic byproducts of relativistic reconnection (in the case of kinetic beaming; Cerutti et al. 2012b; Mehlhaff et al. 2020) or tied to the large-scale evolution (e.g. of the plasmoid chain), which seems rather insensitive to Klein–Nishina radiative physics, as discussed in Sections 4.1 and 4.2.

This, then, makes clear the utility of introducing the individual factors $f_{\text{inj}}$, $f_{\text{nocool}}$, and $f_{\text{noesc}}$: two of them, $f_{\text{noesc}}$ and $f_{\text{inj}}$, are expected to be roughly constant regardless of the exact reconnection parameter values (as long as we are in the Klein–Nishina reconnection regime), and the remaining factor, $f_{\text{nocool}}$, is one that we can estimate based on a substantially reduced set of parameters – the NTPA power-law index, $p$, and cut-off energy, $\gamma_c$. This enables us, as we seek to extrapolate our simulation results to reason about $\xi$ and $\sigma_h^*$ in astrophysical Klein–Nishina reconnection, to skip an exhaustive exploration of

the high-dimensional radiative reconnection parameter space (e.g. all orderings of the energy scales, $\sigma_c$, $\gamma_{\text{max}}$, $\gamma_{\text{rad, T}}$, $\gamma_{\text{cool}}$, $\gamma_{\text{KN}}$, etc.). Instead, we can focus on just $p$ and $\gamma_c$. While it is true that these two are still determined self-consistently by the much larger underlying parameter space, we leave a quantitative characterization of this dependence to future work, having already noted several general trends in Section 4.3.

The established logical framework in terms of $f_{\text{inj}}$, $f_{\text{nocool}}$, and $f_{\text{noesc}}$ – with just the two independent variables $p$ and $\gamma_c$ – equips us to estimate a global upper bound on $\xi$. In the above-mentioned case of extreme NTPA wherein $p \simeq 1$, $\gamma_c \gg 10^3 \gamma_{\text{KN}}$ and, thus, $f_{\text{nocool}} \sim 1$, the combined factor, $\xi = f_{\text{nocool}}/f_{\text{inj}}/f_{\text{noesc}}$, reduces to $\xi \sim f_{\text{inj}}/f_{\text{noesc}} \le f_{\text{inj}} \sim 0.2$. Thus, even in the most favourable circumstances, $\xi$ is still smaller than the most lenient minimum required value, 0.3, found by Mehlhaff et al. (2021) to produce limit cycles. We emphasize the critical role played by simulations in reaching this result, for they provide the necessary bound on the factor, $f_{\text{inj}}$, that caps $\xi$ to beneath the limit cycle value range. In view of this combined input from analytic theory and simulations, we speculate that $\sigma_h$-mediated limit cycles are unlikely, even in astrophysical instances of Klein–Nishina reconnection.

Finally, let us comment on how varying system parameters may impact the pair-loaded magnetization, $\sigma_h^* = 3/2\xi\mathcal{F}$. The discussion thus far can be translated into an expected lower bound on $\sigma_h^*$ as follows. We have already seen that extreme NTPA ($p \simeq 1$ and $\gamma_c \ge 10^3 \gamma_{\text{KN}}$) yields $f_{\text{nocool}} \sim 1$ and, if one also favourably posits $f_{\text{noesc}} = 1$, a global maximum, $\xi \sim 0.1$. In the same NTPA regime, the above-threshold power fraction, $\mathcal{F}$, attains order unity (fig. 10 of Mehlhaff et al. 2021). Thus, with $\xi\mathcal{F}$ as large as possible given numerical measurements and analytic expectations, we arrive at the minimum possible pair-loaded magnetization, $\sigma_h^* \simeq 3/(2 \times 0.1 \times 1) \sim 15$.

If we now relax these extreme assumptions and allow for the (likely) possibility that NTPA is not quite as efficient, we may estimate a more typical value of $\sigma_h^*$. As explained previously, once the power-law index of the particle distribution softens and/or the high-energy cut-off, $\gamma_c$, falls below $\sim 10^3 \gamma_{\text{KN}}$, $f_{\text{nocool}}$ loses dependence on these two parameters and becomes closer to 0.1. In the same regime (again referring to fig. 10 of Mehlhaff et al. 2021), $\mathcal{F}$ becomes more strongly dependent on NTPA, but still attains $\sim 0.5$ when $p \le 2$ for a broad range of $\gamma_c$. This gives a more typical pair-loaded magnetization (realized in our simulation) of $\sigma_h^* \sim 150(0.5/\mathcal{F})(1/f_{\text{noesc}})$.

To summarize, our IC(KN)+PP run allows us to measure the all-important energy recapture efficiency, $\xi$, revealing a value that is too low, in the context of the work by Mehlhaff et al. (2021), to expect violent limit cycles mediated by pair regulation of $\sigma_h$. Such limit cycles are consequently unlikely in real astrophysical systems. The expected evolutionary pathway taken by astrophysical Klein–Nishina reconnection (as long as the newborn pairs remain few in number) is then the same as that of our main IC(KN)+PP run: provided $\sigma_{\text{h,0}} > \sigma_h^*$, the system will self-regulate directly to the $\sigma_{\text{h,0}}$-independent value, $\sigma_h^*$, without overshooting. However, if $\sigma_{\text{h,0}} < \sigma_h^*$, the system remains at its initial magnetization.

<footnote n="9">⁹Mehlhaff et al. (2021) did not include $f_{\text{inj}}$ explicitly in their model, but discuss its effect in their appendix C.</footnote>

<footnote>This paper has been typeset from a TeX/Latex file prepared by the author.</footnote>
et al. 2020) – produce pairs in the downstream region.<fnref n="9" /> Second, $f_{\text{nocool}}$ is that part of the energy deposited into the upstream region as newborn pairs that is not radiated away while those pairs are re-advected toward the layer. Finally, $f_{\text{noesc}}$ is the fraction of newborn upstream particles that do not escape the system (e.g. by travelling along a field line in the $\pm x$-directions) before being swept (in the $\pm y$-directions) into the layer.

To determine how $\xi$ may change with the reconnection parameters, we sketch the dependence of the three governing factors $f_{\text{inj}}$, $f_{\text{nocool}}$, and $f_{\text{noesc}}$, on these parameters as revealed both by analytic theory (Mehlhaff et al. 2021) and by our simulations. Given the periodic boundaries of our setup, particle escape is impossible. Hence, $f_{\text{noesc}}$, perhaps artificially, equals 1 in the simulations, which are then restricted to probing $f_{\text{nocool}}$ and $f_{\text{inj}}$. Of these two, Mehlhaff et al. (2021) provide analytic estimates of $f_{\text{nocool}}$, showing that, similarly to $\mathcal{F}$, it depends only on the shape of the reconnection-energized particle distribution – for a power law, $dN/d\gamma \propto \gamma^{-p}$, on the index, $p$, and on the high-energy cut-off, $\gamma_c$, normalized by $\gamma_{\text{KN}}$. Those authors found that, for $p > 2$, $f_{\text{nocool}}$ is controlled by the low-energy particles, and its value, $f_{\text{nocool}} \sim 0.1$, is therefore independent of $\gamma_c$. For $p < 2$, $f_{\text{nocool}}$ acquires a weak dependence on $\gamma_c$ – signalling the enhanced importance of particles in the high-energy tail – but even then does not reach order unity unless $\gamma_c \gg 10^3 \gamma_{\text{KN}}$ and $p \simeq 1$. In the context of our IC(KN)+PP run, the Mehlhaff et al. (2021) estimates (e.g. their fig. B2) suggest that $f_{\text{nocool}} \sim 0.1$ for $p \simeq 2$ and $\gamma_c \sim 200 \gamma_{\text{KN}}$ (as in Fig. 7). If we adopt this value, we can derive the empirical measurement, $f_{\text{inj}} \sim \xi/f_{\text{nocool}} \sim 0.2$.

Now, while $f_{\text{nocool}}$ depends on the details of reconnection-powered NTPA, $f_{\text{noesc}}$ and $f_{\text{inj}}$ are dictated instead mostly by the large-scale geometry and kinematics of reconnection. For example, $f_{\text{noesc}}$ depends on the free-streaming time, $\le L/c$, for a newborn upstream particle to vacate the system by following an unreconnected magnetic field line (a more complete discussion of the factors influencing $f_{\text{noesc}}$ can be found in the appendix C of Mehlhaff et al. 2021). As another example, the factor, $f_{\text{inj}}$, depends on the effective width of the reconnection layer and on beaming, both that associated with bulk (e.g. plasmoid-chain) motion and that stemming from kinetic effects near reconnection X-points. These remarks suggest that $f_{\text{noesc}}$ and $f_{\text{inj}}$ are constants (at least for $\sigma_{\text{h,0}} \gg 1$), because the processes deciding them are either generic byproducts of relativistic reconnection (in the case of kinetic beaming; Cerutti et al. 2012b; Mehlhaff et al. 2020) or tied to the large-scale evolution (e.g. of the plasmoid chain), which seems rather insensitive to Klein–Nishina radiative physics, as discussed in Sections 4.1 and 4.2.

This, then, makes clear the utility of introducing the individual factors $f_{\text{inj}}$, $f_{\text{nocool}}$, and $f_{\text{noesc}}$: two of them, $f_{\text{noesc}}$ and $f_{\text{inj}}$, are expected to be roughly constant regardless of the exact reconnection parameter values (as long as we are in the Klein–Nishina reconnection regime), and the remaining factor, $f_{\text{nocool}}$, is one that we can estimate based on a substantially reduced set of parameters – the NTPA power-law index, $p$, and cut-off energy, $\gamma_c$. This enables us, as we seek to extrapolate our simulation results to reason about $\xi$ and $\sigma_h^*$ in astrophysical Klein–Nishina reconnection, to skip an exhaustive exploration of

the high-dimensional radiative reconnection parameter space (e.g. all orderings of the energy scales, $\sigma_c$, $\gamma_{\text{max}}$, $\gamma_{\text{rad, T}}$, $\gamma_{\text{cool}}$, $\gamma_{\text{KN}}$, etc.). Instead, we can focus on just $p$ and $\gamma_c$. While it is true that these two are still determined self-consistently by the much larger underlying parameter space, we leave a quantitative characterization of this dependence to future work, having already noted several general trends in Section 4.3.

The established logical framework in terms of $f_{\text{inj}}$, $f_{\text{nocool}}$, and $f_{\text{noesc}}$ – with just the two independent variables $p$ and $\gamma_c$ – equips us to estimate a global upper bound on $\xi$. In the above-mentioned case of extreme NTPA wherein $p \simeq 1$, $\gamma_c \gg 10^3 \gamma_{\text{KN}}$ and, thus, $f_{\text{nocool}} \sim 1$, the combined factor, $\xi = f_{\text{nocool}}/f_{\text{inj}}/f_{\text{noesc}}$, reduces to $\xi \sim f_{\text{inj}}/f_{\text{noesc}} \le f_{\text{inj}} \sim 0.2$. Thus, even in the most favourable circumstances, $\xi$ is still smaller than the most lenient minimum required value, 0.3, found by Mehlhaff et al. (2021) to produce limit cycles. We emphasize the critical role played by simulations in reaching this result, for they provide the necessary bound on the factor, $f_{\text{inj}}$, that caps $\xi$ to beneath the limit cycle value range. In view of this combined input from analytic theory and simulations, we speculate that $\sigma_h$-mediated limit cycles are unlikely, even in astrophysical instances of Klein–Nishina reconnection.

Finally, let us comment on how varying system parameters may impact the pair-loaded magnetization, $\sigma_h^* = 3/2\xi\mathcal{F}$. The discussion thus far can be translated into an expected lower bound on $\sigma_h^*$ as follows. We have already seen that extreme NTPA ($p \simeq 1$ and $\gamma_c \ge 10^3 \gamma_{\text{KN}}$) yields $f_{\text{nocool}} \sim 1$ and, if one also favourably posits $f_{\text{noesc}} = 1$, a global maximum, $\xi \sim 0.1$. In the same NTPA regime, the above-threshold power fraction, $\mathcal{F}$, attains order unity (fig. 10 of Mehlhaff et al. 2021). Thus, with $\xi\mathcal{F}$ as large as possible given numerical measurements and analytic expectations, we arrive at the minimum possible pair-loaded magnetization, $\sigma_h^* \simeq 3/(2 \times 0.1 \times 1) \sim 15$.

If we now relax these extreme assumptions and allow for the (likely) possibility that NTPA is not quite as efficient, we may estimate a more typical value of $\sigma_h^*$. As explained previously, once the power-law index of the particle distribution softens and/or the high-energy cut-off, $\gamma_c$, falls below $\sim 10^3 \gamma_{\text{KN}}$, $f_{\text{nocool}}$ loses dependence on these two parameters and becomes closer to 0.1. In the same regime (again referring to fig. 10 of Mehlhaff et al. 2021), $\mathcal{F}$ becomes more strongly dependent on NTPA, but still attains $\sim 0.5$ when $p \le 2$ for a broad range of $\gamma_c$. This gives a more typical pair-loaded magnetization (realized in our simulation) of $\sigma_h^* \sim 150(0.5/\mathcal{F})(1/f_{\text{noesc}})$.

To summarize, our IC(KN)+PP run allows us to measure the all-important energy recapture efficiency, $\xi$, revealing a value that is too low, in the context of the work by Mehlhaff et al. (2021), to expect violent limit cycles mediated by pair regulation of $\sigma_h$. Such limit cycles are consequently unlikely in real astrophysical systems. The expected evolutionary pathway taken by astrophysical Klein–Nishina reconnection (as long as the newborn pairs remain few in number) is then the same as that of our main IC(KN)+PP run: provided $\sigma_{\text{h,0}} > \sigma_h^*$, the system will self-regulate directly to the $\sigma_{\text{h,0}}$-independent value, $\sigma_h^*$, without overshooting. However, if $\sigma_{\text{h,0}} < \sigma_h^*$, the system remains at its initial magnetization.

<footnote n="9">⁹Mehlhaff et al. (2021) did not include $f_{\text{inj}}$ explicitly in their model, but discuss its effect in their appendix C.</footnote>

<footnote>This paper has been typeset from a TeX/Latex file prepared by the author.</footnote>