High-energy Radiation and Ion Acceleration in Three-dimensional Relativistic Magnetic Reconnection with Strong Synchrotron Cooling

We present the results of 3D particle-in-cell simulations that explore relativistic magnetic reconnection in pair plasma with strong synchrotron cooling and a small mass fraction of nonradiating ions. Our results demonstrate that the structure of the current sheet is highly sensitive to the dynamic efficiency of radiative cooling. Specifically, stronger cooling leads to more significant compression of the plasma and magnetic field within the plasmoids. We demonstrate that ions can be efficiently accelerated to energies exceeding the plasma magnetization parameter, ≫σ, and form a hard power-law energy distribution, f i ∝ γ −1. This conclusion implies a highly efficient proton acceleration in the magnetospheres of young pulsars. Conversely, the energies of pairs are limited to either σ in the strong cooling regime or the radiation burnoff limit, γ syn, when cooling is weak. We find that the high-energy radiation from pairs above the synchrotron burnoff limit, ε c ≈ 16 MeV, is only efficiently produced in the strong cooling regime, γ syn < σ. In this regime, we find that the spectral cutoff scales as ε cut ≈ ε c (σ/γ syn) and the highest energy photons are beamed along the direction of the upstream magnetic field, consistent with the phenomenological models of gamma-ray emission from young pulsars. Furthermore, our results place constraints on the reconnection-driven models of gamma-ray flares in the Crab Nebula.


INTRODUCTION
Over the past two decades, the physics of magnetic reconnection in the magnetospheres of compact astrophysical objects, such as neutron stars and black holes, has been extensively studied both theoretically and numerically.In particular, the key role of reconnection in the process of magnetic energy dissipation and particle acceleration has been widely recognized and accepted (Lyubarskii 1996;Lyubarsky 2005 Sironi et al. 2016;Werner & Uzdensky 2017).In some systems, the presence of large-scale reconnection sites has been demonstrated by global simulations, such as pulsar (e.g., Philippov & Spitkovsky 2014;Chen & Beloborodov 2014) and black hole magnetospheres (e.g., Parfrey et al. 2019;Bransgrove et al. 2021;Ripperda et al. 2022).In other systems, such as the coronae of X-ray binaries or jets of active galactic nuclei (AGNs), the formation of intermittent current layers associated with either turbulent energy cascade or MHD instabilities (e.g., kink instability, Begelman 1998;Tchekhovskoy & Bromberg 2016;Alves et al. 2018;Davelaar et al. 2020;Ortuño-Macías et al. 2022) has been theorized (Beloborodov 2017) or directly demonstrated from intermediate scale simulations (Servidio et al. 2011;Zhdankin et al. 2013Zhdankin et al. , 2017Zhdankin et al. , 2018;;Comisso & Sironi 2018;Davelaar et al. 2020;Chernoglazov et al. 2021;Sironi et al. 2021).Relativistic magnetic reconnection has thus established itself as one of the most important plasma-physical processes powering energy dissipation and non-thermal emission in various contexts of high-energy astrophysics.
Strong radiative cooling, both due to synchrotron emission and inverse-Compton scatterings, has implications for particle acceleration and the reconnection process itself (see, e.g., Uzdensky 2016).The dynamics of the radiative reconnection has so far been mostly studied for 2D current sheets (Jaroschek & Hoshino 2009;Uzdensky & McKinney 2011;Cerutti et al. 2012b,a;Beloborodov 2017;Werner et al. 2019;Sironi & Beloborodov 2020;Mehlhaff et al. 2020;Sridhar et al. 2021).For the case of strong cooling (either synchrotron or inverse-Compton), radiative losses lead to the decrease of the plasma temperature, resulting in a strong compression of plasmoids (e.g., Schoeffler et al. 2019;Hakobyan et al. 2023b).Radiative cooling also substantially limits non-thermal particle acceleration.However, in the case of strong synchrotron cooling, the initial acceleration ("injection") of particles takes place in or close to the X-points, where the radiative cooling is negligible (Cerutti et al. 2014;Kagan et al. 2016a;Hakobyan et al. 2019), allowing the particles to accelerate above the synchrotron "burn-off" limit and emit γ-ray synchrotron photons (Uzdensky et al. 2011;Cerutti et al. 2012a; Kagan et al. 2016b).
The case of 3D relativistic current sheets is substantially less explored (Zenitani & Hoshino 2008;Guo et al. 2014;Sironi & Spitkovsky 2014;Werner & Uzdensky 2017;Zhang et al. 2021a;Werner & Uzdensky 2021;Guo et al. 2021;Schoeffler et al. 2023).It has been shown that there are competing instabilities and acceleration mechanisms specific to 3D.For example, Zenitani & Hoshino (2008), and Cerutti et al. (2014) demonstrated that in the linear stage, the relativistic driftkink instability grows faster than the plasmoid instability.Later, Sironi & Spitkovsky (2014), and Guo et al. (2014) showed that the plasmoid instability still dominates in the non-linear regime, leading to the plasmoiddominated fast reconnection stage in 3D current sheets and efficient particle acceleration.The 3D plasmoids (or flux tubes) resulting from the plasmoid instability are subject to MHD-kink instability, which ultimately leads to turbulence, the nature and properties of which are largely unexplored (Huang & Bhattacharjee 2016;Guo et al. 2021).The mechanisms for particle acceleration are also substantially different between 2D and 3D: while particle "injection," i.e. acceleration to Lorentz factors comparable to the magnetization parameter, σ, proceeds similarly (e.g., Sironi & Spitkovsky 2014;Werner & Uzdensky 2017;Sironi 2022), the acceleration to high energies in 3D is significantly more efficient.In particular, Dahlin et al. (2015), and Zhang et al. (2021aZhang et al. ( , 2023) ) have shown that in 3D, particles are capable of escaping the flux tubes and can be directly accelerated by the reconnection electric field.These studies were mainly focused on the non-radiative regime, and it remains to be seen how radiative cooling affects 3D current sheets (Cerutti et al. 2014;Schoeffler et al. 2023).
Another important aspect of the problem is the composition of current sheets.While most studies of relativistic reconnection focus on pair plasmas, in many scenarios, such as the magnetospheres of pulsars (e.g., Guépin et al. 2020, see Sec. 6.1) and supermassive black holes, we expect a mixture of ions (primarily protons) to be present.Importantly, ions are not susceptible to radiative cooling and have the potential to be accelerated more efficiently by reconnection.Thus far this question, in the context of radiative relativistic reconnection, has not been addressed systematically.
In this paper, we present a set of large 3D PIC simulations of reconnecting current sheets populated by the electron-positron plasma with a small mixture of ions.In our simulations, pairs are subject to dynamically important synchrotron cooling, while ions are unaffected.The goal of our study is to understand the emerging distribution functions of the accelerated ions and pairs, their acceleration mechanisms, and the properties of high-energy radiation.We begin in Sec. 2 with a description of the numerical setup we employ.We then describe the structure of radiatively-cooled 3D current sheets and show that those with the strong cooling form more compressed flux tubes (Sec.3).In Sec. 4 we present our main findings on particle acceleration.In particular, we show that ions are very efficiently accelerated in current sheets with strong radiative cooling, while the Lorentz factors of pairs are limited by the upstream magnetization parameter by the synchrotron losses.In Sec. 5 we discuss the radiation spectra and the beaming of the emission, for varying strength of the synchrotron cooling.We find that 3D current sheets with efficient cooling are capable of emitting high-energy photons above the synchrotron burn-off limit, which are preferentially beamed along the upstream magnetic field.In contrast, current sheets with weak cooling emit photons with energies up to the burn-off limit in the direction perpendicular to the upstream field direction.We summarize our main findings and discuss their implications for magnetospheres of pulsars and black holes in Sec. 6.

NUMERICAL METHODS
For our simulations, we use the Tristan v2 multispecies radiative PIC code (Hakobyan et al. 2023).

