THE EXTENT OF POWER-LAW ENERGY SPECTRA IN COLLISIONLESS RELATIVISTIC MAGNETIC RECONNECTION IN PAIR PLASMAS

Magnetic reconnection is a fundamental plasma physics process in which magnetic field rearrangement and relaxation rapidly converts magnetic energy into particle energy (Zweibel & Yamada 2009). Reconnection is believed to drive many explosive phenomena in the universe, from Earth magnetospheric substorms and solar flares to high-energy X-ray and γ-ray flares in various astrophysical objects. Quite often, the radiation spectra of these flares, and hence the energy distributions of the emitting particles, are observed to be non-thermal (e.g., characterized by power laws). Therefore, understanding the mechanisms of nonthermal particle acceleration and determining the observable characteristics—such as the power-law index and high-energy cutoff—of the resulting particle distribution, is an outstanding problem in modern heliospheric physics and plasma astrophysics.

Of particular interest in high-energy astrophysics is the role of relativistic reconnection—which occurs when the energy density of the reconnecting magnetic field, ${B}_{0}^{2}/8\pi$, exceeds the rest-mass energy density ${n}_{b}{{mc}}^{2}$ of the ambient plasma, leading to relativistic bulk outflows and plasma heating to relativistic temperatures—as a potentially important mechanism for nonthermal particle acceleration to ultra-relativistic energies (with Lorentz factors $\gamma \gg 1$) in various astrophysical sources (Hoshino & Lyubarsky 2012). In particular, this process has been invoked to explain energy dissipation and radiation production in electron–positron (pair) plasmas over multiple scales in pulsar systems—e.g., in the pulsar magnetosphere near the light cylinder, in the striped pulsar wind, and in the pulsar wind nebula (PWN; Coroniti 1990; Lyubarsky 1996; Lyubarsky & Kirk 2001; Sironi & Spitkovsky 2011; Uzdensky et al. 2011; Cerutti et al. 2012a, 2013, 2014a, 2014b; Uzdensky & Spitkovsky 2014). In addition, relativistic reconnection in pair and/or electron–ion plasmas is believed to play an important role in gamma-ray bursts (GRBs; Drenkhahn & Spruit 2002; Giannios & Spruit 2007; McKinney & Uzdensky 2012) and in coronae and jets of accreting black holes, including active galactic nucleus (AGN)/blazar jets, e.g., in the context of TeV blazar flares (Giannios et al. 2009; Nalewajko et al. 2011).

Nonthermal particle acceleration is essentially a kinetic (i.e., non-fluid) phenomenon. Although fluid simulations with test particles have been used to study particle acceleration, particle-in-cell (PIC) simulations include kinetic effects self-consistently. A number of PIC studies have investigated particle acceleration in collisionless relativistic pair-plasma reconnection (Zenitani & Hoshino 2001, 2005, 2007, 2008; Jaroschek et al. 2004; Lyubarsky & Liverts 2008; Liu et al. 2011; Sironi & Spitkovsky 2011; Bessho & Bhattacharjee 2012; Cerutti et al. 2012a, 2013, 2014a, 2014b; Kagan et al. 2013, 2015; Liu et al. 2015); the best evidence for nonthermal particle distributions was provided recently by Sironi & Spitkovsky (2014), Guo et al. (2014). Whereas previous studies have identified power-law slopes of nonthermal spectra, the important question of the energy extent of these power laws has not been systematically addressed.

In this letter we present a comprehensive two-dimensional (2D) PIC investigation of non-thermal particle acceleration in collisionless relativistic reconnection in a pair plasma without guide magnetic field. In particular, we characterize the dependence of the resulting energy distribution function on the system size L and the upstream "cold" magnetization parameter $\sigma \equiv {B}_{0}^{2}/(4\pi {n}_{b}{m}_{e}{c}^{2})$ (relativistic reconnection requires $\sigma \gg 1$). We find empirically that relativistic reconnection produces a high-energy spectrum that is well represented by a power law with exponential and super-exponential cutoffs:

