<figure><img src="image_5.png" /><figcaption>**Figure 5.** Reconnected magnetic flux as a function of time in our four main runs: no rad., IC(Th), IC(KN), and IC(KN)+PP. In this plot, we include information from both the top and bottom reconnection layers in our simulations. Thus, the saturation of each simulation's reconnected flux at ~0.5$B_0L$ indicates that the initially available flux in the box, $B_0L_y/2 = B_0L$, has been approximately half-consumed. The active phase of each run is displayed as a dashed box during which the middle 70 per cent of the simulation's total change in reconnected flux occurs. Each run's average reconnection rate (see Fig. 6) is the slope of the secant line through this box.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 6.** Per-simulation average, median, and $1\sigma$ (16th and 84th) percentile reconnection rates. The reconnection rates between runs are statistically indistinguishable.</figcaption></figure>

alters the late-time flux saturation; it has little impact on the rate of reconnection during the active phases of our simulations.

Additionally, we observe some radiatively dependent skew in our measured reconnection rate distributions. Notably, the run IC(KN)+PP is almost bimodal in its reconnection rate distribution. As visible in Fig. 5, this run reconnects relatively quickly during the first half of its active period and slows down during the second half. As a result, the median reconnection rate falls substantially below the mean. This could be a sign of the pair feedback anticipated by Mehlhaff et al. (2021), where the pairs produced in the upstream region load the upstream plasma, thereby reducing $\sigma_h$ and hence inhibiting reconnection (see also the supporting information of Xie et al. 2023), but it is not statistically significant. We examine the issue of pair feedback again in later sections.

## 4.3 Non-thermal particle acceleration

In our simulations, magnetic reconnection results in efficient energy delivery to the plasma particles in the form of non-thermal particle acceleration (NTPA). To illustrate this, we present the particle energy distributions, time-averaged over each simulation's reconnection active phase (Section 4.2), in Fig. 7. While NTPA is efficient in all cases, differences between the runs now begin to emerge. In particular, the slope and extent of the high-energy power-law tail differ depending on the radiative physics involved.

In the non-radiative run, reconnection promptly energizes particles up to Lorentz factors $\gamma \sim \gamma_X \equiv 4\sigma_{c,0}$ with a hard power law, $dN/d\gamma \propto \gamma^{-p}$, of index $p \simeq 1.2$. This is expected in our weak guide field, extremely highly magnetized ($\sigma_{h,0} \sim 10^3$) regime, in which numerous previous studies have found that $p$ tends to decrease (hardening the power law) with increasing $\sigma_{h,0}$, asymptoting close to unity in the large-$\sigma_{h,0}$ limit (e.g. Zenitani & Hoshino 2001, 2008; Guo et al. 2014, 2015, 2019, 2021; Melzani et al. 2014; Sironi & Spitkovsky 2014; Sironi et al. 2016; Werner et al. 2016, 2018; Werner & Uzdensky 2017; Ball, Sironi & Özel 2018; Sironi 2022; Uzdensky 2022). Because $p < 2$, most of the particle kinetic energy is stored in the high-energy tail of the distribution, forcing a departure from the $p \simeq 1.2$ scaling beyond Lorentz factors a few times $\gamma_X \equiv 4\sigma_{c,0}$. Otherwise, the particles would carry more energy than initially available in the reconnecting magnetic field.

Beyond the steepening around $\gamma = \gamma_X$, the particle distribution declines and then eventually sharply cuts off near $\gamma = \gamma_{max}$. Since the energies $\gamma_X \equiv 4\sigma_{c,0}$ and $\gamma_{max}$ are relatively close together by computational necessity in this study (Section 3.2.2), it is difficult to determine whether the transition near $\gamma = \gamma_X$ leads to a steeper (softer) power law, perhaps scaling roughly as $p \simeq 3$, or to an exponential cut-off that later gives way to an even sharper cut-off near $\gamma \simeq \gamma_{max}$. If, of these two possibilities, it is the former case of a transition to a softer power law that is realized, it would match the picture advanced by Petropoulou & Sironi (2018) and Hakobyan et al. (2021) where particles undergo additional slow acceleration after being processed across the reconnection separatrix and becoming trapped inside of adiabatically compressing plasmons.

The addition of radiative cooling markedly changes the signatures of NTPA in our simulations. In the purely Thomson radiative run, radiative losses impose a decisive cut-off on the maximum particle energy at the Lorentz factor $\gamma_{rad,T}$, as first studied by Werner et al. (2019). This cut-off is well below not only the nominal system-size-limited Lorentz factor, $\gamma_{max}$, but also the intrinsic maximum energy, $\gamma_X$, attainable by particles via non-ideal direct acceleration by the reconnection electric field near X-points. Thus, radiative losses compete with even very rapid particle acceleration. Slower secondary acceleration channels are suppressed altogether, and there is no evidence of a secondary power-law component associated with such channels (consistent with earlier studies, e.g. Sironi & Beloborodov 2020; Mehlhaff et al. 2020, including when strong radiative cooling is facilitated by synchrotron losses instead of IC emission, as observed by Hakobyan et al. 2023b). The suppression of non-impulsive acceleration processes by (synchrotron or IC) radiative cooling has also been observed in the context of relativistic turbulence (e.g. Zhdankin et al. 2020; Comisso & Sironi 2021; Nättilä & Beloborodov 2021; Sobacchi, Nättilä & Sironi 2021; Zhdankin, Uzdensky & Kunz 2021). Had we conducted 3D simulations, the available secondary particle acceleration channels would have probably been faster than in the present 2D setup, with energized particles not indefinitely confined to plasmons, but eventually escaping back into the upstream region to surf along the reconnection electric field on Speiser-like
<figure><img src="image_6.png" /><figcaption>**Figure 5.** Reconnected magnetic flux as a function of time in our four main runs: no rad., IC(Th), IC(KN), and IC(KN)+PP. In this plot, we include information from both the top and bottom reconnection layers in our simulations. Thus, the saturation of each simulation's reconnected flux at ~0.5$B_0L$ indicates that the initially available flux in the box, $B_0L_y/2 = B_0L$, has been approximately half-consumed. The active phase of each run is displayed as a dashed box during which the middle 70 per cent of the simulation's total change in reconnected flux occurs. Each run's average reconnection rate (see Fig. 6) is the slope of the secant line through this box.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Figure 6.** Per-simulation average, median, and $1\sigma$ (16th and 84th) percentile reconnection rates. The reconnection rates between runs are statistically indistinguishable.</figcaption></figure>

alters the late-time flux saturation; it has little impact on the rate of reconnection during the active phases of our simulations.

Additionally, we observe some radiatively dependent skew in our measured reconnection rate distributions. Notably, the run IC(KN)+PP is almost bimodal in its reconnection rate distribution. As visible in Fig. 5, this run reconnects relatively quickly during the first half of its active period and slows down during the second half. As a result, the median reconnection rate falls substantially below the mean. This could be a sign of the pair feedback anticipated by Mehlhaff et al. (2021), where the pairs produced in the upstream region load the upstream plasma, thereby reducing $\sigma_h$ and hence inhibiting reconnection (see also the supporting information of Xie et al. 2023), but it is not statistically significant. We examine the issue of pair feedback again in later sections.

## 4.3 Non-thermal particle acceleration

In our simulations, magnetic reconnection results in efficient energy delivery to the plasma particles in the form of non-thermal particle acceleration (NTPA). To illustrate this, we present the particle energy distributions, time-averaged over each simulation's reconnection active phase (Section 4.2), in Fig. 7. While NTPA is efficient in all cases, differences between the runs now begin to emerge. In particular, the slope and extent of the high-energy power-law tail differ depending on the radiative physics involved.

In the non-radiative run, reconnection promptly energizes particles up to Lorentz factors $\gamma \sim \gamma_X \equiv 4\sigma_{c,0}$ with a hard power law, $dN/d\gamma \propto \gamma^{-p}$, of index $p \simeq 1.2$. This is expected in our weak guide field, extremely highly magnetized ($\sigma_{h,0} \sim 10^3$) regime, in which numerous previous studies have found that $p$ tends to decrease (hardening the power law) with increasing $\sigma_{h,0}$, asymptoting close to unity in the large-$\sigma_{h,0}$ limit (e.g. Zenitani & Hoshino 2001, 2008; Guo et al. 2014, 2015, 2019, 2021; Melzani et al. 2014; Sironi & Spitkovsky 2014; Sironi et al. 2016; Werner et al. 2016, 2018; Werner & Uzdensky 2017; Ball, Sironi & Özel 2018; Sironi 2022; Uzdensky 2022). Because $p < 2$, most of the particle kinetic energy is stored in the high-energy tail of the distribution, forcing a departure from the $p \simeq 1.2$ scaling beyond Lorentz factors a few times $\gamma_X \equiv 4\sigma_{c,0}$. Otherwise, the particles would carry more energy than initially available in the reconnecting magnetic field.

Beyond the steepening around $\gamma = \gamma_X$, the particle distribution declines and then eventually sharply cuts off near $\gamma = \gamma_{max}$. Since the energies $\gamma_X \equiv 4\sigma_{c,0}$ and $\gamma_{max}$ are relatively close together by computational necessity in this study (Section 3.2.2), it is difficult to determine whether the transition near $\gamma = \gamma_X$ leads to a steeper (softer) power law, perhaps scaling roughly as $p \simeq 3$, or to an exponential cut-off that later gives way to an even sharper cut-off near $\gamma \simeq \gamma_{max}$. If, of these two possibilities, it is the former case of a transition to a softer power law that is realized, it would match the picture advanced by Petropoulou & Sironi (2018) and Hakobyan et al. (2021) where particles undergo additional slow acceleration after being processed across the reconnection separatrix and becoming trapped inside of adiabatically compressing plasmons.

The addition of radiative cooling markedly changes the signatures of NTPA in our simulations. In the purely Thomson radiative run, radiative losses impose a decisive cut-off on the maximum particle energy at the Lorentz factor $\gamma_{rad,T}$, as first studied by Werner et al. (2019). This cut-off is well below not only the nominal system-size-limited Lorentz factor, $\gamma_{max}$, but also the intrinsic maximum energy, $\gamma_X$, attainable by particles via non-ideal direct acceleration by the reconnection electric field near X-points. Thus, radiative losses compete with even very rapid particle acceleration. Slower secondary acceleration channels are suppressed altogether, and there is no evidence of a secondary power-law component associated with such channels (consistent with earlier studies, e.g. Sironi & Beloborodov 2020; Mehlhaff et al. 2020, including when strong radiative cooling is facilitated by synchrotron losses instead of IC emission, as observed by Hakobyan et al. 2023b). The suppression of non-impulsive acceleration processes by (synchrotron or IC) radiative cooling has also been observed in the context of relativistic turbulence (e.g. Zhdankin et al. 2020; Comisso & Sironi 2021; Nättilä & Beloborodov 2021; Sobacchi, Nättilä & Sironi 2021; Zhdankin, Uzdensky & Kunz 2021). Had we conducted 3D simulations, the available secondary particle acceleration channels would have probably been faster than in the present 2D setup, with energized particles not indefinitely confined to plasmons, but eventually escaping back into the upstream region to surf along the reconnection electric field on Speiser-like
<figure><img src="image_7.png" /><figcaption>**Figure 7.** Time-averaged particle energy distributions for our four main runs. The averaging interval is the same as Fig. 5: the period during which the middle 70 per cent of the change in the reconnected flux occurs. Error envelopes display the $1\sigma$ (16th and 84th) percentiles during the averaging interval in each energy bin. The normalization is arbitrary, but equal on all panels. The IC(KN)+PP distribution is decomposed into contributions from particles that were originally present at time $t=0$ (orig.) and those produced during the simulation (prod.). Klein-Nishina effects lead to NTPA that is intermediate between the non-radiative and Thomson radiative regimes, with a high-energy cut-off similar to the non-radiative case but a power-law scaling intermediate between those of the non-radiative and Thomson-cooled runs.</figcaption></figure>

trajectories, rapidly accelerating to even higher energies (Zhang, Sironi & Giannios 2021; Zhang et al. 2023). Taking a hint from recent synchrotron radiative simulations (Chernoglazov et al. 2023), we expect that such acceleration would nevertheless remain suppressed in the presence of the intense radiative losses of the IC(Th) run, which would cool particles down before they could escape back into the upstream for further acceleration. These expectations must ultimately be checked by a future 3D study, though.

Besides the radiative cut-off at $\gamma_{\text{rad, T}}$, the high-energy tail of the Thomson-cooled particle energy distribution features a steeper power-law scaling, $dN/d\gamma \propto \gamma^{-p}$, with $p \simeq 2.5$. This is also consistent with earlier work (Werner et al. 2019; Mehlhaff et al. 2020), including the pile-up seen at intermediate energies (yielding a peak in $\gamma^2 dN/d\gamma$) just beyond $\gamma_{\text{cool}}$. This pile-up results from intermittent episodes of explosive particle acceleration ignited at plasmoid mergers. In between such episodes, particles are rapidly cooled inside their host plasmons, reaching a typical energy $\gamma_{\text{cool}} L/c \Delta t_{\text{coll}}$ characterized by the time $\Delta t_{\text{coll}}$ between plasmoid collisions. While, in reality, $\Delta t_{\text{coll}}$ is different for different tiers in the plasmoid hierarchy, the relativistic plasmoid motion in the box of size $L$ dictates that it should be $<L/c$. This is consistent with the pile-up in the high-energy tail (peak in

$\gamma^2 dN/d\gamma$) occurring at a few-to-several times $\gamma_{\text{cool}}$ in the IC(Th) panel (top right) of Fig. 7.

Klein-Nishina effects lead to an NTPA regime that is largely intermediate between the non-radiative and strongly Thomson-cooled cases. This is evidenced by both our runs with general Compton losses, IC(KN) and IC(KN)+PP. The Klein-Nishina reduction in radiative efficiency causes the particle energy distribution's tail to exhibit a flatter (harder) scaling, $p \simeq 2$, and to persist to higher energies (definitively cutting off by the time $\gamma \simeq \gamma_{\text{max}} = \gamma_{\text{rad, IC}}$, though perhaps steepening sooner, near $\gamma = \gamma_X$) than when cooling proceeds purely in the Thomson limit. However, the tail is still steeper than that in the non-radiative run. This intermediate behaviour can be understood from the hierarchy of time-scales, $t_{\text{cool, T}}(\gamma) < t_{\text{cool, IC}}(\gamma) < L/c$, which holds at all Lorentz factors $\gamma < \gamma_{\text{max}}$ in these simulations (Fig. 2, $\gamma_{\text{KN}}^{(\text{med})}$ curves).

The degree to which NTPA more resembles that in the Thomson or non-radiative limits depends on the precise value of $\gamma_{\text{KN}}$. Increasing $\gamma_{\text{KN}}$ tends to bring the generalized cooling time, $t_{\text{cool, IC}}(\gamma)$, closer to its Thomson limit, $t_{\text{cool, T}}(\gamma)$ (Fig. 2, $\gamma_{\text{KN}}^{(\text{hi})}$ curves). Once $\gamma_{\text{KN}}$ surpasses $\gamma_{\text{rad, T}}$, particles are forbidden from experiencing significant Klein-Nishina effects, and the system reverts to purely Thomson radiative
<figure><img src="image_8.png" /><figcaption>**Figure 7.** Time-averaged particle energy distributions for our four main runs. The averaging interval is the same as Fig. 5: the period during which the middle 70 per cent of the change in the reconnected flux occurs. Error envelopes display the $1\sigma$ (16th and 84th) percentiles during the averaging interval in each energy bin. The normalization is arbitrary, but equal on all panels. The IC(KN)+PP distribution is decomposed into contributions from particles that were originally present at time $t=0$ (orig.) and those produced during the simulation (prod.). Klein-Nishina effects lead to NTPA that is intermediate between the non-radiative and Thomson radiative regimes, with a high-energy cut-off similar to the non-radiative case but a power-law scaling intermediate between those of the non-radiative and Thomson-cooled runs.</figcaption></figure>

trajectories, rapidly accelerating to even higher energies (Zhang, Sironi & Giannios 2021; Zhang et al. 2023). Taking a hint from recent synchrotron radiative simulations (Chernoglazov et al. 2023), we expect that such acceleration would nevertheless remain suppressed in the presence of the intense radiative losses of the IC(Th) run, which would cool particles down before they could escape back into the upstream for further acceleration. These expectations must ultimately be checked by a future 3D study, though.

Besides the radiative cut-off at $\gamma_{\text{rad, T}}$, the high-energy tail of the Thomson-cooled particle energy distribution features a steeper power-law scaling, $dN/d\gamma \propto \gamma^{-p}$, with $p \simeq 2.5$. This is also consistent with earlier work (Werner et al. 2019; Mehlhaff et al. 2020), including the pile-up seen at intermediate energies (yielding a peak in $\gamma^2 dN/d\gamma$) just beyond $\gamma_{\text{cool}}$. This pile-up results from intermittent episodes of explosive particle acceleration ignited at plasmoid mergers. In between such episodes, particles are rapidly cooled inside their host plasmons, reaching a typical energy $\gamma_{\text{cool}} L/c \Delta t_{\text{coll}}$ characterized by the time $\Delta t_{\text{coll}}$ between plasmoid collisions. While, in reality, $\Delta t_{\text{coll}}$ is different for different tiers in the plasmoid hierarchy, the relativistic plasmoid motion in the box of size $L$ dictates that it should be $<L/c$. This is consistent with the pile-up in the high-energy tail (peak in

$\gamma^2 dN/d\gamma$) occurring at a few-to-several times $\gamma_{\text{cool}}$ in the IC(Th) panel (top right) of Fig. 7.

Klein-Nishina effects lead to an NTPA regime that is largely intermediate between the non-radiative and strongly Thomson-cooled cases. This is evidenced by both our runs with general Compton losses, IC(KN) and IC(KN)+PP. The Klein-Nishina reduction in radiative efficiency causes the particle energy distribution's tail to exhibit a flatter (harder) scaling, $p \simeq 2$, and to persist to higher energies (definitively cutting off by the time $\gamma \simeq \gamma_{\text{max}} = \gamma_{\text{rad, IC}}$, though perhaps steepening sooner, near $\gamma = \gamma_X$) than when cooling proceeds purely in the Thomson limit. However, the tail is still steeper than that in the non-radiative run. This intermediate behaviour can be understood from the hierarchy of time-scales, $t_{\text{cool, T}}(\gamma) < t_{\text{cool, IC}}(\gamma) < L/c$, which holds at all Lorentz factors $\gamma < \gamma_{\text{max}}$ in these simulations (Fig. 2, $\gamma_{\text{KN}}^{(\text{med})}$ curves).

The degree to which NTPA more resembles that in the Thomson or non-radiative limits depends on the precise value of $\gamma_{\text{KN}}$. Increasing $\gamma_{\text{KN}}$ tends to bring the generalized cooling time, $t_{\text{cool, IC}}(\gamma)$, closer to its Thomson limit, $t_{\text{cool, T}}(\gamma)$ (Fig. 2, $\gamma_{\text{KN}}^{(\text{hi})}$ curves). Once $\gamma_{\text{KN}}$ surpasses $\gamma_{\text{rad, T}}$, particles are forbidden from experiencing significant Klein-Nishina effects, and the system reverts to purely Thomson radiative
<figure><img src="image_8.png" /><figcaption>**Figure 8.** Local average Lorentz factor maps for our four main runs. Each run is pictured at the moment in time when half of its final change in reconnected magnetic flux has occurred. Commensurate with what happens in the particle energy distributions, Klein–Nishina effects result in maps that have properties that are intermediate between the non-radiative and purely Thomson-cooled regimes. Namely, the IC(KN) and IC(KN)+PP runs both feature hot reconnection current sheets like the non-radiative simulation, but plasmons with temperatures between those of the plasmons in the non-radiative and Thomson-cooled cases.</figcaption></figure>

reconnection. On the other hand, reducing $\gamma_{KN}$ lengthens the cooling time, $t_{\text{cool, IC}}(\gamma)$, and, in the extreme case of $\gamma_{KN} \lesssim \gamma_{\text{cool}}$, $t_{\text{cool, IC}}(\gamma)$ exceeds $L/c$ for all $\gamma$: an effectively non-radiative regime (Fig. 2, $\gamma_{KN}^{(\text{lo})}$ curves). We have verified these expectations by running simulations with differing $\gamma_{KN}$ (those outlined in Table 2), though we do not present the NTPA from those simulations in detail. These runs confirm that increasing $\gamma_{KN}$ moves the particle distribution cut-off closer to $\gamma_{\text{rad, T}}$ while also steepening its power-law tail, causing NTPA to resemble the Thomson limit. Conversely, reducing $\gamma_{KN}$ maintains the sharp particle energy cut-off at the system-size limit, $\gamma_{\text{max}}$ (perhaps with an earlier cut-off or power-law transition near $\gamma_X$), while simultaneously hardening the high-energy tail, transitioning the system towards non-radiative NTPA. Thus, the IC(Th) and no rad. runs represent the effective high- and low-$\gamma_{KN}$ extremes, with the IC(KN) and IC(KN)+PP runs falling somewhere in the middle.

In Fig. 7, for the run with pair production, IC(KN)+PP, we decompose the particle energy distribution into its separate contributions from particles that were present originally in the simulation and those produced on the fly. Strikingly, the produced particles dominate not only the high-energy tail, but also the total plasma energy contained in the simulation box. This is in spite of the fact that these particles are far less numerous than their originally present counterparts, which is only possible thanks to their extremely high average energy. That the produced particles should compete with the original particles for energetic dominance despite being fewer in number is in line with the basic predictions for this regime of reconnection advanced by Mehlhaff et al. (2021). We examine more thoroughly the differences between original and produced particles in Section 4.6.

To complement the distributions of Fig. 7 with a spatial view into particle acceleration, we present maps of the local plasma average Lorentz factor for each of our four main runs in Fig. 8. Here again, we see that the non-radiative and Thomson-cooled cases represent two opposite extremes. On the one hand, the non-radiative

