EFFICIENT PRODUCTION OF HIGH-ENERGY NONTHERMAL PARTICLES DURING MAGNETIC RECONNECTION IN A MAGNETICALLY DOMINATED ION–ELECTRON PLASMA

It is thought that the magnetic field is of central importance in pulsar wind nebulae (PWNe) and jets from active galactic nuclei (AGNs) and gamma-ray bursts (GRBs). The initially launched relativistic flows are Poynting-flux dominated (Spruit 2010). Meanwhile, these astrophysical phenomena are observed to have nonthermal spectra that arise from various emission processes of energetic particles. For example, leptonic modeling of TeV blazar emissions requires a power-law electron distribution with energy ${\gamma }_{e}$ ∼ 10³–10⁵ (Böttcher et al. 2013). The spectra of GRBs during the prompt phase suggest that the radiating particles have Lorentz factors ${\gamma }_{e}$ ∼ 10³–10⁷ (Daigne & Mochkovitch 1998; Uhm & Zhang 2014). In hadronic models, protons also need to be accelerated to high energies and may be related to ultra-high-energy cosmic rays (Böttcher et al. 2013). Understanding how the energy stored in magnetic fields being converted into nonthermal particle energy is essential to interpret observed emissions in those high-energy astrophysical processes.

Recent studies of GRBs (Zhang & Yan 2011; McKinney & Uzdensky 2012; Kumar & Crumley 2015), AGNs (Giannios et al. 2009; Zhang et al. 2015a), and PWNe (Uzdensky & Spitkovsky 2014) have shown that models that are magnetically dominant can explain much of the observed data. It is estimated that the magnetization parameter $\sigma \equiv {B}^{2}/(4\pi w)$ can be large and the Alfvén speed approaches the light speed ${v}_{A}\sim c$. Here B is the magnitude of magnetic field and w is the enthalpy. To explain the observed high-energy emissions, an efficient conversion from energies in the magnetized flow into nonthermal particles is required (e.g., Zhang & Yan 2011). Collisionless shocks that are considered to be efficient for particle acceleration in weakly magnetized plasmas are inefficient in dissipating energy and accelerating nonthermal particles in magnetically dominated flows (Sironi et al. 2015; Zhang et al. 2015b).

Magnetic reconnection is a major candidate that unleashes magnetic energy in magnetized flows by reorganizing magnetic topology (Zweibel & Yamada 2009). In the relativistic limit, the magnetic energy gets strongly dissipated in the reconnection region (Blackman & Field 1994; Lyutikov & Uzdensky 2003; Comisso & Asenjo 2014; Liu et al. 2015).

The theory of nonthermal particle acceleration in reconnection layers is still evolving (de Gouveia dal Pino & Lazarian 2005; Drake et al. 2006; Drury 2012; Guo et al. 2014; Zank et al. 2014). Much progress has been made through particle-in-cell (PIC) kinetic simulations that self-consistently model the microscopic physics of collisionless reconnection and nonthermal plasma energization. Earlier studies focused on the low-σ regime ($\sigma \lt 1$) and identified several basic acceleration mechanisms, such as direct acceleration at X-lines (Fu et al. 2006; Huang et al. 2010), Fermi-type acceleration associated with plasmoids (Drake et al. 2006), and further acceleration in plasmoid coalescence regions (Oka et al. 2010). Simulations with mildly relativistic parameters (σ ∼ 1) have uncovered particle acceleration with power-law energy spectra in X-lines (Zenitani & Hoshino 2001; Bessho & Bhattacharjee 2012), but no significant power-law distribution when the spectra are integrated over the whole domain (Zenitani & Hoshino 2007; Liu et al. 2011; Kagan et al. 2013). Recently, there has been an explosion of studies on relativistic magnetic reconnection in the magnetically dominated regime ($\sigma \gg 1$) (Cerutti et al. 2012; Guo et al. 2014, 2015; Melzani et al. 2014b; Sironi & Spitkovsky 2014; Werner et al. 2016; Liu et al. 2015), primarily motivated by the discovery of the spectacular Crab flares (Abdo et al. 2011; Tavani et al. 2011) and the development of the magnetically dominated models. Several simulations have reported hard power-law distributions with spectral index $1\leqslant s\leqslant 2$ when $\sigma \gg 1$ (Guo et al. 2014, 2015; Melzani et al. 2014b; Sironi & Spitkovsky 2014; Werner et al. 2016).⁶  These new simulations claim power-law distributions in the whole reconnection region, suggesting that reconnection in the magnetically dominated regime may be a strong source of nonthermal particles. Guo et al. (2014, 2015) used a force-free current sheet that does not require the hot plasma population in a Harris sheet and showed that the energy spectra of particles within the entire reconnection layer resembles a power law. Through calculating the guiding-center drift motions of all the simulated particles, we have demonstrated that the primary acceleration mechanism is a first-order Fermi mechanism resulting from the particle curvature drift motion along the direction of the $-{\boldsymbol{v}}\times {\boldsymbol{B}}/c$ electric field induced by the relativistic flows. The acceleration leads to the formation of hard power-law spectra $f\propto {(\gamma -1)}^{-s}$ with s ∼ 1 for sufficiently strong magnetic field and large simulation domains. We have developed an analytical model to describe the simulation results and derived a general condition for the formation of the power-law distributions, i.e., the acceleration timescale needs to be shorter than the timescale for particles injected into the reconnection region ${\tau }_{\mathrm{acc}}\lt {\tau }_{\mathrm{inj}}$. The solution may also explain simulations from the Harris current layer (Melzani et al. 2014b; Sironi & Spitkovsky 2014; Werner et al. 2016). The power-law distribution has also been found in nonrelativistic reconnection simulations with a low-β ($\beta =8\pi {nkT}/{B}^{2}\ll 1$) plasma (Li et al. 2015), showing that the mechanism can be applied in a larger parameter range.

