SIMULATIONS OF PARTICLE ACCELERATION BEYOND THE CLASSICAL SYNCHROTRON BURNOFF LIMIT IN MAGNETIC RECONNECTION: AN EXPLANATION OF THE CRAB FLARES

The maximum energy reached by a charged particle in a given astrophysical object is limited by the size of the acceleration region (Hillas 1984). If the relativistic Larmor radius of the particle R is of order the system size L, the particle escapes and is no longer accelerated. The maximum energy is then given by $\mathcal {E}_{\rm max}\lesssim qBL$, where q is the charge of the particle and B is the typical magnetic field strength. Radiative losses within the accelerator decrease this limit (see, e.g., Aharonian et al. 2002; Medvedev 2003). The maximum energy is then set by the balance between the electric acceleration rate and the radiative power lost by the particle. In the case of synchrotron cooling, this balance leads to a remarkable result: the maximum synchrotron photon energy emitted by an electron depends only on the ratio of the electric field to magnetic field perpendicular to the particle's motion, i.e., $\epsilon ^{\rm max}_{\rm sync}=(9mc^2/4\alpha _{\rm F})(E/B_{\perp })$, where α_F is the fine structure constant and mc² is the rest mass energy of the electron. Hence, under ideal magnetohydrodynamic (MHD) conditions where E ⩽ B, the energy of synchrotron radiation should not exceed the fundamental constant 9mc²/4α_F ≈ 160 MeV (Guilbert et al. 1983; de Jager et al. 1996; Lyutikov 2010; Uzdensky et al. 2011). An electron with energy above the radiation reaction limit would lose most of its energy in a fraction of a Larmor gyration. It is then impossible to have electrons radiating synchrotron radiation above 160 MeV with classical models of particle acceleration, all based on ideal MHD (e.g., diffuse shock acceleration), unless the plasma has a relativistic bulk motion with respect to the observer, or the electrons are the by-product of energetic particle decay.

Yet, the gamma-ray space telescopes AGILE and Fermi discovered several powerful gamma-ray flares from the Crab Nebula (Tavani et al. 2011; Abdo et al. 2011; Balbo et al. 2011; Striani et al. 2011, 2013a, 2013b; Buehler et al. 2012; Ojha et al. 2013), presumably during which PeV electrons and positrons emit synchrotron radiation well above the 160 MeV limit. The most powerful flare, recorded in 2011 April, showed clear evidence for synchrotron emission up to 375 MeV (Buehler et al. 2012). This discovery suggests that an extreme and non-conventional particle acceleration mechanism is at work somewhere in the nebula, unless the emission is substantially Doppler-boosted by a factor ≳ 2. However, the typical flow velocity in the nebula (about half the speed of light, e.g., Hester et al. 2002) and its orientation with respect to the observer give a Doppler factor of order unity. The precise location of the flare is unknown because the nebula is not resolved in gamma rays, but the lack of pulsations suggests it does not originate very close to the pulsar. The ≲ 8 hr flux-doubling timescale observed during the flare indicates that a tiny fraction of the nebula is involved. The flaring region radiates about 30 times more than the quiescent emission above 100 MeV, which represents up to 1% of the pulsar spin-down power (Buehler et al. 2012). To explain particle acceleration and emission within the overall duration of the flares, ranging from a few days to a few weeks, the magnetic field should be of order a few milligauss, i.e., much more intense than the average ∼200 μG traditionally inferred from spectral modeling (Horns & Aharonian 2004; Meyer et al. 2010). Simultaneous observations in radio, near-infrared, optical, X-rays, and in TeV gamma-rays were not able to detect a solid counterpart to the flares (see Weisskopf et al. 2013, and references therein), suggesting that the emitting particle spectrum is very hard, perhaps monoenergetic. It is very difficult to reconcile these puzzling features with classical models of particle acceleration and models of pulsar wind nebulae (Rees & Gunn 1974; Kennel & Coroniti 1984; see Kirk et al. 2009; Arons 2012 for recent reviews).

Various models have been proposed to solve the Crab flare mystery. Several studies invoke a relativistic Doppler boosting of the flaring region by a factor of a few. It was proposed that a mildly relativistic flow could be achieved close to the pulsar wind termination shock (Komissarov & Lyutikov 2011; Lyutikov et al. 2012; Bednarek & Idec 2011), in relativistic reconnection events within the nebula (Clausen-Brown & Lyutikov 2012), in a magnetized flow at the base of the Crab jets (Lyubarsky 2012), or in knots of energetic particles (Yuan et al. 2011). The dissipation of the striped pulsar wind structure through the shock (Pétri & Lyubarsky 2007; Sironi & Spitkovsky 2011) could also generate rapidly fluctuating magnetic field on small scales, resulting in synchrotron gamma-ray flares (Bykov et al. 2012). If the magnetic turbulence in the nebula occurs on a length scale shorter than the synchrotron photon formation length mc²/eB, then the particles could emit jitter radiation (Medvedev 2000), with typical energy greater than the classical synchrotron limit (Teraki & Takahara 2013). Alternatively, particle acceleration could occur in regions of strong coherent electric field, in twisted toroidal fields (Sturrock & Aschwanden 2012), or in magnetic reconnection sites within the nebula (Uzdensky et al. 2011; Cerutti et al. 2012a).

Magnetic reconnection offers natural locations (within the diffusion region where the magnetic field is small and reverses) in which the electric field can exceed the magnetic field. In principle, it is possible to accelerate particles above the classical radiation reaction limit at these sites (Kirk 2004). Uzdensky et al. (2011) demonstrated that the highest energy particles are trapped and focused toward the reconnection layer midplane⁴ (see also Contopoulos 2007), following the relativistic analog of Speiser orbits (Speiser 1965; see Figure 1). Once deep inside the layer, the particles are subject to weak radiative losses, but strong coherent electric field. The layer acts as a linear accelerator, and the maximum energy of the particles is then limited just by the total electric potential drop along the layer, i.e., $\mathcal {E}_{\rm max}\sim eEL$, where L is the length of the layer. Using relativistic test-particle simulations, we found in Cerutti et al. (2012a) that magnetic reconnection could generate a quasi-monoenergetic beam of particles above the radiation reaction limit in the Crab Nebula.

Zoom In Zoom Out Reset image size