run contains hot reconnection current sheets and plasmons (except the cold core of the large exhaust plasmoid centred at $x \simeq 20\sigma_{c,0} \rho_0$ – it is composed of initially drifting plasma swept directly into this plasmoid near the beginning of the simulation without ever being energized by reconnection). On the other hand, the simulation with Thomson losses features cold plasmoid cores (cf. Beloborodov 2017; Sironi & Beloborodov 2020; Sridhar et al. 2021, 2023; Groselj et al. 2023) and merely warm current sheets.

Let us see how these differences arise. The plasmons in the non-radiative run accumulate kinetic energy via the hot plasma that is exhausted away from reconnection X-points, thereby containing, collectively, a running tally of the dissipated magnetic energy. While plasmons still collect particles in the IC(Th) case, they no longer amass liberated energy, which instead escapes as IC radiation. Then, instead of the area-filling, particle-accumulating plasmons, it is the quasi-1D current-sheet singularities that host the energetic plasma – that is, where intense acceleration is actively taking place.<fnref n="5" /> And even in these special regions, the local mean energy is radiatively limited to $\sim\gamma_{\text{rad, T}} = 0.2\sigma_{c,0}$, much lower than the intrinsic maximum X-point acceleration Lorentz factor, $\gamma_X = 4\sigma_{c,0}$, reached in the (consequently much hotter) current sheets of the non-radiative run. In the IC(Th) case, once particles vacate rapid acceleration zones near X-points to move into plasmons, they quickly cool down, giving the plasmoid cores a characteristic mean energy of order the pile-up energy, $\langle\gamma\rangle \sim 10^{-2}\sigma_{c,0}$, in the IC(Th) distribution of Fig. 7.

Moving now to the IC(KN) and IC(KN)+PP simulations in Fig. 8, we see again that they are intermediate between the non-radiative and Thomson-cooled extremes. Like the non-radiative case, these two runs contain very hot reconnection current layers, with local Lorentz factors comparable to $\sigma_{c,0}$ and far exceeding $\gamma_{\text{rad, T}}$. This reflects

<footnote n="5">⁵In 3D, the plasmons would be volume-filling instead of area-filling, and the current sheets quasi-2D structures instead of quasi-1D ones.</footnote>
<figure><img src="image_9.png" /><figcaption>**Figure 8.** Local average Lorentz factor maps for our four main runs. Each run is pictured at the moment in time when half of its final change in reconnected magnetic flux has occurred. Commensurate with what happens in the particle energy distributions, Klein–Nishina effects result in maps that have properties that are intermediate between the non-radiative and purely Thomson-cooled regimes. Namely, the IC(KN) and IC(KN)+PP runs both feature hot reconnection current sheets like the non-radiative simulation, but plasmons with temperatures between those of the plasmons in the non-radiative and Thomson-cooled cases.</figcaption></figure>

reconnection. On the other hand, reducing $\gamma_{KN}$ lengthens the cooling time, $t_{\text{cool, IC}}(\gamma)$, and, in the extreme case of $\gamma_{KN} \lesssim \gamma_{\text{cool}}$, $t_{\text{cool, IC}}(\gamma)$ exceeds $L/c$ for all $\gamma$: an effectively non-radiative regime (Fig. 2, $\gamma_{KN}^{(\text{lo})}$ curves). We have verified these expectations by running simulations with differing $\gamma_{KN}$ (those outlined in Table 2), though we do not present the NTPA from those simulations in detail. These runs confirm that increasing $\gamma_{KN}$ moves the particle distribution cut-off closer to $\gamma_{\text{rad, T}}$ while also steepening its power-law tail, causing NTPA to resemble the Thomson limit. Conversely, reducing $\gamma_{KN}$ maintains the sharp particle energy cut-off at the system-size limit, $\gamma_{\text{max}}$ (perhaps with an earlier cut-off or power-law transition near $\gamma_X$), while simultaneously hardening the high-energy tail, transitioning the system towards non-radiative NTPA. Thus, the IC(Th) and no rad. runs represent the effective high- and low-$\gamma_{KN}$ extremes, with the IC(KN) and IC(KN)+PP runs falling somewhere in the middle.

In Fig. 7, for the run with pair production, IC(KN)+PP, we decompose the particle energy distribution into its separate contributions from particles that were present originally in the simulation and those produced on the fly. Strikingly, the produced particles dominate not only the high-energy tail, but also the total plasma energy contained in the simulation box. This is in spite of the fact that these particles are far less numerous than their originally present counterparts, which is only possible thanks to their extremely high average energy. That the produced particles should compete with the original particles for energetic dominance despite being fewer in number is in line with the basic predictions for this regime of reconnection advanced by Mehlhaff et al. (2021). We examine more thoroughly the differences between original and produced particles in Section 4.6.

To complement the distributions of Fig. 7 with a spatial view into particle acceleration, we present maps of the local plasma average Lorentz factor for each of our four main runs in Fig. 8. Here again, we see that the non-radiative and Thomson-cooled cases represent two opposite extremes. On the one hand, the non-radiative

run contains hot reconnection current sheets and plasmons (except the cold core of the large exhaust plasmoid centred at $x \simeq 20\sigma_{c,0} \rho_0$ – it is composed of initially drifting plasma swept directly into this plasmoid near the beginning of the simulation without ever being energized by reconnection). On the other hand, the simulation with Thomson losses features cold plasmoid cores (cf. Beloborodov 2017; Sironi & Beloborodov 2020; Sridhar et al. 2021, 2023; Groselj et al. 2023) and merely warm current sheets.

Let us see how these differences arise. The plasmons in the non-radiative run accumulate kinetic energy via the hot plasma that is exhausted away from reconnection X-points, thereby containing, collectively, a running tally of the dissipated magnetic energy. While plasmons still collect particles in the IC(Th) case, they no longer amass liberated energy, which instead escapes as IC radiation. Then, instead of the area-filling, particle-accumulating plasmons, it is the quasi-1D current-sheet singularities that host the energetic plasma – that is, where intense acceleration is actively taking place.<fnref n="5" /> And even in these special regions, the local mean energy is radiatively limited to $\sim\gamma_{\text{rad, T}} = 0.2\sigma_{c,0}$, much lower than the intrinsic maximum X-point acceleration Lorentz factor, $\gamma_X = 4\sigma_{c,0}$, reached in the (consequently much hotter) current sheets of the non-radiative run. In the IC(Th) case, once particles vacate rapid acceleration zones near X-points to move into plasmons, they quickly cool down, giving the plasmoid cores a characteristic mean energy of order the pile-up energy, $\langle\gamma\rangle \sim 10^{-2}\sigma_{c,0}$, in the IC(Th) distribution of Fig. 7.

Moving now to the IC(KN) and IC(KN)+PP simulations in Fig. 8, we see again that they are intermediate between the non-radiative and Thomson-cooled extremes. Like the non-radiative case, these two runs contain very hot reconnection current layers, with local Lorentz factors comparable to $\sigma_{c,0}$ and far exceeding $\gamma_{\text{rad, T}}$. This reflects

<footnote n="5">⁵In 3D, the plasmons would be volume-filling instead of area-filling, and the current sheets quasi-2D structures instead of quasi-1D ones.</footnote>
the fact that cooling losses do not substantially inhibit acceleration near reconnection X-points in these Klein–Nishina-regime runs. As for plasmons, these are colder than in the non-radiative simulation but warmer than in the IC(Th) case. This stems again from the cooling time-scale hierarchy, $t_{\text{cool, T}}(\gamma) < t_{\text{cool, IC}}(\gamma) < L/c$: particles accelerated at current sheets are not efficiently cooled on the time-scales of their acceleration, but they are still efficiently cooled over one dynamical time, causing plasmons to cool down – just not as quickly as in the IC(Th) run.

To summarize, in this subsection, we have witnessed the first main differences emerge among our simulations, with the differing radiative physics leaving pronounced and distinguishing imprints on NTPA. Even though all runs exhibit non-thermal power-law particle energy distributions, $dN/d\gamma \propto \gamma^{-p}$, in their reconnection active phases, the slopes and extents of their power-laws differ dramatically. The non-radiative regime yields $p$ close to unity, with a departure from this scaling near $\gamma = \gamma_X \equiv 4\sigma_{c,0}$ and a subsequent sharp cut-off at $\gamma = \gamma_{\text{max}}$. The IC(Th) run represents an opposite regime, with $p \simeq 2.5$ followed by an abrupt cut-off at $\gamma \simeq \gamma_{\text{rad, T}} \ll \gamma_{\text{max}}$. Klein–Nishina radiative cooling lies between these two extremes, and for the parameters of our IC(KN) and IC(KN)+PP runs, gives $p \simeq 2$ and a sharp cut-off at $\gamma \simeq \gamma_{\text{max}}$ (perhaps with an earlier steepening near $\gamma = \gamma_X$).

## 4.4 Inverse Compton emission spectra

In this section, we connect the pronounced NTPA in our simulations to their IC emission spectra. Like the underlying particle energy distributions (Section 4.3), these spectra are highly extended and non-thermal, and the distinctions among them reflect the differing radiative physics at play.

We present angle-integrated IC emission spectra for each of our four main runs in Fig. 9. These spectra are computed by summing together the individual spectra from every simulated particle, regardless of the particle’s position or velocity direction. As a result, beaming and light-traveltime effects are neglected. In addition, for the IC(KN)+PP run, pair-production absorption of the emitted gamma-rays is temporarily ignored, but we discuss it briefly at the end of this subsection and in more detail in the next subsection.

In the non-radiative and Thomson-cooled cases, the IC emission spectra are produced in the Thomson regime and are thus given straightforwardly in terms of the underlying particle energy distributions as follows. Particles of Lorentz factor $\gamma$ upscatter photons to energies $\epsilon(\gamma) \sim \gamma^2 \epsilon_{\text{ph}}$. Because the scattering rate, $\sigma_T c n_{\text{ph}}$ [cf. equation (5)], is independent of $\gamma$, the number of photons emitted per unit time into a given energy interval, $d\epsilon dN_{\text{ph}}/d\epsilon d\gamma$, is proportional to the number of particles at the corresponding scattering Lorentz factor, $d\gamma dN/d\gamma|_{\gamma_{\text{alpha}}^{1/2}}$. Thus, if the particle energy distribution is a power law, $dN/d\gamma \propto \gamma^{-p}$, the emitted photon distribution is also a power law: $dN_{\text{ph}}/d\epsilon d\gamma \propto (dN/d\gamma|_{\gamma_{\text{alpha}}^{1/2}})(d\gamma/d\epsilon) \propto \epsilon^{-(p+1)/2}$ (cf. Rybicki & Lightman 1979). In the $\epsilon F_{\text{IC}}(\epsilon)$ representation plotted in Fig. 9, this translates to $\epsilon F_{\text{IC}}(\epsilon) = \epsilon^2 dN_{\text{ph}}/d\epsilon \propto \epsilon^{-(p-3)/2}$.

This result equips us to easily interpret the power-law components in the non-radiative and Thomson-cooled IC emission spectra. The non-radiative run’s particle distribution power law is approximately $dN/d\gamma \propto \gamma^{-1.2}$ (Fig. 7), which yields the expected IC power law, $\epsilon F_{\text{IC}}(\epsilon) \propto \epsilon^{0.9}$. This is in good agreement with the non-radiative IC spectrum in the top left panel of Fig. 9. Additionally, if present, the second/steeper particle power-law component for this same simulation, $dN/d\gamma \propto \gamma^{-3}$ – putatively stemming from slower secondary acceleration channels – should produce a flat spectrum, $\epsilon F_{\text{IC}}(\epsilon) \propto \epsilon^0 = \text{const}$. While this is roughly consistent with the

measured spectrum of Fig. 9, it is difficult to definitively say that such a component truly exists and is not just part of the spectral cut-off. Regarding the Thomson-cooled run, the power law, $dN/d\gamma \propto \gamma^{-2.5}$, should yield the gently increasing spectrum, $\epsilon F_{\text{IC}}(\epsilon) \propto \epsilon^{0.25}$, which agrees with Fig. 9 (top right panel). Finally, all these ideas can be applied not just to connect spectral slopes between the particle and IC spectra, but also to relate their cut-offs. In particular, the cut-offs at $\gamma_{\text{max}}$ and $\gamma_{\text{rad, T}}$ in the respective non-radiative and Thomson-cooled particle distributions (Fig. 7) correspond to the observed cut-offs at photon energies of order $\gamma_{\text{max}}^2 \epsilon_{\text{ph}}$ and $\gamma_{\text{rad, T}}^2 \epsilon_{\text{ph}}$ in these simulations’ IC emission.

This simple framework breaks down in the presence of Klein–Nishina effects. Then, the characteristic scattered photon energy becomes a broken power-law function of the particle’s energy: $\epsilon(\gamma) \sim \gamma^2 \epsilon_{\text{ph}}$ when $\gamma < \gamma_{\text{KN}}$ and $\epsilon(\gamma) \sim \gamma m_e c^2$ otherwise. Furthermore, the scattering rate, $\sigma_T c n_{\text{ph}} \epsilon_{\text{KN}}(\gamma/\gamma_{\text{KN}})$ [equation (5)], becomes a non-trivial, decreasing function of $\gamma$. This suppresses the emission efficiency and breaks the simple correspondence between the particle distribution power-law index and that of the IC emission spectrum. Let us examine how these effects manifest themselves in the IC(KN) and IC(KN)+PP panels (bottom left and bottom right, respectively) of Fig. 9. To begin with, the characteristic emission energies from particles at each of our Lorentz-factor scales are pushed closer together beyond the energy $\epsilon(\gamma_{\text{KN}})$ (because, beyond $\gamma_{\text{KN}}$, $\epsilon$ scales linearly with $\gamma$ instead of quadratically). As an example, even though the scales $\gamma_{\text{max}}$, $\gamma_{\text{rad, T}}$, and $\gamma_{\text{cool}}$ are all equally spaced on a logarithmic scale (because $\gamma_{\text{rad, T}}^2 = \gamma_{\text{max}}^2 \gamma_{\text{cool}}$), the corresponding photon energies, $\epsilon(\gamma_{\text{max}})$, $\epsilon(\gamma_{\text{rad, T}})$, and $\epsilon(\gamma_{\text{cool}})$, are not at all evenly spaced, with $\epsilon(\gamma_{\text{max}})$ and $\epsilon(\gamma_{\text{rad, T}})$ closer together than $\epsilon(\gamma_{\text{rad, T}})$ and $\epsilon(\gamma_{\text{cool}})$ [because the former two lie above the break energy $\epsilon(\gamma_{\text{KN}})$ while the latter lies below it]. Next, even though the power-law scalings of the IC(KN) and IC(KN)+PP particle distributions are both approximately $\gamma^{-2}$ or shallower (i.e. harder), the corresponding $\epsilon F_{\text{IC}}(\epsilon)$ spectra are both steeper (i.e. softer) than the Thomson-limit prediction, $\epsilon^{0.5}$, demonstrating the reduced radiative efficiency in the Klein–Nishina limit.

While the diminished Klein–Nishina cross-section produces a softer emission spectrum for a given particle distribution, it also yields a particle distribution that is harder in the first place (Section 4.3). These two effects somewhat cancel out, and, hence, not much change is observed in the spectral slope from the time-averaged IC(Th) spectrum to those yielded by the IC(KN) and IC(KN)+PP simulations (cf. Moderski et al. 2005). This is even despite the very different shape – the result of different cooling physics – in the particle distributions between these runs.

Finally, we note that the IC(KN) and IC(KN)+PP Compton emission spectra peak far above $\epsilon(\gamma_{\text{KN}})$ – deeply in the Klein–Nishina regime. As a result, most of the radiated energy in the IC(KN)+PP case (grey region in the lower right panel of Fig. 9) is emitted above pair-production threshold, where [unlike the IC(KN) run] it will be recaptured by the system as hot newborn pairs. The peak in the intrinsic emitted IC spectrum is therefore invisible to the observer, who sees only the indirect remnant of this radiation reprocessed to below-threshold energies. We elaborate the observable consequences of this effect in the next subsection.

## 4.5 Light curves and spectral variability

We now complement Section 4.4’s energy-resolved view of the emission from our simulations by discussing the timing of the radiative signatures. When viewed through the lens of timing, the
the fact that cooling losses do not substantially inhibit acceleration near reconnection X-points in these Klein–Nishina-regime runs. As for plasmons, these are colder than in the non-radiative simulation but warmer than in the IC(Th) case. This stems again from the cooling time-scale hierarchy, $t_{\text{cool, T}}(\gamma) < t_{\text{cool, IC}}(\gamma) < L/c$: particles accelerated at current sheets are not efficiently cooled on the time-scales of their acceleration, but they are still efficiently cooled over one dynamical time, causing plasmons to cool down – just not as quickly as in the IC(Th) run.

To summarize, in this subsection, we have witnessed the first main differences emerge among our simulations, with the differing radiative physics leaving pronounced and distinguishing imprints on NTPA. Even though all runs exhibit non-thermal power-law particle energy distributions, $dN/d\gamma \propto \gamma^{-p}$, in their reconnection active phases, the slopes and extents of their power-laws differ dramatically. The non-radiative regime yields $p$ close to unity, with a departure from this scaling near $\gamma = \gamma_X \equiv 4\sigma_{c,0}$ and a subsequent sharp cut-off at $\gamma = \gamma_{\text{max}}$. The IC(Th) run represents an opposite regime, with $p \simeq 2.5$ followed by an abrupt cut-off at $\gamma \simeq \gamma_{\text{rad, T}} \ll \gamma_{\text{max}}$. Klein–Nishina radiative cooling lies between these two extremes, and for the parameters of our IC(KN) and IC(KN)+PP runs, gives $p \simeq 2$ and a sharp cut-off at $\gamma \simeq \gamma_{\text{max}}$ (perhaps with an earlier steepening near $\gamma = \gamma_X$).

## 4.4 Inverse Compton emission spectra

In this section, we connect the pronounced NTPA in our simulations to their IC emission spectra. Like the underlying particle energy distributions (Section 4.3), these spectra are highly extended and non-thermal, and the distinctions among them reflect the differing radiative physics at play.

We present angle-integrated IC emission spectra for each of our four main runs in Fig. 9. These spectra are computed by summing together the individual spectra from every simulated particle, regardless of the particle’s position or velocity direction. As a result, beaming and light-traveltime effects are neglected. In addition, for the IC(KN)+PP run, pair-production absorption of the emitted gamma-rays is temporarily ignored, but we discuss it briefly at the end of this subsection and in more detail in the next subsection.

In the non-radiative and Thomson-cooled cases, the IC emission spectra are produced in the Thomson regime and are thus given straightforwardly in terms of the underlying particle energy distributions as follows. Particles of Lorentz factor $\gamma$ upscatter photons to energies $\epsilon(\gamma) \sim \gamma^2 \epsilon_{\text{ph}}$. Because the scattering rate, $\sigma_T c n_{\text{ph}}$ [cf. equation (5)], is independent of $\gamma$, the number of photons emitted per unit time into a given energy interval, $d\epsilon dN_{\text{ph}}/d\epsilon d\gamma$, is proportional to the number of particles at the corresponding scattering Lorentz factor, $d\gamma dN/d\gamma|_{\gamma_{\text{alpha}}^{1/2}}$. Thus, if the particle energy distribution is a power law, $dN/d\gamma \propto \gamma^{-p}$, the emitted photon distribution is also a power law: $dN_{\text{ph}}/d\epsilon d\gamma \propto (dN/d\gamma|_{\gamma_{\text{alpha}}^{1/2}})(d\gamma/d\epsilon) \propto \epsilon^{-(p+1)/2}$ (cf. Rybicki & Lightman 1979). In the $\epsilon F_{\text{IC}}(\epsilon)$ representation plotted in Fig. 9, this translates to $\epsilon F_{\text{IC}}(\epsilon) = \epsilon^2 dN_{\text{ph}}/d\epsilon \propto \epsilon^{-(p-3)/2}$.

This result equips us to easily interpret the power-law components in the non-radiative and Thomson-cooled IC emission spectra. The non-radiative run’s particle distribution power law is approximately $dN/d\gamma \propto \gamma^{-1.2}$ (Fig. 7), which yields the expected IC power law, $\epsilon F_{\text{IC}}(\epsilon) \propto \epsilon^{0.9}$. This is in good agreement with the non-radiative IC spectrum in the top left panel of Fig. 9. Additionally, if present, the second/steeper particle power-law component for this same simulation, $dN/d\gamma \propto \gamma^{-3}$ – putatively stemming from slower secondary acceleration channels – should produce a flat spectrum, $\epsilon F_{\text{IC}}(\epsilon) \propto \epsilon^0 = \text{const}$. While this is roughly consistent with the

measured spectrum of Fig. 9, it is difficult to definitively say that such a component truly exists and is not just part of the spectral cut-off. Regarding the Thomson-cooled run, the power law, $dN/d\gamma \propto \gamma^{-2.5}$, should yield the gently increasing spectrum, $\epsilon F_{\text{IC}}(\epsilon) \propto \epsilon^{0.25}$, which agrees with Fig. 9 (top right panel). Finally, all these ideas can be applied not just to connect spectral slopes between the particle and IC spectra, but also to relate their cut-offs. In particular, the cut-offs at $\gamma_{\text{max}}$ and $\gamma_{\text{rad, T}}$ in the respective non-radiative and Thomson-cooled particle distributions (Fig. 7) correspond to the observed cut-offs at photon energies of order $\gamma_{\text{max}}^2 \epsilon_{\text{ph}}$ and $\gamma_{\text{rad, T}}^2 \epsilon_{\text{ph}}$ in these simulations’ IC emission.