Most PIC simulations of relativistic magnetic reconnection focused on pair plasmas. In general, the plasma composition is uncertain and can be either electron–positron pairs or electron–ion plasmas. While proton–electron reconnection has been studied by Melzani et al. (2014a, 2014b) with a limited parameter range (${m}_{i}/{m}_{e}=1-50$, ${\sigma }_{0}={B}^{2}/4\pi ({m}_{i}+{m}_{e}){{nc}}^{2}=1-10$), here we show results from two-dimensional PIC simulations for a magnetically dominated ion–electron plasma, the magnetization for each species is ${\sigma }_{i}\sim {\sigma }_{0}$ and ${\sigma }_{e}\sim ({m}_{i}/{m}_{e}){\sigma }_{0}$. We study the magnetically dominated regime for a range of σ₀ with different ion and electron temperatures and mass ratios ${m}_{i}/{m}_{e}$ up to 1836. We demonstrate that reconnection quickly establishes power-law distributions for both electrons and ions within several light-crossing times. The power-law index is 1 < s < 2 for both electrons and ions in the cases with periodic boundaries. The break energies for electrons and ions are ${\gamma }_{{be}}\sim {\sigma }_{e}$ and ${\gamma }_{{bi}}\sim {\sigma }_{i}$. For the antiparallel case we study, the main acceleration mechanism for particles in the power-law is a Fermi-like acceleration through drift motions of charged particles, consistent with earlier work (Guo et al. 2014, 2015). The numerical methods are presented in Section 2. Section 3 discusses the main results of the paper. In Section 4, we discuss the implications of the study in the context of high-energy astrophysics. Our conclusions are summarized in Section 5.

The simulations start from a magnetically dominated force-free current sheet. Previous studies have shown that the intense current sheets can develop through various processes such as striped wind geometry, field-line foot-point motion, and turbulence cascade (Coroniti 1990; Titov et al. 2003; Makwana et al. 2015). The initial condition of our simulations is a force-free current layer with ${\boldsymbol{B}}={B}_{0}{\rm{tanh}}(z/\lambda )\hat{x}+{B}_{0}{\rm{sech}}(z/\lambda )\hat{y}$(Guo et al. 2014, 2015). This corresponds to a rotating magnetic field with a 180° change in direction within a thickness of 2λ. The initial distributions are spatially uniform with density n₀ and relativistic Maxwellian in energy space. The current density ${\boldsymbol{J}}={{en}}_{0}({{\boldsymbol{U}}}_{i}-{{\boldsymbol{U}}}_{e})$ is represented by a drift velocity ${{\boldsymbol{U}}}_{i}=-{{\boldsymbol{U}}}_{e}$ that satisfies Ampere's law.