In this paper, we re-examine particle acceleration in ultra-relativistic pair plasma reconnection, using two-dimensional (2D) particle-in-cell (PIC) simulations. Such simulations capture self-consistently the time evolution of fields and particles at the kinetic level. For this study, we developed a new relativistic PIC code, called Zeltron, that includes the radiation reaction force on the particles. Our first objective is to investigate ab initio whether reconnection can accelerate particles above the radiation reaction limit, under realistic physical conditions. Our second objective is to study the radiative signature of such acceleration following Cerutti et al. (2012b), in order to explain all observed features of the Crab flares. Cerutti et al. (2012b) included optically thin synchrotron radiation as a tracer, but this calculation was not self-consistent because it did not treat radiation reaction effects on particle motion. Section 2 describes the main capabilities of Zeltron, with an emphasis on the radiation reaction force. Section 3 gives the initial setup of the simulations performed in this study, and Section 4 presents the main results on particle acceleration and radiation. We discuss our findings in the context of the Crab gamma-ray flares in Section 5. Section 6 summarizes the main results of this work.

Zeltron is a new 3D, parallel (domain decomposition with MPI), relativistic, electromagnetic PIC code developed independently from other existing codes (for a review about PIC methods, see, e.g., Pritchett 2003; Birdsall & Langdon 2005). Zeltron follows an explicit finite-difference scheme on a Cartesian grid, with a time step Δt set at a fraction of the maximum stable step determined by the Courant–Friedrichs–Lewy condition, i.e., cΔt ⩽ cΔt_CFL = (1/Δx² + 1/Δy² + 1/Δz²)^−1/2, where c is the speed of light and Δx, Δy, and Δz are the minimum grid spacing in the x-, y-, and z-directions. The code uses the Yee algorithm (Yee 1966) to solve the time-dependent Maxwell's equations, given by

$\begin{eqnarray} \frac{\partial \mathbf {E}}{\partial t} &=& c\mathbf {\nabla }\times \mathbf {B}-4\pi \mathbf {J} \end{eqnarray} \tag{ 1 }$

$\begin{eqnarray} \frac{\partial \mathbf {B}}{\partial t} &=& {-}c\mathbf {\nabla }\times \mathbf {E}, \end{eqnarray} \tag{ 2 }$

where E is the electric field, B is the magnetic field, and J is the current density. This algorithm has second-order error in space and time and ensures that ∇ · B = 0 at any instant of the simulation (to the computer round-off accuracy). The numerical scheme does not, however, satisfy the Maxwell-Gauss equation ∇ · E = 4πρ exactly, where ρ is the charge density. The electric field has to be corrected every time step by the small amount δE, obtained by solving Poisson's equation (Pritchett 2003), i.e.,

$\begin{equation} \nabla ^2\left(\delta \phi \right)=-\left(4\pi \rho -\mathbf {\nabla }\cdot \mathbf {E}\right), \end{equation} \tag{ 3 }$

and δE = −∇δϕ. The Poisson solver implemented in the code utilizes an iterative Gauss–Seidel method (with five points in 2D, and seven points in 3D).

The main novelty of Zeltron, compared to most PIC codes,⁵ is the ability to take into account the effect of the radiation reaction force on the motion of the particles. The (non-covariant) equation of motion of a single particle is given by the Lorentz–Abraham–Dirac equation (Landau & Lifshitz 1975)

$\begin{equation} mc\frac{d\mathbf {u}}{dt}=q\left(\mathbf {E_{\rm i}}+\frac{\mathbf {u}\times \mathbf {B_{\rm i}}}{\gamma }\right)+\mathbf {g}, \end{equation} \tag{ 4 }$

where u = γv/c is the four-velocity of the particle, γ = (1 − v²/c²)^−1/2 is the associated Lorentz factor, and g is the radiation reaction force. The fields at the location of the particle E_i and B_i, are linear interpolations of the fields E and B known at the grid nodes. Within the framework of classical electrodynamics, the radiation reaction force is obtained from the Landau–Lifshitz equation (Landau & Lifshitz 1975), valid as long as the product γB is much smaller than the quantum critical magnetic field B_QED = m²c³/eℏ = 4.4 × 10¹³ G, where e is the elementary electric charge, and ℏ = h/2π with h the Planck constant. In the ultra-relativistic regime (γ ≫ 1), the radiation reaction force can simply be expressed as a continuous friction force opposite to the particle's direction of motion (see discussion in Cerutti et al. 2012a), i.e.,

$\begin{equation} \mathbf {g} = -\frac{2}{3}r_{\rm e}^2\gamma \left[\left(\mathbf {E_{\rm i}}+\frac{\mathbf {u}\times \mathbf {B_{\rm i}}}{\gamma }\right)^2-\left(\frac{\mathbf {u}\cdot \mathbf {E_{\rm i}}}{\gamma }\right)^2\right]\mathbf {u}, \end{equation} \tag{ 5 }$

where r_e = e²/mc² is the classical radius of the electron. Following Tamburini et al. (2010), Zeltron uses a modified Boris algorithm (Birdsall & Langdon 2005) to solve the equation of motion with the radiation reaction force given in Equation (5). For an alternative implementation of the radiation reaction force in PIC codes, see also Sokolov et al. (2009) and Capdessus et al. (2012).

The code uses linear interpolation to deposit the charges and currents generated by each particle at the nodes of the computational grid, and computes the charge and current densities for Maxwell's equations. The code assigns variable weights to the macro-particles to model particle density gradients. Zeltron does not strictly conserve the total energy. For the purpose of this study, we ran Zeltron successfully on thousands of cores on the Kraken supercomputer⁶ and on the University of Colorado Janus and Verus supercomputers. We find excellent agreement between Zeltron and the well-tested PIC code Vorpal (Nieter & Cary 2004), in the limit where radiative losses are neglected (the radiation reaction force is currently not implemented in Vorpal).