This simple framework breaks down in the presence of Klein–Nishina effects. Then, the characteristic scattered photon energy becomes a broken power-law function of the particle’s energy: $\epsilon(\gamma) \sim \gamma^2 \epsilon_{\text{ph}}$ when $\gamma < \gamma_{\text{KN}}$ and $\epsilon(\gamma) \sim \gamma m_e c^2$ otherwise. Furthermore, the scattering rate, $\sigma_T c n_{\text{ph}} \epsilon_{\text{KN}}(\gamma/\gamma_{\text{KN}})$ [equation (5)], becomes a non-trivial, decreasing function of $\gamma$. This suppresses the emission efficiency and breaks the simple correspondence between the particle distribution power-law index and that of the IC emission spectrum. Let us examine how these effects manifest themselves in the IC(KN) and IC(KN)+PP panels (bottom left and bottom right, respectively) of Fig. 9. To begin with, the characteristic emission energies from particles at each of our Lorentz-factor scales are pushed closer together beyond the energy $\epsilon(\gamma_{\text{KN}})$ (because, beyond $\gamma_{\text{KN}}$, $\epsilon$ scales linearly with $\gamma$ instead of quadratically). As an example, even though the scales $\gamma_{\text{max}}$, $\gamma_{\text{rad, T}}$, and $\gamma_{\text{cool}}$ are all equally spaced on a logarithmic scale (because $\gamma_{\text{rad, T}}^2 = \gamma_{\text{max}}^2 \gamma_{\text{cool}}$), the corresponding photon energies, $\epsilon(\gamma_{\text{max}})$, $\epsilon(\gamma_{\text{rad, T}})$, and $\epsilon(\gamma_{\text{cool}})$, are not at all evenly spaced, with $\epsilon(\gamma_{\text{max}})$ and $\epsilon(\gamma_{\text{rad, T}})$ closer together than $\epsilon(\gamma_{\text{rad, T}})$ and $\epsilon(\gamma_{\text{cool}})$ [because the former two lie above the break energy $\epsilon(\gamma_{\text{KN}})$ while the latter lies below it]. Next, even though the power-law scalings of the IC(KN) and IC(KN)+PP particle distributions are both approximately $\gamma^{-2}$ or shallower (i.e. harder), the corresponding $\epsilon F_{\text{IC}}(\epsilon)$ spectra are both steeper (i.e. softer) than the Thomson-limit prediction, $\epsilon^{0.5}$, demonstrating the reduced radiative efficiency in the Klein–Nishina limit.

While the diminished Klein–Nishina cross-section produces a softer emission spectrum for a given particle distribution, it also yields a particle distribution that is harder in the first place (Section 4.3). These two effects somewhat cancel out, and, hence, not much change is observed in the spectral slope from the time-averaged IC(Th) spectrum to those yielded by the IC(KN) and IC(KN)+PP simulations (cf. Moderski et al. 2005). This is even despite the very different shape – the result of different cooling physics – in the particle distributions between these runs.

Finally, we note that the IC(KN) and IC(KN)+PP Compton emission spectra peak far above $\epsilon(\gamma_{\text{KN}})$ – deeply in the Klein–Nishina regime. As a result, most of the radiated energy in the IC(KN)+PP case (grey region in the lower right panel of Fig. 9) is emitted above pair-production threshold, where [unlike the IC(KN) run] it will be recaptured by the system as hot newborn pairs. The peak in the intrinsic emitted IC spectrum is therefore invisible to the observer, who sees only the indirect remnant of this radiation reprocessed to below-threshold energies. We elaborate the observable consequences of this effect in the next subsection.

## 4.5 Light curves and spectral variability

We now complement Section 4.4’s energy-resolved view of the emission from our simulations by discussing the timing of the radiative signatures. When viewed through the lens of timing, the
<figure><img src="image_9.png" /><figcaption>**Figure 9.** Time-averaged (over the same interval as in Figs 5–7) IC emission spectra for our four main runs. Instantaneous spectra used for averaging are equal to the box-integrated isotropic IC emissivity, $\epsilon F_{IC}(\epsilon; t) = \epsilon \int d^3x d\Omega j_{IC}(\epsilon, x, \Omega; t)$. Error envelopes display the 1$\sigma$ percentiles over the averaging interval in each energy bin. The normalization is arbitrary but equal across all panels. For the no rad. run, besides the normalization, the luminosity scale itself is arbitrary – not tied to the energetics of the simulation – and is chosen so that the spectral peak is at a similar level to those of the other runs. The characteristic photon energy, $\epsilon(\gamma)$, emitted by a particle of Lorentz factor $\gamma$ is defined piecewise continuously as $4\gamma^2\epsilon_{ph}(\gamma m_e c^2)$ for $\gamma$ less (greater) than $\gamma_{KN}$. The IC(KN)+PP spectrum is decomposed into contributions from originally present (orig.) and produced (prod.) particles. Klein–Nishina effects suppress the IC cross-section, hardening the particle energy distribution (Fig. 7), while softening the emitted spectrum. These competing effects result in IC(KN) spectral scalings that are not very different from that in the IC(Th) simulation. The grey region in the IC(KN)+PP panel indicates photon energies beyond pair-production threshold, $4\gamma_{KN}m_e c^2$. Photons emitted at these energies do not make it to the observer; they are absorbed inside the system to produce new pairs.</figcaption></figure>

differences among the various radiative regimes are accentuated, resulting in highly distinct observable signatures.

We begin by presenting light curves of each run's bolometric luminosity (instantaneous total escaping emitted power) in Fig. 10. For consistency with the flow of energy in the simulations, we only include that part of the luminosity permanently lost by the simulation. This means that, for the light curves of the no rad. and IC(Th) runs, we report the frequency-integrated emitted power as a function of time. In contrast, for the IC(KN)+PP run, we report only the portion of the emission spectrum below pair-production threshold: at photon energies $\epsilon < 4\gamma_{KN}m_e c^2$. By the same reasoning, we are obligated to include all photon energies for the IC(KN) simulation, reporting the total emitted power in that case also – otherwise, since pair production is artificially suppressed in that run, we would not fairly count the energy lost from the simulation. We have checked that this bookkeeping yields the same total energy radiated (integrals of the curves in Fig. 10) by each radiative simulation.

The light curves in Fig. 10 are clearly separated into two main groups: the non-radiative versus the radiative simulations. This dichotomy excellently illustrates a fundamental property of radiative reconnection: prompt emission. That is, in non-radiative reconnection, particles are first accelerated and, then, over much longer timescales than the duration of the reconnection process itself, radiate away their energy as potentially observable emission. In contrast, radiative reconnection features fundamentally prompt emission, where particles radiate their acquired energy on subdynamical timescales, causing radiation to participate in the reconnection dynamics rather than, as in the non-radiative case, passively trace energization that has already occurred. This is reflected in Fig. 10 in that the three radiative light curves broadly track the instantaneous electromagnetic dissipation of their simulations, rising as reconnection gets going and falling again once the reconnected flux saturates (cf. Fig. 5). The light curve of the non-radiative simulation, on the other hand, follows the time integral of the electromagnetic dissipation, growing with the cumulative dissipated energy and reconnected flux.
<figure><img src="image_10.png" /><figcaption>**Figure 9.** Time-averaged (over the same interval as in Figs 5–7) IC emission spectra for our four main runs. Instantaneous spectra used for averaging are equal to the box-integrated isotropic IC emissivity, $\epsilon F_{IC}(\epsilon; t) = \epsilon \int d^3x d\Omega j_{IC}(\epsilon, x, \Omega; t)$. Error envelopes display the 1$\sigma$ percentiles over the averaging interval in each energy bin. The normalization is arbitrary but equal across all panels. For the no rad. run, besides the normalization, the luminosity scale itself is arbitrary – not tied to the energetics of the simulation – and is chosen so that the spectral peak is at a similar level to those of the other runs. The characteristic photon energy, $\epsilon(\gamma)$, emitted by a particle of Lorentz factor $\gamma$ is defined piecewise continuously as $4\gamma^2\epsilon_{ph}(\gamma m_e c^2)$ for $\gamma$ less (greater) than $\gamma_{KN}$. The IC(KN)+PP spectrum is decomposed into contributions from originally present (orig.) and produced (prod.) particles. Klein–Nishina effects suppress the IC cross-section, hardening the particle energy distribution (Fig. 7), while softening the emitted spectrum. These competing effects result in IC(KN) spectral scalings that are not very different from that in the IC(Th) simulation. The grey region in the IC(KN)+PP panel indicates photon energies beyond pair-production threshold, $4\gamma_{KN}m_e c^2$. Photons emitted at these energies do not make it to the observer; they are absorbed inside the system to produce new pairs.</figcaption></figure>

differences among the various radiative regimes are accentuated, resulting in highly distinct observable signatures.

We begin by presenting light curves of each run's bolometric luminosity (instantaneous total escaping emitted power) in Fig. 10. For consistency with the flow of energy in the simulations, we only include that part of the luminosity permanently lost by the simulation. This means that, for the light curves of the no rad. and IC(Th) runs, we report the frequency-integrated emitted power as a function of time. In contrast, for the IC(KN)+PP run, we report only the portion of the emission spectrum below pair-production threshold: at photon energies $\epsilon < 4\gamma_{KN}m_e c^2$. By the same reasoning, we are obligated to include all photon energies for the IC(KN) simulation, reporting the total emitted power in that case also – otherwise, since pair production is artificially suppressed in that run, we would not fairly count the energy lost from the simulation. We have checked that this bookkeeping yields the same total energy radiated (integrals of the curves in Fig. 10) by each radiative simulation.

The light curves in Fig. 10 are clearly separated into two main groups: the non-radiative versus the radiative simulations. This dichotomy excellently illustrates a fundamental property of radiative reconnection: prompt emission. That is, in non-radiative reconnection, particles are first accelerated and, then, over much longer timescales than the duration of the reconnection process itself, radiate away their energy as potentially observable emission. In contrast, radiative reconnection features fundamentally prompt emission, where particles radiate their acquired energy on subdynamical timescales, causing radiation to participate in the reconnection dynamics rather than, as in the non-radiative case, passively trace energization that has already occurred. This is reflected in Fig. 10 in that the three radiative light curves broadly track the instantaneous electromagnetic dissipation of their simulations, rising as reconnection gets going and falling again once the reconnected flux saturates (cf. Fig. 5). The light curve of the non-radiative simulation, on the other hand, follows the time integral of the electromagnetic dissipation, growing with the cumulative dissipated energy and reconnected flux.
<figure><img src="image_10.png" /><figcaption>**Figure 10.** Light curves of the IC luminosity (instantaneous escaping radiated power) for each of our four main runs. The normalization for the non-radiative simulation is not fixed by dynamical radiative cooling and so is scaled to appear on the plot with the other simulations (cf. Fig. 9). Dashed portions of light curves are computed in post-processing (see the text).</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Figure 11.** Cross-correlation of the instantaneous rate of electromagnetic dissipation with the escaping IC luminosity (light curves of Fig. 10). Ambematic of the distinction between radiative and non-radiative reconnection, the non-radiative simulation's luminosity lags its dissipation (proportional to the cumulative dissipated electromagnetic energy), whereas the radiative simulations' luminosities are synchronized with dissipation (proportional to its instantaneous rate). The IC(KN) run features a slightly longer lag than the IC(Th) run owing to its lower radiative efficiency. The IC(KN)+PP run features a still somewhat longer lag with a skewed distribution toward the high-lag end. This stems from the time it takes for radiation injected above pair-production threshold to become reprocessed to below-threshold (and, hence, escaping) energies.</figcaption></figure>

To illustrate these remarks more thoroughly, we also supply Fig. 11, which shows the cross-correlation function, $\propto \int dt L_{\text{IC}}(t + \tau) P_{\text{diss}}(t)$, with $\tau$ the correlation lag, of each $L_{\text{IC}}(t)$ curve from Fig. 10 with the box-integrated electromagnetic dissipation, $P_{\text{diss}}(t) \equiv \int d^3x \mathbf{J}(t) \cdot \mathbf{E}(t)$. As discussed above and confirmed by this figure, there is a significant lag, $\tau > 0$, of the non-radiative IC luminosity behind $P_{\text{diss}}$. In contrast, all three radiative simulations lose their energy promptly, with peak lags close to zero.

There are, however, differences among the radiative runs. The IC(Th) luminosity exhibits the shortest variability time-scales in

Fig. 10, with small bumps atop its overall envelope corresponding to bursts of particle acceleration at plasmoid mergers. In contrast, because the cooling time-scales are slightly longer for the IC(KN) and IC(KN)+PP runs, particles retain enough energy in between these episodic events to smooth out the variations in the resulting light curves. As one might expect from these remarks, the IC(Th) simulation exhibits the smallest (exactly zero) peak lag from its electromagnetic dissipation to its emitted luminosity. In contrast, the IC(KN) and IC(KN)+PP runs have small but finite peak lag. Of these, the IC(KN)+PP case has a slightly longer lag and a more skewed cross-correlation distribution, with more power concentrated at longer lags. This is the result of the reprocessing of gamma-rays (near the peak of the emission spectrum in Fig. 9) to lower energies through gamma-ray radiation and pair production. That is, it takes time for power injected at the high-energy, above-threshold peak of the IC emission spectrum to be processed down to lower energies where it can escape the system.

When combined with spectral information, the lag induced by pair-production and gamma-ray absorption described above leaves a prominent imprint on the observable signatures of Klein–Nishina radiative reconnection – one that could potentially distinguish it from other radiative regimes. To demonstrate this, we present a side-by-side comparison of the luminosities and average photon energies, $\langle\epsilon\rangle \equiv \int d\epsilon \epsilon F_{\text{IC}}(\epsilon)/\int d\epsilon F_{\text{IC}}(\epsilon)$, viewed as a function of time for the IC(Th) and IC(KN)+PP runs in Fig. 12.

To fully appreciate the observational differences between these runs, we extrapolate their spectra beyond the end of each simulation. This is done under the assumption that the simulations are evolved to the point where no further exchange of energy occurs between the particles and the fields. Hence, the particle energy distributions can be passively Compton cooled in 1D energy space. Pair production can be included in this 1D evolution thanks to the homogeneous, isotropic, and static nature of the seed photon bath, which introduces no spatial or velocity-direction dependence. Specifically, the post-processing is implemented by sampling the particle and photon energy distributions near the end of the simulation (at a time-step after the reconnection active period; cf. Section 4.2) with a large number of quasi-particles and quasi-photons: like a PIC code, except that the quasi-particles and quasi-photons only have one coordinate – their energy – rather than six corresponding to their position and momentum vectors. Then, particle cooling is employed after the method described in Section 2.2.3, while pair production is calculated using the procedure from Section 2.2.2. Post-processing, in this way, the luminosity and average photon energy time-series allows us to extend them in Fig. 12 from just shy of $4L/c$ to more than $10L/c$ – a significant gain over the simulations themselves, which would otherwise need to be expensively evolved to more than double their actual duration to reach the same times.

Several facts, each of them accessible to observations, are immediately apparent from Fig. 12. First, the overall luminosity and average photon energy emitted from the IC(Th) simulation are tightly correlated. During the short time-scale variations of the light curve, sporadic magnetic reconnection acceleration yields an extended non-thermal distribution of particles (Fig. 7) with the resulting IC emission spectrum peaking near the high-energy cut-off (Fig. 9): that is, the acceleration mechanism produces a correlation between $\langle\epsilon\rangle$ and $L_{\text{IC}}$ on the rising side of each subpeak in their time-series. In between these reconnection energization episodes, particles emitting at the highest energies – the ones controlling both the overall luminosity and peak photon energy – are also the most rapidly cooled. They thus suddenly plummet to lower energies, inducing
<figure><img src="image_11.png" /><figcaption>**Figure 10.** Light curves of the IC luminosity (instantaneous escaping radiated power) for each of our four main runs. The normalization for the non-radiative simulation is not fixed by dynamical radiative cooling and so is scaled to appear on the plot with the other simulations (cf. Fig. 9). Dashed portions of light curves are computed in post-processing (see the text).</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Figure 11.** Cross-correlation of the instantaneous rate of electromagnetic dissipation with the escaping IC luminosity (light curves of Fig. 10). Ambematic of the distinction between radiative and non-radiative reconnection, the non-radiative simulation's luminosity lags its dissipation (proportional to the cumulative dissipated electromagnetic energy), whereas the radiative simulations' luminosities are synchronized with dissipation (proportional to its instantaneous rate). The IC(KN) run features a slightly longer lag than the IC(Th) run owing to its lower radiative efficiency. The IC(KN)+PP run features a still somewhat longer lag with a skewed distribution toward the high-lag end. This stems from the time it takes for radiation injected above pair-production threshold to become reprocessed to below-threshold (and, hence, escaping) energies.</figcaption></figure>

To illustrate these remarks more thoroughly, we also supply Fig. 11, which shows the cross-correlation function, $\propto \int dt L_{\text{IC}}(t + \tau) P_{\text{diss}}(t)$, with $\tau$ the correlation lag, of each $L_{\text{IC}}(t)$ curve from Fig. 10 with the box-integrated electromagnetic dissipation, $P_{\text{diss}}(t) \equiv \int d^3x \mathbf{J}(t) \cdot \mathbf{E}(t)$. As discussed above and confirmed by this figure, there is a significant lag, $\tau > 0$, of the non-radiative IC luminosity behind $P_{\text{diss}}$. In contrast, all three radiative simulations lose their energy promptly, with peak lags close to zero.

There are, however, differences among the radiative runs. The IC(Th) luminosity exhibits the shortest variability time-scales in

Fig. 10, with small bumps atop its overall envelope corresponding to bursts of particle acceleration at plasmoid mergers. In contrast, because the cooling time-scales are slightly longer for the IC(KN) and IC(KN)+PP runs, particles retain enough energy in between these episodic events to smooth out the variations in the resulting light curves. As one might expect from these remarks, the IC(Th) simulation exhibits the smallest (exactly zero) peak lag from its electromagnetic dissipation to its emitted luminosity. In contrast, the IC(KN) and IC(KN)+PP runs have small but finite peak lag. Of these, the IC(KN)+PP case has a slightly longer lag and a more skewed cross-correlation distribution, with more power concentrated at longer lags. This is the result of the reprocessing of gamma-rays (near the peak of the emission spectrum in Fig. 9) to lower energies through gamma-ray radiation and pair production. That is, it takes time for power injected at the high-energy, above-threshold peak of the IC emission spectrum to be processed down to lower energies where it can escape the system.

When combined with spectral information, the lag induced by pair-production and gamma-ray absorption described above leaves a prominent imprint on the observable signatures of Klein–Nishina radiative reconnection – one that could potentially distinguish it from other radiative regimes. To demonstrate this, we present a side-by-side comparison of the luminosities and average photon energies, $\langle\epsilon\rangle \equiv \int d\epsilon \epsilon F_{\text{IC}}(\epsilon)/\int d\epsilon F_{\text{IC}}(\epsilon)$, viewed as a function of time for the IC(Th) and IC(KN)+PP runs in Fig. 12.

To fully appreciate the observational differences between these runs, we extrapolate their spectra beyond the end of each simulation. This is done under the assumption that the simulations are evolved to the point where no further exchange of energy occurs between the particles and the fields. Hence, the particle energy distributions can be passively Compton cooled in 1D energy space. Pair production can be included in this 1D evolution thanks to the homogeneous, isotropic, and static nature of the seed photon bath, which introduces no spatial or velocity-direction dependence. Specifically, the post-processing is implemented by sampling the particle and photon energy distributions near the end of the simulation (at a time-step after the reconnection active period; cf. Section 4.2) with a large number of quasi-particles and quasi-photons: like a PIC code, except that the quasi-particles and quasi-photons only have one coordinate – their energy – rather than six corresponding to their position and momentum vectors. Then, particle cooling is employed after the method described in Section 2.2.3, while pair production is calculated using the procedure from Section 2.2.2. Post-processing, in this way, the luminosity and average photon energy time-series allows us to extend them in Fig. 12 from just shy of $4L/c$ to more than $10L/c$ – a significant gain over the simulations themselves, which would otherwise need to be expensively evolved to more than double their actual duration to reach the same times.

Several facts, each of them accessible to observations, are immediately apparent from Fig. 12. First, the overall luminosity and average photon energy emitted from the IC(Th) simulation are tightly correlated. During the short time-scale variations of the light curve, sporadic magnetic reconnection acceleration yields an extended non-thermal distribution of particles (Fig. 7) with the resulting IC emission spectrum peaking near the high-energy cut-off (Fig. 9): that is, the acceleration mechanism produces a correlation between $\langle\epsilon\rangle$ and $L_{\text{IC}}$ on the rising side of each subpeak in their time-series. In between these reconnection energization episodes, particles emitting at the highest energies – the ones controlling both the overall luminosity and peak photon energy – are also the most rapidly cooled. They thus suddenly plummet to lower energies, inducing
<figure><img src="image_12.png" /><figcaption>**Figure 12.** Light curves (reproduced from Fig. 10) of escaping IC luminosity and overall (computed using the escaping spectrum) average photon energy for the IC(Th) and IC(KN)+PP simulations. Both quantities are extrapolated beyond the simulation (from time $t_0$ indicated by the vertical black bar) via post-processing (see the text). After $t_0$ (but not before), the luminosities are compensated by $(t/t_0)^6$ – an empirical choice. The time-series of mean photon energy are never compensated. Several robust trends are evident that could distinguish these two types of reconnection in the context of reconnection-powered flares. Thomson radiative reconnection features a tight correlation between IC luminosity and mean photon energy. In contrast, the Klein–Nishina reconnection mean photon energy is ~0.3$\gamma_{KN}m_e c^2$ independently of luminosity. These patterns hold both in the reconnection active phase and in the decaying phase when energization from the fields is either slowing down (as at late simulation times) or completely shut off (as in the post-processing zone of the plot). In the limit of no further energy injection, the $t^{-6}$ brightness decay visible in the Klein–Nishina run sharply contrasts the much more precipitous decline in the Thomson case.</figcaption></figure>

simultaneous drops in both $\langle\epsilon\rangle$ and $L_{IC}$: that is, the radiative cooling mediates the correlation between $\langle\epsilon\rangle$ and $L_{IC}$ on the falling side of each peak in their time-series.

A second observationally pertinent property of the IC(Th) curves in Fig. 12 is that, once energization from the electromagnetic fields is shut off, both the IC photon energy and luminosity drop precipitously. In fact, even though the luminosity time-series are compensated by $(t/t_0)^6$ rightward of the transition time, $t_0$, to the post-processing stage, this does little to stem the fall of the IC(Th) luminosity. Flares from highly radiative Thomson reconnection are thus characterized by *tight correlation between the observed luminosity and photon energy plus a rapid falling phase where both plummet together*.

