# Kinetic simulations and gamma-ray signatures of Klein-Nishina relativistic magnetic reconnection

J Mehlhaff, G Werner, B Cerutti, D Uzdensky, M Begelman

► To cite this version:

J Mehlhaff, G Werner, B Cerutti, D Uzdensky, M Begelman. Kinetic simulations and gamma-ray signatures of Klein-Nishina relativistic magnetic reconnection. Monthly Notices of the Royal Astronomical Society, 2023, 527 (4), pp.11587-11626. 10.1093/mnras/stad3863. hal-04160424

HAL Id: hal-04160424

https://hal.science/hal-04160424

Submitted on 21 Apr 2024

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
ROYAL ASTRONOMICAL SOCIETY

MNRAS 527, 11587–11626 (2024)

https://doi.org/10.1093/mnras/stad3863

Advance Access publication 2023 December 20

# Kinetic simulations and gamma-ray signatures of Klein–Nishina relativistic magnetic reconnection

J. Mehlhaff<sup>1,1,2*</sup> G. Werner<sup>2</sup> B. Cerutti<sup>1</sup> D. Uzdensky<sup>2</sup> and M. Begelman<sup>3,4</sup>

Accepted 2023 December 12. Received 2023 December 11; in original form 2023 June 29

## ABSTRACT

Black hole and neutron star environments often comprise collisionless plasmas immersed in strong magnetic fields and intense baths of low-frequency radiation. In such conditions, relativistic magnetic reconnection can tap the magnetic field energy, accelerating high-energy particles that rapidly cool by inverse Compton (IC) scattering the dense photon background. At the highest particle energies reached in bright gamma-ray sources, IC scattering can stray into the Klein–Nishina regime. Here, the Comptonized photons exceed pair-production threshold with the radiation background and may thus return their energy to the reconnecting plasma as fresh electron–positron pairs. To reliably characterize observable signatures of such Klein–Nishina reconnection, in this work, we present first-principles particle-in-cell simulations of pair-plasma relativistic reconnection coupled to Klein–Nishina and pair-production physics. The simulations show substantial differences between the observable signatures of Klein–Nishina reconnection and reconnection coupled only to low-energy Thomson IC cooling (without pair production). The latter regime exhibits strong harder-when-brighter behaviour; the former involves a stable spectral shape independent of overall brightness. This spectral stability is reminiscent of flat-spectrum radio quasar (FSRQ) GeV high states, furnishing evidence that Klein–Nishina radiative physics operates in FSRQs. The simulated Klein–Nishina reconnection pair yield spans from low to order-unity and follows an exponential scaling law in a single governing parameter. Pushing this parameter beyond its range studied here might give way to a copious pair-creation regime. Besides FSRQs, we discuss potential applications to accreting black hole X-ray binaries, the M87* magnetosphere, and gamma-ray binaries.

**Key words:** acceleration of particles – magnetic reconnection – radiation mechanisms: general – relativistic processes – gamma-rays: general.

## 1 INTRODUCTION

The gamma-ray sky is studded with relativistic compact objects – neutron stars and black holes (of which the most numerous observed varieties are, respectively, pulsars and blazars: Wakely & Horan 2008; Abdollahi et al. 2020). These systems – and connected phenomena including winds, jets, and accretion discs – frequently host collisionless highly magnetized plasmas, with magnetic energy density exceeding not just the pressure (i.e. small plasma beta) but also the rest-mass energy density of the charge-carrying particles. In such plasmas, relativistic magnetic reconnection (Blackman & Field 1994; Lyutikov & Uzdensky 2003; Lyubarsky 2005; Watanabe & Yokoyama 2006) efficiently siphons off the excess magnetic field energy, using it to accelerate relativistic particles and drive relativistic collective motion. The energized particles are then revealed by the light that they shine toward Earth, including in the gamma-ray band. Relativistic magnetic reconnection is, hence, an important candidate

mechanism for powering high-energy phenomena linked to the most compact objects in the Universe.

In some astrophysical situations, there is a well-defined time-scale separation between abrupt reconnection-powered particle acceleration and much slower radiative losses. In this radiatively inefficient regime, observable emission traces particle energization that has occurred in the past. This limit is seldom realized, however, in the plasma environments of compact objects, where intense magnetic and radiation fields lead to rapid synchrotron and inverse Compton (IC) cooling. Then, the problem can no longer be cleanly factorized into a sudden acceleration step followed by a more prolonged cooling stage. Instead, radiative cooling couples in real time to reconnection, tracing active (as opposed to past) particle acceleration, and feeding back on the reconnection process: a qualitatively distinct radiative reconnection regime (Uzdensky 2011, 2016; Mehlhaff et al. 2021).

In the low-energy, optically thin limit where the synchrotron and IC photons freely escape the system – what we might call classical radiative reconnection – radiative losses, while dynamically important, do not change the fundamental flow of energy from the non-radiative case. Magnetic fields serve as the main energy source

<footnote n="1"><sup>1</sup>Univ. Grenoble Alpes, CNRS, IPAG, 38000 Grenoble, France</footnote>

<footnote n="2"><sup>2</sup>Center for Integrated Plasma Studies, Physics Department, University of Colorado, 390 UCB, Boulder, CO 80309, USA</footnote>

<footnote n="3"><sup>3</sup>JILA, University of Colorado and National Institute of Standards and Technology, 440 UCB, Boulder, CO 80309, USA</footnote>

<footnote n="4"><sup>4</sup>Department of Astrophysical and Planetary Sciences, University of Colorado, 391 UCB, Boulder, CO 80309, USA</footnote>

<footnote>*E-mail: john.mehlhaff@univ-grenoble-alpes.fr</footnote>
ROYAL ASTRONOMICAL SOCIETY

MNRAS 527, 11587–11626 (2024)

https://doi.org/10.1093/mnras/stad3863

Advance Access publication 2023 December 20

# Kinetic simulations and gamma-ray signatures of Klein–Nishina relativistic magnetic reconnection

J. Mehlhaff<sup>1,1,2*</sup> G. Werner<sup>2</sup> B. Cerutti<sup>1</sup> D. Uzdensky<sup>2</sup> and M. Begelman<sup>3,4</sup>

Accepted 2023 December 12. Received 2023 December 11; in original form 2023 June 29

## ABSTRACT

Black hole and neutron star environments often comprise collisionless plasmas immersed in strong magnetic fields and intense baths of low-frequency radiation. In such conditions, relativistic magnetic reconnection can tap the magnetic field energy, accelerating high-energy particles that rapidly cool by inverse Compton (IC) scattering the dense photon background. At the highest particle energies reached in bright gamma-ray sources, IC scattering can stray into the Klein–Nishina regime. Here, the Comptonized photons exceed pair-production threshold with the radiation background and may thus return their energy to the reconnecting plasma as fresh electron–positron pairs. To reliably characterize observable signatures of such Klein–Nishina reconnection, in this work, we present first-principles particle-in-cell simulations of pair-plasma relativistic reconnection coupled to Klein–Nishina and pair-production physics. The simulations show substantial differences between the observable signatures of Klein–Nishina reconnection and reconnection coupled only to low-energy Thomson IC cooling (without pair production). The latter regime exhibits strong harder-when-brighter behaviour; the former involves a stable spectral shape independent of overall brightness. This spectral stability is reminiscent of flat-spectrum radio quasar (FSRQ) GeV high states, furnishing evidence that Klein–Nishina radiative physics operates in FSRQs. The simulated Klein–Nishina reconnection pair yield spans from low to order-unity and follows an exponential scaling law in a single governing parameter. Pushing this parameter beyond its range studied here might give way to a copious pair-creation regime. Besides FSRQs, we discuss potential applications to accreting black hole X-ray binaries, the M87* magnetosphere, and gamma-ray binaries.

**Key words:** acceleration of particles – magnetic reconnection – radiation mechanisms: general – relativistic processes – gamma-rays: general.

## 1 INTRODUCTION

The gamma-ray sky is studded with relativistic compact objects – neutron stars and black holes (of which the most numerous observed varieties are, respectively, pulsars and blazars: Wakely & Horan 2008; Abdollahi et al. 2020). These systems – and connected phenomena including winds, jets, and accretion discs – frequently host collisionless highly magnetized plasmas, with magnetic energy density exceeding not just the pressure (i.e. small plasma beta) but also the rest-mass energy density of the charge-carrying particles. In such plasmas, relativistic magnetic reconnection (Blackman & Field 1994; Lyutikov & Uzdensky 2003; Lyubarsky 2005; Watanabe & Yokoyama 2006) efficiently siphons off the excess magnetic field energy, using it to accelerate relativistic particles and drive relativistic collective motion. The energized particles are then revealed by the light that they shine toward Earth, including in the gamma-ray band. Relativistic magnetic reconnection is, hence, an important candidate

mechanism for powering high-energy phenomena linked to the most compact objects in the Universe.

In some astrophysical situations, there is a well-defined time-scale separation between abrupt reconnection-powered particle acceleration and much slower radiative losses. In this radiatively inefficient regime, observable emission traces particle energization that has occurred in the past. This limit is seldom realized, however, in the plasma environments of compact objects, where intense magnetic and radiation fields lead to rapid synchrotron and inverse Compton (IC) cooling. Then, the problem can no longer be cleanly factorized into a sudden acceleration step followed by a more prolonged cooling stage. Instead, radiative cooling couples in real time to reconnection, tracing active (as opposed to past) particle acceleration, and feeding back on the reconnection process: a qualitatively distinct radiative reconnection regime (Uzdensky 2011, 2016; Mehlhaff et al. 2021).

In the low-energy, optically thin limit where the synchrotron and IC photons freely escape the system – what we might call classical radiative reconnection – radiative losses, while dynamically important, do not change the fundamental flow of energy from the non-radiative case. Magnetic fields serve as the main energy source

<footnote n="1"><sup>1</sup>Univ. Grenoble Alpes, CNRS, IPAG, 38000 Grenoble, France</footnote>

<footnote n="2"><sup>2</sup>Center for Integrated Plasma Studies, Physics Department, University of Colorado, 390 UCB, Boulder, CO 80309, USA</footnote>

<footnote n="3"><sup>3</sup>JILA, University of Colorado and National Institute of Standards and Technology, 440 UCB, Boulder, CO 80309, USA</footnote>

<footnote n="4"><sup>4</sup>Department of Astrophysical and Planetary Sciences, University of Colorado, 391 UCB, Boulder, CO 80309, USA</footnote>

<footnote>*E-mail: john.mehlhaff@univ-grenoble-alpes.fr</footnote>
for particle acceleration, while the primary energy sink is the emission mechanisms that efficiently and permanently remove liberated magnetic energy from the system. However, for the gamma-ray-bright relativistic compact objects, the photons emitted at the highest energies are above threshold for various quantum electrodynamical (QED) pair production channels. This fundamentally alters the pathways available to the energy in radiative reconnection, allowing radiation to not just carry energy away from the reconnection site, but also to redistribute it in real time in the form of freshly produced electron-positron pairs. Such QED *radiative reconnection* is thus distinguished from merely classical radiative reconnection in its capacity to alter the plasma material composition and in the key role played by photons as a dynamically important particle species (Uzdensky 2011, 2016; Uzdensky et al. 2019).

To interpret observations of systems where QED reconnection may occur, modelling efforts must employ a self-consistent kinetic plasma description. Such a description is already needed to model collisionless relativistic magnetic reconnection in the non-radiative and classical radiative regimes, because it captures the critical microphysics governing the reconnection rate (i.e. the rate of magnetic energy dissipation) as well as the production of non-thermal particle energy distributions and correspondingly non-thermal emission spectra. In the case of QED reconnection, a kinetic paradigm is even more imperative. The QED cross-sections depend sensitively on the energies of both the emitting particles and the pair-producing photons, placing an even greater importance on capturing energization self-consistently.

All of the necessary kinetic physics can be incorporated by augmenting *ab initio* particle-in-cell (PIC) simulations (Hockney & Eastwood 1988; Birdsall & Langdon 1991) with QED physics. The small number of PIC reconnection studies that have done this have focused on a select few QED interactions. For example, Hakobyan, Philippov & Spitkovsky (2019) present a regime, expected in pulsar magnetospheres, where particles suffer strong synchrotron cooling, and the resulting synchrotron photons, with an emission spectrum peaking at $\gtrsim$ MeV energies, collide with one another, leading to copious pair production in the reconnection inflow region (see also Hakobyan, Ripperda & Philippov 2023b). The same emission and pair-production channels are coupled to reconnection (with an additional IC post-processing step) by Xie et al. (2023) in the context of magnetar magnetospheres. Reconnection in magnetar magnetospheres is also studied by Schoeffler et al. (2019, 2023). However, those authors consider pair production not between colliding synchrotron photons, but from the absorption of single synchrotron photons by the intense magnetar magnetic field. They observe the radiative cooling to remove particle pressure support from the reconnection layer, leading to strong plasma compression there. This locally amplifies synchrotron emission and pair creation, with the end result that pair production in this regime is concentrated not in the plasma fuelling reconnection, but in the heart of the reconnection layer itself. Lastly, Crinquand et al. (2021, 2022) present global models of reconnection in black hole magnetospheres, where the primary radiation mechanism is IC scattering of low-energy (soft) background photons originating from a larger-scale accretion flow. Pair production then occurs when Comptonized photons collide with the soft parent population from which they were first scattered, supplying the plasma to a luminous equatorial reconnection current sheet in the black hole magnetosphere.

The examples above illustrate two general points. First, while QED reconnection is of general high-energy astrophysical interest, the relevant QED interactions depend on the system under study. Second, specializing to certain QED interactions over others not

only decides the applicable astrophysical sources, but can also lead to divergent qualitative dynamics. These remarks underscore the need to understand QED reconnection – in all its astrophysical diversity – as a fundamental physics problem in order to identify its observable signatures in the high-energy Universe.

In this work, we present PIC simulations run using the ZELTRON code (Cerutti et al. 2013; Cerutti & Werner 2019, plus needed auxiliary developments detailed here) of a QED reconnection regime thus far unexplored from first principles as an isolated physics problem. We consider a relativistic magnetic reconnection layer immersed in such an intense bath of soft background radiation that IC scattering strays far into the QED limit. This contrasts the (classical radiative) low-energy, Thomson IC reconnection regime previously studied numerically by Werner, Philippov & Uzdensky (2019), Mehlhaff et al. (2020), Sironi & Beloborodov (2020), and Sridhar, Sironi & Beloborodov (2021, 2023) in two respects. First, we account for quantized gamma-ray emission from the highest energy particles radiating in the Klein–Nishina IC limit (Jones 1968; Blumenthal & Gould 1970). Second, because many of the emitted photons lie above pair-production threshold with the soft background (Mehlhaff et al. 2021), we calculate pair production between the few (low number density) scattered gamma-rays with energies $\epsilon_{scat} \gg m_e c^2$, and the abundant (high number density) soft seed photons with energies $\epsilon_{ph} \ll m_e c^2$. We refer to this realization of QED reconnection as *Klein–Nishina radiative reconnection* (sometimes just *Klein–Nishina reconnection*), omitting explicit reference to pair production since efficient Klein–Nishina IC emission implies pair production in a reconnection context (Mehlhaff et al. 2021).

The QED interactions studied here are the same as those treated by Crinquand et al. (2021, 2022). Here, however, we take a complementary approach, stripping away the global morphology and studying reconnection as a local problem. Computationally, this enables us to concentrate resources toward enhancing the separation among the radiative and plasma microscales. Physically, it permits us to remain more agnostic to the host system, focusing instead on the intrinsic reconnection properties that may be generally applicable.