The simulations assume an electron–ion plasma with mass ratio ${m}_{i}/{m}_{e}=1\to 1836$ and we mainly focus on the case with ${m}_{i}/{m}_{e}\gg 1$. For simplicity, the initial thermal temperatures for ions and electrons are assumed to be the same and varied in a range of ${T}_{i}={T}_{e}=5\to 20{m}_{e}{c}^{2}$ to examine its effect on the resulting energy spectra. This temperature is similar to the temperature inferred for AGN jets (Homan et al. 2006). The full particle simulations are performed using the VPIC code (Bowers et al. 2009), which explicitly solves Maxwell equations and the relativistic equation of motion for the particles. In the simulations, ${\sigma }_{e}$ is adjusted by changing the ratio of the nonrelativistic electron gyrofrequency ${{\rm{\Omega }}}_{{ce}0}={eB}/({m}_{e}c)$ to the nonrelativistic electron plasma frequency ${\omega }_{{pe}0}=\sqrt{4\pi {{ne}}^{2}/{m}_{e}}$, ${\sigma }_{e}\equiv {B}^{2}/{(4\pi {n}_{e}{m}_{e}{c}^{2})=({{\rm{\Omega }}}_{{ce}0}/{\omega }_{{pe}0})}^{2}$ and the total magnetization is ${\sigma }_{0}\sim ({m}_{e}/{m}_{i}){\sigma }_{e}$. We primarily discuss simulations with ${\sigma }_{e}={10}^{3}\to {10}^{5}$. The cell sizes are chosen so ${\rm{\Delta }}x={\rm{\Delta }}z\lt 0.1{d}_{e}$. Here ${d}_{e}=\sqrt{{\gamma }_{0}}{d}_{e0}$, in which ${\gamma }_{0}=1+3{T}_{e}/2{m}_{e}{c}^{2}$ and ${d}_{e0}=c/{\omega }_{{pe}0}$. All the simulations have 100 electron–ion pairs in each cell. Similar to numerous earlier studies, the boundary conditions in the x-direction are periodic for both fields and particles, and in the z-direction the boundaries are conductive for the field and reflect particles when they reach the boundaries. A long-wavelength perturbation with ${B}_{z}=0.05{B}_{0}$ is included to initiate reconnection. To ensure numerical convergence, it is essential to have enough particles per cell and spatial grids to avoid excessive numerical heating (Guo et al. 2015). We have monitored the numerical energy conservation to ensure that its effect does not influence our results. A more comprehensive study will be reported in another publication (F. Guo et al. 2016, in preparation).

Figure 1 shows the current layer evolution in the simulation with ${\sigma }_{e}={10}^{4}$, ${T}_{i}={T}_{e}=10{m}_{e}{c}^{2}$, and ${m}_{i}/{m}_{e}=100$ so ${\sigma }_{0}\sim 100$. The domain size is ${L}_{x}\times {L}_{z}=4000{d}_{e0}\times 2000{d}_{e0}$. This is the representative run that we will discuss in detail. The left panels and right panels show the current density and velocity in the x-direction at ${\omega }_{{pe}0}t=2000$, 3000, and 5000, respectively. The current sheet thins down under the influence of the perturbation and then breaks into many fast-moving secondary plasmoids due to the growth of secondary tearing instability. The thinning phase lasts until ${\omega }_{{pe}0}t\sim 2000$ and the rest of the simulation is dominated by the evolution of plasmoids. These plasmoids coalesce and merge into a single island on the order of the system size. The whole process is similar to the case with pair plasmas (Daughton & Karimabadi 2007) as long as the simulation domain is sufficiently large to contain several secondary tearing modes. The simulation also shows the development of fast bulk plasma flows with bulk Lorentz factors ${\rm{\Gamma }}=1/\sqrt{1-{(V/c)}^{2}}$, associated with plasmoids (Guo et al. 2015; Liu et al. 2015). As we will discuss below, the fast-moving plasmoids contain most of the accelerated energetic particles, which likely generates high-energy nonthermal emissions with strong beaming effects. We also observed magnetosonic waves generated as the plasmoids interact with each other, similar to recent MHD simulations (Yang et al. 2015).

Zoom In Zoom Out Reset image size

Figure 2 summarizes the plasma energetics for the representative run. In Figure 2(a), the magnetic energy in the reconnecting component drops more than $60\%$ within two to three light crossing times and most of it gets converted into plasma energies. Ions gain slightly more energy than electrons ${E}_{{ki}}/{E}_{{ke}}\sim 1.1$. Increasing ${m}_{i}/{m}_{e}$ to 1836, we find ${E}_{{ki}}/{E}_{{ke}}\sim 1.5$, meaning that, in this case, the energy partition does not strongly depend on the ion-to-electron mass ratio. Figure 2(b) shows the evolution of bulk energy, internal energy, and total kinetic energy. The internal energy is defined as the kinetic energy in the local plasma frame. The bulk energy is energy associated with the plasma bulk motion. The figure shows that the bulk flow can carry a nontrivial portion of the total plasma energy during the energy conversion. Initially it exceeds the internal energy and still contains 5%–10% of the total kinetic energy toward the end of the simulation.