Let us now examine how these qualities compare to those of the IC(KN)+PP run. Most importantly, the tight correlation between luminosity and photon energy is broken. Instead, irrespective of the instantaneous IC brightness, the average photon energy remains rock steady, persisting near 0.3$\gamma_{KN}m_e c^2$ – even in the decaying phase of the light curve after electromagnetic energization has ceased. In the following, we argue that this average photon energy is first set during active periods of reconnection-powered particle acceleration and subsequently reinforced, when such acceleration is inactive, by the Klein–Nishina radiative physics, explaining its persistence.

During active acceleration episodes, reconnection produces a hard distribution of radiating particles with corresponding upward-sloped $\epsilon F_{IC}(\epsilon)$ emission spectrum peaking well above pair-production threshold (Fig. 9). Most of the initially emitted energy is, thus, veiled by pair production, and the peak of the apparent/observed spectrum lies instead just before the absorption-induced cut-off: that is, at energies ~$\gamma_{KN}m_e c^2$, as seen in the IC(KN)+PP run. Subsequently, in between reconnection-powered acceleration events, the radiative physics takes over in determining the mean photon energy. As seen during these periods in the Thomson regime, the IC spectrum softens because rapidly cooling particles radiating at the spectral peak energy, ~$\gamma_{rad,T}^2 \epsilon_{ph}$ (Fig. 9), cannot be replenished by electromagnetic energization. The key difference in the IC(KN)+PP

case is that particles emitting at the observed peak energy, ~$\gamma_{KN}m_e c^2$, can still be partially replenished by IC cooling and pair production, which actively reprocess radiation originally emitted at higher, absorbed photon energies down to the observed band. This stabilizes the observed spectrum, even in the falling phase of the light curve when electromagnetic energization is completely absent. Hence, the photon energy 0.3$\gamma_{KN}m_e c^2$ owes its luminosity-independent stability to the fact that the coupling of reconnection-powered NTPA to Klein–Nishina radiative physics results in the same natural photon energy scale as that produced by the radiative physics alone.

We note that radiative reprocessing of initially above-threshold photons also leads, after reconnection has concluded, to the self-similar $t^{-6}$ power-law brightness decay shown in Fig. 12. Although this is much slower than in the case of Thomson radiative cooling, it is still probably too abrupt for gamma-ray instruments to resolve. Thus, what we would like to stress as the main difference between gamma-ray flares powered by Thomson-radiative and Klein–Nishina reconnection is that the latter are characterized by a *constant mean observed photon energy, irrespective of brightness*.

We conduct post-processing experiments in Fig. 13 that suggest that the presented properties of the IC(KN)+PP light curve and spectrum in the absence of particle acceleration are universal. In these experiments, we evolve different initial power-law distributions of particles – $dN/d\gamma \propto \gamma^{-p}$ for $p = 1, 2$, and 3, and $\gamma_{cool} < \gamma < \gamma_{max}$ (dashed lines in Fig. 13, left panel) – solely under the influence of Klein–Nishina emission and pair production (as in the post-processing phase of Fig. 12 except that, there, the initial particle distribution is taken from time $t_0$ of our PIC simulation). We find that, irrespective of the initial power-law slope, the particle distribution always relaxes, in the Thomson regime, $\gamma \lesssim \gamma_{KN}$, to a $\gamma^{-2}$ scaling (solid lines in Fig. 13, left panel).<fnref n="6" /> This corresponds,

<footnote n="6">⁶The slope, $\gamma^{-2}$, can be calculated by considering the trickle of particles from $\gamma > \gamma_{KN}$ as monochromatic particle injection at $\gamma = \gamma_{KN}$ (Mehlhaff et al. 2021).</footnote>
<figure><img src="image_13.png" /><figcaption>**Figure 12.** Light curves (reproduced from Fig. 10) of escaping IC luminosity and overall (computed using the escaping spectrum) average photon energy for the IC(Th) and IC(KN)+PP simulations. Both quantities are extrapolated beyond the simulation (from time $t_0$ indicated by the vertical black bar) via post-processing (see the text). After $t_0$ (but not before), the luminosities are compensated by $(t/t_0)^6$ – an empirical choice. The time-series of mean photon energy are never compensated. Several robust trends are evident that could distinguish these two types of reconnection in the context of reconnection-powered flares. Thomson radiative reconnection features a tight correlation between IC luminosity and mean photon energy. In contrast, the Klein–Nishina reconnection mean photon energy is ~0.3$\gamma_{KN}m_e c^2$ independently of luminosity. These patterns hold both in the reconnection active phase and in the decaying phase when energization from the fields is either slowing down (as at late simulation times) or completely shut off (as in the post-processing zone of the plot). In the limit of no further energy injection, the $t^{-6}$ brightness decay visible in the Klein–Nishina run sharply contrasts the much more precipitous decline in the Thomson case.</figcaption></figure>

simultaneous drops in both $\langle\epsilon\rangle$ and $L_{IC}$: that is, the radiative cooling mediates the correlation between $\langle\epsilon\rangle$ and $L_{IC}$ on the falling side of each peak in their time-series.

A second observationally pertinent property of the IC(Th) curves in Fig. 12 is that, once energization from the electromagnetic fields is shut off, both the IC photon energy and luminosity drop precipitously. In fact, even though the luminosity time-series are compensated by $(t/t_0)^6$ rightward of the transition time, $t_0$, to the post-processing stage, this does little to stem the fall of the IC(Th) luminosity. Flares from highly radiative Thomson reconnection are thus characterized by *tight correlation between the observed luminosity and photon energy plus a rapid falling phase where both plummet together*.

Let us now examine how these qualities compare to those of the IC(KN)+PP run. Most importantly, the tight correlation between luminosity and photon energy is broken. Instead, irrespective of the instantaneous IC brightness, the average photon energy remains rock steady, persisting near 0.3$\gamma_{KN}m_e c^2$ – even in the decaying phase of the light curve after electromagnetic energization has ceased. In the following, we argue that this average photon energy is first set during active periods of reconnection-powered particle acceleration and subsequently reinforced, when such acceleration is inactive, by the Klein–Nishina radiative physics, explaining its persistence.

During active acceleration episodes, reconnection produces a hard distribution of radiating particles with corresponding upward-sloped $\epsilon F_{IC}(\epsilon)$ emission spectrum peaking well above pair-production threshold (Fig. 9). Most of the initially emitted energy is, thus, veiled by pair production, and the peak of the apparent/observed spectrum lies instead just before the absorption-induced cut-off: that is, at energies ~$\gamma_{KN}m_e c^2$, as seen in the IC(KN)+PP run. Subsequently, in between reconnection-powered acceleration events, the radiative physics takes over in determining the mean photon energy. As seen during these periods in the Thomson regime, the IC spectrum softens because rapidly cooling particles radiating at the spectral peak energy, ~$\gamma_{rad,T}^2 \epsilon_{ph}$ (Fig. 9), cannot be replenished by electromagnetic energization. The key difference in the IC(KN)+PP

case is that particles emitting at the observed peak energy, ~$\gamma_{KN}m_e c^2$, can still be partially replenished by IC cooling and pair production, which actively reprocess radiation originally emitted at higher, absorbed photon energies down to the observed band. This stabilizes the observed spectrum, even in the falling phase of the light curve when electromagnetic energization is completely absent. Hence, the photon energy 0.3$\gamma_{KN}m_e c^2$ owes its luminosity-independent stability to the fact that the coupling of reconnection-powered NTPA to Klein–Nishina radiative physics results in the same natural photon energy scale as that produced by the radiative physics alone.

We note that radiative reprocessing of initially above-threshold photons also leads, after reconnection has concluded, to the self-similar $t^{-6}$ power-law brightness decay shown in Fig. 12. Although this is much slower than in the case of Thomson radiative cooling, it is still probably too abrupt for gamma-ray instruments to resolve. Thus, what we would like to stress as the main difference between gamma-ray flares powered by Thomson-radiative and Klein–Nishina reconnection is that the latter are characterized by a *constant mean observed photon energy, irrespective of brightness*.

We conduct post-processing experiments in Fig. 13 that suggest that the presented properties of the IC(KN)+PP light curve and spectrum in the absence of particle acceleration are universal. In these experiments, we evolve different initial power-law distributions of particles – $dN/d\gamma \propto \gamma^{-p}$ for $p = 1, 2$, and 3, and $\gamma_{cool} < \gamma < \gamma_{max}$ (dashed lines in Fig. 13, left panel) – solely under the influence of Klein–Nishina emission and pair production (as in the post-processing phase of Fig. 12 except that, there, the initial particle distribution is taken from time $t_0$ of our PIC simulation). We find that, irrespective of the initial power-law slope, the particle distribution always relaxes, in the Thomson regime, $\gamma \lesssim \gamma_{KN}$, to a $\gamma^{-2}$ scaling (solid lines in Fig. 13, left panel).<fnref n="6" /> This corresponds,

<footnote n="6">⁶The slope, $\gamma^{-2}$, can be calculated by considering the trickle of particles from $\gamma > \gamma_{KN}$ as monochromatic particle injection at $\gamma = \gamma_{KN}$ (Mehlhaff et al. 2021).</footnote>
<figure><img src="image_13.png" /><figcaption>**Figure 13.** Left: numerical experiment wherein different initial power-law particle distributions, $dN/d\gamma \propto \gamma^{-p}$ for $p = 1, 2$, and 3, are passively cooled under IC emission and pair production. All parameters are identical to the IC(KN)+PP simulation. Integration is done from time $t_0$ (dashed curves) to $t_0 + 4.1L/c$ (solid curves). For reference, we present, in blue, the passive cooling of the particle distribution yielded by the IC(KN)+PP simulation (as in Fig. 12). Right: IC emission spectra corresponding to the particle distributions at left. A grey region indicates unobserved (above pair-production threshold) photon energies. After an initial-condition-dependent relaxation period, all distributions converge to an initial-condition-independent universal shape. At observable energies – that is, in the Thomson regime, $\epsilon \lesssim \gamma_{KN}m_e c^2 (\gamma \lesssim \gamma_{KN})$ – this universal shape features an approximate $\gamma^{-2}$ power law, corresponding to $\epsilon F_{IC}(\epsilon) \propto \epsilon^{0.5}$. (The apparent upturn in the particle distributions at low energies is an artefact of short integration times – it results from a bulk population of cooled particles that would gradually become still colder with time, extending the $\gamma^{-2}$ scaling to even lower energies.).</figcaption></figure>

for the Thomson part, $\epsilon \lesssim \gamma_{KN}m_e c^2$, of the emission spectrum, to $\epsilon F_{IC}(\epsilon) \propto \epsilon^{-(2-3)/2} = \epsilon^{1/2}$ (see Section 4.4) – a rising spectrum that continues almost up to the pair-production threshold energy $\epsilon = 4\gamma_{KN}m_e c^2$ (solid lines in Fig. 13, right panel). These universal shapes, once reached, are maintained by the particle distribution and emission spectrum as they fall off, resulting in a self-similar $t^{-6}$ luminosity decay law and constant mean photon energy, $\langle \epsilon \rangle \sim 0.3\gamma_{KN}m_e c^2$: the same as Fig. 12.

The above exercise enables us to reason about the observable signatures of Klein–Nishina reconnection in regimes, not probed by our simulations, where reconnection-powered NTPA is known to yield a steeper particle energy distribution. This occurs, for example, in the presence of a strong guide field (Werner & Uzdensky 2017) or in the transrelativistic regime of electron–proton plasmas (Werner et al. 2018). In such cases, we speculate that one would observe an initially steep emission spectrum corresponding to intrinsic reconnection-powered NTPA, followed, in the passive cooling phase, by a transition – in fact, a hardening! – to the identified universal shape as the flux decays (similar to the $dN/d\gamma \propto \gamma^{-3}$ initial condition in Fig. 13). This might be difficult to observe, however, as the spectrum may dim too much before relaxing to the expected shape, starving gamma-ray instruments of a sufficient number of photons to reconstruct it (cf. the large gap between the initial and final $dN/d\gamma \propto \gamma^{-3}$ curves in Fig. 13).

We remark that the universal shape of the passively evolved (i.e. only via IC losses and pair production but without particle acceleration) particle energy distribution – and the consequent stability of the emission spectrum when energization subsides – depends only on Klein–Nishina effects that are independent of the number of spatial dimensions. This result should therefore generalize to 3D reconnection. However, given the differing acceleration physics in 3D (e.g., as discussed in Section 4.3, particles escaping from the reconnection downstream and subsequently being re-accelerated, Zhang et al. 2021; Chernoglazov et al. 2023; Zhang et al. 2023), the stability of the spectrum on the rising part of the light curve may change. Ultimately, the 3D behaviour needs to be verified by full-fledged 3D simulations.

In this section, we have seen how the different regimes of radiative cooling treated by this study – particularly their influence on particle acceleration and the resulting IC emission spectra – give rise to highly distinct temporal radiative signatures. At the most coarse-grained level, emission from radiative reconnection tracks the electromagnetic dissipation in real time, whereas emission from non-radiative reconnection traces only energization that has occurred in the past. Focusing on more specific observable differences in the context of reconnection-powered flares, the light curve from Thomson radiative reconnection is highly correlated with the observed average photon energy (provided one observes near the spectral peak at $\sim \gamma_{\text{rad,T}}^2 \epsilon_{\text{ph}}$) and features an abrupt decay phase where both drop simultaneously. Klein–Nishina radiative reconnection with pair production, however, yields an exactly opposite trend, with no correlation between the average photon energy, $\sim 0.3\gamma_{KN}m_e c^2$, and the overall brightness, including in the (slower than in the Thomson regime, but still relatively fast) decay phase. These findings can be directly compared with, and tested by, observations of gamma-ray flares from suitable astrophysical systems (see Section 6). Finally, though this section is predicated on simulations with periodic boundaries and, hence, well-defined beginning and end phases of reconnection, the main results are likely to be more general. In particular, the same arguments used to interpret the relative spectral stability of Klein–Nishina reconnection compared to the Thomson case still hold in steady-state reconnection. Whether such effects would be observable in a steady state, however, depends on: (1) whether sufficiently large statistical fluctuations about the time average would occur to yield detectable variability in the observed radiation; and, (2) whether such variability would be distinguishable from that owing to quasi-static changes to the large-scale reconnection layer, as may be caused by macroscopic dynamics of the astrophysical system (cf. Section 6.4).

## 4.6 Newborn pair energy budget and particle count

In the preceding parts of Section 4, we explored consequences of Klein–Nishina and pair-production physics on magnetic reconnec-
<figure><img src="image_14.png" /><figcaption>**Figure 13.** Left: numerical experiment wherein different initial power-law particle distributions, $dN/d\gamma \propto \gamma^{-p}$ for $p = 1, 2$, and 3, are passively cooled under IC emission and pair production. All parameters are identical to the IC(KN)+PP simulation. Integration is done from time $t_0$ (dashed curves) to $t_0 + 4.1L/c$ (solid curves). For reference, we present, in blue, the passive cooling of the particle distribution yielded by the IC(KN)+PP simulation (as in Fig. 12). Right: IC emission spectra corresponding to the particle distributions at left. A grey region indicates unobserved (above pair-production threshold) photon energies. After an initial-condition-dependent relaxation period, all distributions converge to an initial-condition-independent universal shape. At observable energies – that is, in the Thomson regime, $\epsilon \lesssim \gamma_{KN}m_e c^2 (\gamma \lesssim \gamma_{KN})$ – this universal shape features an approximate $\gamma^{-2}$ power law, corresponding to $\epsilon F_{IC}(\epsilon) \propto \epsilon^{0.5}$. (The apparent upturn in the particle distributions at low energies is an artefact of short integration times – it results from a bulk population of cooled particles that would gradually become still colder with time, extending the $\gamma^{-2}$ scaling to even lower energies.).</figcaption></figure>

for the Thomson part, $\epsilon \lesssim \gamma_{KN}m_e c^2$, of the emission spectrum, to $\epsilon F_{IC}(\epsilon) \propto \epsilon^{-(2-3)/2} = \epsilon^{1/2}$ (see Section 4.4) – a rising spectrum that continues almost up to the pair-production threshold energy $\epsilon = 4\gamma_{KN}m_e c^2$ (solid lines in Fig. 13, right panel). These universal shapes, once reached, are maintained by the particle distribution and emission spectrum as they fall off, resulting in a self-similar $t^{-6}$ luminosity decay law and constant mean photon energy, $\langle\epsilon\rangle \sim 0.3\gamma_{KN}m_e c^2$: the same as Fig. 12.

The above exercise enables us to reason about the observable signatures of Klein–Nishina reconnection in regimes, not probed by our simulations, where reconnection-powered NTPA is known to yield a steeper particle energy distribution. This occurs, for example, in the presence of a strong guide field (Werner & Uzdensky 2017) or in the transrelativistic regime of electron–proton plasmas (Werner et al. 2018). In such cases, we speculate that one would observe an initially steep emission spectrum corresponding to intrinsic reconnection-powered NTPA, followed, in the passive cooling phase, by a transition – in fact, a hardening! – to the identified universal shape as the flux decays (similar to the $dN/d\gamma \propto \gamma^{-3}$ initial condition in Fig. 13). This might be difficult to observe, however, as the spectrum may dim too much before relaxing to the expected shape, starving gamma-ray instruments of a sufficient number of photons to reconstruct it (cf. the large gap between the initial and final $dN/d\gamma \propto \gamma^{-3}$ curves in Fig. 13).

We remark that the universal shape of the passively evolved (i.e. only via IC losses and pair production but without particle acceleration) particle energy distribution – and the consequent stability of the emission spectrum when energization subsides – depends only on Klein–Nishina effects that are independent of the number of spatial dimensions. This result should therefore generalize to 3D reconnection. However, given the differing acceleration physics in 3D (e.g., as discussed in Section 4.3, particles escaping from the reconnection downstream and subsequently being re-accelerated, Zhang et al. 2021; Chernoglazov et al. 2023; Zhang et al. 2023), the stability of the spectrum on the rising part of the light curve may change. Ultimately, the 3D behaviour needs to be verified by full-fledged 3D simulations.

In this section, we have seen how the different regimes of radiative cooling treated by this study – particularly their influence on particle acceleration and the resulting IC emission spectra – give rise to highly distinct temporal radiative signatures. At the most coarse-grained level, emission from radiative reconnection tracks the electromagnetic dissipation in real time, whereas emission from non-radiative reconnection traces only energization that has occurred in the past. Focusing on more specific observable differences in the context of reconnection-powered flares, the light curve from Thomson radiative reconnection is highly correlated with the observed average photon energy (provided one observes near the spectral peak at $\sim \gamma_{\text{rad,T}}^2 \epsilon_{\text{ph}}$) and features an abrupt decay phase where both drop simultaneously. Klein–Nishina radiative reconnection with pair production, however, yields an exactly opposite trend, with no correlation between the average photon energy, $\sim 0.3\gamma_{KN}m_e c^2$, and the overall brightness, including in the (slower than in the Thomson regime, but still relatively fast) decay phase. These findings can be directly compared with, and tested by, observations of gamma-ray flares from suitable astrophysical systems (see Section 6). Finally, though this section is predicated on simulations with periodic boundaries and, hence, well-defined beginning and end phases of reconnection, the main results are likely to be more general. In particular, the same arguments used to interpret the relative spectral stability of Klein–Nishina reconnection compared to the Thomson case still hold in steady-state reconnection. Whether such effects would be observable in a steady state, however, depends on: (1) whether sufficiently large statistical fluctuations about the time average would occur to yield detectable variability in the observed radiation; and, (2) whether such variability would be distinguishable from that owing to quasi-static changes to the large-scale reconnection layer, as may be caused by macroscopic dynamics of the astrophysical system (cf. Section 6.4).

## 4.6 Newborn pair energy budget and particle count

In the preceding parts of Section 4, we explored consequences of Klein–Nishina and pair-production physics on magnetic reconnec-
<figure><img src="image_14.png" /><figcaption>**Figure 14.** Left-hand column: snapshots of the relative contributions of original (orig.) and produced (prod.) pairs to (on the top half of each panel) the local particle number density, $(n_{\text{orig.}} - n_{\text{prod.}})/(n_{\text{orig.}} + n_{\text{prod.}})$, and (on the bottom half of each panel) the local particle energy density, $(U_{\text{orig.}} - U_{\text{prod.}})/(U_{\text{orig.}} + U_{\text{prod.}})$. Right-hand column: snapshots of the local cold, $\sigma_c/\sigma_{c,0}$ (top halves of panels), and hot, $\sigma_h/\sigma_{h,0}$ (bottom halves of panels), magnetizations. Contours display values of key quantities along, or just upstream of, the reconnection separatrix. For the $U_{\text{prod.}}$ contour, $U_{B0} = B_0^2/8\pi$. In agreement with Mehlhaff et al. (2021), newborn pairs are everywhere subdominant in terms of their number density, but compete with the original particles for energetic dominance. Moreover, the newborn particles create an energetically dense coat around the reconnection layer where they load the local hot magnetization but not the cold magnetization.</figcaption></figure>

tion, using as control cases the Thomson-radiative and non-radiative regimes. We presented first what is similar to the latter two cases – for example, the overall spatial dynamics and the reconnection rate – and then discussed the main distinctions, culminating with the very different observable signatures of the various radiative regimes. We now go one step farther, leaving behind our control cases in order to address issues that only exist in the context of Klein–Nishina reconnection with pair production. In particular, we comment on the newborn pairs’ contribution to the reconnection system’s energy and particle number budgets, which is presented graphically in Fig. 14.

The left-hand column of Fig. 14 illustrates the local differences in number densities, $(n_{\text{orig.}} - n_{\text{prod.}})/(n_{\text{orig.}} + n_{\text{prod.}})$, and energy densities, $(U_{\text{orig.}} - U_{\text{prod.}})/(U_{\text{orig.}} + U_{\text{prod.}})$, between originally present (orig.) and produced (prod.) particles. In the number density panels,