Indeed, Klein–Nishina reconnection may be realized in a range of astrophysical systems connected to relativistic compact objects, including: flat-spectrum radio quasars (FSRQs), where reconnection occurring in a relativistic jet launched from an active galactic nucleus (AGN) is likely externally illuminated by large-scale circumnuclear structures (cf. Giannios, Uzdensky & Begelman 2009; Nalewajko et al. 2011, 2012; Giannios 2013; Petropoulou, Giannios & Sironi 2016; Sironi, Giannios & Petropoulou 2016; Werner et al. 2018; Christie et al. 2019, 2020; Giannios & Uzdensky 2019; Mehlhaff et al. 2020, 2021; Ortuño-Macías & Nalewajko 2020); the high/soft states of accreting black hole X-ray binaries (BHXRBs), where reconnection in a highly magnetized collisionless coronal region is illuminated by an underlying geometrically thin, optically thick accretion disc (cf. Galeev, Rosner & Vaiana 1979; Di Matteo 1998; Goodman & Uzdensky 2008; Uzdensky & Goodman 2008; Hoshino & Lyubarsky 2012; Uzdensky 2016; Beloborodov 2017; Werne et al. 2019; Sironi & Beloborodov 2020; Mehlhaff et al. 2021; Sridhar et al. 2021, 2023; El Mellah et al. 2022); the magnetospheres of supermassive black holes, particularly the one at the centre of the M87 galaxy, M87*, wherein reconnection may be bathed in photons from a large-scale radiatively inefficient accretion flow (cf. Giannios, Uzdensky & Begelman 2010; Ball et al. 2016; Li, Yuan & Wang 2017; Ripperda, Bacchini & Philippov 2020; Crinquand et al. 2021, 2022; El Mellah et al. 2022; Ripperda et al. 2022; Scepi, Dexter & Begelman 2022; Chen, Uzdensky & Dexter 2023; El Mellah, Cerutti & Crinquand 2023; Galishnikova et al. 2023;
for particle acceleration, while the primary energy sink is the emission mechanisms that efficiently and permanently remove liberated magnetic energy from the system. However, for the gamma-ray-bright relativistic compact objects, the photons emitted at the highest energies are above threshold for various quantum electrodynamical (QED) pair production channels. This fundamentally alters the pathways available to the energy in radiative reconnection, allowing radiation to not just carry energy away from the reconnection site, but also to redistribute it in real time in the form of freshly produced electron-positron pairs. Such QED *radiative reconnection* is thus distinguished from merely classical radiative reconnection in its capacity to alter the plasma material composition and in the key role played by photons as a dynamically important particle species (Uzdensky 2011, 2016; Uzdensky et al. 2019).

To interpret observations of systems where QED reconnection may occur, modelling efforts must employ a self-consistent kinetic plasma description. Such a description is already needed to model collisionless relativistic magnetic reconnection in the non-radiative and classical radiative regimes, because it captures the critical microphysics governing the reconnection rate (i.e. the rate of magnetic energy dissipation) as well as the production of non-thermal particle energy distributions and correspondingly non-thermal emission spectra. In the case of QED reconnection, a kinetic paradigm is even more imperative. The QED cross-sections depend sensitively on the energies of both the emitting particles and the pair-producing photons, placing an even greater importance on capturing energization self-consistently.

All of the necessary kinetic physics can be incorporated by augmenting *ab initio* particle-in-cell (PIC) simulations (Hockney & Eastwood 1988; Birdsall & Langdon 1991) with QED physics. The small number of PIC reconnection studies that have done this have focused on a select few QED interactions. For example, Hakobyan, Philippov & Spitkovsky (2019) present a regime, expected in pulsar magnetospheres, where particles suffer strong synchrotron cooling, and the resulting synchrotron photons, with an emission spectrum peaking at $\gtrsim$ MeV energies, collide with one another, leading to copious pair production in the reconnection inflow region (see also Hakobyan, Ripperda & Philippov 2023b). The same emission and pair-production channels are coupled to reconnection (with an additional IC post-processing step) by Xie et al. (2023) in the context of magnetar magnetospheres. Reconnection in magnetar magnetospheres is also studied by Schoeffler et al. (2019, 2023). However, those authors consider pair production not between colliding synchrotron photons, but from the absorption of single synchrotron photons by the intense magnetar magnetic field. They observe the radiative cooling to remove particle pressure support from the reconnection layer, leading to strong plasma compression there. This locally amplifies synchrotron emission and pair creation, with the end result that pair production in this regime is concentrated not in the plasma fuelling reconnection, but in the heart of the reconnection layer itself. Lastly, Crinquand et al. (2021, 2022) present global models of reconnection in black hole magnetospheres, where the primary radiation mechanism is IC scattering of low-energy (soft) background photons originating from a larger-scale accretion flow. Pair production then occurs when Comptonized photons collide with the soft parent population from which they were first scattered, supplying the plasma to a luminous equatorial reconnection current sheet in the black hole magnetosphere.

The examples above illustrate two general points. First, while QED reconnection is of general high-energy astrophysical interest, the relevant QED interactions depend on the system under study. Second, specializing to certain QED interactions over others not

only decides the applicable astrophysical sources, but can also lead to divergent qualitative dynamics. These remarks underscore the need to understand QED reconnection – in all its astrophysical diversity – as a fundamental physics problem in order to identify its observable signatures in the high-energy Universe.

In this work, we present PIC simulations run using the ZELTRON code (Cerutti et al. 2013; Cerutti & Werner 2019, plus needed auxiliary developments detailed here) of a QED reconnection regime thus far unexplored from first principles as an isolated physics problem. We consider a relativistic magnetic reconnection layer immersed in such an intense bath of soft background radiation that IC scattering strays far into the QED limit. This contrasts the (classical radiative) low-energy, Thomson IC reconnection regime previously studied numerically by Werner, Philippov & Uzdensky (2019), Mehlhaff et al. (2020), Sironi & Beloborodov (2020), and Sridhar, Sironi & Beloborodov (2021, 2023) in two respects. First, we account for quantized gamma-ray emission from the highest energy particles radiating in the Klein–Nishina IC limit (Jones 1968; Blumenthal & Gould 1970). Second, because many of the emitted photons lie above pair-production threshold with the soft background (Mehlhaff et al. 2021), we calculate pair production between the few (low number density) scattered gamma-rays with energies $\epsilon_{scat} \gg m_e c^2$, and the abundant (high number density) soft seed photons with energies $\epsilon_{ph} \ll m_e c^2$. We refer to this realization of QED reconnection as *Klein–Nishina radiative reconnection* (sometimes just *Klein–Nishina reconnection*), omitting explicit reference to pair production since efficient Klein–Nishina IC emission implies pair production in a reconnection context (Mehlhaff et al. 2021).

The QED interactions studied here are the same as those treated by Crinquand et al. (2021, 2022). Here, however, we take a complementary approach, stripping away the global morphology and studying reconnection as a local problem. Computationally, this enables us to concentrate resources toward enhancing the separation among the radiative and plasma microscales. Physically, it permits us to remain more agnostic to the host system, focusing instead on the intrinsic reconnection properties that may be generally applicable.

Indeed, Klein–Nishina reconnection may be realized in a range of astrophysical systems connected to relativistic compact objects, including: flat-spectrum radio quasars (FSRQs), where reconnection occurring in a relativistic jet launched from an active galactic nucleus (AGN) is likely externally illuminated by large-scale circumnuclear structures (cf. Giannios, Uzdensky & Begelman 2009; Nalewajko et al. 2011, 2012; Giannios 2013; Petropoulou, Giannios & Sironi 2016; Sironi, Giannios & Petropoulou 2016; Werner et al. 2018; Christie et al. 2019, 2020; Giannios & Uzdensky 2019; Mehlhaff et al. 2020, 2021; Ortuño-Macías & Nalewajko 2020); the high/soft states of accreting black hole X-ray binaries (BHXRBs), where reconnection in a highly magnetized collisionless coronal region is illuminated by an underlying geometrically thin, optically thick accretion disc (cf. Galeev, Rosner & Vaiana 1979; Di Matteo 1998; Goodman & Uzdensky 2008; Uzdensky & Goodman 2008; Hoshino & Lyubarsky 2012; Uzdensky 2016; Beloborodov 2017; Werne et al. 2019; Sironi & Beloborodov 2020; Mehlhaff et al. 2021; Sridhar et al. 2021, 2023; El Mellah et al. 2022); the magnetospheres of supermassive black holes, particularly the one at the centre of the M87 galaxy, M87*, wherein reconnection may be bathed in photons from a large-scale radiatively inefficient accretion flow (cf. Giannios, Uzdensky & Begelman 2010; Ball et al. 2016; Li, Yuan & Wang 2017; Ripperda, Bacchini & Philippov 2020; Crinquand et al. 2021, 2022; El Mellah et al. 2022; Ripperda et al. 2022; Scepi, Dexter & Begelman 2022; Chen, Uzdensky & Dexter 2023; El Mellah, Cerutti & Crinquand 2023; Galishnikova et al. 2023;
Hakobyan et al. 2023b); and gamma-ray binaries, where a plausible scenario involves a pulsar in tight orbit around a bright type O or Be star, which illuminates reconnection occurring near the pulsar in its magnetosphere and striped wind (i.e. before the pulsar wind shocks with the stellar wind from the companion; cf. Dubus 2006; Cerutti, Dubus & Henri 2008; Dubus 2013; Cerutti & Philippov 2017; Dubus et al. 2017; Philippov & Spitkovsky 2018; Cerutti, Philippov & Dubus 2020). The link between Klein–Nishina reconnection and each of these object classes is, in fact, a major result of this work (Section 6), as further discussed below.

While this study is primarily numerical, analytic, and semi-
analytic modelling are also vital for understanding Klein–Nishina
reconnection and QED reconnection more broadly. Such theoretical
approaches can make targeted, physically motivated arguments for
how results from non-QED reconnection may generalize to the
QED case (e.g. Beloborodov 2017; Mehlhaff et al. 2021; Chen et
al. 2023; Hakobyan et al. 2023b), even if they cannot treat all of
the kinetic physics at play from first principles. This furnishes a
useful interpretive framework for *ab initio* simulations. However,
the reverse is also true: phenomenological models, which sometimes
have the advantage of enhanced physical clarity, can themselves
be refined from the findings of simulations. This study illustrates
both directions of this paradigm. Throughout the text, we make
frequent reference to our earlier work, Mehlhaff et al. (2021), which
analytically considers the setup simulated here. As will be seen,
that study (besides laying much of the theoretical foundation for
this article) serves both as an interpretive lens for our simulations
and as a set of hypotheses that the numerical experiments can
check.

We structure this article as follows. In Section 2, we detail the QED algorithmic developments that enable our PIC simulations. Then, in Section 3, we describe our simulation setup in detail. In Section 4, we present the results of our simulations, comparing and contrasting Klein–Nishina radiative reconnection to two control cases: one of non-radiative reconnection and one of classical radiative reconnection subject to efficient Thomson IC losses. Section 5 then provides a second results section. However, there, instead of delving into a detailed analysis of a few simulations with different radiative physics, we conduct parameter scans with all of the QED physics turned on, characterizing the pair yield of reconnection – a single number computed per simulation – as a function of its main controlling parameters. In Section 6, we survey observations of the four main application systems targeted by this work – FSRQs, BHXRBs, the M87* magnetosphere, and gamma-ray binaries – discussing connections to our simulation results. We conclude with a complete summary of our findings in Section 7. In the remaining part of this Introduction, we preview the three principal astrophysical results of this study.

The first concerns the correlated spectral and temporal signatures of Klein–Nishina reconnection (Section 4.5). As in the non-radiative and classical (Thomson IC) radiative regimes, Klein–Nishina reconnection powers efficient non-thermal particle acceleration (NTPA) and, hence, non-thermal radiative emission. While the time-averaged observable spectrum is similar to that of Thomson IC reconnection, the relationship between the shape of the output spectral energy density and its luminosity are very different. In Thomson IC reconnection, these two are tightly correlated, with a shallower observed spectrum coinciding with a higher luminosity (i.e. ‘harder-when-brighter’). In Klein–Nishina reconnection, this correlation is broken: the spectrum exhibits a constant shape irrespective of overall brightness. These results, potentially observable during gamma-ray flares, represent an important distinguishing property of Klein–

Nishina reconnection and illustrate the importance of collecting
temporally resolved observed spectra.

Our second main astrophysical result concerns the electron–
positron pair yield of Klein–Nishina reconnection (Section 5). Using
simulations, we derive an empirical formula for the Klein–Nishina
reconnection pair yield in terms of a single control parameter. While,
in the regime probed by this study, reconnection generally produces,
at most, order-unity new pairs per processed pair, our derived scaling
law, together with physical arguments for its extrapolation, point
to a potential regime where Klein–Nishina reconnection may be
a copious pair source. However, even in the case of order-unity
pair yield, Klein–Nishina reconnection can still convert an initially
electron-ion plasma into a strongly mixed electron-ion-positron
plasma. This reconnection regime is thus a potentially important
*in situ* antimatter source in astrophysics.

Our final main astrophysical result is a detailed survey of the four main object classes – FSRQs, BHXRBs, the M87* magnetosphere, and gamma-ray binaries – where Klein–Nishina reconnection may occur (Section 6). We find potentially strong observational connections to GeV observations of FSRQs, where observed spectral stability during flares is reminiscent of the anticipated spectral-temporal signatures of Klein–Nishina reconnection. We also sketch a roadmap of the theoretical and instrumental developments necessary to link Klein–Nishina reconnection modelling more rigorously to observations in the other systems. We comment on the potential effect of the Klein–Nishina reconnection pair yield on the global operation and observable aspects of each examined object type.

2 SIMULATION ARCHITECTURE

We here detail the QED capabilities that we added to the radiative electromagnetic PIC code ZELTRON (Cerutti et al. 2013; Cerutti & Werner 2019) to enable the simulations presented in this article. Readers wishing to skip these technical details altogether may proceed directly to Section 3; those interested in only a general description may prefer covering just Section 2.1 before moving on. Excellent additional references on QED methods in PIC simulations can be found in the literature documenting other PIC codes, including: TRISTAN V2 (e.g. Hakobyan et al. 2019, 2023a), OSIRIS (e.g. Fonseca et al. 2002; Del Gaudio et al. 2020), SMILEI (Derouillat et al. 2018), WARPX (Fedeli et al. 2022), and GRZELTRON (e.g. Levinson & Cerutti 2018; Crinquand et al. 2020). With the exception of the null collision method (see Section 2.2.3), the algorithms employed by us are the same as those used in GRZELTRON (cf. the supporting information of Crinquand et al. 2020). However, to our knowledge, many of the specifics – especially the optimization techniques of Section 2.2.3 – have not yet been described to the same level of detail as follows.

## 2.1 The QED PIC method

To provide some context for the modifications we have made to the ZELTRON code, we first review salient general features of electromagnetic PIC codes (though see, e.g., Hockney & Eastwood 1988; Birdsall & Langdon 1991 for a more thorough treatment), discussing how they may be extended to include QED effects. Fig. 1 provides a graphic summary of this discussion.

The electromagnetic PIC (hereafter, simply ‘PIC’) technique is a
computational method for simulating first-principles kinetic plasma
physics. PIC simulations are kinetic in that they self-consistently
describe the full phase space (position + momentum) plasma
distribution function (in contrast to, e.g., fluid plasma frameworks,
Hakobyan et al. 2023b); and gamma-ray binaries, where a plausible scenario involves a pulsar in tight orbit around a bright type O or Be star, which illuminates reconnection occurring near the pulsar in its magnetosphere and striped wind (i.e. before the pulsar wind shocks with the stellar wind from the companion; cf. Dubus 2006; Cerutti, Dubus & Henri 2008; Dubus 2013; Cerutti & Philippov 2017; Dubus et al. 2017; Philippov & Spitkovsky 2018; Cerutti, Philippov & Dubus 2020). The link between Klein–Nishina reconnection and each of these object classes is, in fact, a major result of this work (Section 6), as further discussed below.

While this study is primarily numerical, analytic, and semi-
analytic modelling are also vital for understanding Klein–Nishina
reconnection and QED reconnection more broadly. Such theoretical
approaches can make targeted, physically motivated arguments for
how results from non-QED reconnection may generalize to the
QED case (e.g. Beloborodov 2017; Mehlhaff et al. 2021; Chen et
al. 2023; Hakobyan et al. 2023b), even if they cannot treat all of
the kinetic physics at play from first principles. This furnishes a
useful interpretive framework for *ab initio* simulations. However,
the reverse is also true: phenomenological models, which sometimes
have the advantage of enhanced physical clarity, can themselves
be refined from the findings of simulations. This study illustrates
both directions of this paradigm. Throughout the text, we make
frequent reference to our earlier work, Mehlhaff et al. (2021), which
analytically considers the setup simulated here. As will be seen,
that study (besides laying much of the theoretical foundation for
this article) serves both as an interpretive lens for our simulations
and as a set of hypotheses that the numerical experiments can
check.

We structure this article as follows. In Section 2, we detail the
QED algorithmic developments that enable our PIC simulations.
Then, in Section 3, we describe our simulation setup in detail. In
Section 4, we present the results of our simulations, comparing
and contrasting Klein–Nishina radiative reconnection to two control
cases: one of non-radiative reconnection and one of classical radiative
reconnection subject to efficient Thomson IC losses. Section 5
then provides a second results section. However, there, instead of
delving into a detailed analysis of a few simulations with different
radiative physics, we conduct parameter scans with all of the QED
physics turned on, characterizing the pair yield of reconnection –
a single number computed per simulation – as a function of its
main controlling parameters. In Section 6, we survey observations
of the four main application systems targeted by this work – FSRQs,
BHXRBs, the M87* magnetosphere, and gamma-ray binaries –
discussing connections to our simulation results. We conclude with
a complete summary of our findings in Section 7. In the remaining
part of this Introduction, we preview the three principal astrophysical
results of this study.

The first concerns the correlated spectral and temporal signatures
of Klein–Nishina reconnection (Section 4.5). As in the non-radiative
and classical (Thomson IC) radiative regimes, Klein–Nishina recon-
nection powers efficient non-thermal particle acceleration (NTPA)
and, hence, non-thermal radiative emission. While the time-averaged
observable spectrum is similar to that of Thomson IC reconnection,
the relationship between the shape of the output spectral energy
density and its luminosity are very different. In Thomson IC
reconnection, these two are tightly correlated, with a shallower
observed spectrum coinciding with a higher luminosity (i.e. ‘harder-
when-brighter’). In Klein–Nishina reconnection, this correlation is
broken: the spectrum exhibits a constant shape irrespective of overall
brightness. These results, potentially observable during gamma-ray
flares, represent an important distinguishing property of Klein–

Nishina reconnection and illustrate the importance of collecting
temporally resolved observed spectra.

Our second main astrophysical result concerns the electron–
positron pair yield of Klein–Nishina reconnection (Section 5). Using
simulations, we derive an empirical formula for the Klein–Nishina
reconnection pair yield in terms of a single control parameter. While,
in the regime probed by this study, reconnection generally produces,
at most, order-unity new pairs per processed pair, our derived scaling
law, together with physical arguments for its extrapolation, point
to a potential regime where Klein–Nishina reconnection may be
a copious pair source. However, even in the case of order-unity
pair yield, Klein–Nishina reconnection can still convert an initially
electron-ion plasma into a strongly mixed electron-ion-positron
plasma. This reconnection regime is thus a potentially important
*in situ* antimatter source in astrophysics.

Our final main astrophysical result is a detailed survey of the four
main object classes – FSRQs, BHXRBs, the M87* magnetosphere,
and gamma-ray binaries – where Klein–Nishina reconnection may
occur (Section 6). We find potentially strong observational con-
nections to GeV observations of FSRQs, where observed spectral
stability during flares is reminiscent of the anticipated spectral-
temporal signatures of Klein–Nishina reconnection. We also sketch a
roadmap of the theoretical and instrumental developments necessary
to link Klein–Nishina reconnection modelling more rigorously to
observations in the other systems. We comment on the potential
effect of the Klein–Nishina reconnection pair yield on the global
operation and observable aspects of each examined object type.

## 2 SIMULATION ARCHITECTURE

We here detail the QED capabilities that we added to the radiative
electromagnetic PIC code ZELTRON (Cerutti et al. 2013; Cerutti &
Werner 2019) to enable the simulations presented in this article.
Readers wishing to skip these technical details altogether may
proceed directly to Section 3; those interested in only a general
description may prefer covering just Section 2.1 before moving on.
Excellent additional references on QED methods in PIC simulations
can be found in the literature documenting other PIC codes, includ-
ing: TRISTAN V2 (e.g. Hakobyan et al. 2019, 2023a), OSIRIS (e.g.
Fonseca et al. 2002; Del Gaudio et al. 2020), SMILEI (Derouillat et al.
2018), WARPX (Fedeli et al. 2022), and GRZELTRON (e.g. Levinson &
Cerutti 2018; Crinquand et al. 2020). With the exception of the null
collision method (see Section 2.2.3), the algorithms employed by
us are the same as those used in GRZELTRON (cf. the supporting
information of Crinquand et al. 2020). However, to our knowledge,
many of the specifics – especially the optimization techniques of
Section 2.2.3 – have not yet been described to the same level of
detail as follows.

## 2.1 The QED PIC method

To provide some context for the modifications we have made
to the ZELTRON code, we first review salient general features of
electromagnetic PIC codes (though see, e.g., Hockney & Eastwood
1988; Birdsall & Langdon 1991 for a more thorough treatment),
discussing how they may be extended to include QED effects. Fig. 1
provides a graphic summary of this discussion.

The electromagnetic PIC (hereafter, simply ‘PIC’) technique is a
computational method for simulating first-principles kinetic plasma
physics. PIC simulations are kinetic in that they self-consistently
describe the full phase space (position + momentum) plasma
distribution function (in contrast to, e.g., fluid plasma frameworks,
<figure><img src="image_1.png" /><figcaption>**Figure 1.** Standard PIC loop augmented with steps (dashed outline) to model QED effects. At each PIC time-step, the electromagnetic field – plus an optional continuous radiative drag force to model non-QED cooling – are used to update the particles' positions and momenta according to the Lorentz force law (top left). Photons are also ballistically propagated at the speed of light in this step. Then (top right), QED cross-sections are used to compute Monte Carlo photon emission and pair production. Each time a Monte Carlo photon is emitted, its momentum is self-consistently subtracted from that of its radiating particle. Freshly emitted photons and produced pairs are added to the simulation. In the next step (bottom right), the particle positions and velocities are used to calculate the electromagnetic charge and current densities. These are used in the final step (bottom left) to update the electromagnetic fields via the time-dependent Maxwell's equations.</figcaption></figure>

which track bulk quantities – such as spatial density and local mean velocity – in real space only). The PIC technique is, furthermore, a first-principles method because it evolves physical equations (the Maxwell–Vlasov system) requiring minimal approximations. Owing to these properties, PIC simulations can probe detailed microscopic plasma physical effects, while furnishing vital, self-consistent astro-physical observables such as light curves and spectra.

The PIC method represents the simulated plasma as a large number of discrete charged particles coupled to electric and magnetic fields, **E** and **B**, respectively. The **E** and **B** fields are tracked on a spatial grid, while the particles' positions, **x**, and momenta, **p**, can vary continuously. At each time-step, the gridded field values are interpolated to the positions of the particles, allowing their momenta to be evolved via the Lorentz force law,

$$ \frac{d\mathbf{p}}{dt} = q \left( \mathbf{E} + \frac{\mathbf{v}}{c} \times \mathbf{B} \right) \quad (1) $$

(Fig. 1, red panel; **v** is the 3-velocity for a particle of momentum **p**). Then, to capture the response of the fields to the particles, the particles' positions and momenta are used to calculate the bulk charge, **ρ**, and current, **J**, densities on the computational field grid (Fig. 1, green panel). This allows **E** and **B** to be advanced via the time-dependent Maxwell's equations,

$$ \begin{aligned} \partial_t \mathbf{B} &= -c \nabla \times \mathbf{E} \quad \text{and} \\ \partial_t \mathbf{E} &= c \nabla \times \mathbf{B} - 4\pi \mathbf{J} \end{aligned} \quad (2) $$

(Fig. 1, blue panel). It is worth noting that equations (1) and (2) imply a collisionless PIC method, since the particles only interact collectively through the long-range fields **E** and **B** rather than individually through short-range two- or few-body encounters.

By itself, the procedure described so far does not model the high-energy radiation (with photon frequencies unresolved by the simulation time-step) that is often important in the environments of relativistic compact objects, as in this work. The simplest way that

such radiative effects can be incorporated is by adding a cooling radiative drag term, **f**<sub>rad</sub>, to equation (1) such that it becomes (Fig. 1, red panel)

$$ \frac{d\mathbf{p}}{dt} = q \left( \mathbf{E} + \frac{\mathbf{v}}{c} \times \mathbf{B} \right) + \mathbf{f}_{\text{rad}}. \quad (3) $$

Treating radiative cooling as a continuous drag force is suitable when the radiating particles do not lose a significant fraction of their energy to any single photon emission event. However, when particles begin to emit photons at energies that rival their own, radiation becomes inherently discrete, and a more general approach is necessary. As an additional concern, when the emitted photons are above the threshold energy for one or more pair-production processes, their propagation and absorption must be handled self-consistently.

The needed additional QED operations can be mostly consolidated into one extra step beyond the standard PIC loop (Fig. 1, orange panel). Here, the particle positions and velocities are used – perhaps in conjunction with the electromagnetic fields – to evaluate QED cross-sections for the photon emission and pair-production processes of interest, yielding probabilities for these events to occur. A subset of the possible events are triggered by comparing their probabilities with randomly drawn numbers – a Monte Carlo procedure – and the resulting photons and particles are added to the simulation. Besides this additional QED Monte Carlo step (Fig. 1, orange panel), the positions of photons also need to be evolved. This is typically done alongside the particles (Fig. 1, red panel) and is much easier by comparison because photons follow straight lines (in flat space–time, as in this work). We term the PIC method, expanded to include these extra operations as diagrammed in Fig. 1, the QED PIC method.

To enable this work, we have generalized the Cartesian version of the PIC code ZELTRON to include QED physics as sketched here. The emission and pair-production processes – relevant to the orange panel of Fig. 1 – are: IC scattering (including the high-energy Klein–Nishina limit) of a soft seed photon background by ultrarelativistic particles; and pair production when the high-energy Comptonized photons are absorbed by the soft background. In the following discussion, we provide a sketch of the algorithm used to model these processes.

## 2.2 QED effects implemented in ZELTRON

### 2.2.1 Monte Carlo inverse Compton emission procedure

We begin by describing algorithmic details of the IC emission mechanism. In this section, we present a conceptually simple but computationally expensive implementation. In Section 2.2.3, we discuss the modifications we made to the simpler procedure of this section for the sake of optimization.

Our simulations are immersed in a homogeneous, static, isotropic, and monochromatic photon radiation bath of energy density (per unit energy interval)

$$ u(\epsilon) = U_{\text{ph}}\delta(\epsilon - \epsilon_{\text{ph}}). \quad (4) $$

These background photons are not tracked by the simulation; their energy density is prescribed by hand. Electrons and positrons, through the IC process, can upscatter these background photons. If an upscattered photon attains final energy above the (gamma-ray) threshold for pair production with the background radiation bath, then it is promoted to a tracked photon whose position is evolved explicitly by the simulation.
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Standard PIC loop augmented with steps (dashed outline) to model QED effects. At each PIC time-step, the electromagnetic field – plus an optional continuous radiative drag force to model non-QED cooling – are used to update the particles' positions and momenta according to the Lorentz force law (top left). Photons are also ballistically propagated at the speed of light in this step. Then (top right), QED cross-sections are used to compute Monte Carlo photon emission and pair production. Each time a Monte Carlo photon is emitted, its momentum is self-consistently subtracted from that of its radiating particle. Freshly emitted photons and produced pairs are added to the simulation. In the next step (bottom right), the particle positions and velocities are used to calculate the electromagnetic charge and current densities. These are used in the final step (bottom left) to update the electromagnetic fields via the time-dependent Maxwell's equations.</figcaption></figure>

which track bulk quantities – such as spatial density and local mean velocity – in real space only). The PIC technique is, furthermore, a first-principles method because it evolves physical equations (the Maxwell–Vlasov system) requiring minimal approximations. Owing to these properties, PIC simulations can probe detailed microscopic plasma physical effects, while furnishing vital, self-consistent astro-physical observables such as light curves and spectra.

The PIC method represents the simulated plasma as a large number of discrete charged particles coupled to electric and magnetic fields, **E** and **B**, respectively. The **E** and **B** fields are tracked on a spatial grid, while the particles' positions, **x**, and momenta, **p**, can vary continuously. At each time-step, the gridded field values are interpolated to the positions of the particles, allowing their momenta to be evolved via the Lorentz force law,

$$ \frac{d\mathbf{p}}{dt} = q \left( \mathbf{E} + \frac{\mathbf{v}}{c} \times \mathbf{B} \right) \quad (1) $$

(Fig. 1, red panel; **v** is the 3-velocity for a particle of momentum **p**). Then, to capture the response of the fields to the particles, the particles' positions and momenta are used to calculate the bulk charge, **ρ**, and current, **J**, densities on the computational field grid (Fig. 1, green panel). This allows **E** and **B** to be advanced via the time-dependent Maxwell's equations,

$$ \begin{aligned} \partial_t \mathbf{B} &= -c \nabla \times \mathbf{E} \quad \text{and} \\ \partial_t \mathbf{E} &= c \nabla \times \mathbf{B} - 4\pi \mathbf{J} \end{aligned} \quad (2) $$

(Fig. 1, blue panel). It is worth noting that equations (1) and (2) imply a collisionless PIC method, since the particles only interact collectively through the long-range fields **E** and **B** rather than individually through short-range two- or few-body encounters.

By itself, the procedure described so far does not model the high-energy radiation (with photon frequencies unresolved by the simulation time-step) that is often important in the environments of relativistic compact objects, as in this work. The simplest way that

such radiative effects can be incorporated is by adding a cooling radiative drag term, **f**<sub>rad</sub>, to equation (1) such that it becomes (Fig. 1, red panel)

$$ \frac{d\mathbf{p}}{dt} = q \left( \mathbf{E} + \frac{\mathbf{v}}{c} \times \mathbf{B} \right) + \mathbf{f}_{\text{rad}}. \quad (3) $$

Treating radiative cooling as a continuous drag force is suitable when the radiating particles do not lose a significant fraction of their energy to any single photon emission event. However, when particles begin to emit photons at energies that rival their own, radiation becomes inherently discrete, and a more general approach is necessary. As an additional concern, when the emitted photons are above the threshold energy for one or more pair-production processes, their propagation and absorption must be handled self-consistently.

The needed additional QED operations can be mostly consolidated into one extra step beyond the standard PIC loop (Fig. 1, orange panel). Here, the particle positions and velocities are used – perhaps in conjunction with the electromagnetic fields – to evaluate QED cross-sections for the photon emission and pair-production processes of interest, yielding probabilities for these events to occur. A subset of the possible events are triggered by comparing their probabilities with randomly drawn numbers – a Monte Carlo procedure – and the resulting photons and particles are added to the simulation. Besides this additional QED Monte Carlo step (Fig. 1, orange panel), the positions of photons also need to be evolved. This is typically done alongside the particles (Fig. 1, red panel) and is much easier by comparison because photons follow straight lines (in flat space–time, as in this work). We term the PIC method, expanded to include these extra operations as diagrammed in Fig. 1, the QED PIC method.

To enable this work, we have generalized the Cartesian version of the PIC code ZELTRON to include QED physics as sketched here. The emission and pair-production processes – relevant to the orange panel of Fig. 1 – are: IC scattering (including the high-energy Klein–Nishina limit) of a soft seed photon background by ultrarelativistic particles; and pair production when the high-energy Comptonized photons are absorbed by the soft background. In the following discussion, we provide a sketch of the algorithm used to model these processes.

## 2.2 QED effects implemented in ZELTRON

### 2.2.1 Monte Carlo inverse Compton emission procedure

We begin by describing algorithmic details of the IC emission mechanism. In this section, we present a conceptually simple but computationally expensive implementation. In Section 2.2.3, we discuss the modifications we made to the simpler procedure of this section for the sake of optimization.

Our simulations are immersed in a homogeneous, static, isotropic, and monochromatic photon radiation bath of energy density (per unit energy interval)

$$ u(\epsilon) = U_{\text{ph}}\delta(\epsilon - \epsilon_{\text{ph}}). \quad (4) $$

These background photons are not tracked by the simulation; their energy density is prescribed by hand. Electrons and positrons, through the IC process, can upscatter these background photons. If an upscattered photon attains final energy above the (gamma-ray) threshold for pair production with the background radiation bath, then it is promoted to a tracked photon whose position is evolved explicitly by the simulation.
For a particle with ultrarelativistic Lorentz factor $\gamma \gg 1$ traversing the radiation field (4), the number of photons scattered per unit time is given as

$$ \frac{dN}{dt} = c\sigma_T n_{ph} g_{KN} (\gamma/\gamma_{KN}), \qquad (5) $$

where $\sigma_T$ is the Thomson cross-section, $n_{ph} \equiv U_{ph}/\epsilon_{ph}$ is the background photon number density,

$$ \gamma_{KN} \equiv m_e c^2 / 4 \epsilon_{ph} \qquad (6) $$

is the critical Lorentz factor above which IC scattering transitions to the Klein–Nishina regime, and $g_{KN}(q) \le 1$ is the dimensionless function (calculated in appendix A of Mehlhaff et al. 2021)

$$ g_{KN}(q) = \frac{3}{2q^2} \left[ \left( q + 9 + \frac{8}{q} \right) \ln(1 + q) - \frac{1}{1 + q} \left( \frac{q^2}{2} + 9q + 8 \right) + 4 \mathrm{Li}_2(-q) \right]. \qquad (7) $$

Here, $\mathrm{Li}_2(q)$ is the dilogarithm. As needed to recover the Thomson regime, $g_{KN}(\gamma/\gamma_{KN})$ tends to unity in the low-energy limit, $\gamma \ll \gamma_{KN}$. Over a simulation time-step $\Delta t$, equation (5) yields a probability

$$ p_{\text{emit}} = \Delta t dN/dt \qquad (8) $$

for the particle to scatter a photon. We operate in the regime $p_0 \equiv c\sigma_T n_{ph} \Delta t \ll 1$, implying $p_{\text{emit}} = p_0 g_{KN}(q) \le p_0 \ll 1$.

