ABSTRACT

Title of Dissertation: RADIATIVE PLASMAS IN
PULSAR MAGNETOSPHERES

Alexander Chernoglazov
Doctor of Philosophy, 2024

Dissertation Directed by: Assistant Professor Alexander Philippov
Department of Physics

Pulsars are highly magnetized rotating neutron stars known for their periodic bursts of

radio emission. Decades of astronomical observations revealed that pulsars produce non-thermal

radiation in all energy bands, from radio to gamma rays, covering more than 20 decades in photon

energy. Modern theories consider strongly magnetized relativistic electron-positron plasmas to

be the source of the observed emission. In my Thesis, I investigate physical processes that can

be responsible for plasma production and the observed high-energy emission in the wide range

of photon energies, from eV to TeV.

In the first Chapter of my Thesis, I investigate relativistic magnetic reconnection with

strong synchrotron cooling using three-dimensional particle-in-cell kinetic plasma simulations. I

characterize the spectrum of accelerated particles and emitted synchrotron photons for varying

strengths of synchrotron cooling. I show that the cutoff energy of the synchrotron spectrum can

significantly exceed the theoretical limit of 16 MeV if the plasma magnetization parameter ex-

ceeds the radiation reaction limit. Additionally, I demonstrate that a small fraction of ions present



in the current sheet can be accelerated to the highest energies, making relativistic radiative recon-

nection a promising mechanism for the acceleration of high-energy cosmic rays.

In the second Chapter, I present the first multi-dimensional simulations of the QED pair

production discharge that occurs in the polar region of the neutron star. This process is believed

to be the primary source of the pair plasma in pulsar magnetospheres and also the source of

the radio emission. In this work, I focus on the self-consistently emerging synchronization of

the discharges in different parts of the polar region. I find that pair discharges on neighboring

magnetic field lines synchronize on a scale comparable to the height of the pair production region.

I also demonstrate that the popular “spark” model of pair discharges is incompatible with the

universally adopted force-free magnetospheric model: intermittent discharges fill the entire polar

region that allows pair production, leaving no space for discharge-free regions. My findings

disprove the key assumption of the spark model about the existence of distinct discharge columns.

In the third Chapter, I demonstrate how the key findings of two previous chapters can pro-

vide a self-consistent explanation of the recently discovered very-high-energy, reaching 20TeV,

pulsed emission in Vela pulsar. Motivated by the results of recent global simulations of pulsar

magnetospheres, I propose that this radiation is produced in the magnetospheric current sheet un-

dergoing radiative relativistic reconnection. I show that high-energy synchrotron photons emitted

by reconnection-accelerated particles efficiently produce electron-positron pairs. The density of

secondary pairs exceeds the supply from the polar cap and results in a self-regulated plasma mag-

netization parameter of ∼ 107. Electrons and positrons accelerate to Lorentz factors comparable

to ∼ 107 and emit the observed GeV radiation via the synchrotron process and 10 TeV photons

by Compton scattering of the soft synchrotron photons emitted by secondary pairs. My model

self-consistently accounts for the ratio of the gamma-ray and TeV luminosities.



RADIATIVE PLASMAS IN PULSAR MAGNETOSPHERES

by

Alexander Chernoglazov

Dissertation submitted to the Faculty of the Graduate School of the
University of Maryland, College Park in partial fulfillment

of the requirements for the degree of
Doctor of Philosophy

2024

Advisory Committee:
Assistant Professor Alexander Philippov, Chair/Advisor
Professor James Drake
Professor Michael Coleman Miller
Professor Anatoly Spitkovsky
Dr. Marc Swisdak



Preface

Chapter 2 of this thesis is based on work in collaboration with Prof. Alexander Philippov

and Hayk Hakobyan and published in the Astrophysical Journal as Chernoglazov et al. (2023).

Chapter 3 is based on paper a paper in collaboration with Prof. Andrey Timokhin and Prof.

Alexander Philippov and submitted for publication in the Astrophysical Journal Letters. Chapter

4 presents the preliminary results to be submitted to the Astrophysical Journal Letters.

During my time in graduate school, I had a privilege to work on numerous additional

projects. A significant portion of my time in graduate school was spent working on studying

the dynamics of binary compact objects, and I significantly contributed to a publication Foucart

et al. (2021). During my internship at the Center for Computational Astrophysics, I performed

an extensive study of the relativistic highly-magnetized turbulence, both in strong (Chernoglazov

et al., 2021) and weak (Ripperda et al., 2021; TenBarge et al., 2021) limits, in the fluid approx-

imation. The results and methods that I developed during my work on relativistic reconnection

received further development in studies of neutrinos from coronae of Active Galactic Nuclei

(Mbarek et al., 2024).

ii



Acknowledgments

Professionally and personally, I owe my gratitude to my advisor, Professor Alexander

Philippov, for creating a supportive and creative environment in the research group. His in-

valuable support has helped me work on extremely interesting and challenging problems over the

past many years. I specifically acknowledge the brainstorming we did together during my work

on various projects. His support and guidance are not limited to scientific life but extended to the

times when I asked for his advice at the earlier stage of my career.

I would like to thank all the people I worked with at the Center for Computational Astro-

physics during the first pandemic year: Professor Bart Ripperda, Professor Jens Mahlmann, and

Professor Elias Most. I also want to thank my collaborators on many projects, both included

and not included in this thesis: Dr. Hayk Hakobyan, Professor Andrey Timokhin, and Dr. Ros-

tom Mbarek. I especially want to thank people I met at the early stage of my PhD track at the

University of New Hampshire: Professor Francois Foucart and Professor Benjamin Chandran.

Numerous friendships have enlivened my time in graduate school, including Mariia, Alisa,

Nikita, Anastasia, Alina, Hayk, Tanya, Dima, Olya and Alex. I want to thank them for their

continuous support during my study in graduate school in this turbulent time.

Finally, I must thank my parents for their perpetual support over the years. Their support

of my decisions and encouragement of my interest in science have been integral to my path in

life since my childhood.

iii



Table of Contents

Preface ii

Acknowledgements iii

Table of Contents iv

List of Tables vi

List of Figures vii

List of Abbreviations ix

Chapter 1: Introduction 1
1.1 Radiation from Pulsar Magnetospheres . . . . . . . . . . . . . . . . . . . . . . . 1
1.2 Magnetic Reconnection . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
1.3 Particle-In-Cell Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11

Chapter 2: High-Energy Radiation and Ion Acceleration in Three-dimensional Rela-
tivistic Magnetic Reconnection with Synchrotron Cooling 17

2.1 Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
2.2 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
2.3 Numerical Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21

2.3.1 Initial configuration and boundary conditions . . . . . . . . . . . . . . . 21
2.3.2 Radiation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25

2.4 Structure of the Current Sheet . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27
2.5 Particle acceleration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34

2.5.1 Ion acceleration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38
2.5.2 Pair acceleration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43

2.6 Radiation from the current sheet . . . . . . . . . . . . . . . . . . . . . . . . . . 52
2.7 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57

2.7.1 Acceleration of ions in pulsar magnetospheres . . . . . . . . . . . . . . . 58
2.7.2 High-energy emission from young pulsars and supermassive black holes . 63

Chapter 3: Coherence of Multi-Dimensional Pair Production Discharges in Polar Caps
of Pulsars 66

3.1 Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66
3.2 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67

iv



3.3 Overview of the physical model . . . . . . . . . . . . . . . . . . . . . . . . . . 70
3.3.1 Dynamics of the electromagnetic fields . . . . . . . . . . . . . . . . . . 70
3.3.2 Magnetospheric current distribution . . . . . . . . . . . . . . . . . . . . 72
3.3.3 QED pair production . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74
3.3.4 Atmosphere . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77
3.3.5 Initial plasma state . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78
3.3.6 Numerical details . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79

3.4 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80
3.4.1 Space-Charge-Limited Flow . . . . . . . . . . . . . . . . . . . . . . . . 81
3.4.2 Ruderman-Sutherland Gap . . . . . . . . . . . . . . . . . . . . . . . . . 90
3.4.3 3D Gap, Space-Charge-Limited Flow . . . . . . . . . . . . . . . . . . . 92
3.4.4 Evolution of the Twist of Magnetic Field Lines . . . . . . . . . . . . . . 95

3.5 Conclusions and Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 98
3.5.1 Non-existence of sparks . . . . . . . . . . . . . . . . . . . . . . . . . . 98
3.5.2 Implication for discharges in old pulsars . . . . . . . . . . . . . . . . . . 101
3.5.3 On repetition rate of the cascades . . . . . . . . . . . . . . . . . . . . . . 101

Chapter 4: Origin of the Pulsed 20 TeV Emission in the Vela Pulsar 103
4.1 Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103
4.2 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103
4.3 Magnetospheric Current Sheet . . . . . . . . . . . . . . . . . . . . . . . . . . . 106

4.3.1 Plasma Feeding from the Polar Cap . . . . . . . . . . . . . . . . . . . . 106
4.3.2 Particle Acceleration and Synchrotron Radiation from the Current Sheet . 108

4.4 Model of the HE and VHE Emission . . . . . . . . . . . . . . . . . . . . . . . . 111
4.4.1 Pair Production . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 113
4.4.2 Steady-state Pair Spectrum . . . . . . . . . . . . . . . . . . . . . . . . . 116
4.4.3 Secondary Pair Synchrotron Spectrum . . . . . . . . . . . . . . . . . . . 118
4.4.4 Inverse Compton Signal . . . . . . . . . . . . . . . . . . . . . . . . . . 119

4.5 Conclusion and Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 121

Chapter 5: Conclusion and Future Work 124
5.1 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 124
5.2 Future Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 125