the newborn particles remain subdominant across time and space. However, once reconnection ignites above-threshold gamma-ray emission and pair production, the produced particles build up an energetically dense coat around the reconnection layer. Inside this coat, the newborn pairs completely dominate the energy density of the reconnection upstream region; in the reconnected flux region, while not as clearly dominant, they still vie for control of the energy budget. This demonstrates the main peculiarity of Klein–Nishina reconnection – previously anticipated by Mehlhaff et al. (2021) – that, for a broad range of parameters (made more precise below), the produced particles are hot but tenuous, competing with the original particles for energetic dominance of the system despite being much fewer in number. A parallel view of these effects is presented in the right-hand column of Fig. 14, where we plot the local cold magnetization $\sigma_c \equiv (B_x^2 + B_y^2)/[4\pi(n_{\text{orig.}} + n_{\text{prod.}})m_e c^2]$ and the hot
<figure><img src="image_15.png" /><figcaption>**Figure 14.** Left-hand column: snapshots of the relative contributions of original (orig.) and produced (prod.) pairs to (on the top half of each panel) the local particle number density, $(n_{\text{orig.}} - n_{\text{prod.}})/(n_{\text{orig.}} + n_{\text{prod.}})$, and (on the bottom half of each panel) the local particle energy density, $(U_{\text{orig.}} - U_{\text{prod.}})/(U_{\text{orig.}} + U_{\text{prod.}})$. Right-hand column: snapshots of the local cold, $\sigma_c/\sigma_{c,0}$ (top halves of panels), and hot, $\sigma_h/\sigma_{h,0}$ (bottom halves of panels), magnetizations. Contours display values of key quantities along, or just upstream of, the reconnection separatrix. For the $U_{\text{prod.}}$ contour, $U_{B0} = B_0^2/8\pi$. In agreement with Mehlhaff et al. (2021), newborn pairs are everywhere subdominant in terms of their number density, but compete with the original particles for energetic dominance. Moreover, the newborn particles create an energetically dense coat around the reconnection layer where they load the local hot magnetization but not the cold magnetization.</figcaption></figure>

tion, using as control cases the Thomson-radiative and non-radiative regimes. We presented first what is similar to the latter two cases – for example, the overall spatial dynamics and the reconnection rate – and then discussed the main distinctions, culminating with the very different observable signatures of the various radiative regimes. We now go one step farther, leaving behind our control cases in order to address issues that only exist in the context of Klein–Nishina reconnection with pair production. In particular, we comment on the newborn pairs’ contribution to the reconnection system’s energy and particle number budgets, which is presented graphically in Fig. 14.

The left-hand column of Fig. 14 illustrates the local differences in number densities, $(n_{\text{orig.}} - n_{\text{prod.}})/(n_{\text{orig.}} + n_{\text{prod.}})$, and energy densities, $(U_{\text{orig.}} - U_{\text{prod.}})/(U_{\text{orig.}} + U_{\text{prod.}})$, between originally present (orig.) and produced (prod.) particles. In the number density panels,

the newborn particles remain subdominant across time and space. However, once reconnection ignites above-threshold gamma-ray emission and pair production, the produced particles build up an energetically dense coat around the reconnection layer. Inside this coat, the newborn pairs completely dominate the energy density of the reconnection upstream region; in the reconnected flux region, while not as clearly dominant, they still vie for control of the energy budget. This demonstrates the main peculiarity of Klein–Nishina reconnection – previously anticipated by Mehlhaff et al. (2021) – that, for a broad range of parameters (made more precise below), the produced particles are hot but tenuous, competing with the original particles for energetic dominance of the system despite being much fewer in number. A parallel view of these effects is presented in the right-hand column of Fig. 14, where we plot the local cold magnetization $\sigma_c \equiv (B_x^2 + B_y^2)/[4\pi(n_{\text{orig.}} + n_{\text{prod.}})m_e c^2]$ and the hot
magnetization $\sigma_h \equiv (B_x^2 + B_y^2)/4\pi w$ [note: the enthalpy density, for our ultrarelativistic particle distributions, is $w = 4(U_{\text{orig.}} + U_{\text{prod.}})/3$; see equation (32) and surrounding discussion]. Because the newborn pairs load the upstream plasma energy density without appreciably changing the number density, they reduce $\sigma_h$ from its native/initial value, $\sigma_{h,0}$, while leaving $\sigma_c \simeq \sigma_{c,0}$ essentially untouched.

Even though the newborn pairs are everywhere less numerous than those originally present, they discernibly contribute to the plasma density [e.g. by changing the hue of the $(n_{\text{orig.}} - n_{\text{prod.}})/(n_{\text{orig.}} + n_{\text{prod.}})$ spatial maps] close to and inside of the reconnection separatrix. Mehlhaff et al. (2021) forecasted that the newborn pairs should begin to contribute non-negligibly to the plasma number density passing across the separatrix when $\sigma_{c,0}$ exceeds $\gamma_{\text{KN}}$ by more than a factor of 10 or so. This is consistent with our findings on the Klein–Nishina reconnection pair yield presented in the next section, which suggest that the IC(KN)+PP run ($\sigma_{c,0} \simeq 20\gamma_{\text{KN}}$) is indeed beginning to border on a regime where the *in situ* produced pairs contribute more significantly to the particle count. For the rest of this section, however, we focus on the present regime where the newborn particles are energetically dense but numerically few.

As discussed by Mehlhaff et al. (2021), the energy density injected into the upstream region via the hot newborn pairs leads to a *pair-loaded magnetization*, $\sigma_h^*$, to which the system tends to regulate itself (provided the initial magnetization, $\sigma_{h,0}$, exceeds $\sigma_h^*$). However, if the pair loading is efficient enough, the system may never actually attain a steady state with $\sigma_h = \sigma_h^*$. It will instead flood the upstream energy density via pair production, overshooting to a lower magnetization, $\sigma_h < \sigma_h^*$, and quenching subsequent NTPA until the upstream plasma is vacated and a high $\sigma_h > \sigma_h^*$ is restored. This would restart the process, leading to a limit cycle: the system would indefinitely ricochet between a high and a low magnetization on either side of $\sigma_h^*$ throughout the duration of reconnection.

Mehlhaff et al. (2021) quantified the dependence of the long-term fate of the system – whether it smoothly regulates to, or violently oscillates about, $\sigma_h^*$ – on the efficiency of upstream pair loading. We measure this efficiency from our simulations in Appendix B, finding that it is too low, based on the analysis of Mehlhaff et al. (2021), to trigger $\sigma_h$-mediated limit cycles. We therefore conclude that *limit cycles mediated by pair-loading of $\sigma_h$ are unlikely in astrophysical Klein–Nishina reconnection*.

In this section, we have shown that our simulations probe a Klein–Nishina reconnection regime where, as previously predicted by Mehlhaff et al. (2021), the newborn upstream pairs are hot but tenuous, loading the upstream plasma pressure, energy density, and hot magnetization, but not the upstream plasma density or cold magnetization. The simulations further provide measurements (Appendix B) that disfavour $\sigma_h$-mediated limit cycles. However, if a regime existed featuring a large multiplicity of newborn pairs, the possibility of $\sigma_c$-mediated limit cycles would remain an open question. It is the overall pair yield of Klein–Nishina reconnection, including the potential existence of such a copious pair-production regime, to which we now turn.

## 5 PAIR YIELD

We have already seen how the pairs produced in Klein–Nishina reconnection contribute unique aspects to its observable signatures and self-consistent internal dynamics. In addition to these intrinsic features, another important implication of Klein–Nishina reconnection is its interaction with its environment. Here, pair production opens up a coupling channel that is unique to QED reconnection:

<table><caption>**Table 2.** Parameters for the six-simulation scan over $\gamma_{\text{KN}}$ discussed in section 5. Values that are the same as those of Table 1 are omitted. For computational expediency, we employ a larger time-step in these runs (cf. Tables 1 and 3). The energy error is thus slightly worse, but is still only above 1 per cent for the $\gamma_{\text{KN}} = 1500, 1900$, and 3400 runs, where it is, respectively, 1.2 per cent, 1.3 per cent, and 1.7 per cent.</caption><thead><tr><th>Symbol</th><th>Value</th><th>(= equivalent)</th></tr></thead><tbody><tr><td>σ<sub>c,0</sub></td><td>10<sup>4</sup></td><td></td></tr><tr><td>σ<sub>h,0</sub></td><td>10<sup>2</sup></td><td></td></tr><tr><td>θ<sub>0</sub></td><td>2.5 × 10<sup>-4</sup>σ<sub>c,0</sub></td><td>=25</td></tr><tr><td>L</td><td>210σ<sub>c,0</sub>ρ<sub>0</sub></td><td></td></tr><tr><td>γ<sub>max</sub></td><td>21σ<sub>c,0</sub></td><td>=2.1 × 10<sup>5</sup></td></tr><tr><td>B<sub>g</sub></td><td>0.1B<sub>0</sub></td><td></td></tr><tr><td>Δx, Δy</td><td>σ<sub>c,0</sub>ρ<sub>0</sub>/24</td><td></td></tr><tr><td>Δt</td><td>0.99Δx/√2c</td><td></td></tr><tr><td>N</td><td>5120</td><td></td></tr><tr><td>γ<sub>cool</sub></td><td>2.3 × 10<sup>-2</sup>σ<sub>c,0</sub></td><td>=230</td></tr><tr><td>γ<sub>rad, T</sub></td><td>0.70σ<sub>c,0</sub></td><td>=7.0 × 10<sup>3</sup></td></tr><tr><td>γ<sub>KN</sub></td><td>(0.091, 0.12, 0.15)σ<sub>c,0</sub></td><td>= ( 910, 1200, 1500)</td></tr><tr><td>...</td><td>( 0.19, 0.27, 0.34)σ<sub>c,0</sub></td><td>= (1900, 2700, 3400)</td></tr><tr><td>τ<sub>γγ</sub></td><td>(2.4, 3.1, 4.0, 5.2, 7.2, 9.1)</td><td></td></tr></tbody></table>

<table><caption>**Table 3.** Parameters for the three-simulation scan over L. Values that are the same as those of Table 1 are omitted. The value of $\gamma_{\text{KN}}$ is adjusted slightly between runs to keep $\gamma_{\text{rad, IC}} = \gamma_{\text{max}}$.</caption><thead><tr><th>Symbol</th><th>Value</th><th>(= equivalent)</th></tr></thead><tbody><tr><td>σ<sub>c,0</sub></td><td>400</td><td></td></tr><tr><td>σ<sub>h,0</sub></td><td>400*</td><td></td></tr><tr><td>θ<sub>0</sub></td><td>0.1</td><td></td></tr><tr><td>L</td><td>(34, 67, 100)σ<sub>c,0</sub>ρ<sub>0</sub></td><td></td></tr><tr><td>γ<sub>max</sub></td><td>(3.4, 6.7, 10.)σ<sub>c,0</sub></td><td>= (1300, 2700, 4000)</td></tr><tr><td>B<sub>g</sub></td><td>0.1B<sub>0</sub></td><td></td></tr><tr><td>Δx, Δy</td><td>σ<sub>c,0</sub>ρ<sub>0</sub>/76</td><td></td></tr><tr><td>N</td><td>(2560, 5120, 7680)</td><td></td></tr><tr><td>γ<sub>cool</sub></td><td>(12, 5.9, 4.0) × 10<sup>-3</sup>σ<sub>c,0</sub></td><td>= (4.7, 2.4, 1.6)</td></tr><tr><td>γ<sub>rad, T</sub></td><td>0.2σ<sub>c,0</sub></td><td>= 80</td></tr><tr><td>γ<sub>KN</sub></td><td>(0.065, 0.055, 0.051)σ<sub>c,0</sub></td><td>= (26, 22, 20.)</td></tr><tr><td>τ<sub>γγ</sub></td><td>(3.3, 5.6, 7.8)</td><td></td></tr></tbody></table>

*Note.* * Because the initial plasma temperature, $\theta_0 = 0.1$, is non-relativistically cold, $\sigma_{h,0} = \sigma_{c,0}$.

the possibility to change the ambient plasma composition (e.g. the positron-to-proton ratio). Thus, in this section, we characterize the total pair yield from a Klein–Nishina reconnection event.

To map the dependence of the pair yield on reconnection parameters, we introduce two auxiliary simulation campaigns. Each one explores the pair yield's dependence on one principal variable. The first campaign explores the main new quantity introduced by Klein–Nishina physics, $\gamma_{\text{KN}}$. The second campaign is a system-size (L) scan. It doubles as an opportunity to diagnose how well our results may generalize to astrophysical Klein–Nishina reconnection, where the layer lengths, L, are expected to be much larger compared to the plasma microscales (e.g. $\sigma_{c,0} \rho_0$) than is possible to simulate. The full sets of parameters used for the respective $\gamma_{\text{KN}}$- and L-scans are summarized in Tables 2 and 3.

Besides examining the impact of their respective variables, these additional sets of simulations are separately calibrated to different fiducial parameters. The $\gamma_{\text{KN}}$-scan has a larger normalized system size, $L = 210\sigma_{c,0}\rho_0$, and lower magnetizations, $\sigma_{c,0} = 10^4$ and $\sigma_{h,0} = 10^2$, than our base run, IC(KN)+PP (for which $L \simeq 91\sigma_{c,0}\rho_0$, $\sigma_{c,0} = 1.2 \times 10^5$, and $\sigma_{h,0} = 1250$; Table 1). The L-scan has still
magnetization $\sigma_h \equiv (B_x^2 + B_y^2)/4\pi w$ [note: the enthalpy density, for our ultrarelativistic particle distributions, is $w = 4(U_{\text{orig.}} + U_{\text{prod.}})/3$; see equation (32) and surrounding discussion]. Because the newborn pairs load the upstream plasma energy density without appreciably changing the number density, they reduce $\sigma_h$ from its native/initial value, $\sigma_{h,0}$, while leaving $\sigma_c \simeq \sigma_{c,0}$ essentially untouched.

Even though the newborn pairs are everywhere less numerous than those originally present, they discernibly contribute to the plasma density [e.g. by changing the hue of the $(n_{\text{orig.}} - n_{\text{prod.}})/(n_{\text{orig.}} + n_{\text{prod.}})$ spatial maps] close to and inside of the reconnection separatrix. Mehlhaff et al. (2021) forecasted that the newborn pairs should begin to contribute non-negligibly to the plasma number density passing across the separatrix when $\sigma_{c,0}$ exceeds $\gamma_{\text{KN}}$ by more than a factor of 10 or so. This is consistent with our findings on the Klein–Nishina reconnection pair yield presented in the next section, which suggest that the IC(KN)+PP run ($\sigma_{c,0} \simeq 20\gamma_{\text{KN}}$) is indeed beginning to border on a regime where the *in situ* produced pairs contribute more significantly to the particle count. For the rest of this section, however, we focus on the present regime where the newborn particles are energetically dense but numerically few.

As discussed by Mehlhaff et al. (2021), the energy density injected into the upstream region via the hot newborn pairs leads to a *pair-loaded magnetization*, $\sigma_h^*$, to which the system tends to regulate itself (provided the initial magnetization, $\sigma_{h,0}$, exceeds $\sigma_h^*$). However, if the pair loading is efficient enough, the system may never actually attain a steady state with $\sigma_h = \sigma_h^*$. It will instead flood the upstream energy density via pair production, overshooting to a lower magnetization, $\sigma_h < \sigma_h^*$, and quenching subsequent NTPA until the upstream plasma is vacated and a high $\sigma_h > \sigma_h^*$ is restored. This would restart the process, leading to a limit cycle: the system would indefinitely ricochet between a high and a low magnetization on either side of $\sigma_h^*$ throughout the duration of reconnection.

Mehlhaff et al. (2021) quantified the dependence of the long-term fate of the system – whether it smoothly regulates to, or violently oscillates about, $\sigma_h^*$ – on the efficiency of upstream pair loading. We measure this efficiency from our simulations in Appendix B, finding that it is too low, based on the analysis of Mehlhaff et al. (2021), to trigger $\sigma_h$-mediated limit cycles. We therefore conclude that *limit cycles mediated by pair-loading of $\sigma_h$ are unlikely in astrophysical Klein–Nishina reconnection*.

In this section, we have shown that our simulations probe a Klein–Nishina reconnection regime where, as previously predicted by Mehlhaff et al. (2021), the newborn upstream pairs are hot but tenuous, loading the upstream plasma pressure, energy density, and hot magnetization, but not the upstream plasma density or cold magnetization. The simulations further provide measurements (Appendix B) that disfavour $\sigma_h$-mediated limit cycles. However, if a regime existed featuring a large multiplicity of newborn pairs, the possibility of $\sigma_c$-mediated limit cycles would remain an open question. It is the overall pair yield of Klein–Nishina reconnection, including the potential existence of such a copious pair-production regime, to which we now turn.

## 5 PAIR YIELD

We have already seen how the pairs produced in Klein–Nishina reconnection contribute unique aspects to its observable signatures and self-consistent internal dynamics. In addition to these intrinsic features, another important implication of Klein–Nishina reconnection is its interaction with its environment. Here, pair production opens up a coupling channel that is unique to QED reconnection:

<table><caption>**Table 2.** Parameters for the six-simulation scan over $\gamma_{\text{KN}}$ discussed in section 5. Values that are the same as those of Table 1 are omitted. For computational expediency, we employ a larger time-step in these runs (cf. Tables 1 and 3). The energy error is thus slightly worse, but is still only above 1 per cent for the $\gamma_{\text{KN}} = 1500, 1900$, and 3400 runs, where it is, respectively, 1.2 per cent, 1.3 per cent, and 1.7 per cent.</caption><thead><tr><th>Symbol</th><th>Value</th><th>(= equivalent)</th></tr></thead><tbody><tr><td>σ<sub>c,0</sub></td><td>10<sup>4</sup></td><td></td></tr><tr><td>σ<sub>h,0</sub></td><td>10<sup>2</sup></td><td></td></tr><tr><td>θ<sub>0</sub></td><td>2.5 × 10<sup>-4</sup>σ<sub>c,0</sub></td><td>=25</td></tr><tr><td>L</td><td>210σ<sub>c,0</sub>ρ<sub>0</sub></td><td></td></tr><tr><td>γ<sub>max</sub></td><td>21σ<sub>c,0</sub></td><td>=2.1 × 10<sup>5</sup></td></tr><tr><td>B<sub>g</sub></td><td>0.1B<sub>0</sub></td><td></td></tr><tr><td>Δx, Δy</td><td>σ<sub>c,0</sub>ρ<sub>0</sub>/24</td><td></td></tr><tr><td>Δt</td><td>0.99Δx/√2c</td><td></td></tr><tr><td>N</td><td>5120</td><td></td></tr><tr><td>γ<sub>cool</sub></td><td>2.3 × 10<sup>-2</sup>σ<sub>c,0</sub></td><td>=230</td></tr><tr><td>γ<sub>rad, T</sub></td><td>0.70σ<sub>c,0</sub></td><td>=7.0 × 10<sup>3</sup></td></tr><tr><td>γ<sub>KN</sub></td><td>(0.091, 0.12, 0.15)σ<sub>c,0</sub></td><td>= ( 910, 1200, 1500)</td></tr><tr><td>...</td><td>( 0.19, 0.27, 0.34)σ<sub>c,0</sub></td><td>= (1900, 2700, 3400)</td></tr><tr><td>τ<sub>γγ</sub></td><td>(2.4, 3.1, 4.0, 5.2, 7.2, 9.1)</td><td></td></tr></tbody></table>

<table><caption>**Table 3.** Parameters for the three-simulation scan over L. Values that are the same as those of Table 1 are omitted. The value of $\gamma_{\text{KN}}$ is adjusted slightly between runs to keep $\gamma_{\text{rad, IC}} = \gamma_{\text{max}}$.</caption><thead><tr><th>Symbol</th><th>Value</th><th>(= equivalent)</th></tr></thead><tbody><tr><td>σ<sub>c,0</sub></td><td>400</td><td></td></tr><tr><td>σ<sub>h,0</sub></td><td>400*</td><td></td></tr><tr><td>θ<sub>0</sub></td><td>0.1</td><td></td></tr><tr><td>L</td><td>(34, 67, 100)σ<sub>c,0</sub>ρ<sub>0</sub></td><td></td></tr><tr><td>γ<sub>max</sub></td><td>(3.4, 6.7, 10.)σ<sub>c,0</sub></td><td>= (1300, 2700, 4000)</td></tr><tr><td>B<sub>g</sub></td><td>0.1B<sub>0</sub></td><td></td></tr><tr><td>Δx, Δy</td><td>σ<sub>c,0</sub>ρ<sub>0</sub>/76</td><td></td></tr><tr><td>N</td><td>(2560, 5120, 7680)</td><td></td></tr><tr><td>γ<sub>cool</sub></td><td>(12, 5.9, 4.0) × 10<sup>-3</sup>σ<sub>c,0</sub></td><td>= (4.7, 2.4, 1.6)</td></tr><tr><td>γ<sub>rad, T</sub></td><td>0.2σ<sub>c,0</sub></td><td>= 80</td></tr><tr><td>γ<sub>KN</sub></td><td>(0.065, 0.055, 0.051)σ<sub>c,0</sub></td><td>= (26, 22, 20.)</td></tr><tr><td>τ<sub>γγ</sub></td><td>(3.3, 5.6, 7.8)</td><td></td></tr></tbody></table>

*Note.* * Because the initial plasma temperature, $\theta_0 = 0.1$, is non-relativistically cold, $\sigma_{h,0} = \sigma_{c,0}$.

the possibility to change the ambient plasma composition (e.g. the positron-to-proton ratio). Thus, in this section, we characterize the total pair yield from a Klein–Nishina reconnection event.