Initial configuration and boundary conditions
The simulation domain is initialized with a single current sheet in Harris equilibrium.The magnetic field profile is given by B = B up x tanh (z/δ cs ), i.e., no guide field is present, where x and z are the coordinates along and transverse to the current sheet, correspondingly, and δ cs is the initial half-thickness of the current sheet.We initialize uniform thermal upstream pair plasma with a mixture of ions of overall number density, n up = n i + n e + n e + = 2n e , and a low temperature, T up /(m e c2 ) = 10 −4 .In all simulations presented in the main text, the ions constitute a fraction f i = 0.01 of the total number of particles1 , and we set the ion-toelectron mass ratio to be m i /m e = 1.This is justified because in the process of reconnection at high pair magnetization, σ = B 2 /(4πn e m e c 2 ) ≳ 10 3 (for the realistic mass ratio), ions are quickly accelerated to relativistic energies, ≳ σ(m e /m i ) ≥ 1, similarly to pairs (e.g., Werner et al. 2018).We do not expect a realistic ion-to-electron mass ratio to lead to different results as long as the plasma inertia is dominated by pairs, These assumptions are validated in the Appendix A, where we vary both the mass ratio (up to m i /m e = 25) and the fraction of ions, f i .The magnetization parameter of the upstream plasma, is fixed at σ = 50.The current sheet is initialized as a hot electron-positron-ion plasma of density n cs = N cs n up cosh −2 (z/δ cs ), and temperature T cs /(m e c 2 ) = 2σ/N cs , where N cs = 3 is the maximum overdensity of the layer.Thus, the thermal plasma pressure inside the current sheet balances the pressure of the upstream magnetic field, B 2 up /8π.The initial thickness of the current layer, δ cs = 7d e , is large enough so the tearing instability is not seeded by the numerical noise.
Our simulation domain has outflowing boundary conditions in the ± x-directions, with a thick absorber that gradually damps E and B fields to zero towards the edges of the box.In the ± ẑ-direction we implement a moving plasma injector, similar to Sironi & Spitkovsky (2014); Hakobyan et al. (2019Hakobyan et al. ( , 2021)), which provides a fresh supply of magnetized upstream plasma to replenish the outflowing population.In the third direction, ± ŷ, we impose periodic boundary conditions.The aspect ratio of the domain in the plane of the current sheet is L x /L y = 1, and the size perpendicular to the current sheet is L z ∼ 0.3...0.5L x , large enough to separate the current sheet and the plasma injector.
To trigger the reconnection, we remove thermal pressure in the central part of the current layer as done in Sironi et al. (2016).This procedure launches transients leading to the onset of reconnection and moving outward from the center at approximately the speed of light.To minimize the effects of boundary conditions on the statistical steady state, the radiative cooling and the outflow boundaries are switched on when the transient approaches the edge of the simulation box, at 0.5 L x /c, where L x is the extent of the domain in the x-direction.In all of our analyses, we exclude the particles from the initial current layer, essentially ignoring the evolution of the sheet in its early transient stage.We consider the current sheet to enter the steady-state regime after about t ≳ 1.5 L x /c from the beginning of the simulation.We run all simulations for 3 additional light-crossing times, with a total duration of 4.5 L x /c.All statistical data is averaged over the final 3 light-crossing times, while the current sheet is in statistical steady-state.
To study the numerical convergence of our results, we explore resolutions of d e /∆x = 1...3, where d e = m e c 2 /4πn up e 2 1/2 is the upstream plasma skin depth, and ∆x is the cell size (∆x = ∆y = ∆z); the comparison of these two resolutions is presented in Appendix B, where we find only marginal differences in the reconnection dynamics and particle accel-eration.The choice of σ and d e determines the characteristic resolution of the Larmor orbits of particles, r up L (γ) = γm e c 2 /(|e|B up ), accelerated up to γ ∼ σ: r up L (σ)/∆x = (d e /∆x) √ σ ≈ 7...21, assuming particles gyrate in the perpendicular magnetic field equal to the upstream value.To better characterize the acceleration mechanisms and minimize the boundary effects, we set our fiducial box sizes to be in the range of L x /r up L (σ) ≈ 150...300.Our timestep duration (the CFL number) is fixed at c∆t/∆x = 0.225 for most of our simulations.The summary of all the simulation parameters is given in Table 1; a discussion on how the synchrotron radiation is implemented (last column) follows further (see Sec. 2.2).In our strongest cooling simulation, 1dx1kCool01, we employ a halved CFL number, c∆t/∆x = 0.1125, compared to our default choice of 0.225, in order to properly resolve the plasma oscillation period in regions of dense flux tubes, highly compressed due to strong radiative cooling (see Sec. 3).The upstream plasma is presented with a total concentration of both electrons and positrons averaging to 1 particle per cell3 .To test the convergence, we run a simulation identical to 1dx1kCool02 but with 8 particles per cell in total and obtain identical results.To mitigate the numerical artifacts from the finite number of particles per cell we employ 8 digital filter passes on the deposited currents with a (1/4, 1/2, 1/4) stencil.
To analyze the properties of the current sheet in steady-state, we need to employ a robust method to distinguish the current layer carrying the plasma energized in reconnection from the fresh plasma upstream.The most convenient way is to separately trace particles from the two upstreams (z ≳ 0 and z ≲ 0 regions), and define the mixing factor, M, as suggested by Zhang et al. (2021a), Here n ↑ is the local number density of particles originating from the z > 0 region, and n is the total number density.This definition implies M = 0 in either of the two upstream regions and M ≈ 1 in the mixed region of the current layer.We use the threshold M = 0.7 to separate the upstream from the reconnecting current sheet, however, our results are insensitive to this particular choice, as the upstream-to-downstream transition is quite sharp in terms of the values of M.

Radiation
To model the synchrotron cooling, the code imposes a synchro-curvature radiation reaction force in the Landau & Lifshitz (1975) form, acting exclusively on the electrons and positrons (ions are unaffected by cooling): (3) Here, β = v/c is normalized 3-velocity, and γ = (1 − β 2 ) −1/2 is the Lorentz factor4 .In magnetospheres of compact objects, the radiative drag becomes dynamically important at very high Lorentz factors, γ ≫ 1, depending on the strength of the magnetic field at the location of the reconnecting current sheet.In order to model a similar physical regime in simulations, we rescale the Thomson cross-section σ T , using a dimensionless parameter, γ syn (Uzdensky 2018;Werner et al. 2019;Hakobyan et al. 2019;Sironi & Beloborodov 2020), defined as follows: Thus, the value of γ syn corresponds to the Lorentz factor of particles for which the synchrotron drag force in the perpendicular background field of B up is equal to the acceleration force |e|E rec (where for magnetic reconnection, we can typically estimate the reconnection-driven electric field E rec as ≈ 0.1 B up ).
The relative importance of radiative cooling for particle acceleration is given by the ratio of γ syn and the upstream magnetization parameter, σ.Since the characteristic energy gain in X-points of the reconnecting sheet corresponds to γ ≳ σ (Sironi & Spitkovsky 2014;Werner et al. 2016), we call the cooling regime "strong" if the particles are strongly affected by cooling before their Lorentz factors reach σ around the X-point, i.e., if γ syn ≲ σ.The limit of γ syn > σ is referred to as the "weak" cooling regime.
The second quantity that needs re-scaling is the energy of radiated synchrotron photons, or, equivalently, Planck's constant, ℏ.It is easy to see, that the particles with γ ≈ γ syn , with γ syn being defined as in (4), moving in the perpendicular magnetic field of strength B up radiate photons with the characteristic energy m e c 2 ≈ 16 MeV, (5) where α F ≈ 1/137 is the fine-structure constant (see also Uzdensky et al. 2011;Cerutti et al. 2012a;Uzdensky & Loureiro 2016).Importantly, this value is insensitive to the strength of the magnetic field and is thus a perfect benchmark to calibrate the photon energies.Thus, in our simulations, a particle with Lorentz factor of γ experiencing a magnetic field of B⊥ , defined in equation (3), by definition radiates a spectrum of photons which peaks near (6)