Computationally, the Monte Carlo photon emission process can be accomplished, for each particle at each time-step, by evaluating the probability (8) and drawing a random number to determine whether the emission event occurs. Then, for the subset of particles that actually emit photons, the scattered photon energy $\epsilon_{scat}$ must be assigned. This demands drawing a second random number $\mathcal{R}$ and inverting the probability distribution over $\epsilon_{scat}$: that is, finding $\epsilon_{scat}$ such that

$$ \mathcal{R} = \frac{\int_{0}^{r(\epsilon_{\text{scat}})} K(r', q) dr'}{\int_{0}^{1} K(r', q) dr'} . \qquad (9) $$

Here, $q \equiv \gamma/\gamma_{KN}$ as above, $r(\epsilon_{scat})$ is a proxy for the photon energy defined by

$$ r(\epsilon_{\text{scat}}) \equiv \frac{\epsilon_{\text{scat}} / \gamma m_e c^2}{q (1 - \epsilon_{\text{scat}} / \gamma m_e c^2)}, \qquad (10) $$

and $K(r, q)$ is the single-particle Klein–Nishina scattering kernel as reported by Jones (1968) and Blumenthal & Gould (1970):

$$ K(r, q) = \frac{3}{(1 + qr)^2} \left[ 2r \ln r + (1 + 2r)(1 - r) + \frac{1}{2} \frac{(qr)^2}{1 + qr} (1 - r) \right]. \qquad (11) $$

That is, the number of photons scattered per unit time by a particle with $\gamma = q\gamma_{KN}$ to final proxy photon energy between $r$ and $r + dr$ is

$$ \frac{dN}{dr} dr = c\sigma_T n_{ph} K(r, q) dr. \qquad (12) $$

With the photon energy $\epsilon_{scat}$ known, the momentum of the scattered photon has magnitude $\epsilon_{scat}/c$ and points, for the ultrarelativistic $\gamma \gg 1$ approximation relevant to our simulations, along the direction of the radiating particle's motion. This momentum is subtracted from that of the radiating particle. Furthermore, if $\epsilon_{scat}$ is above pair-production threshold with the background photons – that is, if $\epsilon_{scat}\epsilon_{ph}/(m_e c^2)^2 \ge 1$ – the scattered photon is promoted to a tracked simulation particle, allowing its later potential absorption to be self-consistently calculated.

Though the Monte Carlo emission scheme described so far relies on the ultrarelativistic approximation $\gamma \gg 1$, some of our

simulations contain trans- or even subrelativistic particles. Thus, at lower energies, we need to stitch to a cooling procedure that is non-relativistically valid. To that end, we restrict ourselves to the main regime of astrophysical interest where $\gamma_{KN} = m_e c^2 / 4 \epsilon_{ph} \gg 1$ (i.e. $\epsilon_{ph} \ll 130$ keV). This means that there is a broad range of particle energies $\gamma$ that are ultrarelativistic ($\gamma \gg 1$) but still well below the threshold where quantum Klein–Nishina effects kick in ($\gamma \ll \gamma_{KN}$): a shared applicability range where either the Monte Carlo emission scheme, which demands $\gamma \gg 1$, or a continuous Thomson radiative drag force, which needs $\gamma \ll \gamma_{KN}$, could be used. Within this range, we select a threshold ‘stitching’ particle energy $\gamma_{Th-KN}$. Above $\gamma_{Th-KN}$, we employ the IC Monte Carlo emission procedure. Below $\gamma_{Th-KN}$, we switch to the continuous Thomson radiative drag force (e.g. Blumenthal & Gould 1970; Rybicki & Lightman 1979; Pozdnyakov, Sobol & Syunyaev 1983; for works in the context of relativistic reconnection, cf. Uzdensky 2016; Werner et al. 2019; Mehlhaff et al. 2020, 2021; Sironi & Beloborodov 2020; Sridhar et al. 2021, 2023),

$$ f_{\text{rad}} = -(4/3)\sigma_T \gamma^2 U_{ph}\beta, \qquad (13) $$

that enters (via the method of Tamburini et al. 2010) into the particle push through equation (3) and is, importantly, non-relativistically correct. Regarding the choice of $\gamma_{Th-KN}$, we find that artefacts of the stitching generally disappear when $\gamma_{Th-KN} \le 0.1\gamma_{KN}$ – well inside of the Thomson regime. At the same time, we find that $\gamma_{Th-KN}$ should be at least of order a few, limiting us to $\gamma_{KN}$ values that are above 20 or so.

### 2.2.2 Monte Carlo pair-production procedure

We next describe how pair production is processed in the code. For a collision between a gamma-ray and a background photon with angle $\theta$ between their velocity vectors and respective energies $\epsilon_{hard}$ and $\epsilon_{ph}$, the centre of mass energy is $s = s_0(1 - \cos\theta)/2 \le s_0$, where $s_0 \equiv \epsilon_{hard}\epsilon_{ph}/m_e^2 c^4$ must exceed unity for pair production to be possible. Integrating over all possible collision angles $\theta$ such that $s > 1$ for a given $s_0$ yields an overall cross-section presented to a propagating gamma-ray by the background radiation (Gould & Schréder 1967) of<fnref n="1" />

$$ \sigma_{\gamma\gamma}(s_0) = \frac{3}{8} \frac{\sigma_T}{s_0^2} \left[ \frac{1 + \beta_0^2}{1 - \beta_0^2} \ln W_0 - \beta_0^2 \ln W_0 - \ln^2 W_0 - \frac{4\beta_0}{1 - \beta_0^2} + 2\beta_0 + 4 \ln W_0 \ln (1 + W_0) + 4 \left( \frac{\pi^2}{12} + \mathrm{Li}_2(-W_0) \right) \right], \qquad (14) $$

where $\beta_0$ and $W_0$ are both $s_0$-dependent, reading, respectively,

$$ \beta_0^2(s_0) \equiv 1 - \frac{1}{s_0} \qquad (15) $$

and

$$ W_0(s_0) \equiv \frac{1 + \beta_0(s_0)}{1 - \beta_0(s_0)}. \qquad (16) $$

Thus, the pair-production probability accumulated by a gamma-ray of energy $\epsilon_{hard}$ in one simulation time-step $\Delta t$ is

$$ p_{\text{abs}} = c n_{\text{ph}} \sigma_{\gamma\gamma}(s_0) \Delta t. \qquad (17) $$

<footnote n="1">¹The parenthesized term in equation (14), $4(\pi^2/12 + \mathrm{Li}_2(-W_0))$, is equal to the last term, $-4L(W_0)$, in equation (10) of Gould & Schréder (1967), but is corrected for the missing factor of 4 in that work.</footnote>
For a particle with ultrarelativistic Lorentz factor $\gamma \gg 1$ traversing the radiation field (4), the number of photons scattered per unit time is given as

$$ \frac{dN}{dt} = c\sigma_T n_{ph} g_{KN} (\gamma/\gamma_{KN}), \qquad (5) $$

where $\sigma_T$ is the Thomson cross-section, $n_{ph} \equiv U_{ph}/\epsilon_{ph}$ is the background photon number density,

$$ \gamma_{KN} \equiv m_e c^2 / 4 \epsilon_{ph} \qquad (6) $$

is the critical Lorentz factor above which IC scattering transitions to the Klein–Nishina regime, and $g_{KN}(q) \le 1$ is the dimensionless function (calculated in appendix A of Mehlhaff et al. 2021)

$$ g_{KN}(q) = \frac{3}{2q^2} \left[ \left( q + 9 + \frac{8}{q} \right) \ln(1 + q) - \frac{1}{1 + q} \left( \frac{q^2}{2} + 9q + 8 \right) + 4 \mathrm{Li}_2(-q) \right]. \qquad (7) $$

Here, $\mathrm{Li}_2(q)$ is the dilogarithm. As needed to recover the Thomson regime, $g_{KN}(\gamma/\gamma_{KN})$ tends to unity in the low-energy limit, $\gamma \ll \gamma_{KN}$. Over a simulation time-step $\Delta t$, equation (5) yields a probability

$$ p_{\text{emit}} = \Delta t dN/dt \qquad (8) $$

for the particle to scatter a photon. We operate in the regime $p_0 \equiv c\sigma_T n_{ph} \Delta t \ll 1$, implying $p_{\text{emit}} = p_0 g_{KN}(q) \le p_0 \ll 1$.

Computationally, the Monte Carlo photon emission process can be accomplished, for each particle at each time-step, by evaluating the probability (8) and drawing a random number to determine whether the emission event occurs. Then, for the subset of particles that actually emit photons, the scattered photon energy $\epsilon_{scat}$ must be assigned. This demands drawing a second random number $\mathcal{R}$ and inverting the probability distribution over $\epsilon_{scat}$: that is, finding $\epsilon_{scat}$ such that

$$ \mathcal{R} = \frac{\int_{0}^{r(\epsilon_{\text{scat}})} K(r', q) dr'}{\int_{0}^{1} K(r', q) dr'}. \qquad (9) $$

Here, $q \equiv \gamma/\gamma_{KN}$ as above, $r(\epsilon_{scat})$ is a proxy for the photon energy defined by

$$ r(\epsilon_{\text{scat}}) \equiv \frac{\epsilon_{\text{scat}} / \gamma m_e c^2}{q (1 - \epsilon_{\text{scat}} / \gamma m_e c^2)}, \qquad (10) $$

and $K(r, q)$ is the single-particle Klein–Nishina scattering kernel as reported by Jones (1968) and Blumenthal & Gould (1970):

$$ K(r, q) = \frac{3}{(1 + qr)^2} \left[ 2r \ln r + (1 + 2r)(1 - r) + \frac{1}{2} \frac{(qr)^2}{1 + qr} (1 - r) \right]. \qquad (11) $$

That is, the number of photons scattered per unit time by a particle with $\gamma = q\gamma_{KN}$ to final proxy photon energy between $r$ and $r + dr$ is

$$ \frac{dN}{dr} dr = c\sigma_T n_{ph} K(r, q) dr. \qquad (12) $$

With the photon energy $\epsilon_{scat}$ known, the momentum of the scattered photon has magnitude $\epsilon_{scat}/c$ and points, for the ultrarelativistic $\gamma \gg 1$ approximation relevant to our simulations, along the direction of the radiating particle's motion. This momentum is subtracted from that of the radiating particle. Furthermore, if $\epsilon_{scat}$ is above pair-production threshold with the background photons – that is, if $\epsilon_{scat}\epsilon_{ph}/(m_e c^2)^2 \ge 1$ – the scattered photon is promoted to a tracked simulation particle, allowing its later potential absorption to be self-consistently calculated.

Though the Monte Carlo emission scheme described so far relies on the ultrarelativistic approximation $\gamma \gg 1$, some of our

simulations contain trans- or even subrelativistic particles. Thus, at lower energies, we need to stitch to a cooling procedure that is non-relativistically valid. To that end, we restrict ourselves to the main regime of astrophysical interest where $\gamma_{KN} = m_e c^2/4\epsilon_{ph} \gg 1$ (i.e. $\epsilon_{ph} \ll 130$ keV). This means that there is a broad range of particle energies $\gamma$ that are ultrarelativistic ($\gamma \gg 1$) but still well below the threshold where quantum Klein–Nishina effects kick in ($\gamma \ll \gamma_{KN}$): a shared applicability range where either the Monte Carlo emission scheme, which demands $\gamma \gg 1$, or a continuous Thomson radiative drag force, which needs $\gamma \ll \gamma_{KN}$, could be used. Within this range, we select a threshold ‘stitching’ particle energy $\gamma_{Th-KN}$. Above $\gamma_{Th-KN}$, we employ the IC Monte Carlo emission procedure. Below $\gamma_{Th-KN}$, we switch to the continuous Thomson radiative drag force (e.g. Blumenthal & Gould 1970; Rybicki & Lightman 1979; Pozdnyakov, Sobol & Syunyaev 1983; for works in the context of relativistic reconnection, cf. Uzdensky 2016; Werner et al. 2019; Mehlhaff et al. 2020, 2021; Sironi & Beloborodov 2020; Sridhar et al. 2021, 2023),

$$ f_{\text{rad}} = -(4/3)\sigma_T \gamma^2 U_{ph}\beta, \qquad (13) $$

that enters (via the method of Tamburini et al. 2010) into the particle push through equation (3) and is, importantly, non-relativistically correct. Regarding the choice of $\gamma_{Th-KN}$, we find that artefacts of the stitching generally disappear when $\gamma_{Th-KN} \le 0.1\gamma_{KN}$ – well inside of the Thomson regime. At the same time, we find that $\gamma_{Th-KN}$ should be at least of order a few, limiting us to $\gamma_{KN}$ values that are above 20 or so.

### 2.2.2 Monte Carlo pair-production procedure

We next describe how pair production is processed in the code. For a collision between a gamma-ray and a background photon with angle $\theta$ between their velocity vectors and respective energies $\epsilon_{hard}$ and $\epsilon_{ph}$, the centre of mass energy is $s = s_0(1 - \cos\theta)/2 \le s_0$, where $s_0 \equiv \epsilon_{hard}\epsilon_{ph}/m_e^2 c^4$ must exceed unity for pair production to be possible. Integrating over all possible collision angles $\theta$ such that $s > 1$ for a given $s_0$ yields an overall cross-section presented to a propagating gamma-ray by the background radiation (Gould & Schréder 1967) of<fnref n="1" />

$$ \sigma_{\gamma\gamma}(s_0) = \frac{3}{8} \frac{\sigma_T}{s_0^2} \left[ \frac{1 + \beta_0^2}{1 - \beta_0^2} \ln W_0 - \beta_0^2 \ln W_0 - \ln^2 W_0 - \frac{4\beta_0}{1 - \beta_0^2} + 2\beta_0 + 4 \ln W_0 \ln (1 + W_0) + 4 \left( \frac{\pi^2}{12} + \mathrm{Li}_2(-W_0) \right) \right], \qquad (14) $$

where $\beta_0$ and $W_0$ are both $s_0$-dependent, reading, respectively,

$$ \beta_0^2(s_0) \equiv 1 - \frac{1}{s_0} \qquad (15) $$

and

$$ W_0(s_0) \equiv \frac{1 + \beta_0(s_0)}{1 - \beta_0(s_0)}. \qquad (16) $$

Thus, the pair-production probability accumulated by a gamma-ray of energy $\epsilon_{hard}$ in one simulation time-step $\Delta t$ is

$$ p_{\text{abs}} = c n_{\text{ph}} \sigma_{\gamma\gamma}(s_0) \Delta t. \qquad (17) $$

<footnote n="1">¹The parenthesized term in equation (14), $4(\pi^2/12 + \mathrm{Li}_2(-W_0))$, is equal to the last term, $-4L(W_0)$, in equation (10) of Gould & Schréder (1967), but is corrected for the missing factor of 4 in that work.</footnote>
Because, max ($\sigma_{\gamma\gamma}$) $\simeq \sigma_T/5$, our choice to operate in the regime $p_0 \equiv c\sigma_T n_{ph} \Delta t \ll 1$ implies $p_{abs} \le p_0/5 \ll 1$.

Just as equation (8) determines which potential IC scattering events occur, equation (17) determines which photons produce new pairs. Ordinarily, like the additional steps necessary to determine the final scattered photon energy in the case of IC emission – for example, equation (9) – one would also need to proceed beyond equation (17) to fix the energy of the newborn electron and positron. First, the angle $\theta$, which is integrated out when interested only in the total cross-section (14) and corresponding absorption probability (17), needs to actually be sampled to determine the centre-of-mass energy $s$. Then, one must also sample the angle of one of the newborn particles’ momenta with respect to that of the collision axis in the centre-of-mass frame. We have indeed implemented both steps in ZELTRON, but they are unnecessary when $\epsilon_{hard} \gg m_e c^2 \gg \epsilon_{ph}$ (i.e. when $4\gamma_{KN} \gg 1$), which is all that concerns us in this study. Then, the collision energy budget is entirely dominated by the incoming gamma-ray, and the electron and positron each simply inherit half of the absorbed gamma-ray’s momentum.

### 2.2.3 Optimization methods for inverse Compton emission

The Monte Carlo implementation of IC emission sketched in Section 2.2.1 suffers from two performance bottlenecks. First, the loop computing the IC cross-section, equations (7) and (8), for every simulation particle is costly. Second, within this loop, the inversion of the cumulative distribution function in equation (9) is also expensive (though only necessary for the subset of particles that actually scatters photons). We mitigate these issues using two techniques, which we discuss in turn below. We note that, while similar techniques could be applied to the pair-production procedure, our simulations are optically thick to pair production, and the steady-state number of tracked photons is much smaller than the (always growing) number of particles. Thus, the QED physics involving photons, for us, can tolerate a less rigorously optimized implementation.

First, we speed up the assignment of the scattered photon energy – done in the context of Section 2.2.1 through equation (9). To do this, we consider the IC emission in the rest frame of the scattering particle. As we show below, this requires a larger number of random number draws per particle, but enables the use of simpler expressions. These, in turn, furnish an approximation to the scattering cross-section that effectively replaces equation (9), relieving its performance bottleneck. Our handling of the problem this way, presented in detail below, follows closely Levinson & Cerutti (2018) and Crinquand et al. (2020).

Before transforming to the rest frame of the potentially scattering particle, we isolate to an interaction with a single photon rather than with the entire isotropic radiation bath. This is accomplished by first drawing a random photon angle $\psi$ with respect to the direction of the particle’s velocity $c\beta$ according to the probability density function

$$p(\psi) \propto 1 - \beta \cos \psi, \qquad (18)$$

which accounts for the relative lab-frame rate of encounters of the particle with photons incident from different directions. The selected photon is then boosted to the particle’s rest frame, where it has energy

$$\epsilon_0' = \gamma \epsilon_{ph} (1 - \beta \cos \psi). \qquad (19)$$

Particle rest-frame quantities are primed in our convention.

We now evaluate the spectral (per unit final photon energy) and total (integrated over final photon energies) scattering rate for this interaction in the particle’s rest frame. We discuss after this derivation

how these two quantities are used by the code. In the primed frame, the scattering cross-section is that of ordinary Compton scattering (e.g. Blumenthal & Gould 1970):

$$\frac{d\sigma}{d\epsilon'_{\text{scat}} d\Omega'} = \frac{3\sigma_{\text{T}}}{16\pi} \left(\frac{\epsilon'_{\text{scat}}}{\epsilon'_{0}}\right)^2 \left(\frac{\epsilon'_{0}}{\epsilon'_{\text{scat}}} + \frac{\epsilon'_{\text{scat}}}{\epsilon'_{0}} - \sin^2 \Psi'\right) \times \delta \left[ \epsilon'_{\text{scat}} - \frac{\epsilon'_{0}}{1 + (\epsilon'_{0}/m_e c^2)(1 - \cos \Psi')} \right], \qquad (20)$$

where $\Psi'$ is the angle through which the photon is scattered in the particle frame. The total scattering rate into final photon energies between $\epsilon'_{scat}$ and $\epsilon'_{scat} + d\epsilon'_{scat}$ is then

$$\frac{dN}{dt d\epsilon'_{\text{scat}}} = \frac{1}{\gamma} \frac{dN}{dt' d\epsilon'_{\text{scat}}} = cn_{\text{ph}} \int d\Omega' \frac{d\sigma}{d\epsilon'_{\text{scat}} d\Omega'}, \qquad (21)$$

where we used $n'_{ph} = \gamma n_{ph}$. Defining $x \equiv \epsilon'_{scat}/\epsilon'_{0}$ and $y \equiv m_e c^2/\epsilon'_{0}$, the delta function only activates for $1/(1 + 2/y) \le x \le 1$, in which case the integral evaluates to

$$\frac{dN}{dt dx} = \epsilon'_{0} \frac{dN}{dt d\epsilon'_{\text{scat}}} = \frac{3\sigma_{\text{T}}}{8} \times y \left[ x + \frac{1}{x} + 2y \left(1 - \frac{1}{x}\right) + y^2 \left(\frac{1}{x^2} - \frac{2}{x} + 1\right) \right]. \qquad (22)$$

The total scattering rate is then

$$\frac{dN}{dt} \equiv \int dx \frac{dN}{dt dx} = \frac{3\sigma_{\text{T}}cn_{\text{ph}}}{8} \left[ F(1, y) - F\left(\frac{1}{1+2/y}, y\right) \right], \qquad (23)$$

where

$$F(x, y) \equiv y \left[ -\frac{y^2}{x} + xy(2+y) + \frac{x^2}{2} + (1-2y-2y^2) \ln x \right]. \qquad (24)$$

These results are used by the code as follows. First, equation (23) is evaluated and multiplied by $\Delta t$ to determine the probability $p_{emit} = \Delta t \, dN/dt$ that the scattering occurs during the simulation time-step. This step replaces the evaluation of equation (5). On average, these two procedures are completely equivalent, but in this second method, we have traded the analytic integral over incident photon directions [which yields equation (5)] for a random Monte Carlo sampling over these directions.

Then, for particles that scatter photons, the final photon energy is obtained by drawing a random number $R$ and inverting the cumulative distribution function: that is, finding $x$ such that

$$\mathcal{R} = \frac{3\sigma_{\text{T}}cn_{\text{ph}} [F(x, y) - F(1/(1 + 2/y), y)]}{8 \frac{dN}{dt}} = \frac{F(x, y) - F(1/(1 + 2/y), y)}{F(1, y) - F(1/(1 + 2/y), y)}. \qquad (25)$$

This replaces the inversion of the equivalent lab-frame cumulative distribution function (9) and, as we now show, constitutes the main advantage of this method. The issue with equation (9) is that its solution cannot be expressed analytically. Instead, the cumulative distribution function (or its inverse) must be stored as a table and consulted for each emitting particle – a costly procedure. However, the function $F(x, y)$ can be approximated, when $x \ll 1$, by only its $\ln x$ term. This approximation breaks down as $x$ approaches unity but, as it turns out, does not very strongly disturb the cumulative distribution function on the right-hand side of equation (25) even when used across all $x$. If one then approximates

$$F(x, y) \approx y (1 - 2y - 2y^2) \ln x, \qquad (26)$$
Because, max ($\sigma_{\gamma\gamma}$) $\simeq \sigma_T/5$, our choice to operate in the regime $p_0 \equiv c\sigma_T n_{ph} \Delta t \ll 1$ implies $p_{abs} \le p_0/5 \ll 1$.

Just as equation (8) determines which potential IC scattering events occur, equation (17) determines which photons produce new pairs. Ordinarily, like the additional steps necessary to determine the final scattered photon energy in the case of IC emission – for example, equation (9) – one would also need to proceed beyond equation (17) to fix the energy of the newborn electron and positron. First, the angle $\theta$, which is integrated out when interested only in the total cross-section (14) and corresponding absorption probability (17), needs to actually be sampled to determine the centre-of-mass energy $s$. Then, one must also sample the angle of one of the newborn particles’ momenta with respect to that of the collision axis in the centre-of-mass frame. We have indeed implemented both steps in ZELTRON, but they are unnecessary when $\epsilon_{hard} \gg m_e c^2 \gg \epsilon_{ph}$ (i.e. when $4\gamma_{KN} \gg 1$), which is all that concerns us in this study. Then, the collision energy budget is entirely dominated by the incoming gamma-ray, and the electron and positron each simply inherit half of the absorbed gamma-ray’s momentum.

### 2.2.3 Optimization methods for inverse Compton emission

The Monte Carlo implementation of IC emission sketched in Section 2.2.1 suffers from two performance bottlenecks. First, the loop computing the IC cross-section, equations (7) and (8), for every simulation particle is costly. Second, within this loop, the inversion of the cumulative distribution function in equation (9) is also expensive (though only necessary for the subset of particles that actually scatters photons). We mitigate these issues using two techniques, which we discuss in turn below. We note that, while similar techniques could be applied to the pair-production procedure, our simulations are optically thick to pair production, and the steady-state number of tracked photons is much smaller than the (always growing) number of particles. Thus, the QED physics involving photons, for us, can tolerate a less rigorously optimized implementation.

First, we speed up the assignment of the scattered photon energy – done in the context of Section 2.2.1 through equation (9). To do this, we consider the IC emission in the rest frame of the scattering particle. As we show below, this requires a larger number of random number draws per particle, but enables the use of simpler expressions. These, in turn, furnish an approximation to the scattering cross-section that effectively replaces equation (9), relieving its performance bottleneck. Our handling of the problem this way, presented in detail below, follows closely Levinson & Cerutti (2018) and Crinquand et al. (2020).

Before transforming to the rest frame of the potentially scattering particle, we isolate to an interaction with a single photon rather than with the entire isotropic radiation bath. This is accomplished by first drawing a random photon angle $\psi$ with respect to the direction of the particle’s velocity $c\beta$ according to the probability density function

$$p(\psi) \propto 1 - \beta \cos \psi, \quad (18)$$

which accounts for the relative lab-frame rate of encounters of the particle with photons incident from different directions. The selected photon is then boosted to the particle’s rest frame, where it has energy

$$\epsilon_0' = \gamma \epsilon_{ph} (1 - \beta \cos \psi). \quad (19)$$

Particle rest-frame quantities are primed in our convention.

We now evaluate the spectral (per unit final photon energy) and total (integrated over final photon energies) scattering rate for this interaction in the particle’s rest frame. We discuss after this derivation

how these two quantities are used by the code. In the primed frame, the scattering cross-section is that of ordinary Compton scattering (e.g. Blumenthal & Gould 1970):

$$\frac{d\sigma}{d\epsilon'_{\text{scat}} d\Omega'} = \frac{3\sigma_{\text{T}}}{16\pi} \left(\frac{\epsilon'_{\text{scat}}}{\epsilon'_{0}}\right)^2 \left(\frac{\epsilon'_{0}}{\epsilon'_{\text{scat}}} + \frac{\epsilon'_{\text{scat}}}{\epsilon'_{0}} - \sin^2 \Psi'\right) \times \delta \left[ \epsilon'_{\text{scat}} - \frac{\epsilon'_{0}}{1 + (\epsilon'_{0}/m_e c^2)(1 - \cos \Psi')} \right], \quad (20)$$

where $\Psi'$ is the angle through which the photon is scattered in the particle frame. The total scattering rate into final photon energies between $\epsilon'_{scat}$ and $\epsilon'_{scat} + d\epsilon'_{scat}$ is then

$$\frac{dN}{dt d\epsilon'_{\text{scat}}} = \frac{1}{\gamma} \frac{dN}{dt' d\epsilon'_{\text{scat}}} = cn_{\text{ph}} \int d\Omega' \frac{d\sigma}{d\epsilon'_{\text{scat}} d\Omega'}, \quad (21)$$

where we used $n'_{ph} = \gamma n_{ph}$. Defining $x \equiv \epsilon'_{scat}/\epsilon'_{0}$ and $y \equiv m_e c^2/\epsilon'_{0}$, the delta function only activates for $1/(1 + 2/y) \le x \le 1$, in which case the integral evaluates to

$$\frac{dN}{dt dx} = \epsilon'_{0} \frac{dN}{dt d\epsilon'_{\text{scat}}} = \frac{3\sigma_{\text{T}}}{8} \times y \left[ x + \frac{1}{x} + 2y \left(1 - \frac{1}{x}\right) + y^2 \left(\frac{1}{x^2} - \frac{2}{x} + 1\right) \right]. \quad (22)$$

The total scattering rate is then

$$\frac{dN}{dt} \equiv \int dx \frac{dN}{dt dx} = \frac{3\sigma_{\text{T}}cn_{\text{ph}}}{8} \left[ F(1, y) - F\left(\frac{1}{1+2/y}, y\right) \right], \quad (23)$$

where

$$F(x, y) \equiv y \left[ -\frac{y^2}{x} + xy(2+y) + \frac{x^2}{2} + (1-2y-2y^2) \ln x \right]. \quad (24)$$

These results are used by the code as follows. First, equation (23) is evaluated and multiplied by $\Delta t$ to determine the probability $p_{emit} = \Delta t \, dN/dt$ that the scattering occurs during the simulation time-step. This step replaces the evaluation of equation (5). On average, these two procedures are completely equivalent, but in this second method, we have traded the analytic integral over incident photon directions [which yields equation (5)] for a random Monte Carlo sampling over these directions.

Then, for particles that scatter photons, the final photon energy is obtained by drawing a random number $R$ and inverting the cumulative distribution function: that is, finding $x$ such that

$$\begin{aligned} \mathcal{R} &= \frac{3\sigma_{\text{T}}cn_{\text{ph}} [F(x, y) - F(1/(1+2/y), y)]}{8 \frac{dN}{dt}} \\ &= \frac{F(x, y) - F(1/(1+2/y), y)}{F(1, y) - F(1/(1+2/y), y)}. \end{aligned} \quad (25)$$

This replaces the inversion of the equivalent lab-frame cumulative distribution function (9) and, as we now show, constitutes the main advantage of this method. The issue with equation (9) is that its solution cannot be expressed analytically. Instead, the cumulative distribution function (or its inverse) must be stored as a table and consulted for each emitting particle – a costly procedure. However, the function $F(x, y)$ can be approximated, when $x \ll 1$, by only its $\ln x$ term. This approximation breaks down as $x$ approaches unity but, as it turns out, does not very strongly disturb the cumulative distribution function on the right-hand side of equation (25) even when used across all $x$. If one then approximates

$$F(x, y) \approx y (1 - 2y - 2y^2) \ln x, \quad (26)$$
the solution to equation (25) can be found as

$$x = \exp[-(1 - \mathcal{R}) \ln(1 + 2/y)]. \quad (27)$$

We have checked that the use of this approximation does not change any discernible aspects of our simulations. We have also conducted experiments (Appendix A) to show that it yields nearly the same average emitted photon energy as the exact cumulative distribution function. It therefore provides a powerful speed-up to the code by facilitating the otherwise impossible analytic evaluation of $x$ without compromising the important physics.

Once the rest-frame scattered photon energy $\epsilon'_{scat} = x\epsilon'_0$ is known, it can be boosted back to the lab frame as follows. First, one notes that, in the ultrarelativistic approximation $\gamma \gg 1$ relevant to this work, the incident photon approaches the particle nearly head-on in the primed frame. Thus, $\cos \Psi' = -\cos \psi'_{scat}$, where $\psi'_{scat}$ is the angle between the emerging photon's momentum and the particle's lab velocity. Then, the delta function in equation (20) can be used to write

$$\cos \psi'_{scat} = - \cos \Psi' = \frac{m_e c^2}{\epsilon'_{scat}} - \frac{m_e c^2}{\epsilon'_0} - 1, \quad (28)$$

which yields the lab-frame emerging photon energy via

$$\epsilon_{scat} = \gamma \epsilon'_{scat} (1 + \beta \cos \psi'_{scat}). \quad (29)$$

In addition to considering the IC scattering in the rest frame of each particle, which facilitates the powerful approximations (26) and (27), we also leverage the following second strategy to reduce the cost of the QED module. This standard optimization in the plasma physics literature is known as the *null collision method* (Rees 1968; Lin & Bardsley 1977; Boeuf & Marode 1982; Birdsall 1991). It speeds up the code drastically by avoiding the need to loop over all of the particles during the IC Monte Carlo emission step.

The technique exploits the fact that the probability $p_{emit}$ that an individual particle emits a photon in a given time-step is capped to a global maximum, given by the Thomson limit of equation (5), of $p_0 = c\sigma_T n_{ph} \Delta t$. This is a small number in our simulations: of order $10^{-3}$. Thus, instead of looping over all of the particles to determine whether fewer than 1 in every $1/p_0$ of them emits a photon, we randomly select a small fraction $p_0$ of all the particles, loop over this reduced subset, and exactly compensate the limited particle sample by enhancing the per-particle emission probability by the inverse factor $1/p_0$. The null collision method is analogous to the no-time-counter method (Bird 1989; Del Gaudio et al. 2020) except that the latter applies when the collision in question involves two computational particles rather than, as in our case, one computational particle and a fluid field (for us, the prescribed background seed photon bath).

The two techniques described above – the approximation of the cumulative distribution function on IC emission energies in equation (26) and the use of the null collision method – allow us to run QED-PIC simulations with negligible added cost per time-step (of order 10 per cent) taken by the QED module (orange panel in Fig. 1). The main costs are instead the accumulation of particles and photons in the simulation and the typically larger amount of data that one wishes to dump in QED runs. We note, however, that for the regime of reconnection studied in this work, the total number of simulated particles never grows by more than a factor of a few, and thus we do not need particle merging algorithms (Vranic et al. 2015), as implemented in other codes (e.g. OSIRIS, Grismayer et al. 2016; SMILEI, Derouillat et al. 2018; TRISTAN V2, Hakobyan et al. 2019), to regulate such growth at this stage.

## 3 SIMULATION SETUP

We perform pair-plasma simulations of relativistic collisionless magnetic reconnection using the radiative QED-outfitted (Section 2) electromagnetic PIC code ZELTRON (Cerutti et al. 2013; Cerutti & Werner 2019). The simulations are in a 2D box of physical dimensions $L_x \times L_y = L \times 2L$ and grid size $N_x \times N_y = N \times 2N = 7680 \times 15360$ and are run from $t = 0$ to just over $t = 4L/c$. Spatial dependence is only tracked in the x- and y-directions, but all vectors may have an out-of-plane z-component.

We adopt periodic box boundaries. Each of our simulations can therefore be interpreted as a single reconnection event that completely processes a flux reservoir of macroscopic scale $L$. This facilitates transparent connections to flaring situations in which reconnection suddenly and explosively liberates the magnetic energy contained in such a flux reservoir, and a major application of our results is to gamma-ray-flaring sources (Section 6). In addition, the periodic boundaries enable a straightforward computation of the pair-plasma yield per reconnection-processed particle (Section 5).

Our four main runs presented in Section 4 share the same setup and parameter values but differ in their modelled radiative effects. In particular, we present one case without any radiative cooling; one with purely continuous Thomson IC cooling (similar to, e.g., Werner et al. 2019; Mehlhaff et al. 2020; Sironi & Beloborodov 2020; Sridhar et al. 2021, 2023); one with fully general IC cooling (including the Klein–Nishina regime) but with pair production artificially turned off; and one with general IC cooling and self-consistent pair production. Synchrotron cooling is ignored in all runs. We describe the non-radiative aspects of our setup in Section 3.1 and move on to the radiative details in Section 3.2.

### 3.1 Non-radiative aspects of the setup

Our simulations begin with zero electric field. The initial magnetic field is force-free and undergoes reversals via current sheets located at $y_1 = L/2$ and $y_2 = 3L/2$. It has the form (Harrison & Neukirch 2009; Neukirch, Wilson & Harrison 2009)

$$\mathbf{B} = B_x(y)\hat{\mathbf{x}} + B_z(y)\hat{\mathbf{z}} = \pm B_0 \tanh\left(\frac{y - y_{1,2}}{\delta}\right)\hat{\mathbf{x}} + B_0\sqrt{\text{sech}^2\left(\frac{y - y_{1,2}}{\delta}\right) + \left(\frac{B_g}{B_0}\right)^2}\hat{\mathbf{z}}, \quad (30)$$

where $\hat{\mathbf{x}}$ and $\hat{\mathbf{z}}$ are unit vectors pointing in the respective x- and z-directions, the $+$ (−) sign is taken at $y_1$ ($y_2$), and $\delta$ is the half-thickness of the current sheets (see also Li et al. 2018; Guo et al. 2021; Goodbred & Liu 2022; Zhang et al. 2022; French et al. 2023, for a few relativistic reconnection studies employing the same initial magnetic field). In the reconnection upstream region far away from the current sheets ($|y - y_{1,2}| \gg \delta$), equation (30) reduces to a uniform field with in-plane component $\pm B_0\hat{\mathbf{x}}$ and out-of-plane guide-field component $B_g\hat{\mathbf{z}} = 0.15B_0\hat{\mathbf{z}}$. In addition to this modest guide field, which accompanies the upstream plasma into the reconnection layer, there is also a strong localized ($|y - y_{1,2}| \ll \delta$) initial out-of-plane field of peak strength $\sqrt{B_0^2 + B_g^2}$ that supplies the force-free magnetic field reversal. Both this strong localized out-of-plane field and the upstream guide field provide some pressure support to the plasma energized by reconnection as it radiatively cools down, which helps ensure the Debye length in the simulations remains well resolved.

Our initial fields satisfy the force-free condition, $\mathbf{J} \times \mathbf{B}/c = -\nabla \mathbf{B}^2/8\pi + \mathbf{B} \cdot \nabla \mathbf{B}/4\pi = 0$. Thus, and unlike the case of a Harris
the solution to equation (25) can be found as

$$x = \exp[-(1 - \mathcal{R}) \ln(1 + 2/y)]. \quad (27)$$

We have checked that the use of this approximation does not change any discernible aspects of our simulations. We have also conducted experiments (Appendix A) to show that it yields nearly the same average emitted photon energy as the exact cumulative distribution function. It therefore provides a powerful speed-up to the code by facilitating the otherwise impossible analytic evaluation of $x$ without compromising the important physics.

Once the rest-frame scattered photon energy $\epsilon'_{scat} = x\epsilon'_0$ is known, it can be boosted back to the lab frame as follows. First, one notes that, in the ultrarelativistic approximation $\gamma \gg 1$ relevant to this work, the incident photon approaches the particle nearly head-on in the primed frame. Thus, $\cos \Psi' = -\cos \psi'_{scat}$, where $\psi'_{scat}$ is the angle between the emerging photon's momentum and the particle's lab velocity. Then, the delta function in equation (20) can be used to write

$$\cos \psi'_{scat} = -\cos \Psi' = \frac{m_e c^2}{\epsilon'_{scat}} - \frac{m_e c^2}{\epsilon'_0} - 1, \quad (28)$$

which yields the lab-frame emerging photon energy via

$$\epsilon_{scat} = \gamma \epsilon'_{scat} (1 + \beta \cos \psi'_{scat}). \quad (29)$$

In addition to considering the IC scattering in the rest frame of each particle, which facilitates the powerful approximations (26) and (27), we also leverage the following second strategy to reduce the cost of the QED module. This standard optimization in the plasma physics literature is known as the *null collision method* (Rees 1968; Lin & Bardsley 1977; Boeuf & Marode 1982; Birdsall 1991). It speeds up the code drastically by avoiding the need to loop over all of the particles during the IC Monte Carlo emission step.

The technique exploits the fact that the probability $p_{emit}$ that an individual particle emits a photon in a given time-step is capped to a global maximum, given by the Thomson limit of equation (5), of $p_0 = c\sigma_T n_{ph} \Delta t$. This is a small number in our simulations: of order $10^{-3}$. Thus, instead of looping over all of the particles to determine whether fewer than 1 in every $1/p_0$ of them emits a photon, we randomly select a small fraction $p_0$ of all the particles, loop over this reduced subset, and exactly compensate the limited particle sample by enhancing the per-particle emission probability by the inverse factor $1/p_0$. The null collision method is analogous to the no-time-counter method (Bird 1989; Del Gaudio et al. 2020) except that the latter applies when the collision in question involves two computational particles rather than, as in our case, one computational particle and a fluid field (for us, the prescribed background seed photon bath).

The two techniques described above – the approximation of the cumulative distribution function on IC emission energies in equation (26) and the use of the null collision method – allow us to run QED-PIC simulations with negligible added cost per time-step (of order 10 per cent) taken by the QED module (orange panel in Fig. 1). The main costs are instead the accumulation of particles and photons in the simulation and the typically larger amount of data that one wishes to dump in QED runs. We note, however, that for the regime of reconnection studied in this work, the total number of simulated particles never grows by more than a factor of a few, and thus we do not need particle merging algorithms (Vranic et al. 2015), as implemented in other codes (e.g. OSIRIS, Grismayer et al. 2016; SMILEI, Derouillat et al. 2018; TRISTAN V2, Hakobyan et al. 2019), to regulate such growth at this stage.

## 3 SIMULATION SETUP

We perform pair-plasma simulations of relativistic collisionless magnetic reconnection using the radiative QED-outfitted (Section 2) electromagnetic PIC code ZELTRON (Cerutti et al. 2013; Cerutti & Werner 2019). The simulations are in a 2D box of physical dimensions $L_x \times L_y = L \times 2L$ and grid size $N_x \times N_y = N \times 2N = 7680 \times 15360$ and are run from $t = 0$ to just over $t = 4L/c$. Spatial dependence is only tracked in the $x$- and $y$-directions, but all vectors may have an out-of-plane $z$-component.

We adopt periodic box boundaries. Each of our simulations can therefore be interpreted as a single reconnection event that completely processes a flux reservoir of macroscopic scale $L$. This facilitates transparent connections to flaring situations in which reconnection suddenly and explosively liberates the magnetic energy contained in such a flux reservoir, and a major application of our results is to gamma-ray-flaring sources (Section 6). In addition, the periodic boundaries enable a straightforward computation of the pair-plasma yield per reconnection-processed particle (Section 5).

Our four main runs presented in Section 4 share the same setup and parameter values but differ in their modelled radiative effects. In particular, we present one case without any radiative cooling; one with purely continuous Thomson IC cooling (similar to, e.g., Werner et al. 2019; Mehlhaff et al. 2020; Sironi & Beloborodov 2020; Sridhar et al. 2021, 2023); one with fully general IC cooling (including the Klein–Nishina regime) but with pair production artificially turned off; and one with general IC cooling and self-consistent pair production. Synchrotron cooling is ignored in all runs. We describe the non-radiative aspects of our setup in Section 3.1 and move on to the radiative details in Section 3.2.

### 3.1 Non-radiative aspects of the setup

Our simulations begin with zero electric field. The initial magnetic field is force-free and undergoes reversals via current sheets located at $y_1 = L/2$ and $y_2 = 3L/2$. It has the form (Harrison & Neukirch 2009; Neukirch, Wilson & Harrison 2009)

$$\mathbf{B} = B_x(y)\hat{\mathbf{x}} + B_z(y)\hat{\mathbf{z}} = \pm B_0 \tanh\left(\frac{y - y_{1,2}}{\delta}\right) \hat{\mathbf{x}} + B_0 \sqrt{\sech^2\left(\frac{y - y_{1,2}}{\delta}\right) + \left(\frac{B_g}{B_0}\right)^2} \hat{\mathbf{z}}, \quad (30)$$

where $\hat{\mathbf{x}}$ and $\hat{\mathbf{z}}$ are unit vectors pointing in the respective $x$- and $z$-directions, the $+$ (−) sign is taken at $y_1$ ($y_2$), and $\delta$ is the half-thickness of the current sheets (see also Li et al. 2018; Guo et al. 2021; Goodbred & Liu 2022; Zhang et al. 2022; French et al. 2023, for a few relativistic reconnection studies employing the same initial magnetic field). In the reconnection upstream region far away from the current sheets ($|y - y_{1,2}| \gg \delta$), equation (30) reduces to a uniform field with in-plane component $\pm B_0\hat{\mathbf{x}}$ and out-of-plane guide-field component $B_g\hat{\mathbf{z}} = 0.15 B_0\hat{\mathbf{z}}$. In addition to this modest guide field, which accompanies the upstream plasma into the reconnection layer, there is also a strong localized ($|y - y_{1,2}| \ll \delta$) initial out-of-plane field of peak strength $\sqrt{B_0^2 + B_g^2}$ that supplies the force-free magnetic field reversal. Both this strong localized out-of-plane field and the upstream guide field provide some pressure support to the plasma energized by reconnection as it radiatively cools down, which helps ensure the Debye length in the simulations remains well resolved.

Our initial fields satisfy the force-free condition, $\mathbf{J} \times \mathbf{B}/c = -\nabla \mathbf{B}^2/8\pi + \mathbf{B} \cdot \nabla \mathbf{B}/4\pi = 0$. Thus, and unlike the case of a Harris
equilibrium (Kirk & Skjæraasen 2003), no added plasma pressure is needed inside the initial current layers. This allows us to start the simulations with a plasma of completely homogeneous initial (electron + positron) number density, $n_0$, and temperature, $T_0 = \theta_0 m_e c^2 = 24 m_e c^2$. Specifically, we use a relativistic Maxwell–Jüttner initial plasma distribution function. This setup avoids spurious pair production from radiation emitted by the hot particles that would be needed to supply the initial pressure support in a Harris current sheet.<fnref n="2" />

The number density, $n_0$, and initial reconnecting field strength, $B_0$, together define the cold magnetization,

$$ \sigma_{c,0} \equiv \frac{B_0^2}{4\pi n_0 m_e c^2}, \qquad (31) $$

equal to about twice the reconnecting magnetic field energy per particle. While $\sigma_{c,0}$ defines how much energy a typical reconnection-accelerated particle may acquire, another similar quantity, the hot magnetization,

$$ \sigma_{h,0} \equiv \frac{B_0^2}{4\pi w_0}, \qquad (32) $$

defines the magnetic dominance of the upstream region. Here, the initial plasma enthalpy density, $w_0$, can be written as $w_0 = u_0 + P_0$, where $u_0$ and $P_0$ are, respectively, the initial plasma internal energy density and pressure. For a non-relativistically cold initial temperature, $\theta_0 \ll 1$, the enthalpy density is dominated by rest-mass energy, $w_0 \simeq u_0 = n_0 m_e c^2$, and, hence, $\sigma_{h,0} \simeq \sigma_{c,0}$. For a relativistically hot plasma, $\theta_0 \gg 1$, the thermal kinetic energy and pressure dominate the enthalpy, $w_0 \simeq 4P_0 = 4\theta_0 n_0 m_e c^2$, and, in this case, $\sigma_{h,0} \simeq \sigma_{c,0}/4\theta_0 = 1/(2\beta_{pl})$, where $\beta_{pl}$ is plasma beta. For all regimes, $\sigma_{h,0} \le \sigma_{c,0}$, meaning that $\sigma_{h,0} \gg 1$ is a stronger condition than $\sigma_{c,0} \gg 1$. In fact, because $\sigma_{h,0}$ sets the plasma Alfvén speed $v_A = c\sqrt{\sigma_{h,0}/(1+\sigma_{h,0})}$, a high $\sigma_{h,0}$ means that the energy budget permits not just relativistic individual particles (which merely requires high $\sigma_{c,0}$), but also for the collective bulk motion of the plasma itself to become highly relativistic. In our simulations, we set $\sigma_{c,0} = 1.2 \times 10^5$ and $\sigma_{h,0} = \sigma_{c,0}/4\theta_0 = 1250$.

Though the force-free initial condition does not require any extra plasma density or temperature in the current sheets to balance the upstream magnetic pressure, the field-reversing currents must still be supplied. Correspondingly, we set a local fraction

$$ f_d(y) = \text{sech}(\tilde{y}/\delta) \sqrt{\frac{1 + (B_g/B_0)^2}{1 + [B_g \cosh(\tilde{y}/\delta)/B_0]^2}} \qquad (33) $$

of the positrons in motion at a drift velocity

$$ c\boldsymbol{\beta}_{d,i}(y) = \frac{c\beta_d}{\sqrt{1 + (B_g/B_0)^2}} \times \left[ \mp\sqrt{\text{sech}^2(\tilde{y}/\delta) + (B_g/B_0)^2} \hat{z} - \tanh(\tilde{y}/\delta) \hat{x} \right] \qquad (34) $$

so that they carry half the field-reversing current

$$ \mathbf{J} = c\nabla \times \mathbf{B}/4\pi = \text{ecn}_0 f_d(y) \boldsymbol{\beta}_{d,i}(y), \qquad (35) $$

where the $\mp$ corresponds to equation (30), $e$ is the positron charge, and $\tilde{y} \equiv y - y_{1,2}$. The other half of the current is carried by the initial electrons, which are assigned $\boldsymbol{\beta}_{d,e} = -\boldsymbol{\beta}_{d,i}$. The drifting particles

follow a drifting Maxwell–Jüttner distribution function with initial temperature, $T_0$, defined in their local [boosted by $\boldsymbol{\beta}_{d,i/e}(y)$] rest frame.

The force-free setup ties the current sheet half-thickness $\delta$ to the other length-scales in the problem as follows. The drifting plasma supplies a current density $\text{ecn}_0 \beta_d \le \text{ecn}_0$, whereas the current needed at the heart of the layer is, according to Ampère's law (35), $cB_0/4\pi\delta$. This means that

$$ \delta = \frac{B_0}{4\pi en_0 \beta_d} = \frac{\sigma_{c,0} \rho_0}{\beta_d} \ge \sigma_{c,0} \rho_0, \qquad (36) $$

where we have introduced the nominal gyroradius,

$$ \rho_0 \equiv \frac{m_e c^2}{e B_0}. \qquad (37) $$

Equation (36) demands that the current sheet half-thickness be limited to the typical gyroradius, $\sigma_{c,0} \rho_0$, of reconnection-energized particles. Therefore, in order to achieve a high aspect ratio, $L/\delta$, while also complying with the other demanding radiative constraints described later (Section 3.2), we assign $\beta_d = 0.9$. This corresponds to $\delta = 1.1 \sigma_{c,0} \rho_0$ and to an initial drifting-particle Lorentz factor of 2.3 (still much less than $\theta_0$).

The current sheet aspect ratio, $L/\delta = L\beta_d/\sigma_{c,0} \rho_0 \simeq L/\sigma_{c,0} \rho_0$, represents not just a ratio of length-scales but also one of particle energy scales. This is because, during reconnection, some particles are swept into the vicinity of an X-point (X-line in 3D), which is a region wherein the in-plane magnetic field reconnects. There, they become unmagnetized and linearly accelerated by the out-of-plane reconnection electric field, $E_{\text{rec}} = \beta_{\text{rec}}B_0 \simeq 0.1B_0$, where $\beta_{\text{rec}} \equiv 0.1v_A/c \simeq 0.1$ is the fiducial collisionless reconnection rate. As described by Werner et al. (2016), in systems with modest aspect ratios, for which the reconnection layer is dominated by a single X-point, particles can be accelerated by $E_{\text{rec}}$ until, after travelling an out-of-plane distance $\sim L$, they escape the acceleration region with final system-size-limited Lorentz factor

$$ \gamma_{\text{max}} \equiv \frac{0.1 e B_0 L}{m_e c^2} = \frac{0.1 L}{\rho_0}. \qquad (38) $$

However, for larger systems, the reconnection layer tears into a hierarchical chain of plasma-filled magnetic islands, or plasmons, studded with many X-points. Then, particles cannot travel the whole system size before escaping a given acceleration zone. Instead, at least in 2D, they may travel a distance comparable to that between the smallest-scale plasmons, which subsequently capture the particles, limiting them to energies of order the intrinsic reconnection X-point acceleration Lorentz factor (cf. Sironi et al. 2016; Werner et al. 2016; Uzdensky 2022),

$$ \gamma_X \equiv 4\sigma_{c,0}. \qquad (39) $$

The aspect ratio becomes large enough that the X-point acceleration channel is intrinsically limited by the self-consistent evolution of the plasmoid chain rather than by the modest size of the system when these two energy scales cross each other: $\gamma_{\text{max}}/\gamma_X \ge 1 \Rightarrow L/\sigma_{c,0} \rho_0 \ge 40$ (Werner et al. 2016). Thus, the requirement to have a large aspect ratio $L/\delta \gg 1$ is synonymous with having a healthy separation between $\gamma_{\text{max}}$ and the particle energies, $\sim \gamma_X$, at which reconnection X-point acceleration saturates. We adopt $\gamma_{\text{max}} = 9.1\sigma_{c,0}$, which corresponds to $L = 91\sigma_{c,0} \rho_0$. This meets the fiducial $\gamma_{\text{max}} > \gamma_X$ criterion but, for reasons described below, yields $\gamma_{\text{max}}/\sigma_{c,0}$ smaller than typical for simulations of the numerical size, $N_x \times N_y = 7680 \times 15360$, that we present.

<footnote n="2">²While such pair production could also be avoided in a traditional Harris equilibrium by specifying a high plasma overdensity in the initial current layer, this would reduce the initial Debye length of the layer plasma, and would thus create problems for the Debye length resolution in our simulations.</footnote>
equilibrium (Kirk & Skjæraasen 2003), no added plasma pressure is needed inside the initial current layers. This allows us to start the simulations with a plasma of completely homogeneous initial (electron + positron) number density, $n_0$, and temperature, $T_0 = \theta_0 m_e c^2 = 24 m_e c^2$. Specifically, we use a relativistic Maxwell–Jüttner initial plasma distribution function. This setup avoids spurious pair production from radiation emitted by the hot particles that would be needed to supply the initial pressure support in a Harris current sheet.<fnref n="2" />

The number density, $n_0$, and initial reconnecting field strength, $B_0$, together define the cold magnetization,

$$ \sigma_{c,0} \equiv \frac{B_0^2}{4\pi n_0 m_e c^2}, \qquad (31) $$

equal to about twice the reconnecting magnetic field energy per particle. While $\sigma_{c,0}$ defines how much energy a typical reconnecting-accelerated particle may acquire, another similar quantity, the hot magnetization,

$$ \sigma_{h,0} \equiv \frac{B_0^2}{4\pi w_0}, \qquad (32) $$

defines the magnetic dominance of the upstream region. Here, the initial plasma enthalpy density, $w_0$, can be written as $w_0 = u_0 + P_0$, where $u_0$ and $P_0$ are, respectively, the initial plasma internal energy density and pressure. For a non-relativistically cold initial temperature, $\theta_0 \ll 1$, the enthalpy density is dominated by rest-mass energy, $w_0 \simeq u_0 = n_0 m_e c^2$, and, hence, $\sigma_{h,0} \simeq \sigma_{c,0}$. For a relativistically hot plasma, $\theta_0 \gg 1$, the thermal kinetic energy and pressure dominate the enthalpy, $w_0 \simeq 4P_0 = 4\theta_0 n_0 m_e c^2$, and, in this case, $\sigma_{h,0} \simeq \sigma_{c,0}/4\theta_0 = 1/(2\beta_{pl})$, where $\beta_{pl}$ is plasma beta. For all regimes, $\sigma_{h,0} \le \sigma_{c,0}$, meaning that $\sigma_{h,0} \gg 1$ is a stronger condition than $\sigma_{c,0} \gg 1$. In fact, because $\sigma_{h,0}$ sets the plasma Alfvén speed $v_A = c\sqrt{\sigma_{h,0}/(1+\sigma_{h,0})}$, a high $\sigma_{h,0}$ means that the energy budget permits not just relativistic individual particles (which merely requires high $\sigma_{c,0}$), but also for the collective bulk motion of the plasma itself to become highly relativistic. In our simulations, we set $\sigma_{c,0} = 1.2 \times 10^5$ and $\sigma_{h,0} = \sigma_{c,0}/4\theta_0 = 1250$.

Though the force-free initial condition does not require any extra plasma density or temperature in the current sheets to balance the upstream magnetic pressure, the field-reversing currents must still be supplied. Correspondingly, we set a local fraction

$$ f_d(y) = \text{sech}(\tilde{y}/\delta) \sqrt{\frac{1 + (B_g/B_0)^2}{1 + [B_g \cosh(\tilde{y}/\delta)/B_0]^2}} \qquad (33) $$

of the positrons in motion at a drift velocity

$$ c\boldsymbol{\beta}_{d,i}(y) = \frac{c\beta_d}{\sqrt{1 + (B_g/B_0)^2}} \times \left[ \mp \sqrt{\text{sech}^2(\tilde{y}/\delta) + (B_g/B_0)^2} \hat{z} - \tanh(\tilde{y}/\delta) \hat{x} \right] \qquad (34) $$

so that they carry half the field-reversing current

$$ \mathbf{J} = c\nabla \times \mathbf{B}/4\pi = \text{ecn}_0 f_d(y) \boldsymbol{\beta}_{d,i}(y), \qquad (35) $$

where the $\mp$ corresponds to equation (30), $e$ is the positron charge, and $\tilde{y} \equiv y - y_{1,2}$. The other half of the current is carried by the initial electrons, which are assigned $\boldsymbol{\beta}_{d,e} = -\boldsymbol{\beta}_{d,i}$. The drifting particles

follow a drifting Maxwell–Jüttner distribution function with initial temperature, $T_0$, defined in their local [boosted by $\boldsymbol{\beta}_{d,i/e}(y)$] rest frame.

The force-free setup ties the current sheet half-thickness $\delta$ to the other length-scales in the problem as follows. The drifting plasma supplies a current density $\text{ecn}_0 \beta_d \le \text{ecn}_0$, whereas the current needed at the heart of the layer is, according to Ampère's law (35), $cB_0/4\pi\delta$. This means that

$$ \delta = \frac{B_0}{4\pi en_0 \beta_d} = \frac{\sigma_{c,0} \rho_0}{\beta_d} \ge \sigma_{c,0} \rho_0, \qquad (36) $$

where we have introduced the nominal gyroradius,

$$ \rho_0 \equiv \frac{m_e c^2}{e B_0}. \qquad (37) $$

Equation (36) demands that the current sheet half-thickness be limited to the typical gyroradius, $\sigma_{c,0} \rho_0$, of reconnection-energized particles. Therefore, in order to achieve a high aspect ratio, $L/\delta$, while also complying with the other demanding radiative constraints described later (Section 3.2), we assign $\beta_d = 0.9$. This corresponds to $\delta = 1.1 \sigma_{c,0} \rho_0$ and to an initial drifting-particle Lorentz factor of 2.3 (still much less than $\theta_0$).

The current sheet aspect ratio, $L/\delta = L \beta_d / \sigma_{c,0} \rho_0 \simeq L / \sigma_{c,0} \rho_0$, represents not just a ratio of length-scales but also one of particle energy scales. This is because, during reconnection, some particles are swept into the vicinity of an X-point (X-line in 3D), which is a region wherein the in-plane magnetic field reconnects. There, they become unmagnetized and linearly accelerated by the out-of-plane reconnection electric field, $E_{\text{rec}} = \beta_{\text{rec}} B_0 \simeq 0.1 B_0$, where $\beta_{\text{rec}} \equiv 0.1 v_A / c \simeq 0.1$ is the fiducial collisionless reconnection rate. As described by Werner et al. (2016), in systems with modest aspect ratios, for which the reconnection layer is dominated by a single X-point, particles can be accelerated by $E_{\text{rec}}$ until, after travelling an out-of-plane distance $\sim L$, they escape the acceleration region with final system-size-limited Lorentz factor

$$ \gamma_{\text{max}} \equiv \frac{0.1 e B_0 L}{m_e c^2} = \frac{0.1 L}{\rho_0}. \qquad (38) $$

However, for larger systems, the reconnection layer tears into a hierarchical chain of plasma-filled magnetic islands, or plasmons, studded with many X-points. Then, particles cannot travel the whole system size before escaping a given acceleration zone. Instead, at least in 2D, they may travel a distance comparable to that between the smallest-scale plasmons, which subsequently capture the particles, limiting them to energies of order the intrinsic reconnection X-point acceleration Lorentz factor (cf. Sironi et al. 2016; Werner et al. 2016; Uzdensky 2022),

$$ \gamma_X \equiv 4\sigma_{c,0}. \qquad (39) $$

The aspect ratio becomes large enough that the X-point acceleration channel is intrinsically limited by the self-consistent evolution of the plasmoid chain rather than by the modest size of the system when these two energy scales cross each other: $\gamma_{\text{max}}/\gamma_X \ge 1 \Rightarrow L/\sigma_{c,0} \rho_0 \ge 40$ (Werner et al. 2016). Thus, the requirement to have a large aspect ratio $L/\delta \gg 1$ is synonymous with having a healthy separation between $\gamma_{\text{max}}$ and the particle energies, $\sim \gamma_X$, at which reconnection X-point acceleration saturates. We adopt $\gamma_{\text{max}} = 9.1 \sigma_{c,0}$, which corresponds to $L = 91 \sigma_{c,0} \rho_0$. This meets the fiducial $\gamma_{\text{max}} > \gamma_X$ criterion but, for reasons described below, yields $\gamma_{\text{max}}/\sigma_{c,0}$ smaller than typical for simulations of the numerical size, $N_x \times N_y = 7680 \times 15360$, that we present.

<footnote n="2">²While such pair production could also be avoided in a traditional Harris equilibrium by specifying a high plasma overdensity in the initial current layer, this would reduce the initial Debye length of the layer plasma, and would thus create problems for the Debye length resolution in our simulations.</footnote>
<table><caption>**Table 1.** Simulation parameters used in this study. Non-radiative parameters (described in Section 3.1) are presented in the upper part of the table and radiative parameters (described in Section 3.2) below. The radiative parameters $\gamma_{\text{cool}}$ and $\gamma_{\text{rad,T}}$ apply only to the three simulations with radiative cooling, $\gamma_{\text{KN}}$ applies only to the two simulations with fully general IC cooling (including Klein–Nishina effects), and $\tau_{\gamma\gamma}$ applies only to the one simulation with pair production. The expression for $\sigma_{\text{h,0}}$ assumes $\theta_0 \gg 1$.</caption><thead><tr><th>Parameter</th><th>Symbol (= definition)</th><th>Value (= equivalent)</th></tr></thead><tbody><tr><td>Upstream reconnecting field strength</td><td>$B_0$</td><td></td></tr><tr><td>Nominal gyroradius</td><td>$\rho_0 = m_e c^2 / e B_0$</td><td></td></tr><tr><td>Initial upstream density</td><td>$n_0$</td><td></td></tr><tr><td>Initial cold magnetization</td><td>$\sigma_{c,0} = B_0^2 / (4\pi n_0 m_e c^2)$</td><td>120 000</td></tr><tr><td>Initial hot magnetization</td><td>$\sigma_{h,0} = B_0^2 / (16\pi n_0 \theta_0 m_e c^2)$</td><td>1250</td></tr><tr><td>Initial upstream temperature</td><td>$\theta_0 = T_0 / m_e c^2$</td><td>$2 \times 10^{-4} \sigma_{c,0} = 24$</td></tr><tr><td>System size</td><td>$L_x = L$</td><td>$91 \sigma_{c,0} \rho_0$</td></tr><tr><td>System-size-limited Lorentz factor</td><td>$\gamma_{\text{max}} = 0.1 L / \rho_0$</td><td>$9.1 \sigma_{c,0} = 1.1 \times 10^6$</td></tr><tr><td>Guide field</td><td>$B_g$</td><td>$0.15 B_0$</td></tr><tr><td>Layer drift velocity</td><td>$\beta_d c$</td><td>$0.9 c$</td></tr><tr><td>Layer half-thickness</td><td>$\delta = \sigma_{c,0} \rho_0 / \beta_d$</td><td>$1.1 \sigma_{c,0} \rho_0$</td></tr><tr><td>Cell size</td><td>$\Delta x, \Delta y$</td><td>$\sigma_{c,0} \rho_0 / 85$</td></tr><tr><td>Time-step</td><td>$\Delta t$</td><td>$3 \Delta x / \sqrt{2c}$</td></tr><tr><td>Grid size</td><td>$N_x = N$</td><td>7680</td></tr><tr><td>Computational particles per cell</td><td></td><td>20</td></tr><tr><td>Soft photon energy density</td><td>$U_{\text{ph}}$</td><td></td></tr><tr><td>Soft photon energy</td><td>$\epsilon_{\text{ph}}$</td><td></td></tr><tr><td>Soft photon number density</td><td>$n_{\text{ph}} = U_{\text{ph}} / \epsilon_{\text{ph}}$</td><td></td></tr><tr><td>Nominal cooling Lorentz factor</td><td>$\gamma_{\text{cool}} = 3 m_e c^2 / (4 U_{\text{ph}} \sigma_T L)$</td><td>$3.6 \times 10^{-3} \sigma_{c,0} = 430$</td></tr><tr><td>Thomson IC-limited Lorentz factor</td><td>$\gamma_{\text{rad,T}} = (0.3 e B_0 / 4 \sigma_T U_{\text{ph}})^{1/2} = (\gamma_{\text{max}} \gamma_{\text{cool}})^{1/2}$</td><td>$0.18 \sigma_{c,0} = 2.2 \times 10^4$</td></tr><tr><td>Klein–Nishina Lorentz factor</td><td>$\gamma_{\text{KN}} = m_e c^2 / (4 \epsilon_{\text{ph}})$</td><td>$0.046 \sigma_{c,0} = 5500$</td></tr><tr><td>Box pair-production optical depth</td><td>$\tau_{\gamma\gamma} = n_{\text{ph}} \sigma_T L / 5 = 3 \gamma_{\text{KN}} / 5 \gamma_{\text{cool}}$</td><td>7.7</td></tr></tbody></table>

Our grid resolution $\Delta x = \Delta y$ is set by the need to resolve the initial Debye length,

$$ \lambda_{\text{D,0}} = \sqrt{\frac{\theta_0 m_e c^2}{4\pi e^2 n_0}} = \sqrt{\theta_0 \sigma_{c,0}} \rho_0, \quad (40) $$

and we set $\Delta x = \lambda_{\text{D,0}}/1.2$ in all runs. Because we operate in the highly relativistic regime, $\sigma_{c,0} \ge \sigma_{h,0} \gg 1$, the Debye length is much smaller than the typical energized particle's gyroradius, $\sigma_{c,0} \rho_0 \sim \delta$. This is largely why we cannot afford a larger $\gamma_{\text{max}}/\sigma_{c,0}$—our choice of $\sigma_{\text{h,0}} = 1250$ induces a large separation between the plasma microscales, $\sigma_{c,0} \rho_0 / \lambda_{\text{D,0}} = \sqrt{\sigma_{c,0}/\theta_0} \approx 2\sqrt{\sigma_{\text{h,0}}} \sim 70$, which occupies much of our grid resolution to resolve.<fnref n="3" /> We do, however, under-resolve the gyroradii of particles with energies less than $\gamma_{\Delta x} \equiv (\theta_0 \sigma_{c,0})^{1/2}/1.2$, including the upstream particles (energies $\sim \theta_0 \ll \gamma_{\Delta x}$). We do not observe any strong artefacts of this in our results, and the energy in our simulations is conserved to 1 per cent or better except where noted (i.e. in Table 2).

Given the cell size, $\Delta x$, we employ a time-step $\Delta t = 3\Delta x/\sqrt{2c} \approx 0.5\Delta x/c$: slightly smaller than the maximum allowed Courant–Friedrichs–Lewy step $\Delta x/\sqrt{2c}$. While this does not introduce significant numerical dispersion (see fig. 2 of Greenwood et al. 2004), it typically improves the simulations' energy conservation by about 0.5 to 1 percentage points (see, e.g., the simulations of Table 2 that employ a larger time-step and exhibit slightly poorer energy conservation). We initially place 20 computational particles (electrons + positrons) per cell, except within a few $\delta$ of $\gamma_{1,2}$ where we place 40 per cell. The particle currents are deposited to the grid using the charge-conserving method of Esirkepov (2001).

The principal non-radiative simulation parameters described in this section and their values are summarized in the top part of Table 1.

## 3.2 Radiative aspects of the setup

The two principle radiative parameters in our simulations are the total energy density, $U_{\text{ph}}$, and monochromatic single-photon energy, $\epsilon_{\text{ph}}$, of the IC seed photons [equation (4)]. These photons are not tracked, but provide a static, homogeneous, isotropic target population for the charged particles to scatter. Though only the two numbers $U_{\text{ph}}$ and $\epsilon_{\text{ph}}$ need to be prescribed to fully specify the radiative physics, they introduce a variety of derived energy- and length-scales into the problem, and, hence, divide the parameter space into many different regimes. We first provide a brief overview of these regimes in Section 3.2.1, summarizing the more detailed discussion from our earlier work, Mehlhaff et al. (2021). Afterwards, in Section 3.2.2, we discuss how these parameters are chosen for our simulations (as in Table 1).

### 3.2.1 Reconnection scales introduced by radiative physics

For particles with Lorentz factors $\gamma \ll \gamma_{\text{KN}} \equiv m_e c^2 / 4\epsilon_{\text{ph}}$, IC emission proceeds in the Thomson regime, where typical scattered photons are enhanced to energies $\sim \gamma^2 \epsilon_{\text{ph}}$, up to a maximum of $4\gamma^2 \epsilon_{\text{ph}}$. Because no individual photon robs the particle of a significant fraction of its energy, cooling proceeds continuously and is modelled (Section 2.1) as a radiative drag force, $f_{\text{rad}}$. The total power radiated by a given particle in the Thomson regime is [cf. equation (13)]

$$ P_{\text{T}}(\gamma) = |c f_{\text{rad}} \cdot \boldsymbol{\beta}| = \frac{4}{3} \sigma_{\text{T}} c \gamma^2 \beta^2 U_{\text{ph}}, \quad (41) $$

where $c\beta$ is the 3-velocity of the particle.

However, at sufficiently high $\gamma$, the maximum Thomson emission energy, $4\gamma^2 \epsilon_{\text{ph}}$, rivals the scattering particle's energy $\gamma m_e c^2$ (the two

<footnote n="3">³While $\sigma_{c,0} \rho_0 / \lambda_{\text{D,0}} = \sqrt{\sigma_{c,0}/\theta_0}$ is generally true, the approximation, $\sqrt{\sigma_{c,0}/\theta_0} \approx 2\sqrt{\sigma_{\text{h,0}}}$, assumes $\theta_0 \gg 1$.</footnote>
<table><caption>**Table 1.** Simulation parameters used in this study. Non-radiative parameters (described in Section 3.1) are presented in the upper part of the table and radiative parameters (described in Section 3.2) below. The radiative parameters $\gamma_{\text{cool}}$ and $\gamma_{\text{rad, T}}$ apply only to the three simulations with radiative cooling, $\gamma_{\text{KN}}$ applies only to the two simulations with fully general IC cooling (including Klein–Nishina effects), and $\tau_{\gamma\gamma}$ applies only to the one simulation with pair production. The expression for $\sigma_{\text{h,0}}$ assumes $\theta_0 \gg 1$.</caption><thead><tr><th>Parameter</th><th>Symbol (= definition)</th><th>Value (= equivalent)</th></tr></thead><tbody><tr><td>Upstream reconnecting field strength</td><td>$B_0$</td><td></td></tr><tr><td>Nominal gyroradius</td><td>$\rho_0 = m_e c^2 / e B_0$</td><td></td></tr><tr><td>Initial upstream density</td><td>$n_0$</td><td></td></tr><tr><td>Initial cold magnetization</td><td>$\sigma_{\text{c,0}} = B_0^2 / (4\pi n_0 m_e c^2)$</td><td>120 000</td></tr><tr><td>Initial hot magnetization</td><td>$\sigma_{\text{h,0}} = B_0^2 / (16\pi n_0 \theta_0 m_e c^2)$</td><td>1250</td></tr><tr><td>Initial upstream temperature</td><td>$\theta_0 = T_0 / m_e c^2$</td><td>$2 \times 10^{-4} \sigma_{\text{c,0}} = 24$</td></tr><tr><td>System size</td><td>$L_x = L$</td><td>$91 \sigma_{\text{c,0}} \rho_0$</td></tr><tr><td>System-size-limited Lorentz factor</td><td>$\gamma_{\text{max}} = 0.1 L / \rho_0$</td><td>$9.1 \sigma_{\text{c,0}} = 1.1 \times 10^6$</td></tr><tr><td>Guide field</td><td>$B_g$</td><td>$0.15 B_0$</td></tr><tr><td>Layer drift velocity</td><td>$\beta_d c$</td><td>$0.9 c$</td></tr><tr><td>Layer half-thickness</td><td>$\delta = \sigma_{\text{c,0}} \rho_0 / \beta_d$</td><td>$1.1 \sigma_{\text{c,0}} \rho_0$</td></tr><tr><td>Cell size</td><td>$\Delta x, \Delta y$</td><td>$\sigma_{\text{c,0}} \rho_0 / 85$</td></tr><tr><td>Time-step</td><td>$\Delta t$</td><td>$3 \Delta x / \sqrt{2c}$</td></tr><tr><td>Grid size</td><td>$N_x = N$</td><td>7680</td></tr><tr><td>Computational particles per cell</td><td></td><td>20</td></tr><tr><td>Soft photon energy density</td><td>$U_{\text{ph}}$</td><td></td></tr><tr><td>Soft photon energy</td><td>$\epsilon_{\text{ph}}$</td><td></td></tr><tr><td>Soft photon number density</td><td>$n_{\text{ph}} = U_{\text{ph}} / \epsilon_{\text{ph}}$</td><td>$3.6 \times 10^{-3} \sigma_{\text{c,0}} = 430$</td></tr><tr><td>Nominal cooling Lorentz factor</td><td>$\gamma_{\text{cool}} = 3 m_e c^2 / (4 U_{\text{ph}} \sigma_T L)$</td><td>$0.18 \sigma_{\text{c,0}} = 2.2 \times 10^4$</td></tr><tr><td>Thomson IC-limited Lorentz factor</td><td>$\gamma_{\text{rad, T}} = (0.3 e B_0 / 4 \sigma_T U_{\text{ph}})^{1/2} = (\gamma_{\text{max}} \gamma_{\text{cool}})^{1/2}$</td><td>$0.046 \sigma_{\text{c,0}} = 5500$</td></tr><tr><td>Klein–Nishina Lorentz factor</td><td>$\gamma_{\text{KN}} = m_e c^2 / (4 \epsilon_{\text{ph}})$</td><td></td></tr><tr><td>Box pair-production optical depth</td><td>$\tau_{\gamma\gamma} = n_{\text{ph}} \sigma_T L / 5 = 3 \gamma_{\text{KN}} / 5 \gamma_{\text{cool}}$</td><td>7.7</td></tr></tbody></table>

Our grid resolution $\Delta x = \Delta y$ is set by the need to resolve the initial Debye length,

$$ \lambda_{\text{D,0}} = \sqrt{\frac{\theta_0 m_e c^2}{4\pi e^2 n_0}} = \sqrt{\theta_0 \sigma_{\text{c,0}}} \rho_0, \qquad (40) $$

and we set $\Delta x = \lambda_{\text{D,0}}/1.2$ in all runs. Because we operate in the highly relativistic regime, $\sigma_{\text{c,0}} \ge \sigma_{\text{h,0}} \gg 1$, the Debye length is much smaller than the typical energized particle’s gyroradius, $\sigma_{\text{c,0}}\rho_0 \sim \delta$. This is largely why we cannot afford a larger $\gamma_{\text{max}}/\sigma_{\text{c,0}}$ – our choice of $\sigma_{\text{h,0}} = 1250$ induces a large separation between the plasma microscales, $\sigma_{\text{c,0}}\rho_0/\lambda_{\text{D,0}} = \sqrt{\sigma_{\text{c,0}}/\theta_0} \approx 2\sqrt{\sigma_{\text{h,0}}} \sim 70$, which occupies much of our grid resolution to resolve.<fnref n="3" /> We do, however, under-resolve the gyroradii of particles with energies less than $\gamma_{\Delta x} \equiv (\theta_0 \sigma_{\text{c,0}})^{1/2}/1.2$, including the upstream particles (energies $\sim \theta_0 \ll \gamma_{\Delta x}$). We do not observe any strong artefacts of this in our results, and the energy in our simulations is conserved to 1 per cent or better except where noted (i.e. in Table 2).

Given the cell size, $\Delta x$, we employ a time-step $\Delta t = 3\Delta x/\sqrt{2c} \approx 0.5\Delta x/c$: slightly smaller than the maximum allowed Courant–Friedrichs–Lewy step $\Delta x/\sqrt{2c}$. While this does not introduce significant numerical dispersion (see fig. 2 of Greenwood et al. 2004), it typically improves the simulations’ energy conservation by about 0.5 to 1 percentage points (see, e.g., the simulations of Table 2 that employ a larger time-step and exhibit slightly poorer energy conservation). We initially place 20 computational particles (electrons + positrons) per cell, except within a few $\delta$ of $\gamma_{1,2}$ where we place 40 per cell. The particle currents are deposited to the grid using the charge-conserving method of Esirkepov (2001).

The principal non-radiative simulation parameters described in this section and their values are summarized in the top part of Table 1.

## 3.2 Radiative aspects of the setup

The two principle radiative parameters in our simulations are the total energy density, $U_{\text{ph}}$, and monochromatic single-photon energy, $\epsilon_{\text{ph}}$, of the IC seed photons [equation (4)]. These photons are not tracked, but provide a static, homogeneous, isotropic target population for the charged particles to scatter. Though only the two numbers $U_{\text{ph}}$ and $\epsilon_{\text{ph}}$ need to be prescribed to fully specify the radiative physics, they introduce a variety of derived energy- and length-scales into the problem, and, hence, divide the parameter space into many different regimes. We first provide a brief overview of these regimes in Section 3.2.1, summarizing the more detailed discussion from our earlier work, Mehlhaff et al. (2021). Afterwards, in Section 3.2.2, we discuss how these parameters are chosen for our simulations (as in Table 1).

### 3.2.1 Reconnection scales introduced by radiative physics

For particles with Lorentz factors $\gamma \ll \gamma_{\text{KN}} \equiv m_e c^2 / 4\epsilon_{\text{ph}}$, IC emission proceeds in the Thomson regime, where typical scattered photons are enhanced to energies $\sim \gamma^2 \epsilon_{\text{ph}}$, up to a maximum of $4\gamma^2 \epsilon_{\text{ph}}$. Because no individual photon robs the particle of a significant fraction of its energy, cooling proceeds continuously and is modelled (Section 2.1) as a radiative drag force, $f_{\text{rad}}$. The total power radiated by a given particle in the Thomson regime is [cf. equation (13)]

$$ P_{\text{T}}(\gamma) = |c f_{\text{rad}} \cdot \boldsymbol{\beta}| = \frac{4}{3} \sigma_{\text{T}} c \gamma^2 \beta^2 U_{\text{ph}}, \qquad (41) $$

where $c\beta$ is the 3-velocity of the particle.

However, at sufficiently high $\gamma$, the maximum Thomson emission energy, $4\gamma^2 \epsilon_{\text{ph}}$, rivals the scattering particle’s energy $\gamma m_e c^2$ (the two

<footnote n="3">³While $\sigma_{\text{c,0}}\rho_0/\lambda_{\text{D,0}} = \sqrt{\sigma_{\text{c,0}}/\theta_0}$ is generally true, the approximation, $\sqrt{\sigma_{\text{c,0}}/\theta_0} \approx 2\sqrt{\sigma_{\text{h,0}}}$, assumes $\theta_0 \gg 1$.</footnote>
are formally equal at $\gamma_{KN} = m_e c^2/4\epsilon_{ph})$. Then, particles pass into the Klein–Nishina regime where they lose energy in discrete quanta, ~$\gamma m_e c^2$. Here, the cross-section [e.g. equations (7) and (20)] needs to be described by QED and gives rise to an average power radiated per particle of

$$P_{\text{IC}}(\gamma) = P_{\text{T}}(\gamma) f_{\text{KN}}(\gamma/\gamma_{\text{KN}}), \quad (42)$$

where $f_{KN}(q)$, with $q = \gamma/\gamma_{KN}$, is the dimensionless function (cf. Jones 1968; Nalewajko, Yuan & Chruścińska 2018; Mehlhaff et al. 2021)

$$f_{\text{KN}}(q) = \frac{9}{q^3} \left[ \left( \frac{q}{2} + 6 + \frac{6}{q} \right) \ln(1+q) - \frac{1}{(1+q)^2} \left( \frac{11}{12}q^3 + 6q^2 + 9q + 4 \right) - 2 + 2 \text{Li}_2(-q) \right] \quad (43)$$

and $\text{Li}_2(q)$ is the dilogarithm. As necessary to recover the Thomson limit, $f_{KN}(q \ll 1) \to 1$. In the opposite, deep Klein–Nishina regime, $f_{KN}(q \gg 1) \simeq (9/2q^2)[\ln(q) - 11/6]$. Equations (41) and (42) define the respective Thomson-limit and general-case IC cooling times,

$$\begin{aligned} t_{\text{cool},T}(\gamma) &= \frac{\gamma m_e c^2}{P_T(\gamma)} \approx \frac{3m_e c}{4\sigma_T U_{\text{ph}}\gamma} \\ &= \frac{L}{c} \frac{\gamma_{\text{cool}}}{\gamma} \end{aligned} \quad (44)$$

and

$$\begin{aligned} t_{\text{cool},\text{IC}}(\gamma) &= \frac{\gamma m_e c^2}{P_{\text{IC}}(\gamma)} \approx \frac{3m_e c}{4\sigma_T U_{\text{ph}}\gamma f_{\text{KN}}(\gamma/\gamma_{\text{KN}})} \\ &= \frac{L}{c} \frac{\gamma_{\text{cool}}}{\gamma f_{\text{KN}}(\gamma/\gamma_{\text{KN}})}, \end{aligned} \quad (45)$$

where we make the relativistic $\beta \sim 1$ approximation in both cases. The above expressions also invoke the *nominal efficient-cooling* Lorentz factor

$$\gamma_{\text{cool}} \equiv \frac{3m_e c^2}{4\sigma_T U_{\text{ph}} L}, \quad (46)$$

defined by $t_{\text{cool},T}(\gamma_{\text{cool}}) \equiv L/c$. Thus, in the Thomson regime, $\gamma_{\text{cool}}$ corresponds to the minimum Lorentz factor for a particle to cool on time-scales shorter than the system light/Alfvén-crossing time, $L/c$. One may express $\gamma_{\text{cool}}$ in terms of the radiative compactness, $\ell \equiv U_{\text{ph}}\sigma_T L/m_e c^2$, as $\gamma_{\text{cool}} = 3/4\ell$. Unlike a real particle Lorentz factor, the formal parameter $\gamma_{\text{cool}}$ can be less than 1. This corresponds to the highly compact case, $\ell > 1$, and signals that all particles cool to non-relativistic energies in less than $L/c$.

In addition to $\gamma_{\text{cool}}$ – and ignoring Klein–Nishina effects for the moment – another radiative Lorentz factor scale may be defined by equating the radiative cooling time, $t_{\text{cool},T}(\gamma)$, to the linear acceleration time for particles experiencing the reconnection electric field $E_{\text{rec}} = 0.1B_0$ near reconnection X-points,

$$t_{\text{X}}(\gamma) = \frac{\gamma m_e c^2}{0.1 e c B_0} = \frac{10\gamma\rho_0}{c}. \quad (47)$$

Putting $t_X(\gamma_{\text{rad},T}) \equiv t_{\text{cool},T}(\gamma_{\text{rad},T})$ yields the *nominal radiatively-limited Lorentz factor* (cf. Nalewajko 2016; Uzdensky 2016; Werner et al. 2019; Mehlhaff et al. 2020, 2021; Sironi & Beloborodov 2020; Sridhar et al. 2021, 2023),

$$\gamma_{\text{rad},T} \equiv \sqrt{\frac{0.3e B_0}{4\sigma_T U_{\text{ph}}}}. \quad (48)$$

A particle with Lorentz factor $\gamma > \gamma_{\text{rad},T}$ experiences a Thomson radiative drag force stronger than the acceleration force from the reconnection electric field. Exceeding this energy is impossible in the Thomson regime in the absence of other, faster (i.e. operating on time-scales $<t_X$) acceleration mechanisms (Mehlhaff et al. 2021).

The Lorentz factors $\gamma_{\text{cool}}$ and $\gamma_{\text{rad},T}$ are not independent; they are both set by the same underlying radiative parameter $U_{\text{ph}}$. Hence, $\gamma_{\text{rad},T}$ is fixed by $\gamma_{\text{cool}}$ and $\gamma_{\text{max}}$ (Mehlhaff et al. 2021):

$$\gamma_{\text{rad},T}^2 = \gamma_{\text{cool}}\gamma_{\text{max}}. \quad (49)$$

Equation (49) reflects the fact that, since the X-point acceleration time – which defines $\gamma_{\text{rad},T}$ – is faster than the system light crossing time – which defines $\gamma_{\text{cool}}$ – one generally has $\gamma_{\text{cool}} < \gamma_{\text{rad},T} < \gamma_{\text{max}}$. The exception is the non-radiative regime, where $\gamma_{\text{cool}} > \gamma_{\text{rad},T} > \gamma_{\text{max}}$ and all particles cool on time-scales exceeding $L/c$.

The Lorentz factors $\gamma_{\text{cool}}$, $\gamma_{\text{rad},T}$, and $\gamma_{\text{max}}$ characterize the topology of the time-scales $L/c$, $t_X(\gamma)$, and $t_{\text{cool},T}(\gamma)$ by defining all three pairwise intersection points in $\gamma$-space: $t_{\text{cool},T}(\gamma_{\text{cool}}) \equiv L/c$, $t_X(\gamma_{\text{rad},T}) \equiv t_{\text{cool},T}(\gamma_{\text{rad},T})$, and [using equations (38) and (47)] $t_X(\gamma_{\text{max}}) \equiv L/c$. We depict this graphically in Fig. 2. Such a topological view helps illustrate a noteworthy point of conceptual consistency. Namely, because $L/c$, $t_X(\gamma)$, and $t_{\text{cool},T}(\gamma)$ depend, respectively, on $L$, $B_0$, and $U_{\text{ph}}$, the following hold: $\gamma_{\text{cool}}$, at the intersection of $t_{\text{cool},T}(\gamma)$ and $L/c$, depends on $U_{\text{ph}}$ and $L$ but not on $B_0$; $\gamma_{\text{rad},T}$, at the intersection of $t_{\text{cool},T}(\gamma)$ and $t_X(\gamma)$, depends on $U_{\text{ph}}$ and $B_0$ but not on $L$; and $\gamma_{\text{max}}$, at the intersection of $t_X(\gamma)$ and $L/c$, depends on $B_0$ and $L$ but not on $U_{\text{ph}}$ (it is the only non-radiative parameter of the three).

Let us now summarize the situation for Thomson IC cooling. In this regime, according to equation (41), cooling depends only on $U_{\text{ph}}$. The seed photon energy, $\epsilon_{\text{ph}}$, completely drops out of the dynamics (though it still influences the observed photon energies). The parameter $U_{\text{ph}}$ can then be recast in terms of the particle energy scales $\gamma_{\text{cool}}$ and $\gamma_{\text{rad},T}$, which characterize intersection points of important time-scales in the reconnection problem (Fig. 2).

To illuminate the influence of the seed photon energy, $\epsilon_{\text{ph}}$, in the general IC case (i.e. including the Klein–Nishina regime), a similar procedure can be employed as for the Thomson limit. Here, the relevant energy scale in terms of which $\epsilon_{\text{ph}}$ is recast, $\gamma_{\text{KN}}$, has already been introduced; it is the Lorentz factor above which particles lose their energies in discrete photon quanta. We then only need to repeat the comparison of time-scales as done above, but now replacing the Thomson-regime $t_{\text{cool},T}(\gamma)$ with the more general $t_{\text{cool},\text{IC}}(\gamma)$. This is done graphically in Fig. 2, where, in addition to the time-scales $L/c$, $t_X(\gamma)$, and $t_{\text{cool},T}(\gamma)$ discussed above, we also plot $t_{\text{cool},\text{IC}}(\gamma)$ for three representative values of $\gamma_{\text{KN}}$: $\gamma_{\text{KN}}^{(\text{lo})}$, $\gamma_{\text{KN}}^{(\text{med})}$, and $\gamma_{\text{KN}}^{(\text{hi})}$. In the figure, we define the new auxiliary energy scale $\gamma_{\text{rad},\text{IC}}$ as the generalization of $\gamma_{\text{rad},T}$ to include Klein–Nishina effects: $t_{\text{cool},\text{IC}}(\gamma_{\text{rad},\text{IC}}) \equiv t_X(\gamma_{\text{rad},\text{IC}})$. As is the case for $\gamma_{\text{rad},T}$, exceeding $\gamma_{\text{rad},\text{IC}}$ is impossible barring acceleration channels faster than direct acceleration near reconnection X-points (Mehlhaff et al. 2021).

The three sets of curves in Fig. 2 corresponding to the different $\gamma_{\text{KN}}$ values illustrate the three time-scale topologies that may be realized as $\gamma_{\text{KN}}$ is brought in from infinity (equivalently, as the seed photon energy, $\epsilon_{\text{ph}}$, is increased from zero). First, when $\gamma_{\text{KN}} > \gamma_{\text{rad},T}$ (Fig. 2, green $\gamma_{\text{KN}}^{(\text{hi})}$ curves) the generalized radiative Lorentz factor, $\gamma_{\text{rad},\text{IC}}$, remains approximately equal to $\gamma_{\text{rad},T}$. Thus, radiative losses inhibit particles from accessing energies $\ge \gamma_{\text{KN}}$ where they would experience Klein–Nishina effects. This regime therefore essentially reduces to the Thomson limit of radiative reconnection.

The next regime occurs once $\gamma_{\text{KN}}$ falls below $\gamma_{\text{rad},T}$ (Fig. 2, purple $\gamma_{\text{KN}}^{(\text{med})}$ curves). Then, Klein–Nishina effects [entering through
are formally equal at $\gamma_{KN} = m_e c^2/4\epsilon_{ph})$. Then, particles pass into the Klein–Nishina regime where they lose energy in discrete quanta, ~$\gamma m_e c^2$. Here, the cross-section [e.g. equations (7) and (20)] needs to be described by QED and gives rise to an average power radiated per particle of

$$P_{\text{IC}}(\gamma) = P_{\text{T}}(\gamma) f_{\text{KN}}(\gamma/\gamma_{\text{KN}}), \quad (42)$$

where $f_{KN}(q)$, with $q = \gamma/\gamma_{KN}$, is the dimensionless function (cf. Jones 1968; Nalewajko, Yuan & Chruścińska 2018; Mehlhaff et al. 2021)

$$f_{\text{KN}}(q) = \frac{9}{q^3} \left[ \left( \frac{q}{2} + 6 + \frac{6}{q} \right) \ln(1+q) - \frac{1}{(1+q)^2} \left( \frac{11}{12}q^3 + 6q^2 + 9q + 4 \right) - 2 + 2 \text{Li}_2(-q) \right] \quad (43)$$

and $\text{Li}_2(q)$ is the dilogarithm. As necessary to recover the Thomson limit, $f_{KN}(q \ll 1) \to 1$. In the opposite, deep Klein–Nishina regime, $f_{KN}(q \gg 1) \simeq (9/2q^2)[\ln(q) - 11/6]$. Equations (41) and (42) define the respective Thomson-limit and general-case IC cooling times,

$$\begin{aligned} t_{\text{cool,T}}(\gamma) &= \frac{\gamma m_e c^2}{P_{\text{T}}(\gamma)} \approx \frac{3m_e c}{4\sigma_{\text{T}} U_{\text{ph}} \gamma} \\ &= \frac{L}{c} \frac{\gamma_{\text{cool}}}{\gamma} \end{aligned} \quad (44)$$

and

$$\begin{aligned} t_{\text{cool,IC}}(\gamma) &= \frac{\gamma m_e c^2}{P_{\text{IC}}(\gamma)} \approx \frac{3m_e c}{4\sigma_{\text{T}} U_{\text{ph}} \gamma f_{\text{KN}}(\gamma/\gamma_{\text{KN}})} \\ &= \frac{L}{c} \frac{\gamma_{\text{cool}}}{\gamma f_{\text{KN}}(\gamma/\gamma_{\text{KN}})}, \end{aligned} \quad (45)$$

where we make the relativistic $\beta \sim 1$ approximation in both cases. The above expressions also invoke the *nominal efficient-cooling* Lorentz factor

$$\gamma_{\text{cool}} \equiv \frac{3m_e c^2}{4\sigma_{\text{T}} U_{\text{ph}} L}, \quad (46)$$

defined by $t_{\text{cool,T}}(\gamma_{\text{cool}}) \equiv L/c$. Thus, in the Thomson regime, $\gamma_{\text{cool}}$ corresponds to the minimum Lorentz factor for a particle to cool on time-scales shorter than the system light/Alfvén-crossing time, $L/c$. One may express $\gamma_{\text{cool}}$ in terms of the radiative compactness, $\ell \equiv U_{\text{ph}}\sigma_{\text{T}}L/m_e c^2$, as $\gamma_{\text{cool}} = 3/4\ell$. Unlike a real particle Lorentz factor, the formal parameter $\gamma_{\text{cool}}$ can be less than 1. This corresponds to the highly compact case, $\ell > 1$, and signals that all particles cool to non-relativistic energies in less than $L/c$.

In addition to $\gamma_{\text{cool}}$ – and ignoring Klein–Nishina effects for the moment – another radiative Lorentz factor scale may be defined by equating the radiative cooling time, $t_{\text{cool,T}}(\gamma)$, to the linear acceleration time for particles experiencing the reconnection electric field $E_{\text{rec}} = 0.1B_0$ near reconnection X-points,

$$t_{\text{X}}(\gamma) = \frac{\gamma m_e c^2}{0.1 e c B_0} = \frac{10\gamma\rho_0}{c}. \quad (47)$$

Putting $t_X(\gamma_{\text{rad,T}}) \equiv t_{\text{cool,T}}(\gamma_{\text{rad,T}})$ yields the *nominal radiatively-limited Lorentz factor* (cf. Nalewajko 2016; Uzdensky 2016; Werner et al. 2019; Mehlhaff et al. 2020, 2021; Sironi & Beloborodov 2020; Sridhar et al. 2021, 2023),

$$\gamma_{\text{rad,T}} \equiv \sqrt{\frac{0.3eB_0}{4\sigma_{\text{T}}U_{\text{ph}}}}. \quad (48)$$

A particle with Lorentz factor $\gamma > \gamma_{\text{rad,T}}$ experiences a Thomson radiative drag force stronger than the acceleration force from the reconnection electric field. Exceeding this energy is impossible in the Thomson regime in the absence of other, faster (i.e. operating on time-scales $<t_X$) acceleration mechanisms (Mehlhaff et al. 2021).

The Lorentz factors $\gamma_{\text{cool}}$ and $\gamma_{\text{rad,T}}$ are not independent; they are both set by the same underlying radiative parameter $U_{\text{ph}}$. Hence, $\gamma_{\text{rad,T}}$ is fixed by $\gamma_{\text{cool}}$ and $\gamma_{\text{max}}$ (Mehlhaff et al. 2021):

$$\gamma_{\text{rad,T}}^2 = \gamma_{\text{cool}}\gamma_{\text{max}}. \quad (49)$$

Equation (49) reflects the fact that, since the X-point acceleration time – which defines $\gamma_{\text{rad,T}}$ – is faster than the system light crossing time – which defines $\gamma_{\text{cool}}$ – one generally has $\gamma_{\text{cool}} < \gamma_{\text{rad,T}} < \gamma_{\text{max}}$. The exception is the non-radiative regime, where $\gamma_{\text{cool}} > \gamma_{\text{rad,T}} > \gamma_{\text{max}}$ and all particles cool on time-scales exceeding $L/c$.

The Lorentz factors $\gamma_{\text{cool}}$, $\gamma_{\text{rad,T}}$, and $\gamma_{\text{max}}$ characterize the topology of the time-scales $L/c$, $t_X(\gamma)$, and $t_{\text{cool,T}}(\gamma)$ by defining all three pairwise intersection points in $\gamma$-space: $t_{\text{cool,T}}(\gamma_{\text{cool}}) \equiv L/c$, $t_X(\gamma_{\text{rad,T}}) \equiv t_{\text{cool,T}}(\gamma_{\text{rad,T}})$, and [using equations (38) and (47)] $t_X(\gamma_{\text{max}}) \equiv L/c$. We depict this graphically in Fig. 2. Such a topological view helps illustrate a noteworthy point of conceptual consistency. Namely, because $L/c$, $t_X(\gamma)$, and $t_{\text{cool,T}}(\gamma)$ depend, respectively, on $L$, $B_0$, and $U_{\text{ph}}$, the following hold: $\gamma_{\text{cool}}$, at the intersection of $t_{\text{cool,T}}(\gamma)$ and $L/c$, depends on $U_{\text{ph}}$ and $L$ but not on $B_0$; $\gamma_{\text{rad,T}}$, at the intersection of $t_{\text{cool,T}}(\gamma)$ and $t_X(\gamma)$, depends on $U_{\text{ph}}$ and $B_0$ but not on $L$; and $\gamma_{\text{max}}$, at the intersection of $t_X(\gamma)$ and $L/c$, depends on $B_0$ and $L$ but not on $U_{\text{ph}}$ (it is the only non-radiative parameter of the three).

Let us now summarize the situation for Thomson IC cooling. In this regime, according to equation (41), cooling depends only on $U_{\text{ph}}$. The seed photon energy, $\epsilon_{\text{ph}}$, completely drops out of the dynamics (though it still influences the observed photon energies). The parameter $U_{\text{ph}}$ can then be recast in terms of the particle energy scales $\gamma_{\text{cool}}$ and $\gamma_{\text{rad,T}}$, which characterize intersection points of important time-scales in the reconnection problem (Fig. 2).

To illuminate the influence of the seed photon energy, $\epsilon_{\text{ph}}$, in the general IC case (i.e. including the Klein–Nishina regime), a similar procedure can be employed as for the Thomson limit. Here, the relevant energy scale in terms of which $\epsilon_{\text{ph}}$ is recast, $\gamma_{\text{KN}}$, has already been introduced; it is the Lorentz factor above which particles lose their energies in discrete photon quanta. We then only need to repeat the comparison of time-scales as done above, but now replacing the Thomson-regime $t_{\text{cool,T}}(\gamma)$ with the more general $t_{\text{cool,IC}}(\gamma)$. This is done graphically in Fig. 2, where, in addition to the time-scales $L/c$, $t_X(\gamma)$, and $t_{\text{cool,T}}(\gamma)$ discussed above, we also plot $t_{\text{cool,IC}}(\gamma)$ for three representative values of $\gamma_{\text{KN}}$: $\gamma_{\text{KN}}^{(\text{lo})}$, $\gamma_{\text{KN}}^{(\text{med})}$, and $\gamma_{\text{KN}}^{(\text{hi})}$. In the figure, we define the new auxiliary energy scale $\gamma_{\text{rad,IC}}$ as the generalization of $\gamma_{\text{rad,T}}$ to include Klein–Nishina effects: $t_{\text{cool,IC}}(\gamma_{\text{rad,IC}}) \equiv t_X(\gamma_{\text{rad,IC}})$. As is the case for $\gamma_{\text{rad,T}}$, exceeding $\gamma_{\text{rad,IC}}$ is impossible barring acceleration channels faster than direct acceleration near reconnection X-points (Mehlhaff et al. 2021).

The three sets of curves in Fig. 2 corresponding to the different $\gamma_{\text{KN}}$ values illustrate the three time-scale topologies that may be realized as $\gamma_{\text{KN}}$ is brought in from infinity (equivalently, as the seed photon energy, $\epsilon_{\text{ph}}$, is increased from zero). First, when $\gamma_{\text{KN}} > \gamma_{\text{rad,T}}$ (Fig. 2, green $\gamma_{\text{KN}}^{(\text{hi})}$ curves) the generalized radiative Lorentz factor, $\gamma_{\text{rad,IC}}$, remains approximately equal to $\gamma_{\text{rad,T}}$. Thus, radiative losses inhibit particles from accessing energies $\ge \gamma_{\text{KN}}$ where they would experience Klein–Nishina effects. This regime therefore essentially reduces to the Thomson limit of radiative reconnection.

The next regime occurs once $\gamma_{\text{KN}}$ falls below $\gamma_{\text{rad,T}}$ (Fig. 2, purple $\gamma_{\text{KN}}^{(\text{med})}$ curves). Then, Klein–Nishina effects [entering through
<figure><img src="image_2.png" /><figcaption>**Figure 2.** Important reconnection time-scales plotted as functions of $\gamma$. Intersection points define characteristic Lorentz factors. The isosceles triangle formed by the non-radiative time-scales, $L/c$ (horizontal solid black line) and $t_X(\gamma)$ [equation (47); red line], with the Thomson IC cooling time, $t_{\text{cool},T}(\gamma)$ [equation (44), blue line], yields the condition $\gamma_{\text{rad},T}^2 = \gamma_{\text{cool}}\gamma_{\text{max}}$ [equation (49)]. Klein–Nishina effects enter through finite $\gamma_{\text{KN}}$ [equation (7)] and are captured by replacing $t_{\text{cool},T}(\gamma)$ with the general IC cooling time, $t_{\text{cool},\text{IC}}(\gamma)$ [equation (45), where $t_{\text{cool},T}(\gamma) = \lim_{\gamma_{\text{KN}} \to \infty} t_{\text{cool},\text{IC}}(\gamma)$]. We depict three illustrative values of $\gamma_{\text{KN}}$, $\gamma_{\text{KN}}^{(\text{lo})} < \gamma_{\text{KN}}^{(\text{med})} < \gamma_{\text{KN}}^{(\text{hi})}$, corresponding, respectively, to the orange, purple, and green sets of curves. These exemplify the three main topologies, described in Section 3.2.1, that can be formed by $t_{\text{cool},\text{IC}}(\gamma)$ with the non-radiative time-scales, $L/c$ and $t_X(\gamma)$. The depicted $\gamma_{\text{cool}}$, $\gamma_{\text{rad},T}$, and $\gamma_{\text{max}}$ equal those in the runs of Section 4 with radiative cooling: IC(Th), IC(KN), and IC(KN)+PP. In addition, while $\gamma_{\text{KN}}^{(\text{lo})}$ and $\gamma_{\text{KN}}^{(\text{hi})}$ are chosen arbitrarily for illustrative purposes, $\gamma_{\text{KN}}^{(\text{med})} \simeq 13\gamma_{\text{cool}}$ matches the runs of Section 4 with Klein–Nishina effects: IC(KN) and IC(KN)+PP.</figcaption></figure>

the function $f_{\text{KN}}(\gamma/\gamma_{\text{KN}})$ falling appreciably below unity] create a pronounced departure of $t_{\text{cool},\text{IC}}(\gamma)$ from the Thomson limit. As a result, $\gamma_{\text{rad},\text{IC}}$ rises rapidly (it depends super-exponentially on $\gamma_{\text{rad},T}/\gamma_{\text{KN}} > 1$; Mehlhaff et al. 2021) until it eventually crosses $\gamma_{\text{max}}$. The $\gamma_{\text{KN}}^{(\text{med})}$ curves of Fig. 2 illustrate the particular case where $\gamma_{\text{rad},\text{IC}} = \gamma_{\text{max}}$, which already occurs when $\gamma_{\text{KN}}$ is smaller than $\gamma_{\text{rad},T}$ by just a factor of a few. As $\gamma_{\text{KN}}$ continues to diminish from this point, $\gamma_{\text{rad},\text{IC}}$ comes to exceed $\gamma_{\text{max}}$, meaning that Klein–Nishina effects effectively remove the radiative limit on direct X-point acceleration (though of course the intrinsic and system-size limits, $\gamma_X$ and $\gamma_{\text{max}}$, are still present). Nevertheless, as long as $\gamma_{\text{KN}} \gtrsim \gamma_{\text{cool}}$, there remains a range of particle energies characterized by relatively fast cooling, with $t_{\text{cool},\text{IC}}(\gamma) < L/c$.

Finally, as $\gamma_{\text{KN}}$ approaches $\gamma_{\text{cool}}$ (Fig. 2, orange $\gamma_{\text{KN}}^{(\text{lo})}$ curves), the $t_{\text{cool},\text{IC}}(\gamma)$ curve is lifted above the line $L/c$ for all $\gamma$. Here, Klein–Nishina effects suppress cooling to the point that all particles cool on times longer than $L/c$, an effectively non-radiative regime.

Of these three cases, the one where Klein–Nishina effects are expected to influence the reconnection dynamics is the second; only there can particles access energies $>\gamma_{\text{KN}}$ while maintaining a relatively rapid cooling time. We hence call this the Klein–Nishina radiative regime, characterized by the scale hierarchy

$$ \gamma_{\text{cool}} < \gamma_{\text{KN}} < \gamma_{\text{rad},T} < \gamma_{\text{max}}. \qquad (50) $$

As shown next, this hierarchy is not just important from the standpoint of radiative cooling; it is also where pair production from the emitted photons may play an important role.

Particles with Lorentz factor $\gamma_{\text{KN}}$ scatter photons to typical energies $\epsilon_{\text{scat}} \sim \gamma_{\text{KN}}m_e c^2 = (m_e c^2)^2/4\epsilon_{\text{ph}}$: that is, close to pair-production threshold with the seed-photon background. This means that particles must necessarily be accelerated above $\gamma_{\text{KN}}$ in order for a significant fraction of their radiated energy to be recaptured as fresh pairs. In addition, the fiducial pair-production optical depth is [using the peak cross-section $\sigma_{\gamma\gamma} \simeq \sigma_T/5$; equation (14); see also Mehlhaff et al. (2021)]

$$ \tau_{\gamma\gamma} = \frac{U_{\text{ph}}\sigma_T L}{5\epsilon_{\text{ph}}} = \frac{3\gamma_{\text{KN}}}{5\gamma_{\text{cool}}}, \qquad (51) $$

meaning that placing $\gamma_{\text{KN}}$ between $\gamma_{\text{cool}}$ and $\gamma_{\text{rad},T}$ enables both the emission of above-threshold gamma-rays and their absorption inside the system on time-scales $L/c\tau_{\gamma\gamma} < L/c$.

In summary, the Klein–Nishina radiative scale hierarchy, $\gamma_{\text{cool}} < \gamma_{\text{KN}} < \gamma_{\text{rad},T} < \gamma_{\text{max}}$, triggers three simultaneous and important QED effects:

(i) it permits particles to reach energies $>\gamma_{\text{KN}}$, where their IC cooling transitions to the quantized Klein–Nishina regime;
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Important reconnection time-scales plotted as functions of $\gamma$. Intersection points define characteristic Lorentz factors. The isosceles triangle formed by the non-radiative time-scales, $L/c$ (horizontal solid black line) and $t_X(\gamma)$ [equation (47); red line], with the Thomson IC cooling time, $t_{\text{cool},T}(\gamma)$ [equation (44), blue line], yields the condition $\gamma_{\text{rad},T}^2 = \gamma_{\text{cool}}\gamma_{\text{max}}$ [equation (49)]. Klein–Nishina effects enter through finite $\gamma_{\text{KN}}$ [equation (7)] and are captured by replacing $t_{\text{cool},T}(\gamma)$ with the general IC cooling time, $t_{\text{cool},\text{IC}}(\gamma)$ [equation (45), where $t_{\text{cool},T}(\gamma) = \lim_{\gamma_{\text{KN}} \to \infty} t_{\text{cool},\text{IC}}(\gamma)$]. We depict three illustrative values of $\gamma_{\text{KN}}$, $\gamma_{\text{KN}}^{(\text{lo})} < \gamma_{\text{KN}}^{(\text{med})} < \gamma_{\text{KN}}^{(\text{hi})}$, corresponding, respectively, to the orange, purple, and green sets of curves. These exemplify the three main topologies, described in Section 3.2.1, that can be formed by $t_{\text{cool},\text{IC}}(\gamma)$ with the non-radiative time-scales, $L/c$ and $t_X(\gamma)$. The depicted $\gamma_{\text{cool}}$, $\gamma_{\text{rad},T}$, and $\gamma_{\text{max}}$ equal those in the runs of Section 4 with radiative cooling: IC(Th), IC(KN), and IC(KN)+PP. In addition, while $\gamma_{\text{KN}}^{(\text{lo})}$ and $\gamma_{\text{KN}}^{(\text{hi})}$ are chosen arbitrarily for illustrative purposes, $\gamma_{\text{KN}}^{(\text{med})} \simeq 13\gamma_{\text{cool}}$ matches the runs of Section 4 with Klein–Nishina effects: IC(KN) and IC(KN)+PP.</figcaption></figure>

the function $f_{\text{KN}}(\gamma/\gamma_{\text{KN}})$ falling appreciably below unity] create a pronounced departure of $t_{\text{cool},\text{IC}}(\gamma)$ from the Thomson limit. As a result, $\gamma_{\text{rad},\text{IC}}$ rises rapidly (it depends super-exponentially on $\gamma_{\text{rad},T}/\gamma_{\text{KN}} > 1$; Mehlhaff et al. 2021) until it eventually crosses $\gamma_{\text{max}}$. The $\gamma_{\text{KN}}^{(\text{med})}$ curves of Fig. 2 illustrate the particular case where $\gamma_{\text{rad},\text{IC}} = \gamma_{\text{max}}$, which already occurs when $\gamma_{\text{KN}}$ is smaller than $\gamma_{\text{rad},T}$ by just a factor of a few. As $\gamma_{\text{KN}}$ continues to diminish from this point, $\gamma_{\text{rad},\text{IC}}$ comes to exceed $\gamma_{\text{max}}$, meaning that Klein–Nishina effects effectively remove the radiative limit on direct X-point acceleration (though of course the intrinsic and system-size limits, $\gamma_X$ and $\gamma_{\text{max}}$, are still present). Nevertheless, as long as $\gamma_{\text{KN}} \gtrsim \gamma_{\text{cool}}$, there remains a range of particle energies characterized by relatively fast cooling, with $t_{\text{cool},\text{IC}}(\gamma) < L/c$.

Finally, as $\gamma_{\text{KN}}$ approaches $\gamma_{\text{cool}}$ (Fig. 2, orange $\gamma_{\text{KN}}^{(\text{lo})}$ curves), the $t_{\text{cool},\text{IC}}(\gamma)$ curve is lifted above the line $L/c$ for all $\gamma$. Here, Klein–Nishina effects suppress cooling to the point that all particles cool on times longer than $L/c$, an effectively non-radiative regime.

Of these three cases, the one where Klein–Nishina effects are expected to influence the reconnection dynamics is the second; only there can particles access energies $>\gamma_{\text{KN}}$ while maintaining a relatively rapid cooling time. We hence call this the Klein–Nishina radiative regime, characterized by the scale hierarchy

$$ \gamma_{\text{cool}} < \gamma_{\text{KN}} < \gamma_{\text{rad},T} < \gamma_{\text{max}}. \qquad (50) $$

As shown next, this hierarchy is not just important from the standpoint of radiative cooling; it is also where pair production from the emitted photons may play an important role.

Particles with Lorentz factor $\gamma_{\text{KN}}$ scatter photons to typical energies $\epsilon_{\text{scat}} \sim \gamma_{\text{KN}}m_e c^2 = (m_e c^2)^2/4\epsilon_{\text{ph}}$: that is, close to pair-production threshold with the seed-photon background. This means that particles must necessarily be accelerated above $\gamma_{\text{KN}}$ in order for a significant fraction of their radiated energy to be recaptured as fresh pairs. In addition, the fiducial pair-production optical depth is [using the peak cross-section $\sigma_{\gamma\gamma} \simeq \sigma_T/5$; equation (14); see also Mehlhaff et al. (2021)]

$$ \tau_{\gamma\gamma} = \frac{U_{\text{ph}}\sigma_T L}{5\epsilon_{\text{ph}}} = \frac{3\gamma_{\text{KN}}}{5\gamma_{\text{cool}}}, \qquad (51) $$

meaning that placing $\gamma_{\text{KN}}$ between $\gamma_{\text{cool}}$ and $\gamma_{\text{rad},T}$ enables both the emission of above-threshold gamma-rays and their absorption inside the system on time-scales $L/c\tau_{\gamma\gamma} < L/c$.

In summary, the Klein–Nishina radiative scale hierarchy, $\gamma_{\text{cool}} < \gamma_{\text{KN}} < \gamma_{\text{rad},T} < \gamma_{\text{max}}$, triggers three simultaneous and important QED effects:

(i) it permits particles to reach energies $>\gamma_{\text{KN}}$, where their IC cooling transitions to the quantized Klein–Nishina regime;
(ii) it permits particles to reach energies >γ<sub>KN</sub>, where their IC-emitted gamma-rays exceed pair-production threshold with the background radiation bath; and

(iii) it guarantees almost all of these gamma-rays to be absorbed inside the system on subdynamical (<L/c) time-scales, allowing the resulting pairs to feed back on reconnection.

This is the target regime of this study, and how we realize it in
numerical simulations is the topic of the next section.

3.2.2 *Selection of the radiative reconnection parameters*

We now discuss our choices of the radiative reconnection parameters, U<sub>ph</sub> and ϵ<sub>ph</sub>, recast, as described in the preceding Section 3.2.1, in terms of the energy scales γ<sub>cool</sub>, γ<sub>rad, T</sub>, and γ<sub>KN</sub>, plus the pair-production optical depth τ<sub>γγ</sub>. We set as a first goal the Klein–Nishina scale hierarchy, γ<sub>cool</sub> < γ<sub>KN</sub> < γ<sub>rad, T</sub> < γ<sub>max</sub>, discussing afterward the necessary placement of the remaining non-radiative energy scale, σ<sub>c,0</sub> (or, equivalently, γ<sub>X</sub> ≡ 4σ<sub>c,0</sub>) within this base ordering. A key issue is that as much space as possible needs to be opened up between each successive energy scale in the Klein–Nishina hierarchy. The reason for this is that a high optical depth, τ<sub>γγ</sub> = 3γ<sub>KN</sub>/5γ<sub>cool</sub>, demands that γ<sub>KN</sub> substantially exceed γ<sub>cool</sub>. At the same time, γ<sub>KN</sub> cannot be as large as γ<sub>rad, T</sub> without, as previously described, leading to a prohibitively small radiatively limited Lorentz factor, γ<sub>rad, IC</sub>. Thus, the gaps from γ<sub>cool</sub> to γ<sub>KN</sub> and, then, from γ<sub>KN</sub> to γ<sub>rad, T</sub> both need to be as wide as resolution requirements permit. In effect, the separation from γ<sub>cool</sub> to γ<sub>rad, T</sub> needs to be maximized, which is equivalent, through equation (49), to maximizing the ratios γ<sub>max</sub>/γ<sub>rad, T</sub> and γ<sub>max</sub>/γ<sub>cool</sub>. Hence, resolving the necessary scale hierarchy demands pushing the outer scales, γ<sub>cool</sub> and γ<sub>max</sub>, as far away from each other as possible.

Let us now examine the implications this has on numerical
cost. We begin by noting that γ<sub>cool</sub> cannot be made arbitrarily
small. Otherwise, the ambient upstream plasma becomes efficiently
radiative, appreciably cooling down over the course of the simulation.
Physically, this renders the initial background plasma temperature,
T<sub>0</sub>, meaningless and, hence, makes interpreting the simulation
results more difficult; it also means, numerically, that the upstream
Debye length, λ<sub>D,0</sub>, quickly becomes unresolved (we only resolve
it marginally to begin with: λ<sub>D,0</sub> = 1.2Δx), leading to spurious
numerical heating. If, to avoid this, one requires the upstream plasma
cooling time to be at least some factor M > 1 longer than L/c, then
equation (44) implies γ<sub>cool</sub> = M⟨γ⟩ = 3Mθ<sub>0</sub>. Here, ⟨γ⟩ = 3θ<sub>0</sub> = 72
is the initial mean upstream Lorentz factor, and t<sub>cool, T</sub> can be used
instead of t<sub>cool, IC</sub> because upstream particles cool in the Thomson
regime. We find empirically that our simulations take ≃ 3L/c for the
reconnected magnetic flux to saturate, and thus we conservatively set
M = 6, giving γ<sub>cool</sub> = 430 = 3.6 × 10<sup>-3</sup>σ<sub>c,0</sub> = γ<sub>max</sub>/2500.

What the preceding paragraph shows, importantly, is that, regardless of the particular choice for *M*, the minimum value of γ<sub>cool</sub> is inevitably tied, through θ<sub>0</sub>, to the Debye-length resolution requirement. The ratio γ<sub>max</sub>/γ<sub>cool</sub> thus becomes a proxy for L/λ<sub>D,0</sub> ≃ L/Δx ≡ N, the number of cells (in the x-direction) across the simulation box. Specifically

$$
\begin{align}
N = L/\Delta x &= 1.2L/\lambda_{\mathrm{D},0} = 1.2L/\sqrt{\theta_0 \sigma_{\mathrm{c},0}} \rho_0 \nonumber \\
&\approx 12\gamma_{\mathrm{max}}/2\theta_0 \sqrt{\sigma_{\mathrm{h},0}} = 18M\gamma_{\mathrm{max}}/\gamma_{\mathrm{cool}}\sqrt{\sigma_{\mathrm{h},0}} = 7680, \tag{52}
\end{align}
$$

where we used the θ₀ ≫ 1 approximation, σ<sub>h,0</sub> ≃ σ<sub>c,0</sub>/4θ₀. Our
available resources limit us to N = 7680, and, hence, for σ<sub>h,0</sub> =

1250, to γ<sub>max</sub>/γ<sub>cool</sub> ≃ 2500. This ratio is, nevertheless, sufficient to
realize a healthy Klein–Nishina energy-scale hierarchy.

With γ<sub>max</sub>/γ<sub>cool</sub> = √γ<sub>max</sub>/γ<sub>rad, T</sub> = √γ<sub>rad, T</sub>/γ<sub>cool</sub> set, we must now decide where to place γ<sub>KN</sub>. We empirically find that a good fiducial choice is to set γ<sub>KN</sub>/γ<sub>rad, T</sub> such that γ<sub>rad, IC</sub> = γ<sub>max</sub>. Raising γ<sub>KN</sub> from here lowers γ<sub>rad, IC</sub>, somewhat inhibiting reconnection-powered NTPA and, hence, the production of high-energy gamma-rays and pairs. On the other hand, lowering γ<sub>KN</sub> limits the pair-production optical depth, τ<sub>γγ</sub> ∝ γ<sub>KN</sub>/γ<sub>cool</sub>, and, simultaneously, the range of energies where particles are efficiently cooled, t<sub>cool, IC</sub>(γ) < L/c (an extreme case of which are the γ<sub>KN</sub><sup>(lo)</sup> curves of Fig. 2). Overall, this slows the pair-production response of the system: high-energy particles take longer to radiate pair-producing photons, and those photons travel farther before being absorbed to create pairs. Choosing, then, t<sub>cool, IC</sub>(γ<sub>rad, IC</sub>) = t<sub>X</sub>(γ<sub>rad, IC</sub>) = L/c to set γ<sub>rad, IC</sub> = γ<sub>max</sub>, gives γ<sub>KN</sub> = 5500 = 0.046σ<sub>c,0</sub> = 0.26γ<sub>rad, T</sub> = 13γ<sub>cool</sub>, and τ<sub>γγ</sub> = 3γ<sub>KN</sub>/5γ<sub>cool</sub> = 7.7.

Although our fiducial γ<sub>KN</sub> is finely tuned, this parameter is
much more flexible in real astrophysical systems where numerical
requirements do not limit the scale separation. With much larger
values of γ<sub>max</sub>/γ<sub>cool</sub>, astrophysical reconnection may have γ<sub>KN</sub>
different from our numerical sweet spot while still preserving the
important features: copious particle acceleration above γ<sub>KN</sub>, efficient
cooling of high-energy particles, and short mean-free paths of emitted
gamma-rays. A much more detailed discussion of this point for two
concrete astrophysical systems – FSRQ jets and black hole accretion
disc coronae – is presented by Mehlhaff et al. (2021).

Having set γ<sub>KN</sub>, γ<sub>cool</sub>, γ<sub>max</sub>, and τ<sub>γγ</sub>, we are at last equipped to
justify our choice of σ<sub>c,0</sub>. Even though σ<sub>c,0</sub> is technically a non-
radiative reconnection parameter (that we have already discussed
in Section 3.1), it has unique consequences in the presence of
Klein–Nishina radiative physics. Namely, σ<sub>c,0</sub> needs to be at least
a factor of several higher than γ<sub>KN</sub> in order for the intrinsic X-
point particle acceleration limit, γ<sub>X</sub> = 4σ<sub>c,0</sub>, to lie deeply in the
regime where reconnection-energized particles emit pair-producing
photons. As discussed by Mehlhaff et al. (2021), since the pair-
production cross-section, σ<sub>γγ</sub> [equation (14)], peaks at gamma-ray
energies ≃ 3.6(m<sub>e</sub>c<sup>2</sup>)<sup>2</sup>/ϵ<sub>ph</sub> ≃ 14γ<sub>KN</sub>m<sub>e</sub>c<sup>2</sup>, and because a particle
in the Klein–Nishina limit tends to donate about half its energy
to a Comptonized photon, particles emitting photons at peak pair-
production cross-section have typical Lorentz factors γ ≃ 30γ<sub>KN</sub>.
To permit reconnection-energized particles to easily pass this limit,
we have set σ<sub>c,0</sub> = 22γ<sub>KN</sub> (γ<sub>X</sub> = 88γ<sub>KN</sub>). This is why our
simulations feature such a high σ<sub>c,0</sub> and, consequently, a relatively
low γ<sub>max</sub>/σ<sub>c,0</sub> = 0.1L/σ<sub>c,0</sub>ρ<sub>0</sub> = 9.1 given their numerical size. The
issue is that σ<sub>c,0</sub> needs to be on the upper end of our energy-scale
hierarchy, which is already strapped to maximize γ<sub>max</sub>/γ<sub>cool</sub>, and so
σ<sub>c,0</sub> ends up being somewhat close to γ<sub>max</sub>.

A summary of all our radiative and non-radiative simulation
parameters is given in Table 1. In addition, we present a graph-
ical description of how the various reconnection energy scales
relate to key physical and numerical quantities in Fig. 3. The
figure also illustrates how changing one quantity in this high-
dimensional parameter space affects the others. For example, a
simple system-size scan requires added care in the presence of all
of these radiative effects, because changing L (e.g. by changing N),
displaces not only γmax higher, but also γcool lower, reducing the
cooling time of the upstream particles with respect to L/c. Thus,
to conduct such a scan, one would need to take care to set θ0
low enough such that, even at the end of the scan (highest L), the
upstream particles are still sufficiently cold that they do not radiate
appreciably.
(ii) it permits particles to reach energies >γ<sub>KN</sub>, where their IC-emitted gamma-rays exceed pair-production threshold with the background radiation bath; and

(iii) it guarantees almost all of these gamma-rays to be absorbed inside the system on subdynamical (<L/c) time-scales, allowing the resulting pairs to feed back on reconnection.

This is the target regime of this study, and how we realize it in
numerical simulations is the topic of the next section.

3.2.2 *Selection of the radiative reconnection parameters*

We now discuss our choices of the radiative reconnection parameters, $U_{\text{ph}}$ and $\epsilon_{\text{ph}}$, recast, as described in the preceding Section 3.2.1, in terms of the energy scales $\gamma_{\text{cool}}$, $\gamma_{\text{rad, T}}$, and $\gamma_{\text{KN}}$, plus the pair-production optical depth $\tau_{\gamma\gamma}$. We set as a first goal the Klein-Nishina scale hierarchy, $\gamma_{\text{cool}} < \gamma_{\text{KN}} < \gamma_{\text{rad, T}} < \gamma_{\text{max}}$, discussing afterward the necessary placement of the remaining non-radiative energy scale, $\sigma_{c,0}$ (or, equivalently, $\gamma_X \equiv 4\sigma_{c,0}$) within this base ordering. A key issue is that as much space as possible needs to be opened up between each successive energy scale in the Klein-Nishina hierarchy. The reason for this is that a high optical depth, $\tau_{\gamma\gamma} = 3\gamma_{\text{KN}}/5\gamma_{\text{cool}}$, demands that $\gamma_{\text{KN}}$ substantially exceed $\gamma_{\text{cool}}$. At the same time, $\gamma_{\text{KN}}$ cannot be as large as $\gamma_{\text{rad, T}}$ without, as previously described, leading to a prohibitively small radiatively limited Lorentz factor, $\gamma_{\text{rad, IC}}$. Thus, the gaps from $\gamma_{\text{cool}}$ to $\gamma_{\text{KN}}$ and, then, from $\gamma_{\text{KN}}$ to $\gamma_{\text{rad, T}}$ both need to be as wide as resolution requirements permit. In effect, the separation from $\gamma_{\text{cool}}$ to $\gamma_{\text{rad, T}}$ needs to be maximized, which is equivalent, through equation (49), to maximizing the ratios $\gamma_{\text{max}}/\gamma_{\text{rad, T}}$ and $\gamma_{\text{max}}/\gamma_{\text{cool}}$. Hence, resolving the necessary scale hierarchy demands pushing the outer scales, $\gamma_{\text{cool}}$ and $\gamma_{\text{max}}$, as far away from each other as possible.

Let us now examine the implications this has on numerical
cost. We begin by noting that γ<sub>cool</sub> cannot be made arbitrarily
small. Otherwise, the ambient upstream plasma becomes efficiently
radiative, appreciably cooling down over the course of the simulation.
Physically, this renders the initial background plasma temperature,
T<sub>0</sub>, meaningless and, hence, makes interpreting the simulation
results more difficult; it also means, numerically, that the upstream
Debye length, λ<sub>D,0</sub>, quickly becomes unresolved (we only resolve
it marginally to begin with: λ<sub>D,0</sub> = 1.2Δx), leading to spurious
numerical heating. If, to avoid this, one requires the upstream plasma
cooling time to be at least some factor M > 1 longer than L/c, then
equation (44) implies γ<sub>cool</sub> = M⟨γ⟩ = 3Mθ<sub>0</sub>. Here, ⟨γ⟩ = 3θ<sub>0</sub> = 72
is the initial mean upstream Lorentz factor, and t<sub>cool, T</sub> can be used
instead of t<sub>cool, IC</sub> because upstream particles cool in the Thomson
regime. We find empirically that our simulations take ≃ 3L/c for the
reconnected magnetic flux to saturate, and thus we conservatively set
M = 6, giving γ<sub>cool</sub> = 430 = 3.6 × 10<sup>-3</sup>σ<sub>c,0</sub> = γ<sub>max</sub>/2500.

What the preceding paragraph shows, importantly, is that, re-
gardless of the particular choice for *M*, the minimum value of
γ<sub>cool</sub> is inevitably tied, through θ<sub>0</sub>, to the Debye-length resolution
requirement. The ratio γ<sub>max</sub>/γ<sub>cool</sub> thus becomes a proxy for L/λ<sub>D,0</sub>
≈ L/Δx ≡ N, the number of cells (in the x-direction) across the
simulation box. Specifically

$$
\begin{align}
N = L/\Delta x &= 1.2L/\lambda_{\mathrm{D},0} = 1.2L/\sqrt{\theta_0 \sigma_{\mathrm{c},0} \rho_0} \nonumber \\
&\approx 12\gamma_{\mathrm{max}}/2\theta_0 \sqrt{\sigma_{\mathrm{h},0}} = 18M\gamma_{\mathrm{max}}/\gamma_{\mathrm{cool}}\sqrt{\sigma_{\mathrm{h},0}} = 7680, \tag{52}
\end{align}
$$

where we used the $\theta_0 \gg 1$ approximation, $\sigma_{h,0} \simeq \sigma_{c,0}/4\theta_0$. Our
available resources limit us to $N = 7680$, and, hence, for $\sigma_{h,0} =$

1250, to $\gamma_{\text{max}}/\gamma_{\text{cool}} \simeq 2500$. This ratio is, nevertheless, sufficient to
realize a healthy Klein–Nishina energy-scale hierarchy.

With $\gamma_{\text{max}}/\gamma_{\text{cool}} = \sqrt{\gamma_{\text{max}}/\gamma_{\text{rad,T}}} = \sqrt{\gamma_{\text{rad,T}}/\gamma_{\text{cool}}}$ set, we must now decide where to place $\gamma_{\text{KN}}$. We empirically find that a good fiducial choice is to set $\gamma_{\text{KN}}/\gamma_{\text{rad,T}}$ such that $\gamma_{\text{rad,IC}} = \gamma_{\text{max}}$. Raising $\gamma_{\text{KN}}$ from here lowers $\gamma_{\text{rad,IC}}$, somewhat inhibiting reconnection-powered NTPA and, hence, the production of high-energy gamma-rays and pairs. On the other hand, lowering $\gamma_{\text{KN}}$ limits the pair-production optical depth, $\tau_{\gamma\gamma} \propto \gamma_{\text{KN}}/\gamma_{\text{cool}}$, and, simultaneously, the range of energies where particles are efficiently cooled, $t_{\text{cool,IC}}(\gamma) < L/c$ (an extreme case of which are the $\gamma_{\text{KN}}^{(\text{o})}$ curves of Fig. 2). Overall, this slows the pair-production response of the system: high-energy particles take longer to radiate pair-producing photons, and those photons travel farther before being absorbed to create pairs. Choosing, then, $t_{\text{cool,IC}}(\gamma_{\text{rad,IC}}) = t_X(\gamma_{\text{rad,IC}}) = L/c$ to set $\gamma_{\text{rad,IC}} = \gamma_{\text{max}}$, gives $\gamma_{\text{KN}} = 5500 = 0.046\sigma_{c,0} = 0.26\gamma_{\text{rad,T}} = 13\gamma_{\text{cool}}$, and $\tau_{\gamma\gamma} = 3\gamma_{\text{KN}}/5\gamma_{\text{cool}} = 7.7$.

Although our fiducial $\gamma_{\text{KN}}$ is finely tuned, this parameter is much more flexible in real astrophysical systems where numerical requirements do not limit the scale separation. With much larger values of $\gamma_{\text{max}}/\gamma_{\text{cool}}$, astrophysical reconnection may have $\gamma_{\text{KN}}$ different from our numerical sweet spot while still preserving the important features: copious particle acceleration above $\gamma_{\text{KN}}$, efficient cooling of high-energy particles, and short mean-free paths of emitted gamma-rays. A much more detailed discussion of this point for two concrete astrophysical systems – FSRQ jets and black hole accretion disc coronae – is presented by Mehlhaff et al. (2021).

Having set $\gamma_{\text{KN}}$, $\gamma_{\text{cool}}$, $\gamma_{\text{max}}$, and $\tau_{\gamma\gamma}$, we are at last equipped to
justify our choice of $\sigma_{c,0}$. Even though $\sigma_{c,0}$ is technically a non-
radiative reconnection parameter (that we have already discussed
in Section 3.1), it has unique consequences in the presence of
Klein–Nishina radiative physics. Namely, $\sigma_{c,0}$ needs to be at least
a factor of several higher than $\gamma_{\text{KN}}$ in order for the intrinsic X-
point particle acceleration limit, $\gamma_X = 4\sigma_{c,0}$, to lie deeply in the
regime where reconnection-energized particles emit pair-producing
photons. As discussed by Mehlhaff et al. (2021), since the pair-
production cross-section, $\sigma_{\gamma\gamma}$ [equation (14)], peaks at gamma-ray
energies $\simeq 3.6(m_e c^2)^2 / \epsilon_{\text{ph}} \simeq 14\gamma_{\text{KN}} m_e c^2$, and because a particle
in the Klein–Nishina limit tends to donate about half its energy
to a Comptonized photon, particles emitting photons at peak pair-
production cross-section have typical Lorentz factors $\gamma \simeq 30\gamma_{\text{KN}}$.
To permit reconnection-energized particles to easily pass this limit,
we have set $\sigma_{c,0} = 22\gamma_{\text{KN}}$ ($\gamma_X = 88\gamma_{\text{KN}}$). This is why our
simulations feature such a high $\sigma_{c,0}$ and, consequently, a relatively
low $\gamma_{\text{max}}/\sigma_{c,0} = 0.1L/\sigma_{c,0}\rho_0 = 9.1$ given their numerical size. The
issue is that $\sigma_{c,0}$ needs to be on the upper end of our energy-scale
hierarchy, which is already strapped to maximize $\gamma_{\text{max}}/\gamma_{\text{cool}}$, and so
$\sigma_{c,0}$ ends up being somewhat close to $\gamma_{\text{max}}$.

A summary of all our radiative and non-radiative simulation parameters is given in Table 1. In addition, we present a graphical description of how the various reconnection energy scales relate to key physical and numerical quantities in Fig. 3. The figure also illustrates how changing one quantity in this high-dimensional parameter space affects the others. For example, a simple system-size scan requires added care in the presence of all of these radiative effects, because changing *L* (e.g. by changing *N*), displaces not only *γ*<sub>max</sub> higher, but also *γ*<sub>cool</sub> lower, reducing the cooling time of the upstream particles with respect to *L/c*. Thus, to conduct such a scan, one would need to take care to set *θ*<sub>0</sub> low enough such that, even at the end of the scan (highest *L*), the upstream particles are still sufficiently cold that they do not radiate appreciably.
<figure><img src="image_3.png" /><figcaption>**Figure 3.** Important radiative ($\gamma_{\text{cool}}$ and $\gamma_{\text{rad,T}}$) and non-radiative ($\theta_0$, $\sigma_{c,0}$, and $\gamma_{\text{max}}$) energy scales in reconnection connected with lines indicating how their relative magnitudes depend on: the system size, *L*; the number of simulation cells in the *x*-direction, *N*; the characteristic cooling time, $ct_{\text{cool,T}}(3\theta_0)/L$, of the initial upstream particles (*M* in the text); and the hot magnetization, $\sigma_{h,0}$. The labels $\sigma_{h,0}$ and $ct_{\text{cool,T}}(3\theta_0)/L$ follow only in the relativistically hot case, $\theta_0 \gg 1$. The Klein–Nishina Lorentz factor, $\gamma_{\text{KN}}$, and, consequently, $\tau_{\gamma\gamma} \propto \gamma_{\text{KN}}/\gamma_{\text{cool}}$ [equation (51)], can be freely selected independently of all other diagrammed quantities (since $\gamma_{\text{KN}}$ depends only on $\epsilon_{\text{ph}}$, but not on any of the other parameters, like $U_{\text{ph}}$ and $L$, that the other scales depend on). They are placed to illustrate the scale hierarchy of equation (50), $\gamma_{\text{cool}} < \gamma_{\text{KN}} < \gamma_{\text{rad,T}} < \gamma_{\text{max}}$, which is realized in all the simulations of this work that include Klein–Nishina effects.</figcaption></figure>

# 4 KLEIN–NISHINA IMPACT ON RECONNECTION

Here, we present and compare the four main simulations whose setup is discussed in Section 3. We stress that, all parameters being the same (Table 1), the simulations differ only in their modelled physics. In our non-radiative run, ‘no rad.’, radiative cooling and pair production are completely turned off. In our Thomson-radiative run, ‘IC(Th)’, we apply continuous Thomson radiative cooling, but pair production remains absent. In the runs ‘IC(KN)’ and ‘IC(KN)+PP’, we employ our general IC cooling scheme (including the high-energy Klein–Nishina regime; Sections 2.2.1 and 2.2.3), but only self-consistently calculate pair production (Section 2.2.2) in the IC(KN)+PP run. In this sense, the IC(KN) run is artificial. In it, we simply pretend that all radiation emitted by particles, even the part above pair threshold, is permanently lost from the simulation.

In each of the following subsections, we analyse one aspect of the simulations, starting with those that are more similar among the four and moving to those that are more different. We finish by addressing unique properties of the run IC(KN)+PP that only exist in the context of pair production.

## 4.1 Global spatial evolution

We present the large-scale temporal evolution of the simulations in Fig. 4. For generality, the figure depicts the run, IC(KN)+PP, with general IC cooling and pair production, but the temporal evolution in terms of the spatial plasma number density (left-hand column) is similar in all runs.

At the simulation onset, we apply a 1 per cent perturbation to the in-plane magnetic field. This seeds the tearing instability in the initial current sheet, rapidly disrupting it into a chain of plasmons separated by smaller current sheets, themselves tearing-unstable. The successive tearing of these current sheets down to smaller and smaller scales yields a self-similar hierarchy of plasmons and inter-plasmoid current layers (Uzdensky, Loureiro & Schekochihin 2010), as in, for example, the snapshots of Fig. 4 at $t = 0.6L/c$ and $t = 1.2L/c$. As the plasmons merge with one another, they also create separate miniature reconnection sites with current sheets oriented perpendicular to those of the main plasmoid chain (parallel to the $yz$-plane instead of to the $xz$-plane). These reconnection sites – for

example at $x \simeq 35\sigma_{c,0}\rho_0$ in the $t = 1.2L/c$ snapshot and at $x \simeq 5\sigma_{c,0}\rho_0$ in the $t = 1.9L/c$ snapshot – beget their own plasmoid hierarchies.<fnref n="4" />

In both cases – whether in the mid-plane plasmoid chain or at vertical reconnection regions between colliding plasmons – reconnecting current sheets are prominent sites of particle acceleration. This fact is underscored in the depicted IC(KN)+PP run by the copious emission of gamma-ray (above pair-production threshold) radiation from these regions in the middle column of Fig. 4.

Given the periodic boundary conditions of our setup, the large plasmoid with centre at $x \simeq 25\sigma_{c,0}\rho_0$ in our simulations serves as an exhaust for the plasma processed by reconnection. Reconnection slows down and eventually stalls once about half of the initial magnetic flux in the box is reconnected. This corresponds to a state where the *separatrix* – the topological boundary between the domains of reconnected and unreconnected flux – crosses the mid-plane at an angle of about 45° (crosses itself at 90° angles) at the dominant X-point and opens up around the large exhaust plasmoid, which is then the only plasmoid remaining in the layer.

Finally, even though the runs are quite similar in their global spatial evolution, one unique aspect of the IC(KN)+PP run is the difference in spatial coherence between the original particles (Fig. 4, left-hand column) and those produced *in situ* (Fig. 4, right-hand column). Original particles are bound to magnetic field lines (except at reconnection X-points), and thus the striations in the original-particle number density follow closely the wrapping of magnetic field lines around large plasmons. This forms a tree-ring pattern: density striations along field lines trace the history of magnetic flux accumulation onto each plasmoid – by merging with and absorbing smaller plasmons – in the same way that rings on a tree stump trace the felled tree’s lifecycle. In contrast, for the produced particles, such an effect, while still discernible, is much less pronounced. This owes to the added channel through which produced particles can take up residence in plasmons: they can be born there directly. They are not constrained, like the original particles, to essentially follow the

<footnote n="4">⁴In a real instance of astrophysical reconnection, the outer scale, *L*, would dwarf the plasma microscales (e.g. $\sigma_{c,0} \rho_0$) by many orders of magnitude, leading to a deep self-similar hierarchy both in the main plasmoid chain and in the recursive ones birthed between merging plasmons. However, in our simulations with limited computational resources, we only witness the primary chain and the first secondary plasmoid-merger chains.</footnote>
<figure><img src="image_4.png" /><figcaption>**Figure 3.** Important radiative ($\gamma_{\text{cool}}$ and $\gamma_{\text{rad,T}}$) and non-radiative ($\theta_0$, $\sigma_{c,0}$, and $\gamma_{\text{max}}$) energy scales in reconnection connected with lines indicating how their relative magnitudes depend on: the system size, *L*; the number of simulation cells in the *x*-direction, *N*; the characteristic cooling time, $ct_{\text{cool,T}}(3\theta_0)/L$, of the initial upstream particles (*M* in the text); and the hot magnetization, $\sigma_{h,0}$. The labels $\sigma_{h,0}$ and $ct_{\text{cool,T}}(3\theta_0)/L$ follow only in the relativistically hot case, $\theta_0 \gg 1$. The Klein–Nishina Lorentz factor, $\gamma_{\text{KN}}$, and, consequently, $\tau_{\gamma\gamma} \propto \gamma_{\text{KN}}/\gamma_{\text{cool}}$ [equation (51)], can be freely selected independently of all other diagrammed quantities (since $\gamma_{\text{KN}}$ depends only on $\epsilon_{\text{ph}}$, but not on any of the other parameters, like $U_{\text{ph}}$ and $L$, that the other scales depend on). They are placed to illustrate the scale hierarchy of equation (50), $\gamma_{\text{cool}} < \gamma_{\text{KN}} < \gamma_{\text{rad,T}} < \gamma_{\text{max}}$, which is realized in all the simulations of this work that include Klein–Nishina effects.</figcaption></figure>

# 4 KLEIN–NISHINA IMPACT ON RECONNECTION

Here, we present and compare the four main simulations whose setup is discussed in Section 3. We stress that, all parameters being the same (Table 1), the simulations differ only in their modelled physics. In our non-radiative run, ‘no rad.’, radiative cooling and pair production are completely turned off. In our Thomson-radiative run, ‘IC(Th)’, we apply continuous Thomson radiative cooling, but pair production remains absent. In the runs ‘IC(KN)’ and ‘IC(KN)+PP’, we employ our general IC cooling scheme (including the high-energy Klein–Nishina regime; Sections 2.2.1 and 2.2.3), but only self-consistently calculate pair production (Section 2.2.2) in the IC(KN)+PP run. In this sense, the IC(KN) run is artificial. In it, we simply pretend that all radiation emitted by particles, even the part above pair threshold, is permanently lost from the simulation.

In each of the following subsections, we analyse one aspect of the simulations, starting with those that are more similar among the four and moving to those that are more different. We finish by addressing unique properties of the run IC(KN)+PP that only exist in the context of pair production.

## 4.1 Global spatial evolution

We present the large-scale temporal evolution of the simulations in Fig. 4. For generality, the figure depicts the run, IC(KN)+PP, with general IC cooling and pair production, but the temporal evolution in terms of the spatial plasma number density (left-hand column) is similar in all runs.

At the simulation onset, we apply a 1 per cent perturbation to the in-plane magnetic field. This seeds the tearing instability in the initial current sheet, rapidly disrupting it into a chain of plasmons separated by smaller current sheets, themselves tearing-unstable. The successive tearing of these current sheets down to smaller and smaller scales yields a self-similar hierarchy of plasmons and inter-plasmoid current layers (Uzdensky, Loureiro & Schekochihin 2010), as in, for example, the snapshots of Fig. 4 at $t = 0.6L/c$ and $t = 1.2L/c$. As the plasmons merge with one another, they also create separate miniature reconnection sites with current sheets oriented perpendicular to those of the main plasmoid chain (parallel to the $yz$-plane instead of to the $xz$-plane). These reconnection sites – for

example at $x \simeq 35\sigma_{c,0}\rho_0$ in the $t = 1.2L/c$ snapshot and at $x \simeq 5\sigma_{c,0}\rho_0$ in the $t = 1.9L/c$ snapshot – beget their own plasmoid hierarchies.<fnref n="4" />

In both cases – whether in the mid-plane plasmoid chain or at vertical reconnection regions between colliding plasmons – reconnecting current sheets are prominent sites of particle acceleration. This fact is underscored in the depicted IC(KN)+PP run by the copious emission of gamma-ray (above pair-production threshold) radiation from these regions in the middle column of Fig. 4.

Given the periodic boundary conditions of our setup, the large plasmoid with centre at $x \simeq 25\sigma_{c,0}\rho_0$ in our simulations serves as an exhaust for the plasma processed by reconnection. Reconnection slows down and eventually stalls once about half of the initial magnetic flux in the box is reconnected. This corresponds to a state where the *separatrix* – the topological boundary between the domains of reconnected and unreconnected flux – crosses the mid-plane at an angle of about 45° (crosses itself at 90° angles) at the dominant X-point and opens up around the large exhaust plasmoid, which is then the only plasmoid remaining in the layer.

Finally, even though the runs are quite similar in their global spatial evolution, one unique aspect of the IC(KN)+PP run is the difference in spatial coherence between the original particles (Fig. 4, left-hand column) and those produced *in situ* (Fig. 4, right-hand column). Original particles are bound to magnetic field lines (except at reconnection X-points), and thus the striations in the original-particle number density follow closely the wrapping of magnetic field lines around large plasmons. This forms a tree-ring pattern: density striations along field lines trace the history of magnetic flux accumulation onto each plasmoid – by merging with and absorbing smaller plasmons – in the same way that rings on a tree stump trace the felled tree’s lifecycle. In contrast, for the produced particles, such an effect, while still discernible, is much less pronounced. This owes to the added channel through which produced particles can take up residence in plasmons: they can be born there directly. They are not constrained, like the original particles, to essentially follow the

<footnote n="4">⁴In a real instance of astrophysical reconnection, the outer scale, *L*, would dwarf the plasma microscales (e.g. $\sigma_{c,0} \rho_0$) by many orders of magnitude, leading to a deep self-similar hierarchy both in the main plasmoid chain and in the recursive ones birthed between merging plasmons. However, in our simulations with limited computational resources, we only witness the primary chain and the first secondary plasmoid-merger chains.</footnote>
<figure><img src="image_4.png" /><figcaption>**Figure 4.** Time snapshots of the original pair number density (left-hand column), above-threshold photon energy density (middle column), and accumulated produced pair number density (right-hand column) for the main Klein–Nishina reconnection simulation with pair production, IC(KN)+PP. In each panel, the reconnection separatrix is drawn in white and intersects at the main reconnection X-point. All runs exhibit similar large-scale temporal evolution to that depicted: the initial current sheet quickly tears into a hierarchical chain of merging plasmons (e.g. left-hand column). In addition to this familiar picture, in the IC(KN)+PP run, regions of active particle acceleration are signalled by flashes of gamma-rays (middle column). These photons are eventually absorbed to build up the population of particles born *in situ* (right-hand column).</figcaption></figure>

reconnection of a given field line onto a plasmoid. This pollutes, in the number density maps of the produced particles, what are, on the maps of original-particle density, fairly pristine tree-ring-like density striations around plasmons. These remarks are specific to 2D reconnection; in 3D, out-of-plane variation would likely obscure, even in the original particles, the tree-ring structures observed here.

## 4.2 Reconnection rate

Consistent with their similar global evolution, all four of our principal runs exhibit statistically indistinguishable reconnection rates, $\beta_{rec}$. We measure these reconnection rates in Figs 5 and 6. Fig. 5 shows the reconnected flux, $\Phi(t)$, in each simulation as a function of time, $t$. To determine the characteristic reconnection rate, we consider the *active phase* of each simulation, defined as the period, $t \in [t_{start}, t_{end}]$, during which the middle 70 per cent of the change in reconnected flux (difference between initial and final points on the Fig. 5 curves) occurs. This insulates the measurement from artificially slow values during reconnection onset while also reducing sensitivity to the late-

time slowdown during reconnection saturation (see discussion in Section 4.1). We define the average reconnection rate as $[\Phi(t_{end}) - \Phi(t_{start})]/B_{0c}(t_{end} - t_{start})$. Our measurements are summarized in Fig. 6.

Even in the active phase, however, reconnection does not proceed at a precisely uniform rate. To provide some idea of the temporal variation, we also compute instantaneous reconnection rates by differentiating the curves in Fig. 5 with respect to time. We report the median and one standard-deviation (16th and 84th) percentiles of the resulting reconnection rate distributions in Fig. 6.

Concerning the differences between the simulations, the radiative runs generally reconnect about 20 per cent more magnetic flux than the non-radiative simulation. This effect was previously noted by Werner et al. (2019), who interpreted it as resulting from radiative cooling tending to reduce plasma pressure support inside plasmons, enhancing plasmoid compression (cf. Hakobyan et al. 2019, 2023b; Schoeffler et al. 2019, 2023; Chernoglazov, Hakobyan & Philippov 2023). As a result, more reconnected magnetic flux is needed to achieve the saturated geometry where the separatrix forms 90° angles with itself at the principal X-point (Section 4.1). This effect primarily
<figure><img src="image_5.png" /><figcaption>**Figure 4.** Time snapshots of the original pair number density (left-hand column), above-threshold photon energy density (middle column), and accumulated produced pair number density (right-hand column) for the main Klein–Nishina reconnection simulation with pair production, IC(KN)+PP. In each panel, the reconnection separatrix is drawn in white and intersects at the main reconnection X-point. All runs exhibit similar large-scale temporal evolution to that depicted: the initial current sheet quickly tears into a hierarchical chain of merging plasmons (e.g. left-hand column). In addition to this familiar picture, in the IC(KN)+PP run, regions of active particle acceleration are signalled by flashes of gamma-rays (middle column). These photons are eventually absorbed to build up the population of particles born *in situ* (right-hand column).</figcaption></figure>

reconnection of a given field line onto a plasmoid. This pollutes, in the number density maps of the produced particles, what are, on the maps of original-particle density, fairly pristine tree-ring-like density striations around plasmons. These remarks are specific to 2D reconnection; in 3D, out-of-plane variation would likely obscure, even in the original particles, the tree-ring structures observed here.

## 4.2 Reconnection rate

Consistent with their similar global evolution, all four of our principal runs exhibit statistically indistinguishable reconnection rates, $\beta_{rec}$. We measure these reconnection rates in Figs 5 and 6. Fig. 5 shows the reconnected flux, $\Phi(t)$, in each simulation as a function of time, $t$. To determine the characteristic reconnection rate, we consider the *active phase* of each simulation, defined as the period, $t \in [t_{start}, t_{end}]$, during which the middle 70 per cent of the change in reconnected flux (difference between initial and final points on the Fig. 5 curves) occurs. This insulates the measurement from artificially slow values during reconnection onset while also reducing sensitivity to the late-

time slowdown during reconnection saturation (see discussion in Section 4.1). We define the average reconnection rate as $[\Phi(t_{end}) - \Phi(t_{start})]/B_{0c}(t_{end} - t_{start})$. Our measurements are summarized in Fig. 6.

Even in the active phase, however, reconnection does not proceed at a precisely uniform rate. To provide some idea of the temporal variation, we also compute instantaneous reconnection rates by differentiating the curves in Fig. 5 with respect to time. We report the median and one standard-deviation (16th and 84th) percentiles of the resulting reconnection rate distributions in Fig. 6.

Concerning the differences between the simulations, the radiative runs generally reconnect about 20 per cent more magnetic flux than the non-radiative simulation. This effect was previously noted by Werner et al. (2019), who interpreted it as resulting from radiative cooling tending to reduce plasma pressure support inside plasmons, enhancing plasmoid compression (cf. Hakobyan et al. 2019, 2023b; Schoeffler et al. 2019, 2023; Chernoglazov, Hakobyan & Philippov 2023). As a result, more reconnected magnetic flux is needed to achieve the saturated geometry where the separatrix forms 90° angles with itself at the principal X-point (Section 4.1). This effect primarily