We present a series of simulations of ultra-relativistic electron–positron pair plasma reconnection with radiation reaction force in 2.5D (the fields depend on two coordinates, but the motion of particles is 3D), using Zeltron. The initial setup is very similar to our previous study (Cerutti et al. 2012b), which is standard in such reconnection simulations (Zenitani & Hoshino 2001, 2007, 2008; Jaroschek et al. 2004a, 2004b; Bessho & Bhattacharjee 2007, 2012; Daughton & Karimabadi 2007; Pétri & Lyubarsky 2007; Lyubarsky & Liverts 2008; Yin et al. 2008; Jaroschek & Hoshino 2009; Liu et al. 2011; Sironi & Spitkovsky 2011; Kagan et al. 2012). The computational domain is a rectangle of size L_x × L_y, where L_{x, y} is the length of the box in each direction. It contains two anti-parallel relativistic, flat, Harris current sheets (Kirk & Skjæraasen 2003) in the xz-plane, which enables us to set periodic boundary conditions in the x- and y-directions. The reconnecting magnetic field, B_x, reverses across each layer along the y-direction, and is given by

$\begin{equation} B_{\rm x} = \left\lbrace \begin{array}{@{}lcl}-B_0\tanh \left(\frac{y-L_{\rm y}/4}{\delta }\right) &{/rm if} & y<L_{\rm y}/2\\ \ms B_0\tanh \left(\frac{y-3L_{\rm y}/4}{\delta }\right) & {\rm if} & y>L_{\rm y}/2\end{array} \right. \!\!, \end{equation} \tag{ 6 }$

where B₀ is the initial upstream reconnecting field strength and δ is the initial layer half-thickness (Figure 1). The strength of the uniform guide field component (Bz) is a free parameter that varies from 0 to B₀ in our simulations (see Section 4.5). There is no electric field initially, E = 0. We apply an initial perturbation to the magnetic field (10% maximum amplitude in the magnetic flux function; see Figure 2), in order to force the onset of reconnection at the beginning of the simulation, thereby reducing the computing time. If there is no perturbation, we find that reconnection eventually happens spontaneously in the simulation because of the high numerical noise inherent to PIC codes.

Zoom In Zoom Out Reset image size

A population of relativistic thermal pairs, following a Maxwell–Jüttner distribution and concentrated in the current layers, balances the upstream magnetic pressure and carries the initial current in the ±z-directions given by the bulk motion of these particles, drifting at a velocity v_drift/c = β_drift = 0.6. The density profile of the drifting particles is

$\begin{equation} n_{\rm drift} = \left\lbrace \begin{array}{@{}lcl}n_0\left[\cosh \left(\frac{y-L_{\rm y}/4}{\delta }\right)\right]^{-2} &{\rm if} & y<L_{\rm y}/2\\ \ms n_0\left[\cosh \left(\frac{y-3L_{\rm y}/4}{\delta }\right)\right]^{-2} & {\rm if} & y>L_{\rm y}/2\end{array} \right. \!\!, \end{equation} \tag{ 7 }$

with

$\begin{equation} n_0 = \frac{k T_{\rm drift}\left(1-\beta _{\rm drift}^2\right)^{1/2}}{4\pi e^2 \beta _{\rm drift}^2\delta ^2}, \end{equation} \tag{ 8 }$

where k is the Boltzmann constant and T_drift is the temperature of the drifting particles defined in their co-moving frame. In the simulation, the drifting particles are injected uniformly throughout the box, but with variable weight according to Equation (7). The purpose of this procedure is to represent low-density distributions with less noise. In addition, we fill the box with a uniform density n_bg = 0.1n₀ of isotropic non-thermal ultra-relativistic particles (see explanation below), such that

$\begin{equation} \frac{d{\rm N}_{\rm bg}}{d\gamma } = \left\lbrace \begin{array}{@{}lcl}K\gamma ^{-p} &{\rm if} & \gamma _1<\gamma <\gamma _2\\ 0 & {\rm otherwise} &\end{array} \right.\!\!, \end{equation} \tag{ 9 }$

where p is the index of the power law and K is a normalization constant.

We initialize the simulations with physical parameters consistent with the conditions thought to exist in the flaring region of the Crab Nebula, although our results are general and scalable. The energy and spatial scales of the problem are obtained from the reconnecting field strength B₀, of order a few milligauss during Crab flares. Following Cerutti et al. (2012a), we set B₀ = 5 mG. The fiducial radiation-reaction-limited energy of the particles, given by the balance of the electric force with the radiation reaction force for E = B₀, is then equal to

$\begin{equation} \gamma _{\rm rad}=\left(\frac{3 e}{2 r_{\rm e}^2 B_0}\right)^{1/2}\approx 1.3\times 10^9 B_{\rm 5\,mG}^{-1/2}, \end{equation} \tag{ 10 }$

where B_5 mG = B₀/5 mG, allowing electrons with energies up to in the PeV range. The Larmor radius associated with such a particle is R_rad = γ_radmc²/eB₀ ≈ 4.6 × 10¹⁴ cm, or 4.25 light-hours. The maximum energy expected to be reached by the particles in the simulation (in the absence of radiative losses) is limited by the electric potential drop along the layer, which is proportional to the system size L and the dimensionless reconnection rate β_rec < 1, such that γ_max ≲ β_receB₀L/mc². To observe particle acceleration above the radiation reaction limit (Equation (10)), the system size L must be bigger than R_rad. R_rad must also be large compared with the smallest spatial scales in the system, given by the Larmor radius of the lowest energy particles injected initially. We set the minimum Lorentz factor of the background particles to γ₁ = 4 × 10⁷ with a power-law index of p = 2 extending up to γ₂ = 4 × 10⁸ < γ_rad, and an ultra-relativistic temperature kT_drift/mc² = 4 × 10⁷ for the drifting particles. In reality, the typical energy of the background particles in the Crab Nebula is probably much lower, of order γ₁ ∼ 10⁴–10⁶. However, such a broad range in energy cannot be reached with the current computing power available, even for 2D simulations. Our simulations model only the high-energy end of the Crab Nebula spectrum. The initial −2 power law of the electrons approximately represents the purely non-thermal background particles emitting the quiescent emission of the Crab Nebula, inferred from spectral modeling (Horns & Aharonian 2004; Meyer et al. 2010). This assumes pre-acceleration of the particles in the nebula, possibly by shock acceleration or by other reconnection events.