To map the dependence of the pair yield on reconnection parameters, we introduce two auxiliary simulation campaigns. Each one explores the pair yield's dependence on one principal variable. The first campaign explores the main new quantity introduced by Klein–Nishina physics, $\gamma_{\text{KN}}$. The second campaign is a system-size (L) scan. It doubles as an opportunity to diagnose how well our results may generalize to astrophysical Klein–Nishina reconnection, where the layer lengths, L, are expected to be much larger compared to the plasma microscales (e.g. $\sigma_{c,0} \rho_0$) than is possible to simulate. The full sets of parameters used for the respective $\gamma_{\text{KN}}$- and L-scans are summarized in Tables 2 and 3.

Besides examining the impact of their respective variables, these additional sets of simulations are separately calibrated to different fiducial parameters. The $\gamma_{\text{KN}}$-scan has a larger normalized system size, $L = 210\sigma_{c,0}\rho_0$, and lower magnetizations, $\sigma_{c,0} = 10^4$ and $\sigma_{h,0} = 10^2$, than our base run, IC(KN)+PP (for which $L \simeq 91\sigma_{c,0}\rho_0$, $\sigma_{c,0} = 1.2 \times 10^5$, and $\sigma_{h,0} = 1250$; Table 1). The L-scan has still
<figure><img src="image_15.png" /><figcaption>**Figure 15.** Pair yield for the simulations of the γKN parameter scan (blue; Table 2), the L parameter scan (green; Table 3), and the IC(KN)+PP run (red; Table 1). The smallest-L simulation (open circle) is omitted when fitting for the exponential scaling law (solid line with error envelope indicating reported 1σ uncertainty on fit parameters). Despite the large differences in their fiducial and scanned parameters, all simulations' pair yields (with the exception of the smallest one in the L-scan; see the text) fall onto a scaling law with a single control parameter. This parameter is the ratio of two Lorentz factors: that, γ<sub>pp,max</sub> ≃ 30γ<sub>KN</sub>, of a particle that typically scatters photons up to the peak energy of the pair-production cross-section; and the characteristic maximum Lorentz factor, γ<sub>X</sub> ≡ 4σ<sub>c,0</sub>, attained near reconnection X-points. This suggests that the pair yield is controlled mainly by the efficiency of particle acceleration up to and beyond the optimal pair-producing particle energy.</figcaption></figure>

different initial magnetizations, σ<sub>c,0</sub> = σ<sub>h,0</sub> = 400, and features a non-relativistic initial upstream plasma, θ<sub>0</sub> = 0.1. Spreading out, in this way, our auxiliary campaigns around the Klein–Nishina reconnection parameter space helps us identify a reduced set of control parameters (in fact, one single parameter; Fig. 15) that decides the pair yield. These control parameters, in turn, shed light on the main physical mechanisms responsible for the pair-production efficiency while also providing a potential method for estimating this efficiency in astrophysical systems.

We define the pair yield in our simulations as the ratio of the total number, N<sub>prod.</sub>, of leptons (electrons and positrons) produced on the fly to the cumulative count, N<sub>rec.</sub>, of originally present leptons processed by reconnection (i.e. swept across the separatrix). Generally, about 60 per cent of the initial upstream particles cross the separatrix before reconnection saturates, the same as the percentage of the initial magnetic flux that is reconnected (see Section 4.2 and Fig. 5). Though the processing of upstream magnetic field and original particles essentially finishes by the end of our simulations, pair production does not; there remains a prominent population of high-energy particles that have yet to cool down by emitting pair-producing gamma-rays. In Section 4.5, we post-processed the passive cooling associated with these particles to diagnose the decaying phase of the reconnection-powered light curve. Here, we use the same technique to continue evolving the total newborn pair count past the end of each simulation. Once the count saturates (typically by 6L/c or so), we record it as N<sub>prod.</sub>. As long as we begin the post-processing after the energy transfer from fields to particles is mostly complete, the final N<sub>prod.</sub> does not depend much (less than 10 per cent) on the exact moment in time when the post-processing starts.

Let us now comment on the generality of results obtained via the above pair yield measurement recipe. We note that the measurement procedure is easily formulated and performed in the context of our

periodic simulation box. The fact that the processing of magnetic flux, advection of upstream particles across the separatrix, and particle energization all eventually slow down and cease enables a straightforward calculation of the pairs created per reconnection-processed original particle. This does not necessarily limit the applicability of the pair yield of this section to flaring sources with discrete reconnection episodes, however. On the contrary, if the pair yield is truly local, scaling linearly – for a fixed set of parameters (those outlined in Section 3) – with the number of processed upstream particles, then the results would also hold in the context of steady-state reconnection.

As a potential caveat, one should note that, like our periodic simulation boundaries, our pair yield measurement technique ignores the potential for particles and photons to escape the ambient radiation field before pair production is complete. This issue is less important when the pair-production optical depth, τ<sub>γγ</sub>, is large, corresponding to shorter mean-free paths of photons and more rapid particle cooling times (Fig. 2 and Section 3.2.1). When τ<sub>γγ</sub> becomes small, the pair yield calculated in this way still has meaning as long as the extent of the ambient radiation field is much larger than the size of the reconnection system (and provided, once particles exit the reconnection system, they are no longer significantly energized).

Having discussed the general applicability and limitations of our pair yield measurements, we now present the pair yield calculated for our auxiliary γKN and L simulation campaigns, as well as for the single IC(KN)+PP run discussed earlier (Section 4), in Fig. 15. Remarkably, when plotted as a function of just the single control parameter, γ<sub>pp,max</sub>/4σ<sub>c,0</sub>, the pair yields from all simulations – despite the very different fiducial and scanned parameters – collapse onto the same exponential scaling law,

$$N_{\text{prod.}}/N_{\text{rec.}} = 1.5 \exp \left( -1.7 \frac{\gamma_{\text{pp,max}}}{4\sigma_{c,0}} \right). \quad (53)$$

Here, γ<sub>pp,max</sub> ≡ 3.6 × 8 × γKN ≃ 30γKN is the characteristic Lorentz factor of particles that scatter background photons to energies at peak pair-production cross-section, ϵ<sub>scat</sub> ∼ γ<sub>pp,max</sub>/2 ∼ 3.6(m<sub>e</sub>c<sup>2</sup>)<sup>2</sup>/ϵ<sub>ph</sub>. Meanwhile, 4σ<sub>c,0</sub> is the characteristic maximum energy, γ<sub>X</sub> ≡ 4σ<sub>c,0</sub>, that particles acquire near reconnection X-points. Equation (53) suggests that what controls the final pair yield of reconnection is how broad a distribution of high-energy particles can be energized near reconnection X-points to radiate photons close to or above the peak pair-production cross-section: that is, by how much γ<sub>X</sub> exceeds γ<sub>pp,max</sub>.

Mehlhaff et al. (2021) predicted that the ratio of the newborn-to-original upstream pair densities flowing across the reconnection separatrix should be proportional to σ<sub>c,0</sub>/γKN (times a non-trivial function depending on NTPA in the reconnection layer). Although this number density ratio is not the same as the global ratio of newborn-to-reconnection-processed particles measured here – the latter also includes the non-negligible number of pairs born on the exhaust side of the reconnection separatrix – both results share the same main controlling parameter, σ<sub>c,0</sub>/γKN ∝ γX/γ<sub>pp,max</sub>.

Now let us discuss where the scaling (53) may break down. First, all of the simulations in our campaigns have an effective radiative cut-off Lorentz factor, γ<sub>rad, IC</sub> ≥ 4σ<sub>c,0</sub>. This means that we need to be vigilant as γKN increases (coinciding with larger γ<sub>pp,max</sub>/4σ<sub>c,0</sub>), because γ<sub>rad, IC</sub> may then fall below 4σ<sub>c,0</sub>. In that case, our present interpretation of equation (53) leads us to expect another cut-off in the pair yield, for example, of the form, exp(−γ<sub>pp,max</sub>/γ<sub>rad, IC</sub>), to kick in, signalling that X-point particle acceleration has become limited by γ<sub>rad, IC</sub> instead of by σ<sub>c,0</sub>. Second, in the opposite limit, where γ<sub>pp,max</sub>/4σ<sub>c,0</sub> becomes small, we are likely to transition to a regime in which extremely efficient particle acceleration (giving a particle
<figure><img src="image_16.png" /><figcaption>**Figure 15.** Pair yield for the simulations of the γKN parameter scan (blue; Table 2), the L parameter scan (green; Table 3), and the IC(KN)+PP run (red; Table 1). The smallest-L simulation (open circle) is omitted when fitting for the exponential scaling law (solid line with error envelope indicating reported 1σ uncertainty on fit parameters). Despite the large differences in their fiducial and scanned parameters, all simulations' pair yields (with the exception of the smallest one in the L-scan; see the text) fall onto a scaling law with a single control parameter. This parameter is the ratio of two Lorentz factors: that, γpp,max ≃ 30γKN, of a particle that typically scatters photons up to the peak energy of the pair-production cross-section; and the characteristic maximum Lorentz factor, γX ≡ 4σc,0, attained near reconnection X-points. This suggests that the pair yield is controlled mainly by the efficiency of particle acceleration up to and beyond the optimal pair-producing particle energy.</figcaption></figure>

different initial magnetizations, σ<sub>c,0</sub> = σ<sub>h,0</sub> = 400, and features a non-relativistic initial upstream plasma, θ<sub>0</sub> = 0.1. Spreading out, in this way, our auxiliary campaigns around the Klein–Nishina reconnection parameter space helps us identify a reduced set of control parameters (in fact, one single parameter; Fig. 15) that decides the pair yield. These control parameters, in turn, shed light on the main physical mechanisms responsible for the pair-production efficiency while also providing a potential method for estimating this efficiency in astrophysical systems.

We define the pair yield in our simulations as the ratio of the total number, N<sub>prod.</sub>, of leptons (electrons and positrons) produced on the fly to the cumulative count, N<sub>rec.</sub>, of originally present leptons processed by reconnection (i.e. swept across the separatrix). Generally, about 60 per cent of the initial upstream particles cross the separatrix before reconnection saturates, the same as the percentage of the initial magnetic flux that is reconnected (see Section 4.2 and Fig. 5). Though the processing of upstream magnetic field and original particles essentially finishes by the end of our simulations, pair production does not; there remains a prominent population of high-energy particles that have yet to cool down by emitting pair-producing gamma-rays. In Section 4.5, we post-processed the passive cooling associated with these particles to diagnose the decaying phase of the reconnection-powered light curve. Here, we use the same technique to continue evolving the total newborn pair count past the end of each simulation. Once the count saturates (typically by 6L/c or so), we record it as N<sub>prod.</sub>. As long as we begin the post-processing after the energy transfer from fields to particles is mostly complete, the final N<sub>prod.</sub> does not depend much (less than 10 per cent) on the exact moment in time when the post-processing starts.

Let us now comment on the generality of results obtained via the above pair yield measurement recipe. We note that the measurement procedure is easily formulated and performed in the context of our

periodic simulation box. The fact that the processing of magnetic flux, advection of upstream particles across the separatrix, and particle energization all eventually slow down and cease enables a straightforward calculation of the pairs created per reconnection-processed original particle. This does not necessarily limit the applicability of the pair yield of this section to flaring sources with discrete reconnection episodes, however. On the contrary, if the pair yield is truly local, scaling linearly – for a fixed set of parameters (those outlined in Section 3) – with the number of processed upstream particles, then the results would also hold in the context of steady-state reconnection.

As a potential caveat, one should note that, like our periodic simulation boundaries, our pair yield measurement technique ignores the potential for particles and photons to escape the ambient radiation field before pair production is complete. This issue is less important when the pair-production optical depth, τ<sub>γγ</sub>, is large, corresponding to shorter mean-free paths of photons and more rapid particle cooling times (Fig. 2 and Section 3.2.1). When τ<sub>γγ</sub> becomes small, the pair yield calculated in this way still has meaning as long as the extent of the ambient radiation field is much larger than the size of the reconnection system (and provided, once particles exit the reconnection system, they are no longer significantly energized).

Having discussed the general applicability and limitations of our pair yield measurements, we now present the pair yield calculated for our auxiliary γKN and L simulation campaigns, as well as for the single IC(KN)+PP run discussed earlier (Section 4), in Fig. 15. Remarkably, when plotted as a function of just the single control parameter, γ<sub>pp,max</sub>/4σ<sub>c,0</sub>, the pair yields from all simulations – despite the very different fiducial and scanned parameters – collapse onto the same exponential scaling law,

$$N_{\text{prod.}}/N_{\text{rec.}} = 1.5 \exp \left( -1.7 \frac{\gamma_{\text{pp,max}}}{4\sigma_{c,0}} \right). \quad (53)$$

Here, γ<sub>pp,max</sub> ≡ 3.6 × 8 × γKN ≃ 30γKN is the characteristic Lorentz factor of particles that scatter background photons to energies at peak pair-production cross-section, ϵ<sub>scat</sub> ∼ γpp,max/2 ∼ 3.6(m<sub>e</sub>c<sup>2</sup>)<sup>2</sup>/ϵ<sub>ph</sub>. Meanwhile, 4σ<sub>c,0</sub> is the characteristic maximum energy, γX ≡ 4σ<sub>c,0</sub>, that particles acquire near reconnection X-points. Equation (53) suggests that what controls the final pair yield of reconnection is how broad a distribution of high-energy particles can be energized near reconnection X-points to radiate photons close to or above the peak pair-production cross-section: that is, by how much γX exceeds γpp,max.

Mehlhaff et al. (2021) predicted that the ratio of the newborn-to-original upstream pair densities flowing across the reconnection separatrix should be proportional to σ<sub>c,0</sub>/γKN (times a non-trivial function depending on NTPA in the reconnection layer). Although this number density ratio is not the same as the global ratio of newborn-to-reconnection-processed particles measured here – the latter also includes the non-negligible number of pairs born on the exhaust side of the reconnection separatrix – both results share the same main controlling parameter, σ<sub>c,0</sub>/γKN ∝ γX/γpp,max.

Now let us discuss where the scaling (53) may break down. First, all of the simulations in our campaigns have an effective radiative cut-off Lorentz factor, γ<sub>rad, IC</sub> ≥ 4σ<sub>c,0</sub>. This means that we need to be vigilant as γKN increases (coinciding with larger γpp,max/4σ<sub>c,0</sub>), because γ<sub>rad, IC</sub> may then fall below 4σ<sub>c,0</sub>. In that case, our present interpretation of equation (53) leads us to expect another cut-off in the pair yield, for example, of the form, exp(−γpp,max/γ<sub>rad, IC</sub>), to kick in, signalling that X-point particle acceleration has become limited by γ<sub>rad, IC</sub> instead of by σ<sub>c,0</sub>. Second, in the opposite limit, where γpp,max/4σ<sub>c,0</sub> becomes small, we are likely to transition to a regime in which extremely efficient particle acceleration (giving a particle
distribution power-law index approaching the non-radiative limit, $p \simeq 1$) causes the pair yield to continue to grow with $4\sigma_{c,0}/\gamma_{pp,\text{max}}$. Though this is not what a naive extrapolation of (53) produces, it is more coherent with our present interpretation of that formula, in which more particles being accelerated into the optimal pair-producing energy range enhances the pair yield.

Let us examine, for the sake of argument, what might happen if pushing $\gamma_{\text{pp,max}}/4\sigma_{c,0}$ to smaller values than studied here indeed led to a copious pair yield. Then Klein–Nishina reconnection would begin to move out of the regime of our simulations (discussed in Section 4.6 and Appendix B) where the produced particles are hot but tenuous. Instead, the newborn particles would become hot and abundant, potentially even loading the upstream cold magnetization to a reduced value, $\sigma_c^* < \sigma_{c,0}$, whereas in our simulations they modify only the hot magnetization (Section 4.6). If a steady state were to exist in such a regime, the cold magnetization would need to self-adjust until the pair yield as defined in equation (53) attains unity (times some efficiency factor to account for the fact that not all of the produced particles would be born into the upstream region; cf. Appendix B). According to Fig. 15, this occurs when $\gamma_{\text{pp,max}}/4\sigma_c^* \simeq 0.25$, or, equivalently, $\sigma_c^* \simeq \gamma_{\text{pp,max}} \simeq 30\gamma_{\text{KN}}$. Such a result would be highly attractive, for it would open up the possibility of measuring the lepton material composition (via $\sigma_c^*$) in terms of the seed photon energy, a much easier quantity to infer observationally.

Lastly, we note that the trend identified in Fig. 15 is broken at smaller system sizes. Indeed, we excluded the pair yield measurement from our $L = 34\sigma_{c,0}\rho_0$ simulation when fitting for equation (53). The break from the formula in smaller systems reflects findings of earlier simulation studies that collisionless reconnection transitions to an asymptotically large-system limit, corresponding to the multiple X-point, plasmoid-mediated regime, only once $L \gtrsim 40\sigma_{c,0}\rho_0$ (Werner et al. 2016). When respecting this limit, our simulations overlap the identified pair yield trend. However, we cannot rule out larger systems yielding even more efficient pair production than (53).

In this section, we have characterized the pair yield of Klein–Nishina reconnection in terms of a 1D exponential scaling law, equation (53) – despite the high-dimensional parameter space of this problem. This scaling law appears to be robust across an order of magnitude or more in $\gamma_{\text{KN}}$, $\sigma_{\text{h},0}$, $\sigma_{c,0}$, and $\theta_0$, while being respected across a factor of 3 or 4 in system size (the most that we can afford to probe in the large-system, $L > 40\sigma_{c,0}\rho_0$, regime given the stringent parameter constraints of the problem; cf. Section 3). Equation (53) may need to be modified in the transition region between Klein–Nishina ($\gamma_{\text{KN}} < \gamma_{\text{rad,T}}$) and Thomson ($\gamma_{\text{KN}} > \gamma_{\text{rad,T}}$) radiative reconnection, and it may also give way to a new regime of efficient (much greater than order-unity) pair yield when $\gamma_{\text{pp,max}}/4\sigma_{c,0}$ becomes much smaller than the values we test. Nevertheless, we are able to capture an order of magnitude in the control parameter, $\gamma_{\text{pp,max}}/4\sigma_{c,0}$, including near the point, $\gamma_{\text{pp,max}}/4\sigma_{c,0} \simeq 0.25$, where the pair yield reaches 1.

# 6 DISCUSSION

In this section, we discuss the relevance of our findings to gamma-ray observations of selected astrophysical systems: FSRQs, black hole accretion disc coronae, the M87* magnetosphere, and gamma-ray binaries. In each case, we argue why the operation of Klein–Nishina reconnection in these systems is expected on theoretical grounds. We further discuss consequences, for each system, of the results of Sections 4.5 and 5. For reference, we briefly recapitulate those results here in the context of potential links to observations.

The main finding of Section 4.5 is the marked departure of the observable signatures of radiative reconnection in the Klein–Nishina regime (with pair production) from the Thomson-cooled regime (without pair production). Namely, while Thomson-cooled reconnection features a tight correlation between the mean observed photon energy and the system’s total luminosity (i.e. ‘harder-when-brighter’), Klein–Nishina reconnection breaks this correlation, featuring a steady mean photon energy irrespective of the luminosity. This includes the decaying phase of a flaring event, wherein the average photon energy is preserved even as the luminosity drops. These findings serve therefore as an observational diagnostic. They can be directly compared to gamma-ray observations to build a case (or not) for the operation of Klein–Nishina reconnection in a given object. In contrast, the results of Section 5 function as an inference criterion, providing a method to estimate a quantity that is difficult to constrain from observations – the emitting region’s material composition – using quantities that may be more readily measured or estimated (specifically, the magnetization, $\sigma_{c,0}$, and the Lorentz factor, $\gamma_{\text{pp,max}}$, of particles whose photons are at maximum pair-production cross-section with the radiation bath).

As an aid to navigating the astrophysical diversity treated by the following discussion, each system-specific section (6.1–6.4) is further subdivided into three parts: an overview of pertinent background information, the principal discussion of connections to this work, and a brief object-specific synopsis. Following our system-by-system presentation, we summarize our broad conclusions across all systems in Section 6.5 and Table 4. Readers interested in a particular object may skip directly to the corresponding subsection; those seeking a more general overview may wish to skip first to the end summary material (Section 6.5).

## 6.1 Flat-spectrum radio quasars and other blazars

### 6.1.1 Background information

Blazars are AGNs that launch bipolar relativistic jets, one of which (hereafter, the singular ‘jet’) travels toward the Earth. The jet’s relativistic motion Doppler boosts its emission, leading to dramatic observable consequences. For example, blazars dominate the discrete sources on the extragalactic gamma-ray sky (e.g. Wakely & Horan 2008; Abdollahi et al. 2020) and, in the optical band, they routinely outshine the cumulative starlight of their host galaxies (Olguín-Iglesias et al. 2016). Blazar jet emission is also exceptionally broad, extending from radio frequencies up to gamma-rays in a characteristic non-thermal double-humped spectrum (Fossati et al. 1998; Ghisellini 2011; Ghisellini et al. 2017). The lower energy spectral hump originates from synchrotron radiation by relativistic electrons and positrons (henceforth ‘leptons’) spiraling around magnetic field lines in the jet. The higher-energy peak is frequently attributed to IC radiation also by relativistic jet leptons (e.g. Prandini & Ghisellini 2022).

Blazars are phenomenologically decomposed into two main subdivisions: FSRQs and BL Lacs.<fnref n="7" /> FSRQs exhibit lower energy spectra, with the synchrotron component peaking at infrared energies and the higher-energy IC component peaking in the MeV-to-GeV gamma-rays. In contrast, the maxima of the synchrotron and IC spectra

<footnote n="7">⁷For our purposes, we fold the extreme/ultrahigh-frequency-peaked BL Lacs (EHBLs/UHBLs) high-frequency-peaked BL Lacs (HBLs), intermediate-frequency-peaked BL Lacs (IBLs), and low-frequency-peaked BL Lacs (LBLs) into the BL Lac class.</footnote>
distribution power-law index approaching the non-radiative limit, $p \simeq 1$) causes the pair yield to continue to grow with $4\sigma_{c,0}/\gamma_{pp,\text{max}}$. Though this is not what a naive extrapolation of (53) produces, it is more coherent with our present interpretation of that formula, in which more particles being accelerated into the optimal pair-producing energy range enhances the pair yield.