Appendix A: Appendices for Chapter 2 128
A.1 Varying the fraction of ions and the mass ratio . . . . . . . . . . . . . . . . . . . 128
A.2 Varying the resolution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 130

Appendix B: Appendices for Chapter 3 132
B.1 Additional Details on radiation and QED physics . . . . . . . . . . . . . . . . . 132

v



List of Tables

2.1 Summary of Simulation Parameters for Current Sheets . . . . . . . . . . . . . . 25

3.1 Summary of Simulation Parameters for Discharge Simulations. . . . . . . . . . . 77

vi



List of Figures

1.1 Schematic of the pulsar magnetosphere . . . . . . . . . . . . . . . . . . . . . . . 3
1.2 Schematic of the magnetic reconnection . . . . . . . . . . . . . . . . . . . . . . 8
1.3 Schematic of the PIC algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . 14

2.1 Volume Rendering of the Current Sheet . . . . . . . . . . . . . . . . . . . . . . 28
2.2 Distribution of Perpendicular Sizes of the Current Sheet . . . . . . . . . . . . . . 29
2.3 Time Evolution of the Reconnection Rate . . . . . . . . . . . . . . . . . . . . . 31
2.4 The Internal Structure of a Plasmoid . . . . . . . . . . . . . . . . . . . . . . . . 33
2.5 Distribution Functions of Pairs and Ions . . . . . . . . . . . . . . . . . . . . . . 36
2.6 Acceleration of Particles in X-points . . . . . . . . . . . . . . . . . . . . . . . . 37
2.7 Trajectories and Energy Gain of Ions . . . . . . . . . . . . . . . . . . . . . . . . 39
2.8 Time Evolution of Sample of Ions . . . . . . . . . . . . . . . . . . . . . . . . . 41
2.9 Time Evolution of a Typical High-Energy Electron . . . . . . . . . . . . . . . . 45
2.10 Close-up View of an Electron Acceleration Dynamics . . . . . . . . . . . . . . . 47
2.11 Distribution of the Effective Perpendicular Magnetic Field . . . . . . . . . . . . 50
2.12 Distribution of the Lifetimes of Strongly Cooled Pairs . . . . . . . . . . . . . . . 51
2.13 Radiation Spectra from Current Sheets . . . . . . . . . . . . . . . . . . . . . . . 53
2.14 Energy-dependent Angular Anisotropy of the Radiation . . . . . . . . . . . . . . 55
2.15 Spatial Distribution of Radiating Regions . . . . . . . . . . . . . . . . . . . . . 55
2.16 Estimates for Burnoff Limit of Protons in Pulsar Magnetospheres . . . . . . . . . 61

3.1 Distribution of the magnetospheric currents . . . . . . . . . . . . . . . . . . . . 73
3.2 SCLF discharge with small gap . . . . . . . . . . . . . . . . . . . . . . . . . . . 83
3.3 SCLF discharge, different gap models . . . . . . . . . . . . . . . . . . . . . . . 87
3.4 RS discharge, different gap sizes . . . . . . . . . . . . . . . . . . . . . . . . . . 91
3.5 3D SCLF discharge . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 94

4.1 Plasma supply to the current sheet . . . . . . . . . . . . . . . . . . . . . . . . . 107
4.2 Particle and photon spectra from a current sheet . . . . . . . . . . . . . . . . . . 110
4.3 Particle and photon spectra in Vela . . . . . . . . . . . . . . . . . . . . . . . . . 117

5.1 Vlasov simulation of the polar discharge . . . . . . . . . . . . . . . . . . . . . . 127

A.1 Energy Spectra for Different Fractions of Ions . . . . . . . . . . . . . . . . . . . 129
A.2 Comparison of Different Resolutions for Current Sheets . . . . . . . . . . . . . . 131

vii



B.1 1D model of SCLF discharge . . . . . . . . . . . . . . . . . . . . . . . . . . . . 136

viii



List of Abbreviations

CSS Current Sheet Synchrotron
FFE Force-Free Electrodynamics
GJ Goldreich-Julian
HE High Energy
H.E.S.S. High Energy Stereoscopic System
IC Inverse Compton
LC Light Cylinder
NS Neutron star
PC Polar Cap
PIC Particle-In-Cell
QED Quantum Electrodynamics
RS Ruderman-Sutherland
SCLF Space-Charge-Limited Flow
SPS Secondary Pairs Synchrotron
VHE Very High Energy

ix



Chapter 1: Introduction

1.1 Radiation from Pulsar Magnetospheres

Pulsars are lighthouses of the Universe. Being initially discovered as sources of extremely

stable periodic bursts of bright radio emission (Hewish et al., 1968), they were quickly linked to

rapidly rotating neutron stars: compact astrophysical objects with characteristic radii of 12 km

and masses of 1.4 − 2.2M⊙, and typical rotational periods P ≈ 10 ms − 10 sec. The next

≳ 55 years of observations show that pulsars are broad-band emitters of non-thermal radiation,

covering about 20 orders of magnitude in photon energy (see, e.g., Philippov & Kramer, 2022).

The studies of pulsar magnetospheres were revolutionized by the detection of their γ-ray emission

using Fermi, with more than 200 pulsars detected until now (Abdollahi et al., 2020): for the first

time, observations opened a window into an energetically significant process operating in pulsar

magnetospheres.

Early theoretical models of pulsar magnetospheres were based on a model of a vacuum

magnetic dipole, µ, rotating with angular velocity Ω (Pacini, 1967), which radiates its rotational

energy through the Poynting flux equivalent to a luminosity L = (2/3)(µ2Ω4)/(c3) sin2 χ (Lan-

dau & Lifshitz, 1975), where χ is an angle between the rotational axis and the magnetic dipole,

c is the speed of light. Poynting flux is supported by the rotation-induced electric field, which is

quadrupolar for the dipolar magnetic field. Such an electric field has a parallel to the magnetic

1



field component that is able to extract charged particles from the stellar surface and accelerate

them to enormous energies. The fatal flaw of the vacuum model is its inconsistency with obser-

vations, requiring charged particles to produce the observed emission. The presence of plasma

leads to modifications of the vacuum fields; in particular, it leads to a screening of the parallel

component of the electric field. Under the assumption of a plasma-filled magnetosphere, mag-

netic field lines loaded with plasma cannot corotate with the neutron star beyond a certain radius.

This distance is called the light cylinder radius and is determined as a distance where the coro-

tational speed is equal to the speed of light, RLC = c/Ω. While the magnetic field lines that do

not cross the light cylinder remain mostly not distorted, the ones that cross have to charge their

geometry. This condition separates the pulsar magnetosphere into two regions: the zone of the

closed field lines and the zone of the open field lines. These two zones and the light cylinder are

shown in the sketch 1.1.

Open field lines originate in a small region on the stellar surface called the polar cap. Its

characteristic size is RPC ≈ R⋆

√
R⋆/RLC, where R⋆ is a stellar radius. The exact size and shape

of the polar cap depends on the mutual orientation of the magnetic µ and rotational Ω axes. The

field-aligned flow of charged particles extracted from the polar cap produces an electric current,

I . It generates the twist of the magnetic field, resulting in a toroidal component of the magnetic

field,Bφ = 2I(r⊥)/cr⊥. This toroidal magnetic field together with corotational electric field Ecor

create an outward component of the Poynting flux, (c/4π)(Ecor ×Bφ). The corotational electric

field, Ecor = −Ω×r×B/c, is realized by a minimally required charge density, ∇·Ecor ≈ 4πρGJ,

where, ρGJ = −Ω ·B/(2πc) is the Goldreich-Julian charge density (Goldreich & Julian, 1969).

Integration of the Poynting flux through the light cylinder results in an estimate of the energy

2



Ω

μ

cu
rre

nt 
sh

ee
t

separatrix

γ
γ

e+
e−

⟹

⟹

jmag/jGJ > 0

jmag/jGJ < 0

RLC

closed 
field lines

NS

open field lines
γ

B

e−

e−e+
⟹

E||

Figure 1.1: Sketch of the magnetosphere of the 45◦-inclined pulsar. The main elements of the
magnetosphere are shown: zones of open and closed field lines, magnetospheric current sheet,
and separatrix. The shadowed region of the open field lines zone highlights the field lines carrying
volume return current. The zoom-in shows the scheme of the polar discharge.

3



losses

L ≈ µ2Ω4

c3
∼ 1035erg/s, (1.1)

for typical surface magnetic field strength, B = 1012G, and rotational period, P = 0.1s. This

estimate turns out to be close to the vacuum dipolar losses, with the main difference being that the

losses are now non-zero for the aligned rotator as well. Since the Poynting flux carries away the

rotational energy of the star, it produces a spindown, Ė = −IΩΩ̇ = 4π2IṖ /P 3. For a moment

of inertia of the star, I = 1045g × cm2, the increase of the rotational period is small, Ṗ ∼ 10−15,

but observable and provides an important constraint on the strength of the magnetic field. The

luminosity of the radio emission, which pulsars are known for, accounts only for a small fraction

of the total spindown luminosity: Lradio ∼ 10−7 − 10−2Lspindown (Lorimer & Kramer, 2005).

Its brightness temperature1, however, is incredibly high, Tb ∼ 1025 − 1040K (Hankins & Eilek,

2007), requiring a coherent emission mechanism.

In the Goldreich-Julian model, the accelerating electric field is screened by the plasma ex-

tracted from the stellar surface, and no particle acceleration occurs. The model was modified by

Sturrock (1971), who assumed that the non-zero accelerating electric field exists in the magneto-

sphere. It leads to the production of additional pair plasma via quantum electrodynamics (QED)