In the following, we express all spatial scales in terms of the initial minimum Larmor radius ρ₁ = γ₁mc²/eB₀ ≈ 1.4 × 10¹³ cm, and timescales with respect to the corresponding inverse Larmor frequency $\omega _1^{-1}=\rho _1/c\approx 455$ s. We perform all the simulations with L_x = L_y = 500ρ₁, which corresponds to 7 × 10¹⁵ cm or 2.7 lt-day. For a typical value of β_rec = 0.1–0.2, we have γ_max ≲ 50–100γ₁ ≈ (2–4) × 10⁹, allowing us to see particles with gamma above γ_rad ≈ 1.3 × 10⁹. The other relevant quantities of the plasmas are the electron skin-depth in the layer d_e = (kT_drift/4πn₀e²)^1/2 ≈ 1.8ρ₁, the initial layer thickness $\delta =2kT_{\rm drift}(1-\beta _{\rm drift}^2)^{1/2}/\beta _{\rm drift}eB_0\approx 2.7\rho _1$, and the upstream magnetization parameter $\sigma =B_0^2/4\pi n_{\rm bg} \gamma _1 mc^2\approx 16$. The computational grid is Cartesian and uniform, composed of 1440² cells, giving a spatial resolution of ≈3 cells per ρ₁ or ≈8 cells per δ. We note that the simulation becomes unstable at late times if Δt = Δt_CFL, only when radiative cooling is strong and if there is no guide field. In particular, the electric field oscillates between adjacent time steps, leading to an artificial increase of radiative losses. We find that setting the time step to Δt = 0.3Δt_CFL is good enough to quench the development of this numerical instability. The time resolution is then $\Delta t\approx 0.07\omega _1^{-1}$. The total initial number of particles per cell is 100 (25 electrons and 25 positrons for each population). Table 1 gives all the physical and numerical parameters and their values chosen for this study.

Table 1. Physical and Numerical Parameters Common to All the Simulations Reported in This Work

Physical Parameters	Set Values
B0	5 mG
γrad	1.3 × 109
ρ1	1.4 × 1013 cm
de/ρ1	1.8
δ/ρ1	2.7
$\omega _1^{-1}$	455 s
kTdrift/mc2	4 × 107
βdrift	0.6
nbg/n0	0.1
σ	16
γ1	4 × 107
γ2	4 × 108
p	2
Numerical Parameters	Set Values
ρ1/Δx	3
ρ1/Δy	3
Δt × ω1	0.07
Particles/cell	100
Lx/ρ1	500
Ly/ρ1	500

Download table as:  ASCIITypeset image

This section presents the results based on a set of simulations of size (L_x × L_y) = (500ρ₁)² to investigate the general properties of reconnection under strong synchrotron cooling conditions and particle acceleration above the radiation reaction limit (Sections 4.1–4.3); we compare with an identical simulation without the radiation reaction force. In Section 4.4, we investigate the variability patterns (spectra, light curves, power-spectra) of the high-energy radiation emitted by the layer, and we show that these features are robust and reproducible, using 10 identical simulations with the radiation reaction force. The statistical variations in this sample of simulations originate only from the initial random positions and velocities of the particles in the box. We present also a set of four simulations to study the effect of the guide field on particle acceleration above the radiation reaction limit (Section 4.5).

4.1. Time Evolution of Reconnection

At the beginning of the simulation, the system remains quasi-static⁷ for about tω₁ ≲ 110, before the layers become unstable to tearing modes and break up into several islands of closed magnetic field loops filled with plasma ("magnetic islands" or "plasmoids"; see Figure 2). This instability produces multiple secondary X-points between the plasmoids, where the magnetic field reconnects and where the electric field, mostly along the ±z-directions (E ≈ E_z), is intense, of order B₀, leading to strong particle acceleration (see Figure 3, bottom-left panel). In addition, the formation of islands induces a strong bipolar reconnected magnetic field B_y ∼ B₀ concentrated at the edges of islands between two X-points (see Figure 3, middle-right panel). The magnetic tension of this field drives the reconnection outflow by pushing the plasma away from X-points into the direction of magnetic islands. Later on, magnetic islands become rounder and merge with each other to form bigger islands, until there is only a single big island and a single X-point left (per layer) at the end of reconnection (at tω₁ ≳ 450, Figure 2; the end of the simulation does not reach the full saturated state). This peculiar symmetric final state arises from the choice of double-periodic boundary conditions.

Zoom In Zoom Out Reset image size

We find that the reconnection rate increases by a factor of two when the radiation reaction is self-consistently included, β_rec ≈ 0.3. Because the initial pressure balance is approximately maintained during the simulation, the layer compresses due to high synchrotron energy losses inside the layer, leading to an increase in the reconnection rate (Uzdensky & McKinney 2011). The magnetic islands tend to be smaller and denser with radiative losses. Apart from this, we do not find any qualitative difference in the overall time evolution of the reconnection process, with or without the radiation reaction force. Figure 4 shows the time evolution of the different energy components in the system, i.e., the energy of the fields, of the particles and of the radiation. About 63% of the initial magnetic energy is dissipated and entirely radiated away by the particles at the end of the simulation. The total energy is conserved to within about 5% throughout the simulation. This moderate error comes from the overestimation of the synchrotron energy losses due to the high-frequency fluctuating electric field described in Section 3. If the radiation reaction force is neglected, or if there is a non-zero guide field, the total energy is very well conserved with less than 0.2% error.

Zoom In Zoom Out Reset image size

To summarize, the overall time evolution of reconnection can be divided schematically into three main phases shown in Figure 4:

• 1.  
  Phase 1. Quasi-static state where the electric field builds up smoothly with little magnetic reconnection. About 12% of the total magnetic energy is dissipated during this period.
• 2.  
  Phase 2. Plasmoid-dominated reconnection. This is the most active and bursty period of magnetic dissipation, during which about 40% of the total magnetic energy is dissipated. This phase is characterized by the strong competition between particle acceleration and cooling.
• 3.  
  Phase 3. Saturated state, where particle cooling dominates. The particles and the magnetic field are approximately in equipartition. The reconnection rate decreases so that only 11% of the total magnetic energy is dissipated during this period.

4.2. Particle and Synchrotron Radiation Spectra and Anisotropies