STRUCTURE OF THE CURRENT SHEET
The dynamics of a current sheet, undergoing a 3D magnetic reconnection, is generally similar to that in 2D, but there are also significant differences specific to 3D that must be considered.The relativistic drift-kink instability is the fastest growing in 3D, and dominates the initial evolution.As the current sheet evolves, the plasmoid instability in the nonlinear phase becomes dominant (e.g., Sironi & Spitkovsky 2014).The dynamical state of the current sheet at late stages is thus dominated by a number of large flux tubes (plasmoids), separated by kinetic-scale X-points.Similar to 2D, the plasmoids are continuously forming and undergoing merging processes.However, contrary to 2D, in 3D they are also being continuously "dissipated" by the MHD-kink instability (Werner & Uzdensky 2021).
To visualize the complex 3D structure, Fig. 1 shows volume renderings of the plasma density in the fullydeveloped non-linear reconnection stage of 3D current sheets from two different simulations.Fig. 1a shows a density snapshot of the simulation with weak cooling (γ syn /σ = 5; 3dx3kCool5), while Fig. 1b corresponds to the strong cooling case (γ syn /σ = 0.2; 3dx3kCool02).The small ripples, visible in the perpendicular-to-field direction (y-z plane), are caused by the relativistic driftkink instability.The imprint of the ideal kink instability is visible as the overall twist and deformation of the large flux tubes in both panels of Fig. 1.The lines in the figure represent the magnetic field lines and highlight the main features of the current sheet, including the Xpoints where field lines reconnect crossing the z = 0 plane, and the plasmoids where magnetic field lines are helically wrapped around high-density structures.The color of the field lines shows the strength of the magnetic field |B|, particularly demonstrating that it falls to zero in the X-points.
To systematically quantify the variations of the current sheet structure due to synchrotron cooling, we evaluate the current sheet width, w (along the z axis), for individual points in the horizontal (x-y) plane using the definition of the mixed region in Sec.2.1 (this includes both the thin current sheet and plasmoids of various sizes).Fig. 2 shows the distribution of widths, P(w), for a large collection of points for four simulations with varying cooling strength (in the order of increasing cooling strength: 3dx3kCoolInf, 3dx3kCool5, 3dx3kCool1, 3dx3kCool02).There are two important conclusions one can draw from this plot.First of all, the peak of the P(w) is slightly smaller than w peak ∼ r up the Lorentz factor close to σ in the upstream magnetic field, B up .The presence of this peak is associated with the X-points, where particles are accelerated to characteristic Lorentz factors of γ ≲ σ (see Sec. 4.1 and insets of Fig. 6).More importantly, current sheets with stronger cooling display a steeper slope and shorter distribution span for P(w), indicating that the plasmoids in the current layer become smaller and more concentrated as cooling strength increases (more detailed analysis of this issue was recently performed by Schoeffler et al. 2023).
Variation of the current sheet structure due to synchrotron cooling is reflected in one of the most important macroscopic parameters of the reconnection -its rate, β rec .The evolution of the rate, measured as the inflow velocity of the plasma, β rec = (E × B) z /|B| 2 , averaged in the thin layer far upstream parallel to the current sheet, for different cooling regimes is shown in Fig. 3.As mentioned before, we focus our analysis in the time range t ≳ 1.5 L x /c, when the simulation reaches a dynamic steady state.We find that the reconnection rate in the weak cooling case (γ syn ≳ σ) is β rec ≈ 0.06...0.1 (cf.Zhang et al. 2021a).Notably, in the simulations with stronger cooling (γ syn ≲ σ), the reconnection rate is systematically larger, with values ranging around β rec ≈ 0.12...0.14.We associate this correlation with the change in the filling fraction of X-points in the current layer, which for the stronger cooled runs is approximately twice as large, compared to the run without cooling as seen from the peak values in Fig. 2.This conclusion is also consistent with the characteristically smaller plasmoid sizes in the strong cooling regime, as discussed above.
Plasmoids accumulate the hot plasma energized in the X-points.In 2D the energetic particles are well magnetized and confined, resulting in plasmoids being relatively coherent and clearly distinguishable from other relativistic outflows present in the sheet.In 3D the structure of plasmoids is richer and less coherent.Nonetheless, in both cases, their equilibrium is governed by an adiabatic balance between the gradient of thermal pressure of the hot plasma, P , and the Lorentz force: j × B = c∇P , with j being the total current density.In general, plasma pressure is anisotropic and needs to be computed self-consistently as a moment of the particle distribution function.To compute it, we separate the bulk component of the particle momentum from the thermal one by defining the fluid (Eckart) frame for species s, moving with a velocity U µ s = N µ s / N ν,s N ν s , where N µ s ≡ f s (u µ /u 0 )d 3 u, and f s and u µ are the distribution and the four-velocity of particle species s.The plasma pressure tensor for species s, P µν s , is then computed as the projection of the stress-energy tensor: s is the projection operator, and η µν is the Minkowski tensor (p µ = m s u µ is the fourmomentum of the particle of species s).We find that the non-diagonal elements of the pressure tensor compared to diagonal elements are typically at the level of 5% inside plasmoids in the weak-to-no cooling regime, with an increase to non-negligible ≈ 20% for the strongest cooling regime, γ syn /σ = 0.2.In addition, inside plasmoids, the diagonal component perpendicular to the magnetic field is suppressed by a factor of ∼ 5 in the case of  , d), where we plot the magnetic field, B 2 /8π, and gas pressure, P , with the corresponding magnetic tension force, (j × B)z, and the gradient of the pressure tensor, c∇iP iz .In the uncooled case, the plasmoids are in equilibrium in terms of magnetic and thermal pressure, with the thermal pressure typically dominating inside the plasmoids.On the other hand, the equilibrium inside the cooled plasmoids is more dynamic, with the magnetic field pressure being larger than the upstream value, and dominating the plasmoid cores.
strong cooling.5Guided by these findings, we consider the gradient of the full pressure tensor to analyze the force balance of flux tubes.
In typical 2D simulations (with or without cooling), the force balance inside plasmoids is achieved by strong compression evident in stronger magnetic fields and higher densities (Sironi et al. 2016).The presence of strong cooling exacerbates this effect6 , by removing the pressure support, which results in extremely compressed structures (see, e.g., Hakobyan et al. 2019).In 3D, the picture is quantitatively different.The forcebalance condition is demonstrated in Fig. 4, where in panels (a) and (c) we show the cross-section of the plasma density in x-z plane for two simulations: one without cooling (3dx3kCoolInf), and the other with strong cooling (3dx3kCool02).We pick intermediatesize plasmoids in both of these simulations (indicated by the dotted vertical lines in left panels) and in panels (b) and (d) plot the magnetic pressure, B 2 /8π, the scalar plasma pressure, defined as P = 1/3 s P i s, i , as well as the z-components of the two forces: the pressure gradient, c∇ i P iz , and the Lorentz force, (j × B) z , as functions of z.In the no-cooling simulation (Fig. 4a, b), large plasmoids are typically weakly magnetized, while in the strong cooling regime (Fig. 4c, d), the strength of the magnetic field inside the plasmoids is several times larger than the upstream value, B up , (compare blue solid lines in Fig. 4b, d).Despite this, even in the strongest cooling regime, the magnetic field compression inside the plasmoids is much smaller than what is observed in 2D (Hakobyan et al. 2021(Hakobyan et al. , 2023b;;Schoeffler et al. 2023).Notably, the strongest B-field component inside the plasmoid (in the strong cooling regime) is the outof-plane one, B y .The generation of this out-of-plane magnetic field component, as well as the dissipation of the in-plane B xz components (evident in the no-cooling regime), are both consequences of the MHD kink instability of 3D flux tubes.This instability dissipates the x-z component of the magnetic field, converting it to the plasma thermal energy (e.g., Werner & Uzdensky 2021), while simultaneously acting as a dynamo by generating an out-of-plane field component, B y , until the flux tube becomes stable to kinking (Werner & Uzdensky 2017;Alves et al. 2018;Davelaar et al. 2020;Zhang et al. 2021b).Importantly, this interplay is inherent specifically to 3D reconnection and is thus completely overlooked in 2D.
The plasma inside flux tubes both in uncooled and strongly cooled simulations is compressed and heated to a rough equipartition with the upstream magnetic field pressure: P ≈ B 2 up /8π.In the latter case, the pressure is mainly delivered by higher plasma density (compare densities in panels (a) and (c) in Fig. 4), since strong synchrotron cooling prevents plasmoids from heating up by effectively radiating away the components of particle momenta perpendicular to the magnetic field.Orange lines in Fig. 4b and d show the two force terms: the gradient of the plasma pressure tensor, and the magnetic force, both plotted against the vertical axis z.The balance between these two forces is well satisfied for the flux tube in the uncooled simulation (Fig. 4b), and is reasonable in the strongly cooled simulation (Fig. 4d).
In this section, we largely neglected the presence of ions (insusceptible to synchrotron cooling) in plasmoids.While their contribution to the total pressure was selfconsistently included in our calculations, their relative contribution was typically small for the timescales considered, e.g., in Fig. 4, due to their smaller number density.However, at later times ions can accelerate to much higher Lorentz factors than the pairs (see Sec. 4.1), at which point their lack in the number density is compensated by the higher energy.When the ion pressure dominates in flux tubes, we expect the structures of plasmoids to be similar to the case of no-cooling (Fig. 4a, b and Appendix A).The applicability of this regime for the current sheets in pulsar magnetospheres is studied in detail in Sec.6.1.

PARTICLE ACCELERATION
During the non-linear stage of reconnection, cold particles from upstream -both pairs and ions -are advected into the current sheet, where they are energized in either the X-points or the relativistic outflows from the X-points (e.g., French et al. 2022).After this initial "injection" stage most of these particles enter into the plasmoids.In 2D, strong compression of the magnetic field inside plasmoids leads to adiabatically slow secondary acceleration of these particles to high energies limited by the system size (Petropoulou & Sironi 2018;Hakobyan et al. 2019).However, in 3D, the finite size of plasmoids in the out-of-plane direction allows the most energetic particles to freely escape from their bounds.Moreover, as described above, compression of the magnetic field inside plasmoids is significantly smaller in 3D.Thus, the same secondary acceleration mechanism as in 2D cannot operate efficiently in 3D.On the other hand, in 3D, another acceleration channel has been observed for uncooled pair plasmas by Zhang et al. (2021aZhang et al. ( , 2023)), For reference, we also show the corresponding Larmor radii, r up L .
with the large-scale reconnection electric field in the upstream, E rec , being its main driver.The dynamics of plasma in X-points is insensitive to the presence of synchrotron cooling, as the magnetic field vanishes in these regions.However, as soon as the energized pairs leave the X-points, they lose their energy shortly after being exposed to a strong perpendicular magnetic field component.As a result, the secondary acceleration channel in the reconnection upstream is strongly suppressed for pairs.Ions, on the other hand, are unaffected by cooling, and can thus tap the upstream electric field and get accelerated to potentially large Lorentz factors, ≫ σ.This can clearly be seen in Fig. 5, where we plot the energy spectra of both pairs and ions.For the simulations with strong-to-marginal cooling (γ syn /σ = 0.2...1), the distribution of pairs (solid lines) typically spans up to ∼ σ, with the slope being close to γ −2 at highest energies.Ions (dashed lines), on the other hand, regardless of the cooling strength, extend to Lorentz factors significantly above σ, with the cutoff energy itself increasing in time.In the uncooled simulation (red lines), the cutoff in the distribution of pairs is exactly the same as for the ions, suggesting that the underlying secondary acceleration mechanisms are identical.Surprisingly, in simulations with stronger cooling, ions form a much harder power law distribution, with the typical slope close to γ −1 for the strongest cooling case, and steepening to γ −1.7 for the uncooled case (this effect is explored in Sec 4.1).
We further quantify the contributions of different acceleration channels in Fig. 6, where we plot the energy spectra of both ions and pairs in the uncooled and strongly cooled simulations.To test the dominant acceleration mechanism up to γ ∼ σ (the injection stage), we choose the particles that experienced the non-ideal electric field, E > B, during their lifetime (dashed lines in Fig. 6).In the inset panels of Fig. 6 we also plot the distribution of energies, ∆γ, gained during the passage of E > B regions.Notably, the particles that have passed through an E > B region dominate at energies ≳ σ (cf.Sironi 2022).However, the net energy gain of particles in X-points in both the uncooled and strongly cooled cases does not exceed a fraction of σ, and thus cannot account for the high-energy power-law tail of either the pairs or the ions (e.g., Werner et al. 2016;Uzdensky 2022).The main energy gain, thus, occurs due to the acceleration by ideal electric fields.For pairs, as we argue further in Sec.4.2, "pick-up" acceleration by outflows from the X-points allows them to gain energies up to σ even in the strong cooling case.

Ion acceleration
In our simulations, ions are modeled as separate particle species with a mass equal (or comparable) to that of pairs.Additionally, ions do not undergo synchrotron cooling.In Fig. 7, we present 2D projections of the ion trajectories obtained from the 3D simulation with strong sycnhrotron cooling, γ rad = 0.2σ, 3dx3kCool02.Each point on the trajectory is color-coded according to the Lorentz factor of the ion, this information is also displayed on the inset axes.The thickness of the line shows whether the ion is upstream (thick) or inside the current sheet (thin).We initiate the tracking by picking two representative ions with energies γ ∼ σ at t = 2.25 L x /c (corresponding to the same moment shown in Fig. 5, when the spectrum shape is in the steady-state), and follow their evolution until the end of the simulation t = 4.5 L x /c.Interestingly, the behavior of the two ions differs significantly.One ion reaches a high Lorentz factor, ∼ 30σ, while the energy of the other ion oscillates around a few σ.Notably, the ion that gains large energy experiences most of its acceleration while moving outside of the current layer (Zhang et al. 2021a(Zhang et al. , 2023)), as evident from the thickness of the line representing the trajectory in Fig. 7.These trajectories correspond to bouncing of accelerated ions between the two converging upstream flows7 .In contrast, the low-energy ion displays a tangled trajectory due to its entrapment inside a plasmoid or frequent scattering on magnetic field inhomogeneities inside the current sheet.This behavior is likely associated with self-generated turbulence (e.g., Guo et al. 2021).
The acceleration of the more energetic ion takes place linearly in time, with its trajectory showing almost no scatterings.Notably, the preferred direction of motion for this ion is in y, while the acceleration rate is almost constant.Both of these facts point towards the large-scale reconnection-driven electric field, E rec ≈ E y , being the dominant acceleration mechanism.The rate of acceleration for these trajectories can be written as .Archetypal trajectories (projected to the x-y and the x-z plane) and the energy evolution of two ions.The color of each dot corresponds to the Lorentz factor of the ion at that particular time, while the size of the dot indicates whether the ion is inside the current sheet (small dots), or upstream (large dots).The ion on the right accelerates to high energies (≫ σ) by the large-scale electric field, Erec, in the region upstream of the layer, while the one on the left gains only a marginal amount of energy (∼ σ).The more energetic ion on the right primarily moves upstream of the current layer along the reconnection electric field (y direction) and accelerates linearly in time.The slower particle, on the other hand, is trapped inside the current sheet and experiences constant scatterings on local inhomogeneities.
where we used E rec = β rec B up (Zhang et al. 2021a).This slope is overplotted with a dashed line in the inset of the left panel of Fig. 7, showing excellent agreement with the acceleration history.Since ions are by construction not susceptible to synchrotron cooling, we find that this acceleration mechanism operates similarly for both the strong cooling and no-cooling regimes.
Despite the ion acceleration mechanism being identical for different cooling regimes, the distribution of ions is much harder in the stronger cooling cases, as highlighted with dashed lines in Fig. 5. Two parameters from Eq. ( 7) that could potentially contribute to different ac-celeration rates are the reconnection rate, β rec ,8 and the beaming of ions towards the direction of the reconnection electric field, β y .To inspect the relative importance of the two effects, in Fig. 8 we plot time-dependent statistical distributions of particle Lorentz factors, γ (panels a1...a4), their velocity component along the reconnection electric field in y, β y (panels b1...b4), and the electric field strength in y (panels c1...c4).The first two columns represent ions from the uncooled simula- tion, while the last two -from the simulation with the strongest cooling, γ syn /σ = 0.2.Since we are primarily interested in the post-injection phase, in the oddnumbered columns we select ions that have reached energies γ in ≈ σ (the core of the distribution), while in the even-numbered ones, we pick ions with γ in ≈ 5σ (the middle of the tail of the distribution).
If we consider the population of particles with energies γ ≈ σ (consider a slice at t = 0 in Fig. 8b1, b3), we see that the motion of ions is either not beamed at all in the y-direction (the case with no cooling), or has only marginal beaming (the case with strong cooling).By this time, the energy gain of these ions, γ in ≈ σ, is dominated by the Fermi acceleration in the current sheet (either via a Fermi kick or the "pick-up" mech-anism).This prelude is consistent with the fact that the ions at these energies are primarily moving in the plane of the reconnecting magnetic field, x-z, i.e., perpendicular to the reconnection electric field, β y ≈ 0, (cf.French et al. 2022).At a later time, a fraction of these particles begin accelerating upstream of the layer along the y direction, i.e., along E rec ; that fraction of "lucky" particles ultimately determines the power-law slope in the distribution.As it is evident from the first column of Fig. 8, in the uncooled run only a small fraction of ions follows the γ ∝ E rec • v linear acceleration path, with the bulk beaming in y-direction remaining marginal.This is due to the fact, that the cross-section of typical plasmoids in the simulation with weak-to-no cooling is significantly larger (see Fig. 2), making them more capable of trapping the freely accelerating ions.On the other hand, in the strong cooling case almost all of the particles that reach γ in ≈ σ are "lucky" to land on the efficiently accelerating trajectories, beamed in the direction of the electric field in the upstream (see the third column in Fig. 8; these trajectories are similar to the one shown in the right column of Fig. 7).As the "lucky" particles become more energized, they are essentially bound to remain on these accelerating trajectories.This can be clearly seen in the second and the fourth panels of Fig. 8, where we show the same tracking statistics, but only for the particles that have already reached γ in ≈ 5σ.For this selection, the fraction of "lucky" particles is substantially larger.This suggests that even in the uncooled case the ion energization mechanism to γ ≫ σ is the same large-scale acceleration in the upstream, albeit a smaller amount of accelerated particles, free-streaming in the upstream, results in a steeper power-law slope in their distribution.
We observe a general tendency that the acceleration rate for individual "lucky" ions is larger for the strongly cooled simulations.This can be clearly seen from the typical strength of the electric field experienced by the ions in Fig. 8c1...c4: for stronger cooling the electric field is systematically larger than in the uncooled case, even if we consider particles beyond γ in ≈ 5σ.We associate this difference with the fact that the global reconnection rate, β rec , is larger for the strong cooling case, as shown in Fig. 3, which in turn means larger E y = E rec ≈ β rec B up .Another notable feature is the evolution of the distribution of electric field values, E y , felt by ions (especially for γ in ≈ σ, panels c1 and c3).At early times, while particles still have relatively low energies, they move primarily in the current sheet, where the electric field is mostly chaotic.When accelerated to large-enough energies, most of these ions are moving primarily in the upstream (cf.Fig. 7), where the elec-tric field is coherent and almost exclusively points in the y-direction.

Pair acceleration
Acceleration mechanisms for pairs in the weak-to-no cooling regimes are identical to that for ions: pairs are first accelerated in the X-points and their outflows by non-ideal and ideal electric fields up to Lorentz factors comparable to σ, after which a fraction of them escape to the upstream where they are linearly accelerated further by the large-scale E rec .However, when the cooling is enabled, even if it is dynamically weak, γ syn ≳ σ, this linear acceleration cannot energize the pairs to arbitrarily high energies.In this case, the typical value of the effective perpendicular magnetic field, B, defined in (3), for particles accelerating upstream is close to B up .This implies that the upstream acceleration in the weak cooling regime is limited by the burnoff limit, γ ≲ γ syn .In Fig. 5, this effect can clearly be seen by comparing the uncooled simulation with the one with γ syn = 5σ.While in the uncooled case, the distributions of both pairs and ions extend to several tens of σ in energy, in the weakly cooled simulation the distribution of pairs essentially cuts off at ≈ γ syn = 5σ, with ions being relatively unaffected.
In the strong cooling regime, the upstream acceleration is essentially prohibited, as the strong synchrotron cooling disallows particles to leave the X-points by crossing magnetic field lines.As shown in Fig. 5, the spectrum of accelerated leptons does not extend beyond ≈ σ.In this regime, two acceleration mechanisms are particularly efficient: the X-point acceleration by the non-ideal electric field, and the "pick-up" mechanism.In Fig. 6b, we show that particle energy gain in the E > B zones is not sufficient to extend the power-law to γ ≈ σ; this implies that pairs with the highest energies are additionally accelerated by the "pick-up" mechanism.
An example of a typical high-energy electron trajectory in the strong cooling regime (3dx2kCool02) is shown in Fig. 9.In panels (a)...(e) we plot the x-z slices of the plasma density, and the magnetic field lines at different timesteps, with the current positions of the particle indicated with colored circles.The instantaneous four-velocities of the particle are indicated with colored arrows.We also over-plot in blue the magnetic field line at the position of the particle, with the width of the line indicating its strength.On the right of Fig. 9 we show the time evolution of the Lorentz factor γ, with the dashed lines indicating the times of the snapshots on the left.In panel (a) the particle is upstream from the layer, moving towards the current sheet with a small Lorentz factor, γ ∼ 1.In panel (b), the particle enters show the density slice and the magnetic field lines in the x-z plane at a particular time of the simulation (corresponding to the right panel).The position and the velocity of the electron are shown with a colored circle and an arrow (color corresponds to its energy at that particular time).The magnetic field line passing through the position of the particle is shown as a blue line, with its thickness corresponding to the strength of the magnetic field.The electron starts cold upstream of the current layer (a) and is accelerated to the Lorentz factor up to ≈ σ in the current sheet (b), after which it is picked up by the outflow from the X-point and gains energy further (c).Ultimately, the electron experiences an intensive slow-down when it collides with the magnetic field inside the plasmoid (d), where it remains throughout the rest of its lifetime (e).
the X-point and is accelerated by the non-ideal electric field, E > B, to γ ≈ 0.5σ.After this rapid acceleration event, the particle moves along with the relativistic outflow from the X-point constantly crossing the upstream/downstream boundary, which is perturbed vertically due to RDKI (see below).As the particle hits the upstream, it experiences subtle energy losses indicated as oscillations in its energy, γ.After that, the particle is picked up by a contracting magnetic field line performing a slingshot-like motion.After that, as shown in panel (c), for an extended period of time particle moves along with a rapidly receding field line while being slowly accelerated to γ ≳ σ.Once the slingshot-like motion of the field line comes to a halt, the particle is injected into the forming plasmoid, and experiences a catastrophic cooling event, being exposed to a large B⊥ (panel (d); this event is studied in more details later).During its interaction with the plasmoid, the particle radiates most of its momentum and is then trapped inside the plasmoid until the end of the simulation, with a low Lorentz factor, γ ∼ few, as shown in panel (e).
We further describe the 3D motion of the same accelerated electron in Fig. 10a...d, where we show 3D volume renderings of the plasma density and the magnetic field lines in the vicinity of the electron, also shown in Fig. 9. Panels on the right display the energization history of the particle, as well as the interesting quantities zoomed at two distinct moments of the particle's evolution (roughly corresponding to panels (a)/(b), and (c)/(d) of the same figure).In the upper half of Fig. 10, we specifically focus on the time when the particle experiences successive scattering events while moving in the RDKI-perturbed current layer,9 after being pre-energized in the X-point (corresponding to Fig. 9b).As indicated in the upper right panels, just before the scattering event, the particle crosses the null line, |B| ≈ 0, and approaches the peak of the RDKI perturbation, where both B⊥ and |B| grow rapidly (time t ≈ 1.8 in units indicated on the plot; field strengths at particle position are shown with blue and orange lines).B⊥ reaches the critical value, ≈ B up γ syn /γ (see below), after which deceleration occurs due to the radiative drag force, f drag (shown with an orange line), with the particle being slowed down until the field drops below critical (time t ≈ 1.9).
The catastrophic cooling event, shown in the lower half of Fig. 10 (panels c,d; see also Fig. 9d),10 is different from the subtle scattering event discussed above.In this case, the particle is accelerated by the "pick-up" mechanism to energy γ ≳ σ, while moving along a slingshotlike trajectory.Before t ≈ 3.5, although the velocity of the particle is perpendicular to the local magnetic field (as indicated in Fig. 10c, and Fig. 9c, and d), the effective perpendicular magnetic field, B⊥ is almost zero, allowing the particle to avoid being radiatively cooled.In particular, we plot the strength of the magnetic field, and the B⊥ , experienced by the electron, as well as the critical magnetic field, Bupγsyn/γ.We also show the contributions of two forces, the acceleration by the electric field, eE • v, and the cooling, f drag • v, to the energy evolution of the particle, γ.Finally, we show the curvature radius of the magnetic field at the location of the particle.At the time t ≈ 3.55, when the particle enters into the dense flux tube, B⊥ increases dramatically (above the critical value, B up (γ syn /γ)), and the particle experience a strong radiative drag force (as shown in the second panel on the lower-right of Fig. 10 with an orange line, f drag ).In this catastrophic cooling episode, the particle loses most of its energy and is further trapped inside the flux tube (also shown in Fig. 9e).
In both the scattering and the catastrophic cooling examples, the full radiation reaction force, given by Eq. 3, consists of both synchrotron and curvature contributions.As shown in the top and the bottom panels on the right in Fig. 10, in both cases the local radius of curvature of the field line11 during the deceleration is comparable to the characteristic Larmor radius of the particle with the Lorentz factor γ ≈ σ, i.e., it is close to the gyro-radius of the particle.Under these conditions, the synchrotron and curvature radiation mechanisms are indistinguishable, with the radius of curvature of field lines being microscopic and set by the local field perturbations inside the current layer, i.e., the curvature radius is comparable to the characteristic sizes of the flux tubes.
After the acceleration phase, pairs experience a highly inhomogeneous magnetic field while propagating inside the current sheet.Synchrotron cooling strength for a given particle with a Lorentz factor γ is controlled by the effective perpendicular magnetic field strength at the position of the particle, B⊥ .Assuming a balance between the energy gain from the electric field of characteristic strength, |e|E rec , and the radiative drag force, |e|E rec ( B⊥ /B up ) 2 (γ/γ syn ) 2 , one can define the critical value for the effective magnetic field: B⊥ ≈ B up (γ syn /γ).As discussed above, pairs can experience short episodes of stronger B⊥ , when they rapidly cool down on the timescales of gyration.This condition thus implies, that in simulations with a strong cooling, the highest-energy pairs (with γ ≈ σ) that sustain their motion for long periods of time experience a significantly weaker effective perpendicular magnetic field, either biasing towards local weak-field regions or outflows from the X-points.
To directly measure the effect of energy-dependent field inhomogeneities, in Fig. 11 we plot the distribution of energies, γ, and effective magnetic field strengths experienced by particles, B⊥ , for pairs in simulations with varying cooling strength.The dashed line indicates the position of the critical magnetic field, B⊥ ≈ B up (γ syn /γ).Evidently, for all the cases particle distributions do not cross the critical threshold, with the highest-energy particles, γ ≈ σ, following trajectories with B⊥ ≲ B up (γ syn /γ) (cf.Hakobyan et al. 2023b).Importantly, the combination of particle energy, γ, and the effective magnetic field it experiences, B⊥ , determines the peak energy of synchrotron emission, ε p ∝ γ 2 B⊥ , and in the context of the highest energy pairs, this determines the cutoff of the overall radiation spectrum.The results shown in Fig. 11 suggest that the naive prediction of the radiation spectrum from the pair distribution (or vice-versa) typically used in one-zone models can be significantly modified by this non-trivial dependence of B⊥ on the particle energy, especially in systems with dynamically strong synchrotron cooling.As we have seen above, in strongly cooled simulations inhomogeneities of the magnetic field in the current layer do not allow pairs to accelerate indefinitely.In particular, we can estimate the characteristic timescale of catastrophic energy losses for the highest-energy pairs as where ω L = |e|B up /(m e c). 12 We directly measure the continuous residence time, ∆t, at high energies, γ ≳ γ thr , for all the pairs in our strongly cooled runs (γ syn = 0.2σ).The distribution of residence times, normalized to r up L (σ)/c, is shown in Fig. 12 for two values of γ thr : 2γ syn , and σ.To ensure our results are independent of the size of the box, we also over-plot the results with varying domain sizes and separation of scales.Evidently, most of the particles spend a very short amount of time, ∆t ≲ t c , at high energies, after which they rapidly cool down, as we have discussed in the example above.Notice, that this timescale is typically much shorter than the light-crossing time of our simulation domain, L x /c, which implies two important points.First of all, this indicates that our simulation domain is large enough so that the pairs are free to experience many accelerationdeceleration events during their lifetime inside the current layer.Despite this, pairs are unable to gain energy indefinitely in successive acceleration events, as each of these episodes is followed by a rapid energy loss.Second, the statistical balance between these two episodes essentially enforces the critical field criterion, B⊥ γ ≈ B up γ syn , allowing us to extrapolate the upper ⊥ γ 2 , and substitute γsyn, and B⊥ .
Figure 13.Normalized flux density of the emitted radiation for different cooling regimes.The additional panel shows the position of the exponential cutoff of the spectra as a function of the cooling strength.Current sheets with stronger cooling produce a more extended emission spectrum, with significant power above the burnoff limit, εc, and a cutoff being in good agreement with the prediction of Eq. 8: εcut ∼ εc(σ/γsyn).
Current sheets with weak cooling, on the other hand, produce almost no emission above the burnoff limit.
bound on pair acceleration to the current layers with realistic scale separations.

RADIATION FROM THE CURRENT SHEET
The primary probe of the underlying acceleration physics from the reconnecting sheets in astrophysical systems is the properties of the outgoing high-energy radiation.So far we have only focused on the acceleration dynamics of ions and pairs under the influence of the synchrotron drag force, in this section we study the properties of the produced emission.
In Sec.4.2, we demonstrated that in the weak-tomoderate cooling regime, γ syn ≳ σ, pair acceleration is limited to γ max ≲ γ syn .We expect the spectrum of the high-energy radiation to cut off at energies comparable to the burnoff limit, ε c , from (5).In contrast, in the strong cooling regime, γ syn ≲ σ, pairs can accelerate to γ max ∼ σ (see, e.g., Fig. 5), and we expect the radiation spectrum to extend above the burnoff limit, as follows from Eq. ( 6).Following the discussion in Sec.2.2, the estimate of the cutoff energy for the strong cooling regime can be obtained from the critical value of the effective perpendicular magnetic field, B⊥ ≈ B up (γ syn /γ): To verify our expectations, in Fig. 13 we show photon spectra, ε 2 dn ε /dε ≡ ε 2 f ε , from simulations with differ- γ syn = 5σ γ syn = 0.1σ In the weak cooling case, γsyn = 5σ, emission above the burnoff limit is significantly suppressed, and most of the high-energy photons are emitted perpendicular to the upstream magnetic field.On the other hand, in the strong cooling regime, most of the high-energy photons above the burnoff limit are emitted along the upstream magnetic field.
ent cooling regimes.Each of these spectra is averaged in the interval of 1.5 < tc/L < 4.5.To construct the spectrum, we assume that every particle emits a continuous synchrotron spectrum (e.g., Rybicki & Lightman 1979): where K 5/3 is the modified Bessel function, and the energy ε p is defined according to Eq. ( 6).All of the spectra rise linearly at low energies, ε 2 f ε ∝ ε, which is expected from the hard power-law of reconnectionaccelerated particles, f ∼ γ −1 .The rise is then followed by a broad peak, and a near-exponential decay at high energies.To find the position of the exponential cutoff, we fit the high-energy parts of the spectra as ε 2 f ε ∝ exp (−ε/ε cut ).The inferred values, ε cut , are shown in Fig. 13 with colored vertical lines, as well as in the inset panel at the top right.Evidently, our results are consistent both with the prediction of Eq. ( 8) for the strong cooling, ε cut ≈ ε c (σ/γ syn ), as well as for the weak cooling, where the cutoff approaches ε cut ≈ ε c .The angular distribution of the produced emission is another milestone toward applying our results to realistic astrophysical systems.Previous analysis of the emission anisotropy was done theoretically (e.g., Uzdensky et al. 2011;Cerutti et al. 2012a), in 2D simulations (Cerutti et al. 2012b(Cerutti et al. ,a, 2013;;Kagan et al. 2016a;Mehlhaff et al. 2020) or in 3D simulations with a limited separation of scales (Cerutti et al. 2014).In Fig. 14 we show the distribution of photon energies, f (θ, ε), with respect to θ, the angle between the photon (emitted along the instantaneous motion of the electron/positron) momentum and the upstream magnetic field, x, for different values of the cooling strength.In the weakest cooling case, γ syn = 5σ (Fig. 14a), the highest energy pairs are accelerated in the upstream, along the reconnection electric field, E rec ≈ E y .This results in highenergy photons being emitting preferentially perpendicular to the upstream magnetic field (Zhang et al. 2021a, and Sec. 4.2).For the strong cooling (Fig. 14c, d), the emission of photons with energies above the burnoff limit, ε ≳ few • ε c , is strongly anisotropic and is beamed towards the direction of the upstream field.The bulk of these photons is emitted by the mechanism described in Sec.4.2 and particularly demonstrated in the lower half of Fig. 10: the particle is accelerated by the "pick-up" mechanism in the X-point outflow, in the direction of the upstream magnetic field (Fig. 10c), and later produces the short radiation burst, after colliding face-on with the plasmoid.
To understand the spatial distribution of emitting regions, in Fig. 15 we plot a volumetric map of the radiated synchrotron power, P sync = e ± F rad • v, with the F rad defined in Eq. ( 3) for simulations with weak (Fig. 15a), and strong cooling (Fig. 15b), (the summation, e ± , is done over all the pairs in the given cell).We separately plot the total radiated power (green) and the power carried by the photons above the burnoff limit, ε > ε c (blue).For visual guidance, we use the same snapshot and viewing angle as in Fig. 1, and also plot the total plasma density.To properly quantify the radiated power, we normalize it by the incoming Poynting flux, β rec cB 2 up /(4π).
Figure 15.Spatial distribution of regions producing the synchrotron radiation shown together with the total plasma density, for the same snapshots as in Fig. 1. Green volume rendering corresponds to the total synchrotron intensity, while the blue one shows the emission intensity above the burnoff limit, εc.Panel (a) corresponds to the weak cooling, while panel (b) is for the strong cooling.The emission in the strongly cooled current sheet is mainly associated with the peaks of the RDKI and the outer shells of the plasmoids (both for the total emission and the emission above the burnoff).The much fainter emission from the weakly cooled current sheet (note the difference in color scales between the two regimes) is significantly more uniform, with the emission above the burnoff being essentially absent.
In the case of weak cooling (Fig. 15a), the density of radiation losses is very weak, i.e., significantly below the Poynting flux inflow (indicated by the colorbar scales for panels (a) and (b)), substantially uniform, and marginally biased towards plasmoids.This is due to the fact that a substantial fraction of the highest energy particles are accelerated by E rec in the upstream of the current sheet, and trapped by plasmoids.On the other hand, for the strong cooling regime (Fig. 15b) a large fraction (∼ 10%) of the incoming Poynting flux is being radiated away, with most of the power going into photons comparable to and above the burnoff limit, ε c .There are two main zones producing the radiation in this case: the outflows from the X-points, and the boundaries of plasmoids.The former is associated with the emission of the X-point-accelerated particles, which interact with the small-scale field inhomogeneities in the outflows of the current sheet (see also Fig. 10a, b).Radiation near the plasmoid boundaries is associated with catastrophic cooling events when the energetic particles collide face-on with the strong magnetic field surrounding the plasmoids (Fig. 10c, d).These particles lose energy very intensively at a very short timescale and are ultimately trapped by plasmoids while having very low Lorentz factors, γ ≪ σ (also visible in Fig. 9e).Because of this, in the strong cooling regime, plasmoid interiors have almost no contribution to the total radiated power, in contrast to the weak cooling regime.
Emission above the burnoff limit (blue regions) is almost negligible in the case of weak cooling (Fig. 15a), which is consistent with our discussion above, and is also highlighted in the emission spectra in Fig. 13.For the strong cooling (Fig. 15b), the power above the burnoff limit (blue) is similar to the total power (green), with the main difference being that the highest-energy emission is suppressed in plasmoid interiors.Thus, in the strong cooling regime, the inner volume of plasmoids only emit soft photons, with the photons above the burnoff being emitted primarily near their exteriors, in the outflows from the X-points, and at the local field inhomogeneities (such as the RDKI ripples, see also Fig. 10).

DISCUSSION
In this paper, we study the 3D dynamics of magnetic reconnection in radiatively cooled pair-dominated plasma with a small mixture of ions.Our study focuses on different regimes of synchrotron cooling, parametrized by the burnoff Lorentz factor, γ syn , which quantifies the dynamical importance of cooling compared to acceleration by the reconnection electric field.In the weak cooling regime, γ syn ≳ σ, particle acceleration is effective both for pairs and ions, both of which form a power-law distribution function with an index f (γ) ∝ ∼ γ −1.7 that extends to energies substantially above σ.Here, the highest energy particles (ions and pairs) are mainly accelerated by the large-scale reconnection electric field upstream of the current sheet.In contrast to that, in the strong cooling regime, γ syn ≲ σ, the acceleration of pairs is severely limited by the radiative losses, and their maximum Lorentz factor reaches only γ ≈ σ. Surprisingly, the acceleration of uncooled ions in this regime becomes more efficient as they form a very hard power-law distribution, with the distribution approaching γ −1 .As we uncovered in Sec.4.1, this hardening of the power law in ion distribution is due to plasmoids being compressed by radiative cooling, which limits their ability to trap the accelerated particles and prevent acceleration in the reconnection electric field.Simultaneously, the strong cooling increases the filling fraction of the X-points in the current layer, effectively enhancing the reconnection rate and, as a result, the strength of the electric field.
The difference in acceleration mechanisms results in a number of important implications for the properties of the high-energy radiation.In the case of weak cooling, the spectral cutoff is close to the synchrotron burnoff limit, ε c ≈ 16 MeV, as the particles accelerated in the upstream have large pitch angles and, thus, their energy cannot substantially exceed γ syn .Conversely, in the case of strong cooling particles are accelerated up to γ ≈ σ ≳ γ syn in X-points and relativistic outflows, which results in significant emission beyond the burnoff limit, with the cutoff being close to ε cut ≈ ε c (σ/γ syn ).In addition, the difference in acceleration mechanisms results in a significant difference in the anisotropy of the high-energy radiation.For strong cooling, the most energetic photons are emitted by particles accelerated by the "pick-up" mechanism, which primarily move in the direction of the upstream magnetic field.Thus, the highest energy radiation is emitted in the direction along the upstream magnetic field.For the weak cooling, the most energetic pairs propagate along the reconnection electric field, and thus in this regime the most energetic photons are emitted along E rec and perpendicular to the upstream magnetic field B up .

Acceleration of ions in pulsar magnetospheres
Current sheets in the magnetospheres of pulsars are predominantly filled with the electron-positron plasma, however, a small mixture of ions is also expected (Arons 2003;Amato & Arons 2006;Fang et al. 2012;Guépin et al. 2020).Ions can be extracted from the atmosphere of the neutron star near the polar cap region and can carry volumetric and separatrix return current (Timokhin 2006;Philippov & Spitkovsky 2018).Depending on the inclination angle between the magnetic and the rotational axes of the pulsar, a significant fraction of these ions can be carried toward the reconnecting current sheet.Due to their large mass, ions are relatively unaffected by the synchrotron cooling, and may thus be accelerated more efficiently.
The assumption of negligible cooling for ions can be justified for the observed population of pulsars.Using Estimates for burnoff limit of protons, γ i syn , (blue lines) and the Lorentz factor of protons accelerated in the full potential drop near the light cylinder, γ i max , (red lines) for the pulsars in the ATNF catalog.Axes correspond to the rotation period of the star, P⋆, and the magnetic field strength near the surface of the star.The shaded region corresponds to the regime, when the synchrotron cooling of protons becomes dynamically important: γ i max > γ i syn .Evidently, for all of the observed pulsars, the radiation drag of ions is negligible in the context of ion acceleration.
the definition of the synchrotron burnoff Lorentz factor, Eq. ( 4), taking the Thomson cross-section of ions, σ i T = σ T Z 4 /A 2 (m e /m p ) 2 , m i = Am p (m p is the mass of a proton), and q i = Z|e|, we find that for the ions γ i syn = AZ −3/2 (m p /m e )γ syn , with γ syn being the burnoff Lorentz factor for pairs (see Sec. 2.2).Since most of the high-energy emission takes place in regions of the current layer close to the light cylinder, we may further employ B LC ≈ B ⋆ (R ⋆ /R LC ) 3 as a proxy for B up , with R ⋆ , and B ⋆ being the radius of the neutron star, and the magnetic field strength near its surface, calculated as B ⋆ = 10 12 G (P ⋆ /1 s)( Ṗ⋆ /10 −15 ), and R LC = cP ⋆ /2π being the light cylinder radius (P ⋆ is the rotation period of the star).We then find which can be compared with the total potential drop in the current layer, q i β rec B LC R LC , normalized to m i c 2 : Then the assumption of negligible cooling for ions holds when γ i syn ≳ γ i max .We compare these two quantities for protons (A = Z = 1) for the population of about 3000 radio pulsars from the ATNF catalog (Manchester et al. 2005; The Australia Telescope National Facility (ATNF) 2023) by overplotting the isocontours on the P ⋆ − B ⋆ diagram in Fig. 16.The gray-shaded region in the parameter space indicates the strong cooling regime for ions, γ i syn < γ i max , and evidently, none of the observed radio pulsar fall into that category.This justifies our assumption made throughout the paper on neglecting the radiation drag for the ions.As seen from Fig. 16, a small population of the most energetic pulsars with small-enough rotation period and large surface magnetic fields (e.g., Crab) are able to accelerate protons to Lorentz factors ≳ 10 6 , i.e., above PeV energies.For heavier nuclei with Z ≈ A > 1, the cooling is slightly more efficient, as γ i syn /γ i max ∝ Z −1/2 , and may thus inhibit their acceleration even for the most energetic pulsars.Acceleration near the light cylinder is also severely radiation-limited for the parameters of fast new-born magnetars (see, e.g., Arons 2003;Fang et al. 2012).In this case, however, if the reconnection continues to be efficient further from the light cylinder, additional acceleration is possible in the wind, where the magnetic field is weaker and the radiative losses are less severe.
The mass density fraction of ions in pulsar current sheets is low, f i ≈ (m i /m e )α/λ ≈ 0.01 (which is the value we employed in our simulations), where λ = n ± /n GJ ∼ 10 4 is the pair multiplicity, and α = n i /n GJ ∼ 0.1 is the ion density relative to the Goldreich-Julian density 13 , n GJ = |Ω • B|/2πce, where Ω is the angular velocity of the pulsar.The energy budget of even this small fraction of ions, n i m i c 2 ⟨γ i ⟩, can become substantial as they accelerate to high energies, which will inevitably lead to steepening in their distribution function.In the parameter regime where synchrotron cooling is strong for pairs (γ syn < σ LC ), and weak for ions (γ i syn ≳ γ i max ), the energy of pairs is limited to γ max ≈ σ LC , while the distribution of ions can extend to the full potential drop, γ i max , Eq. ( 11).Assuming a power-law index of −1 for ions and ⟨γ ± ⟩ ≈ γ syn for pairs (cf.our simulation 3dx3kCool02), we may rewrite the total energy budget for both species as: 13 Low fractions of the extracted ions are due to the fact that the bulk of the return magnetospheric current is carried by produced electron-positron pairs (e.g., Timokhin & Arons 2013;Chen & Beloborodov 2014;Philippov & Spitkovsky 2018).
On the other hand, the magnetization near the light cylinder, σ LC , can be expressed as σ LC = (P ⋆ /4πλ)(|e|B LC /m e c) = 5γ ± max /λ, where γ ± max = β rec |e|B LC R LC /m e c 2 is defined as the Lorentz factor of an electron or a positron, which taps β rec fraction of the total voltage drop.We can thus find E ± /E i = (5/α)(γ syn /σ LC ), where we used the relation m e γ ± max = m i γ i max .For a Crab-like pulsar with γ syn /σ LC ≈ 0.1, we find that pairs dominate in terms of the total energy budget, i.e., E ± ≫ E i , even when the ions are accelerated to the full voltage.This implies that the feedback of ions on the reconnection dynamics is negligible, and the ion spectrum extends to γ i max without significant deviations from the f (γ) ∝ γ −1 predicted in our simulations.In fact, we do not observe the spectrum to steepen even when E ± /E i ≈ 1/3, as happens in the simulation 3dx3kCool02 at late times.
Another channel of ion acceleration in pulsar magnetospheres is the energy gain by the parallel electric field in the discharge region close to the polar cap.The amplitude of the unscreened voltage is limited by pair cascade, which in young pulsars is triggered when the energetic pairs reach Lorentz factors, γ th ∼ 10 6 ...10 7 .The same potential drop will correspond to the maximum Lorentz factor of protons γ i th ∼ 10 3 ...10 4 (smaller by a factor of m i /m e ).These values are much smaller than the ones corresponding to the acceleration in the current sheets of Crab-like pulsars, γ i max ≳ 10 6 , reinforcing the notion that the acceleration in the current sheet beyond the light cylinder is the most efficient channel of ion energization in magnetospheres of young pulsars.

High-energy emission from young pulsars and supermassive black holes
In Sec. 5, we discussed the predictions for the observed high energy spectrum from the accelerated pairs.To summarize, in the strong cooling regime applicable in the context of the magnetospheric current sheets of young pulsars, the flux density of the synchrotron emission rises linearly νF ν ∝ ν (with ν being the photon frequency) up to the energy comparable to the synchrotron burnoff limit, 16 MeV.The broad peak at energies higher than the burnoff is then followed by a steep cutoff at energies close to 16 • (σ/γ syn ) MeV (cf.Hakobyan et al. 2023a).For the young pulsars, such as the Crab, where the values for the gamma-ray cutoff are between one-to-few GeV (Abdo et al. 2013), our prediction yields σ LC /γ syn ∼ 100 (i.e., the cooling regime is strong).For Crab pulsar, where γ syn ≈ 4 • 10 4 close to the light cylinder, this implies that the magnetization parameter near the current layer is around σ LC ≈ 10 6 .Notably, the presence of pairs with energies ≳ few • 10 6 m e c 2 is a strict requirement to explain the pulsed emission above TeV energies in the Crab pulsar (Ansoldi et al. 2016).
The beaming of the high energy emission is another central outcome of our work.Modeling of gamma-ray lightcurves from global simulations of pulsar magnetospheres (Bai & Spitkovsky 2010;Cerutti et al. 2016;Kalapotharakos et al. 2018) discovered that in order to explain the observed lightcurves the emission has to be beamed along the magnetic field lines, in the corotating frame of the pulsar.It has thus been not entirely clear how magnetic reconnection in the outer-magnetospheric current sheet could be responsible for such energization, as the naive expectation implied that particles were beamed in the direction of the reconnecting electric field, perpendicular to the direction of the magnetic field.Our findings outlined in Sec. 5 and in Fig. 14, however, indicate that in the regime of strong synchrotron cooling, γ syn < σ LC , most of the high-energy emission is indeed beamed along the upstream magnetic field, as a result of particles being reaccelerated by the "pick-up" mechanism, while moving along the relativistic outflows from the X-points in the direction of the upstream field (see also Cerutti et al. 2012b for the 2D study).
Magnetospheres of low-luminosity accreting supermassive black holes have also been predicted to host intermittently formed reconnecting current layers with sizes comparable to few-to-ten black-hole gravitational radii (e.g., Ripperda et al. 2020Ripperda et al. , 2022)).In particular, for the black hole at the center of the M87 galaxy, due to abundant pair-loading, the regime of reconnection is predicted to be somewhat similar to the reconnection in Crab pulsar (Ripperda et al. 2022), with σ ∼ 10 7 , and γ syn ∼ 10 6 ...10 7 .In that regard, our simulations are largely consistent with the results of 2D simulations by Hakobyan et al. (2023c), where the authors predict that most of the power is radiated at 10...100 MeV range, with the most energetic pairs, γ ∼ max (σ, γ syn ) ∼ 10 7 , also producing the observed TeV signal (via inverse Compton scattering of soft disk photons).Remarkably, the luminosity ratio between the TeV signal, and the jet power, 0.1%, is close to the luminosity ratio between the TeV and GeV signal from the Crab pulsar, suggesting that the Compton amplification parameters for both of these vastly different systems are somewhat similar.
Finally, in most of the discussion above, we focused on the strong cooling regime, where significant emission is expected to be observed above the burnoff limit.One important instance where it is unclear whether this regime applies is the Crab Nebula, where reconnection has been proposed (see, e.g., Cerutti et al. 2014;Lyutikov et al. 2016) as a possible mechanism for the gamma-ray flares at energies ≳ 160 MeV observed by the Fermi satellite (Tavani et al. 2011;Buehler et al. 2012).Here, the upstream pair plasma is relativistically hot, ⟨γ⟩ ≫ 1, and the characteristic energy gain per particle is given by the hot magnetization parameter, σ h = B 2 /(4πnm e c 2 h), where h ∼ ⟨γ⟩ corresponds to the plasma enthalpy per particle.The strong cooling regime then applies only if σ h ⟨γ⟩ ≳ γ syn ∼ 10 9 , where the estimate for γ syn is made assuming the strength of the magnetic field is close to milligauss.Our findings suggest that a significant flux at energies above ≳ 160 MeV can only be produced when σ h ≳ 10 • (γ syn /⟨γ⟩) ≈ 1000(γ syn /10 9 )/(⟨γ⟩/10 7 ) −1 .It is currently unclear whether such conditions can be realized in the Nebula.As long as the mass density fraction of ions is relatively small, nimi/(neme) ≪ 1 (which is not the case in panel c), the power-law index of their distribution is insensitive to either fi, or the mi/me.The x-axis of panel (d) is normalized by the mass of the particles.

APPENDIX
A. VARYING THE FRACTION OF IONS AND THE MASS RATIO In this appendix, we compare the efficiency of ion acceleration in the current sheet with a varying number density fraction of ions, f i = n i /n up , as well as the mass ratio, m i /m e .For this study, we use simulations with a spatial resolution of d e = ∆x and a box size L x /d e = 1000.In Sec.4.1, we have demonstrated that in simulations with no cooling and mass ratio m i /m e = 1 the evolution of positrons and ions is identical, and, obviously, there is no difference for different fractions of ions in the simulation.The difference is most pronounced when pairs are strongly cooled.In this section, we compare the evolution of the distribution of ions in the strong cooling regime, γ syn /σ = 0.2.
In Fig. 17a, b, c we show a comparison of the distribution of ions measured in the steady state of the reconnection, t = 4L x /c, for varying number density fractions, f i = 1%, 10%, 49%, while keeping m i /m e = 1 (simulation shown in Fig. 17a is identical to the previously discussed simulation 1dx1kCool02).From these results, we see that as long as the number density contribution of ions is small, f i ≪ 1, their distribution remains very hard, dn i /dγ ∝ γ −1 for energies γ ≳ σ.In case when the number of positrons is low, and ions constitute almost half of the particle population (Fig. 17c), we observe a significant deviation, with the resulting ion distribution, dn i /dγ ∝ γ −1.7 , being similar to that observed in the simulations without synchrotron cooling.In this case, the pressure contribution of hot ions to the dynamics of the reconnection layer and plasmoids is no longer negligible.As a result, plasmoids in this case are much larger than in strongly cooled simulations with a lower number of ions, and this leads to more efficient trapping of the accelerated ions and, ultimately, to a steeper power-law index.In all cases, the spectrum of the strongly cooled pairs is not affected, as their acceleration primarily occurs in the X-point region.
Another limitation of our simulations is the reduced mass ratio, m i /m e , which we have chosen to be 1 for most of the discussions throughout the paper.The reason for this choice is to maximize the separation of scales between the Larmor radius of the most energetic ions, r L , and the size of the box, and to minimize the effects of boundary conditions on particle acceleration.To confirm that the reduced mass ratio m i /m e = 1 adequately captures the general properties of the acceleration of ions, we carried out two simulations with the parameters described above, but with larger mass ratios, m i /m e = 5 and m i /m e = 25.To keep the mass fraction of ions constant, ρ i /ρ e = 2%, we also decreased their number density fraction, f i .The results of these simulations are shown in Fig. 17d.Evidently, the powerlaw slope of the distribution is only marginally affected by the mass ratio, with the power-law tail extending to energies ≫ m e σ (where σ takes into account the corresponding mass ratio as shown in Eq. ( 1)).

B. VARYING THE RESOLUTION
To ensure that our results are robust, and do not depend on the resolution, i.e., the number of cells per upstream plasma skin depth, in this section we vary this parameter, d e = 1...3∆x, while keeping m i /m e = 1, and f i = 1%.For a direct comparison, we will keep other parameters fixed in the simulation: the ratio of the size of the simulation box to the skin depth, L x /d e = 1000, and the magnetization parameter, σ = 50.
We compare two particular cases for the strong cooling, γ syn /σ = 0.2 (simulations 3dx3kCool02 and 1dx1kCool02), where the results are expected to be the most sensitive to the resolution, as these simulations demonstrate a strong compression of plasmoids, making the local skin depth potentially difficult to resolve.We choose the efficiency of particle acceleration, i.e., the total distribution function as well as the distribution of energy gains in E > B regions, as the main criteria for comparing different numerical resolutions.The results of this comparison are shown in Fig. 18.Panels (a) and (b) show slices in y = 0 plane of the plasma density at time t = 2.7L x /c.Panel (c) shows the normalized distribution of both pairs and ions, as well as the contributions by particles that encountered E > B during their lifetime.Panel (d) shows the normalized distribution functions of the energy gained by both pairs and ions during their flythrough in the E > B region.Evidently, the cross-sections of the current sheet are similar in both cases (Fig. 18a, b), with the lower-resolution simulation showing less high-density regions.The power-law slopes of the distributions, especially at energies γ ≳ σ, are also similar both for the high and the low-resolution runs.The energy gain in the regions of non-ideal field E > B is identical for both simulations, showing that the dynamics around X-points is equally resolved for both simulations.

Figure 1 .
Figure 1.Volume rendering of the plasma density in the current sheet with weak (a), and strong (b) cooling at t ∼ 2.2Lx/c.Magnetic field lines are also displayed, with the color representing the strength, |B|.From the field lines, one can easily identify the non-ideal regions -the X-points (selected ones are shown by blue arrows)-where the magnetic field vanishes (this corresponds to regions where the field lines are blue).The sizes of plasmoids (selected ones shown by green arrows) are significantly different in these two simulations, with the strong cooling run displaying smaller structures and larger compression (note that the colorbars for the two plots have different scales).

Figure 2 .
Figure2.Distribution of transverse sizes (along z) of the current layer for different cooling regimes.The scale is zoomed in for w < r up L (σ).Stronger cooled current sheets form on average smaller structures due to significant compression.The peak of the distribution, i.e., most of the current sheet, has a thickness below the Larmor radius of particles accelerated up to σ: w ∼ r up L (σ).

Figure 3 .
Figure 3.Time evolution of the reconnection rate in simulations with different cooling regimes.Measured rates are of order βrec ≈ 0.06...0.14, and are systematically larger for the case of strongest cooling, γsyn/σ = 0.2.

Figure 4 .
Figure 4.The cross-section of typical plasmoids for uncooled (a, b) and strongly cooled (c, d) current sheets.The left column (a, c) shows density slices and the magnetic field lines.The region highlighted with a cyan-dotted boundary is shown in 1D cuts in (b, d), where we plot the magnetic field, B 2 /8π, and gas pressure, P , with the corresponding magnetic tension force, (j × B)z, and the gradient of the pressure tensor, c∇iP iz .In the uncooled case, the plasmoids are in equilibrium in terms of magnetic and thermal pressure, with the thermal pressure typically dominating inside the plasmoids.On the other hand, the equilibrium inside the cooled plasmoids is more dynamic, with the magnetic field pressure being larger than the upstream value, and dominating the plasmoid cores.

Figure 5 .
Figure5.Energy distributions of the pairs (solid lines) and the ions (dotted lines) at t = 2.25 Lx/c for different cooling regimes: purple lines correspond to strong cooling γsyn/σ = 0.2, and the red ones for the uncooled runs, γsyn/σ = ∞.Dynamically strong cooling effectively limits the maximum energy of pairs to ∼ σ, but at the same time facilitates the formation of harder distribution for the ions.For reference, we also show the corresponding Larmor radii, r up L .

Figure 6 .
Figure6.Distribution functions of the pairs (blue) and the ions (orange) for both the uncooled (a) and the strongly cooled (b) simulations.Dashed lines show parts of the distributions contributed by the particles having experienced E > B during their lifetime in the sheet.Inset panels show the distribution of energy gains, ∆γ, of particles during their contiguous presence in these regions.The non-ideal electric field in the X-points only accelerates particles to Lorentz factors of a fraction of σ.At the same time, this process is essential to promote the particles (pairs and ions) for further acceleration up to the highest energies: in relativistic outflows along the layer and the large-scale reconnection electric field upstream.As a result, the high energy tail of the spectrum is fully formed by the particles passed through E > B.
Figure7.Archetypal trajectories (projected to the x-y and the x-z plane) and the energy evolution of two ions.The color of each dot corresponds to the Lorentz factor of the ion at that particular time, while the size of the dot indicates whether the ion is inside the current sheet (small dots), or upstream (large dots).The ion on the right accelerates to high energies (≫ σ) by the large-scale electric field, Erec, in the region upstream of the layer, while the one on the left gains only a marginal amount of energy (∼ σ).The more energetic ion on the right primarily moves upstream of the current layer along the reconnection electric field (y direction) and accelerates linearly in time.The slower particle, on the other hand, is trapped inside the current sheet and experiences constant scatterings on local inhomogeneities.

Figure 8 .
Figure8.Time evolution of a sample of ions that reach energies, correspondingly, γin = σ (odd columns), and γin = 5σ (even columns).Results from two simulations are shown: with no cooling (first two columns), and strong cooling current (last two columns).The upper row (a1...a4) shows the evolution of the Lorentz factor of the ions.The middle row (b1...b4) shows the evolution of the y component of their 3-velocities (along the reconnecting electric field Erec).The lower row (c1...c4) shows the value of the reconnection electric field Ey experienced by these ions.In the strong cooling regime, most of the ions successfully accelerate in the reconnection electric field and experience little-to-no scatterings (evident from the narrow spread in their energy evolution, a3-a4).In the uncooled regime, the energy gain is more chaotic as ions constantly scatter inside the current sheet.

Figure 9 .
Figure 9.Time evolution of a typical accelerated electron in the strongly cooled simulation.The right panel shows the energy evolution of the particle with time.Panels (a)...(e)show the density slice and the magnetic field lines in the x-z plane at a particular time of the simulation (corresponding to the right panel).The position and the velocity of the electron are shown with a colored circle and an arrow (color corresponds to its energy at that particular time).The magnetic field line passing through the position of the particle is shown as a blue line, with its thickness corresponding to the strength of the magnetic field.The electron starts cold upstream of the current layer (a) and is accelerated to the Lorentz factor up to ≈ σ in the current sheet (b), after which it is picked up by the outflow from the X-point and gains energy further (c).Ultimately, the electron experiences an intensive slow-down when it collides with the magnetic field inside the plasmoid (d), where it remains throughout the rest of its lifetime (e).

Figure 10 .
Figure10.Close-up view of the acceleration dynamics of the electron shown in Fig.9.Panels (a...d) show 3D snapshots of the density zoomed in at the position of the particle (shown with a black sphere).We also show the vector of its velocity, and the magnetic field lines close to the location of the particle.Line colors represent the amplitude of the magnetic field, with the colorbar being identical to the one in Fig.1.The middle panel on the right shows the energy evolution of the electron, with the moments of the 3D snapshots indicated by vertical lines.Three extra panels at the top and the bottom of the right column show the time evolution of several important quantities, limited in time to the intervals of the scattering (top), and the catastrophic cooling events (bottom).In particular, we plot the strength of the magnetic field, and the B⊥ , experienced by the electron, as well as the critical magnetic field, Bupγsyn/γ.We also show the contributions of two forces, the acceleration by the electric field, eE • v, and the cooling, f drag • v, to the energy evolution of the particle, γ.Finally, we show the curvature radius of the magnetic field at the location of the particle.

Figure 11 .
Figure 11.The effective perpendicular magnetic field, B⊥ , experienced by the pairs of different energies, γ.The 2D histograms are shown for all the cooling regimes, from weak, γsyn/σ = 5, to strong, γsyn/σ = 0.1.The red dashed lines show the estimated critical magnetic field, B⊥ γ = Bupγsyn.Notably, the acceleration of the highest-energy pairs is inhibited by the synchrotron cooling, and, as a result, most of the pairs do not cross the critical line.

Figure 12 .
Figure12.Distribution of the lifetimes of strongly cooled pairs at high energies (γ ≳ 2γsyn and γ ≳ σ) for different sizes of the simulation domain.For both thresholds, the lifetime is close to one gyro-period of the most energetic particles, r up L (σ)/c and, importantly, is independent of the domain size.

Figure 14 .
Figure 14.Angular anisotropy of the radiation at different photon energies for different cooling regimes.The angle is measured with respect to the direction of the upstream magnetic field Bup ≡ Bx x.In the weak cooling case, γsyn = 5σ, emission above the burnoff limit is significantly suppressed, and most of the high-energy photons are emitted perpendicular to the upstream magnetic field.On the other hand, in the strong cooling regime, most of the high-energy photons above the burnoff limit are emitted along the upstream magnetic field.

Figure 17 .
Figure17.Distribution functions of pairs (dashed) and ions (solid) in simulations containing a different fraction of ions, fi = ni/nup = 1%, 10%, 49% (panels a...c), and different mass ratio mi/me = 5, 25 (panel d).As long as the mass density fraction of ions is relatively small, nimi/(neme) ≪ 1 (which is not the case in panel c), the power-law index of their distribution is insensitive to either fi, or the mi/me.The x-axis of panel (d) is normalized by the mass of the particles.

Figure 18 .
Figure 18.Comparison of 2D plasma density slices (a, b) and normalized particle distributions (c, d) for different numerical resolutions de/∆x = 1, 3. Dark and bright lines correspond to ions and pairs, while blue and orange correspond to high-and low-resolution runs.The acceleration mechanisms, the statistics of the energy gains in E > B regions, and thus the resulting distributions are largely insensitive to numerical resolution.The compression of plasmoids, on the other hand, is somewhat smaller for the low-resolution run.

Table 1 .
Summary of simulation parameters.