effects. The driver behind this is a huge electric potential drop across the polar cap, V = ErotRPC.

Electric fields lead to efficient particle acceleration along the curved magnetic field lines and sub-

sequent emission of high-energy curvature photons with the energies εcurv = (3/2)ℏ(c/ρ)γ3,

where ρ is the curvature radius of the magnetic field line. These photons, emitted initially along

the magnetic field lines, propagate along a straight line, and the angle between the wave vector

1The brightness temperature is the temperature needed for a blackbody radiator to produce the same intensity as
the observed source.

4



and the magnetic field lines increases. As a result, the probability of the photon decay into an

electron-positron pair via a single photon conversion process increases, γ+B → e−+e+ (Erber,

1966). The sketch of this process is shown in the inset of Figure 1.1. The produced pairs have

a non-zero pitch angle and intensively radiate transverse energy into synchrotron photons which

are again converted into new electron-positron pairs. During the acceleration process, a primary

charged particle emits multiple energetic photons, and each of them is converted into pairs. As

a result, the magnetosphere is filled with a high-density plasma with number density n = λnGJ,

which screens the parallel electric field. Here, nGJ = ρGJ/e is the characteristic density of parti-

cles extracted from the stellar atmosphere (Beloborodov, 2008; Timokhin & Arons, 2013). The

dimensionless parameter λ≫ 1 is called the multiplicity and characterizes how many secondary

particles can be produced by a single primary particle accelerated in the gap – the zone where the

unscreened electric field could exist. The multiplicity can be estimated from energy conservation:

estimates for Lorentz-factors of the primary particles, γth ∼ 106−107, and of the secondary ones,

γsec ∼ 102−103 (Timokhin & Harding, 2015), result in the multiplicity λ ∼ γth/γsec ∼ 104−105.

The process of particle acceleration with subsequent secondary plasma creation and screening of

the electric field is called a discharge. It is considered the primary source of pair plasma in the

magnetosphere.

Ruderman & Sutherland (1975) developed these ideas further and produced the first quan-

titative model of the pair discharge. They demonstrated that the position of the pulsars’ “death

line” corresponds to the potential drop across the polar cap, Φ ≈ 1012V. At lower electric po-

tential, primary particles are not accelerated to high enough energy to ignite the pair-production

discharge, and the pulsar is no longer able to produce radio emission. Initially, the discharge was

considered as a stationary process (e.g., Ruderman & Sutherland, 1975; Arons & Scharlemann,

5



1979). Later, using first-principles plasma simulations, Timokhin (2010) demonstrated that the

discharge proceeds intermittently: intense bursts of pair-production are followed by quiet mo-

ments when the QED-produced plasma screens an electric field. After the plasma leaves the gap

zone, the unscreened electric field forms again, and the process restarts. The cyclic behavior of

the discharge provides a natural explanation for the observed microstructure of individual radio

pulses. Recently, the process of cyclic screening of the electric field has been directly connected

to the formation of the radio emission (e.g., Philippov et al., 2020; Melrose et al., 2021a). In this

Thesis, I develop the first self-consistent multi-dimensional models of the pair discharge.

In the plasma-filled magnetosphere, the bundles of the open field lines from two opposite

polar caps cross the light cylinder and approach each other. Since their magnetic field lines have

opposite polarities, a discontinuity appears in the region of contact of these two fields - a mag-

netospheric current sheet (shown in blue in Figure 1.1). This current sheet undergoes magnetic

reconnection and is the site of the efficient particle acceleration and subsequent γ-ray emission

(Lyubarskii, 1996; Cerutti et al., 2015; Philippov & Spitkovsky, 2018). Other possible sites of

the γ-ray emission are possible in the magnetosphere: the escaping high-energy photons from the

gap region close to the stellar surface, and the curvature radiation of the energetic particles propa-

gating along the boundary of open and closed field lines, the separatrix, carrying intense currents

(Contopoulos & Kalapotharakos, 2010). Global magnetospheric simulations show that the energy

release in the inner magnetosphere is insignificant (Hu & Beloborodov, 2022; Hakobyan et al.,

2023a). Observations rule out high-energy radiation from near the surface region: data show no

super-exponential cutoff of the spectrum at high energies expected from single-photon conver-

sion. The mismatch of rotational phases at which the radio and γ-ray light curves peak (e.g.,

Kuiper & Hermsen, 2015) also disfavors the inner magnetosphere as the source of the observed

6



high-energy radiation (the Crab pulsar and some other very energetic pulsars are an exception,

Eilek & Hankins 2016). At the same time, Fermi observations demonstrate that the high-energy

emission contains ∼ 10% of the total pulsar spindown luminosity (Abdo et al., 2010a,b, 2013;

Ansoldi et al., 2016a), pointing out to the magnetic reconnection as the particle energization pro-

cess. Emission produced near and beyond the light cylinder accumulates in caustics: all of the

photons, emitted by relativistic particles propagating to the observer, travel together with an emit-

ting particle and arrive to an observer simultaneously, leading to the pile-up of the photons in the

skymap (Bai & Spitkovsky, 2010; Kalapotharakos et al., 2014). This caustic results in narrow

peaks in γ-ray lightcurves, consistent with observations.

1.2 Magnetic Reconnection

Magnetic reconnection is believed to be the most efficient mechanism of converting mag-

netic field energy into plasma energy and is accompanied by the change of the magnetic field

topology schematically shown in Fig. 1.2a. Magnetic reconnection events were directly observed

near the surface of the Sun, in the Solar corona, and also in the Earth magnetotail (Yamada

et al., 2010). It appears plausible that the reconnection can occur in remote astrophysical systems

as well, although its presence can not be demonstrated directly. Around compact astrophysical

objects, neutron stars, and black holes, magnetic reconnection occurs in highly magnetized en-

vironments with βp ≡ 8πPth/B
2 ≪ 1, and transfer of the magnetic field energy into particles

significantly affects their dynamics.

The switch of the field polarity necessitates the presence of a large current confined in the

current sheet. The characteristic width of the current sheet in magnetized pair plasma is given by

7



Bup

vin

vout

Bup

current sheet plasmoid

(a) (b)

Figure 1.2: Sketch of the current sheet formation. Subplot (a) shows the formation of a short
Sweet-Parker-like current sheet. Arrows indicate the directions of the plasma inflow and outflow
and the magnetic field orientation. Subplot (b) shows a long plasmoid-unstable current sheet with
a developed plasmoid chain.

a skin depth, de = c/ωpl,e where ωpl,e =
√
4πne2/me is the plasma frequency. At smaller scales,

the magnetic field is no longer frozen into the plasma, and the rearrangement of the magnetic

field and consecutive energy release are spontaneously triggered.

The first analytical steady-state model of reconnection was proposed by Parker (1957) and

Sweet (1958), the Sweet-Parker model. In their model, two opposite polarities of the magnetic

field are advected into the current sheet with a constant speed, vin. In the steady-state, the as-

sociated motional electric field, E = −vin ×B/c, is constant everywhere, including inside the

current sheet, due to the continuity of the tangential electric field at the boundaries of the up-

stream zones and the current sheet. The inflow and outflow mass fluxes are equal, resulting in

nδvout = nLvin. Here, n is the plasma density, δ is the half-width of the current sheet, and L is

the length of the interface of the opposite magnetic field regions. The dimensionless parameter

r = vin/vout = δ/L is called the “reconnection rate” and characterizes the efficiency of advec-

tion and dissipation of the magnetic field. The outflow velocity, vout, can be found from energy

conservation under the assumption that all of the magnetic energy is converted into the kinetic

8



energy of plasma, B2/8π = nmev
2
out/2. Hence, the outflow velocity is equal to the Alfven ve-

locity, vout = B/
√
4πρ ≡ vA. For the plasma resistivity η, Ohm’s law determines the thickness

of the current sheet, J ≈ B/δ ≈ vinB/(ηc). These results determine the reconnection rate as

r =
1√
S
, S =

vAL

η
− Lundquist number. (1.2)

This fact reveals the main drawback of the Sweet-Parker theory. For astrophysical systems,

the Lundquist number is huge, S ∼ 1020, leading to an extremely low reconnection rate and,

hence, inefficient particle acceleration. However, the current sheets of such high Lundquist num-

bers are always unstable. Furth et al. (1963) showed that the thinning of the current sheet leads

to the growth of the tearing instability. Biskamp (1986); Loureiro et al. (2007) demonstrated

that long enough current sheets, with S ≳ 104, have to split into smaller ones on the dynami-

cal timescales. As a result of continuous fragmentation into the magnetic islands, or plasmoids,

current sheets maintain multiple evolving reconnection sites (a sketch of the plasmoid-unstable

current sheet is shown in Fig. 1.2b). These plasmoids were directly observed in flares in the Solar

corona (Lin et al., 2005) and in the Earth magnetotails (Zong et al., 2004). The fragmentation of

a long current sheet marks the “fast reconnection” regime when the reconnection rate becomes

independent of the length of the current sheet and is about 10%, r ≈ 0.1, for collisionless plasmas

(e.g., Zenitani & Hoshino, 2008; Bessho & Bhattacharjee, 2012).

Similar conclusions hold for reconnection in collisionless relativistic plasmas (Lyubarsky,

2005). Here, the ratio of the supplied magnetic energy and plasma energy is characterized by a

9



dimensionless magnetization parameter:

σ =
B2

4πnmec2
. (1.3)

For magnetized environments of black holes and neutron stars, σ ≫ 1, making reconnection

proceed in the relativistic regime. For these conditions, the Alfven velocity reaches the speed of