Figure 5 presents the energy distributions of all the background particles (γ²dN/dγ, top panel) and their instantaneous, transparent synchrotron radiation (νF_ν, bottom panel) at tω₁ = 0, 110, 220, 397, and 662. The contribution from the drifting particles is not shown here because of their small number compared with the background particles. We assume that the radiation is emitted continuously (valid if γB ≪ B_QED) and tangentially to the particle's orbit (valid if γ ≫ 1). Photons do not interact with the plasma (optically thin approximation), hence there is no need to solve the full radiative transfer equation. Figure 5 shows the emergence of a high-energy tail of particles, whose maximum energy increases with time during the most active phase of reconnection, tω₁ ≲ 450 (i.e., during the plasmoid-dominated phase). During the early stage of reconnection, the initial −2 power law extends to higher energies and cuts off exponentially. The spectrum marginally extends to lower energies. At tω₁ ≳ 220, there are some particles that are accelerated above the radiation-reaction limit energy γ_rad ≈ 1.3 × 10⁹. We find that these extremely energetic particles are accelerated at X-points, surfing multiple times across the reconnection layer, following relativistic Speiser orbits. The details of the Speiser acceleration mechanism are described below in Section 4.3. At tω₁ = 397, the spectrum above γ = 3 × 10⁸ is well described by a steep power law, i.e., dN/dγ∝γ^−3.8, followed by a sharp cut-off beyond γ = γ_rad (Figure 7, top panel). We find that about 0.03% of all the particles are above γ_rad, and represent about 0.39% of the total kinetic energy. These particles are responsible for the excess of synchrotron radiation above 160 MeV (see Figure 5), which represents about 4.3% of the total isotropic radiative power. After tω₁ = 450, the high end of the spectrum contracts to lower energies, and very few particles above γ_rad survive at tω₁ ≳ 662. We attribute the disappearance of the most energetic particles to synchrotron cooling. In the final saturated state, most particles are located within the big islands where they cool progressively. There is little acceleration and the magnetic field remains strong, of order B₀. If the radiation reaction force is neglected, the high-energy component of the spectrum does not evolve once established (for tω₁ ≳ 550; see Figure 5), and extends up to γ_max ≈ 3 × 10⁹. An identical simulation performed without radiative losses with Vorpal gives very similar results.

Zoom In Zoom Out Reset image size

We investigate the angular distribution of the particles' velocities, as a function of their energy. In agreement with our previous study (Cerutti et al. 2012b), we find a pronounced energy-dependent anisotropy of the particles and their synchrotron emission, increasing with energy. Figure 6 illustrates the strong anisotropy of the expected synchrotron radiation, as a function of the photon energy. Following Cerutti et al. (2012b), we use the spherical angles ϕ (latitude) and λ (longitude) to study the angular distributions. The latitude varies between −90° and +90°, and the longitude varies between −180° and +180°. A radial unit vector has the coordinates x = cos ϕsin λ, y = sin ϕ, z = cos ϕcos λ. At tω₁ = 397, we find that half of the >160 MeV radiative flux is concentrated into less than 4% of the total solid angle 4π. The high-energy beam of radiation is concentrated in the mid-plane (xz-plane, ϕ = 0°), preferentially toward the ±x-directions (ϕ = 0°, λ = ±90°). This result can be explained by the deflection of the particles' trajectories from the ±z-directions (ϕ = 0°, λ = 0°, ±180°, along which the particles are accelerated by Ez) to the ±x-directions by the reconnected field. The direction of the beam is changing with time, wiggling around the plane of the reconnection layer during the active phases of reconnection (tω₁ ≲ 450). The beam broadens and stabilizes along the z-direction at later times.

Zoom In Zoom Out Reset image size

The strong anisotropy of the emitted radiation leads to an apparent boosting of the flux seen by an observer looking in the direction of the beam (the so-called kinetic beaming; see Cerutti et al. 2012b). The energy distribution of the particles pointing in the direction λ = +70°, ϕ = 0° (indicated by the black box in Figure 6) within the solid angle ΔΩ/4π ≈ 3 × 10⁻³ is substantially harder than the isotropic one. A power-law fit yields dN/dγ∝γ^−2.5 above γ = 3 × 10⁸ up to γ_max ≈ γ_rad (Figure 7, top panel). In this case, the particles with Lorentz factors above γ_rad account for about 0.5% of the particles and 8.5% of the energy. Similarly, the apparent spectral energy distribution is dominated by the high-energy radiation, peaking at around 100 MeV with about 20% of the radiative power above 160 MeV. The spatial distribution of high-energy particles is strongly inhomogeneous (see also Cerutti et al. 2012b). Figure 8 shows the spatial distribution of a sample of high-energy particles with γ > 5 × 10⁸, at tω₁ = 309. The energetic particles are clustered into compact bunches within the reconnection layer and magnetic islands. We note that the high-energy particles are preferentially located at the periphery of the big islands or inside small, newly formed islands. The particles near the centers of big islands are not very energetic because they are no longer accelerated and have cooled radiatively over the time their host island has grown to a large size.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The pronounced inhomogeneity and anisotropy of the extremely energetic particles above γ_rad, and the associated radiation above 160 MeV, are key elements in explaining the Crab gamma-ray flares (see Section 5).

4.3. Relativistic Speiser Orbits

In this section, we examine in detail the acceleration mechanism of the most energetic particles with γ > γ_rad. We follow the trajectories of 20, 000 particles, picked randomly and uniformly throughout the box at t = 0. In this sample, there are about a dozen particles reaching a maximum energy above γ_rad. Figure 9 shows three typical particle trajectories projected in the yz-plane and in the xy-plane, as well as the time history along the particle trajectory of the particle Lorentz factor γ, the relative strength of the electric force and the radiation reaction force, and the synchrotron critical energy _sync = 3heB_⊥γ²/4πmc.

Zoom In Zoom Out Reset image size

These extremely energetic particles systematically follow the same simple time evolution, which can be decomposed into four main phases:

• 1.  
  Drifting. The particle is initially located upstream and is well magnetized where ideal MHD holds (E < B_⊥). The particle moves together with the field lines toward the layer, gains little energy, and does not radiate much. This phase ends at the reference point "in" in Figure 9 when the particle gets inside the layer.