$\begin{eqnarray}&&f(\gamma )=\displaystyle \frac{{dN}}{d\gamma }\propto {\gamma }^{-\alpha }\mathrm{exp}(-\gamma /{\gamma }_{c1}-{\gamma }^{2}/{\gamma }_{c2}^{2}).\end{eqnarray} \tag{ 1 }$

The different cutoffs serendipitously allow us to distinguish different scalings with σ and L: γ_c1 ∼ 4σ depends on σ, while ${\gamma }_{c2}\sim 0.1\;L/{\rho }_{0}$ depends on L (here ${\rho }_{0}\equiv {m}_{e}{c}^{2}/{{eB}}_{0}$ is the nominal Larmor radius).

Equality of the two cutoffs, ${\gamma }_{c1}\simeq {\gamma }_{c2}$, defines a critical size ${L}_{c}\simeq 40\sigma {\rho }_{0}$ separating the small- and large-system regimes. We find that for large systems ($L/\sigma {\rho }_{0}\gg 40$), the energy spectrum of accelerated particles (hence ${\gamma }_{c1}$) is essentially independent of L. Importantly (as we discuss later), L_c is approximately the length at which a current layer, with thickness equal to the average Larmor radius ${\rho }_{e}=\bar{\gamma }{\rho }_{0}$, becomes tearing-unstable and breaks up into multiple plasmoids and secondary current sheets. (Here, $\bar{\gamma }{m}_{e}{c}^{2}$ is the average dissipated energy per background particle, $\bar{\gamma }\simeq \kappa ({B}_{0}^{2}/8\pi )/({n}_{b}{m}_{e}{c}^{2})=\kappa \sigma /2;$ in our simulations $\kappa \simeq 0.6$, so $\bar{\gamma }\simeq 0.3\sigma$.) Therefore, we propose that (at least in 2D with an initially cold background plasma) reconnection in the plasmoid-dominated regime yields a high-energy particle spectrum that is predominantly independent of system size $L\gg {L}_{c}$. Consequently, nonthermal particle acceleration in huge, astrophysically relevant systems may be studied via merely large simulations, i.e., with $L\gtrsim {L}_{c}$.

This study focuses on reconnection in 2D without guide field (${B}_{z}=0$). Although some important differences in the reconnection dynamics emerge between 2D and 3D, such as the development of the drift-kink instability (Zenitani & Hoshino 2008), the dimensionality is not believed to affect the particle energy spectra at late stages (W. Daughton 2014, private communication; J. F. Drake 2014, private communication; Guo et al. 2014; Sironi & Spitkovsky 2014). Working in 2D (much less costly than 3D) enabled investigation of large system sizes.

We simulate systems of size ${L}_{x}={L}_{y}=L$ with periodic boundary conditions and two antiparallel reconnection layers. The two layers begin as relativistic Harris current sheets (Kirk & Skjæraasen 2003) with upstream magnetic field ${B}_{x}={B}_{0}$ and a peak drifting plasma simulation-frame-density n_d (at the layer centers) that is 10 times the uniform background density n_b. A small (1%) initial magnetic-flux perturbation facilitates reconnection onset. Electrons and positrons in each Harris layer drift (in opposite directions) with average velocity ${\beta }_{d}c=0.6c$, and initial Maxwell–Jüttner temperature ${\theta }_{d}\equiv {k}_{B}{T}_{d}/{m}_{e}{c}^{2}=\sigma /16;$ the initial layer half-thickness is $\delta =(8/3){\theta }_{d}\;{\rho }_{0}=\sigma {\rho }_{0}/6$. The background plasma is initially at rest, with temperature ${\theta }_{b}\ll \sigma ;$ however, due to the finite grid instability (Birdsall & Maron 1980), the background plasma is expected to heat until its Debye length is resolved, which occurs at a temperature of order ${\theta }_{D}\sim \sigma /512$ for ${\rm{\Delta }}x=\sigma {\rho }_{0}/32$.