light, vA ≈ c
√
σ/(σ + 1) ≈ c, and the reconnection rate can be defined as r ≈ βrec ≡ vin/c.

Following the conversion of the magnetic field energy into plasma energy at reconnection sites,

the characteristic energies of the accelerated particles become comparable to the magnetization

parameter, E ∼ mec
2σ, and in the relativistic regime particles form a hard non-thermal power-

law distribution (e.g., Guo et al., 2014; Sironi & Spitkovsky, 2014): f(E) ∝ E−1. Subsequent

acceleration to energies exceeding the value of the magnetization parameter happens only for a

small fraction of particles, and three-dimensional simulations show steepening of the distribution

function, f(E) ∝ E−2 (Zhang et al., 2021a).

The formation of current sheets is ubiquitous in astrophysical systems. In some of them,

it is determined by the global geometries of the magnetic field, e.g., magnetospheres of pulsars

(e.g., Spitkovsky, 2006; Philippov et al., 2015; Cerutti et al., 2020; Hakobyan et al., 2023a) and

accreting black holes (e.g., Ripperda et al., 2020; Galishnikova et al., 2023). Additionally to

global current sheets, localized reconnection events can occur in magnetized relativistic turbu-

lence in both collisional (Zhdankin et al., 2013; Chernoglazov et al., 2021; Galishnikova et al.,

2022) and collisionless (e.g., Comisso & Sironi, 2018, 2019) regimes, due to the scale-dependent

intermittency of turbulence (Boldyrev, 2006; Mallet et al., 2015; Chandran et al., 2015). In some

astrophysical systems, the system-sized current sheet forms intermittently and is believed to drive

10



powerful flares in erupting magnetospheres of magnetars (e.g., Beloborodov, 2017a; Mahlmann

et al., 2023), kink-unstable relativistic jets (e.g., Davelaar et al., 2020), black holes surrounded by

magnetically-arrested disks (Ripperda et al., 2022), and in coronae of x-ray binaries and Active

Galactic Nuclei (Beloborodov, 2017b). In this Thesis, I investigate three-dimensional reconnec-

tion in conditions appropriate for pulsar magnetospheres.

In the vicinity of compact astrophysical objects, electromagnetic radiation emitted by ac-

celerated particles is so powerful that the radiation reaction force significantly alters the dynamics

of particles and, thus, of reconnection itself (Uzdensky, 2016). In pulsar magnetospheres and,

likely, in magnetospheres of black holes, electron-positron pairs created by high-energy photons

significantly affect the mass loading. Thus, self-consistent predictions of the high-energy emis-

sion from reconnecting current sheets require the inclusion of radiative effects into kinetic simu-

lations of relativistic reconnection. The 2D regime of radiative reconnection has been thoroughly

explored in recent years (Jaroschek & Hoshino, 2009; Uzdensky & McKinney, 2011; Cerutti

et al., 2012b,a; Beloborodov, 2017b; Werner et al., 2019; Sironi & Beloborodov, 2020; Mehlhaff

et al., 2020; Sridhar et al., 2021). However, little progress has been made in studying three-

dimensional reconnection, incorporating at least some of the radiative effects self-consistently,

which motivates my explorations in Chapter 2.

1.3 Particle-In-Cell Method

The dynamics of collisionless plasmas in many astrophysical systems are determined by the

self-consistent interplay between the charged particles and electromagnetic fields. This evolution

11



is described by the Vlasov equation (e.g., Lifshitz & Pitaevskii, 1981)

∂fs
∂t

+ v · ∇fs −
qs
ms

(
E +

v ×B

c

)
· ∂fs
∂v

= 0, (1.4)

which is a kinetic equation for the 6D distribution function f(x,v, t) of charged particles of

specie s, with charge qs and mass ms. E and B are the coarse-grained electric and magnetic

fields consistent with the charge density and electric currents created by the charged particles in

the system. The evolution of the electromagnetic field is described by Maxwell’s equations:

1

c

∂B

∂t
= −∇×E,

1

c

∂E

∂t
= ∇×B − 4π

c
j, (1.5)

and satisfies the constraints ∇ · E = 4πρ and ∇ · B = 0. Here, the charge, ρ, and current, j,

densities are given by the distribution function fs as

ρ =
∑
s

qs

∫
fsd

3v, j =
∑
s

qs

∫
vfsd

3v. (1.6)

Equations (1.4-1.5) for a collisionless plasma form a set of hyperbolic evolution equations with

a closure (1.6). These equations can rarely be solved analytically in the non-linear regime, re-

quiring the adoption of numerical techniques. The straightforward numerical solution of the

kinetic equation (1.4) requires evolution of the distribution function in 6-dimensional space and

is extremely computationally expensive. The most widely applied method for solving the Vlasov

equations is the Particle-In-Cell (PIC) algorithm, which I employ in this Thesis.

The PIC algorithm (Birdsall & Langdon, 1991) discretizes the distribution function by

sampling it with a finite number of macro-particles. These particles can be considered as a large

12



collection of real particles with mass ms and charge qs. At any given moment t, the distribution

function is expressed as a sum over the sampling macro-particles

fs(x,v, t) =
∑
i

S(x− xi(t))δ(v − vi(t)), (1.7)

where xi(t), vi(t) are position and velocity of a macro-particle i at the moment of time t. While

the macro-particle has a determined velocity vi(t), its spatial distribution is extended and is de-

scribed by a shape-function S. For the relativistic plasmas, further discussed in the Thesis, the

equation of motion of a separate particle can be conveniently described by the evolution of the

4-velocity ui = viγi. Here, γi = 1/
√

1− v2i /c
2 is the Lorentz factor of a particle. Thus, the

equation of motion of each particle is given by

dxs,i

dt
=

us,i

γs,i
(t), (1.8)

dus,i

dt
=

1

ms

Fs,i. (1.9)

Consistency with the Vlasov equation (1.4) requires spatial averaging of the Lorentz force, Fs,i,

using the same shape function, S:

Fs,i = qs

(
E(rs,i) +

us,i ×B(rs,i)

γs,ic

)
. (1.10)

While the macroparticles sample the plasma distribution function, the electric and magnetic

fields are sampled on a spatial grid. The scheme of such a grid is shown in Fig.1.3a. The staggered

positions of the field components in the cell are chosen to satisfy the divergence-free constraint

13



Computing the EM force on 
particles


∇ × E = − 1
c

∂B
∂t

∇ × B = 1
c

∂E
∂t

+ j

                             


                                          
dpi

dt
= q (E + v × B

c ) − Frad

Charge weighting on the grid

Solving Maxwell’s equations on the 
grid

Pushing particles

Δt

ρ, jE, B

γβ, r

dxi

dt
= 1

mΓ pi

Plasma injection

Ez, Jz

Ex, Jx

Ey, Jy

By

Bz

Bx

(i, j, k) (i + 1,j, k)

(i, j, k + 1)

(i + 1,j + 1,k)

(a)

(b)

Figure 1.3: Schematic representation of the PIC algorithm: Subplot (a) shows the sampling of the
fields and electric currents on the grid (Yee, 1966). Subplot (b) shows the key steps performed
by a PIC algorithm at each timestep.

14



to the machine precision and discretize the spatial derivatives with the second-order numerical

accuracy (Yee, 1966). The particle current is collected from macroparticles via the current deposit

scheme, which can be done to satisfy the charge-conservation equation to the machine’s precision.

Further description of the method will be based on PIC code Tristan v2 (Hakobyan et al.,

2023), which I used for simulations in this Thesis. The scheme of the algorithm is presented in

Fig. 1.3b.

The main steps of the PIC algorithm are:

1. Computation of the change of particle 4-velocity during the timestep due to the Lorentz

force, Eqn. (1.10), and the corresponding coordinate update, using a leap-frog algorithm

(Boris, 1970).

2. Deposit of the particle electric currents j onto the grid. This process needs to preserve the

charge conservation constraint, ∇·E = 4πρ. Tristan v2 code uses the zigzag method,

which conserves charge with machine precision (Umeda et al., 2003).

3. Computation of the time evolution of the electromagnetic field, E,B, governed by Maxwell

equations, using the grid-based values of the fields and the electric current. As with most

other PIC codes, Tristan v2 adopts the leap-frog scheme for the time update. In Chap-

ter 3, I adopt higher-order spatial discretization (Blinne et al., 2018) in order to improve

the dispersive properties of the field solver.

4. Interpolation of the updated electromagnetic field from the grid into the position of each

particle. Here, similar to the current deposit, the spatial non-locality of the particles is used.

Sampling of the continuous distribution function by the finite number of macroparticles leads

15



to the introduction of the numerical shot noise, which often results in parasitic heating of the

plasma (e.g., Shalaby et al., 2017). There are two main ways to reduce it: to increase the num-

ber of computational particles or to use a smoother shape function (Birdsall & Langdon, 1991).

Since the first method is numerically expensive and only slowly reduces the noise, Tristan

v2 employs the second method by applying digital filter passes on the deposited currents with a

(1/4, 1/2, 1/4) stencil.

In order to model relativistic radiative plasmas in pulsar magnetospheres, I have developed

a synchrotron emission module, which takes into account the synchrotron kernel, and the com-

plete QED pair production scheme that models the single-photon decay process, γ+B → e−+e+,

into Tristan v2. I describe details of these numerical implementations in Chapters 2 and 3

and corresponding Appendices.

16



Chapter 2: High-Energy Radiation and Ion Acceleration in Three-dimensional

Relativistic Magnetic Reconnection with Synchrotron Cooling

2.1 Abstract

We present the results of 3D particle-in-cell (PIC) simulations that explore relativistic mag-

netic reconnection in pair plasma with strong synchrotron cooling and a small mass fraction of