• 2.  
  Linear acceleration. This phase starts when the particle reaches the reconnection layer, where it is no longer magnetized. The strong electric field Ez accelerates the particle almost linearly along the z-direction. In addition, the magnetic field B_x reversing across the layer confines the particle toward the layer mid-plane. The particle follows a relativistic Speiser orbit (Speiser 1965; Kirk 2004; Contopoulos 2007; Uzdensky et al. 2011), whose meandering width y_m (maximum distance of the particle from the neutral sheet) decreases with time as the particle energy increases. The particle is then confined closer and closer to the layer mid-plane, where the perpendicular magnetic field strength and hence the radiation reaction force decrease, while the electric force keeps on accelerating the particle. During this process, the particle's trajectory is also significantly deflected toward the ±x-directions by the reconnected field B_y. This phase begins at the reference point "in" and ends when the particle reaches its maximum energy, marked by the red cross in Figure 9.
• 3.  
  Ejection and emission. While it remains deep inside the layer, the particle is accelerated above the radiation reaction limit γ_rad. The particle is then ejected from the layer, in most cases when it encounters the final big magnetic island where it feels a sharp increase of the radiation reaction force. The big magnetic island acts effectively as a "beam dump." The particle loses most of its energy in a fraction of a Larmor cycle and emits synchrotron photons above 160 MeV. This phase happens just before and after the particle exits the layer (around the reference point "out" in Figure 9).
• 4.  
  Cooling. The particle is back in a region where ideal MHD conditions apply and cools progressively. The particle does not experience any significant acceleration once the saturated state of reconnection is reached.

This simple picture highlights the distinction between the acceleration zone (inside the layer) and the >160 MeV synchrotron radiating zone (upstream, or inside magnetic islands). The acceleration zone is of order the system size ∼L_x, while the radiating zone is a fraction of the Larmor radius ≪γ_radmc²/eB₀. The acceleration and focusing mechanisms described above in phase 2 agree surprisingly well with our previous test-particle simulations (Uzdensky et al. 2011; Cerutti et al. 2012a), despite the simplistic assumptions on the fields used in these studies (prescribed and static). Uzdensky et al. (2011) predicted a relationship between γ and the angle θ₀ between the particle's velocity vector and the layer mid-plane defined at each crossing, in two extreme regimes. If the particle's meandering width y_m is much greater than the layer thickness δ, and if radiative losses are negligible (i.e., during the first Speiser cycles), then |θ₀|∝γ^−2/3. In contrast, if the particle is deep inside the layer and reaches the local radiation reaction limit energy $\gamma ^{\prime }_{\rm rad}$ (defined with the perpendicular field at the location of the particle B_⊥ < B₀, so that $\gamma ^{\prime }_{\rm rad}>\gamma _{\rm rad}$) within each cycle, then |θ₀|∝γ^−3/2. Figure 10 shows the tracks followed by a representative sample of eight high-energy particles in the θ₀–γ plane (which are not necessarily accelerated above γ_rad). The mid-plane crossing angle is given by $\theta _0=\pi /2-\arccos (v_{\rm y}/\sqrt{\mathbf {v}\cdot \mathbf {v}})$, where v is the three-velocity vector of the particle. The agreement with the analytical expectations is very good: the particles remain between these two power laws, tending to a −2/3 index at low energies and to a −3/2 index at the highest energies. This is a robust and clean feature of the most energetic particles accelerated and focused through the Speiser mechanism.

Zoom In Zoom Out Reset image size

4.4. 100 MeV Emission">Variability Pattern of the >100 MeV Emission

In this section, we investigate the time-dependent radiation escaping in the +x-direction where most of the high-energy radiation is expected (Figure 6). Figure 11 presents the expected synchrotron flux integrated above 100 MeV as a function of time, taking into account the time delay due to the light crossing time through the box. In the case of radiation into the +x-direction, the propagation time is given by t_propag = (L_x − x_e)/c, where x_e is the location of the emitting electron/positron. In agreement with Cerutti et al. (2012b), the high-energy radiation is highly variable on timescales much shorter than the light crossing time of the layer (≲ 0.1 L_x/c, or ≲ 6 hr). The light curve is composed of multiple intense spikes that are nearly symmetric in time. This result is a direct consequence of the strong focusing of the energetic particles accelerated through the Speiser mechanism. The beam of energetic particles wiggles around the reconnection layer and crosses the line of sight several times. The bunching of the high-energy particles into compact blobs within the layer and within the magnetic islands also contributes to the multiple, powerful sub-flares in the light curve (Cerutti et al. 2012b). This dramatic variability disappears if, instead of considering one particular direction, the emission is averaged over all directions. Figure 11 also shows the energy dependence of the light curve. The amplitude of the spikes increases with the energy of the radiation considered, because of the increasing emission anisotropy (Figure 6). Figure 12 presents the resulting power density spectrum (PDS) of the light curve (given by the squared modulus of the fast Fourier transform), in the three energy bands defined in Figure 11. The observed PDS above 100 MeV is well fit by a hard power law of index ≈ − 0.5. At lower energies, the best-fit indexes are ≈ − 1.0 in the 10 MeV<₁ < 100 MeV band and ≈ − 1.2 in the 1 MeV<₁ < 10 MeV band. As expected, the PDS slope hardens with increasing photon energy, indicating that the highest energy radiation is also the most rapidly variable.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The received spectrum is also highly time variable. We decompose the light curve into 12 blocks of duration 12 hr each (see Figure 11). Within each period of time, we compute the time-averaged synchrotron spectrum received by the observer. We fit the high-energy component above 100 MeV only with a power law times an exponential cut-off, $\nu F_{\nu }=K_{\epsilon }\epsilon _1^{\alpha }\exp \left(-\epsilon _1/\epsilon _{\rm cut}\right)$, where the free parameters are the normalization constant K, the spectral index α, and the cut-off energy _cut. This analysis reveals a strong correlation between the cut-off energy and the total flux above 100 MeV (see Figure 13). A power-law fit gives νF_ν(₁ > 100 MeV)$\equiv \int _{100\, {\rm MeV}}^{+\infty }F_{\nu }d\epsilon _1\propto \epsilon _{\rm cut}^{+3.8\pm 0.6}$. To test the robustness of this correlation, we performed a series of 10 identical simulations, where only the initial positions and velocities of particles vary. We find a strong and positive flux/cut-off energy correlation in all the simulations, using 50 time bins to refine our analysis. The power-law index varies from about corr ≈ + 2 to corr ≈ + 4, with a mean index corr ≈ + 3 (see Table 2). We note that the slope depends on the low-energy cut-off (100 MeV). However, because the PDS varies significantly from one simulation to another, we are only able to see a hint of a reproducible pattern. On average, the PDSs are best fitted by power laws with the following indices: ≈ − 1.4 in the 1 MeV<₁ < 10 MeV band, ≈ − 1.3 in the 10 MeV<₁ < 100 MeV band, and ≈ − 1.0 above 100 MeV. The deviations from the mean indices within the sample are large, of order ±0.5.