Zoom In Zoom Out Reset image size

Figure 3 presents the energy spectra for both electrons and ions at different simulation times for the representative case. In energy space, electron distribution shows a power-law distribution with s ∼ 1.35, whereas the ion distribution shows a power law with another slope s ∼ 1.2, starting from several times the initial thermal energy. The break energies correspond to ${\gamma }_{{be}}\sim {\sigma }_{e}$ and ${\gamma }_{{bi}}\sim {\sigma }_{i}$ and above that the spectra resemble an exponential cutoff. During the initial phase (${\omega }_{{pe}0}t\lt 3000$), the flat spectrum is associated with the acceleration by the parallel electric field during the thinning process and later the acceleration is dominated by Fermi acceleration associated with plasmoid motion (see discussion below). Figure 3(c) shows the momentum distribution $f(p){p}^{2}$ for both of the species. We find that if the two distributions are plotted in momentum space, they give almost the same slope $f(p)\propto {p}^{-{s}_{p}}$ and ${s}_{p}\sim 3.35$. This is a general feature of Fermi acceleration, but not obviously predicted by other mechanism proposed to explain power laws generated in reconnection regions (Sironi & Spitkovsky 2014; Werner et al. 2016).

Zoom In Zoom Out Reset image size

Figure 4 compares electron energy spectra for simulations with different (a) σ_e, (b) initial temperature, (c) box size, and (d) mass ratio. Results show that for a sufficiently large domain, both electrons and ions get efficiently accelerated and form a power-law distribution. For different σ_e, we observe slightly harder spectra, as the available magnetic energy increases compare to the plasma energy. The trend is similar as we reduce the thermal energy. Panel (c) shows that for a larger simulation domain, the power-law distribution becomes fully developed. As the box size increased to ${L}_{x}=4000{d}_{e0}$ for ${m}_{i}/{m}_{e}=100$, the energy spectra are similar for larger systems. We also find that as long as the simulation domain is sufficiently large, the mass ratio does not significantly influence the electron spectrum. As for the spectral index, we find that the simulations give 1 < s < 2 for the periodic closed system typically used for reconnection studies. This, however, may be changed for the case with open boundary conditions (Daughton et al. 2006), where the particle escape from the boundaries could lead to a softer spectrum (Guo et al. 2014, 2015).

Zoom In Zoom Out Reset image size

Now we briefly discuss the acceleration of protons and electrons using a particle tracking module newly implemented in VPIC. We find that particles initially in the current sheet mainly experience strong direct acceleration through the parallel electric field. Typical examples for both electrons and ions are shown in Figure 5(a). For particles that get accelerated to the highest energies γ_max in the simulation, this is the main type of acceleration. However, for the accelerated particles that form the power-law distribution $\gamma \lt {\gamma }_{b}$, the main acceleration pattern is a Fermi-like process where particles gain energy by bouncing back and forth a few times associated with the motion of plasmoids as the plasmoids evolve and coalesce. In Figure 5(b), we show typical particle trajectories that resemble the Fermi mechanism. The acceleration process for both electrons and ions are very similar. We verified that the size of the largest island is well above the gyroradius of power-law ions with the real mass ratio. In a forthcoming paper, we will discuss the statistics of particle acceleration in detail, similar to previous studies (Guo et al. 2014; Li et al. 2015).

Zoom In Zoom Out Reset image size

It is thought that as relativistic jets or pulsar winds are launched, initially, the energy flow is magnetically dominated and then the magnetic energy gets converted into plasma energy (Spruit 2010). During this process, strong sheared current sheets may form (McKinney & Uzdensky 2012) and magnetic energy gets converted into plasma energy through magnetic reconnection and produces nonthermal emission. Our results imply that the observed nonthermal emission from AGNs, GRBs, and PWNe could come from strong particle acceleration in magnetic reconnection regions. For example, when modeling the emissions from TeV blazars, it is required that the electron spectrum is a power-law distribution with an energy range of ${\gamma }_{e}={10}^{3}-{10}^{5}$ (Böttcher et al. 2013). Recent observations show that strong multiwavelength flares are sometimes accompanied by strong polarization variations, indicating the active participation of the magnetic field during flares (Zhang et al. 2015a). Our calculation shows the development of power-law distribution with break energy ${\gamma }_{{be}}\sim {\sigma }_{e}$ regardless of the ion-to-electron mass ratio, suggesting that the nonthermal electrons can be produced when $1\lesssim {\sigma }_{0}\lesssim 100$ ($1800\lesssim {\sigma }_{e}\lesssim 1.8\times {10}^{5}$) for the real mass ratio. Our simulations show that the spectral index is $1\lt s\lt 2$, which is sometimes reported in TeV blazars (Hayashida et al. 2015; Yan et al. 2015). Most of the energetic particles are contained in plasmoids and the bulk flow Lorentz factor Γ can be a few so we expect that radiation from the nonthermal particles has strong beaming effects that may help explain the fast variability observed in TeV blazars (Giannios et al. 2009).