non-radiating ions. Our results demonstrate that the structure of the current sheet is highly sen-

sitive 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 ≈ 16 MeV, is only efficiently produced in the strong cooling regime, γsyn < σ.

In this regime, we find that the spectral cutoff scales as εcut ≈ εc(σ/γsyn), and the highest en-

ergy 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

17



place constraints on the reconnection-driven models of gamma-ray flares in the Crab Nebula.

2.2 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 pro-

cess 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 & Loureiro, 2016; Uzdensky, 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., Philippov & Spitkovsky, 2014;

Chen & Beloborodov, 2014) and black hole magnetospheres (e.g., Parfrey et al., 2019; Brans-

grove et al., 2021; Ripperda et al., 2022). In other systems, such as the coronae of X-ray binaries

or jets of active galactic nuclei (AGNs), the formation of intermittent current layers associated

with either turbulent energy cascade or MHD instabilities (e.g., kink instability, Begelman, 1998;

Tchekhovskoy & Bromberg, 2016; Alves et al., 2018; Davelaar et al., 2020; Ortuño-Macı́as et al.,

2022) has been theorized (Beloborodov, 2017b) 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 re-

connection has thus established itself as one of the most important plasma-physical processes

powering energy dissipation and non-thermal emission in various contexts of high-energy astro-

physics.

18



Strong radiative cooling, both due to synchrotron emission and inverse-Compton scatter-

ings, has implications for particle acceleration and the reconnection process itself (see, e.g., Uz-

densky, 2016). The dynamics of the radiative reconnection has so far been mostly studied for 2D

current sheets (Jaroschek & Hoshino, 2009; Uzdensky & McKinney, 2011; Cerutti et al., 2012b,a;

Beloborodov, 2017b; Werner et al., 2019; Sironi & Beloborodov, 2020; Mehlhaff et al., 2020;

Sridhar et al., 2021). For the case of strong cooling (either synchrotron or inverse-Compton),

radiative losses lead to the decrease of the plasma temperature, resulting in a strong compression

of plasmoids (e.g., Schoeffler et al., 2019; Hakobyan et al., 2019). Radiative cooling also sub-

stantially limits non-thermal particle acceleration. However, in the case of strong synchrotron

cooling, the initial acceleration (“injection”) of particles takes place in or close to the X-points,

where the radiative cooling is negligible (Cerutti et al., 2014; Kagan et al., 2016a; Hakobyan

et al., 2019), allowing the particles to accelerate above the synchrotron “burn-off” limit and emit

γ-ray synchrotron photons (Uzdensky et al., 2011; Cerutti et al., 2012a; Kagan et al., 2016b).

The case of 3D relativistic current sheets is substantially less explored (Zenitani & Hoshino,

2008; Guo et al., 2014; Sironi & Spitkovsky, 2014; Werner & Uzdensky, 2017; Zhang et al.,

2021a; Werner & Uzdensky, 2021; Guo et al., 2021; Schoeffler et al., 2023). It has been shown

that there are competing instabilities and acceleration mechanisms specific to 3D. For exam-

ple, Zenitani & Hoshino (2008), and Cerutti et al. (2014) demonstrated that in the linear stage,

the relativistic drift-kink instability grows faster than the plasmoid instability. Later, Sironi &

Spitkovsky (2014), and Guo et al. (2014) showed that the plasmoid instability still dominates in

the non-linear regime, leading to the 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 ultimately leads to turbulence,

19



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 magne-

tization 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, particles

are capable of escaping the flux tubes and can be directly accelerated by the reconnection electric

field. These studies were mainly focused on the non-radiative regime, and it remains to be seen

how radiative cooling affects 3D current sheets (Cerutti et al., 2014; Schoeffler et al., 2023).

Another important aspect of the problem is the composition of current sheets. While most

studies of relativistic reconnection focus on pair plasmas, in many scenarios, such as the mag-

netospheres of pulsars (e.g., Guépin et al., 2020, see Sec. 2.7.1) and supermassive black holes,

we expect a mixture of ions (primarily protons) to be present. Importantly, ions are not suscepti-

ble 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 chapter, we present a set of large 3D PIC simulations of reconnecting current sheets

populated by the electron-positron plasma with a small mixture of ions. In our simulations, pairs

are subject to dynamically important synchrotron cooling, while ions are unaffected. The goal

of our study is to understand the emerging distribution functions of the accelerated ions and

pairs, their acceleration mechanisms, and the properties of high-energy radiation. We begin in

Sec. 2.3 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

20



compressed flux tubes (Sec. 2.4). In Sec. 2.5, we present our main findings on particle accelera-

tion. In particular, we show that ions are very efficiently accelerated in current sheets with strong

radiative cooling, while the Lorentz factors of pairs are limited by the upstream magnetization

parameter by the synchrotron losses. In Sec. 2.6 we discuss the radiation spectra and the beaming

of the emission, for varying strength of the synchrotron cooling. We find that 3D current sheets

with efficient cooling are capable of emitting high-energy photons above the synchrotron burn-

off limit, which are preferentially beamed along the upstream magnetic field. In contrast, current

sheets with weak cooling emit photons with energies up to the burn-off limit in the direction

perpendicular to the upstream field direction. We summarize our main findings and discuss their

implications for magnetospheres of pulsars and black holes in Sec. 2.7.

2.3 Numerical Methods

For our simulations, we use the Tristan v2multi-species radiative PIC code (Hakobyan

et al., 2023).

2.3.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 = Bupx̂ tanh (z/δcs), i.e., no guide field is present, where

x and z are the coordinates along and transverse to the current sheet, correspondingly, and δcs is

the initial half-thickness of the current sheet, and the thermal pressure in the current sheet is set

to compensate the gradient of magnetic pressure B2/8π. We initialize uniform thermal upstream

pair plasma with a mixture of ions of overall number density, nup = ni + ne + ne+ = 2ne, and

21



a low temperature, Tup/(mec
2) = 10−4. In all simulations presented in the main text, the ions

constitute a fraction fi = 0.01 of the total number of particles1, 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.2 These assumptions

are validated in the Appendix A.1, 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,

σ ≡ B2
up/4π

menupc2
[
(1− fi) + fi

mi

me

] ≈ B2
up/4π

menupc2
, (2.1)

is fixed at σ = 50. The current sheet is initialized as a hot electron-positron-ion plasma of

density ncs = Ncsnup cosh
−2 (z/δcs), 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, B2
up/8π. The initial thickness of the current

layer, δcs = 7de, is large enough so the tearing instability is not seeded by the numerical noise.

Our simulation domain has outflowing boundary conditions in the ±x̂-directions, with a

thick absorber that gradually damps E and B fields to zero towards the edges of the box. In the

±ẑ-direction we implement a moving plasma injector, similar to Sironi & Spitkovsky (2014);

1A simulation with no ions, but otherwise identical to 1dx1kCool02, demonstrates similar overall dynamics
of the reconnection process.

2This situation applies to the current sheets in pulsar magnetospheres and can be relevant to the pair-dominated
coronae of accreting black holes.

22



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, 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 min-

imize the effects of boundary conditions on the statistical steady state, the radiative cooling and

the outflow boundaries are switched on when the transient approaches the edge of the simulation

box, at 0.5 Lx/c, where Lx is the extent of the domain in the x-direction. In all of our analyses,

we exclude the particles from the initial current layer, essentially ignoring the evolution of the

sheet in its early transient stage. We consider the current sheet to enter the steady-state regime

after about t ≳ 1.5 Lx/c from the beginning of the simulation. We run all simulations for 3

additional light-crossing times, with a total duration of 4.5 Lx/c. All statistical data is averaged

over the final 3 light-crossing times, while the current sheet is in statistical steady-state.

To study the numerical convergence of our results, we explore resolutions of de/∆x =

1...3, where de = (mec
2/4πnupe

2)
1/2 is the upstream plasma skin depth, and ∆x is the cell size

(∆x = ∆y = ∆z); the comparison of these two resolutions is presented in Appendix A.2,

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,

rupL (γ) = γmec
2/(|e|Bup), accelerated up to γ ∼ σ: rupL (σ)/∆x = (de/∆x)

√
σ ≈ 7...21,

assuming particles gyrate in the perpendicular magnetic field equal to the upstream value. To

23



better characterize the acceleration mechanisms and minimize the boundary effects, we set our

fiducial box sizes to be in the range of Lx/r
up
L (σ) ≈ 150...300. Our timestep duration (the

CFL number) is fixed at c∆t/∆x = 0.225 for most of our simulations. The summary of all the

simulation parameters is given in Table 2.1; a discussion on how the synchrotron radiation is

implemented (last column) follows further (see Sec. 2.3.2). In our strongest cooling simulation,

1dx1kCool01, we employ a halved CFL number, c∆t/∆x = 0.1125, compared to our default

choice of 0.225, in order to properly resolve the plasma oscillation period in regions of dense flux

tubes, highly compressed due to strong radiative cooling (see Sec. 2.4). The upstream plasma

is presented with a total concentration of both electrons and positrons averaging to 1 particle

per cell3. To test the convergence, we run a simulation identical to 1dx1kCool02 but with 8

particles per cell in total and obtain identical results. To mitigate the numerical artifacts from the

finite number of particles per cell we employ 8 digital filter passes on the deposited currents with

a (1/4, 1/2, 1/4) stencil.

To analyze the properties of the current sheet in steady-state, we need to employ a robust

method to distinguish the current layer carrying the plasma energized in reconnection from the

fresh plasma upstream. The most convenient way is to separately trace particles from the two

upstreams (z ≳ 0 and z ≲ 0 regions), and define the mixing factor, M, as suggested by Zhang