The simulations begin with ${N}_{p}=128$ (macro)particles per grid cell with cell sizes ${\rm{\Delta }}x={\rm{\Delta }}y={\theta }_{d}\;{\rho }_{0}/2=\sigma {\rho }_{0}/32\approx 0.2\delta$ (except for $\sigma =3$, where the mildly relativistic particles allowed ${\rm{\Delta }}x={\theta }_{d}{\rho }_{0}$ without sacrificing accuracy). The total energy is conserved within 1% during each simulation. Convergence tests with respect to ${\rm{\Delta }}x$ and N_p indicate that our simulations are well resolved and, in particular, that the high-energy parts of the particle distributions are robust.

The Vorpal code (Nieter & Cary 2004), employed for this study, uses an explicit electromagnetic PIC time advance, with Yee electromagnetics and a relativistic Boris particle push.

To determine the power-law index α and the energy cutoffs ${\gamma }_{c1}$ and ${\gamma }_{c2}$ as functions of the upstream magnetization σ and the system size L, we ran simulations with σ = 3, 10, 30, 100, 300, 1000, and, for each σ, a range of system sizes up to $L/\sigma {\rho }_{0}=100$ for $\sigma =1000,300$, up to $L/\sigma {\rho }_{0}=200$ for $\sigma =100,30,10$, and up to $L/\sigma {\rho }_{0}=400$ for σ = 3.

We focus on the energy distribution of background particles, excluding the initially drifting particles, which contribute negligibly to the overall distribution for large L. Evolution to a nonthermal distribution proceeds rapidly (Figure 1); and the shape of the high-energy spectrum, as characterized by α and ${\gamma }_{c1/2}$, ceases to evolve well before all available flux has reconnected, especially for large systems.

Zoom In Zoom Out Reset image size

We find that the late-time high-energy spectrum is a power law with a high-energy cutoff significantly above the average particle energy, in agreement with Sironi & Spitkovsky (2014). We further observe (Figure 1) that spectra for large systems have exponential cutoffs, $\mathrm{exp}(-\gamma /{\gamma }_{c1})$, while small systems have sharper cutoffs, which we empirically model with a super-exponential $\mathrm{exp}(-{\gamma }^{2}/{\gamma }_{c2}^{2})$. We therefore fit all spectra with the universal form of Equation (1) to determine the power-law index α and the cutoffs ${\gamma }_{c1}$, ${\gamma }_{c2};$ for small systems, the best-fit ${\gamma }_{c1}$ is typically much larger than ${\gamma }_{c2}$ (hence irrelevant and highly uncertain), while for large L, ${\gamma }_{c2}$ is larger and uncertain.

Each spectrum is fit to Equation (1) over an interval $[{\gamma }_{f1},{\gamma }_{f2}]$, chosen as large as possible while maintaining a good fit. Because spectra depart from a power law at lowest energies, and because of increased noise at highest energies, larger fitting intervals yield unacceptably poor fits. Noise is reduced (and fit improved) by averaging over short time intervals and, if available, over multiple simulations (identical except for randomized initial particle velocities). Because the choices of acceptable fit quality and the durations of averaging intervals are somewhat subjective, we perform many fits using different choices, and finally report the median values with "error" bars encompassing the middle 68% of the fits (i.e., ±1 standard deviation if the data were Gaussian-distributed); small error bars thus demonstrate insensitivity to the fitting process. Very uncertain and large (hence irrelevant) cutoff values are discarded.