Zoom In Zoom Out Reset image size

Table 2. Correlation between the Energy Flux and the High-energy Photon Spectral Cut-off Energy, for a Statistical Sample of 10 Identical Simulations Where Only the Initial Positions and Velocities of Particles Vary

Simulation	Best-fit Index Corr
1	3.14 ± 0.45
2	2.57 ± 0.44
3	3.26 ± 0.41
4	2.90 ± 0.50
5	1.96 ± 0.38
6	4.07 ± 0.43
7	3.33 ± 0.52
8	2.89 ± 0.39
9	3.05 ± 0.37
10	3.44 ± 0.56

Notes. The correlation coefficient is the best-fit power-law index as shown in Figure 13, i.e., νF_ν(₁ > 100 MeV) $\propto \epsilon _{\rm cut}^{\rm corr}$, using 50 bins in time.

Download table as:  ASCIITypeset image

4.5. Effect of the Guide Field

All the simulations presented above assumed a zero magnetic guide field component B_z = 0. In the general case, however, reconnection is not purely anti-parallel and there is always a finite guide field. In this section, we examine the effect of a moderate and uniform guide magnetic field on particle acceleration and radiation, by performing a set of four simulations with B_z/B₀ = 0.1, 0.25, 0.5, and 1, with the physical parameters given in Table 1. The general time evolution of the reconnection dynamics described in Section 4.1 is still valid with a guide field. However, we find that the morphology of the particle distribution inside the plasmoids is more complex. In the presence of a finite guide field, two streams of particles form symmetrically on both sides of the layer, in the opposite direction to the island motion (Figure 14, left panels). Each arm is produced by the deflection of the particles by the guide field in the ±y-directions, caused by the Lorentz force, and is composed of electrons or positrons only. As a result, the guide field creates a separation of charges within each island across the layer, which in return induces a strong electric field E_y of order B₀ across the layer (see Figure 14, right panels). We also note the induction of an E_x electric field of order 0.1B₀, reversing across the layer, due to the motion of the incoming plasma toward the layer at a velocity of order β_recc in the ±y-directions. We also find that the strong energy-dependent anisotropy of the particles and radiation is preserved with a guide field, although the structure of the beam is more complex (high-energy particles occasionally point toward high latitudes |ϕ| > 0°).

Zoom In Zoom Out Reset image size

The deflection of the particles away from the layer does not limit significantly the efficiency of particle acceleration through the Speiser mechanism. The particle energy distribution extends up to γ ≈ 2 × 10⁹. Analysis of a sample of particles tracked throughout the simulation shows that the energetic particles are inevitably deflected away from the layer, but only after several crossings of the layer. Figure 15 shows three typical high-energy particle trajectories, projected in the yz-plane. The particles are progressively carried outside the layer by the guide field, and have time to cool over a Larmor timescale. The radiation reaction force gradually overcomes the electric force because the pitch angle of the particle to the magnetic field line (and hence the perpendicular magnetic field) increases progressively as it exits the layer. In contrast, the ejection of the energetic particles into strong B_⊥ in the weak guide field case is abrupt: the radiation reaction force jumps rapidly over a short period of time and exceeds the electric force, leading to a rapid (sub-Larmor) cooling and the emission of >160 MeV synchrotron photons (Figure 9). Hence, the presence of a strong guide field (B_z ≳ B₀) effectively inhibits the emission of >160 MeV synchrotron radiation. Figure 16 shows the fraction of the total radiative flux emitted above 160 MeV, as a function of the guide field strength, from B_z = 0 to B_z = B₀. For B_z = B₀, the >160 MeV flux is about 50 times smaller than the zero-guide field case. We note that the >160 MeV flux varies by a factor 2–3 within our statistical sample of 10 simulations with no guide field introduced in Section 4.4, as illustrated in Figure 16.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

A moderate guide field may be beneficial for particle acceleration in 3D pair plasma reconnection. It tends to suppress the development of the drift kink instability that occurs only in 3D, which may quench non-thermal particle acceleration (Zenitani & Hoshino 2008; see also the discussion in Liu et al. 2011; Sironi & Spitkovsky 2011; Kagan et al. 2012). A large 3D simulation such as those presented here in 2D is currently beyond our computational resources. We leave this issue to a future study.

In this section, we discuss the implications of our findings in the context of the Crab gamma-ray flares. Below, we review the observational features of the flares, and show how reconnection can naturally explain them:

• 1.  
  Synchrotron radiation >160 MeV. All flares show synchrotron radiation well above the 160 MeV burnoff limit, up to 375 MeV during the 2011 April super-flare (Buehler et al. 2012). We showed in Section 4.3 that particles can be accelerated well above the classical radiation reaction limit, via the relativistic Speiser mechanism deep inside the layer, and radiate synchrotron radiation up to ≳ 400 MeV. The most energetic particles go through a substantial fraction of the electric potential drop available in the simulation, i.e., $\mathcal {E}_{\rm max}\sim e\beta _{\rm rec}B_0 L_{\rm x}\approx 3$ PeV, with β_rec = E_z/B₀ ≈ 0.3, B₀ = 5 mG, and L_x = 2.7 lt-day.
• 2.  
  Hard spectrum. The 2011 April flare unambiguously shows an extra synchrotron component on top of the quiescent emission. The spectrum averaged over the flare is hard, such that F_ν∝ν^{−0.27 ± 0.12} (Buehler et al. 2012), i.e., consistent with an emitting population of pairs distributed with a power law of index p ≈ 1.6. The energy of the distribution is then dominated by the highest energy particles. We find that a hard spectrum (although not as hard as observed) is naturally expected in our simulations during the brightest periods of high-energy emission (periods 8, 9, and 10 in Figure 11), when F_ν∝ν^−0.4, ν^−0.5 above 100 MeV. A hard spectrum could explain why there is no detectable counterpart at lower energies (X-rays, optical, IR, radio).