et al. (2021a),

M = 1− 2

∣∣∣∣n↑
n

− 1

2

∣∣∣∣ . (2.2)

Here n↑ is the local number density of particles originating from the z > 0 region, and n is the

total number density. This definition implies M = 0 in either of the two upstream regions and

3I.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.

24



M ≈ 1 in the mixed region of the current layer. We use the threshold M = 0.7 to separate

the upstream from the reconnecting current sheet, however, our results are insensitive to this

particular choice, as the upstream-to-downstream transition is quite sharp in terms of the values

of M.

Table 2.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

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

Frad =
σT
4π

(κR − γ2B̃2
⊥β),

κR = (E + β ×B)×B + (β ·E)E,

B̃2
⊥ = (E + β ×B)2 − (β ·E)2.

(2.3)

25



Here, β = v/c is normalized 3-velocity, and γ = (1− β2)−1/2 is the Lorentz factor4. In mag-

netospheres 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 re-

connecting current sheet. In order to model a similar physical regime in simulations, we rescale

the Thomson cross-section σT , using a dimensionless parameter, γsyn (Uzdensky, 2018; Werner

et al., 2019; Hakobyan et al., 2019; Sironi & Beloborodov, 2020), defined as follows:

|e|Erec ≡
1

4π
σTB

2
upγ

2
syn. (2.4)

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.

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 (2.4), moving in the perpendicular magnetic field of strength Bup radiate

4The total radiated power, Psync ∝ γ2B̃2
⊥, corresponds to a square of the electric field in the rest frame of the

moving particle, E′ = γ(E + β ×B)− (γ − 1)(β ·E)β/β2.

26



photons with the characteristic energy

εc =
3

2
ℏ
|e|Bup

mec
γ2syn =

9

4αF

Erec

Bup

mec
2 ≈ 16 MeV, (2.5)

where αF ≈ 1/137 is the fine-structure constant (see also Uzdensky et al., 2011; Cerutti et al.,

2012a; Uzdensky & Loureiro, 2016). Importantly, this value is insensitive to the strength of the

magnetic field and is thus a perfect benchmark to calibrate the photon energies. Thus, in our

simulations, a particle with Lorentz factor of γ experiencing a magnetic field of B̃⊥, defined in

equation (2.3), by definition radiates a spectrum of photons which peaks near

εp = 16 MeV
B̃⊥
Bup

(
γ

γsyn

)2

. (2.6)

2.4 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 con-

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

27



To visualize the complex 3D structure, Fig. 2.1 shows volume renderings of the plasma

density in the fully-developed non-linear reconnection stage of 3D current sheets from two dif-

ferent simulations. Fig. 2.1a shows a density snapshot of the simulation with weak cooling

(γsyn/σ = 5; 3dx3kCool5), while Fig. 2.1b corresponds to the strong cooling case (γsyn/σ =

0.2; 3dx3kCool02). The small ripples, visible in the perpendicular-to-field direction (y-z

plane), are caused by the relativistic drift-kink instability. The imprint of the ideal kink insta-

bility is visible as the overall twist and deformation of the large flux tubes in both panels of

Fig. 2.1. The lines in the figure represent the magnetic field lines and highlight the main fea-

tures 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.

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

28



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 hori-

zontal (x-y) plane using the definition of the mixed region in Sec. 2.3.1 (this includes both the thin

current sheet and plasmoids of various sizes). Fig. 2.2 shows the distribution of widths, P(w), for

a large collection of points for four simulations with varying cooling strength (in the order of in-

creasing cooling strength: 3dx3kCoolInf, 3dx3kCool5, 3dx3kCool1, 3dx3kCool02).

There are two important conclusions one can draw from this plot. First of all, the peak of the

P(w) is slightly smaller than wpeak ∼ rupL (σ) = de
√
σ, which corresponds to the Larmor radii of

particles with the Lorentz factor close to σ in the upstream magnetic field, Bup. The presence of

0 1 10 20 30 40

w [rup
L (σ)]

10−5

10−4

10−3

10−2

P
(w

)

γsyn/σ

∞
5

1

0.2

Figure 2.2: Distribution of transverse sizes (along z) of the current layer for different cooling
regimes. The scale is zoomed in for w < rupL (σ). 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 ∼
rupL (σ).

29



this peak is associated with the X-points, where particles are accelerated to characteristic Lorentz

factors of γ ≲ σ (see Sec. 2.5.1). More importantly, current sheets with stronger cooling dis-

play a steeper slope and shorter distribution span for P(w), indicating that the plasmoids in the

current layer become smaller and more concentrated as cooling strength increases (more detailed

analysis of this issue was recently performed by Schoeffler et al. 2023).

Variation of the current sheet structure due to synchrotron cooling is reflected in one of

the most important macroscopic parameters of the reconnection – its rate, βrec. The evolution

of the rate, measured as the inflow velocity of the plasma, βrec = (E ×B)z /|B|2, averaged in

the thin layer far upstream parallel to the current sheet, for different cooling regimes is shown in

Fig. 2.3. As mentioned before, 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 (cf. Zhang et al. 2021a). Notably, in the simulations with

stronger cooling (γsyn ≲ σ), the reconnection rate is systematically larger, with values ranging

around βrec ≈ 0.12...0.14. We associate this correlation with the change in the filling fraction of

X-points in the current layer, which for the stronger cooled runs is approximately twice as large,

compared to the run without cooling as seen from the peak values in Fig. 2.2. This conclusion is

also consistent with the characteristically smaller plasmoid sizes in the strong cooling regime, as

discussed above.

Plasmoids accumulate the hot plasma energized in the X-points. In 2D the energetic par-

ticles 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 struc-

ture 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,

30



0 1 2 3 4

t [Lx/c]

0.00

0.05

0.10

0.15

0.20

β
re

c

initial
transient

γsyn/σ

∞
5

1

0.2

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

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 mo-

mentum from the thermal one by defining the fluid (Eckart) frame for species s, moving with

a velocity Uµ
s = Nµ

s /
√
Nν,sNν

s , where Nµ
s ≡

∫
fs(u

µ/u0)d3u, and fs and uµ are the distri-

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

P µν
s , is then computed as the projection of the stress-energy tensor: P µν

s = ∆µ
α∆

ν
βT

αβ
s , where

T µν
s ≡

∫
(pµpν/p0)fsd

3p, ∆µν
s ≡ ηµν+Uµ

s U
ν
s 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 non-

31



diagonal elements of the pressure tensor compared to diagonal elements are typically at the level

of 5% inside plasmoids in the weak-to-no cooling regime, with an increase to non-negligible

≈ 20% for the strongest cooling regime, γsyn/σ = 0.2. In addition, inside plasmoids, the diago-

nal component perpendicular to the magnetic field is suppressed by a factor of ∼ 5 in the case of

strong cooling.5 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 6, by removing the pressure

support, which results in extremely compressed structures (see, e.g., Hakobyan et al., 2019).

In 3D, the picture is quantitatively different. The force-balance condition is demonstrated in

Fig. 2.4, where in panels (a) and (c) we show the cross-section of the plasma density in x-z

plane for two simulations: one without cooling (3dx3kCoolInf), and the other with strong

cooling (3dx3kCool02). We pick intermediate-size plasmoids in both of these simulations

(indicated by the dotted vertical lines in left panels) and in panels (b) and (d) plot the magnetic

pressure, B2/8π, the scalar plasma pressure, defined as P = 1/3
∑

s P
i
s, i, as well as the z-

components of the two forces: the pressure gradient, c∇iP
iz, and the Lorentz force, (j ×B)z,

as functions of z. In the no-cooling simulation (Fig. 2.4a, b), large plasmoids are typically weakly

magnetized, while in the strong cooling regime (Fig. 2.4c, d), the strength of the magnetic field

inside the plasmoids is several times larger than the upstream value, Bup, (compare blue solid

lines in Fig. 2.4b, d). Despite this, even in the strongest cooling regime, the magnetic field

5The role of the pressure anisotropy due to synchrotron cooling is discussed by Zhdankin et al. (2022).
6This effect is 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.

32



