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

Alexander Chernoglazov1,2 , Hayk Hakobyan3,4 , and Alexander Philippov1,2
1 Department of Physics, University of Maryland, College Park, MD 20742, USA; achernog@umd.edu

2 Institute for Research in Electronics and Applied Physics, University of Maryland, College Park, MD 20742, USA
3 Computational Sciences Department, Princeton Plasma Physics Laboratory (PPPL), Princeton, NJ 08540, USA
4 Physics Department & Columbia Astrophysics Laboratory, Columbia University, New York, NY 10027, USA
Received 2023 May 4; revised 2023 September 26; accepted 2023 September 27; published 2023 December 13

Abstract

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, fi∝ γ−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≈ 16MeV, 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.

Unified Astronomy Thesaurus concepts: Compact radiation sources (289); Plasma astrophysics (1261); High
energy astrophysics (739)

Supporting material: animation

1. 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; Uzdensky et al. 2010;
Sironi & Spitkovsky 2014; Uzdensky & Spitkovsky 2014;
Uzdensky 2016; Uzdensky & Loureiro 2016; 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., Chen & Beloborodov
2014; Philippov & Spitkovsky 2014) and black hole magneto-
spheres (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 orMHD 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. 2013,
2017, 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 nonthermal 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, 2012a; Beloborodov 2017; Werner et al. 2019;
Mehlhaff et al. 2020; Sironi & Beloborodov 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 nonthermal 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 “burnoff” 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; Guo et al. 2021; Werner & Uzdensky 2021; Schoeffler
et al. 2023). It has been shown that there are competing
instabilities and acceleration mechanisms specific to 3D. For

The Astrophysical Journal, 959:122 (20pp), 2023 December 20 https://doi.org/10.3847/1538-4357/acffc6
© 2023. The Author(s). Published by the American Astronomical Society.

Original content from this work may be used under the terms
of the Creative Commons Attribution 4.0 licence. Any further

distribution of this work must maintain attribution to the author(s) and the title
of the work, journal citation and DOI.

1

https://orcid.org/0000-0001-5121-1594
https://orcid.org/0000-0001-5121-1594
https://orcid.org/0000-0001-5121-1594
https://orcid.org/0000-0001-8939-6862
https://orcid.org/0000-0001-8939-6862
https://orcid.org/0000-0001-8939-6862
https://orcid.org/0000-0001-7801-0362
https://orcid.org/0000-0001-7801-0362
https://orcid.org/0000-0001-7801-0362
mailto:achernog@umd.edu
http://astrothesaurus.org/uat/289
http://astrothesaurus.org/uat/1261
http://astrothesaurus.org/uat/739
http://astrothesaurus.org/uat/739
https://doi.org/10.3847/1538-4357/acffc6
https://crossmark.crossref.org/dialog/?doi=10.3847/1538-4357/acffc6&domain=pdf&date_stamp=2023-12-13
https://crossmark.crossref.org/dialog/?doi=10.3847/1538-4357/acffc6&domain=pdf&date_stamp=2023-12-13
http://creativecommons.org/licenses/by/4.0/