In GRBs, the spectra of prompt emission suggest that the radiating electrons can reach ${\gamma }_{e}\sim 100$ and form a nonthermal energy distribution. In our simulation, the power-law distribution starts from several times the thermal energy (${\gamma }_{1}$ ∼ 50–200) and extends to high energies (${\gamma }_{2}\sim {\sigma }_{e}$) up to 10⁵. As the small-scale reconnection events continuously occur, electrons can be continually accelerated so the so-called "slow cooling" regime is likely to be valid. The highly efficient acceleration may help to solve the efficiency problem in GRB models and the relativistic bulk motion may explain the fast variability observed in some GRB events (Zhang & Yan 2011; Deng et al. 2015).

Theoretical analysis has shown that the strong particle acceleration and/or the beaming effect in the reconnection region can explain the gamma-ray Crab flares (Abdo et al. 2011; Tavani et al. 2011). Our simulation shows that reconnection is a credible process that drives both strong particle acceleration and relativistic bulk motion. The strong energy conversion from magnetic energy into the plasma kinetic energy during relativistic magnetic reconnection may also help us understand the so-called σ-problem, where the energy conversion is needed for dissipating the Poynting-flux-dominated flow (e.g., Coroniti 1990).

In some hadronic models, a proton spectrum with s < 2 may be required to model TeV blazars (Böttcher et al. 2013). Our results show that ion energy gain is comparable to the electron energy gain and most of the energy goes into a nonthermal power-law spectrum. The maximum energy of protons γ_imax depends on the system size and the largest gyroradius can be comparable to the system size. This suggests that magnetic reconnection can be a viable mechanism for providing energetic protons in hadronic models.

Most observed energy spectra are considerably softer than the simulated results. Moreover, the hard power laws (s < 2) obtained in this and other works are only in a limited energy range. While the maximum energy of particles grows with time, the power-law break energy does not and is limited by the amount of dissipated magnetic energy ${\gamma }_{{be}}\sim ({m}_{i}/{m}_{e}){\gamma }_{{bi}}\sim {[\delta {\sigma }_{e}(2-s)]}^{1/(2-s)}$, where $\delta \sim 0.2$ is the fraction of the dissipated magnetic energy channeled into each species. This poses a serious problem that must be addressed in the future. In our earlier papers (Guo et al. 2014, 2015), we have shown that particle escape from the acceleration region could lead to a softer spectrum. Since our current simulations employ periodic systems, the acceleration tends to give the hardest spectrum. This suggests more realistic open boundary conditions (e.g., Daughton et al. 2006) that include the effect of particle escape from the simulation domain are necessary. Also, our results are tested for a domain size up to $L\sim {10}^{4}{d}_{e0}$, but for larger systems the spectrum could still evolve as acceleration continues.

Astrophysical reconnection sites have long been expected to be strong sources of nonthermal particles. Our current study has demonstrated that relativistic magnetic reconnection is efficient at generating nonthermal particles though the underlying energy conversion and acceleration process in a magnetically dominated proton–electron plasma. During reconnection, the energy distributions for both ions and electrons develop power-law spectra, which is important to understand the nonthermal emission from objects like PWNe, AGNs, and GRBs. We also find that during the active reconnection phase the outflow speeds approach the light speed and the bulk Lorentz factor can reach a few, which may help explain the observed fast variability and high luminosity emissions in several systems. These findings on plasma energization and particle acceleration in reconnection regions further demonstrate the important effect of magnetic reconnection in high-energy astrophysical systems.

We are grateful to Jonathan Jara-Almonte, who developed an initial version of the particle tracking module. This work is supported by the DOE through the LDRD program at LANL and DOE/OFES support to LANL in collaboration with CMSO. X.L. is supported by NASA Headquarters under the NASA Earth and Space Science Fellowship Program-Grant NNX13AM30H. Contributions from W.D. were supported by NASA from the Heliophysics Theory Program. N.L.-R. is supported by the M. Hildred Blewett Fellowship of the APS. Simulation resources are provided by LANL institutional computing.