compression inside the plasmoids is much smaller than what is observed in 2D (Hakobyan et al.,

2021; 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 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.

−50

0

50

z
/d
e

(a) no cooling

400 500 600 700 800 900

x/de

−50

0

50

z
/
d
e

(c) strong cooling

100 101 102
total density [nup]

0

1

2

3

4

5

−0.4

−0.2

0.0

0.2

0.4(b)

B2/8π

P

(j ×B)z

c∇iP iz

−50 0 50

z/de

0

1

2

3

4

5

p
re

ss
u

re
[B

2 u
p
/
8
π

]

−0.4

−0.2

0.0

0.2

0.4

force
[B

2u
p /

(8
π
d
e )]

(d)

Figure 2.4: The cross-section of typical plasmoids for uncooled (a, b) and strongly cooled (c,
d) current sheets. The left column (a, c) shows density slices and the magnetic field lines. The
region highlighted with a cyan-dotted boundary is shown in 1D cuts in (b, d), where we plot
the magnetic field, 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. 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.

33



The plasma inside flux tubes both in uncooled and strongly cooled simulations is com-

pressed and heated to a rough equipartition with the upstream magnetic field pressure: P ≈

B2
up/8π. In the latter case, the pressure is mainly delivered by higher plasma density (compare

densities in panels (a) and (c) in Fig. 2.4), since strong synchrotron cooling prevents plasmoids

from heating up by effectively radiating away the components of particle momenta perpendicular

to the magnetic field. Orange lines in Fig. 2.4b and d show the two force terms: the gradient of

the plasma pressure tensor, and the magnetic force, both plotted against the vertical axis z. The

balance between these two forces is well satisfied for the flux tube in the uncooled simulation

(Fig. 2.4b), and is reasonable in the strongly cooled simulation (Fig. 2.4d).

In this section, we largely neglected the presence of ions (insusceptible to synchrotron cool-

ing) 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 Fig. 2.4, due to their smaller number density. However, at later times ions can accelerate to

much higher Lorentz factors than the pairs (see Sec. 2.5.1), at which point their lack in the num-

ber density is compensated by the higher energy. When the ion pressure dominates in flux tubes,

we expect the structures of plasmoids to be similar to the case of no-cooling (Fig. 2.4a, b and

Appendix A.1). The applicability of this regime for the current sheets in pulsar magnetospheres

is studied in detail in Sec. 2.7.1.

2.5 Particle acceleration

During the non-linear stage of reconnection, cold particles from upstream – both pairs and

ions – are advected into the current sheet, where they are energized in either the X-points or the

34



relativistic outflows from the X-points (e.g., French et al., 2022). After this initial “injection”

stage most of these particles enter into the plasmoids. In 2D, strong compression of the magnetic

field inside plasmoids leads to adiabatically slow secondary acceleration of these particles to

high energies limited by the system size (Petropoulou & Sironi, 2018; Hakobyan et al., 2019).

However, in 3D, the finite size of plasmoids in the out-of-plane direction allows the most energetic

particles to freely escape from their bounds. Moreover, as described above, compression of

the magnetic field inside plasmoids is significantly smaller in 3D. Thus, the same secondary

acceleration mechanism as in 2D cannot operate efficiently in 3D. On the other hand, in 3D,

another acceleration channel has been observed for uncooled pair plasmas by Zhang et al. (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,

as the magnetic field vanishes in these regions. However, as soon as the energized pairs leave the

X-points, they lose their energy shortly after being exposed to a strong perpendicular magnetic

field component. As a result, the secondary acceleration channel in the reconnection upstream is

strongly suppressed for pairs. Ions, on the other hand, are unaffected by cooling, and can thus

tap the upstream electric field and get accelerated to potentially large Lorentz factors, ≫ σ. This

can clearly be seen in Fig. 2.5, where we plot the energy spectra of both pairs and ions. For the

simulations with strong-to-marginal cooling (γsyn/σ = 0.2...1), the distribution of pairs (solid

lines) typically spans up to ∼ σ, with the slope being close to γ−2 at highest energies. Ions

(dashed lines), on the other hand, regardless of the cooling strength, extend to Lorentz factors

significantly above σ, with the cutoff energy itself increasing in time. In the uncooled simulation

(red lines), the cutoff in the distribution of pairs is exactly the same as for the ions, suggesting

that the underlying secondary acceleration mechanisms are identical. Surprisingly, in simulations

35



0.1σ σ 10σ 100σ

γ

103

104

105

106

107

108

γ
d
n
/
d
γ

[a
rb

.
u

n
it

s]

f(γ) ∝ γ−1.7

γ−2

γsyn/σ

∞
5

1

0.5

0.2

de 10de 100de
rup
L (γ)

pairs

ions

Figure 2.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 ones for the uncooled runs, γsyn/σ = ∞. Dynamically strong cooling effectively
limits the maximum energy of pairs to ∼ σ, but at the same time facilitates the formation of
harder distribution for the ions. For reference, we also show the corresponding Larmor radii, rupL .

with stronger cooling, ions form a much harder power law distribution, with the typical slope

close to γ−1 for the strongest cooling case, and steepening to γ−1.7 for the uncooled case (this

effect is explored in Sec 2.5.1).

We further quantify the contributions of different acceleration channels in Fig. 2.6, where

we plot the energy spectra of both ions and pairs in the uncooled and strongly cooled simulations.

To test the dominant acceleration mechanism up to γ ∼ σ (the injection stage), we choose the

particles that experienced the non-ideal electric field, E > B, during their lifetime (dashed lines

in Fig. 2.6). In the inset panels of Fig. 2.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

36



0.1σ σ 10σ 100σ

γ

103

104

105

106

107

108

γ
d
n
/d
γ

[a
rb

.
u

n
it

s]

γsyn

pairs ions

0.1σ σ 10σ10σ 100σ

∆γ

100

104

108
f

(∆
γ

)

σ/4

103

104

105

106

107

108

γ
d
n
/d
γ

[a
rb

.
u

n
it

s]

particles that experienced E > B

0.1σ σ 10σ10σ 100σ

∆γ

100

104

108

f
(∆
γ

)

σ/4

acceleration in E > B region

(b) strong cooling

(a) no cooling

Figure 2.6: Distribution functions of the pairs (blue) and the ions (orange) for both the uncooled
(a) and the strongly cooled (b) simulations. Dashed lines show parts of the distributions con-
tributed by the particles having experienced E > B during their lifetime in the sheet. Inset
panels show the distribution of energy gains, ∆γ, of particles during their contiguous presence
in these regions. The non-ideal electric field in the X-points only accelerates particles to Lorentz
factors of a fraction of σ. At the same time, this process is essential to promote the particles
(pairs and ions) for further acceleration up to the highest energies: in relativistic outflows along
the layer and the large-scale reconnection electric field upstream. As a result, the high energy tail
of the spectrum is fully formed by the particles passed through E > B.

37



region dominate at energies ≳ σ (cf. Sironi, 2022). However, the net energy gain of particles

in X-points in both the uncooled and strongly cooled cases does not exceed a fraction of σ, and

thus cannot account for the high-energy power-law tail of either the pairs or the ions (e.g., Werner

et al., 2016; Uzdensky, 2022). The main energy gain, thus, occurs due to the acceleration by ideal

electric fields. For pairs, as we argue further in Sec. 2.5.2, “pick-up” acceleration by outflows

from the X-points allows them to gain energies up to σ even in the strong cooling case.

2.5.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 Fig. 2.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 Fig. 2.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 ex-

periences most of its acceleration while moving outside of the current layer (Zhang et al., 2021a,

2023), as evident from the thickness of the line representing the trajectory in Fig. 2.7. These

trajectories correspond to bouncing of accelerated ions between the two converging upstream

38



flows7. In contrast, the low-energy ion displays a tangled trajectory due to its entrapment inside a

plasmoid or frequent scattering on magnetic field inhomogeneities inside the current sheet. This

behavior is likely associated with self-generated turbulence (e.g., Guo et al., 2021).

0

200

400

600

800

1000

y
/d
e

Erec

Bup

particle location

upstream

current sheet

400 500 600 700 800 900 1000 1100 1200

x/de

−50

0

50

z
/d
e

Bup

Erec

100 200 300 400 500 600 700 800 900

x/de

0 1 2

ct/Lx

σ

10σ

20σ

30σ

γ

0 1 2

ct/Lx

σ

10σ

20σ

30σ

γ

γ̇ = (qi/mic
2)Erecvy

10de

100de

200de

rup
L (γ)

σ

10σ

20σ

30σ

γ

Figure 2.7: Archetypal trajectories (projected to the x-y and the x-z plane) and the energy evo-
lution of two ions. The color of each dot corresponds to the Lorentz factor of the ion at that
particular time, while the size of the dot indicates whether the ion is inside the current sheet
(small dots), or upstream (large dots). The ion on the right accelerates to high energies (≫ σ) by
the large-scale electric field, Erec, in the region upstream of the layer, while the one on the left
gains only a marginal amount of energy (∼ σ). The more energetic ion on the right primarily
moves upstream of the current layer along the reconnection electric field (y direction) and accel-
erates linearly in time. The slower particle, on the other hand, is trapped inside the current sheet
and experiences constant scatterings on local inhomogeneities.

The acceleration of the more energetic ion takes place linearly in time, with its trajectory
7Note that the classical Speiser orbits correspond to the acceleration inside the current sheet, as opposed to the

acceleration upstream (for more discussion, see Zhang et al. 2021a).

39



showing almost no scatterings. Notably, the preferred direction of motion for this ion is in y,

while the acceleration rate is almost constant. Both of these facts point towards the large-scale

reconnection-driven electric field, Erec ≈ Ey, being the dominant acceleration mechanism. The

rate of acceleration for these trajectories can be written as

γ̇ =
qi
mic

(
Erec ·

v

c

)
=⇒ ∆γ ≈ βrecβy

qiBup

mic
∆t, (2.7)

where we used Erec = βrecBup (Zhang et al., 2021a). This slope is overplotted with a dashed

line in the inset of the left panel of Fig. 2.7, showing excellent agreement with the acceleration

history. Since ions are by construction not susceptible to synchrotron cooling, we find that this

acceleration mechanism operates similarly for both the strong cooling and no-cooling regimes.

Despite the ion acceleration mechanism being identical for different cooling regimes, the

distribution of ions is much harder in the stronger cooling cases, as highlighted with dashed

lines in Fig. 2.5. Two parameters from Eq. (2.7) that could potentially contribute to different

acceleration rates are the reconnection rate, βrec,8 and the beaming of ions towards the direction

of the reconnection electric field, βy. To inspect the relative importance of the two effects, in

Fig. 2.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 – 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

8As it was already mentioned in Sec. 2.4, 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.

40



σ

10σ

20σ

30σ

40σ

γ

(a1)

γin = σ

(a2)

γin = 5σ

(a3)

γ̇
=
q i
E re

c
/m
ic

γin = σ

(a4)

γin = 5σ

−1.0

−0.5

0.0

0.5

1.0

β
y

(b1) (b2) (b3) (b4)

0 1 2

time [Lx/c]

−0.1

0.0

0.1

E
y
/B

u
p

(c1)

0 1 2

time [Lx/c]

(c2)

0 1 2

time [Lx/c]

(c3)

0 1 2

time [Lx/c]

(c4)
〈Ey〉

101

102

103

101

102

103

101

102

103

no cooling strong cooling

Figure 2.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 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 fieldEy 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 en-
ergy evolution, a3-a4). In the uncooled regime, the energy gain is more chaotic as ions constantly
scatter inside the current sheet.

41



the even-numbered ones, we pick ions with γin ≈ 5σ (the middle of the tail of the distribution).

If we consider the population of particles with energies γ ≈ σ (consider a slice at t = 0

in Fig. 2.8b1, b3), we see that the motion of ions is either not beamed at all in the y-direction

(the case with no cooling), or has only marginal beaming (the case with strong cooling). By

this time, the energy gain of these ions, γin ≈ σ, is dominated by the Fermi acceleration in the

current sheet (either via a Fermi kick or the “pick-up” mechanism). This prelude is consistent

with the fact that the ions at these energies are primarily moving in the plane of the reconnecting

magnetic field, x-z, i.e., perpendicular to the reconnection electric field, βy ≈ 0, (cf. French

et al., 2022). At a later time, a fraction of these particles begin accelerating upstream of the layer

along the y direction, i.e., along Erec; that fraction of “lucky” particles ultimately determines

the power-law slope in the distribution. As it is evident from the first column of Fig. 2.8, in

the uncooled run only a small fraction of ions follows the γ̇ ∝ Erec · v linear acceleration path,

with the bulk beaming in y-direction remaining marginal. This is due to the fact, that the cross-

section of typical plasmoids in the simulation with weak-to-no cooling is significantly larger (see

Fig. 2.2), making them more capable of trapping the freely accelerating ions. On the other hand,

in the strong cooling case almost all of the particles that reach γin ≈ σ are “lucky” to land on the

efficiently accelerating trajectories, beamed in the direction of the electric field in the upstream

(see the third column in Fig. 2.8; these trajectories are similar to the one shown in the right

column of Fig. 2.7). As the “lucky” particles become more energized, they are essentially bound

to remain on these accelerating trajectories. This can be clearly seen in the second and the fourth

panels of Fig. 2.8, where we show the same tracking statistics, but only for the particles that have

already reached γin ≈ 5σ. For this selection, the fraction of “lucky” particles is substantially

larger. This suggests that even in the uncooled case the ion energization mechanism to γ ≫ σ is

42



the same large-scale acceleration in the upstream, albeit a smaller amount of accelerated particles,

free-streaming in the upstream, results in a steeper power-law slope in their distribution.

We observe a general tendency that the acceleration rate for individual “lucky” ions is

larger for the strongly cooled simulations. This can be clearly seen from the typical strength

of the electric field experienced by the ions in Fig. 2.8c1...c4: for stronger cooling the electric

field is systematically larger than in the uncooled case, even if we consider particles beyond

γin ≈ 5σ. We associate this difference with the fact that the global reconnection rate, βrec, is

larger for the strong cooling case, as shown in Fig. 2.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 are moving primarily in

the upstream region (cf. Fig. 2.7), where the electric field is coherent and almost exclusively

points perpendicular to the magnetic field direction y.

2.5.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 non-ideal and ideal

electric fields up to Lorentz factors comparable to σ, after which a fraction of them escape to

the upstream where they are linearly accelerated further by the large-scale 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

43



perpendicular magnetic field, B̃, defined in (2.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 Fig. 2.5, this effect can clearly be seen by comparing the uncooled

simulation with the one with γsyn = 5σ. While in the uncooled case, the distributions of both pairs

and ions extend to several tens of σ in energy, in the weakly cooled simulation the distribution of

pairs essentially cuts off at ≈ γsyn = 5σ, with ions being relatively unaffected.

In the strong cooling regime, the upstream acceleration is essentially prohibited, as the

strong synchrotron cooling disallows particles to leave the X-points by crossing magnetic field

lines. As shown in Fig. 2.5, the spectrum of accelerated leptons does not extend beyond ≈ σ.

In this regime, two acceleration mechanisms are particularly efficient: the X-point accelera-

tion by the non-ideal electric field, and the “pick-up” mechanism when particle is injected non-

adiabatically into a high velocity reconnection exhaust. In Fig. 2.6b, we show that particle energy

gain in the E > B zones is not sufficient to extend the power-law to γ ≈ σ; this implies that pairs

with the highest energies are additionally accelerated by the “pick-up” mechanism.

An example of a typical high-energy electron trajectory in the strong cooling regime is

shown in Fig. 2.9 (3dx2kCool02). In panels (a)...(e) we plot the x-z slices of the plasma den-

sity, and the magnetic field lines at different timesteps, with the current positions of the particle

indicated with colored circles. The instantaneous four velocities of the particle are indicated with

colored arrows. We also over-plot in blue the magnetic field line at the position of the particle,

with the width of the line indicating its strength. On the right of Fig. 2.9 we show the time

evolution of the Lorentz factor γ, with the dashed lines indicating the times of the snapshots on

the left. In panel (a) the particle is upstream from the layer, moving towards the current sheet

with a small Lorentz factor, γ ∼ 1. In panel (b), the particle enters the X-point and is acceler-

44



350 400 450 500 550 600

−50

0

50

z
/
d
e

(a) y/de = 286

100 101 102
total density [nup]

350 400 450 500 550 600

−50

0

50

z
/
d
e

(b) y/de = 302

0.1 2.5
B/Bup

300 350 400 450 500 550

−50

0

50

z
/
d
e

(c) y/de = 330

1
30

60

γ

150 200 250 300 350 400

−50

0

50

z
/
d
e

(d) y/de = 366

100 150 200 250 300 350

x/de

−50

0

50

z
/
d
e

(e) y/de = 400

0 3de 6de 9de

rup
L (γ)

0 25 50
γ

0

1

2

3

4

5

tim
e

[10
3
m
e c/|e|B

u
p ]

σ

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

45



ated by the non-ideal electric field, E > B, to γ ≈ 0.5σ. After this rapid acceleration event,

the particle moves along with the relativistic outflow from the X-point, constantly crossing the

upstream/downstream boundary, which is perturbed vertically due to RDKI (see below). As the

particle hits the upstream, it experiences subtle energy losses indicated as oscillations in its en-

ergy, γ. After that, the particle is picked up by a contracting magnetic field line performing a

slingshot-like motion. After that, as shown in panel (c), for an extended period of time particle

moves along with a rapidly receding field line while being slowly accelerated to γ ≳ σ. Once

the slingshot-like motion of the field line comes to a halt, the particle is injected into the forming

plasmoid, and experiences a catastrophic cooling event, being exposed to a large B̃⊥ (panel (d);

this event is studied in more details later). During its interaction with the plasmoid, the parti-

cle radiates most of its momentum and is then trapped inside the plasmoid until the end of the

simulation, with a low Lorentz factor, γ ∼ few, as shown in panel (e).

We further describe the 3D motion of the same accelerated electron in Fig. 2.10a...d, where

we show 3D volume renderings of the plasma density and the magnetic field lines in the vicinity

of the electron, also shown in Fig. 2.9. Panels on the right display the energization history of

the particle, as well as the interesting quantities zoomed at two distinct moments of the particle’s

evolution (roughly corresponding to panels (a)/(b), and (c)/(d) of the same figure). In the upper

half of Fig. 2.10, we specifically focus on the time when the particle experiences successive

scattering events while moving in the RDKI-perturbed current layer, after being pre-energized in

the X-point (corresponding to Fig. 2.9b). As indicated in the upper right panels, just before the

scattering event, the particle crosses the null line, |B| ≈ 0, and approaches the peak of the RDKI

perturbation, where both B̃⊥ and |B| grow rapidly (time t ≈ 1.8 in units indicated on the plot;

field strengths at particle position are shown with blue and orange lines). B̃⊥ reaches the critical

46



(a)

(b)

(c)

(d)

0 2 4

time [103 mec/|e|Bup]

0

20

40

60

γ

(a) (b)

(c) (d)

0

5

10

R
c
/
ru

p
L

(σ
)

−1

0

p
ow

er
[|e
|B

u
p
c]

eE · v
fdrag · v
γ̇mec2

1.7 1.8 1.9
0

1

2

B
-fi

el
d

[B
u
p
]

|B|
γsyn/γ

B̃⊥

0

1

2

B
-fi

el
d

[B
u
p
]

|B|
γsyn/γ

B̃⊥

−1

0

p
ow

er
[|e
|B

u
p
c]

eE · v
fdrag · v
γ̇mec2

3.4 3.6
0

5

10

R
c
/r

u
p
L

(σ
)

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

47



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 Fig. 2.10 (panels c,d; see also

Fig. 2.9d), 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 Fig. 2.10c, and Fig. 2.9c, and d), the effective perpendicular

magnetic field, B̃⊥ is almost zero, allowing the particle to avoid being radiatively cooled. At

the 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 of Fig. 2.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 Fig. 2.9e).

In both the scattering and the catastrophic cooling examples, the full radiation reaction

force, given by Eq. 2.3, consists of both synchrotron and curvature contributions. As shown in

the top and the bottom panels on the right in Fig. 2.10, in both cases the local radius of curvature

of the field line9 during the deceleration is comparable to the characteristic Larmor radius of the

particle with the Lorentz factor γ ≈ σ, i.e., it is close to the gyro-radius of the particle. Under

these conditions, the synchrotron and curvature radiation mechanisms are indistinguishable, with

the radius of curvature of field lines being microscopic and set by the local field perturbations

inside the current layer, i.e., the curvatu