example, Zenitani & Hoshino (2008) and Cerutti et al. (2014)
demonstrated that the relativistic drift-kink instability grows
faster than the plasmoid instability in the linear stage. Later,
Sironi & Spitkovsky (2014) and Guo et al. (2014) showed that
the plasmoid instability still dominates in the nonlinear regime,
leading to the plasmoid-dominated 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 ulti-
mately 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.
(2021a, 2023) have shown that in 3D the 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 nonradiative 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 reconnec-
tion focus on pair plasmas, in many scenarios, such as the
magnetospheres of pulsars (e.g., Guépin et al. 2020, see
Section 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 particle-in-cell
(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
Section 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 (Section 3). In Section 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 Section 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 burnoff limit, which are preferentially
beamed along the upstream magnetic field. In contrast, current
sheets with weak cooling emit photons with energies up to the
burnoff 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 Section 6.

2. Numerical Methods

For our simulations, we use the Tristan v2 multi-species
radiative PIC code (Hakobyan et al. 2023).

2.1. 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 xB ztanhup csd= , 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 n n n n2i e e eup = + + =+ , and a low temperature,
Tup/(mec

2)= 10−4. In all of the simulations presented in the
main text, the ions constitute a fraction fi= 0.01 of the total
number of particles,5 and we set the ion-to-electron mass ratio
to be mi/me= 1. This is justified because in the process of
reconnection at high pair magnetization, σ= B2/(4πnemec

2)
103 (for the realistic mass ratio), ions are quickly accelerated to
relativistic energies, σ(me/mi)� 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, nimi/(2neme)= fi · (mi/me) 1.6

These assumptions are validated in Appendix A, where we
vary both the mass ratio (up to mi/me= 25) and the fraction of
ions, fi.
The magnetization parameter of the upstream plasma,

( )
( )

B

m n c f f

B

m n c

4

1

4
, 1

e i i
m

m
e

up
2

up
2

up
2

up
2i

e

s
p p

º
- +

»
⎡
⎣

⎤
⎦

is fixed at σ= 50. The current sheet is initialized as a hot electron-
positron-ion plasma of density ( )n N n zcoshcs cs up

2
csd= - , and

temperature Tcs/(mec
2)= 2σ/Ncs, where Ncs= 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 8up

2 p. The initial thickness of the current layer, δcs= 7de,
is large enough for the tearing instability to not be 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 toward the edges of the box. In
the ẑ -direction, we implement a moving plasma injector,
similar to Sironi & Spitkovsky (2014) and Hakobyan et al.
(2019, 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 Lx/Ly= 1, and the size perpendicular to the current sheet is
Lz∼ 0.3...0.5Lx, which is 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

5 A simulation with no ions, but otherwise identical to 1dx1kCool02,
demonstrates similar overall dynamics of the reconnection process.
6 This situation applies to the current sheets in pulsar magnetospheres and can
be relevant to the pair-dominated coronae of accreting black holes.

2

The Astrophysical Journal, 959:122 (20pp), 2023 December 20 Chernoglazov, Hakobyan, & Philippov



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 Lx/c, where Lx is the extent of the domain in the x-
direction. In all of our analyzes, 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 Lx/c from
the beginning of the simulation. We run all of the simulations
for three additional light-crossing times, with a total duration of
4.5 Lx/c. All statistical data is averaged over the final three
light-crossing times, while the current sheet is in statistical
steady-state.

To study the numerical convergence of our results, we explore
resolutions of de/Δx= 1...3, where ( )d m c n e4e e

2
up

2 1 2p= is
the upstream plasma skin depth, and Δx is the cell size
(Δx=Δy=Δz); a comparison of these two resolutions is
presented in Appendix B, where we find only marginal
differences in the reconnection dynamics and particle
acceleration. The choice of σ and de determines the
characteristic resolution of the Larmor orbits of particles,

( )rL
up g = (∣ ∣ )m c e Be

2
upg , accelerated up to γ∼ σ: ( )r xL

up s D =
( )d x 7 ... 21e sD » , assuming particles gyrate in the perpend-
icular 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 r 150 ... 300x L

up s » . Our timestep duration (the CFL
number) is fixed at cΔt/Δx= 0.225 for most of our simulations.
A summary of all the simulation parameters is given in Table 1; a
discussion on how the synchrotron radiation is implemented (last
column) follows later on (see Section 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, which are highly compressed due
to strong radiative cooling (see Section 3). The upstream plasma
is presented with a total concentration of both electrons and
positrons averaging to 1 particle per cell.7 To test the
convergence, we run a simulation identical to 1dx1kCool02
but with eight particles per cell in total and obtain identical
results. To mitigate the numerical artifacts from the finite
number of particles per cell, we employ eight 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,, as suggested
by Zhang et al. (2021a),

( )
n

n
1 2

1

2
, 2= - -

where n↑ is the local number density of particles originating
from the z> 0 region and n is the total number density. This
definition implies 0= in either of the two upstream regions
and 1» in the mixed region of the current layer. We use
the threshold  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.

2.2. 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):

( ˜ )
( ) ( · )

˜ ( ) ( · ) ( )

F

E B B E E

E B E

B

B

,

,

, 3

R

R

rad 4
2 2

2 2 2

T k b

k b b

b b

g= -

= + ´ ´ +

= + ´ -

s
p ^

^

where β= v/c is normalized 3-velocity, and ( )1 2 1 2g b= - -

is the Lorentz factor.8 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.
To model a similar physical regime in simulations, we rescale
the Thomson cross-section σT, using a dimensionless para-
meter, γsyn (Uzdensky 2018; Hakobyan et al. 2019; Werner
et al. 2019; Sironi & Beloborodov 2020), defined as follows:

∣ ∣ ( )e E B
1

4
. 4Trec up

2
syn
2

p
s gº

Thus, the value of γsyn corresponds to the Lorentz factor of
particles for which the synchrotron drag force in the
perpendicular background field of Bup is equal to the
acceleration force |e|Erec (where for magnetic reconnection,
we can typically estimate the reconnection-driven electric field
Erec as ≈0.1 Bup).
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.

Table 1
Summary of Simulation Parameters

Simulation Lx × Ly × Lz de/Δx cΔt/Δx γsyn/σ

3dx2kCool02 20002 × 1000 3 0.225 0.2
3dx2kCool5 20002 × 1000 3 0.225 5.0
3dx3kCool02 30002 × 1000 3 0.225 0.2
3dx3kCool05 30002 × 1000 3 0.225 0.5
3dx3kCool1 30002 × 1000 3 0.225 1.0
3dx3kCool5 30002 × 1600 3 0.225 5.0
3dx3kCoolInf 30002 × 1600 3 0.225 ∞

1dx1kCool02 10002 × 1000 1 0.225 0.2
1dx2kCool02 20002 × 1000 1 0.225 0.2
1dx1kCool01 10002 × 1000 1 0.1125 0.1

7 I.e., each cell either contains an electron and a positron injected at the same
location to satisfy zero numerical charge or does not contain any simulation
particles.

8 The total radiated power, ˜P Bsync
2 2gµ ^, corresponds to a square of the

electric field in the rest frame of the moving particle, E¢ =
( ) ( )( · )E B E1 2b b bg g b+ ´ - - .

3

The Astrophysical Journal, 959:122 (20pp), 2023 December 20 Chernoglazov, Hakobyan, & Philippov



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 Bup radiate photons with the
characteristic energy

∣ ∣
( )

e B

m c

E

B
m c

3

2

9

4
16 MeV, 5c

e F
e

up
syn
2 rec

up

2e g
a

= = »

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 that peaks near

˜
( )B

B
16 MeV . 6p

up syn

2

e
g
g

= ^
⎜ ⎟
⎛

⎝

⎞

⎠

3. 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, Figure 1 shows volume
renderings of the plasma density in the fully developed nonlinear
reconnection stage of 3D current sheets from two different

simulations. Figure 1(a) shows a density snapshot of the
simulation with weak cooling (γsyn/σ= 5; 3dx3kCool5), while
Figure 1(b) corresponds to the strong cooling case (γsyn/σ= 0.2;
3dx3kCool02). The small ripples that are visible in the
perpendicular-to-field direction (y-z plane) are caused by the
relativistic drift-kink 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 Figure 1. The lines in the figure
represent the magnetic field lines and highlight the main features
of the current sheet, including the X-points 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 Section 2.1 (this includes both the thin current sheet and
plasmoids of various sizes). Figure 2 shows the distribution of
widths, ( ) w , for a large collection of points for four
simulations with varying cooling strengths (in the order of
increasing cooling strength: 3dx3kCoolInf, 3dx3kCool5,
3dx3kCool1, and 3dx3kCool02). Two important conclu-
sions can be drawn from this plot. First, the peak of the ( ) w is
slightly smaller than ( )w r dL epeak

up s s~ = , which corre-
sponds to the Larmor radii of particles with the Lorentz factor
close to σ in the upstream magnetic field, Bup. The presence of
this peak is associated with the X-points, where particles are
accelerated to characteristic Lorentz factors of γ σ (see
Section 4.1). More importantly, current sheets with stronger
cooling display a steeper slope and shorter distribution span for

( ) w , indicating that the plasmoids in the current layer become
smaller and more concentrated as cooling strength increases (a
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, i.e., its rate, βrec. The evolution

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 nonideal 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 the 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).

4

The Astrophysical Journal, 959:122 (20pp), 2023 December 20 Chernoglazov, Hakobyan, & Philippov

MAC004054
Pencil



of the rate, measured as the inflow velocity of the plasma,
( ) ∣ ∣E B Bzrec

2b = ´ , averaged in the thin layer far upstream
parallel to the current sheet, for different cooling regimes is
shown in Figure 3. As mentioned earlier, we focus our analysis
in the time range t 1.5 Lx/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 (see 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 is approximately twice as large for the
stronger cooled runs when compared to the run without
cooling, as seen from the peak values in Figure 2. This
conclusion is also consistent with the characteristically smaller
plasmoid sizes in the strong cooling regime, as discussed
earlier.

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 velocityU N N Ns s s s,=m m

n
n ,

where ( )N f u u d us s
0 3òºm m , and fs and uμ are the distribution

and the four-velocity of particle species s. The plasma pressure
tensor for species s, Ps

mn , is then computed as the projection
of the stress-energy tensor: P Ts s= D Dmn

a
m

b
n ab, where Ts ºmn

( )p p p f d ps
0 3ò m n , U Us s shD º +mn mn m n is the projection

operator, and ημν is the Minkowski tensor (pμ=msu
μ is the

four-momentum of the particle of species s). We find that the
nondiagonal 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 nonnegligible ≈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 strong cooling.9 Guided 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
effect by removing the pressure support,10 which results in
extremely compressed structures (see, e.g., Hakobyan et al.
2019). In 3D, the picture is quantitatively different. The force-
balance condition is demonstrated in Figure 4, where we show
the cross-section of the plasma density in x-z plane for two
simulations in panels (a) and (c): one without cooling
(3dx3kCoolInf), and the other with strong cooling
(3dx3kCool02). We pick intermediate-size plasmoids in
both of these simulations (indicated by the dotted vertical lines
in the left-hand panels), and in panels (b) and (d) plot the
magnetic pressure, B2/8π, the scalar plasma pressure, defined
as P P1 3 s s i

i
,= å , as well as the z-components of the two

forces: the pressure gradient, c∇iP
iz, and the Lorentz force,

( )j B z´ , as functions of z. In the no-cooling simulation
(Figure 4(a), (b)), large plasmoids are typically weakly
magnetized, while in the strong cooling regime (Figures 4(c),
(d)), the strength of the magnetic field inside the plasmoids is
several times larger than the upstream value, Bup, (compare the
blue-solid lines in Figures 4(b), (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, 2023b; Schoeffler et al. 2023). Notably,
the strongest B-field component inside the plasmoid (in the
strong cooling regime) is the out-of-plane one, By. The
generation of this out-of-plane magnetic field component, as
well as the dissipation of the in-plane Bxz components (evident
in the no-cooling regime), are both consequences of the MHD-

Figure 2. Distribution of transverse sizes (along z) of the current layer for
different cooling regimes. The scale is zoomed in for ( )w rL

up s< . 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 rL

up s~ .

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.

9 The role of the pressure anisotropy due to synchrotron cooling is discussed
by Zhdankin et al. (2023).
10 This effect was investigated in detail by Schoeffler et al. (2019), who
observed that this compression provides a positive feedback loop: magnetic
field amplification by compression inside plasmoids leads to stronger
synchrotron cooling, which in turn promotes further compression.

5

The Astrophysical Journal, 959:122 (20pp), 2023 December 20 Chernoglazov, Hakobyan, & Philippov



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, By, 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 in both uncooled and strongly
cooled simulations is compressed and heated to a rough
equipartition with the upstream magnetic field pressure:
P B 8up

2 p» . In the latter case, the pressure is mainly delivered
by higher plasma density (compare densities in panels (a) and
(c) in Figure 4) because strong synchrotron cooling prevents
plasmoids from heating up by effectively radiating away the
components of particle momenta perpendicular to the magnetic
field. The orange lines in Figure 4(b) 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 (Figure 4(b)) and is reasonable
in the strongly cooled simulation (Figure 4(d)).

In this section, we largely neglected the presence of ions
(insusceptible to synchrotron cooling) in plasmoids. While their
contribution to the total pressure was self-consistently included
in our calculations, their relative contribution was typically
small for the timescales considered, e.g., in Figure 4, due to
their smaller number density. However, at later times the ions
can accelerate to much higher Lorentz factors than the pairs
(see Section 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 (Figures 4(a),
(b) and Appendix A). The applicability of this regime for the
current sheets in pulsar magnetospheres is studied in detail in
Section 6.1.

4. Particle Acceleration

During the nonlinear 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.
2023). 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. More-
over, as described earlier, 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. Meanwhile, in 3D, another acceleration
channel has been observed for uncooled pair plasmas by Zhang
et al. (2021a, 2023), with the large-scale reconnection electric
field in the upstream, Erec, being its main driver.
The dynamics of plasma in X-points is insensitive to the

presence of synchrotron cooling because 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.
Meanwhile, ions 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 Figure 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.
Furthermore, ions (dashed lines), regardless of the cooling
strength, extend to Lorentz factors significantly above σ, with
the cutoff energy itself increasing in time. In the uncooled

Figure 4. A cross-section of typical plasmoids for uncooled ((a), (b)) and strongly cooled ((c), (d)) current sheets. The left-hand 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, B2/8π, and gas
pressure, P, with the corresponding magnetic tension force, (j × B)z, and the gradient of the pressure tensor, c∇iP

iz. On the one hand, 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.

6

The Astrophysical Journal, 959:122 (20pp), 2023 December 20 Chernoglazov, Hakobyan, & Philippov



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 Section 4.1).

We further quantify the contributions of different accelera-
tion channels in Figure 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 nonideal electric field, E> B, during their
lifetime (dashed lines in Figure 6). In the inset panels of
Figure 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 σ
(see 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). Thus, the main energy
gain occurs due to the acceleration by ideal electric fields. For
pairs, as we argue further in Section 4.2, the “pick-up”
acceleration by outflows from the X-points allows them to gain
energies up to σ even in the strong cooling case.

4.1. 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
Figure 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 Lx/c (corresponding to the same moment shown in
Figure 5, when the spectrum shape is in the steady-state), and

follow their evolution until the end of the simulation
t= 4.5 Lx/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, 2023), which is evident from the thickness of the
line representing the trajectory in Figure 7. These trajectories
correspond to bouncing of accelerated ions between the two
converging upstream flows.11 In contrast, the low-energy ion
displays a tangled trajectory due to its entrapment inside a
plasmoid or frequent scattering on magnetic field inhomogene-
ities 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 toward the large-scale reconnection-driven
electric field, Erec≈Ey, being the dominant acceleration

Figure 5. 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 lines correspond to 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, rL

up.

Figure 6. Distribution functions of the pairs (blue) and the ions (orange) for
both the uncooled (a) and the strongly cooled (b) simulations. The dashed lines
show parts of the distributions contributed by the particles having experienced
E > B during their lifetime in the sheet. The inset panels show the distribution
of energy gains, Δγ, of particles during their contiguous presence in these
regions. The nonideal 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.

11 Note that the classical Speiser orbits correspond to the acceleration inside
the current sheet, as opposed to the acceleration upstream (for further
discussion, see Zhang et al. 2021a).

7

The Astrophysical Journal, 959:122 (20pp), 2023 December 20 Chernoglazov, Hakobyan, & Philippov



mechanism. The rate of acceleration for these trajectories can
be written as

· ⟹ ( ) E
vq

m c c

q B

m c
t, 7i

i
y

i

i
rec rec

upg g b b= D » D⎛
⎝

⎞
⎠

where we used Erec= βrecBup (Zhang et al. 2021a). This slope
is overplotted with a dashed line in the inset of the left-hand
panel of Figure 7, showing excellent agreement with the
acceleration history. Since ions are by construction not
susceptible to synchrotron cooling, we find that this accelera-
tion 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, which is highlighted with
dashed lines in Figure 5. Two parameters from Equation (7)
that could potentially contribute to different acceleration rates
are the reconnection rate, βrec,

12 and the beaming of ions

toward the direction of the reconnection electric field, βy. To
inspect the relative importance of these two effects, in Figure 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 simulation,
while the last two columns represent ions from the simulation
with the strongest cooling, γsyn/σ= 0.2. Since we are primarily
interested in the post-injection phase, in the odd-numbered
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 Figures 8(b1) and (b3)), then 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 mechanism). This prelude is consistent with the fact
that the ions at these energies are primarily moving in the plane

Figure 7. 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-
hand side accelerates to high energies (?σ) by the large-scale electric field, Erec, in the region upstream of the layer, while the ion on the left-hand side 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 is trapped inside the current sheet and experiences constant scatterings on local inhomogeneities.

12 As was already mentioned in Section 3, the global reconnection rate is
larger for the strongest cooling regime, which in turn means that the large-scale
electric field, Ey, is typically stronger.

8

The Astrophysical Journal, 959:122 (20pp), 2023 December 20 Chernoglazov, Hakobyan, & Philippov



of the reconnecting magnetic field, x-z, i.e., perpendicular to the
reconnection electric field, βy≈ 0, (see French et al. 2023). At a
later time, a fraction of these particles begin accelerating
upstream of the layer along the y-direction, i.e., along Erec; that
fraction of “lucky” particles ultimately determines the power-
law slope in the distribution. As is evident from the first column
of Figure 8, in the uncooled run only a small fraction of ions
follows the · E vrecg µ 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 Figure 2),
making them more capable of trapping the freely accelerating