Let us examine, for the sake of argument, what might happen if pushing $\gamma_{pp,\text{max}}/4\sigma_{c,0}$ to smaller values than studied here indeed led to a copious pair yield. Then Klein–Nishina reconnection would begin to move out of the regime of our simulations (discussed in Section 4.6 and Appendix B) where the produced particles are hot but tenuous. Instead, the newborn particles would become hot and abundant, potentially even loading the upstream cold magnetization to a reduced value, $\sigma_c^* < \sigma_{c,0}$, whereas in our simulations they modify only the hot magnetization (Section 4.6). If a steady state were to exist in such a regime, the cold magnetization would need to self-adjust until the pair yield as defined in equation (53) attains unity (times some efficiency factor to account for the fact that not all of the produced particles would be born into the upstream region; cf. Appendix B). According to Fig. 15, this occurs when $\gamma_{pp,\text{max}}/4\sigma_c^* \simeq 0.25$, or, equivalently, $\sigma_c^* \simeq \gamma_{pp,\text{max}} \simeq 30\gamma_{\text{KN}}$. Such a result would be highly attractive, for it would open up the possibility of measuring the lepton material composition (via $\sigma_c^*$) in terms of the seed photon energy, a much easier quantity to infer observationally.

Lastly, we note that the trend identified in Fig. 15 is broken at smaller system sizes. Indeed, we excluded the pair yield measurement from our $L = 34\sigma_{c,0}\rho_0$ simulation when fitting for equation (53). The break from the formula in smaller systems reflects findings of earlier simulation studies that collisionless reconnection transitions to an asymptotically large-system limit, corresponding to the multiple X-point, plasmoid-mediated regime, only once $L \gtrsim 40\sigma_{c,0}\rho_0$ (Werner et al. 2016). When respecting this limit, our simulations overlap the identified pair yield trend. However, we cannot rule out larger systems yielding even more efficient pair production than (53).

In this section, we have characterized the pair yield of Klein–Nishina reconnection in terms of a 1D exponential scaling law, equation (53) – despite the high-dimensional parameter space of this problem. This scaling law appears to be robust across an order of magnitude or more in $\gamma_{\text{KN}}$, $\sigma_{h,0}$, $\sigma_{c,0}$, and $\theta_0$, while being respected across a factor of 3 or 4 in system size (the most that we can afford to probe in the large-system, $L > 40\sigma_{c,0}\rho_0$, regime given the stringent parameter constraints of the problem; cf. Section 3). Equation (53) may need to be modified in the transition region between Klein–Nishina ($\gamma_{\text{KN}} < \gamma_{\text{rad, T}}$) and Thomson ($\gamma_{\text{KN}} > \gamma_{\text{rad, T}}$) radiative reconnection, and it may also give way to a new regime of efficient (much greater than order-unity) pair yield when $\gamma_{\text{pp, max}}/4\sigma_{c,0}$ becomes much smaller than the values we test. Nevertheless, we are able to capture an order of magnitude in the control parameter, $\gamma_{\text{pp, max}}/4\sigma_{c,0}$, including near the point, $\gamma_{\text{pp, max}}/4\sigma_{c,0} \simeq 0.25$, where the pair yield reaches 1.

# 6 DISCUSSION

In this section, we discuss the relevance of our findings to gamma-ray observations of selected astrophysical systems: FSRQs, black hole accretion disc coronae, the M87* magnetosphere, and gamma-ray binaries. In each case, we argue why the operation of Klein–Nishina reconnection in these systems is expected on theoretical grounds. We further discuss consequences, for each system, of the results of Sections 4.5 and 5. For reference, we briefly recapitulate those results here in the context of potential links to observations.

The main finding of Section 4.5 is the marked departure of the observable signatures of radiative reconnection in the Klein–Nishina regime (with pair production) from the Thomson-cooled regime (without pair production). Namely, while Thomson-cooled reconnection features a tight correlation between the mean observed photon energy and the system’s total luminosity (i.e. ‘harder-when-brighter’), Klein–Nishina reconnection breaks this correlation, featuring a steady mean photon energy irrespective of the luminosity. This includes the decaying phase of a flaring event, wherein the average photon energy is preserved even as the luminosity drops. These findings serve therefore as an observational diagnostic. They can be directly compared to gamma-ray observations to build a case (or not) for the operation of Klein–Nishina reconnection in a given object. In contrast, the results of Section 5 function as an inference criterion, providing a method to estimate a quantity that is difficult to constrain from observations – the emitting region’s material composition – using quantities that may be more readily measured or estimated (specifically, the magnetization, $\sigma_{c,0}$, and the Lorentz factor, $\gamma_{\text{pp, max}}$, of particles whose photons are at maximum pair-production cross-section with the radiation bath).

As an aid to navigating the astrophysical diversity treated by the following discussion, each system-specific section (6.1–6.4) is further subdivided into three parts: an overview of pertinent background information, the principal discussion of connections to this work, and a brief object-specific synopsis. Following our system-by-system presentation, we summarize our broad conclusions across all systems in Section 6.5 and Table 4. Readers interested in a particular object may skip directly to the corresponding subsection; those seeking a more general overview may wish to skip first to the end summary material (Section 6.5).

## 6.1 Flat-spectrum radio quasars and other blazars

### 6.1.1 Background information

Blazars are AGNs that launch bipolar relativistic jets, one of which (hereafter, the singular ‘jet’) travels toward the Earth. The jet’s relativistic motion Doppler boosts its emission, leading to dramatic observable consequences. For example, blazars dominate the discrete sources on the extragalactic gamma-ray sky (e.g. Wakely & Horan 2008; Abdollahi et al. 2020) and, in the optical band, they routinely outshine the cumulative starlight of their host galaxies (Olguín-Iglesias et al. 2016). Blazar jet emission is also exceptionally broad, extending from radio frequencies up to gamma-rays in a characteristic non-thermal double-humped spectrum (Fossati et al. 1998; Ghisellini 2011; Ghisellini et al. 2017). The lower energy spectral hump originates from synchrotron radiation by relativistic electrons and positrons (henceforth ‘leptons’) spiraling around magnetic field lines in the jet. The higher-energy peak is frequently attributed to IC radiation also by relativistic jet leptons (e.g. Prandini & Ghisellini 2022).

Blazars are phenomenologically decomposed into two main subdivisions: FSRQs and BL Lacs.<fnref n="7" /> FSRQs exhibit lower energy spectra, with the synchrotron component peaking at infrared energies and the higher-energy IC component peaking in the MeV-to-GeV gamma-rays. In contrast, the maxima of the synchrotron and IC spectra

<footnote n="7">⁷For our purposes, we fold the extreme/ultrahigh-frequency-peaked BL Lacs (EHBLs/UHBLs) high-frequency-peaked BL Lacs (HBLs), intermediate-frequency-peaked BL Lacs (IBLs), and low-frequency-peaked BL Lacs (LBLs) into the BL Lac class.</footnote>
<table><caption>**Table 4.** A graphic summary of Section 6.</caption><thead><tr><th>(1) Object class</th><th>(2) KN hierarchy realized?</th><th>(3) Observational connection</th><th>(4) Limiting factor</th><th>(5) Explanation</th></tr></thead><tbody><tr><td>FSRQs</td><td></td><td></td><td></td><td></td></tr><tr><td>IC(BLR) scenario</td><td>Yes</td><td>Strong</td><td>–</td><td>Fermi-LAT observations are broadly consistent with anticipated temporal-spectral signatures of KN reconnection (Section 4.5; Figs 16–18)</td></tr><tr><td>IC(HDR) scenario</td><td>Yes</td><td>Limited</td><td>Instrumental</td><td>The upcoming CTA will provide enhanced temporal resolution at the relevant TeV energies, making possible more explicit comparisons with this work</td></tr><tr><td>Black hole accretion disc coronae</td><td>Yes</td><td>Limited</td><td>Sources</td><td>Sources are not bright enough in the relevant 10 + MeV range to probe temporal variability. However, the spectral cut-off in Cyg X-1 is potentially consistent with a KN reconnection model</td></tr><tr><td>M87* magnetosphere</td><td>Yes</td><td>–</td><td>Modelling</td><td>Need to account for synchrotron cooling in order to make relevant predictions</td></tr><tr><td>Gamma-ray binaries</td><td>Yes</td><td>Limited</td><td>Modelling</td><td>Orbital modulations of light curves at the relevant GeV energies necessitate global modelling</td></tr></tbody></table>

Notes. Column (1) indicates the object class (roughly one per Sections 6.1–6.4). Column (2) indicates whether basic theoretical estimates suggest that the Klein–Nishina reconnection scale hierarchy (50), $\gamma_{\text{cool}} < \gamma_{\text{KN}} < \gamma_{\text{rad, T}} < \gamma_{\text{max}}$, is achieved in the given system. Column (3) indicates our judgment of how firmly the results of this study – in particular the observable signatures of Klein–Nishina reconnection discussed in Section 4.5 – can be connected to presently available observational data. If strong observational connections cannot be made, column (4) indicates what, in our view, is the primary reason for this. Finally, column (5) gives a short explanation of the judgments in columns (3) and (4).

in BL Lacs typically lie in the UV/X-ray and GeV-to-TeV bands, respectively. Despite their lower photon energies, FSRQs are more luminous and exhibit much larger ratios of IC-to-synchrotron power. Finally, where the norm for BL Lacs is featureless non-thermal spectra, FSRQs usually exhibit prominent broad emission lines or quasi-thermal radiation at lower energies. These are usually attributed to emission by the underlying AGN accretion disc and to reprocessing of the accretion disc light by circumnuclear material. (Illustrative references pertaining to this entire paragraph include: Fossati et al. 1998; Ghisellini 2011; Madejski & Sikora 2016; Ghisellini et al. 2017; Blandford, Meier & Readhead 2019; Prandini & Ghisellini 2022.)

The circumnuclear regions that are observed at lower energies in FSRQs can provide intense sources of seed photons for IC emission in the jet (Begelman & Sikora 1987; Melia & Königl 1989; Sikora, Begelman & Rees 1994). Particularly bright are the broad emission line region (BLR) and the hot dust region (HDR). Of these two, the BLR is smaller, occupying an inner zone (up to roughly 0.1 pc from the nucleus) where irradiation from the accretion disc ionizes the ambient gas, and subsequent recombination emits line emission, broadened by rapid orbital motion, of characteristic UV energy

$$ \epsilon_{\text{BLR}} = 10 \text{ eV}, \qquad (54) $$

onto the jet (Tavecchio & Ghisellini 2008; Sikora et al. 2009; Nalewajko et al. 2012). The outer circumnuclear reprocessing region (up to roughly 4 pc from the nucleus) is the HDR, which comprises dust radiatively heated by the accretion disc light up to a temperature of about $T_{\text{HDR}} \sim 1200$ K. The hot dust shines a quasi-thermal spectrum onto the jet of characteristic energy (Nenkova et al. 2008a, b; Sikora et al. 2009; Nalewajko et al. 2012)

$$ \epsilon_{\text{HDR}} \sim 3k_B T_{\text{HDR}} = 0.3 \text{ eV}. \qquad (55) $$

The radiation fields from the BLR and the HDR furnish excellent conditions for comparing with our simulations. First, they are energetically dense at distances far from the central engine such that the magnetic field energy density is small compared to that of the seed photons, $B_0^2/8\pi \ll U_{\text{ph}}$, a necessary condition for neglecting synchrotron losses, as we do in our simulations. Second, the resulting

radiation field is expected to be homogeneous not just across the reconnection region, but also across the whole jet width. This creates a direct opportunity for applying our pair yield law found in Section 5, which ignores the possibility of above-threshold photons escaping the ambient radiation field before being absorbed to produce electron–positron pairs.

We conduct a detailed analysis of scenarios where reconnection powers high-energy IC emission in FSRQ jets in our previous analytic work (Mehlhaff et al. 2021). There, we estimate the Lorentz factor energy scales $\gamma_{\text{cool}}$, $\gamma_{\text{KN}}$, $\gamma_{\text{rad, T}}$, and $\gamma_{\text{max}}$ either for the case where the reconnection region lies within the more energetically dense BLR, and thus leptons scatter primarily BLR photons, or for the case where the reconnection region is outside the BLR but inside the HDR such that the BLR radiation field is diluted and the HDR supplies the dominant seed photons. Our estimates in both scenarios yield fiducial energy scales that are in the required order, $\gamma_{\text{cool}} < \gamma_{\text{KN}} < \gamma_{\text{rad, T}} < \gamma_{\text{max}}$ [equation (50)], to realize Klein–Nishina reconnection.

In Mehlhaff et al. (2021), we also pointed out that the BLR and HDR are optically thick to gamma-rays above the pair-production threshold energies

$$ \epsilon_{\text{th,BLR}} = m_e^2 c^4 / \epsilon_{\text{BLR}} = 30 \text{ GeV} \qquad (56) $$

and

$$ \epsilon_{\text{th,HDR}} = m_e^2 c^4 / \epsilon_{\text{HDR}} = 0.9 \text{ TeV}, \qquad (57) $$

respectively. The corresponding characteristic energies radiated by a particle of energy $\gamma_{\text{KN}}$ in each case are $\gamma_{\text{KN,BLR}}m_e c^2/2 = \epsilon_{\text{th,BLR}}/8 = 3$ GeV and $\gamma_{\text{KN,HDR}}m_e c^2/2 = \epsilon_{\text{th,HDR}}/8 = 0.1$ TeV. This means that observations by the Fermi Large Area Telescope (LAT), which is sensitive roughly to energies in the 0.1 – 100 GeV range (Atwood et al. 2009), are able to probe emission by particles at $\gamma_{\text{KN,BLR}}$ up through the BLR gamma-ray absorption cut-off. At the same time, Imaging Atmospheric Cherenkov Telescopes (IACTs), typically sensitive in the 0.1 – 10 TeV band (Cherenkov Telescope Array Consortium et al. 2019), stand best to capture the analogous physics for the IC(HDR) scenario. This is fortuitous because particles with energies near $\gamma_{\text{KN}}$ are precisely those responsible for the characteristic spectral and timing signatures of Klein–Nishina reconnection uncovered in
<table><caption>**Table 4.** A graphic summary of Section 6.</caption><thead><tr><th>(1) Object class</th><th>(2) KN hierarchy realized?</th><th>(3) Observational connection</th><th>(4) Limiting factor</th><th>(5) Explanation</th></tr></thead><tbody><tr><td>FSRQs</td><td></td><td></td><td></td><td></td></tr><tr><td>IC(BLR) scenario</td><td>Yes</td><td>Strong</td><td>–</td><td>Fermi-LAT observations are broadly consistent with anticipated temporal-spectral signatures of KN reconnection (Section 4.5; Figs 16–18)</td></tr><tr><td>IC(HDR) scenario</td><td>Yes</td><td>Limited</td><td>Instrumental</td><td>The upcoming CTA will provide enhanced temporal resolution at the relevant TeV energies, making possible more explicit comparisons with this work</td></tr><tr><td>Black hole accretion disc coronae</td><td>Yes</td><td>Limited</td><td>Sources</td><td>Sources are not bright enough in the relevant 10 + MeV range to probe temporal variability. However, the spectral cut-off in Cyg X-1 is potentially consistent with a KN reconnection model</td></tr><tr><td>M87* magnetosphere</td><td>Yes</td><td>–</td><td>Modelling</td><td>Need to account for synchrotron cooling in order to make relevant predictions</td></tr><tr><td>Gamma-ray binaries</td><td>Yes</td><td>Limited</td><td>Modelling</td><td>Orbital modulations of light curves at the relevant GeV energies necessitate global modelling</td></tr></tbody></table>

Notes. Column (1) indicates the object class (roughly one per Sections 6.1–6.4). Column (2) indicates whether basic theoretical estimates suggest that the Klein–Nishina reconnection scale hierarchy (50), $\gamma_{\text{cool}} < \gamma_{\text{KN}} < \gamma_{\text{rad, T}} < \gamma_{\text{max}}$, is achieved in the given system. Column (3) indicates our judgment of how firmly the results of this study – in particular the observable signatures of Klein–Nishina reconnection discussed in Section 4.5 – can be connected to presently available observational data. If strong observational connections cannot be made, column (4) indicates what, in our view, is the primary reason for this. Finally, column (5) gives a short explanation of the judgments in columns (3) and (4).

in BL Lacs typically lie in the UV/X-ray and GeV-to-TeV bands, respectively. Despite their lower photon energies, FSRQs are more luminous and exhibit much larger ratios of IC-to-synchrotron power. Finally, where the norm for BL Lacs is featureless non-thermal spectra, FSRQs usually exhibit prominent broad emission lines or quasi-thermal radiation at lower energies. These are usually attributed to emission by the underlying AGN accretion disc and to reprocessing of the accretion disc light by circumnuclear material. (Illustrative references pertaining to this entire paragraph include: Fossati et al. 1998; Ghisellini 2011; Madejski & Sikora 2016; Ghisellini et al. 2017; Blandford, Meier & Readhead 2019; Prandini & Ghisellini 2022.)

The circumnuclear regions that are observed at lower energies in FSRQs can provide intense sources of seed photons for IC emission in the jet (Begelman & Sikora 1987; Melia & Königl 1989; Sikora, Begelman & Rees 1994). Particularly bright are the broad emission line region (BLR) and the hot dust region (HDR). Of these two, the BLR is smaller, occupying an inner zone (up to roughly 0.1 pc from the nucleus) where irradiation from the accretion disc ionizes the ambient gas, and subsequent recombination emits line emission, broadened by rapid orbital motion, of characteristic UV energy

$$ \epsilon_{\text{BLR}} = 10 \, \text{eV}, \qquad (54) $$

onto the jet (Tavecchio & Ghisellini 2008; Sikora et al. 2009; Nalewajko et al. 2012). The outer circumnuclear reprocessing region (up to roughly 4 pc from the nucleus) is the HDR, which comprises dust radiatively heated by the accretion disc light up to a temperature of about $T_{\text{HDR}} \sim 1200$ K. The hot dust shines a quasi-thermal spectrum onto the jet of characteristic energy (Nenkova et al. 2008a, b; Sikora et al. 2009; Nalewajko et al. 2012)

$$ \epsilon_{\text{HDR}} \sim 3k_B T_{\text{HDR}} = 0.3 \, \text{eV}. \qquad (55) $$

The radiation fields from the BLR and the HDR furnish excellent conditions for comparing with our simulations. First, they are energetically dense at distances far from the central engine such that the magnetic field energy density is small compared to that of the seed photons, $B_0^2/8\pi \ll U_{\text{ph}}$, a necessary condition for neglecting synchrotron losses, as we do in our simulations. Second, the resulting

radiation field is expected to be homogeneous not just across the reconnection region, but also across the whole jet width. This creates a direct opportunity for applying our pair yield law found in Section 5, which ignores the possibility of above-threshold photons escaping the ambient radiation field before being absorbed to produce electron–positron pairs.

We conduct a detailed analysis of scenarios where reconnection powers high-energy IC emission in FSRQ jets in our previous analytic work (Mehlhaff et al. 2021). There, we estimate the Lorentz factor energy scales $\gamma_{\text{cool}}$, $\gamma_{\text{KN}}$, $\gamma_{\text{rad, T}}$, and $\gamma_{\text{max}}$ either for the case where the reconnection region lies within the more energetically dense BLR, and thus leptons scatter primarily BLR photons, or for the case where the reconnection region is outside the BLR but inside the HDR such that the BLR radiation field is diluted and the HDR supplies the dominant seed photons. Our estimates in both scenarios yield fiducial energy scales that are in the required order, $\gamma_{\text{cool}} < \gamma_{\text{KN}} < \gamma_{\text{rad, T}} < \gamma_{\text{max}}$ [equation (50)], to realize Klein–Nishina reconnection.

In Mehlhaff et al. (2021), we also pointed out that the BLR and HDR are optically thick to gamma-rays above the pair-production threshold energies

$$ \epsilon_{\text{th,BLR}} = m_e^2 c^4 / \epsilon_{\text{BLR}} = 30 \, \text{GeV} \qquad (56) $$

and

$$ \epsilon_{\text{th,HDR}} = m_e^2 c^4 / \epsilon_{\text{HDR}} = 0.9 \, \text{TeV}, \qquad (57) $$

respectively. The corresponding characteristic energies radiated by a particle of energy $\gamma_{\text{KN}}$ in each case are $\gamma_{\text{KN,BLR}}m_e c^2/2 = \epsilon_{\text{th,BLR}}/8 = 3$ GeV and $\gamma_{\text{KN,HDR}}m_e c^2/2 = \epsilon_{\text{th,HDR}}/8 = 0.1$ TeV. This means that observations by the Fermi Large Area Telescope (LAT), which is sensitive roughly to energies in the 0.1–100 GeV range (Atwood et al. 2009), are able to probe emission by particles at $\gamma_{\text{KN,BLR}}$ up through the BLR gamma-ray absorption cut-off. At the same time, Imaging Atmospheric Cherenkov Telescopes (IACTs), typically sensitive in the 0.1–10 TeV band (Cherenkov Telescope Array Consortium et al. 2019), stand best to capture the analogous physics for the IC(HDR) scenario. This is fortuitous because particles with energies near $\gamma_{\text{KN}}$ are precisely those responsible for the characteristic spectral and timing signatures of Klein–Nishina reconnection uncovered in
<figure><img src="image_16.png" /><figcaption>**Figure 16.** Fermi-LAT photon flux (top) and spectral index (bottom), both binned into one-week time intervals, detected from 3C 279 during the period presented by Hayashida et al. (2012). Error bars indicate symmetric Gaussian equivalent $1\sigma$ error. Top and bottom panels, correspond, respectively, to panels (c) and (f) of those authors’ fig. 1. The data were retrieved from the Fermi-LAT Light Curve Repository (Abdollahi et al. 2023). The photon index, $\Gamma$, is defined such that the flux of photons between energies $\epsilon$ and $\epsilon + d\epsilon$ is proportional to $\epsilon^{-\Gamma}$. This is connected to the $\epsilon F(\epsilon)$ representation, for example, of Fig. 9, in that $\epsilon F(\epsilon) \propto \epsilon^{-\Gamma + 2}$.</figcaption></figure>