By applying this fitting procedure to the background particle spectrum for each different value of $(L,\sigma )$, we mapped out α, ${\gamma }_{c1}$, and ${\gamma }_{c2}$ as functions of σ and L, up to sufficiently large L to estimate the asymptotic values ${\alpha }_{*}(\sigma )={\mathrm{lim}}_{L\to \infty }\alpha (\sigma ,L)$ (Figure 2). We find that ${\alpha }_{*}(\sigma )$ starts above 2 for modest σ, and decreases to ${\alpha }_{*}(\sigma )\approx 1.2$ in the ultra-relativistic limit of $\sigma \gg 1$ (Figure 3), a result that is broadly consistent with previous studies (Zenitani & Hoshino 2001; Jaroschek et al. 2004; Lyubarsky & Liverts 2008; Guo et al. 2014; Melzani et al. 2014; Sironi & Spitkovsky 2014); while our measurement is closer to 1.2 than 1, the uncertainty is too large to rule out ${\alpha }_{*}\to 1$, predicted by some (Larrabee et al. 2003; Guo et al. 2014).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In contrast to the power-law index α, the energy extent of the power law has received relatively little attention in relativistic reconnection literature (Larrabee et al. 2003; Lyubarsky & Liverts 2008). We find that the high-energy cutoffs scale as ${\gamma }_{c1}\sim 4\sigma$ (independent of L) and ${\gamma }_{c2}\sim 0.1\;L/{\rho }_{0}$ (independent of σ) (Figures 4 and 5). Thus $L/\sigma {\rho }_{0}\ll 40$ implies ${\gamma }_{c2}\ll {\gamma }_{c1}$, and a super-exponential cuts off the power-law at an energy determined by the system size. Larger system sizes $L/\sigma {\rho }_{0}\gg 40$ have ${\gamma }_{c1}\ll {\gamma }_{c2}$, and so ${\gamma }_{c1}$ determines where the power law ends, independent of L.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The scaling of the high-energy cutoffs can be explained in terms of the distance a particle could travel within the reconnection field ${E}_{z}\sim {\beta }_{r}{B}_{0}$ (where ${\beta }_{r}\sim 0.1$ is the reconnection rate). By calculating analytic trajectories in fields around a single X-point, (Larrabee et al. 2003) concluded that $f(\gamma )\propto {\gamma }^{-1}\mathrm{exp}(-\gamma /{{\rm{\Gamma }}}_{0})$ with ${{\rm{\Gamma }}}_{0}=12e{\beta }_{r}^{2}{B}_{0}{{\ell }}_{x}/$ ${m}_{e}{c}^{2}\sim {{eE}}_{z}{{\ell }}_{x}/$ ${m}_{e}{c}^{2}\sim 0.1{{\ell }}_{x}/{\rho }_{0}$, with ${{\ell }}_{x}$ being the size of the reconnection region in x,⁸ a result that was supported by 2D PIC simulation in Lyubarsky & Liverts (2008).