ions. Meanwhile, 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 Figure 8; these
trajectories are similar to the one shown in the right-hand
column of Figure 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 Figure 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

Figure 8. 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 with 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.

9

The Astrophysical Journal, 959:122 (20pp), 2023 December 20 Chernoglazov, Hakobyan, & Philippov



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
Figures 8(c1)...(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 Figure 3, which in turn
means larger Ey= Erec≈ βrecBup. Another notable feature is the
evolution of the distribution of electric field values, Ey, 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 primarily move in the upstream (see
Figure 7), where the electric field is coherent and almost
exclusively points in the y-direction.

4.2. 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 nonideal 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 Erec. 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 Bup. This implies that
the upstream acceleration in the weak cooling regime is limited
by the burnoff limit, γ γsyn. In Figure 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 because the strong synchrotron cooling
disallows particles to leave the X-points by crossing magnetic
field lines. As shown in Figure 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 nonideal electric field and the pick-up
mechanism. In Figure 6(b), 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
Figure 9. In panels (a)–(e) we plot the x-z slices of the plasma
density, and we plot 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-hand
side of Figure 9, we show the time evolution of the Lorentz
factor γ, with the dashed lines indicating the times of the
snapshots on the left-hand side. In panel (a), the particle is
upstream from the layer and moving toward the current sheet
with a small Lorentz factor, γ∼ 1. In panel (b), the particle
enters the X-point and is accelerated by the nonideal 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 bound-
ary, which is perturbed vertically due to RDKI (see below). As
the particle hits the upstream, it experiences subtle energy
losses, which are indicated as oscillations in its energy, γ. After
that, the particle is picked up by a contracting magnetic field
line performing a slingshot-like motion. Then, as shown in
panel (c), for an extended period of time the particle moves
along with a rapidly receding field line while being slowly
accelerated to γ σ. Once the slingshot-like motion of the field

Figure 9. Time evolution of a typical accelerated electron in the strongly
cooled simulation. The right-hand 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-hand 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 further energy (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).

10

The Astrophysical Journal, 959:122 (20pp), 2023 December 20 Chernoglazov, Hakobyan, & Philippov



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
detail later on). 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 Figures 10(a)–(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 Figure 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
Figure 10, we specifically focus on the time when the particle
experiences successive scattering events while moving in the
RDKI-perturbed current layer,13 after being pre-energized in
the X-point (corresponding to Figure 9(b)). 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, ≈Bupγsyn/γ (see below), after which deceleration
occurs due to the radiative drag force, fdrag (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
Figure 10 (panels (c) and (d); see also Figure 9(d)),14 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 slingshot-like trajectory.
Before t≈ 3.5, although the velocity of the particle is
perpendicular to the local magnetic field (as indicated in
Figure 10(c), and Figures 9(c) and (d)), the effective
perpendicular magnetic field, B̃̂ is almost zero, which allows
the particle to avoid being radiatively cooled. At time t≈ 3.55,
when the particle enters into the dense flux tube, B̃̂ increases
dramatically (above the critical value, Bup(γsyn/γ)) and the
particle experience a strong radiative drag force (as shown in
the second panel on the lower right-hand side of Figure 10 with
an orange line, fdrag). In this catastrophic cooling episode, the
particle loses most of its energy and is further trapped inside the
flux tube (also shown in Figure 9(e)).

In both the scattering and the catastrophic cooling examples,
the full radiation reaction force, given by Equation (3), consists
of both synchrotron and curvature contributions. As shown in
the top and the bottom panels on the right-hand side of
Figure 10, in both cases the local radius of curvature of the field
line15 during the deceleration is comparable to the character-
istic 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 mechan-
isms 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, the pairs experience a highly

inhomogeneous magnetic field while propagating inside the
current sheet. The 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|Erec, and the
radiative drag force, ∣ ∣ ( ˜ ) ( )e E B Brec up

2
syn

2g g^ , one can define

the critical value for the effective magnetic field: B̃ »^
( )Bup syng g . As discussed earlier, 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 toward local weak-field regions
or outflows from the X-points.
To directly measure the effect of energy-dependent field

inhomogeneities, in Figure 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 Bup syng g»^ . Evidently, for all the cases
particle distributions do not cross the critical threshold, with the
highest-energy particles, γ≈ σ, following trajectories with
˜ ( )B Bup syng g^ (see 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, B̃p

2e gµ ^, and in the context of the
highest-energy pairs this determines the cutoff of the overall
radiation spectrum. The results shown in Figure 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 nontrivial dependence
of B̃̂ on the particle energy, especially in systems with
dynamically strong synchrotron cooling.
As we have seen earlier, in strongly cooled simulations the

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 ( )tc Lrecg b w »

( ( ) )( )r cLrec
1 upb s g s- , where ωL= |e|Bup/(mec).

16 We directly
measure the continuous residence time, Δt, at high energies,
γ γthr, for all of the pairs in our strongly cooled runs
(γsyn= 0.2σ). The distribution of residence times, normalized
to ( )r cL

up s , is shown in Figure 12 for two values of γthr: 2γsyn,
and σ. To ensure that 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 tc, 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,
Lx/c, which implies two important points. First, this indicates
that our simulation domain is large enough for the pairs to be
free to experience many acceleration-deceleration events

13 A 3D movie of this moment is available via the following link: https://
youtu.be/AILs-clOaPY.
14 A 3D movie focusing on this event can be found via the link: https://youtu.
be/zWn4Co1mLJU.
15 (( ˆ · ) ˆ)b b 1 - , where ˆ ∣ ∣b B Bº is a unit vector along the local magnetic field.

16 Here we used the definition of γsyn from (4) and approximated B̃ »^
( )Bup syng g . Then one can write ( ) ( )t f m c m c B4c e e Tdrag

1 2 2g p s g= = -
^ , and

substitute γsyn, and B̃̂ .

11

The Astrophysical Journal, 959:122 (20pp), 2023 December 20 Chernoglazov, Hakobyan, & Philippov

https://youtu.be/AILs-clOaPY
https://youtu.be/AILs-clOaPY
https://youtu.be/zWn4Co1mLJU
https://youtu.be/zWn4Co1mLJU


Figure 10. Close-up view of the acceleration dynamics of the electron shown in Figure 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 Figure 1. The middle panel on the right-hand side 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, fdrag · v, to the energy evolution of the particle, g . Finally, we show the curvature radius of the
magnetic field at the location of the particle. An animation of these acceleration dynamics is available in the online journal. The animation consists of two separate
sequences, corresponding to panels (a) and (b) and panels (c) and (d) in this figure. Each animation includes a 3D visualization of particle evolution and two different
line plots showing the energy evolution of the particle.

(An animation of this figure is available.)

12

The Astrophysical Journal, 959:122 (20pp), 2023 December 20 Chernoglazov, Hakobyan, & Philippov



during their lifetime inside the current layer. Despite this, the
pairs are unable to gain energy indefinitely in successive
acceleration events because 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̃ Bup syng g»^ , allowing us to extrapolate the upper
bound on pair acceleration to the current layers with realistic
scale separations.

5. 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 Section 4.2, we demonstrated that in the weak-to-
moderate cooling regime, γsyn σ, pair acceleration is limited
to max syng g . 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 maxg s~ (see, e.g., Figure 5), and we
expect the radiation spectrum to extend above the burnoff limit,
as follows from Equation (6). Following the discussion in

Section 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 Bup syng g»^ :

˜
( )B

B
. 8c ccut

up

max

syn

2

syn
e e

g
g

e
s
g

= =^
⎜ ⎟⎜ ⎟ ⎜ ⎟
⎛

⎝

⎞

⎠

⎛

⎝

⎞

⎠

⎛

⎝

⎞

⎠
To verify our expectations, in Figure 13 we show photon

spectra, dn d f2 2e e eºe e, from simulations with different
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):

( ) ( )dn

d
K x dx, 9

p
5 3

p
òe

e
e
e

µe

e e

¥

where K5/3 is the modified Bessel function, and the energy εp is
defined according to Equation (6). All of the spectra rise
linearly at low energies, ε2fε∝ ε, which is expected from the

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 syng g=^ . 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. 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 cL

up s and, importantly, is independent of the domain size.

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 Equation (8): εcut ∼ εc(σ/γsyn). However, current sheets with
weak cooling produce almost no emission above the burnoff limit.

13

The Astrophysical Journal, 959:122 (20pp), 2023 December 20 Chernoglazov, Hakobyan, & Philippov



hard power law of reconnection-accelerated 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

( )f exp2
cute e eµ -e . The inferred values, εcut, are shown in

Figure 13 with colored vertical lines, as well as in the inset
panel at the top right-hand side. Evidently, our results are
consistent both with the prediction of Equation (8) for the
strong cooling, εcut≈ εc(σ/γsyn) and 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, 2012a, 2013; Kagan et al.
2016a; Mehlhaff et al. 2020) or in 3D simulations with a
limited separation of scales (Cerutti et al. 2014). In Figure 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σ
(Figure 14(a)), the highest-energy pairs are accelerated in the
upstream, along the reconnection electric field, Erec≈ Ey. This
results in high-energy photons being emitting preferentially
perpendicular to the upstream magnetic field (Zhang et al.
2021a, and Section 4.2). For the strong cooling (Figure 14(c),
(d)), the emission of photons with energies above the burnoff
limit, ε few · εc, is strongly anisotropic and is beamed toward
the direction of the upstream field. The bulk of these photons is
emitted by the mechanism described in Section 4.2 and
particularly demonstrated in the lower half of Figure 10: the
particle is accelerated by the pick-up mechanism in the X-point
outflow, in the direction of the upstream magnetic field
(Figure 10(c)), and later produces the short radiation burst
after colliding face-on with the plasmoid.

To understand the spatial distribution of emitting regions, in
Figure 15 we plot a volumetric map of the radiated synchrotron
power, ·F vP esync rad= å  , with the Frad defined in Equation (3)

for simulations with weak (Figure 15(a)), and strong cooling
(Figure 15(b)); 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 Figure 1 and also plot the total plasma density.
To properly quantify the radiated power, we normalize it by the
incoming Poynting flux, ( )cB 4rec up

2b p .
In the case of weak cooling (Figure 15(a)), 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 toward
plasmoids. This is due to the fact that a substantial fraction of
the highest-energy particles are accelerated by Erec in the upstream
of the current sheet and trapped by plasmoids. Meanwhile, for the
strong cooling regime (Figure 15(b)), 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 Figures 10(a), (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 (Figures 10(c), (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 Figure 9(e)). Because of
this, plasmoid interiors in the strong cooling regime 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 (Figure 15(a)), which is
consistent with our earlier discussion, and is also highlighted in
the emission spectra in Figure 13. For the strong cooling
(Figure 15(b)), the power above the burnoff limit (blue) is
similar to the total power (green), the main difference is that the
highest-energy emission is suppressed in plasmoid interiors.

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 ˆB xBxup º . 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. However, in the strong cooling regime, most of the high-energy photons above the burnoff limit are
emitted along the upstream magnetic field.

14

The Astrophysical Journal, 959:122 (20pp), 2023 December 20 Chernoglazov, Hakobyan, & Philippov



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 inhomogene-
ities (such as the RDKI ripples, see also Figure 10).

6. 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 accelera-
tion is effective for both 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, 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 Section 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 reconnec-
tion 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≈ 16MeV,
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 the 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 Erec and perpendicular to
the upstream magnetic field Bup.

6.1. Acceleration of Ions in Pulsar Magnetospheres

Although current sheets in the magnetospheres of pulsars are
predominantly filled with the electron-positron plasma, 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 the definition of the
synchrotron burnoff Lorentz factor, Equation (4), taking the
Thomson cross-section of ions, ( )Z A m mT

i
T e p

4 2 2s s= ,
mi= Amp (mp is the mass of a proton), and qi= Z|e|, we find
that for the ions ( )AZ m mi

p esyn
3 2

syng g= - , with γsyn being the
burnoff Lorentz factor for pairs (see Section 2.2). Since most of
the high-energy emission takes place in regions of the current

Figure 15. Spatial distribution of regions producing the synchrotron radiation shown together with the total plasma density for the same snapshots as in Figure 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) corresponds to 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.

15

The Astrophysical Journal, 959:122 (20pp), 2023 December 20 Chernoglazov, Hakobyan, & Philippov



layer close to the light cylinder, we may further employ
( ) B B R RLC LC

3» as a proxy for Bup, with Rå, and Bå being
the radius of the neutron star, and the magnetic field strength near
its surface, calculated as ( )( )  

B P P10 G 1 s 1012 15= - , and
RLC= cPå/2π being the light cylinder radius (På is the rotation
period of the star). We then find

( ) A

Z

P B
10

1 s 10 G
, 10i

syn
10

3 2

3 2

12

1 2

g »
-

⎛
⎝

⎞
⎠

⎛
⎝

⎞
⎠

which can be compared with the total potential drop in the
current layer, qiβrecBLCRLC, normalized to mic

2:

· ( ) Z

A

P B
4 10

1 s 10 G
. 11i

max
3

2

12
g »

-
⎛
⎝

⎞
⎠

⎛
⎝

⎞
⎠

The assumption of negligible cooling for ions then holds when
i i

syn maxg g . 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 Figure 16. The gray-
shaded region in the parameter space indicates the strong
cooling regime for ions, i i

syn maxg g< , 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 Figure 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 106, i.e., above
PeV energies. For heavier nuclei with Z≈ A> 1, the cooling is
slightly more efficient, as Zi i

syn max
1 2g g µ - , 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, then 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, fi≈ (mi/me)α/λ≈ 0.01 (which is the value that we
employed in our simulations), where λ= n±/nGJ∼ 104 is the
pair multiplicity, and α= ni/nGJ∼ 0.1 is the ion density relative
to the Goldreich–Julian density,17nGJ= |Ω ·B|/2πce, where Ω
is the angular velocity of the pulsar. The energy budget of even
this small fraction of ions, nimic

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 i

syn maxg g ), the energy of
pairs is limited to max LCg s» , while the distribution of ions can
extend to the full potential drop, i

maxg , Equation (11).
Assuming a power-law index of −1 for ions and 〈γ±〉≈ γsyn
for pairs (see our simulation 3dx3kCool02), we may rewrite
the total energy budget for both species as:

( )




n m c

n m c

,

. 12

e

i i
i

GJ
2

syn

GJ
2

max

l g

a g

»

»


Meanwhile, the magnetization near the light cylinder, σLC, can be
expressed as ( )(∣ ∣ )P e B m c4 5eLC LC maxs pl g l= =  , where

∣ ∣e B R m cemax rec LC LC
2g b= 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 ( )( )  5i syn LCa g s= , where

we used the relation m me i
i

max maxg g= . For a Crab-like pulsar
with γsyn/σLC≈ 0.1, we find that pairs dominate in terms of the
total energy budget, i.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

maxg without significant deviations from the
f (γ)∝ γ−1 predicted in our simulations. In fact, we do not
observe the spectrum to steepen, even when   1 3i » , as
happens in the simulation 3dx3kCool02 at late times.
Another channel of ion acceleration in pulsar magneto-

spheres 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∼ 106...107. The same potential drop will corre-
spond to the maximum Lorentz factor of protons

10 ... 10i
th

3 4g ~ (smaller by a factor of mi/me). These values
are much smaller than the ones corresponding to the
acceleration in the current sheets of Crab-like pulsars,

 10i
max

6g , 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.

Figure 16. Estimates for burnoff limit of protons, i
syng , (blue lines) and the

Lorentz factor of protons accelerated in the full potential drop near the light
cylinder, i

maxg , (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 i
max syng g> . Evidently, for all of the observed pulsars, the radiation drag of

ions is negligible in the context of ion acceleration.

17 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).

16

The Astrophysical Journal, 959:122 (20pp), 2023 December 20 Chernoglazov, Hakobyan, & Philippov



6.2. High-energy Emission from Young Pulsars and
Supermassive Black Holes

In Section 5, we discussed the predictions for the observed
high-energy spectrum from the accelerated pairs. To summar-
ize, 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 (see Hakobyan
et al. 2023a). For 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 the Crab pulsar, where
γsyn≈ 4 · 104 close to the light cylinder, this implies that the
magnetization parameter near the current layer is around
σLC≈ 106. Notably, the presence of pairs with energies
few · 106mec

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 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
because the naive expectation implied that particles were
beamed in the direction of the reconnecting electric field,
perpendicular to the direction of the magnetic field. However,
our findings outlined in Section 5 and in Figure 14 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, which is 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.
2020, 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
σ∼ 107, and γsyn∼ 106...107. In that regard, our simulations
are largely consistent with the results of 2D simulations by
Hakobyan et al. (2023b), where the authors predict that most of
the power is radiated at 10...100 MeV range, with the most
energetic pairs, ( )max , 10syn

7g s g~ ~ , 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,
which suggests 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. 2018) as a possible mechanism for the gamma-ray
flares at energies 160MeV 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=B2/(4πnmec

2h), where h∼ 〈γ〉 corresponds to the plasma
enthalpy per particle. The strong cooling regime then applies only
if σh〈γ〉 γsyn∼ 109, where the estimate for γsyn is made
assuming that the strength of the magnetic field is close to
milligauss. Our findings suggest that a significant flux at energies
above 160MeV can only be produced when · 10hs
( ) ( ) ( )1000 10 10syn syn