Section 4.5. These particles radiate just below pair-production threshold, producing the highest energy observable photons, and they are actively replenished by radiative reprocessing from higher (above-threshold and, hence, invisible) photon energies, which stabilizes the observed spectral energy density. Thus, FSRQ flares in the GeV and TeV bands are ideally suited to probe the expected observational signatures of reconnection in the regime studied in this work.

### 6.1.2 Connection to this work

Let us examine a few instructive flares observed by the Fermi-LAT [which, again, probes the IC(BLR) scenario] from the FSRQ, 3C 279. Hayashida et al. (2012) report 3C 279 flaring periods over the first two years of Fermi operations. For reference, we reproduce the light curve and photon index time-series presented by those authors in their fig. 1, as well as the correlation between the total gamma-ray flux and spectral shape shown in their fig. 3, in our respective Figs 16 and 17 using data retrieved from the Fermi LAT Light Curve Repository (Abdollahi et al. 2023). Hayashida et al. (2012) note mild ‘harder-than-brighter’ behaviour over the entire observation period, but this is somewhat quenched during the brightest periods (Fig. 17), during which the Fermi-measured spectral index becomes remarkably flux-independent (e.g. between MJD 54 800 and 54 900 as well as between MJD 55 000 and 55 100 in Fig. 16). Roughly similar behaviour is often, but not uniformly, seen in later observations of the same object. In two even brighter outbursts from 3C 279 reported by Hayashida et al. (2015) and Ackermann et al. (2016), the gamma-ray flux reached high-enough levels to reconstruct spectra for individual orbits of the Fermi satellite. In the latter event, the flaring individual-orbit spectra reveal a photon energy index that remains between about 1.9 and 2.1 while the flux varies across about a factor of 3 (table 1 of Ackermann et al. 2016). The pre- and post-outburst phases of the event also feature rather stable photon indices (though this appears more statistically significant for the pre-outburst phase; see fig. 1 of Ackermann et al. 2016). As an example of when such spectral stability is not seen, the first flare analysed by Hayashida

<figure><img src="image_17.png" /><figcaption>**Figure 17.** The data from Fig. 16 re-arranged into a scatter plot of gamma-ray spectral index versus photon flux. This figure roughly corresponds to fig. 3, panel (A) of Hayashida et al. (2012). The stabilization of the spectrum during high-flux periods, notably from MJD 54 800 to 54 900 and from MJD 55 000 to 55 100 in Fig. 16, is quite reminiscent of that found in Section 4.5 (cf. Fig. 12).</figcaption></figure>

et al. (2015) exhibits an extreme hardening of the photon index that then softens on the trailing edge of the flare.

To provide a complete view of the correlation between spectral hardness and luminosity of 3C 279 in the GeV band, we also supply Fig. 18. This is similar to Fig. 17 except that the Fermi-LAT spectral index is plotted against the gamma-ray flux level for all 14+ yr of archived data. The full set of Fermi observations demonstrates broad consistency with, for example, the individual flaring period presented by Hayashida et al. (2012) and displayed in Figs 16 and 17: as the object brightens, the variation in photon index appears to decrease, suggesting a stabler spectrum during flaring periods.

At an even more general level, Meyer, Scargle & Blandford (2019) conducted a statistical analysis of the brightest flares from 6 of the most luminous FSRQs detected by Fermi. They find hints of ‘harder-
<figure><img src="image_17.png" /><figcaption>**Figure 16.** Fermi-LAT photon flux (top) and spectral index (bottom), both binned into one-week time intervals, detected from 3C 279 during the period presented by Hayashida et al. (2012). Error bars indicate symmetric Gaussian equivalent $1\sigma$ error. Top and bottom panels, correspond, respectively, to panels (c) and (f) of those authors’ fig. 1. The data were retrieved from the Fermi-LAT Light Curve Repository (Abdollahi et al. 2023). The photon index, $\Gamma$, is defined such that the flux of photons between energies $\epsilon$ and $\epsilon + d\epsilon$ is proportional to $\epsilon^{-\Gamma}$. This is connected to the $\epsilon F(\epsilon)$ representation, for example, of Fig. 9, in that $\epsilon F(\epsilon) \propto \epsilon^{-\Gamma + 2}$.</figcaption></figure>

Section 4.5. These particles radiate just below pair-production threshold, producing the highest energy observable photons, and they are actively replenished by radiative reprocessing from higher (above-threshold and, hence, invisible) photon energies, which stabilizes the observed spectral energy density. Thus, FSRQ flares in the GeV and TeV bands are ideally suited to probe the expected observational signatures of reconnection in the regime studied in this work.

### 6.1.2 Connection to this work

Let us examine a few instructive flares observed by the Fermi-LAT [which, again, probes the IC(BLR) scenario] from the FSRQ, 3C 279. Hayashida et al. (2012) report 3C 279 flaring periods over the first two years of Fermi operations. For reference, we reproduce the light curve and photon index time-series presented by those authors in their fig. 1, as well as the correlation between the total gamma-ray flux and spectral shape shown in their fig. 3, in our respective Figs 16 and 17 using data retrieved from the Fermi LAT Light Curve Repository (Abdollahi et al. 2023). Hayashida et al. (2012) note mild ‘harder-than-brighter’ behaviour over the entire observation period, but this is somewhat quenched during the brightest periods (Fig. 17), during which the Fermi-measured spectral index becomes remarkably flux-independent (e.g. between MJD 54 800 and 54 900 as well as between MJD 55 000 and 55 100 in Fig. 16). Roughly similar behaviour is often, but not uniformly, seen in later observations of the same object. In two even brighter outbursts from 3C 279 reported by Hayashida et al. (2015) and Ackermann et al. (2016), the gamma-ray flux reached high-enough levels to reconstruct spectra for individual orbits of the Fermi satellite. In the latter event, the flaring individual-orbit spectra reveal a photon energy index that remains between about 1.9 and 2.1 while the flux varies across about a factor of 3 (table 1 of Ackermann et al. 2016). The pre- and post-outburst phases of the event also feature rather stable photon indices (though this appears more statistically significant for the pre-outburst phase; see fig. 1 of Ackermann et al. 2016). As an example of when such spectral stability is not seen, the first flare analysed by Hayashida

<figure><img src="image_18.png" /><figcaption>**Figure 17.** The data from Fig. 16 re-arranged into a scatter plot of gamma-ray spectral index versus photon flux. This figure roughly corresponds to fig. 3, panel (A) of Hayashida et al. (2012). The stabilization of the spectrum during high-flux periods, notably from MJD 54 800 to 54 900 and from MJD 55 000 to 55 100 in Fig. 16, is quite reminiscent of that found in Section 4.5 (cf. Fig. 12).</figcaption></figure>

et al. (2015) exhibits an extreme hardening of the photon index that then softens on the trailing edge of the flare.

To provide a complete view of the correlation between spectral hardness and luminosity of 3C 279 in the GeV band, we also supply Fig. 18. This is similar to Fig. 17 except that the Fermi-LAT spectral index is plotted against the gamma-ray flux level for all 14+ yr of archived data. The full set of Fermi observations demonstrates broad consistency with, for example, the individual flaring period presented by Hayashida et al. (2012) and displayed in Figs 16 and 17: as the object brightens, the variation in photon index appears to decrease, suggesting a stabler spectrum during flaring periods.

At an even more general level, Meyer, Scargle & Blandford (2019) conducted a statistical analysis of the brightest flares from 6 of the most luminous FSRQs detected by Fermi. They find hints of ‘harder-
<figure><img src="image_18.png" /><figcaption>**Figure 18.** The same as Fig. 17, but plotting all (blue background crosses) of the one-week-binned LAT flux levels against corresponding spectral indices for over 14 yr of 3C 279 observations in the *Fermi*-LAT Light Curve Repository (Abdollahi et al. 2023). Additionally, we average the spectral indices (red foreground crosses) in flux bins of width $5 \times 10^7$ photons cm$^{-2}$ s$^{-1}$. This bin width is indicated by horizontal red error bars. Vertical red error bars show $1\sigma$ uncertainty levels after bin averaging.</figcaption></figure>

when-brighter' behaviour in some flares from some objects, but no statistical significance. At the same time, they report that higher flaring flux tends to coincide with reduced spectral variability. These remarks seem to be representative of other individual *Fermi*-detected FSRQ outbursts that are (quasi-)contemporaneous with TeV flares seen by one or more IACTs, including from PKS 1222+216 (Aleksić et al. 2011), PKS 1441+25 (Abeysekara et al. 2015), PKS 1510–089 (Ahnen et al. 2017), PKS 0736+017 (H. E. S. S. Collaboration 2020), and QSO B1420+326 (MAGIC Collaboration 2021). In these events, the GeV (sub-)flare often coincides with a relative stabilization of the GeV spectral index: a seeming decoupling between the flux and the spectral shape – including on the falling part of the flare (albeit the spectral index error bars are larger there). Nevertheless, one does see hints of ‘harder-when-brighter’ trends in some of the GeV (sub-)flares.

Thus, FSRQ flares in the *Fermi*-LAT sensitivity band exhibit broad, though perhaps imperfect, consistency with Klein–Nishina reconnection coupled to BLR Compton seed photons. The trends seem to persist even into the flare decay: if radiative losses reverted entirely to the Thomson regime, this would induce a strong (and unobserved) correlation between the gamma-ray spectral index and the flux in the decaying part of the flares. This suggests the role of Klein–Nishina effects, as unveiled in this work, in stabilizing the spectrum even as the light curve declines. For the cases where a harder-when-brighter trend appears more evident, this could be induced by weak synchrotron losses perturbing the dominant temporal-spectral behaviour induced by Klein–Nishina and pair-production physics.

One caveat to our association of GeV FSRQ flares with IC(BLR) Klein–Nishina reconnection is the occasional (quasi-)contemporaneous detection of TeV gamma-rays (a few examples of which are listed two paragraphs prior). For these cases, a simple one-zone emission model for both the GeV and TeV outburst places the emission region outside the BLR, which would otherwise absorb the TeV photons [equation (56)]. This caveat is made more severe by population studies of *Fermi*-detected FSRQs purely in the GeV band, which find no evidence for gamma-ray absorption by the BLR seed photons in the vast majority of objects (Costamante et al. 2018; Meyer et al. 2019). As noted by Costamante et al. (2018), however,

even within a single zone framework, these constraints become less severe if one attributes the flaring emission to a structure (in our case, a reconnection layer) that comoves with the jet rather than a stationary feature (e.g. a standing shock). Then, over an observing period, $\Delta t_{\text{obs}}$, of just 1 d (typically comparable to or shorter than variability time-scales identified by TeV FSRQ observations, with two notable exceptions: Aleksić et al. 2011 and H. E. S. S. Collaboration 2021), Doppler time compression of the light curve allows the emitting zone to travel a distance, $d \sim \Gamma_j^2 c \Delta t_{\text{obs}} \sim 10^{17}$ cm, from the central engine for a fiducial jet Lorentz factor, $\Gamma_j = 10$. This is at the edge of the BLR, which only extends up to roughly 0.1 pc $\sim 10^{17}$ cm from the nucleus (Tavecchio & Ghisellini 2008; Sikora et al. 2009; Nalewajko et al. 2012; Mehlhaff et al. 2021), reducing the importance of absorption for the (potentially up to TeV) part of the emission produced at larger distances. In this view, attributing emission from longer flaring periods (e.g. as in Fig. 16) to IC(BLR)-coupled reconnection demands positing the ejection of multiple reconnecting structures from the central engine, as, for example, in striped-jet models (Giannios & Uzdensky 2019).

We next consider FSRQ flares at TeV energies observed by IACTs. Such events probe the the IC(HDR) scenario wherein reconnection couples to Klein–Nishina and pair-production physics mediated by the HDR seed photons. Suitable observations are much more difficult to obtain here. Not only do FSRQs shine intrinsically weakly in the TeV band (recall that their quiescent IC spectral hump peaks already at MeV or, sometimes, GeV energies), but they suffer both intrinsic and external absorption at these energies. Intrinsically, TeV emission produced inside the BLR will be absorbed, rendering invisible TeV flaring regions too close to the central engine. Furthermore, because FSRQs (unlike BL Lacs) are distributed in the Universe preferentially at higher redshift (Ackermann et al. 2015; Ajello et al. 2020), their TeV gamma-rays may be absorbed by the extragalactic background light while en route toward Earth. Owing to these combined effects, only a handful of FSRQs have even been detected at TeV energies (9 at the time of writing, Wakely & Horan 2008). Of these, quiescent emission has only been seen from one object, PKS 1510–089 (MAGIC Collaboration 2018). The rest are detected exclusively in high or flaring states (e.g. 3C 279, MAGIC Collaboration 2008; PKS 1222+216, Aleksić et al. 2011; PKS 1441+25, Abeysekara et al. 2015; QSO B0218+257, Sitarek et al. 2015; Ahnen et al. 2016; PKS 0736+017, H. E. S. S. Collaboration 2020; QSO B1420+326, MAGIC Collaboration 2021). Even during flares, the gamma-ray flux is rarely sufficient to provide detailed temporal information for the total luminosity itself, much less for the spectrum.

Considering mostly time-integrated spectra, what is generally seen for FSRQs in the TeV band is a much steeper spectrum, even after correcting for absorption by the extragalactic background light, than in the *Fermi* range. TeV spectral indices are almost always greater than 2.5 and commonly larger than 3: for 3C 279, $\Gamma \simeq 4.2$ (H. E. S. S. Collaboration 2019); for PKS 1510–089, $\Gamma \simeq 2.5$ (Aleksić et al. 2014), $\Gamma \simeq 3.2$ and 4.3 (Ahnen et al. 2017), $\Gamma \simeq 2.9$ and 3.4 (Zacharias et al. 2017), and $\Gamma \simeq 3.3$ (MAGIC Collaboration 2018); for PKS 1222+216, $\Gamma \simeq 2.7$ (Aleksić et al. 2011); for PKS 1441+25, $\Gamma \simeq 3.4$ (Abeysekara et al. 2015); for QSO B0218+257, $\Gamma \simeq 2.4$ (Ahnen et al. 2016); and for QSO B1420+326, $\Gamma \simeq 2.9$ (MAGIC Collaboration 2021). In the IC(HDR) scenario, these steep spectral indices suggest that reconnection proceeds in a regime, unlike that probed in detail in this study, where its intrinsic particle acceleration index is steeper (e.g. with a strong guide field: Werner & Uzdensky 2017). Even in this case, Klein–Nishina radiative physics should still harden the spectrum in the decaying part of a flare (Section 4.5). However, such detailed temporal behaviour has not
<figure><img src="image_19.png" /><figcaption>**Figure 18.** The same as Fig. 17, but plotting all (blue background crosses) of the one-week-binned LAT flux levels against corresponding spectral indices for over 14 yr of 3C 279 observations in the *Fermi*–LAT Light Curve Repository (Abdollahi et al. 2023). Additionally, we average the spectral indices (red foreground crosses) in flux bins of width $5 \times 10^7$ photons cm$^{-2}$ s$^{-1}$. This bin width is indicated by horizontal red error bars. Vertical red error bars show $1\sigma$ uncertainty levels after bin averaging.</figcaption></figure>

when-brighter' behaviour in some flares from some objects, but no statistical significance. At the same time, they report that higher flaring flux tends to coincide with reduced spectral variability. These remarks seem to be representative of other individual *Fermi*-detected FSRQ outbursts that are (quasi-)contemporaneous with TeV flares seen by one or more IACTs, including from PKS 1222+216 (Aleksić et al. 2011), PKS 1441+25 (Abeysekara et al. 2015), PKS 1510–089 (Ahnen et al. 2017), PKS 0736+017 (H. E. S. S. Collaboration 2020), and QSO B1420+326 (MAGIC Collaboration 2021). In these events, the GeV (sub-)flare often coincides with a relative stabilization of the GeV spectral index: a seeming decoupling between the flux and the spectral shape – including on the falling part of the flare (albeit the spectral index error bars are larger there). Nevertheless, one does see hints of ‘harder-when-brighter’ trends in some of the GeV (sub-)flares.

Thus, FSRQ flares in the *Fermi*–LAT sensitivity band exhibit broad, though perhaps imperfect, consistency with Klein–Nishina reconnection coupled to BLR Compton seed photons. The trends seem to persist even into the flare decay: if radiative losses reverted entirely to the Thomson regime, this would induce a strong (and unobserved) correlation between the gamma-ray spectral index and the flux in the decaying part of the flares. This suggests the role of Klein–Nishina effects, as unveiled in this work, in stabilizing the spectrum even as the light curve declines. For the cases where a harder-when-brighter trend appears more evident, this could be induced by weak synchrotron losses perturbing the dominant temporal-spectral behaviour induced by Klein–Nishina and pair-production physics.

One caveat to our association of GeV FSRQ flares with IC(BLR) Klein–Nishina reconnection is the occasional (quasi-)contemporaneous detection of TeV gamma-rays (a few examples of which are listed two paragraphs prior). For these cases, a simple one-zone emission model for both the GeV and TeV outburst places the emission region outside the BLR, which would otherwise absorb the TeV photons [equation (56)]. This caveat is made more severe by population studies of *Fermi*-detected FSRQs purely in the GeV band, which find no evidence for gamma-ray absorption by the BLR seed photons in the vast majority of objects (Costamante et al. 2018; Meyer et al. 2019). As noted by Costamante et al. (2018), however,

even within a single zone framework, these constraints become less severe if one attributes the flaring emission to a structure (in our case, a reconnection layer) that comoves with the jet rather than a stationary feature (e.g. a standing shock). Then, over an observing period, $\Delta t_{\text{obs}}$, of just 1 d (typically comparable to or shorter than variability time-scales identified by TeV FSRQ observations, with two notable exceptions: Aleksić et al. 2011 and H. E. S. S. Collaboration 2021), Doppler time compression of the light curve allows the emitting zone to travel a distance, $d \sim \Gamma_j^2 c \Delta t_{\text{obs}} \sim 10^{17}$ cm, from the central engine for a fiducial jet Lorentz factor, $\Gamma_j = 10$. This is at the edge of the BLR, which only extends up to roughly 0.1 pc $\sim 10^{17}$ cm from the nucleus (Tavecchio & Ghisellini 2008; Sikora et al. 2009; Nalewajko et al. 2012; Mehlhaff et al. 2021), reducing the importance of absorption for the (potentially up to TeV) part of the emission produced at larger distances. In this view, attributing emission from longer flaring periods (e.g. as in Fig. 16) to IC(BLR)-coupled reconnection demands positing the ejection of multiple reconnecting structures from the central engine, as, for example, in striped-jet models (Giannios & Uzdensky 2019).

We next consider FSRQ flares at TeV energies observed by IACTs. Such events probe the the IC(HDR) scenario wherein reconnection couples to Klein–Nishina and pair-production physics mediated by the HDR seed photons. Suitable observations are much more difficult to obtain here. Not only do FSRQs shine intrinsically weakly in the TeV band (recall that their quiescent IC spectral hump peaks already at MeV or, sometimes, GeV energies), but they suffer both intrinsic and external absorption at these energies. Intrinsically, TeV emission produced inside the BLR will be absorbed, rendering invisible TeV flaring regions too close to the central engine. Furthermore, because FSRQs (unlike BL Lacs) are distributed in the Universe preferentially at higher redshift (Ackermann et al. 2015; Ajello et al. 2020), their TeV gamma-rays may be absorbed by the extragalactic background light while en route toward Earth. Owing to these combined effects, only a handful of FSRQs have even been detected at TeV energies (9 at the time of writing, Wakely & Horan 2008). Of these, quiescent emission has only been seen from one object, PKS 1510–089 (MAGIC Collaboration 2018). The rest are detected exclusively in high or flaring states (e.g. 3C 279, MAGIC Collaboration 2008; PKS 1222+216, Aleksić et al. 2011; PKS 1441+25, Abeysekara et al. 2015; QSO B0218+257, Sitarek et al. 2015; Ahnen et al. 2016; PKS 0736+017, H. E. S. S. Collaboration 2020; QSO B1420+326, MAGIC Collaboration 2021). Even during flares, the gamma-ray flux is rarely sufficient to provide detailed temporal information for the total luminosity itself, much less for the spectrum.

Considering mostly time-integrated spectra, what is generally seen for FSRQs in the TeV band is a much steeper spectrum, even after correcting for absorption by the extragalactic background light, than in the *Fermi* range. TeV spectral indices are almost always greater than 2.5 and commonly larger than 3: for 3C 279, $\Gamma \simeq 4.2$ (H. E. S. S. Collaboration 2019); for PKS 1510–089, $\Gamma \simeq 2.5$ (Aleksić et al. 2014), $\Gamma \simeq 3.2$ and 4.3 (Ahnen et al. 2017), $\Gamma \simeq 2.9$ and 3.4 (Zacharias et al. 2017), and $\Gamma \simeq 3.3$ (MAGIC Collaboration 2018); for PKS 1222+216, $\Gamma \simeq 2.7$ (Aleksić et al. 2011); for PKS 1441+25, $\Gamma \simeq 3.4$ (Abeysekara et al. 2015); for QSO B0218+257, $\Gamma \simeq 2.4$ (Ahnen et al. 2016); and for QSO B1420+326, $\Gamma \simeq 2.9$ (MAGIC Collaboration 2021). In the IC(HDR) scenario, these steep spectral indices suggest that reconnection proceeds in a regime, unlike that probed in detail in this study, where its intrinsic particle acceleration index is steeper (e.g. with a strong guide field: Werner & Uzdensky 2017). Even in this case, Klein–Nishina radiative physics should still harden the spectrum in the decaying part of a flare (Section 4.5). However, such detailed temporal behaviour has not
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_{KN}$ particle is $\gamma_{KN, disc}m_e c^2/2 = \epsilon_{th, disc}/8 = 30 \text{ MeV}$, where $\gamma_{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_{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_{KN, disc}m_e c^2$, where Klein–Nishina effects are likely most prominent, are not promising. At lower, X-ray energies,