Particle Acceleration and Plasma Dynamics during Magnetic Reconnection in the Magnetically-dominated Regime
DDraft version September 12, 2018
Preprint typeset using L A TEX style emulateapj v. 5/2/11
PARTICLE ACCELERATION AND PLASMA DYNAMICS DURING MAGNETIC RECONNECTION IN THEMAGNETICALLY-DOMINATED REGIME
Fan Guo , Yi-Hsin Liu , William Daughton , and Hui Li Draft version September 12, 2018
ABSTRACTMagnetic reconnection is thought to be the driver for many explosive phenomena in the universe.The energy release and particle acceleration during reconnection have been proposed as a mechanismfor producing high-energy emissions and cosmic rays. We carry out two- and three-dimensional kineticsimulations to investigate relativistic magnetic reconnection and the associated particle acceleration.The simulations focus on electron-positron plasmas starting with a magnetically dominated, force-freecurrent sheet ( σ ≡ B / (4 πn e m e c ) (cid:29) f ∝ ( γ − − p and approaches p = 1 forsufficiently large σ and system size. Eventually most of the available magnetic free energy is convertedinto nonthermal particle kinetic energy. An analytic model is presented to explain the key resultsand predict a general condition for the formation of power-law distributions. The development ofreconnection in these regimes leads to relativistic inflow and outflow speeds and enhanced reconnectionrates relative to non-relativistic regimes. In the three-dimensional simulation, the interplay betweensecondary kink and tearing instabilities leads to strong magnetic turbulence, but does not significantlychange the energy conversion, reconnection rate, or particle acceleration. This study suggests thatrelativistic reconnection sites are strong sources of nonthermal particles, which may have importantimplications to a variety of high-energy astrophysical problems. Subject headings: acceleration of particles – magnetic reconnection – relativistic process – gamma-raybursts: general – galaxies: jets – pulsars: general INTRODUCTIONMagnetic reconnection is a fundamental plasma pro-cess that rapidly rearranges magnetic topology and con-verts magnetic energy into various forms of plasma ki-netic energy, including bulk plasma flow, thermal andnonthermal plasma distributions (Kulsrud 1998; Priest& Forbes 2000). It is thought to play an important roleduring explosive energy release processes of a wide va-riety of laboratory, space and astrophysical systems in-cluding tokamak, planetary magnetospheres, solar flares,and high-energy astrophysical objects. Relativistic mag-netic reconnection is often invoked to explain high-energyemissions and ultra-high-energy cosmic rays from ob-jects such as pulsar wind nebulae (PWNe; Kirk 2004;Arons 2012; Uzdensky & Spitkovsky 2014), jets from ac-tive galactic nuclei (AGN; de Gouveia dal Pino & Lazar-ian 2005; Giannios et al. 2009), and gamma-ray bursts(GRBs; Thompson 1994; Zhang & Yan 2011; McKinney& Uzdensky 2012). In those systems, the magnetizationparameter σ ≡ B / (4 πn e m e c ), is often estimated to bemuch larger than unity σ (cid:29) e n speedapproaches the speed of light v A ∼ c . To explain theobserved high-energy emissions, often an efficient energyconversion mechanism is required (e.g., Zhang et al. 2007;Celotti & Ghisellini 2008; Zhang & Yan 2011; Zhang et al.2013). Collisionless shocks, which can efficiently convert [email protected] Los Alamos National Laboratory, Los Alamos, NM 87545,USA NASA Goddard Space Flight Center, Greenbelt, MD 20771,USA plasma flow energy into thermal and nonthermal ener-gies in low- σ flows, are inefficient in dissipating magneti-cally dominated flows, where most of the energy is storedin magnetic fields. In these regimes, magnetic reconnec-tion is the primary candidate for dissipating and convert-ing magnetic energy into relativistic particles and subse-quent radiation. Understanding magnetic reconnectionis also important to solve the so-called σ -problem (Coro-niti 1990; Lyubarsky & Kirk 2001; Kirk & Skjæraasen2003; Porth et al. 2013), where strong magnetic dissipa-tion may be required to convert the magnetically domi-nated flow ( σ (cid:29)
1) to a matter-dominated flow ( σ (cid:28) σ (cid:29) u out approaches the speed of light, and therate of relativistic magnetic reconnection and inflow ve-locity u in may increase compared to the nonrelativisticcase. This is because of the enhanced outflow densityarising from the Lorentz contraction of plasma passingthrough the diffusion region u in ∝ u out Γ out / Γ in , whereΓ out and Γ in are Lorentz factors of outflows and in-flows, respectively. However, later analysis (Lyubarsky2005) showed that for a pressure-balanced current layer( B / π ∼ nk ( T i + T e )), the thermal pressure constrainsthe outflow speed to be mildly relativistic and hence the a r X i v : . [ a s t r o - ph . H E ] A p r Guo et al.effect of Lorentz contraction is negligible. Although therate of relativistic magnetic reconnection is reported toenhance in a number of studies using different numeri-cal models (Zenitani et al. 2009; Bessho & Bhattacharjee2012; Takamoto 2013; Sironi & Spitkovsky 2014; Guoet al. 2014; Melzani et al. 2014a), its nature is not clear.This issue is recently revisited by carefully analyzing re-sults from fully kinetic two-dimensional (2D) simulations(Liu et al. 2015), which shows that the plasma densityand pressure around the X-line drop significantly as theinitial high-pressure region is depleted during reconnec-tion. This results in a reconnection region with σ (cid:29) v in ∼ c . The local recon-nection rate across the diffusion region is well predictedby a simple model that includes the Lorentz contrac-tion. However, the extension of these results to three-dimensional (3D) kinetic simulations was not considered.Plasma energization during magnetic reconnection hasbeen extensively discussed in the literature. However,the primary acceleration mechanism is still unclear. Ro-manova & Lovelace (1992) analyzed particle motions in alarge-scale reconnection region and predicted a spectrum dN/dγ = γ − p with p = 1 . p = 2 when particles are accelerated in a direct elec-tric field associated with magnetic reconnection. Using amodel for the motions of particles in a steady magneticreconnection region, Larrabee et al. (2003) have foundstrong particle acceleration in the reconnection layer andobtained a hard energy spectrum with a spectral indexof about p = 1. The first-order Fermi acceleration inconverging reconnection inflows has been discussed (deGouveia dal Pino & Lazarian 2005; Lazarian & Opher2009; Kowal et al. 2012). Drury (2012) studied the accel-eration in a reconnection layer including energy changein both inflows and outflows and show fluid compressionis crucial for efficient particle acceleration. Test-particleapproach has been applied to interpret the strong parti-cle acceleration responsible for γ -ray flares from the Crabpulsar (Cerutti et al. 2012). The results show that mag-netic reconnection may be the site of extreme particleacceleration required to explain high energy emissionsfrom the Crab flares (see also Cerutti et al. 2013). Self-consistent kinetic simulations have been widely used tostudy plasma dynamics and particle energization duringmagnetic reconnection. Most of previous kinetic studieshave focused on the regime with σ (cid:46)
1, and found a num-ber of acceleration mechanisms such as direct accelera-tion at X-line regions (Drake et al. 2005; Fu et al. 2006;Pritchett 2006; Oka et al. 2010; Huang et al. 2010) andFermi-type acceleration in reconnection induced plasmaflows within magnetic islands (Drake et al. 2006, 2010;Oka et al. 2010; Huang et al. 2010). The high- σ regime( σ >
1) has been explored in a number of papers usingthe Harris equilibrium (Zenitani & Hoshino 2001, 2007;Liu et al. 2011; Bessho & Bhattacharjee 2012; Ceruttiet al. 2013; Sironi & Spitkovsky 2014; Melzani et al.2014b; Werner et al. 2014). However, the initial con-dition employed in these studies requires a hot plasmacomponent inside the current sheet to maintain forcebalance, which may not be justified for high- σ plasmas.Recently, several studies have reported hard power-lawdistributions 1 ≤ p ≤ σ (cid:29) f ∝ ( γ − − p withspectral index approaching p = 1 for a sufficiently high σ and a large system size. An analytical model was de-veloped to describe the main feature of the simulationsand it gives a general condition for the formation of thepower-law particle energy distribution. The solution alsoappears to explain simulations from the Harris currentlayer, in which the particles initially in the current layerform a heated thermal distribution and particles injectedfrom the upstream region are accelerated into a power-law distribution (Sironi & Spitkovsky 2014; Melzani et al.2014b; Werner et al. 2014).Another important issue is the influence of 3D dynam-ics that may significantly modify the reconnection rate,energy release, and particle acceleration process. Re-cently, the rate of 3D nonrelativistic magnetic reconnec-tion has been explored and compared with 2D simula-tions in a number of non-relativistic studies (Liu et al.2013; Daughton et al. 2014), which showed only mod-est differences between 2D and 3D simulations althoughstrong 3D effects emerge as the tearing mode developsover a range of oblique angles (Daughton et al. 2011). Forrelativistic magnetic reconnection with a pair plasma,Sironi & Spitkovsky (2014) reported a factor of fourdecrease of reconnection rate for 3D simulations com-pared to 2D simulations. This is in contrast to Guoet al. (2014), who observed similar reconnection rate andenergy conversion between 2D and 3D simulations, al-though the kink mode (Daughton 1999) strongly inter-acts with the tearing mode leading to a turbulent recon-nection layer (Yin et al. 2008). For particle accelerationin 3D reconnection simulations, early studies reportedthat the drift kink instability can modify the electric andmagnetic field structures in an antiparallel reconnectionlayer and prohibit nonthermal acceleration (Zenitani &Hoshino 2005b, 2007, 2008). However, recent large-scale3D simulations have found strong nonthermal particlespectra are produced even when the kink mode is active(Liu et al. 2011; Sironi & Spitkovsky 2014; Guo et al.2014). Earlier large-scale 3D studies have shown the de-velopment of turbulence in the reconnection layer (Yinet al. 2008), but its effect to energetic particle acceler-ation is unknown. It is therefore important to furtherstudy particle acceleration in relativistic regimes usinglarge-scale 3D kinetic simulations.In this paper, we perform 2D and 3D fully kinetic simu-lations starting from a force-free current sheet with uni-form plasma density and temperature to model recon-nection over a broad range in the magnetization param-eter σ = 0 .
25 - 1600. This paper builds upon earlierwork (Guo et al. 2014) and gives further details regard-article Acceleration in Relativistic Magnetic Reconnection 3ing the plasma dynamics and particle acceleration duringrelativistic magnetic reconnection in the high- σ regime.We also present detailed results from a 3D simulationthat shows a turbulent reconnection layer arising fromthe interaction between the secondary tearing and kinkmodes. In Section 2, we describe the numerical methodsand parameters. Section 3 discusses the main results ofthe paper. In section 4, we present an analytical modelthat explains the main feature of particle acceleration inthe simulations. The implications from this work for arange of astrophysical problems are discussed in Section5 and our conclusions are summarized in Section 6. Inaddition, we have also explicitly examined the numericalconvergence for this problem and the effect of numeri-cal heating in our simulations, which is discussed in theAppendix. NUMERICAL METHODSWe envision a situation where intense current sheetsare developed within a magnetically dominated plasma.Earlier work in non-relativistic low- β plasmas has shownthat the gradual evolution of the magnetic field can leadto formation of intense nearly force-free current layerswhere magnetic reconnection may be triggered (Titovet al. 2003; Galsgaard et al. 2003). In the presentstudy, the critical parameter is the magnetization pa-rameter defined as σ ≡ B / (4 πn e m e c ), which roughlycorresponds to the available magnetic energy per par-ticle. The numerical simulations presented in this pa-per are initialized from a force-free current layer with B = B tanh( z/λ )ˆ x + B sech( z/λ )ˆ y (Che et al. 2011;Liu et al. 2013, 2014), which corresponds to a mag-netic field with magnitude B rotating by 180 ◦ acrossthe central layer with a half-thickness of λ . No exter-nal guide field is included in this study but there isan intrinsic guide field B y associated with the centralsheet. The plasma consists of electron-positron pairswith mass ratio m i /m e = 1. The initial distributionsare Maxwellian with a spatially uniform density n anda thermal temperature ( kT i = kT e = 0 . m e c ). Par-ticles in the central sheet have a net drift U i = − U e to represent a current density J = en ( U i − U e ) thatis consistent with ∇ × B = 4 π J /c . Since the force-freecurrent sheet does not require a hot plasma componentto balance the Lorentz force, this initial setup is moresuitable to study reconnection in low β and/or high- σ plasmas. The full particle simulations are performed us-ing the VPIC code (Bowers et al. 2009) and NPIC code(Daughton et al. 2006; Daughton & Karimabadi 2007),both of which solve Maxwell equations and push particlesusing relativistic approaches. The VPIC code directlysolves electric and magnetic fields in Maxwell equations,whereas in the NPIC code the fields are advanced us-ing the scalar and vector potentials. Therefore the twocodes have different algorithms and can be used to cross-check numerical results for this problem. We find thatthey give consistent results for this problem. In addition,we have developed a particle-tracking module to analyzethe detailed physics of the particle energization process.In the simulations, we define and adjust σ by changingthe ratio of the electron gyrofrequency Ω ce = eB/ ( m e c )to the electron plasma frequency ω pe = (cid:112) πne /m e , σ ≡ B / (4 πn e m e c ) = (Ω ce /ω pe ) . For 2D simulations, we have performed simulations with σ = 0 . → L x × L z = 300 d i × d i , 600 d i × d i ,and 1200 d i × d i , where d i is the inertial length c/ω pe .For 3D simulations, the largest case is L x × L y × L z =300 d i × d i × d i with σ = 100. For high- σ cases( σ > x = ∆ y = 1 . / √ σd i and ∆ z = 0 . / √ σd i , so the particle gyromotion scale ∼ v the / √ σd i is resolved. The time step is chosen tocorrespond to a Courant number C r = c ∆ t/ ∆ r = 0 . r = ∆ x ∆ y ∆ z/ (∆ x ∆ y + ∆ y ∆ z + ∆ x ∆ z ). Thehalf-thickness of the current sheet is λ = 6 d i for σ ≤ d i for σ = 400, and 24 d i for σ = 1600 in order to sat-isfy the drift velocity U i < c . For both 2D and 3D sim-ulations, we have more than 100 electron-positron pairsin each cell. The boundary conditions for 2D simula-tions are periodic for both fields and particles in the x -direction, while in the z -direction the boundaries areconducting for the field and reflecting for the particles.In the 3D simulations, the boundary conditions are peri-odic for both fields and particles in the y -direction, whilethe boundary conditions in the x and z directions are thesame as the 2D cases. A weak long-wavelength pertur-bation (Birn et al. 2001) with B z = 0 . B is included toinitiate reconnection. The parameters for different runsare summarized in Table 1, which also lists key resultssuch as maximum energy of particles, spectral index, thefraction of kinetic energy converted from the magneticenergy and the portion of energy gain arising from theperpendicular electric fields.Using the set of numerical parameters described above,all of the simulations show excellent energy conservationwith violation of energy conservation less than 10 − ofthe total energy in all cases. However, we note that toaccurately determine the particle energy spectra, the vi-olation in energy conservation should be smaller than theinitial plasma kinetic energy, which is only a small frac-tion of the total energy for the problem we study. Cau-tion is needed when using a small number of particlesper cell and a small initial plasma kinetic energy in thesimulations (Sironi & Spitkovsky 2014), since numericalheating may significantly modify the particle distribu-tion. In the Appendix, we have extensively tested howthe numerical convergence varies with the initial plasmatemperature, cell size, number of particles per cell, andtime step. For all the cases we present in the main paper,the violation of energy conservation is a few percent ofthe initial kinetic energy in the system, meaning effectssuch as numerical heating have a negligible influence onthe simulated energy spectra. SIMULATION RESULTS3.1.
General feature and energy conversion
Figure 1 gives an overview of the evolution of the cur-rent layer in the case with σ = 100 and domain size L x × L z = 300 d i × d i ( L y = 300 d i for the 3D sim-ulation) from runs 2D-7 and 3D-7. Panel (a) shows thecolor-coded current density from the 2D simulation andPanel (b) shows a 2D cut of the current density anda 3D isosurface of plasma density colored by the cur-rent density from the 3D simulation at ω pe t = 175 and ω pe t = 375, respectively. Starting from the initial per-turbation, the current sheet gradually narrows as thecurrent density is concentrated in the central region. Guo et al.In the 2D simulation, the extended thin current sheetbreaks into a number of fast-moving secondary plasmoids( ω pe t ∼ E B , electric field energy E E , kinetic energy E k , and energy carried by relativisticparticles with γ > − of the initial value.The evolutions of different forms of energies between 2Dand 3D simulations are very similar. In both the 2Dand 3D simulations, about 25% of the magnetic energyis converted into plasma kinetic energy, most of whichis carried by relativistic particles. Figure 2 (b) showsthe time-integrated energy conversion from magnetic en-ergy into plasma energy in the simulation (cid:82) t dt (cid:82) dV J · E and its contribution from parallel and perpendicular elec-tric field terms J (cid:107) · E (cid:107) and J ⊥ · E ⊥ , respectively. Here (cid:82) dV = (cid:82) dxdydz . The difference in energy conversionbetween the 2D and 3D simulations can be as large as afactor of two at ω pe t = 300, but at the end of the simula-tions both cases have converted about the same amountof magnetic energy. This shows that the kink instabil-ity that may modify the magnetic field does not signif-icantly change the overall energy conversion. While theenergy conversion through parallel electric field is impor-tant when the thin current layer initially develops, mostof the energy conversion is due to perpendicular electricfields induced by relativistic flows as the system is dom-inated by secondary plasmoids/flux ropes. This analysishas been done in all the cases and summarized in Table1, which shows that in most of cases, the perpendicularelectric field plays a dominant role in converting mag-netic energy into plasma kinetic energy. This can alsobe seen in Figure 3, which shows the color-coded inten-sities of J · E , J ⊥ · E ⊥ and J (cid:107) · E (cid:107) from the 2D and 3Dsimulations at ω pe t = 175 and ω pe t = 375, respectively.Figure 2 (c) compares the energy spectra from the 2Dand 3D simulations at various times. The most strikingfeature is that a hard power-law spectrum f ∝ ( γ − − p with a spectral index p ∼ .
35 forms in both 2D and3D runs. Although a fraction of particles are acceler-ated in the early phase when the parallel electric fieldis important, most of the particles in the power-law dis-tribution are accelerated when the system is dominatedby plasmoids/flux ropes. As we will discuss below, theformation of power law is closely related to the motionalelectric field induced by the fast moving plasmoids. Inthe subpanel, the energy spectrum for all particles in the 3D simulation at ω pe t = 700 is shown by the red line.The low-energy portion can be fitted by a Maxwelliandistribution (black) and the nonthermal part resemblesa power-law distribution (blue) starting at γ ∼ γ (cid:38) ∼
25% of particles and ∼
95% of the kineticenergy. The maximum particle energy of the system canbe predicted approximately using the reconnecting elec-tric field m e c ( γ max −
1) = (cid:82) | qE rec | cdt until the gyro-radius is comparable to the system size (see also Figure6b). Although we observe a strong kink instability inthe 3D simulations, the energy conversion and particleenergy spectra are remarkably similar to the 2D results,indicating the 3D effects are not crucial for the particleacceleration. The fast acceleration is distinct from thatof nonrelativistic magnetic reconnection, where particlesare at most accelerated to mildly relativistic energy (e.g.,Fu et al. 2006; Drake et al. 2006; Pritchett 2006; Okaet al. 2010). The nonthermal-dominated distribution inthe simulations is also quite different from distributionsin the relativistic shock regions (e.g., Spitkovsky 2008),where the particles are heated at the shock front and forman extended thermal distribution containing most of thedissipated energy. The power-law spectral index p ∼ p ∼ Particle Acceleration
We now discuss the details of particle acceleration. Wewill first present some analysis of particle trajectoriesto show the acceleration mechanism. Then the domi-nant acceleration mechanism is distinguished by track-ing all the particles and calculating the energy gain us-ing the guiding-center drift approximation. The resultsdemonstrate that the dominant acceleration mechanismis a first-order Fermi acceleration through curvature driftmotion along the motional electric field induced by therelativistic reconnection flows. We calculate the acceler-ation rate α = ∆ ε/ ( ε ∆ t ) and its time integral for caseswith σ = 6 − ε is the averaged energy gainfor particles of energy ε over a period ∆ t . Finally, wesummarize the character of the energy spectra. Thesemain results will be discussed and interpreted in detailin Section 4, where we present the acceleration model.Figure 4 and Figure 5 present the trajectory analysisfor the motions of accelerated particles in the 2D casewith σ = 100 and L x × L z = 600 d i × d i . Theseparticles are selected to show the common accelerationpattern of accelerated particles, which is consistent withthe results of statistical analysis for the acceleration ofparticles in Figure 6. The first three panels of Figure 4show (a) the trajectory of a representative particle closeto the central sheet between ω pe t = 30 - 300 togetherwith E (cid:107) at ω pe t = 180, (b) the trajectory of the sameparticle between ω pe t = 310 - 510 together with E y at ω pe t = 400, and (c) the trajectory of the particle between ω pe t = 510 - 720 together with E y at ω pe t = 640, respec-tively. The starting and ending locations of the particleare labeled by ‘+’ and ‘ × ’ signs, respectively. Note thatthe field is highly variable in time and the location ofthe particle at the same time step as the field contouris drawn by the ‘ ∗ ’ sign. The two bottom panels showarticle Acceleration in Relativistic Magnetic Reconnection 5 (a)(b) i i xz i i i xy i i i i xzy z |J|/J ω pe t=375 ω pe t=175 ω pe t=175 ω pe t=375 Fig. 1.—
Evolution of 2D and 3D simulations with σ = 100 and domain size L x × L z = 300 d i × d i ( L y = 300 d i for 3D); (a) Color-codedcurrent density from the 2D simulation at ω pe t = 175 and ω pe t = 375, respectively; (b) 2D cut of current density and a 3D isosurface ofthe plasma density colored by the current density at ω pe t = 175 and ω pe t = 375, respectively. the evolution of the particle energy as a function of time(d) and energy as a function of the x position (e), re-spectively. Each period corresponding to that in ( a )-( c )is labeled by the same color. The green curve representsthe energy gain in the parallel electric field integratedfrom t = 0. Initially the particle is close to the centrallayer and gains energy by the parallel electric field. It isthen strongly accelerated by perpendicular electric fieldwhen the reconnection region breaks into multiple islandsand the electric field is mostly the motional electric field E = − V × B /c generated by relativistic plasma outflows.The figure also shows that the acceleration by E ⊥ resem-bles a Fermi process by bouncing back and forth withina magnetic island.Figure 5 presents another view of the particle accel-eration physics. It is similar to Figure 4, but the fieldcontours show the outflow speed to highlight the role of V x in the particle’s energization. This clearly illustrates arelativistic first-order Fermi process by bouncing in out-flow regions of the reconnection layer. Note the energygain from the parallel electric field for this sample par-ticle is negligible since it entered the reconnection layerlonger after the development of multiple plasmoids.In Figure 6, we present more analysis for the mecha-nism of particle acceleration. Panel (a) shows the energyas a function of the x -position of four accelerated parti-cles. Similar to Figure 5, the electrons gain energy bybouncing back and forth within the reconnection layer. We have analyzed trajectories of a large number of par-ticles and found the energy gain for each cycle is ∆ ε ∼ ε ,which demonstrates that the acceleration mechanism is afirst-order Fermi process (Drake et al. 2006, 2010; Kowalet al. 2011). Panel (b) shows the maximum particle en-ergy in the system as a function of time. This is plot-ted using different count level from the 1-particle levelto the 1000-particle level. Also plotted is the estimatedmaximum energy resulting from the reconnecting electricfield by assuming particles moving along the electric fieldat the speed of light (cid:82) | qE rec | cdt . This shows that themaximum possible energy occurs for a small number ofparticles that continuously sample the reconnection elec-tric field m e c γ max = (cid:82) | qE rec | cdt . At late time, as theparticle gyroradius becomes large and comparable to thesystem size, the maximum energy saturates. To showthe Fermi process more rigorously, we have tracked theenergy change for all the particles in the simulation andthe relative contributions arising from the parallel elec-tric field ( m e c ∆ γ = (cid:82) qv (cid:107) E (cid:107) dt ) and curvature drift ac-celeration ( m e c ∆ γ = (cid:82) q v curv · E ⊥ dt ) similar to (Dahlinet al. 2014), where v curv = γv (cid:107) ( b × ( b · ∇ ) b ) / Ω ce , v (cid:107) isthe particle velocity parallel to the magnetic field, and b = B / | B | . Panel (c) shows the averaged energy gainand the contribution from parallel electric field and cur-vature drift acceleration over an interval of 25 ω − pe as afunction of energy starting at ω pe t = 350. The energygain follows ∆ ε ∼ αε , confirming the first-order Fermi Guo et al. (a)(b) ε / ε t o t a l ω pe t (c) γ -1 f -4 -2 -8 thermal all particlesnonthermal -1 -2 -3 -4 ε / ε t o t a l ω pe t -2 -4 -6 -8 -2 -1 ∫ dt ∫ dVJ · E2D ∫ dt ∫ dV(J · E) ⊥ ∫ dt ∫ dV(J · E) || ∫ dt ∫ dVJ · E3D ∫ dt ∫ dV(J · E) ⊥ ∫ dt ∫ dV(J · E) || Fig. 2.—
Plasma energetics in 2D and 3D simulations with σ = 100 and domain size L x × L z = 300 d i × d i ( L y = 300 d i for3D); (a) Evolution of magnetic energy E B , electric field energy E E ,plasma kinetic energy E k and energy carried by relativistic parti-cles with Lorentz factor γ >
4; (b) Energy conversion from mag-netic energy into plasma energy integrated over time (cid:82) t dt (cid:82) dV J · E and its contribution from parallel and perpendicular electric field J (cid:107) · E (cid:107) and J ⊥ · E ⊥ ; (c) Evolution of particle energy spectra from2D and 3D simulations. Subpanel: Energy spectrum from the 3Dsimulations at ω pe t = 700. The low energy is fitted with a thermaldistribution and rest of the distribution is a nonthermal power lawwith an exponential cutoff. process identified from particle trajectories. The energygain from the parallel motion depends weakly on energy,whereas the energy gain from the curvature drift accel-eration is roughly proportional to energy. In the earlyphase, the parallel electric field is strong but only accel-erates a small portion of particles, and the curvature driftdominates the acceleration starting at about ω pe t = 250.The contribution from the gradient drift was also evalu-ated and found to be negligible in comparison. Panel (d)shows α = < ∆ ε > / ( ε ∆ t ) measured directly from theenergy gain of the particles in the perpendicular elec-tric field ( m e c ∆ γ = (cid:82) q v ⊥ · E ⊥ dt ) and estimated fromthe expression for the curvature drift acceleration. Theclose agreement demonstrates that curvature drift termdominates the particle energization.For higher σ and larger domains, the acceleration isstronger and reconnection is sustained over a longer du-ration. In Figure 7(a), we present the energy spectra atthe end of simulation for a number of cases with different σ and system size L x × L z = 600 d i × d i . A summaryfor the spectral index can be found in Table 1. In Fig-ure 7(b), a summary for the measured spectral index forthe power-law ranges of all the 2D runs shows that thespectrum is harder for higher σ and larger domain sizes,and approaches the limit p = 1. Note that the spec-tral indexes appear systematically harder than in otherrecent papers (Sironi & Spitkovsky 2014; Melzani et al.2014b; Werner et al. 2014). However, the energy spectrain these studies are plotted using total relativistic en-ergy γmc and here we use kinetic energy ( γ − mc .Using total relativistic energy in the energy spectra sig-nificantly distorts the spectral index in the energy rangeof 0 < γ − <
10, which may alter the interpretationof the results (Sironi & Spitkovsky 2014; Melzani et al.2014b; Werner et al. 2014) ‡ .3.3. Reconnection rate and relativistic flows
Figure 8 (a) shows the time-dependent reconnectionrates normalized using the initial asymptotic magneticfield B in 2D and 3D simulation with σ = 100 (Run2D-7 and 3D-7). The 2D reconnection rate is computedfrom R = E rec B = 1 B V A < ∂ψ∂t >, where ψ = max( A y ) − min( A y ) along the central layer z = 0, A y is the vector potential along the y direction, <> represents a time average over δtω pe = 25 (Liu et al.2014), V A = v A / (cid:112) v A /c ) = (cid:112) σ/ (2 + σ ) c is therelativistic Alfven speed in the cold-plasma limit. Here v A = B / (cid:112) πn ( m i + m e ) is the non-relativistic Alfvenspeed based on B . The 3D reconnection rate is esti-mated by using the mixing of plasma across the sepa-ratrix surfaces (Daughton et al. 2014). The rate in the ‡ In fact, our simulation results show that the “ −
1” spectra canbe obtained as long as the magnetic energy dominates over theinitial plasma kinetic energy 8 πnkT /B = β (cid:28)
1. An examplecan be seen in the Appendix (Figure 13), which robustly shows the p = 1 spectrum can be obtained when σ = 25 and kT = 0 . mc .The same spectrum gives a “ p ∼
2” slope when it is plotted as afunction of γ , which may explain the different conclusions reportedby other papers (Sironi & Spitkovsky 2014; Melzani et al. 2014b;Werner et al. 2014) article Acceleration in Relativistic Magnetic Reconnection 7 (a)(c) i i xz i i i i i i i J·EJ·E (J·E) ⊥ (J·E) || (J·E) ⊥ (J·E) || i i i i i i i i i i i i i i (b) i i i xz ω pe t=175 ω pe t=375(d) 2D2D3D ω pe t=175 ω pe t=3753D J · E ( J · E ) || ( J · E ) ⊥ ( J · E ) || J · E ( J · E ) ⊥ Fig. 3.—
Color-coded intensity of energy conversion rate J · E normalized using n m e c ω pe and contributions from J ⊥ · E ⊥ and J (cid:107) · E (cid:107) for the 2D and 3D simulations with σ = 100 at ω pe t = 175 and ω pe t = 375, respectively. In the early stage the conversion by parallelelectric field is important and the perpendicular electric field plays a dominant role when multiple-plasmoids (flux ropes in 3D) developdue to the secondary tearing instability.
2D simulation is quite variable but the range is within afactor of two times of the 3D results, meaning that the2D and 3D simulations give roughly the same reconnec-tion rate. Figure 8 (b) shows the peak reconnection ratefor a number of 2D cases with σ from 0 .
25 to 1600 andbox size 1200 d i × d i . The rate is observed to increasewith σ from E rec ∼ . B for σ = 1 to E rec ∼ . B for σ = 1600. It shows that the peak reconnection fieldincreases with σ and starts to saturate around σ = 1000.For low- σ cases with σ <
1, the reconnecting electricfield is consistent with previous work for nonrelativisticreconnection (e.g., Daughton & Karimabadi 2007). Moredetailed analyses have shown that for high- σ cases, thereconnection rate normalized using the magnetic field B u upstream of the diffusion region E rec /B u is close to 1 for σ (cid:38)
100 (Liu et al. 2015).In Figure 9 (a) and (b) we plot the maximum flowvelocity in the x direction (outflow direction) and thecorresponding Lorentz factor Γ x . The 2D results arerepresented by blue symbols and the 3D results are inred symbols, respectively. Although we have only used a small simulation domain that may be affected by counter-streaming particles, a relativistic outflow still developswith Γ x of a few. In Figure 9 (c) and (d) we plot the max-imum flow velocity in the z direction (inflow direction)and the corresponding Lorentz factor Γ z , respectively.Interestingly, the inflow speed can also be relativistic forhigh- σ cases. Detailed analysis for the diffusion regionhas been discussed in Liu et al. (2015), which shows thatthe inflow speed can be predicted by a model based onthe Lorentz contraction of the plasma passing throughthe diffusion region.The enhanced reconnection rate and development ofrelativistic inflow/outflow structures are in contrast tothe results reported earlier (Sironi & Spitkovsky 2014),where the outflow can only be mildly relativistic and theinflow speed remains nonrelativistic. Note that Liu et al.(2015) has also reported the development of relativis-tic inflow for both Harris and force-free current sheets,indicating that this property of relativistic magnetic re-connection does not strongly depend on the initial setup. Guo et al. -0.800-0.533-0.2670.0000.2670.5330.800 x/d i
360 400 440 x/d i ω pe t γ - z / d i z / d i z / d i -10.00-6.67-3.330.003.336.6710.00 -25.00-16.67-8.330.008.3316.6725.00 (a) E || ω pe t = 180Particle ω pe t = 30 - 310(b) E y ω pe t = 400Particle ω pe t = 310 - 510(c) E y ω pe t = 640Particle ω pe t = 510 - 720 (d) (e) -25.00-16.67-8.330.008.3316.6725.00-25.00-16.67-8.330.008.3316.6725.00 Fig. 4.—
Panels (a)-(c) show a particle trajectory in the x - z plane together with the color-coded electric field (a) E (cid:107) , (b) E y , and (c) E y . Panels (d) and (e) show the particle energy as a function of time and energy as a function of the x position, respectively. In (d) and(e), curves with different colors represent the energy evolution during time periods in (a)-(c). The green curve shows the integrated energygain from the parallel electric field.
3D Dynamics
In our three-dimensional simulation, we also findstrong bulk Γ x ∼ y direction and a volume rendering of thecurrent density in the 3D simulation with σ = 100 at ω pe t = 708. The power spectrum shows a clear iner-tial range with a slope of “ −
2” and steeper slope forhigher wave numbers k ⊥ d i (cid:38)
1. As we have discussed,the 3D simulation allows the development and interac-tion of secondary tearing instability and kink instability,leading to a turbulent magnetic field in the reconnectionlayer. For the whole simulation, the range of scales for the 2D magnetic islands is similar to the observed 3Dflux ropes. The maximum energy in both 2D and 3Dagrees well with the time integral of energy gain fromreconnecting electric field. This is in contrast to the ear-lier kinetic simulations (Zenitani & Hoshino 2005a, 2007,2008). The energy distributions reported in this paperare remarkably similar in 2D and 3D, suggesting that theunderlying Fermi acceleration is rather robust and doesnot depend on the existence of well-defined magnetic is-lands. A SIMPLE MODELIt is often argued that some loss mechanism is neededto form a power-law distribution (Zenitani & Hoshino2001; Drake et al. 2010, 2013; Hoshino 2012). How-article Acceleration in Relativistic Magnetic Reconnection 9 -1.300-0.867-0.4330.0000.4330.8671.300 x/d i
400 450 500 x/d i ω pe t γ - z / d i z / d i z / d i -1.250-0.833-0.4170.0000.4170.8331.250 -1.000-0.667-0.3330.0000.3330.6671.000 -1.250-0.833-0.4170.0000.4170.8331.250 -1.000-0.667-0.3330.0000.3330.6671.000-1.000-0.667-0.3330.0000.3330.6671.000-1.000-0.667-0.3330.0000.3330.6671.000 z / d i (a) V x ω pe t=520 -25.00-16.67-8.330.008.3316.6725.00-25.00-16.67-8.330.008.3316.6725.00 Particle ω pe t = 480 - 600(b) V x ω pe t=670 Particle ω pe t = 600 - 770(c) V x ω pe t=870 Particle ω pe t = 770 - 960(d) V x ω pe t=1020 Particle ω pe t = 960 - 1120 (e) (f) Fig. 5.—
Panels (a)-(d) show a particle trajectory in the x - z plane together with the fluid velocity in the x direction V x . Panels (e) and(f) show the particle energy as a function of time and energy as a function of the x position, respectively. Different colored curves representthe energy evolution during time periods in (a)-(d), showing that the particle gains energy by bouncing in the relativistic flow generatedby reconnection. ever, the simulation results reported in this paper clearlyshow power-law distributions in a closed periodic sys-tem. Here we present a simple model to explain thepower-law energy spectrum observed in our PIC simu-lations. The model is illustrated by Figure 11(a). Asreconnection proceeds, cold plasma in the upstream re-gion advects into the acceleration zone at a constant ve-locity that is determined by reconnection electric field V in = c E rec × B /B . The process lasts τ ∼ L z / V in ,where L z is the size of the simulation box along the z direction. In the acceleration region, our analysis hasshown that a first-order Fermi process dominates theenergy gain during reconnection. We solve the energy-continuity equation for the energy distribution function f ( ε, t ) within the acceleration region ∂f∂t + ∂∂ε (cid:18) ∂ε∂t f (cid:19) = 0 , (1)with acceleration ∂ε/∂t = αε , where ε = m e c ( γ − /kT is the normalized kinetic energy and α is the constantacceleration rate in the first-order Fermi process. Weassume that the initial distribution within the layer f is Maxwellian with initial temperature kT < m e c , suchthat f ∝ γ ( γ − / exp( − ε ) (2) ≈ √ ε (cid:18) kT m e c ε + .... (cid:19) exp( − ε ) . Fermi acc. from Eq. 8energy gain from E ⊥ (c)(a)
420 460 540 x/d i γ - γ -1 Δ γ ω pe t
600 800 (d) α ( ω p e ) total energy gaincurvature drift acc.E // acc. curvature drift acc. (b) γ m a x ω pe t
200 400 600 800 1000 rec
Fig. 6.— (a) Energy as a function of x -position of four accelerated particles; (b) The maximum energy of particles in the system as afunction of time from the 1-particle count level to the 1000-particle count level. The red dashed line shows the maximum energy estimatedfor a particle moving in the direction of the reconnecting electric field at the speed of light m e c γ max = (cid:82) | qE rec | cdt ; (c) Averaged energygain and contributions from parallel electric fields and curvature drift acceleration over a time interval of 25 ω − pe as a function of particleenergy starting at ω pe t = 350;(d) α = < ∆ ε > / ( ε ∆ t ) from energy gain in perpendicular electric field and by curvature drift acceleration,and from the Equation (6) using the averaged flow speed and island size.Run σ system size λ p γ max E kin % ( J · E ) ⊥ % ατ inj d i × d i d i d i × d i d i d i × d i d i d i × d i d i d i × d i d i d i × d i d i d i × d i d i d i × d i × d i d i d i × d i d i d i × d i d i d i × d i d i d i × d i d i d i × d i d i d i × d i d i d i × d i d i d i × d i d i TABLE 1List of simulation runs with σ (cid:62) . The spectral index p , the maximum energy ( -particle level) at the end of thesimulation γ max , the percentage of magnetic energy that is converted into kinetic energy E kin % , the conversion ofmagnetic energy caused by perpendicular electric field ( J · E ) ⊥ and ατ inj estimated by tracking particles in the system. For simplicity, we consider the lowest order (nonrela-tivistic) term in this expansion and normalize f = N √ π √ ε exp( − ε ) by the number of particles N within theinitial layer. The distribution after time t is: f ( ε, t ) = 2 N √ π √ εe − αt/ exp( − εe − αt ) , (3) which remains a thermal distribution with a tempera-ture e αt T , consistent with that obtained by Drake et al.(2010). However, since upstream particles enter contin-uously into the acceleration region, the number of par-ticles in the acceleration zone increases with time. Weconsider a particle distribution f inj = N inj √ π √ ε exp( − ε )with number of particles N inj ∝ V in τ inj injected fromupstream, where τ inj is the time scale for particle injec-article Acceleration in Relativistic Magnetic Reconnection 11 f ( a r b it a r y un it s ) γ -1 -2 α τ i n j (c) σ i X i i X i i X i ∝ σ / (b) S p ec t r a l I nd e x p σ (a) i X i i X i i X i σ =625100400 Fig. 7.—
Energy spectra at the end of simulations for a series of2D runs with system size L x × L z = 600 d i × d i and different σ from 6 to 1600. (b) Spectra index for all 2D simulations with σ from 6 to 1600. (c) Time integrated ατ inj for cases with σ = 6-400and different system sizes. tion. To highlight the key role that time-dependent injec-tion plays in setting up the power-law, we first consider aquick heuristic derivation of the main result. To proceed,we split f inj into N groups and release the j th group intothe acceleration region at time t = j ∆ t . Since each groupwill satisfy the Equation 3 for a different initial time, af-ter we have injected the final group at t = τ inj , the total ω pe t E r ec / B E r ec / B σ (a)(b) σ = 100 3D Fig. 8.— (a) Time-dependent 2D and 3D reconnection electricfield normalized by the initial magnetic field E rec /B . (b) Nor-malized peak electric field E rec /B as a function of σ in 2D simu-lations. distribution (in the limit N → ∞ ) is f ( ε, t ) ∼ N inj √ πτ inj (cid:90) τ inj √ εe − αt/ exp ( − εe − αt ) dt (4)= N inj ατ inj [ erf( ε / ) − erf( ε / e − ατ inj / ) ε + 2 √ π e − ατ inj / exp( − εe − ατ inj ) − e − ε ε / ]In the limit of ατ (cid:29)
1, this gives the relation f ∝ /ε in the energy range 1 < ε < e ατ inj . Figure 11(b) shows(4) for different ατ inj . A power-law spectrum with p = 1emerges as ατ inj increases ατ inj >
1. Note that for aclosed system, since the averaged magnetic energy perparticle is only σm e c / γ max ∼ σ/ ∂f∂t + ∂∂ε (cid:18) ∂ε∂t f (cid:19) = f inj τ inj − fτ esc , (5)2 Guo et al. (a) σ
10 100 1000 (c) σ
10 100 1000 m a x ( Γ z ) σ
10 100 10001 σ
10 100 1000 (b)(d) m a x ( V x ) m a x ( V z ) m a x ( Γ x ) Fig. 9.— (a) The maximum flow velocity in the x direction V x as a function of σ ; (b) The maximum flow Lorentz factor in the x directionΓ x = 1 / (1 − V x /c ) as a function of σ ; (c) The maximum flow velocity in the z direction V z as a function of σ ; (d) The maximum flowLorentz factor in the z direction Γ z = 1 / (1 − V z /c ) as a function of σ . where τ esc is the escape time for particles. For the initialcurrent-layer distribution f and injected particle distri-bution f inj considered above, the solution can be writtenas f ( ε, t ) = 2 N √ π √ εe − (3 / β ) αt exp( − εe − αt ) (6)+ 2 N inj √ π ( ατ inj ) ε β (cid:2) Γ (3 / β ) ( εe − αt ) − Γ (3 / β ) ( ε ) (cid:3) , where β = 1 / ( ατ esc ) and Γ s ( x ) is the incomplete Gammafunction. The first term accounts for particles initially inthe acceleration region while the second term describesthe evolution of injected particles. In the limit of no in-jection or escape ( τ esc → ∞ and τ inj → ∞ ), the firstterm in (6) remains a thermal distribution the same as(3). However, as reconnection proceeds new particles en-ter continuously into the acceleration region and due tothe periodic boundary conditions there is no particle es-cape. Thus considering the case τ esc → ∞ and assuming N (cid:28) N inj , at the time t = τ inj when reconnection satu-rates the second term in (6) simplifies to (4). Thus in thelimit N ∼ N inj the first term in (6) should be retainedand the power-law produced is sub-thermal relative tothis population. While it is straightforward to obtain therelativistic corrections arising from the injected distribu-tion (2), we emphasize that these terms do not alter thespectral index. This solution explains results from oursimulations, and also appears to explain the results fromseveral recent papers, which obtained power-law distri- butions by subtracting the initial hot plasma componentin the current layer (Sironi & Spitkovsky 2014; Melzaniet al. 2014b; Werner et al. 2014). In particular, Melzaniet al. (2014b) explicitly discussed the evolution of parti-cle distribution initially in the current layer and reportedit as a heated Maxwellian distribution.In order to estimate the acceleration rate α , the en-ergy change of each particle can be approximated by arelativistic collision formula (e.g., Longair 1994)∆ ε = (cid:18) Γ V (1 + 2 V v x c + V c ) − (cid:19) ε, (7)where V is the outflow speed, Γ V = 1 / (1 − V /c ), and v x is the particle velocity in the x direction. The time be-tween two collisions is about L is /v x , where L is is the typ-ical size of the magnetic islands (or flux ropes in 3D). As-suming that relativistic particles have a nearly isotropicdistribution v x ∼ c/
2, then α = ∆ εε ∆ t ∼ c (Γ V (1 + Vc + V c ) − L is . (8)Using this expression, we measure the averaged V and L is from the simulations and estimate the time-dependent acceleration rate α ( t ). An example is shownin Figure 6 (d). This agrees reasonably well with thatobtained from perpendicular acceleration and curvaturedrift acceleration. Figure 7(c) shows the time-integratedvalue of ατ inj = (cid:82) τ inj α ( t ) dt for various simulations witharticle Acceleration in Relativistic Magnetic Reconnection 13 k ? d e / k ? E B ( k ? ) Evidence for Turbulence in Magnetic Fluctuation Spectrum y z x ! pe t = 708 J/J k ? d i Fig. 10.—
The evidence for turbulence in the 3D simulation. Left: power spectrum of magnetic fluctuations with wave numbersperpendicular to the y direction. Right: volume rendering of the current density J/J in the 3D simulation at ω pe t = 708. σ = 6 − ατ inj >
1, a hard power-law distribution with spectral index p ∼ σ and larger system size, the magnitude of ατ inj increases approximately as ∝ σ / .Better agreement between the simple model and thePIC simulations can be reached by considering the time-dependent acceleration rate α ( t ). As the magnetic recon-nection saturates, the acceleration rate decreases. Figure11 (c) shows the solution that uses the time-dependentacceleration rate α ( t ) in Figure 6(d) using a stochastic in-tegration technique described by Guo et al. (2010). Thefinal spectral index is about p = 1 .
25, similar to thatfrom the PIC simulation shown in Figure 7(a). IMPLICATIONSWe discuss the implication of the above conclusions forunderstanding the role of magnetic reconnection in mag-netically dominated astrophysical systems. Based on thecurrent understanding of magnetic reconnection, mul-tiple X-line reconnection develops when the secondarytearing instability is active in large-scale collisionlessplasma system. This process may also be importantwhen a hierarchy of collisional plasmoids (Loureiro et al. 2007; Bhattacharjee et al. 2009; Uzdensky et al. 2010)develops kinetic scale current layers that may triggercollisionless reconnection (Daughton et al. 2009; Ji &Daughton 2011). Therefore the collisionless reconnec-tion process discussed here is relevant to a range of high-energy astrophysical problems below (see Ji & Daughton2011, for a comprehensive summary of astrophysicalproblems with relevant physics).5.1.
Pulsar Wind Nebulae
In PWN models, magnetic reconnection has beenproposed as a mechanism for dissipating magnetic en-ergy in Poynting-flux dominated flows (Coroniti 1990;Lyubarsky & Kirk 2001; Kirk & Skjæraasen 2003; Porthet al. 2013) and accelerating particles to high energies(Kirk 2004). In PWNe, the emission flux usually hasspectral indices α ν = 0 − . dN/dγ ∝ γ − p with p = 1 − . p = 2 α ν + 1), too hard to be ex-plained by diffusive shock acceleration (Atoyan 1999).The recently detected > (a) V in ~V A E rec /B L is τ inj ~L z /(2V in ) L z (b) f ( γ ) / N i n j -2 -4 -6 -8 -2 γ -1 (c)10 -2 -4 -6 -2 γ -1 f ( γ ) / N i n j Thermal
Fig. 11.— (a) Illustration of the acceleration model for the forma-tion of power-law distributions; (b) Analytical results for different ατ inj obtained from Eq. (4). (c) The solution of Eq. (1) usingtime-dependent α ( t ) from Figure 6 (d). celeration theory (Abdo et al. 2011; Tavani et al. 2011;B¨uhler & Blandford 2014). There are two main possibil-ities for explaining the photon energies, (1) a relativisticDoppler boosting of the emitting region (Clausen-Brown& Lyutikov 2012) and/or (2) a strong particle acceler-ation in a nonideal electric field where E > B ⊥ , where B ⊥ is the magnetic field perpendicular to particle veloc-ity (Cerutti et al. 2013; Lyutikov et al. 2014).These observations suggest that relativistic magneticreconnection may occur in the Crab nebula. The powerlaw index revealed in this study is p = 1 −
2, consis-tent with the inferred spectra in the radio range (Atoyan1999) and in high-energy during the Crab γ − ray flares(Tavani et al. 2011). Explaining these observations re-quires a fast and efficient dissipation mechanism thatconverts a substantial fraction of magnetic energy intorelativistic particles (Lyutikov et al. 2014). In the Crab pulsar, magnetic reconnection is estimated to be in theplasmoid dominated regime and can dissipate a nontriv-ial fraction of the pulsar spin-down power (Uzdensky &Spitkovsky 2014). Our simulations have shown that for amagnetically dominated reconnection layer with σ (cid:29) max (cid:38)
10) associated with reconnection,which may boost the emission photon energy and helpto explain the observed Crab flares (Clausen-Brown &Lyutikov 2012). It is interesting to note that the recon-nection acceleration may also explain the pulsed γ -rayemission, although observations at higher energies is re-quired to further constrain the model (Mochol & P´etri2015). 5.2. AGN Jets
In AGN jets, a number of γ -ray sources have flat radiospectra with indices around α ν = 0, meaning the electronenergy distribution index may be close to p = 1 (Abdoet al. 2010; Hayashida et al. 2015). Several blazars haveshown extremely fast variability in TeV range on the or-der of several minutes (Aharonian et al. 2007; Albertet al. 2007). Hard power laws p ∼ σ e which ismeasured as magnetic energy power to the electron en-ergy power is very high up to the order of 100 (Zhanget al. 2014).Explaining the fast variability requires the relativis-tic beaming effect possibly arising from relativistic re-connection outflows (Giannios et al. 2009; Deng et al.2015). Our kinetic simulations have shown that theLorentz factor of the maximum outflow speed Γ x ∼ σ ∼ σ inferred from the model fittingcan significantly exceed unity σ (cid:29) Gamma-ray bursts
In gamma-ray bursts (GRBs), the traditional internalshock model of prompt emission is difficult to reconcilearticle Acceleration in Relativistic Magnetic Reconnection 15with observations (see Zhang & Yan 2011, and referencestherein). Magnetic reconnection and associated particleacceleration have been proposed as a key process in GRBmodels such as ICMART model (Zhang & Yan 2011)and reconnection-switch model (McKinney & Uzdensky2012). The efficient magnetic dissipation and particleacceleration during reconnection may be important tounderstand the emission mechanism in GRBs (Kumar1999; Spruit et al. 2001; Drenkhahn & Spruit 2002). Gru-ber et al. (2014) have shown a series of features in GRBprompt emission that are not consistent with the simplesynchrotron shock model. For example, the hard low-energy spectra, where the particle energy spectral indexis close to p = 1 assuming synchrotron radiation (Ghis-ellini et al. 2000; Preece et al. 2002) and the thermal emis-sion component predicted in the fireball-internal-shockmodel has been rarely seen in GRBs (Zhang & Pe’er2009; Abdo et al. 2009).From our simulation results and analytical model, theparticle energy spectral index is close to p = 1, con-sistent with low-energy photon spectra observed in mostGRBs (Band et al. 1993; Preece et al. 2000; Gruber et al.2014). The acceleration in reconnection layers is muchfaster than the radiation cooling and can maintain thehard spectrum. Using PIC simulation, Spitkovsky (2008)found that in the downstream region of highly relativis-tic shocks the number of particles in the nonthermal tailis ∼
1% of the entire downstream population, and theycarry 10% of the kinetic energy in the downstream re-gion. In our simulations of relativistic reconnection, thenumber of nonthermal relativistic particles is ∼
25% ofthe total number particles in the simulation and theycarry ∼
95% of kinetic energy in the system, meaningrelativistic reconnection is much more efficient in pro-ducing nonthermal relativistic particles. This efficientconversion from magnetic energy into kinetic energy ofnonthermal particles may help solve the efficiency prob-lem in GRB models (Zhang et al. 2007; Deng et al. 2015).5.4.
Nonrelativistic reconnection sites
While the primary focus of this paper is relativisticmagnetic reconnection, the physics of Fermi accelerationand the formation of power-law distribution is also appli-cable to the nonrelativistic regimes previously discussed(Drake et al. 2006, 2010, 2013). Based on our analyti-cal model, the power-law distribution forms only when ατ inj >
1. The results in this paper demonstrate thatthis condition is more easily achieved in regimes with σ (cid:29)
1, but it may also occur with σ < β = 8 πnkT /B ∼ .
001 - 0 .
01 ( σ < m i /m e , strong trapping at X-line region, and particleescape from the system need to be investigated further(Egedal & Daughton 2015, in preparation). DISCUSSION AND CONCLUSIONThe dissipation of magnetic field and particle energiza-tion in the magnetically dominated systems is of stronginterest in high energy astrophysics. In this study, we use2D and 3D fully kinetic simulations that resolve the fullrange of plasma physics to investigate the particle accel-eration and plasma dynamics during collisionless mag-netic reconnection in a pair plasma with magnetizationparameter σ varying from 0 .
25 to 1600. A force-free cur-rent layer, which does not require a hot plasma popu-lation in the current layer, is implemented as the initialcondition.We find that the evolution of the current sheet and ac-celeration of particles has two stages. In the early stage,an extended reconnection region forms and generates aparallel electric field that accelerates particles in the cur-rent layer. As time proceeds, the layer breaks into mul-tiple plasmoids (flux ropes in 3D) due to the secondarytearing instability. The motional electric field in the re-connection layer strongly accelerates energetic particlesvia a first-order relativistic Fermi process leading to theconversion of most of the free energy in the system. Alarge fraction of the magnetic energy is quickly convertedinto the kinetic energy of nonthermal relativistic parti-cles (within a few light-crossing times) and the eventualenergy spectra show a power law f ∝ ( γ − − p , withthe spectral index p decreasing with σ and system sizeand approaching p = 1. The formation of the power-law distribution can be described by a simple model thatincludes both inflow and the Fermi acceleration. Thismodel also appears to explain recent PIC simulations(Sironi & Spitkovsky 2014; Melzani et al. 2014b; Werneret al. 2014), which reported hard power-law distribu-tions after subtracting the initial hot plasma popula-tion inside the current layer. For the more realistic limitwith both particle loss and injection, the spectral index p = 1+1 / ( ατ esc ), recovering the classical Fermi solution.If the escape is caused by convection out of the recon-nection region τ esc = L x /V x , the spectral index shouldapproach p = 1 when ατ esc (cid:29) σ regime. Inpreliminary 2D simulations using open boundary condi-tions, we have confirmed this trend and will report else-where. For the nonrelativistic limit, the reconnectionneeds to be sustained over a longer time to form a powerlaw.We have also shown that in sufficiently high- σ regimesthe magnetic reconnection rate is enhanced and relativis-tic inflow and outflow structures develop. The scalingfollows the prediction based on the Lorentz contractionof plasma passing through the diffusion region. Although3D magnetic turbulence is generated as a consequence ofthe growth of the secondary tearing instability and kinkinstability, the particle acceleration, energy release andreconnection rate in the 3D simulation are comparableto the corresponding 2D simulation.Our study has demonstrated that relativistic mag-netic reconnection is a highly efficient energy-dissipationmechanism in the magnetically dominated regimes. Theplasma distribution in the reconnection layer features6 Guo et al.power-law energy spectra, which may be important inunderstanding the nonthermal emissions from objectslike pulsars, jets from black holes, and gamma-ray bursts.Both the inflow and outflow speeds approach the speedof light and have Lorentz factors of a few, which may ex-plain the fast variability and high luminosity observed inthose high-energy astrophysical systems. These findingson particle acceleration and plasma dynamics during rel-ativistic reconnection substantiate the important role ofmagnetic reconnection in high-energy astrophysical sys-tems. ACKNOWLEDGEMENTWe gratefully acknowledge useful discussions with andcomments from Andrey Beresnyak, Xuhui Chen, WeiCui, Wei Deng, Brenda Dingus, Jim Drake, Joe Gi-acalone, Dimitrios Giannios, Serguei Komissarov, PawanKumar, Xiaocan Li, Maxim Lyutikov, Rob Preece, MarcSwisdak, Alexander Tchekhovskoy, Dmitri Uzdensky,Yajie Yuan, Gary Zank, Bing Zhang, and HaochengZhang. This work is supported by the DOE throughthe LDRD program at LANL and DOE/OFES supportto LANL in collaboration with CMSO. The research ispart of the Blue Waters sustained-petascale computingproject, which is supported by the NSF (Grand No. OCI07-25070) and the state of Illinois. Additional simula-tions were performed at the National Center for Compu-tational Sciences at ORNL and with LANL institutionalcomputing.APPENDIX: NUMERICAL CONVERGENCE The accuracy of particle-in-cell (PIC) kinetic simula-tions depends on a series of numerical parameters suchas cell size, time step, and the number of macro-particlesper cell (e.g., Birdsall & Langdon 1991). The numericalconvergence of simulation results has been rarely explic-itly checked when modeling astrophysical problems usingPIC simulations, and often a small number of macro-particles are used. Here we examine the numerical con-vergence of our results on these numerical parametersusing VPIC code for different initial temperatures from kT = 0 .
01 to 0 . m e c . Our test case has σ = 25 withbox size L x × L z = 600 d i × d i and simulation time ω pe t = 3000. We find that numerical heating can becomeunacceptably large when a small number of particles percell is used. In Table 2 we list the key parameters for thetest. Although for most cases, the violation in energyconservation is small ( E err /E tot within 1%), the numeri-cal heating can significantly modify the particle distribu-tion since the initial kinetic energy is a small fraction ofthe total energy. Therefore to obtain trustworthy resultsthat are numerically converged, the violation of energyconservation should be much less than the initial kineticenergy E err /E k (cid:28)
1. Figure 13 shows several caseswith kT = 0 . m e c with grid number 2048 × C r = 0 .
7, and different numbers of par-ticles per cell from 8 to 512. Figure 14 shows several casesfor kT = 0 . m e c with grid number 4096 × C r = 0 . E err /E k (cid:28)
1, the numerical heat-ing has a negligible effect on the distribution function.
REFERENCESAbdo, A. A., et al. 2009, Science, 323, 1688—. 2010, ApJ, 716, 30—. 2011, Science, 331, 739Achterberg, A., Gallant, Y. A., Kirk, J. G., & Guthmann, A. W.2001, MNRAS, 328, 393Aharonian, F., Akhperjanian, A. G., Bazer-Bachi, A. R., Behera,B., Beilicke, M., Benbow, W., & Berge, D. 2007, ApJ, 664, L71Aharonian, F., Akhperjanian, A. G., Bazer-Bachi, A. R., Beilicke,M., Benbow, W., Berge, D., & Bernl¨ohr, K. 2006, Nature, 440,1018Albert, J., et al. 2007, ApJ, 669, 862Arons, J. 2012, Space Sci. Rev., 173, 341Atoyan, A. M. 1999, A&A, 346, L49Band, D., et al. 1993, ApJ, 413, 281Bessho, N., & Bhattacharjee, A. 2012, ApJ, 750, 129Bhattacharjee, A., Huang, Y.-M., Yang, H., & Rogers, B. 2009,Physics of Plasmas, 16, 112102Birdsall, C. K., & Langdon, A. B. 1991, Plasma Physics viaComputer SimulationBirn, J., et al. 2001, J. Geophys. Res., 106, 3715Blackman, E. G., & Field, G. B. 1994, Physical Review Letters,72, 494Blandford, R., & Eichler, D. 1987, Phys. Rep., 154, 1Bowers, K. J., Albright, B. J., Yin, L., Daughton, W.,Roytershteyn, V., Bergen, B., & Kwan, T. J. T. 2009, Journalof Physics Conference Series, 180, 012055B¨uhler, R., & Blandford, R. 2014, Reports on Progress inPhysics, 77, 066901Celotti, A., & Ghisellini, G. 2008, MNRAS, 385, 283Cerutti, B., Uzdensky, D. A., & Begelman, M. C. 2012, ApJ, 746,148Cerutti, B., Werner, G. R., Uzdensky, D. A., & Begelman, M. C.2013, ApJ, 770, 147Che, H., Drake, J. F., & Swisdak, M. 2011, Nature, 474, 184Chen, X., et al. 2014, MNRAS, 441, 2188Clausen-Brown, E., & Lyutikov, M. 2012, MNRAS, 426, 1374Coroniti, F. V. 1990, ApJ, 349, 538Dahlin, J. T., Drake, J. F., & Swisdak, M. 2014, Physics ofPlasmas, 21, 092304Daughton, W. 1999, Physics of Plasmas, 6, 1329 Daughton, W., & Karimabadi, H. 2007, Physics of Plasmas, 14,072303Daughton, W., Nakamura, T. K. M., Karimabadi, H.,Roytershteyn, V., & Loring, B. 2014, Physics of Plasmas, 21,052307Daughton, W., Roytershteyn, V., Albright, B. J., Karimabadi, H.,Yin, L., & Bowers, K. J. 2009, Physical Review Letters, 103,065004Daughton, W., Roytershteyn, V., Karimabadi, H., Yin, L.,Albright, B. J., Bergen, B., & Bowers, K. J. 2011, NaturePhysics, 7, 539Daughton, W., Scudder, J., & Karimabadi, H. 2006, Physics ofPlasmas, 13, 072101de Gouveia dal Pino, E. M., & Lazarian, A. 2005, A&A, 441, 845Deng, W., Li, H., Zhang, B., & Li, S. 2015, ApJ, submitted,arXiv 1501.07595Drake, J. F., Opher, M., Swisdak, M., & Chamoun, J. N. 2010,ApJ, 709, 963Drake, J. F., Shay, M. A., Thongthai, W., & Swisdak, M. 2005,Physical Review Letters, 94, 095001Drake, J. F., Swisdak, M., Che, H., & Shay, M. A. 2006, Nature,443, 553Drake, J. F., Swisdak, M., & Fermo, R. 2013, ApJ, 763, L5Drenkhahn, G., & Spruit, H. C. 2002, A&A, 391, 1141Drury, L. O. 2012, MNRAS, 422, 2474Fu, X. R., Lu, Q. M., & Wang, S. 2006, Physics of Plasmas, 13,012309Galsgaard, K., Titov, V. S., & Neukirch, T. 2003, ApJ, 595, 506Ghisellini, G., Celotti, A., & Lazzati, D. 2000, MNRAS, 313, L1Giannios, D., Uzdensky, D. A., & Begelman, M. C. 2009,MNRAS, 395, L29Gruber, D., et al. 2014, ApJS, 211, 12Guo, F., Jokipii, J. R., & Kota, J. 2010, ApJ, 725, 128Guo, F., Li, H., Daughton, W., & Liu, Y.-H. 2014, PhysicalReview Letters, 113, 155005Hayashida, M., et al. 2015, ArXiv 1502.04699Hoshino, M. 2012, Physical Review Letters, 108, 135003Huang, C., Lu, Q., & Wang, S. 2010, Physics of Plasmas, 17,072306Ji, H., & Daughton, W. 2011, Physics of Plasmas, 18, 111207 article Acceleration in Relativistic Magnetic Reconnection 17
Run kT /m e c Grid numbers Time step (Cr) NPC E err /E total E err /E k A-1 0 .
36 4096 × .
36 4096 × .
4% 38%A-3 0 .
36 4096 × .
56% 9%A-4 0 .
36 4096 × . . .
36 2048 × .
36 2048 × .
2% 20%A-7 0 .
36 2048 × .
3% 5%A-8 0 .
36 2048 × .
9% 30%A-9 0 .
36 2048 × .
45% 7%A-10 0 .
36 2048 × .
12% 1 . .
36 2048 × .
04% 0 . .
36 2048 × .
75% 12%A-13 0 .
36 2048 × .
19% 3%A-14 0 .
36 2048 × .
05% 0 . .
09 4096 × .
09 4096 × .
9% 100%B-3 0 .
09 4096 × .
42% 22%B-4 0 .
09 4096 × .
05% 2 . .
09 2048 × .
2% 212%B-6 0 .
09 2048 × .
9% 45%B-7 0 .
09 2048 × .
24% 12%B-8 0 .
09 2048 × .
6% 80%B-9 0 .
09 2048 × .
37% 19%B-10 0 .
09 2048 × .
1% 5%B-11 0 .
09 2048 × .
6% 30%B-12 0 .
09 2048 × .
13% 7%B-13 0 .
09 2048 × .
04% 2%C-1 0 .
01 4096 × .
5% 3000%C-2 0 .
01 4096 × .
7% 595%C-3 0 .
01 4096 × .
36% 126%C-4 0 .
01 4096 × .
09% 32%C-5 0 .
01 4096 × .
03% 10%C-6 0 .
01 2048 × .
75% 265%C-7 0 .
01 2048 × .
19% 67%C-8 0 .
01 2048 × .
08% 29%C-9 0 .
01 2048 × .
28% 102%C-10 0 .
01 2048 × .
08% 28%C-11 0 .
01 2048 × . .
01 2048 × .
11% 39%C-13 0 .
01 2048 × . .
01 2048 × . TABLE 2List of simulation runs used to test numerical convergence. All the runs are for σ = 25 and L x × L z = 600 d i × d i andwere performed over a duration ω pe t = 3000 . Note kT /m e c is the initial plasma temperature normalized by rest energy m e c . Time step is represented by the dimensionless Courant number C r = c ∆ t/ ∆ r , where ∆ r = ∆ x ∆ y ∆ z/ (∆ x ∆ y + ∆ y ∆ z + ∆ x ∆ z ) . NPC represents the number of particle pairs per cell. E err /E total represents theratio between change of total energy compare to the initial total energy. E err /E k represents the ratio between changeof total energy compare to the initial plasma kinetic energy. Kirk, J. G. 2004, Physical Review Letters, 92, 181101Kirk, J. G., & Skjæraasen, O. 2003, ApJ, 591, 366Kowal, G., de Gouveia Dal Pino, E. M., & Lazarian, A. 2011,ApJ, 735, 102—. 2012, Physical Review Letters, 108, 241102Krennrich, F., Dwek, E., & Imran, A. 2008, ApJ, 689, L93Krucker, S., & Battaglia, M. 2014, ApJ, 780, 107Krucker, S., Hudson, H. S., Glesener, L., White, S. M., Masuda,S., Wuelser, J.-P., & Lin, R. P. 2010, ApJ, 714, 1108Krucker, S., et al. 2008, A&A Rev., 16, 155Kulsrud, R. M. 1998, Physics of Plasmas, 5, 1599Kumar, P. 1999, ApJ, 523, L113Larrabee, D. A., Lovelace, R. V. E., & Romanova, M. M. 2003,ApJ, 586, 72Lazarian, A., & Opher, M. 2009, ApJ, 703, 8Litvinenko, Y. E. 1999, A&A, 349, 685Liu, W., Li, H., Yin, L., Albright, B. J., Bowers, K. J., & Liang,E. P. 2011, Physics of Plasmas, 18, 052105Liu, Y.-H., Daughton, W., Karimabadi, H., Li, H., & Peter Gary,S. 2014, Physics of Plasmas, 21, 022113Liu, Y.-H., Daughton, W., Karimabadi, H., Li, H., &Roytershteyn, V. 2013, Physical Review Letters, 110, 265004 Liu, Y.-H., Guo, F., Daughton, W., Li, H., & Hesse, M. 2015,Physical Review Letters, 114, 095002Longair, M. S. 1994, High energy astrophysics. Volume 2. Stars,the Galaxy and the interstellar medium.Loureiro, N. F., Schekochihin, A. A., & Cowley, S. C. 2007,Physics of Plasmas, 14, 100703Lyubarsky, Y., & Kirk, J. G. 2001, ApJ, 547, 437Lyubarsky, Y. E. 2005, MNRAS, 358, 113Lyutikov, M., Komissarov, S., & Sironi, L. 2014, to be submittedLyutikov, M., & Uzdensky, D. 2003, ApJ, 589, 893McKinney, J. C., & Uzdensky, D. A. 2012, MNRAS, 419, 573Melzani, M., Walder, R., Folini, D., Winisdoerffer, C., & Favre,J. M. 2014a, A&A, 570, A111—. 2014b, A&A, 570, A112Mochol, I., & P´etri, J. 2015, MNRAS, 449, L51Oka, M., Phan, T.-D., Krucker, S., Fujimoto, M., & Shinohara, I.2010, ApJ, 714, 915Porth, O., Komissarov, S. S., & Keppens, R. 2013, MNRAS, 431,L48Preece, R. D., Briggs, M. S., Giblin, T. W., Mallozzi, R. S.,Pendleton, G. N., Paciesas, W. S., & Band, D. L. 2002, ApJ,581, 1248 E e rr / E k0 ( % ) -2 E e rr / E t o t ( % ) NPC=832 10.1 -3 ω pe t (a)(b) f ( a r b it a r y un it s ) -2 γ -1 NPC=8 32 128 512 p=1.3
Fig. 12.—
Several cases with kT = 0 . m e c with grid number2048 × C r = 0 . E e rr / E k0 ( % ) -2 E e rr / E t o t ( % ) NPC=28 10 -3 ω pe t (a)(b) f ( a r b it a r y un it s ) -4 γ -1 NPC=2832 p=1.0 -2 p=2.3 p=2 Fig. 13.—
Several cases for kT = 0 . m e c with grid number4096 × C r = 0 ..