9 7 1g g g gá ñ » á ñ - . It is currently unclear
whether such conditions can be realized in the Nebula.

Acknowledgments

We thank James Drake, Lorenzo Sironi, and Dmitri
Uzdensky for the enlightening conversations, and the anon-
ymous referee for the useful comments. This work was
supported by NSF grant No. PHY-2231698, and facilitated
by Multimessenger Plasma Physics Center (MPPC), NSF grant
No. PHY-2206607. This work was supported by a grant from
the Simons Foundation (00001470, AP). Computing resources
were provided and supported by the Division of Information
Technology at the University of Maryland (Zaratan
cluster).18 This research is part of the Frontera computing
project at the Texas Advanced Computing Center (LRAC-
AST21006). Frontera is made possible by NSF award OAC-
1818253. The computing resources for this research were
partially provided and supported by Princeton Institute for
Computational Science and Engineering and Princeton
Research Computing. H.H. was partially supported by the
U.S. Department of Energy under contract No. DE-AC02-
09CH11466.

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, fi= ni/nup, as well as the mass ratio, mi/me.
For this study, we use simulations with a spatial resolution of
de=Δx and a box size Lx/de= 1000. In Section 4.1, we have
demonstrated that in simulations with no cooling and mass
ratio mi/me= 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 panels (a)–(c) of Figure 17 we show a comparison of the