• 3.  
  Ultra-rapid time variability. The first short intra-flare variability was detected in the 4-day-long 2010 September flare, over a 12 hr timescale (Balbo et al. 2011). During the 9-day-long 2011 April flare, there are significant variations of the flux over periods ≲ 8 hr (Buehler et al. 2012). Hence, there are even smaller structures in the flaring region emitting each spike of high-energy radiation. Along the lines of Cerutti et al. (2012b), we find in Section 4.4 that super-fast time variability of the observed high-energy flux is expected, due to the strong inhomogeneity and anisotropy of the most energetic particles responsible for the >100 MeV emission. We see a flare on Earth only when the beam of high-energy radiation crosses our line of sight. The symmetry of the sub-flare time profile is also consistent with observations. We attribute the overall duration of the Crab flares (from a few days to a couple of weeks) to the reconnection timescale, i.e., of order the light crossing time of the reconnecting region, 6 days in our simulations (Figure 11). We associate the intra-flare time variability (hours to days) with the light crossing time of magnetic islands generated by the tearing instability, ≲ 6 hr in the simulations (see also Giannios 2013 for a similar interpretation in the context of super-fast TeV flares in blazars). Our study also reveals that the PDS of the >100 MeV light curve is expected to be a hard power law of index ≈ − 1 ± 0.5 (Figure 12). This is consistent with observations that show a power-law index ≈ − 1 above the noise floor (Buehler et al. 2012).
• 4.  
  Flux/energy correlation. One of the most remarkable features of the 2011 April flare is the discovery of a clear correlation between the >100 MeV flux and the cut-off energy, such that νF_ν(> 100 MeV)$\,{\propto}\, \epsilon _{\rm cut}^{+3.42\pm 0.86}$ (Buehler et al. 2012). This result is interpreted by Buehler et al. (2012) as the signature of the rapid variations of relativistic Doppler-boosted emission (Lind & Blandford 1985). We discovered in Section 4.4 that reconnection can reproduce such a correlation as well, with a power-law index corr ≈ + 3, and mimic the effect of Doppler beaming. Relativistic bulk outflows with Lorentz factor of order $\Gamma \sim \sqrt{\sigma }=4$ are expected in relativistic Petschek reconnection (Lyubarsky 2005), and lead to a Doppler amplification of the observed emission (Giannios et al. 2009, 2010; Nalewajko et al. 2011; Clausen-Brown & Lyutikov 2012; Giannios 2013). However, a Petschek-like reconnection configuration is traditionally associated with either the Hall effect or a highly localized anomalous resistivity due to current-driven plasma microinstabilities, both of which are absent in our 2D pair plasma simulations.⁸ Hence, the scaling $\Gamma \sim \sqrt{\sigma }$ may not be valid in the plasmoid-dominated reconnection regime explored by the simulations. Our results suggest that the reconnection outflows are not strongly relativistic, i.e., Γ ∼ 1, as expected in relativistic Sweet–Parker reconnection (Lyubarsky 2005). Here, we attribute the beaming of the radiation to the strong energy-dependent anisotropy of the particles only (Section 4.2; see also Cerutti et al. 2012b). The beaming pattern obtained in the simulations is not consistent with a relativistic Doppler effect that beams all the emission by the same factor, regardless of the particle energy. This is a possible observational test of the kinetic beaming.

The main result of this paper is the first discovery of particle acceleration above the classical synchrotron radiation reaction limit, and the associated emission of >160 MeV synchrotron radiation, in numerical simulations of collisionless pair plasma reconnection. For this purpose, we developed a new, parallel, relativistic PIC code, called Zeltron, that includes the effects of the radiation reaction force on particle dynamics self-consistently. We confirm in every detail the expectations of earlier (semi-)analytical works (Kirk 2004; Contopoulos 2007; Uzdensky et al. 2011; Cerutti et al. 2012a): the most energetic particles are almost linearly accelerated and confined deep inside the reconnection layer, where radiative losses are small while the accelerating electric field is strong. If the particle stays long enough within the layer, it can be accelerated above the radiation reaction limit defined in the upstream plasma. Then, once the particle eventually escapes the layer, it radiates >160 MeV synchrotron radiation within a fraction of a Larmor gyration. The acceleration mechanism that emerges from this study is remarkably simple and robust.

Although this work addresses a fundamental question in particle acceleration, it is mainly motivated by the mystery of the Crab Nebula gamma-ray flares. We have shown in this study that all the puzzling aspects of the flares are consistent with a reconnection event in the nebula. In addition to the >160 MeV synchrotron radiation, our simulations can explain the ultra-rapid time variability of the observed high-energy radiation. The emission originates mostly from extremely anisotropic and compact bunches of energetic particles. An external observer sees a flare when the beam of energetic particles crosses the line of sight. We found that the power spectrum of the light curve is well fit by a power law, which shows hints of hardening with the energy band of the radiation. In addition, we discovered that there is a strong positive correlation between the emitted radiative flux and the high-energy spectral cut-off, mimicking the effect of a relativistic Doppler boost. This correlation is also consistent with the observations of the Crab flares. A strong guide field, i.e., ≳ B₀, deflects particles out of the reconnection layer and tends to suppress the emission >160 MeV synchrotron radiation. Our results support the magnetic reconnection scenario at the origin of the Crab flares.

We thank Rolf Buehler, Dimitrios Giannios, Geoffroy Lesur, Krzysztof Nalewajko, Anatoly Spitkovsky, and the referee for helpful discussions. This research was supported by an allocation of advanced computing resources provided by the National Science Foundation, by NSF grant PHY-0903851, NSF grant AST-0907872, DoE grant DE-SC0008409, and NASA Astrophysics Theory Program grant NNX09AG02G. Numerical simulations were performed on the local CIPS computer cluster Verus, on Kraken at the National Institute for Computational Sciences (www.nics.tennessee.edu/). This work also utilized the Janus supercomputer, which is supported by the National Science Foundation (award number CNS-0821794), the University of Colorado Boulder, the University of Colorado Denver, and the National Center for Atmospheric Research. The Janus supercomputer is operated by the University of Colorado Boulder.