In general, small systems reconnect mainly with one X-point, so ${{\ell }}_{x}\sim L$ and ${{\rm{\Gamma }}}_{0}\sim 0.1\;L/{\rho }_{0}$, which equals our ${\gamma }_{c2}$. (The observed super-exponential form presumably results from the simulation's boundary conditions.)

In large systems, however, the tearing instability breaks up current layers with full-length greater than ${{\ell }}_{{\rm{tear}}}\sim 100\bar{\delta }$, where $\bar{\delta }$ is the layer half-thickness (Loureiro et al. 2005; Ji & Daughton 2011), resulting in a hierarchy of layers ending with elementary layers, which are marginally stable against tearing (Shibata & Tanuma 2001; Loureiro et al. 2007; Uzdensky et al. 2010). The half-thickness of elementary (single X-point, laminar) layers should be about the average Larmor radius $\bar{\delta }\sim {\rho }_{e}=\bar{\gamma }{\rho }_{0}$ (Kirk & Skjæraasen 2003). Although Larrabee et al. (2003) considered single X-point reconnection, we propose that their formula for ${{\rm{\Gamma }}}_{0}$ can also be used in the context of plasmoid-dominated reconnection in large systems if applied to elementary layers (instead of the entire global layer): ${{\ell }}_{x}\sim {{\ell }}_{{\rm{tear}}}\sim 100\bar{\gamma }{\rho }_{0}\sim 30\sigma {\rho }_{0}$ (instead of ${{\ell }}_{x}\sim L$). Then, ${{\rm{\Gamma }}}_{0}\sim 0.1{{\ell }}_{{\rm{tear}}}/{\rho }_{0}\sim 3\sigma$, which is essentially our ${\gamma }_{c1}$ (and consistent with the measurement of ${{\rm{\Gamma }}}_{0}=35$ for $\sigma =9$ in Lyubarsky & Liverts 2008).

This explanation of high-energy-cutoff scaling in terms of elementary layer lengths may be robust despite the potentially important roles played by other acceleration mechanisms (Hoshino & Lyubarsky 2012). For example, significant additional acceleration may occur within contracting plasmoids (Drake et al. 2006; Dahlin et al. 2014; Guo et al. 2014, 2015) or—especially for the highest-energy particles—in the (anti-)reconnection electric field of secondary plasmoid mergers (Oka et al. 2010; Sironi & Spitkovsky 2014; K. Nalewajko et al. 2015, in preparation).

It is interesting to compare our high-energy cutoffs to the upper bound imposed on a power-law distribution by a finite energy budget. When $1\lt \alpha \lt 2$, most of the kinetic energy resides in high-energy particles, so the available energy per particle $\bar{\gamma }\sim 0.3\sigma$ limits the extent of the power law. If $f(\gamma )\sim {\gamma }^{-\alpha }$ extends from ${\gamma }_{{\rm{min}}}$ to ${\gamma }_{{\rm{max}}}\gg {\gamma }_{{\rm{min}}}$, then $\bar{\gamma }\approx [(\alpha -1)/(2-\alpha )]{\gamma }_{{\rm{min}}}^{\alpha -1}{\gamma }_{{\rm{max}}}^{2-\alpha }$ (Sironi & Spitkovsky 2014). For $\alpha \approx 1$, ${\gamma }_{{\rm{max}}}$ can extend well beyond $\bar{\gamma }$, but ${\gamma }_{{\rm{max}}}/\bar{\gamma }$ depends weakly on system parameters, consistent with our finding ${\gamma }_{c1}\sim \bar{\gamma }\sim \sigma$. E.g., for $\alpha =1.2$, ${\gamma }_{{\rm{max}}}/\bar{\gamma }\approx {({10}^{3}\bar{\gamma }/{\gamma }_{{\rm{min}}})}^{1/4}$. However, when $\alpha \gt 2$ (e.g., for low σ), the energy budget imposes no upper bound, since ${\int }_{{\gamma }_{{\rm{min}}}}^{\infty }\gamma {\gamma }^{-\alpha }d\gamma$ is finite. Nevertheless, for $\sigma =3$ where ${\alpha }_{*}\gt 2$, we observe ${\gamma }_{c1}\sim 4\sigma$, the same as for smaller ${\alpha }_{*}$.

The exponential cutoff at energies above ${\gamma }_{c1}\sim 4\sigma \sim 10\bar{\gamma }$ has important astrophysical implications for particle acceleration in systems such as pulsar magnetospheres, winds, PWN, and relativistic jets in GRBs and AGNs. Our results (insofar as they are ultra-relativistic) can be generalized to relativistically hot upstream plasmas by scaling all the energies by ${\bar{\gamma }}_{b}$, the average Lorentz factor of background particles. The "hot" magnetization ${\sigma }^{({\rm{hot}})}\equiv {B}_{0}^{2}/(4\pi {nw})$ therefore parameterizes similar simulations, since the relativistic specific enthalpy w also scales with ${\bar{\gamma }}_{b}$ (i.e., $w={\bar{\gamma }}_{b}{m}_{e}{c}^{2}+{p}_{b}/{n}_{b}$, where p_b is the background plasma pressure; for ${\bar{\gamma }}_{b}\gg 1$, $w\approx (4/3){\bar{\gamma }}_{b}{m}_{e}{c}^{2}$).⁹ For example, our reconnection-based model (Uzdensky et al. 2011; Cerutti et al. 2012a, 2013, 2014a, 2014b) for high-energy γ-ray flares in the Crab PWN (Abdo et al. 2011; Tavani et al. 2011) relies upon acceleration of a significant number of particles from ${\bar{\gamma }}_{b}\sim 3\times {10}^{6}$ to $\gamma \gtrsim {10}^{9}$. If, to achieve this, we need ${\gamma }_{c1}\gt {10}^{9}$, then direct extrapolation of the results from this letter would require ${\sigma }^{({\rm{hot}})}\gtrsim (1/4){\gamma }_{c1}/(w/{m}_{e}{c}^{2})\approx 60;$ this should be comparable (via scaling equivalence) to simulations presented in this work with $\sigma \sim 60$ (corresponding to a power-law index ${\alpha }_{*}\sim 1.3$). This required ${\sigma }^{({\rm{hot}})}$ is significantly higher than what is expected in the Crab Nebula. However, here we analyzed the entire spectrum of background particles, while (Cerutti et al. 2012b) suggested that bright flares observed in the Crab Nebula result from preferential focusing of the highest-energy particles into tight beams with energy spectra that differ from the entire spectrum. We also note that our present simulations are initialized with a Maxwellian plasma, whereas the ambient plasma filling the Crab Nebula has a power-law distribution, which may result in a higher high-energy cutoff.

We ran a series of collisionless relativistic pair-plasma magnetic reconnection simulations with no guide field, covering a wide range of system sizes L and upstream magnetizations $\sigma \geqslant 3$. We observed acceleration of the background plasma particles to a nonthermal energy distribution $f(\gamma )\sim {\gamma }^{-\alpha (L,\sigma )}$ with a high-energy cutoff. The cutoff energy is proportional to the maximum length of elementary, single X-point layers, which is limited by L in small systems, and by the secondary tearing instability in large systems. For small systems ($L\ll 40\sigma {\rho }_{0}$) we observe $f(\gamma )\sim {\gamma }^{-\alpha }\mathrm{exp}(-{\gamma }^{2}/{\gamma }_{c2}^{2})$ with ${\gamma }_{c2}\sim 0.1\;L/{\rho }_{0}$, and for large systems, $f(\gamma )\sim {\gamma }^{-\alpha }\mathrm{exp}(-\gamma /{\gamma }_{c1})$ with ${\gamma }_{c1}\sim 4\sigma$. As L becomes large, the power-law index $\alpha (L,\sigma )$ asymptotically approaches ${\alpha }_{*}(\sigma )$, which in turn decreases to $\approx 1.2$ as $\sigma \to \infty$. This characterization of power-law slope and high-energy cutoffs can be used to link ambient plasma conditions (i.e., σ) with observed radiation from high-energy particles, to investigate the role that reconnection plays in high-energy particle acceleration in the universe.

This work was supported by DOE grants DE-SC0008409 and DE-SC0008655, NASA grant NNX12AP17G, and NSF grant AST-1411879. Numerical simulations were made possible by the Extreme Science and Engineering Discovery Environment (XSEDE), which is supported by National Science Foundation (NSF) grant number ACI-1053575—and in particular by the NSF under Grant numbers 0171134, 0933959, 1041709, and 1041710 and the University of Tennessee through the use of the Kraken computing resource at the National Institute for Computational Sciences (www.nics.tennessee.edu/). This work also used the Janus supercomputer, which is supported by the NSF (award number CNS-0821794) and the University of Colorado Boulder; the Janus supercomputer is a joint effort of the University of Colorado Boulder, the University of Colorado Denver, and the National Center for Atmospheric Research. We gratefully acknowledge the developers responsible for the Vorpal code.