distribution of ions measured in the steady state of the
reconnection, t= 4Lx/c, for varying number density fractions,
fi= 1%, 10%, and 49%, while keeping mi/me= 1 (the
simulation shown in Figure 17(a) 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, fi= 1, their distribution remains very hard,

18 http://hpcc.umd.edu

17

The Astrophysical Journal, 959:122 (20pp), 2023 December 20 Chernoglazov, Hakobyan, & Philippov

http://hpcc.umd.edu


dn di
1g gµ - for energies γ σ. In a case where the number

of positrons is low and ions constitute almost half of the
particle population (Figure 17c), we observe a significant
deviation, with the resulting ion distribution, dn di

1.7g gµ - ,
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 because their acceleration primarily occurs in the
X-point region.

Another limitation of our simulations is the reduced mass
ratio, mi/me, 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, rL, and the size of the box, and to
minimize the effects of boundary conditions on particle
acceleration. To confirm that the reduced mass ratio
mi/me= 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:
mi/me= 5 and mi/me= 25. To keep the mass fraction of ions
constant, ρi/ρe= 2%, we also decreased their number density
fraction, fi. The results of these simulations are shown in
Figure 17(d). Evidently, the power-law slope of the distribution
is only marginally affected by the mass ratio, with the
power-law tail extending to energies ?meσ (where σ takes
into account the corresponding mass ratio as shown in
Equation (1)).

Appendix 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, de= 1...3Δx,
while keeping mi/me= 1, and fi= 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,
Lx/de= 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 because 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 Figure 18. Panels (a) and (b) show slices
in y= 0 plane of the plasma density at time t= 2.7Lx/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 (Figures 18(a), (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 nonideal field E>B is identical for
both simulations, showing that the dynamics around X-points is
equally resolved for both simulations.

Figure 17. 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.

18

The Astrophysical Journal, 959:122 (20pp), 2023 December 20 Chernoglazov, Hakobyan, & Philippov



ORCID iDs

Alexander Chernoglazov https://orcid.org/0000-0001-
5121-1594
Hayk Hakobyan https://orcid.org/0000-0001-8939-6862
Alexander Philippov https://orcid.org/0000-0001-
7801-0362

References

Abdo, A. A., Ajello, M., Allafort, A., et al. 2013, ApJS, 208, 17
Alves, E. P., Zrake, J., & Fiuza, F. 2018, PhRvL, 121, 245101
Amato, E., & Arons, J. 2006, ApJ, 653, 325
Ansoldi, S., Antonelli, L. A., Antoranz, P., et al. 2016, A&A, 585, A133
Arons, J. 2003, ApJ, 589, 871
Bai, X.-N., & Spitkovsky, A. 2010, ApJ, 715, 1282
Begelman, M. C. 1998, ApJ, 493, 291
Beloborodov, A. M. 2017, ApJ, 850, 141
Bransgrove, A., Ripperda, B., & Philippov, A. 2021, PhRvL, 127, 055101
Buehler, R., Scargle, J. D., Blandford, R. D., et al. 2012, ApJ, 749, 26
Cerutti, B., Philippov, A. A., & Spitkovsky, A. 2016, MNRAS, 457, 2401
Cerutti, B., Uzdensky, D. A., & Begelman, M. C. 2012a, ApJ, 746, 148
Cerutti, B., Werner, G. R., Uzdensky, D. A., & Begelman, M. C. 2012b, ApJL,

754, L33
Cerutti, B., Werner, G. R., Uzdensky, D. A., & Begelman, M. C. 2013, ApJ,

770, 147
Cerutti, B., Werner, G. R., Uzdensky, D. A., & Begelman, M. C. 2014, ApJ,

782, 104
Chen, A. Y., & Beloborodov, A. M. 2014, ApJL, 795, L22
Chernoglazov, A., Ripperda, B., & Philippov, A. 2021, ApJL, 923, L13
Comisso, L., & Sironi, L. 2018, PhRvL, 121, 255101
Dahlin, J. T., Drake, J. F., & Swisdak, M. 2015, PhPl, 22, 100704
Davelaar, J., Philippov, A. A., Bromberg, O., & Singh, C. B. 2020, ApJL,

896, L31
Fang, K., Kotera, K., & Olinto, A. V. 2012, ApJ, 750, 118
French, O., Guo, F., Zhang, Q., & Uzdensky, D. 2023, ApJ, 948, 19
Guépin, C., Cerutti, B., & Kotera, K. 2020, A&A, 635, A138
Guo, F., Li, H., Daughton, W., & Liu, Y.-H. 2014, PhRvL, 113, 155005
Guo, F., Li, X., Daughton, W., et al. 2021, ApJ, 919, 111

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

Hakobyan, H., Philippov, A., & Spitkovsky, A. 2019, ApJ, 877, 53
Hakobyan, H., Philippov, A., & Spitkovsky, A. 2023a, ApJ, 943, 105
Hakobyan, H., Ripperda, B., & Philippov, A. A. 2023b, ApJL, 943, L29
Hakobyan, H., Spitkovsky, A., Chernoglazov, A., et al., 2023

PrincetonUniversity/tristan-mp-v2: v2.6, v2.6, Zenodo, doi:10.5281/
zenodo.75667256

Huang, Y.-M., & Bhattacharjee, A. 2016, ApJ, 818, 20
Jaroschek, C. H., & Hoshino, M. 2009, PhRvL, 103, 075002
Kagan, D., Nakar, E., & Piran, T. 2016a, ApJ, 826, 221
Kagan, D., Nakar, E., & Piran, T. 2016b, ApJ, 833, 155
Kalapotharakos, C., Brambilla, G., Timokhin, A., Harding, A. K., &

Kazanas, D. 2018, ApJ, 857, 44
Landau, L. D., & Lifshitz, E. M. 1975, The Classical Theory of Fields (4th ed.;

Amsterdam: Elsevier)
Lyubarskii, Y. E. 1996, A&A, 311, 172
Lyubarsky, Y. E. 2005, MNRAS, 358, 113
Lyutikov, M., Sironi, L., Komissarov, S., & Porth, O. 2018, JPlPh, 84,

635840201
Manchester, R. N., Hobbs, G. B., Teoh, A., & Hobbs, M. 2005, AJ, 129, 1993
Mehlhaff, J. M., Werner, G. R., Uzdensky, D. A., & Begelman, M. C. 2020,

MNRAS, 498, 799
Ortuño-Macías, J., Nalewajko, K., Uzdensky, D. A., et al. 2022, ApJ, 931, 137
Parfrey, K., Philippov, A., & Cerutti, B. 2019, PhRvL, 122, 035101
Petropoulou, M., & Sironi, L. 2018, MNRAS, 481, 5687
Philippov, A. A., & Spitkovsky, A. 2014, ApJL, 785, L33
Philippov, A. A., & Spitkovsky, A. 2018, ApJ, 855, 94
Ripperda, B., Bacchini, F., & Philippov, A. A. 2020, ApJ, 900, 100
Ripperda, B., Liska, M., Chatterjee, K., et al. 2022, ApJL, 924, L32
Rybicki, G. B., & Lightman, A. P. 1979, Radiative Processes in Astrophysics

(New York: Interscience (Wiley-Interscience))
Schoeffler, K. M., Grismayer, T., Uzdensky, D., Fonseca, R. A., & Silva, L. O.

2019, ApJ, 870, 49
Schoeffler, K. M., Grismayer, T., Uzdensky, D., & Silva, L. O. 2023, MNRAS,

523, 3812
Servidio, S., Dmitruk, P., Greco, A., et al. 2011, NPGeo, 18, 675
Sironi, L. 2022, PhRvL, 128, 145102
Sironi, L., & Beloborodov, A. M. 2020, ApJ, 899, 52
Sironi, L., Giannios, D., & Petropoulou, M. 2016, MNRAS, 462, 48
Sironi, L., Rowan, M. E., & Narayan, R. 2021, ApJL, 907, L44

Figure 18. Comparison of 2D plasma density slices ((a) and (b)) and normalized particle distributions ((c), and(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. However, the compression of plasmoids is
somewhat smaller for the low-resolution run.

19

The Astrophysical Journal, 959:122 (20pp), 2023 December 20 Chernoglazov, Hakobyan, & Philippov

https://orcid.org/0000-0001-5121-1594
https://orcid.org/0000-0001-5121-1594
https://orcid.org/0000-0001-5121-1594
https://orcid.org/0000-0001-5121-1594
https://orcid.org/0000-0001-5121-1594
https://orcid.org/0000-0001-5121-1594
https://orcid.org/0000-0001-5121-1594
https://orcid.org/0000-0001-5121-1594
https://orcid.org/0000-0001-5121-1594
https://orcid.org/0000-0001-8939-6862
https://orcid.org/0000-0001-8939-6862
https://orcid.org/0000-0001-8939-6862
https://orcid.org/0000-0001