The interplay of magnetically-dominated turbulence and magnetic reconnection in producing nonthermal particles
DDraft version September 5, 2019
Typeset using L A TEX twocolumn style in AASTeX62
The interplay of magnetically-dominated turbulence and magnetic reconnection in producing nonthermal particles
Luca Comisso and Lorenzo Sironi Department of Astronomy and Columbia Astrophysics Laboratory, Columbia University, New York, NY 10027, USA
ABSTRACTMagnetized turbulence and magnetic reconnection are often invoked to explain the nonthermal emis-sion observed from a wide variety of astrophysical sources. By means of fully-kinetic 2D and 3Dparticle-in-cell simulations, we investigate the interplay between turbulence and reconnection in gen-erating nonthermal particles in magnetically-dominated (or, equivalently, “relativistic”) pair plasmas.A generic by-product of the turbulence evolution is the generation of a nonthermal particle spec-trum with a power-law energy range. The power-law slope p is harder for larger magnetizations andstronger turbulence fluctuations, and it can be as hard as p (cid:46)
2. The Larmor radius of particles atthe high-energy cutoff is comparable to the size l of the largest turbulent eddies. Plasmoid-mediatedreconnection, which self-consistently occurs in the turbulent plasma, controls the physics of particleinjection. Then, particles are further accelerated by stochastic scattering off turbulent fluctuations.The work done by parallel electric fields — naturally expected in reconnection layers — is responsi-ble for most of the initial energy increase, and is proportional to the magnetization σ of the system,while the subsequent energy gain, which dominates the overall energization of high-energy particles,is powered by the perpendicular electric fields of turbulent fluctuations. The two-stage accelerationprocess leaves an imprint in the particle pitch-angle distribution: low-energy particles are aligned withthe field, while the highest energy particles move preferentially orthogonal to it. The energy diffusioncoefficient of stochastic acceleration scales as D γ ∼ . σ ( c/l ) γ , where γ is the particle Lorentz factor.This results in fast acceleration timescales t acc ∼ (3 /σ ) l/c . Our findings have important implicationsfor understanding the generation of nonthermal particles in high-energy astrophysical sources. INTRODUCTIONGeneration of energetic particles far exceeding ther-mal energies is ubiquitous in the collisionless plasmasfound in space and astrophysical environments. Thus,it is not surprising that over the last several decades,significant efforts have been made to understand themechanisms of particle acceleration. Among such mech-anisms, plasma turbulence has been often invoked toexplain nonthermal particles in a variety of astrophys-ical systems (e.g. Melrose 1980; Petrosian 2012; Lazar-ian et al. 2012). Indeed, turbulence is ubiquitous in as-trophysics, in systems as diverse as stellar coronae andwinds (Matthaeus et al. 1999; Cranmer et al. 2007), theinterstellar medium (Armstrong et al. 1995; Lithwick &Goldreich 2001), supernova remnants (Weiler & Sramek1988; Roy et al. 2009), pulsar wind nebulae (Porth et al.
Corresponding author: Luca Comisso, Lorenzo [email protected], [email protected] a r X i v : . [ a s t r o - ph . H E ] S e p erate in synergy, and a comprehensive understandingof the particle acceleration physics in a turbulent en-vironment will require a detailed investigation of theirinterplay.Here, we want to study the physics of the generationof energetic particles in magnetically-dominated turbu-lence (Thompson & Blaes 1998; Cho 2005; Inoue et al.2011; Zrake & MacFadyen 2012; Cho & Lazarian 2014;Zrake 2014). In this case, the magnetic energy densityexceeds not only the pressure, but also the rest massenergy of the plasma, and the Alfv´en speed approachesthe speed of light. Understanding the process of parti-cle acceleration in this turbulence regime is importantto shed light on the bright nonthermal synchrotron andinverse Compton signatures that are routinely observedfrom high-energy astrophysical sources such as pulsarmagnetospheres and winds (B¨uhler & Blandford 2014),jets from active galactic nuclei (AGNs) (Begelman et al.1984), or coronae of accretion disks (Yuan & Narayan2014). In particular, there are several crucial questionsthat need to be answered: (i) how efficient is the tur-bulence acceleration process in these systems? (ii) whatis the slope of a (potential) power-law high-energy tailgenerated by turbulence? (iii) what is the maximum at-tainable particle energy? (iv) which physical mechanismgoverns the injection of particles from the thermal poolto higher energies? and (v) on what timescales particleacceleration proceeds?Given the complexity of the problem, an analytictreatment is often insufficient, and one must rely on nu-merical simulations. In this case, most of the previousworks have used test particle simulations, where turbu-lence was represented by prescribed fields (e.g. Micha(cid:32)lek& Ostrowsky 1996; Arzner et al. 2006; Fraschetti &Melia 2008; O’Sullivan et al. 2009; Teraki & Asano 2019)or it was provided by turbulent fields obtained fromMHD simulations (e.g. Ambrosiano et al. 1988; Dmitruket al. 2004; Kowal et al. 2012; Dalena et al. 2014; Lynnet al. 2014; Kimura et al. 2016; Beresnyak & Li 2016;Isliker et al. 2017; Gonz´alez et al. 2017; Kimura et al.2019). These approaches offer a useful strategy to studythe problem of particle acceleration with relatively inex-pensive computational simulations. On the other hand,they have also some limitations, e.g., the absence of backreaction to the imposed electromagnetic fields and ad-hoc particle injection prescriptions. These limitationsare overcome by recent hybrid (kinetic ions and fluidelectrons) (Servidio et al. 2012; Kunz et al. 2016; Pec-ora et al. 2018) and fully-kinetic PIC simulations (Zh-dankin et al. 2017; Comisso & Sironi 2018; Zhdankinet al. 2018, 2019a; Wong et al. 2019; N¨attil¨a 2019; Zh-dankin et al. 2019b), where the particle acceleration pro- cess can be followed self-consistently durying the turbu-lence evolution. These simulations have confirmed in aself-consistent way that in a collisionless plasma, turbu-lence can drive particles out of thermal equilibrium.In our earlier work (Comisso & Sironi 2018), we per-formed large-scale fully-kinetic simulations to show thatdecaying turbulence in magnetically-dominated plas-mas can generate a large fraction of nonthermal par-ticles with a power-law distribution that extends tovery high energies. The simulation domains were largeenough to capture both the MHD cascade at largescales and the kinetic cascade at small scales, andin this astrophysically-relevant setting we found thatthe power-law slope attains an asymptotic, system-size-independent value, while the high-energy cutoff in-creases linearly with the system size. Zhdankin et al.(2017, 2018) found that driven plasma turbulence isalso a viable astrophysical particle accelerator. Indeed,they showed that nonthermal energy distributions pro-duced by driven turbulence converge to a system-size-independent power-law slope for sufficiently large do-mains. In order to explain the formation of nonther-mal particle populations in magnetically-dominated tur-bulence, in Comisso & Sironi (2018) we analyzed self-consistent particle trajectories from one of the PIC sim-ulations, finding that most of the particles enter intothe acceleration process through an injection phase thatoccurs at reconnecting current sheets which form self-consistently in the turbulent system. However, we alsofound that this initial energy gain, mediated by recon-nection, is relatively small. At higher energies, particleswere stochastically accelerated by scattering off the tur-bulent fluctuations, thereby experiencing a biased ran-dom walk in momentum space.In this paper, we extend our previous analysis of theparticle acceleration process to a suite of large-scale PICsimulations. In particular, we analyze in a more ex-tended way the impact of magnetic reconnection on theinitial stage of particle acceleration, the properties of theparticle diffusion process in energy space due to stochas-tic scattering off turbulence fluctuations, and the sig-natures of these acceleration processes on the particledistribution. We show that elongated current sheets areprone to the rapid development of the plasmoid instabil-ity and break up into plasmoids/flux ropes separated bysecondary current sheets, which gives rise to fast recon-nection and efficient particle injection. Plasmoids/fluxropes are ubiquitous in both 2D and 3D simulations,as a consequence of the large scale separation betwenthe energy-containing eddies and the plasma skin depth.The initial energization of particles (i.e., at injection) iscontrolled by the work done by the electric field parallelto the local magnetic field, which is nonzero at recon-necting current sheets. On the other hand, after the firstenergization phase, the work done by the perpendicularelectric field takes over and eventually dominates theoverall energization for high-energy particles. Indeed,also the slope of the power-law high-energy tail is con-trolled by energization via perpendicular electric fields.We show that the particle pitch-angle distribution bearsmemory of the different energization processes, showingthat particle velocities are preferentially aligned withthe magnetic field at low energies, while they are pref-erentially oriented in the direction perpendicular to themagnetic field at high particle energies. We also deter-mine the diffusion coefficient in energy space that char-acterizes the physics of stochastic acceleration by tur-bulent fluctuations. In both 2D and 3D simulations, inthe energy interval pertaining to the nonthermal power-law tail, the energy diffusion coefficient increases linearlywith the plasma magnetization and quadratically withthe particle energy. For high plasma magnetizations,this yields a fast rate of particle energy gain, which canbe comparable or even higher then the particle energygain rate from fast magnetic reconnection.This paper is organized as follows. In Section 2, we de-scribe our computational method and simulations setup.This is followed, in Section 3, by a description of thefully developed turbulence state and the resulting parti-cle energy spectra for different plasma conditions. Thefollowing sections are mostly devoted to the analysis ofthe acceleration mechanisms and their signature on theparticle distribution function. In particular, in Section 4we investigate the role of magnetic reconnection in pro-viding an efficient particle injection mechanism. In Sec-tion 5 we study the different contributions of the parallel vs perpendicular electric field in driving the energizationof particles. The properties of pitch angle particle dis-tributions, four-velocity distribution functions, and mix-ing of the energized particles, are presented in Section6. Then, in Section 7, we study the properties of diffu-sion in energy space of the particles that are acceleratedby stochastic scattering off the turbulent fluctuations.Finally, in Section 8 we summarize our findings. NUMERICAL METHOD AND SETUPIn order to study the particle acceleration processfrom first principles, we solve the full Vlasov-Maxwellsystem of equations through the particle-in-cell (PIC)method (Birdsall & Langdon 1985), which evolves elec-tromagnetic fields via Maxwell’s equations and parti-cle trajectories via the Lorentz force. To this purpose,we employ the electromagnetic fully-relativistic PICcode TRISTAN-MP (Buneman 1993; Spitkovsky 2005), which allows us to perform large-scale two-dimensional(2D) and three-dimensional (3D) simulations of plasmaturbulence. In 2D our computational domain is a squareof size L in the xy -plane, while in 3D it is a cube ofsize L . We use periodic boundary conditions in all di-rections. For both 2D and 3D domains, all three com-ponents of particle momenta and electromagnetic fieldsare evolved in time.We initialize a uniform electron-positron plasma withtotal particle density n according to a Maxwell-J¨uttnerdistribution f ( p ) = 14 πm c θ K (1 /θ ) exp (cid:18) − γ ( p ) θ (cid:19) , (1)where γ ( p ) = (cid:112) p /mc ) is the particle Lorentz fac-tor, θ = k B T /mc is the dimensionless temperature,and K ( z ) is the modified Bessel function of the secondkind. Here, as usual, k B indicates the Boltzmann con-stant, T is the initial plasma temperature, m denotesthe particle mass, p is the particle momentum and c isthe speed of light in vacuum. In all the simulations, weset up a uniform mean magnetic field along the z direc-tion, B = B ˆ z . The initial equilibrium is perturbed bymagnetic fluctuations of the form δ B ( x ) = (cid:88) k δB ( k ) ˆ ξ ( k ) exp [ i ( k · x + φ k )] , (2)where δB ( k ) is the Fourier amplitude of the mode withwave vector k , ˆ ξ ( k ) = i k × B / | k × B | are Alfv´enicpolarization unit vectors, and φ k are random phases. Bysetting δB ( − k ) = δB ( k ) and φ − k = − φ k we ensure that δ B ( x ) is a real function. We adopt equal amplitude permode and wave vector components k j = 2 πn j /L withmode numbers in the interval n j ∈ { , . . . , N j } . We set N x = N y = 8 in 2D simulations, while N x = N y = 4and N z = 2 in 3D simulations. The choice of perturbinglower mode numbers in 3D simulations is due to thesmaller domain size affordable in 3D and the desire tomaximize the inertial range of the turbulent cascade.With these choices, the initial magnetic energy spectrumpeaks near k N = 2 πN max /L (where N max = 8 in 2D and N max = 4 in 3D). In the following, we will use l = 2 π/k N as our unit length, which we also refer to as the energy-carrying scale.The strength of the initial magnetic field fluctuationsis parameterized by the magnetization σ = δB πn w mc , (3)where δB rms0 = (cid:104) δB ( t = 0) (cid:105) / is the space-averagedroot-mean-square value of the initial magnetic field fluc-tuations and w mc = [ K (1 /θ ) /K (1 /θ )] mc is the Table 1.
Simulation parameters.Sim
L/d e σ δB rms0 /B θ N max . . . . . . . . . . . . . . . . . . . . . . . Note —We mark the reference simulations withan asterisk (*). The magnetization parame-ter σ is defined with the initial magnetic fieldfluctuations, σ = δB / πn w mc , where δB rms0 = (cid:104) δB ( t = 0) (cid:105) / . In this paper we usealso the instantaneous magnetization parameter σ = δB / πn wmc , where δB rms = (cid:104) δB (cid:105) / (and wmc is the instantaneous enthalpy per par-ticle), and the magnetization associated with themean magnetic field, σ z = B / πn w mc = σ ( B /δB rms0 ) . initial enthalpy per particle, with K n ( z ) indicating themodified Bessel function of the second kind of order n .Since we are interested in magnetically-dominated envi-ronments, we present results from simulations with dif-ferent values of σ (from 2 . σ (cid:29) v A = c (cid:112) σ / (1 + σ ) ∼ c . We find thatwith our definition of σ , our results do not depend onthe choice of the initial dimensionless temperature θ ,apart from an overall energy rescaling (see Section 3). We resolve the initial plasma skin depth d e = c/ω p with 10 cells in 2D and 3 cells in 3D (in 2D we havechecked that runs with d e = 3 or 10 cells give identicalresults, including the development of turbulent struc-tures, as can be seen in the Appendix). Note thatthe initial plasma skin depth is defined with the rela-tivistic plasma frequency ω p = (cid:112) πn e /γ th m , where γ th = w − θ is the initial mean particle Lorentz factor.In order to capture the full plasma turbulence cas-cade from macroscopic MHD scales to kinetic scales, wesolve the kinetic system of equations on large compu-tational domains. This is achieved by adopting a boxof 2460 cells in 3D simulations and 16400 cells in 2Dsimulations. For the 2D analysis, we also present resultsfrom three simulations with 32800 cells and one simu-lation with 65600 cells. In our reference 2D simulationwe employ 64 particles per cell, while 16 particles percell are adopted for our reference 3D simulation. Forthe other runs, we employ 16 particles per cell in 2Dand 4 particles per cell in 3D. We have tested that inthe magnetically-dominated regime of interest here, thediscussed results are the same when using up to 256 par-ticles per cell (see a particle spectrum comparison in theAppendix).The simulation timestep is controlled by the numeri-cal speed of light of 0.45 cells/timestep. The simulationsare run for ct/l = 12 −
15, at which point most of the tur-bulent magnetic energy has been transferred to the par-ticles. Our study is focused on magnetically-dominatedturbulence, and for this purpose we have performed sev-eral simulations at different magnetizations σ . In 2Dwe have investigated σ ∈ { . , , , , , } . In 3Dwe have explored σ ∈ { , , , } . If not otherwisespecified, the simulations start with δB rms0 /B = 1 and θ = 0 .
3. Cases with different δB rms0 /B and θ havealso been performed in 2D. For convenience, we havesummarized the physical parameters of the presentedsimulations in Table 1. Our reference 2D and 3D simu-lations are indicated with an asterisk. PLASMA TURBULENCE AND PARTICLESPECTRUMIn this section, we give an overview of the plasma tur-bulence state in 2D and 3D PIC simulations, with aparticular focus on the particle energy spectrum that de-velops self-consistently. We first present the character-istic fluid structures of the magnetized turbulence state,and the time evolution of the magnetic power spectrum.Then we show the time evolution of the particle energyspectrum and we discuss its dependence on the mainphysical parameters.3.1.
Plasma turbulence y / l x/l δ B y / l x/l -2.0-1.00.01.02.0 J z y / l x/l -1.7-1.1-0.50.10.7 l og ( Γ β ) y / l x/l -1.2-0.60.00.61.2 l og ( n / n ) Figure 1.
2D plots of different fluid structures in fully developed 2D turbulence (at ct/l = 4 .
6) with σ = 10, δB rms0 /B = 1,and L/d e = 1640 (with l = L/ B / π , the current density J z along the mean magnetic field in units of en c , the bulk dimensionlessfour-velocity Γ β , and the particle density ratio n/n . Note that the color bars for Γ β and n/n are in logarithmic scale. Turbulence structures from our reference 2D simula-tion are illustrated in Fig. 1. We plot the magneticfield squared fluctuations δB , the out-of-plane elec-tric current density J z , the bulk dimensionless four-velocity Γ β , and the particle density ratio n/n . Here,Γ = 1 / (cid:113) − ( V /c ) indicates the plasma bulk Lorentzfactor and β = | V | /c is the dimensionless plasma bulkspeed obtained by averaging the velocities of individ-ual particles. We can see that the fluctuations δB aregenerally stronger in large-scale flux tubes (see the circu-lar structures of size comparable to the energy-carryingscale l ), but high values of δB are also obtained insmall-scale structures identified with reconnection plas-moids (see the circular structures with size (cid:28) l ). Theseare “secondary” magnetic islands (flux ropes in 3D) that are produced by magnetic reconnection (Biskamp2000). In such plasmoids, the particle number density n exhibits strong enhancements in excess of n ∼ n .High values of particle number density occur in large-scale flux tubes as well. In general, the particle densitydisplays strong compressions in coherent quasi-circularstructures spanning a range of scales.In between flux tubes, reconnection layers reveal theformation of plasmoids within narrow current sheets.Indeed, current sheets with high aspect ratio tend tofragment into plasmoids and secondary current sheetsas a result of magnetic reconnection. Smaller-size cur-rent sheets are also ubiquitous, spanning a wide range ofscales. We will see in the following sections that recon-necting current sheets, which are a natural by-productof turbulent cascades in magnetized plasmas (e.g. Ser- Figure 2.
3D plots of different fluid structures in fully developed 3D turbulence (at ct/l = 2 .
7) with σ = 10, δB rms0 /B = 1,and L/d e = 820 (with l = L/ B / π , the current density J z along the mean magnetic field in units of en c , the bulk dimensionlessfour-velocity Γ β , and the particle density ratio n/n . Note that the color bars for Γ β and n/n are in logarithmic scale. Ananimation showing the current density J z in different xy slices can be found at https://doi.org/10.7916/d8-prt9-kn88. vidio et al. 2009; Wan et al. 2013; Cerri & Califano 2017;Franci et al. 2017; Haggerty et al. 2017; Dong et al. 2018;Comisso & Sironi 2018; Papini et al. 2019), play an im-portant role for particle injection into the accelerationprocess (Comisso & Sironi 2018). Finally, we also pointout that in the strongly magnetized regime of plasmaturbulence investigated here, the plasma bulk speed canreach very high values. In particular, we observe ultra-relativistic flows with bulk Lorentz factor as high asΓ ∼
5. Such high speeds develop predominantly in be- tween the large-scale flux tubes, although high-velocityfluctuations occur all over the spatial domain.We now consider the fluid structures that develop in3D plasma turbulence. Our reference 3D simulation has
L/d e = 820, which is half the size of the reference 2Dsimulation. However, since in 3D we adopt perturbationnumbers up to N max = 4 (as compared to N max = 8 in2D), we still have a well-extended turbulence inertialrange. In fact, the ratio of the initial energy-carryingscale l = 2 π/k N to the plasma skin depth d e remainsthe same between our reference 2D and 3D simulations, -3 -2 -1 kd e -5 -4 -3 -2 -1 P B ( k ) k - k - k - ct/l ( δ B r m s / B ) Figure 3.
Power spectrum of the magnetic field for the 2Dsimulation in Fig. 1, showing a well-developed inertial rangeand a kinetic range scaling roughly as P B ( k ) ∝ k − . . Theinset shows the time evolution of δB = (cid:104) δB (cid:105) normalizedto B , with vertical dashed lines indicating the times whenthe magnetic power spectra presented in the main panel arecomputed (same color coding). leading to the same high-energy cutoff of the particleenergy spectrum (see Comisso & Sironi (2018) as wellas Eq. (9) in the following subsection).The turbulence structures from our reference 3D sim-ulation are displayed in Fig. 2. The magnetic fieldsquared fluctuations δB present both large-scale andsmall-scale structures. However, in this case, the large-scale fluctuations are not organized in coherent fluxtubes (as it was in 2D, where they were a result of theconstrained 2D dynamics). Despite differences in thelarge-scale structure of the magnetic field, there is still acopious presence of current sheets (current ribbons whenconsidering the third direction). Due to the presence ofthe mean magnetic field B = B ˆ z , current ribbons aremostly elongated along ˆ z . We can see that the z com-ponent of the electric current, J z , displays a variety ofcurrent sheets of different sizes. Some of these currentlayers break into plasmoids (see also Sec. 4), as highlyelongated layers cannot be stable against the plasmoidinstability, also in 3D geometry (e.g. Daughton et al.2011; Sironi & Spitkovsky 2014; Guo et al. 2015; Huang& Bhattacharjee 2016; Ebrahimi 2017; Werner & Uzden-sky 2017; Baalrud et al. 2018; Stanier et al. 2019). Here,we show that plasmoids/flux ropes are self-consistentlycreated in fully 3D plasma turbulence (see Sec. 4), wherecurrent sheets are self-consistently generated by the tur-bulence itself. As for 2D plasma turbulence, we will seethat these current sheets play an important role in theinitial stages of particle acceleration (Sections 4-6). Locations characterized by strong electric current den-sities are typically accompanied by strong gradients inparticle density. In localized regions, the particle den-sity can exceed n ∼ n , similar to the 2D case. Onthe other hand, large-scale structures like the overdenseregions at the core of 2D large-scale flux tubes are miss-ing in 3D. Finally, we observe that also in 3D, due tothe high magnetization of the system, the plasma flowspeed is generally very high. We can see regions withultra-relativistic flow speeds having bulk Lorentz factoras high as Γ ∼ P B ( k ) for the reference2D simulation, where P B ( k ) dk = (cid:88) k ∈ dk δ B k · δ B ∗ k B (4)is computed from the discrete Fourier transform δ B k ofthe fluctuating magnetic field. Each curve refers to a dif-ferent time (from brown to orange), as indicated by thecorresponding vertical dashed lines in the inset, wherewe present the temporal decay of the energy in turbulentfluctuations δB /B . We can see that at MHD scales( kd e (cid:46) .
5) the magnetic power spectrum is consistentwith a Kolmogorov scaling P B ( k ) ∝ k − / (Biskamp2003) (compare with the dot-dashed line), while theIroshnikov-Kraichnan scaling P B ( k ) ∝ k − / (Irosh-nikov 1963; Kraichnan 1965) (triple-dot-dashed line) ispossibly approached at late times. At kinetic scales( kd e (cid:38) . P B ( k ) ∝ k − . (compare with thedashed line). A similar slope was proposed for magne-tized turbulence at sub-inertial scales in a cold plasma(Abdelhamid et al. 2016; Passot et al. 2017). We finallyobserve that the turbulence integral scale (cid:96) ( t ) = 2 πk I ( t ) = 2 π (cid:82) ∞ k − P B ( k, t ) dk (cid:82) ∞ P B ( k, t ) dk , (5)which is close to the energy-carrying scale associatedto the wavenumber where P B ( k, t ) peaks, increases asthe magnetic energy decays in time. This is due to themerging of the large-scale magnetic flux tubes, whichdrives an inverse energy transfer to scales larger thanthe initial integral scale (e.g. Biskamp & Schwarz 2001). -2 -1 k ⊥ d e -5 -4 -3 -2 -1 P B ( k ⊥ ) k ⊥ - k ⊥ - k ⊥ - ct/l ( δ B r m s / B ) Figure 4.
Power spectrum of the magnetic field for the 3Dsimulation in Fig. 2, showing a well-developed inertial rangeand a kinetic range scaling roughly as P B ( k ) ∝ k − . . Theinset shows the time evolution of δB = (cid:104) δB (cid:105) normalizedto B , with vertical dashed lines indicating the times whenthe magnetic power spectra presented in the main panel arecomputed (same color coding). We now consider the time evolution of the magneticpower spectrum in 3D. Due to the presence of the large-scale mean magnetic field B , turbulence becomes in-creasingly anisotropic toward small scales, within the in-ertial range. To account for this global anisotropy withrespect to B , we consider the magnetic power spectrumwith respect to the wavenumber k ⊥ = ( k x + k y ) / per-pendicular to the mean field, obtained from the discreteFourier transform of the fluctuating magnetic field as P B ( k ⊥ ) dk ⊥ = (cid:88) k ∈ dk ⊥ δ B k · δ B ∗ k B . (6)Figure 4 shows the time evolution of the magnetic powerspectrum P B ( k ⊥ ), which does exhibit inertial and ki-netic ranges of the turbulence cascade at times ct/l (cid:38) k ⊥ d e (cid:46) .
5) of the magnetic power spectrumtends to flatten from P B ( k ⊥ ) ∝ k − / ⊥ (Goldreich & Srid-har 1995; Thompson & Blaes 1998) (dot-dashed line)to P B ( k ⊥ ) ∝ k − / ⊥ (Boldyrev 2006) (triple-dot-dashedline). At kinetic scales ( k ⊥ d e (cid:38) . P B ( k ⊥ ) ∝ k − . ⊥ (dashed line),similar to the 2D result and in agreement with theoreti-cal predictions for magnetized turbulence at sub-inertialscales in cold plasmas (Abdelhamid et al. 2016; Passotet al. 2017). Note also that in the 3D case, the magneticenergy decays faster than in the 2D case (compare in-sets of Figs. 3 and 4). We will show that this leads to areduced particle acceleration rate at late times. -1 γ - -6 -5 -4 -3 -2 -1 d N / d l n ( γ - ) p=2.8 ct/l σ p ( δ B rms0 /B ) = δ B rms0 /B ) = Figure 5.
Time evolution of the particle spectrum dN/d ln( γ −
1) for the simulation in Fig. 1. At late times,the particle spectrum displays a power-law tail with index p = − d log N/d log( γ − ∼ .
8. About 17% of the particleshave γ ≥
12 at ct/l = 12 (twice the peak of the particleenergy spectrum at that time), which gives an indication ofthe percentage of nonthermal particles. The inset shows thepower-law index p as a function of the magnetization σ fortwo values of δB rms0 /B . Particle spectrum
The most interesting outcome of the turbulent cascadeis the generation of a large population of nonthermalparticles. This is shown in Fig. 5 (for the 2D setup),where the time evolution of the particle energy spec-trum dN/d ln( γ −
1) is presented ( γ − E k /mc is thenormalized particle kinetic energy). As a result of tur-bulent field dissipation, the spectrum shifts to energiesmuch larger than the initial Maxwellian, which is shownby the blue line peaking at γ − ∼ γ th − (cid:39) .
6. At latetimes, when most of the turbulent energy has decayed,the spectrum stops evolving (orange and red lines): itpeaks at γ − ∼
5, and extends well beyond the peakinto a nonthermal tail of ultra-relativistic particles thatcan be described by a power-law dNdγ = N (cid:18) γ − γ st − (cid:19) − p , for γ st < γ < γ c , (7)and a sharp cutoff for γ ≥ γ c . Here, N is the normal-ization of the power-law and p is the power-law index,which is about 2 . dN/d ln( γ −
1) to emphasize the particle content, propor-tional to ( γ − − p +1 dγ for the distribution in Eq. (7)).The percentage of the particles in the nonthermal tail(measured as the number of particles with Lorentz factorexceeding twice the thermal peak) is high, ζ nt ∼ N , which is close to the thermal peak. The startingpoint of the power-law, γ st , is roughly only a factor oftwo larger than the peak of the particle energy spec-trum at late times. Therefore, dropping O (1) factors,the starting point of the power-law can be estimated as γ st ∼ γ σ = (cid:16) σ (cid:17) γ th , (8)since most of the magnetic energy is converted to par-ticle energy by the time the particle energy spectrumhas saturated (see inset of Fig. 3). On the other hand,the high-energy cutoff γ c depends on the system size.As discussed in the following Sections, stochastic accel-eration by turbulent fluctuations dominates the energygain of the most energetic particles. High-energy par-ticles cease to be efficiently scattered by turbulent fluc-tuations when their Larmor radius ρ L = ( γmc/eB ) v ⊥ exceeds the integral length scale (cid:96) = 2 π/k I , implying anupper limit to their Lorentz factor of γ c ∼ e (cid:112) (cid:104) B (cid:105) (cid:96)mc ∼ πk I d e √ σ z γ th , (9)where (cid:104) B (cid:105) is the space-averaged mean-square value ofthe magnetic field, and σ z = B πn w mc = σ (cid:18) B δB rms0 (cid:19) , (10)This argument assumes that the turbulence surviveslong enough to allow the particles to reach this upperlimit (we also assumed B /δB rms (cid:38) k I ∼ k N , was presentedin Comisso & Sironi (2018) by performing simulationswith different domain sizes. We point out that inversemagnetic energy transfer (e.g. Biskamp & Schwarz 2001;Zrake 2014; Brandenburg et al. 2015) can possibly drivea substantial decrease in time of k I , which, in turn, canallow the most energetic particles to reach even higherenergies.We observe that the slope of the power-law is not uni-versal, but it depends on the magnetization σ and theratio δB rms0 /B (Comisso & Sironi 2018). The inset ofFig. 5 shows how the power-law index changes with themagnetization σ from two series of simulations having δB rms0 /B = 1 and δB rms0 /B = 2. We can see that theslope of the power-law becomes harder for larger magne-tization, and that for fixed σ it is harder when increas-ing δB rms0 /B (see also Fig. 7). The decrease of thepower-law index p for increasing magnetization σ (seealso Zhdankin et al. (2017); Comisso & Sironi (2018))is in analogy with the results of PIC simulations ofrelativistic magnetic reconnection (Sironi & Spitkovsky2014; Guo et al. 2014; Werner et al. 2016; Lyutikov et al. -1 γ - -6 -5 -4 -3 -2 -1 d N / d l n ( γ - ) p=2.9 ct/l p σ γ c Figure 6.
Time evolution of the particle spectrum dN/d ln( γ −
1) for the simulation in Fig. 2. At late times,the spectrum displays a power-law tail with index p = − d log N/d log( γ − ∼ .
9. About 16% of the particles have γ ≥
15 at ct/l = 12 (twice the peak of the particle energyspectrum), which gives an indication of the percentage ofnonthermal particles. The inset shows the power-law index p and the cutoff Lorentz factor γ c as a function of the magne-tization σ . The dashed line indicates the scaling γ c ∝ σ / expected for a σ -independent domain size L/d e = 820. dN/d ln( γ −
1) starting from the initial Maxwellianpeaked at γ − ∼ γ th − (cid:39) .
6. As time progresses, theparticle energy spectrum shifts to higher energies anddevelops a high-energy tail containing a large fractionof particles. At late times, when most of the turbulentenergy has decayed, the particle energy spectrum stopsevolving (orange and red lines) and it peaks at γ − ∼ p ∼ . ct/l = 12 we find thatabout 16% of particles have or exceed twice the energy0 γ - -4 -3 -2 -1 d N / d l n ( γ - ) p=1.9 δ B rms0 /B = δ B rms0 /B = δ B rms0 /B = Figure 7.
Particle spectra dN/d ln( γ −
1) at late timesfor simulations with magnetization σ = 40, system size L/d e = 3280 (with l = L/ δB rms0 /B ∈ { , , } . For the case with largerinitial fluctuations, the late time particle spectrum displays apower-law tail with index p = − d log N/d log( γ − ∼ . γ ≥
25 at ct/l = 12 (twicethe peak of the particle energy spectrum at that time), whichwhich gives an indication of the percentage of particles in thenonthermal tail. of the spectral peak, which provides an indication of thepercentage of particles in the nonthermal tail ζ nt .In order to understand the dependence of the high-energy power-law slope on the initial magnetization in3D, we performed four large-scale 3D simulations with σ ∈ { , , , } and same δB rms0 /B = 1, L/d e =820. The power-law index p decreases for increasing σ (see top inset in Fig. 6), with values that are close tothe ones from the corresponding 2D simulations with δB rms0 /B = 1 (blue curve from the inset in Fig. 5).Here we show also the scaling of the high-energy cutoff γ c (bottom inset in Fig. 6), defined as the Lorentz fac-tor where the spectrum drops one order of magnitudebelow the power-law best fit. The high-energy cutoff γ c increases as γ c ∝ σ / (compare with dashed linein the inset), which is consistent with the expectationfrom Eqs. (9) and (10) for a σ -independent domainsize L/d e and fixed δB rms0 /B .Several astrophysical systems are thought to have δB rms /B larger than unity (e.g., δB /B ∼ δB rms0 /B = 1 , ,
4, withfixed initial magnetization σ = 40, and a larger domainsize L/d e = 3280. Fig. 7 shows that the power-law be-comes harder with increasing δB rms0 /B , with p < ( γ - / γ th -6 -5 -4 -3 -2 -1 d N / d l n ( γ - ) p=2.8 θ = θ = θ = θ = θ = Figure 8.
Particle spectra dN/d ln( γ −
1) at ct/l = 12for simulations with fixed σ = 10, δB rms0 /B = 1, and L/d e = 1640 (with l = L/ θ = k B T /mc ∈ { . , . , , , } . The x -axishas been normalized to the initial thermal Lorentz factor γ th to facilitate comparison among the different cases. subject to energy constraints, as we now discuss. Thestarting point of the power-law tail, γ st , could be lowerthan indicated in Eq. (8), if only a minor fraction ofthe available energy goes into thermal particles, whilemost of the energy goes into the nonthermal tail. Also,while in the case p > γ c → ∞ as k I d e →
0, the case 1 < p < γ c , since the mean energyper particle in the power-law tail has to be (Sironi &Spitkovsky 2014)1 − p − p ( γ c − − p − ( γ st − − p ( γ c − − p − ( γ st − − p = (cid:16) χ σ (cid:17) γ th , (11)where χ is the fraction of turbulent magnetic energyconverted into particles belonging to the power-law tail.We conclude this section with the results of 2D sim-ulations having different initial plasma temperature θ .From Fig. 8, we can see that the slope p , the fraction ofparticles in the nonthermal tail, and the extent of thenonthermal tail γ c /γ st do not depend on θ . Indeed, thisplot shows that spectra obtained from simulations withdifferent θ nearly overlap, when shifted by an amountequal to the initial thermal Lorentz factor γ th . The roleof the initial choice of temperature is only to produce anenergy rescaling, since both γ st and γ c are proportionalto γ th , as can be seen from the relations (8) and (9),and the definitions of σ and √ σ z /d e already take intoaccount relativistic thermal effects.Up to this point, we have discussed general features ofthe particle spectrum generated as a by-product of theplasma turbulence. We have found that despite some1differences between 2D and 3D settings, the producedparticle spectrum does not depend on the dimensional-ity of the simulation domain (see also Comisso & Sironi(2018)). In both cases, the high-energy power-law rangeextends from (about) the thermal peak to a maximumenergy set by the energy-containing scale of turbulence.These common features, combined with the fact thatthe slope of the power-law is also similar, yield a simi-lar percentage of particles in the power-law tail. In thenext Sections, we will shed light on the particle acceler-ation mechanisms that produce the nonthermal particlespectrum. PARTICLE INJECTION AND FASTRECONNECTIONIn this section, we investigate the physics behind theinitial rapid acceleration of particles from low energies( γmc ∼ γ th mc ), to energies well above the thermalpeak ( γmc (cid:29) γ th mc ), which is usually referred to asthe injection mechanism. The investigation of the in-jection mechanism will not be limited to this section,but it will be pursued also in parts of Sections 5 and6. Here, specifically, as a continuation of our earlieranalysis (Comisso & Sironi 2018), we want to examinethe spatial locations where the injection process occurs,and understand what is special about these locations.To this aim, we have tracked the time evolution of alarge sub-sample of particles that were randomly se-lected from our reference PIC simulations. Followingin time their trajectory and energy evolution, we cananalyze, for the fraction of particles that experience aninjection process, the physical conditions at the momentof their rapid initial acceleration phase. Then we calcu-late the conditions for having efficient particle injectionby reconnection, which are linked to the onset of fastmagnetic reconnection mediated by the plasmoid insta-bility. Indeed, despite their small filling fraction, weshow that reconnecting current sheets can inject a largefraction of particles in a few outer-scale eddy turnovertimes.4.1. Particle injection at reconnecting current sheets
We begin our analysis from the reference 2D case, andthen we extend the analysis to the reference 3D case.For the injection analysis presented in this section, weemployed a sub-sample of ∼ tracked particles for the2D case, and a sub-sample of ∼ tracked particles forthe 3D case.We show in Fig. 9(a) the time evolution of the Lorentzfactor for 10 representative particles that eventuallypopulate the nonthermal tail at ct/l = 12 (see parti-cle spectrum in Fig. 5). These particles have a distinct ct/l γ (a) 0 5 10 15 20 25 30 | J z,p | /J z, rms -4 -3 -2 -1 P D F high-energy particles at t inj all particles at ct/l=3.5 (b) Figure 9.
Relation between particle injection and elec-tric current density from the 2D simulation with σ = 10, δB rms0 /B = 1, and L/d e = 1640. Top frame: Time evo-lution of the Lorentz factor for 10 representative particlesselected to end up in different energy bins at ct/l = 12(matching the different colors in the color bar on the right).Bottom frame: Probability density functions of | J z,p | /J z, rms experienced by high-energy particles at their injection time t inj (red circles) and by all our tracked particles at ct/l = 3 . | J z,p | ≥ J z, rms . moment in which they are “extracted” from the thermalpool at γ ∼ γ th and injected to higher Lorentz factors γ (cid:29) γ th . To identify this moment, that we call injec-tion time t inj , we evaluate when the rate of increase ofthe particle Lorentz factor (averaged over c ∆ t/d e = 45)satisfies ∆ γ/ ∆ t ≥ ˙ γ thr , and prior to this time the par-ticle Lorentz factor was γ ≤ γ th ∼
6. We take thethreshold ˙ γ thr (cid:39) . √ σ γ th ω p , but we have verifiedthat our identification of t inj is nearly the same whenvarying ˙ γ thr around this value by up to a factor of three(the factor 0 .
01 is much lower than the typical colli-sionless reconnection rate [ ∼ .
1, in units of the Alfv´enspeed], which is the appropriate reference scaling here,as showed in Comisso & Sironi (2018) and below).2 y / l x/l | J z | ≥ J z, rms ct/l = ct/l = ct/l = y / l x/l x/l x/l Figure 10.
Spatial correlation between particle injectionand reconnecting current sheets for the same simulation asin Fig. 9. Top frame: Regions of space with | J z | ≥ (cid:10) J z (cid:11) / (shown in black) at ct/l = 3 .
5, with red circles indicating thepositions of the particles undergoing injection around thistime. Bottom frames: Shaded isocontours of J z in the spa-tial domain ( x/l, y/l ) ∈ [2 . , . × [2 . , .
9] (correspondingto the area within the rectangular blue contour in the topframe) at times ct/l = 3 . ct/l = 3 . ct/l = 3 . J z <
0, while red indicates regions with J z > Once t inj is determined for the population of parti-cles at hand, it is possible to explore the properties ofthe electromagnetic fields at the injection location. Inthis case, by analyzing the fields at injection, we findthat the out-of-plane current density J z is particularlyrevealing. In particular, J z has, in general, high valuesat injection locations. To provide a statistical measureof the likelihood of this occurrence, we can construct the probability density function (PDF) of the magnitude ofthe out-of-plane electric current density experienced bythe particles at their injection time, | J z,p | , normalized by J z, rms , i.e., the standard deviation of the current den-sity J z, rms = (cid:10) J z (cid:11) / in the whole domain at that time.The outcome of this analysis is shown in Fig. 9(b) bythe red circles. The PDF of the high-energy particlesat injection should be contrasted with the PDF of theentire population of particles at a representative time(here, ct/l = 3 . | J z,p | ∼ J z, rms . In par-ticular, approximately ∼
95% of the high-energy parti-cles are injected at locations with | J z,p | ≥ J z, rms . Onthe other hand, by taking all the particles at the rep-resentative time ct/l = 3 .
5, only ∼
9% of them happento be in regions where | J z,p | ≥ J z, rms . Note also thatthe PDF of the overall particle population does not fol-low Gaussian statistics due to the intermittent natureof current sheets in turbulence (e.g. Servidio et al. 2009;Cerri et al. 2017; Haggerty et al. 2017; Dong et al. 2018),which is therefore reflected in the PDF of the particlesthat sample the entire domain.To obtain further insight, we look now at the mor-phology of regions with out-of-plane current density | J z | ≥ J z, rms , and we correlate it with the spatial loca-tions of the particles undergoing injection at t inj . This isshown in Fig. 10(a), where we can see that the vast ma-jority of the structures with | J z | ≥ J z, rms are sheet-likestructures, namely current sheets, and the overwhelm-ing majority of particles at injection resides in these re-gions. A large fraction of these current sheets are ac-tive reconnection layers, fragmenting into plasmoids. Atypical case of such reconnecting current sheets is il-lustrated in Fig. 10(b), where we show a small por-tion of the domain, corresponding to the area withinthe rectangular blue contour in Fig. 10(a), at differenttimes ct/l = 3 . , . , .
5. The reconnecting current sheetevolves in time and breaks up in shorter sheets due tothe formation of plasmoids. During this period of time,particles are constantly injected up to nonthermal ener-gies, as shown by the red circles in Fig. 10(b).These results are also robust in 3D, for which we haveperformed the same type of analysis. Fig. 11(a) showsthe time evolution of the Lorentz factor for 10 repre-sentative particles selected to end up in different energybins of the nonthermal tail at ct/l = 12 (see particlespectrum in Fig. 6). As in 2D, we can see a sudden ac-celeration episode with particles that are extracted from3 ct/l γ (a) | J z,p | /J z, rms -4 -3 -2 -1 P D F high-energy particles at t inj all particles at ct/l=2.5 (b) Figure 11.
Relation between particle injection and elec-tric current density from the 3D simulation with σ = 10, δB rms0 /B = 1, and L/d e = 820. Top frame: Time evo-lution of the Lorentz factor for 10 representative particlesselected to end up in different energy bins at ct/l = 12(matching the different colors in the color bar on the right).Bottom frame: Probability density functions of | J z,p | /J z, rms experienced by the high-energy particles at their t inj (redcircles) and by all our tracked particles at ct/l = 2 . | J z,p | ≥ J z, rms . the thermal pool and injected into the acceleration pro-cess. We identify the injection time t inj as for the 2Dcase, by evaluating when the Lorentz factor increases ata rate exceeding the same threshold ˙ γ thr adopted for 2D,starting from a value that is γ ≤ γ th ∼ | J z,p | /J z, rms for the high-energy particles at injection and for all particles at arepresentative time (taken at ct/l = 2 . y / l x/l | J z | ≥ J z, rms y / l x/l | J z | ≥ J z, rms Figure 12.
Spatial correlation between particle injectionand reconnecting current sheets for the same 3D simulationas in Fig. 11. In black, we show regions of space with strongcurrent density | J z | ≥ (cid:10) J z (cid:11) / at ct/l = 2 .
5, for two repre-sentative planes of the 3D domain, taken at z/l = 0 . z/l = 3 . B is in the out-of-plane direction. The redcircles indicate the positions of particles undergoing injectionaround this time. case. Fig. 11(b) indeed shows that the PDF of the parti-cles at injection (red circles) peaks at | J z,p | /J z, rms ∼ . ct/l = 2 . | J z,p | /J z, rms ∼
0. Again, parti-cles at injection feel a substantial electric current den-4
Figure 13.
Chain of flux ropes formed in a reconnectingcurrent sheet that self-consistently develops in 3D turbulence(with σ = 10, δB rms0 /B = 1, and L/d e = 820). Isosur-faces of the current density J z are shown in blue color in thezoomed region, highlighting four flux ropes (3D plasmoids)elongated along ˆ z , i.e., the direction of the mean magneticfield. The color scheme for the shaded isocontours is suchthat blue indicates regions with J z <
0, while red indicatesregions with J z > sity in the direction of the mean magnetic field. Thepeak of the PDF for the particles at injection is at alower value of | J z,p | /J z, rms than in 2D, and in generalthere are weaker | J z,p | /J z, rms wings for both the PDFof all particles and the PDF of particles experiencinginjection. This can be attributed to the lower levels ofintermittency that characterize 3D magnetized turbu-lence with respect to its 2D counterpart (e.g. Biskamp2003). Nevertheless, about 80% of the particles are in-jected in regions with | J z,p | ≥ J z, rms . On the otherhand, only approximately 11% of the entire populationof particles (at the representative time ct/l = 2 .
5) re-side at | J z,p | ≥ J z, rms . Therefore, also in 3D, special locations of high electric current density are associatedwith particle injection.The spatial locations with | J z | ≥ J z, rms are associ-ated with current ribbons that are predominantly elon-gated along the mean magnetic field B . In Fig. 12, weshow the morphology of these regions for two represen-tative planes perpendicular to B (taken at ct/l = 2 . Plasmoid-mediated disruption of the current sheetsand efficiency of reconnection-mediated injection
Reconnecting current sheets are a viable source of par-ticle injection in typical astrophysical systems ( (cid:96) ≫ d e )only if the injection efficiency (i.e., the fraction of par-ticles going through the injection phase) is large andindependent of system size. Here we show that this isindeed expected for our turbulence studies.The rate at which a reconnecting current sheet canprocess particles is proportional to the normalized re-connection speed β R = v R /c , which essentially quan-tifies the speed of the reconnection process. This ratewould be low for very elongated current sheets, as thelarge aspect ratio has the effect of throttling the recon-nection rate. Indeed, a stable current sheet would beable to reach an asymptotic width determined by the mi-crophysics of the plasma. For a relativistic pair plasma,the steady-state solution for the half-width of a recon-necting current sheet is (Comisso & Asenjo 2014) λ ∞ (cid:39) d w = (cid:114) mc πne w , (12)where w = K (1 /θ ) /K (1 /θ ) is the enthalpy per particlein units of mc . For a thinning current sheet, λ ∞ is theasymptotic limit of its half-width. For θ = k B T /mc (cid:29) d w = (cid:112) γ th mc / πne ∼ d e . Then, for a currentsheet of half-length ξ (cid:29) λ ∞ , and a compression ratiobetween inflow and outflow of order unity, the steady-5 y / l x/l Figure 14.
Chains of plasmoids in plasma turbulence from a 2D simulation with
L/d e = 6560 ( σ = 10, δB rms0 /B = 1).The shaded isocontours represent the electric current density J z in a portion of the spatial domain given by ( x/l, y/l ) ∈ [2 . , . × [1 . , .
0] at time ct/l = 4 .
5. The color scheme is such that blue represents the most negative value, and red the mostpositive value. Zoomed-in subdomains are used to reveal one plasmoid chain. ct/l = ct/l = ct/l = y / l x/l x/l x/l Figure 15.
Plasmoid formation and development from a2D simulation with
L/d e = 6560 ( σ = 10, δB rms0 /B = 1).The shaded isocontours represent the electric current density J z in a portion of the spatial domain given by ( x/l, y/l ) ∈ [7 . , . × [2 . , .
2] at times ct/l = 4 . ct/l = 4 . ct/l = 4 . J z <
0) to red ( J z > state reconnection rate is β R ∼ d w ξ (cid:28) . (13)Since current sheets generated by outer-scale eddies(which, as we discuss below, are the ones that dominateparticle injection) have half-length ξ ∼ (cid:96) larger than d w by many orders of magnitude, the reconnection rate, aswell as the injection efficiency, would be extremely lowin this scenario.However, plasmoids (which form copiously in our sim-ulations) can break the reconnection layer into shorterelements, consequently leading to a regime of fast non-linear reconnection (Daughton et al. 2006; Daughton &Karimabadi 2007; Daughton et al. 2009; Bhattachar-jee et al. 2009; Huang & Bhattacharjee 2010; Uzden-sky et al. 2010). This can happen if the plasmoids dis-rupt the current sheet within its characteristic lifetime,i.e., within one nonlinear eddy turnover time (Carboneet al. 1990; Mallet et al. 2017; Loureiro & Boldyrev2017; Boldyrev & Loureiro 2017; Comisso et al. 2018;Dong et al. 2018; Walker et al. 2018). Fast magneticreconnection essentially begins when plasmoids becomenonlinear, namely when the current density fluctuationscaused by the growing plasmoids are of the same orderof the current density of the reconnection layer (see Fig.4 in Huang et al. (2017)). Therefore, understanding theplasmoid formation in the context of a forming currentsheet is essential to understand the onset of fast mag-netic reconnection and ensuing particle injection.In order to evaluate the conditions for plasmoid for-mation and current sheet disruption, we need to analyzethe growth rate of tearing (or “reconnecting”) modes insuch current sheet. The tearing mode dispersion rela-tion for a relativistic pair plasma can be obtained fromthe relativistic pair-plasma fluid equations (e.g. Koide2009) by applying the standard tearing mode analysis(Furth et al. 1963; Coppi et al. 1976; Ara et al. 1978).6In this way, one can obtain, for arbitrary values of thetearing stability parameter ∆ (cid:48) (Furth et al. 1963), thedispersion relation γ / τ / H (cid:18) λd w (cid:19) / Γ [(Υ − / /
4] = − π ∆ (cid:48) , (14)where τ H = 1 k ξ v Aλ , Υ = γ τ H λd w , (15) γ is the growth rate, k ξ is the wavenumber in the ξ -direction, λ is the current sheet half-width, v Aλ is theAlfv´en speed based on the reconnecting magnetic field,and Γ( z ) indicates the gamma function. This dispersionrelation matches the non-relativistic one (Porcelli 1991)when w →
1, i.e., when the plasma is cold. Eq. (14)can be further simplified for short-wavelength modes(small-∆ (cid:48) ) and long-wavelength modes (large-∆ (cid:48) ), whichis convenient in order to derive analytically the condi-tions for current sheet disruption. For Υ (cid:28)
1, the small-∆ (cid:48) regime, the growth rate and the inner tearing layerhalf-width (where the ideal MHD approximation breaksdown due to the finite electron and positron inertia) ofthe instability are γ s = (cid:20) Γ( )2 π Γ( ) (cid:21) (∆ (cid:48) λ ) τ H (cid:18) d w λ (cid:19) , δ in = d w ∆ (cid:48) . (16)On the other hand, for Υ → − , in the large-∆ (cid:48) regime,the growth rate and the inner tearing layer half-widthare γ l = 1 τ H (cid:18) d w λ (cid:19) , δ in = d w . (17)Using these relations, together with the tearing stabilityindex ∆ (cid:48) = 2 λ (cid:18) k ξ λ − k ξ λ (cid:19) , (18)it can be shown that the dominant tearing mode at cur-rent sheet disruption scales like the fastest growing mode(see Comisso et al. (2016, 2017); Huang et al. (2017)),so that the instability wavenumber at current sheet dis-ruption turns out to be simply k ξ,d ∼ d w λ d , (19) For the purpose of this study we have not considered obliquetearing modes, which can be included in a more general dispersionrelation. Here we assume a Harris-type current sheet (Harris 1962),which is a reasonably good approximation of current sheets oc-curring in magnetized turbulence (e.g. Servidio et al. 2010) andcoalescing magnetic islands (e.g. Huang et al. 2017)). where the subscript “ d ” denotes current sheet disrup-tion. This implies also that the growth rate and theinner tearing layer half-width at current sheet disrup-tion are γ d ξv Aλ ∼ d w λ d ξ , δ in ,d ∼ d w . (20)From these expressions, one still needs to determine thecurrent sheet half-width at disruption, λ d , in order toknow the wavenumber k ξ,d and the growth rate γ d . Wecalculate the width of the current sheet at disruption byusing the principle of least time introduced in Comissoet al. (2016, 2017), substituting the resistive tearingmode dispersion relation with the collisionless disper-sion relation discussed above. Then, for a rapid currentsheet that forms on the Alfv´enic timescale, the modethat becomes nonlinear in the shortest time disrupts thecurrent sheet when λ d ξ (cid:34) ln (cid:32) (cid:15) / (cid:18) d w ξ (cid:19) α/ (cid:18) ξλ d (cid:19) / α (cid:33)(cid:35) / (cid:39) c λ (cid:18) d w ξ (cid:19) / , (21)where c λ is an O (1) constant, ˆ (cid:15) = (cid:15)/ ( δB λ ξ ) is a nor-malized amplitude of the noise that seeds the instability(evaluated at the disruption scale), δB λ is the character-istic magnetic field fluctuation at scale λ , and α is an in-dex that depends on the spectrum of the noise, which isrelated to the turbulence spectrum as P B ( k ξ ) ∝ k ξ − α (Comisso et al. 2018). Eq. (21) can be solved exactly interms of the Lambert W function, but here we prefer toconsider an asymptotic solution that yields more trans-parent results. Therefore, we solve Eq. (21) by iterationobtaining, at the first order, the solution λ d ∼ d / w ξ / (cid:34) ln (cid:32) (cid:15) / (cid:18) d w ξ (cid:19) − α (cid:33)(cid:35) − / , (22)which gives us the critical current sheet width that de-termines the layer disruption and the onset of fast re-connection. Finally, the growth rate of the instabilitywhen the current sheet reaches this ratio is γ d ∼ v Aλ ξ ln (cid:32) (cid:15) / (cid:18) d w ξ (cid:19) − α (cid:33) , (23)while the wavenumber of the dominant mode becomes k ξ,d ∼ d − / w ξ − / (cid:34) ln (cid:32) (cid:15) / (cid:18) d w ξ (cid:19) − α (cid:33)(cid:35) / . (24)We obtain from Eq. (22) that λ d (cid:29) λ ∞ (cid:39) d w forouter-scale current sheets with ξ ∼ (cid:96) (cid:29) d w , as it is7expected under typical astrophysical conditions. There-fore, an outer-scale current sheet disrupts in a chain ofplasmoids before reaching the kinetic scale d w , whileinter-plasmoid layers, being shorter, can reach the thick-ness d w . Eq. (22) tells us also that current sheets disruptat a larger thickness for larger noise levels. However, thedependence is only logarithmic. From the other two re-lations, Eq. (23) and Eq. (24), we have that the growthrate of the plasmoid instability is γ d ξ/v Aλ (cid:29) ∝ k ξ,d (cid:96) , increases as (cid:96)/d w increases and the noise of the system decreases. As anexample, we show in Figs. 14 and 15 that a larger num-ber of plasmoids forms when the domain is increased bya factor 4 with respect to the reference 2D simulation.In this simulation, as we argue below in this section, ef-ficient plasmoid formation keeps the reconnection speedand the injection efficiency high when increasing systemsize. As a result, the fraction of nonthermal particles re-mains about the same when moving from the referencebox size L/d e = 1640 up to L/d e = 6560 (see Fig.2(b) in Comisso & Sironi (2018)).When the reconnection layer becomes dominated bythe presence of plasmoids, soon after the condition λ ∼ λ d is met, the complexity of the dynamics givesrise to a strongly time-dependent process. Nevertheless,in a statistical steady-state, we may expect that the re-connection layer containing the main X -point, which isthe one that determines the global reconnection rate,has a bounded aspect ratio ξ X /λ X . If ξ X /λ X ∼
1, thereconnection process would choke itself off, since thiswould imply β R ∼ ξ X /λ X (cid:29) X -point cannot be longer than the marginallystable sheet. Indeed, the fractal-like process of currentsheet disruption due to the plasmoid instability termi-nates when the length of the innermost local currentlayer of the chain is shorter than the critical length ξ c (Huang & Bhattacharjee 2010; Uzdensky et al. 2010;Comisso et al. 2015; Comisso & Grasso 2016). There-fore, ξ X (cid:46) ξ c is also expected. At present there are noanalytical estimates for the aspect ratio ξ c /λ X , whichmight also depend on the noise level (e.g. Ni et al. 2010;Huang et al. 2017; Shi et al. 2018). However, numer-ical simulations have found ξ c /λ X ∼
50 in the colli-sionless regime (e.g. Daughton et al. 2006; Daughton& Karimabadi 2007; Ji & Daughton 2011). As a con-sequence, for a compression ratio between inflow and outflow of order unity, the reconnection rate is boundedfrom above and below as1 / (cid:46) β R (cid:28) , (25)which classify it as a fast reconnection rate. More pre-cisely, numerical simulations consistently indicates that β R is an O (0 .
1) quantity (for relativistic pair plasmas,see, e.g, Zenitani et al. (2009); Bessho & Bhattachar-jee (2012); Cerutti et al. (2012); Guo et al. (2014); Ka-gan et al. (2015); Liu et al. (2015); Sironi et al. (2016);Werner & Uzdensky (2017); Liu et al. (2017)).The aforementioned properties of reconnecting cur-rent sheets are important in regulating the particle in-jection efficiency. Here we show that the fraction ofparticles processed by reconnecting current sheets is in-dependent of the system size and is quite large (despitethe small filling fraction of current sheets) as long as thereconnection rate is high. To this purpose, let us con-sider a generic current sheet of characteristic length 2 ξ and thickness 2 λ , whose lifetime is approximately givenby the local eddy turnover time τ nl ∼ τ Aξ = ξ/v Aλ ,assuming critical balance (Goldreich & Sridhar 1995;Boldyrev 2006). If fast reconnection occurs for a timeclose to the eddy turnover time (see, e.g, Fig. 15), asingle reconnecting current sheet can “process” the up-stream plasma up to a distance λ R,j = β R,j c τ nl ,j ∼ β R,j c ξ j v Aλ,j , (26)where j labels the j -th current sheet among the popu-lation of current sheets present at a given time, and thesubscript “ R ” stands for reconnection. Since the surfaceprocessed by the entire population of reconnecting cur-rent sheets is in good approximation the one processedby the largest-scale ones, whose length scale correspondsto the turbulence integral length (cid:96) (e.g. Servidio et al.2009), we have that magnetic reconnection can processa plasma surface A R = (cid:88) j λ R,j ξ j ∼ β R L (27)in one large-eddy turnover time. Here, we have used n cs ∼ ( L/(cid:96) ) as an estimate for the number of outer-scale current sheets, and β R is the average reconnectionrate. Furthermore, if we consider that current sheets in3D are sheet-like structures with 2 λ (cid:28) ξ (cid:46) l (cid:107) , with2 l (cid:107) indicating the direction along the magnetic field, wecan obtain that, in one large-eddy turnover time, thereconnecting current sheets process a plasma volume V R = (cid:88) j λ R,j ξ j l (cid:107) ,j ∼ β R L . (28)Therefore, according to Eqs. (27)/(28), the plasmasurface/volume processed by the reconnecting current8sheets is a fixed fraction of the domain if β R is indepen-dent of the system size, as discussed above. Moreover,since β R is an O (0 .
1) quantity, magnetic reconnectioncan process large volumes of magnetic energy in fewouter-scale eddy turnover times.In the next sections, we will address how particles areenergized both in the injection phase and in the subse-quent stochastic acceleration phase, and we will analyzethe signatures of the acceleration process on the particledistribution. MECHANISMS OF PARTICLE ENERGIZATIONIn order to distinguish the relative roles of differentenergization mechanisms, it is convenient to computethe work done by the parallel electric field, W (cid:107) ( t ) = q (cid:82) t E (cid:107) ( t (cid:48) ) · v ( t (cid:48) ) dt (cid:48) , as well as the work done by the per-pendicular electric field, W ⊥ ( t ) = q (cid:82) t E ⊥ ( t (cid:48) ) · v ( t (cid:48) ) dt (cid:48) ,for a statistically significant sample of particles (here, asusual, q is the electric charge, E is the electric field, and v is the particle velocity). To this aim, we tracked a sam-ple of ∼ particles randomly selected from each of ourPIC simulations. Note that in this section, parallel ( (cid:107) )and perpendicular ( ⊥ ) components are defined with re-spect to the local magnetic field, i.e. E (cid:107) = ( E · B ) B /B and E ⊥ = E − E (cid:107) . The main results of our analysis, forthe reference 2D and 3D simulations (see Table 1), arepresented in Fig. 16 (left column for 2D and right col-umn for 3D). We first discuss the energization process ofrepresentative particles that end up in the high-energytail, and then we present a statistical analysis that al-lows us to quantify the contributions of parallel and per-pendicular electric fields for the overall acceleration ofnonthermal particles.Figs. 16(a) and 16(b) show the particle energy gainnormalized to rest mass energy, ∆ γ ( t ) = γ ( t ) − γ (0),as well as the relative contributions W (cid:107) ( t ) /mc and W ⊥ ( t ) /mc , for representative high-energy particles in2D and 3D turbulence. The total work done by the elec-tric field is not plotted here, since q (cid:82) t E ( t (cid:48) ) · v ( t (cid:48) ) dt (cid:48) = mc ∆ γ ( t ) is satisfied to high accuracy and is essentiallyindistinguishable from the black solid line representing∆ γ ( t ). Both figures indicate that the work done by E (cid:107) is responsible for the initial energy gain, while the workdone by E ⊥ takes over at relatively low energies andpropels the particle to the highest energies. Alternativeplots that provide similar information, but can be moreeasily generalized to analyze a large population of par-ticles (as we do below), are shown in Figs. 16(c) and16(d), for 2D and 3D, respectively. In this case, the rel- W (cid:107) ( t ) and W ⊥ ( t ) are computed on the fly in order to achievehigh accuracy, regardless of the time sampling of particle outputs. ative contributions W (cid:107) /mc and W ⊥ /mc are plotted asa function of ∆ γ , and the black solid line indicates theexpected sum of the two terms. The plots show thatthe low ∆ γ -range is dominated by W (cid:107) , while W ⊥ (cid:29) W (cid:107) when particles reach high energies.Figs. 16(c) and 16(d) are generalized in Figs. 16(e)and 16(f), respectively, to account for a statistical as-sessment of the energization of a sample of particles.We consider all tracked particles that end up well intothe nonthermal tail at late times, more precisely alltracked particles for which γ ≥ σ at ct/l = 12. Thefigures show the distribution f (cid:0) ∆ γ, W (cid:107) /mc (cid:1) of parti-cles with respect to ∆ γ and W (cid:107) /mc . We normalize f (cid:0) ∆ γ, W (cid:107) /mc (cid:1) such that ∞ (cid:90) −∞ f (cid:18) ∆ γ, W (cid:107) mc (cid:19) d (∆ γ ) = 1 . (29)The distribution f (cid:0) ∆ γ, W ⊥ /mc (cid:1) is not plotted heresince it conveys the same message. We can see thatthe peak of the distribution for a given ∆ γ is around W (cid:107) /mc ∼
40, for all ∆ γ >
50. This occurs both in 2Dand 3D simulations. We also calculated the median ofthe histogram as a function of ∆ γ , which is shown asa black dashed line in Figs. 16(e) and 16(f). The me-dian also approaches a constant value W (cid:107) /mc ∼ γ > This confirms for a statistically-significant sample of particles the same conclusions pre-sented above: high-energy particles are first energizedvia v · E (cid:107) , which brings them up to ∆ γ ∼ W (cid:107) /mc ∼ W ⊥ (cid:29) W (cid:107) for the highest energyparticles.In summary, we find both for individual particles andfor a large sample of particles, that the initial stagesof acceleration are controlled by parallel electric fields.This is consistent with the fact that strong parallel elec-tric fields are expected at active reconnection layers,where we have indeed shown that particle injection (i.e.,the first stage of acceleration) occurs.The initial energy gain due to v · E (cid:107) is dependent onmagnetization. It can be seen from Fig. 17 that thetypical energy gain provided by v · E (cid:107) increases with σ .From our simulations we find that the typical increasein Lorentz factor during the injection process (which is For low values of ∆ γ , the mode and the median of f (cid:0) W (cid:107) /mc | ∆ γ (cid:1) are independent of the final particle energy onlyif the selected threshold satisfies γ (cid:29) ( σ / γ th at late times,i.e., for particles that end up well into the nonthermal tail (seealso Sec. 6). ct/l ∆γ W / mc || W / mc ⊥ ct/l ∆γ W / mc || W / mc ⊥ ∆γ ∆γ W / mc || W / mc ⊥ ∆γ ∆γ W / mc || W / mc ⊥ ∆γ median[ f ( W / mc | ∆γ )] || ∆γ ∆γ median[ f ( W / mc | ∆γ )] || ∆γ f ( ∆γ , W / mc ) || ∆γ median[ f ( W / mc | ∆γ )] || ∆γ ∆γ median[ f ( W / mc | ∆γ )] || ∆γ f ( ∆γ , W / mc ) || Figure 16.
Relative contributions of E (cid:107) = ( E · B ) B /B and E ⊥ = E − E (cid:107) to the particle energization in 2D (left) and 3D(right) simulations with σ = 10 and δB rms0 /B = 1. The 2D simulation has domain size L/d e = 1640 (with l = L/ L/d e = 820 (with l = L/ γ (black solid line), normalized work done by the parallel electric field, W (cid:107) /mc (red solid line),and normalized work done by the perpendicular electric field, W ⊥ /mc (blue solid line). Middle row: scatter plot of W (cid:107) /mc versus ∆ γ (red triangles) and W ⊥ /mc versus ∆ γ (blue diamonds), for the same particle displayed in the top frame. The solidblack line indicates the expected sum W (cid:107) /mc + W ⊥ /mc = ∆ γ . Bottom row: distribution of particles with respect to ∆ γ and W (cid:107) /mc , for particles ending up with γ ≥ σ at ct/l = 12. The median of the conditional PDF at given ∆ γ , f (cid:0) W (cid:107) /mc | ∆ γ (cid:1) ,is shown with a dashed black line. Again, the solid black line indicates the expected sum W (cid:107) /mc + W ⊥ /mc = ∆ γ . ∆γ / σ σ = σ = σ = σ =
400 5 10 15 ∆γ / σ m e d i a n [ f ( W || / m c | ∆ γ )] / σ ∆γ / σ σ = σ = σ = σ =
400 5 10 15 ∆γ / σ m e d i a n [ f ( W || / m c | ∆ γ )] / σ Figure 17.
Median of f (cid:0) W (cid:107) /mc | ∆ γ (cid:1) , divided by σ , forhigh-energy particles from different 3D simulations having σ = 5 (blue), σ = 10 (green), σ = 20 (orange), and σ =40 (red). We considered all tracked particles with γ ≥ σ at ct/l = 12. All the simulations have δB rms0 /B = 1 and L/d e = 820. The solid black line indicates the expectedsum W (cid:107) /mc + W ⊥ /mc = ∆ γ . governed by parallel fields at reconnection layers) is∆ γ inj ∼ W (cid:107) /mc ∼ κσ γ th , σ (cid:29) κ is a numerical factor of order unity ( κ ∼ σ = δB / πn wmc decreases with time in decayingturbulence, implying that the time-dependent ∆ γ inj = κσγ th also decreases with time ( γ th is the instantaneousmean Lorentz factor).The length l (cid:107) along B (which in reconnection layers isdominated by the mean field B ) required to attain theenergy gain ∆ γ inj can be obtained from the reconnectionelectric field E R by assuming particles moving along E (cid:107) at the speed | v | ∼ c . If E (cid:107) ∼ E R ≈ const during theacceleration time, then∆ γ inj = β R δB rms e l (cid:107) mc . (31)Here we have used the typical reconnection electric field E R = β R δB rms v A /c ∼ β R δB rms . Therefore, the lengthscale l inj required to attain the increase ∆ γ inj can beexpressed as l inj = κβ R (cid:114) σw γ th d e , (32)where the different physical quantities have to be eval-uated at the injection time. This expression indicatesthat a sufficient length for particle injection is alwaysguaranteed for a large enough system, i.e., l (cid:29) l inj . Sim-ilarly, as most of injection happens at outer-scale current -1 γ - -6 -5 -4 -3 -2 -1 d N / d l n ( γ - ) p=2.8 self-consistent particlestest particles with no E || Figure 18.
Particle spectra dN/d ln( γ −
1) at ct/l = 12for normal particles (red solid line) and test-particles thatare evolved with E → E − E (cid:107) (blue solid line) from a 2Dsimulation with σ = 10, δB rms0 /B = 1, and L/d e = 1640.The two spectra display similar power-law index but differentpower-law normalization. The high-energy tail with γ ≥ .
2% of thetest-particles are contained in the tail with γ ≥ sheets, the time τ inj required for accelerating particlesup to this energy is always granted if τ nl = l/δV rms (cid:29) τ inj = l inj /c , where τ nl is the outer-scale nonlinear timeand δV rms = (cid:104) δV (cid:105) / is the space-averaged root-mean-square value of the velocity field fluctuations. The tworequirements coincide for δV rms → c .Even though W ⊥ (cid:29) W (cid:107) for high-energy particles, theinitial v · E (cid:107) energization process is important to pro-mote the particles to energies large enough such thatthey can experience the subsequent v · E ⊥ acceleration.Hence, energization by v · E (cid:107) controls the number ofparticles that have the possibility to reach nonthermalenergies. This point, which was already discussed inSection 4, can be probed in a direct way by compar-ing the self-consistent PIC particles with a populationof test-particles for which we artificially exclude accel-eration by E (cid:107) , assuming that the electric field they feelis E → E − E (cid:107) . To this aim, we performed a 2D PICsimulation where we added a population of ∼ × such test-particles. The resulting particle spectra at latetime ( ct/l = 12) are shown in Fig. 18. The particlespectrum of normal particles has a much larger fractionof particles contained in the high-energy tail (17% vs0 . Equivalently, the normalization of the power-law tail in the test-particle spectrum is much lower than In both cases we consider a nonthermal tail starting at γ ≥
12, since the power-law range starts at γ ∼
12 for both particlespectra. Note that this value is quite larger than the termal peak p = − d log N/d log( γ −
1) of the power-law tailis similar (see dashed black lines), indicating that the v · E ⊥ energization is the crucial process responsible forsetting the power-law slope. Also, the cutoff energy isabout the same for the two population of particles, indi-cating that it is not controlled by parallel electric fields.In the next section we will see that the two differ-ent energization processes, which dominate in differentenergy ranges ( v · E (cid:107) for ∆ γ inj (cid:46) κσγ th and v · E ⊥ for∆ γ inj (cid:38) κσγ th ), also affect the anisotropy of the particledistribution. ANISOTROPY AND PARTICLE MIXINGAnisotropic features of the particle distribution canhave a significant impact on the observed synchrotronemission (e.g. Tavecchio & Sobacchi 2019). Here weshow that even if the initial velocity distribution isisotropic, the particle energization process drives a sig-nificant energy-dependent anisotropy, as the pitch anglescattering rate is not sufficient to keep the particle dis-tribution close to isotropy. In order to characterize theanisotropy of the particle distribution, we first examinethe pitch-angle distribution of the particles, namely thestatistics of the pitch-angle cosine cos α = v · B / ( | v | | B | ).Then, we analyze the anisotropy of the four-velocity dis-tribution of the particles. We perform these analyses ona statistically significant sample of ∼ particles, bothin 2D and 3D. These particles were randomly selectedand tracked over time for each of the simulations. Fi-nally, we also look at the spatial mixing of particles.6.1. Pitch-angle distribution
The time evolution of the overall particle distributionwith respect to cos α is shown in Fig. 19(a) for the ref-erence 2D simulation and in Fig. 19(b) for the reference3D simulation. As turbulence evolves, the pitch-angledistribution becomes anisotropic with strong peaks atcos α = ±
1, i.e., for particles moving along the mag-netic field lines. Pronounced peaks of the pitch-angledistribution near cos α = ± β p (e.g. Pecoraet al. 2018), with β p = n k B T / ( (cid:104) B (cid:105) / π ) indicating theplasma beta, i.e., the ratio of thermal pressure to mag-netic pressure. Indeed, the low- β p regime is similar tothe high- σ regime investigated here, in the sense that inboth cases the magnetic energy density dominates overthe thermal energy density (in our simulations the initialplasma beta is β p = 2 θ / [ w ( σ z + σ )]). The PDFs of of the test-particle population, which is consistent with the lownormalization N of its nonthermal power-law tail.
2D and 3D simulations are similar; nevertheless, the 3Dcase exhibits higher probability peaks near cos α = ± E (cid:107) . As we have shown, reconnecting current sheets canprocess a large fraction of particles in just a few c/l (seeEqs. (27) and (28)).The PDFs illustrated in Fig. 19 are dominated bylow-energy particles (i.e., near the spectral peak), sincethey control the number census (see Figs. 5 and 6).In order to characterize the anisotropy of particles ofhigher energy, we construct PDFs of cos α for differ-ent populations of particles depending on their Lorentzfactor. We collected particle data from a time range ct/l ∈ [3 ,
12] in order to increase statistics. However, wehave verified that decreasing this range (up to a singletime snapshot taken at late times) does not modify theresults. These results are shown in Fig. 20(a) and Fig.20(b) for 2D and 3D turbulence, respectively. At verylow energies ( γ ∼ γ ∼ α = ±
1, in analogywith the results shown in Fig. 19. At higher Lorentzfactors, the pitch-angle distribution evolves into a “but-terfly distribution” with minima at both cos α = ± α = 0. This phenomenon occurs at Lorentz factors γ ∼
50 for the simulations with σ = 10 shown in Fig.19. At even higher energies ( γ (cid:29) α = 0, i.e.for particles moving in the plane perpendicular to the lo-cal magnetic field. This trend can be displayed using adistribution f (cos α, γ ) of particles with respect to cos α and γ . This distribution, shown in Fig. 20(c) and Fig.20(d) for 2D and 3D turbulence, respectively, has beennormalized such that (cid:90) − f (cos α, γ ) d (cos α ) = 1 . (33)In these plots, the peaks of f (cos α, γ ) are located atcos α = ± α = 0 until γ ∼
80. At higher energies, the peakof the distribution remains located at cos α = 0, withparticles that lie progressively more perpendicular to thelocal magnetic field as their energy increases.2 -1.0 -0.5 0.0 0.5 1.00.00.51.01.52.02.5-1.0 -0.5 0.0 0.5 1.0cos α P D F ct/l (a) -1.0 -0.5 0.0 0.5 1.00.00.51.01.52.02.53.0-1.0 -0.5 0.0 0.5 1.0cos α P D F ct/l (b) Figure 19.
Probability density functions of the pitch-angle cosine cos α = v · B / ( | v | | B | ) at different times, obtained from2D (left) and 3D (right) simulations. Both simulations have σ = 10 and δB rms0 /B = 1. The 2D simulation has domain size L/d e = 1640 (with l = L/ L/d e = 820 (with l = L/ -1.0 -0.5 0.0 0.5 1.00.00.51.01.52.02.5-1.0 -0.5 0.0 0.5 1.0cos α P D F [1,1.5] [15,20] [35,50] [60,80] [150,200] γ -intervals (a) -1.0 -0.5 0.0 0.5 1.00.00.51.01.52.02.5-1.0 -0.5 0.0 0.5 1.0cos α P D F [1,1.5] [15,20] [35,50] [60,80] [150,200] γ -intervals (b) -1.0 -0.5 0.0 0.5 1.0cos α γ -1.0 -0.5 0.0 0.5 1.050100150200250300 γ cos α , γ )( f -1.0 -0.5 0.0 0.5 1.0cos α γ -1.0 -0.5 0.0 0.5 1.050100150200250300 γ cos α , γ )( f Figure 20.
Particle distributions obtained from 2D (left) and 3D (right) simulations with σ = 10 and δB rms0 /B = 1. The2D simulation has domain size L/d e = 1640 (with l = L/ L/d e = 820 (with l = L/ α = v · B / ( | v | | B | ) for particle Lorentz factors inthe intervals γ ∈ [1 , .
5] ( × symbol), γ ∈ [15 ,
20] ( ◦ symbol), γ ∈ [35 ,
50] ( (cid:3) symbol), γ ∈ [60 ,
80] ( (cid:5) symbol), and γ ∈ [150 , (cid:77) symbol). Bottom row: particle distribution with respect to the pitch-angle cosine cos α and the Lorentz factor γ . The plotsare obtained from data in the time range ct/l ∈ [3 ,
12] to increase statistics. -1.0 -0.5 0.0 0.5 1.001234-1.0 -0.5 0.0 0.5 1.0cos α P D F σ = σ = σ = σ = γ ∈ [0.8 σ , σ ] γ ∈ [4 σ , σ ] γ ∈ [16 σ , σ ] Figure 21.
Probability density functions of the pitch-anglecosine cos α = v · B / ( | v | | B | ) for particles with Lorentz fac-tors γ ∈ [0 . σ , . σ ] (solid lines), γ ∈ [4 σ , σ ] (long-dashed lines), and γ ∈ [16 σ , σ ] (dashed lines). Differentcolors refer to different 3D simulations having σ = 5 (blue), σ = 10 (green), σ = 20 (orange), and σ = 40 (red). All3D simulations have δB rms0 /B = 1 and L/d e = 820. Wealso recall that γ th ≈ .
58. Data is collected from a timerange ct/l ∈ [3 , The energy-dependent anisotropy illustrated in Fig.20 reflects the different acceleration mechanisms thatoperate at different energies (see Section 5). At lowenergies, the contribution of the v (cid:107) · E (cid:107) energizationis dominant, so that particles end up being stronglyaligned/antialigned with the magnetic field (cos α ∼± v ⊥ · E ⊥ energization takes over and propels the parti-cles in the direction perpendicular to the local magneticfield. The time scale of this acceleration is fast comparedto the pitch-angle scattering timescale, so that particlesretain their orientation cos α ∼ σ =10 hold also for the other magnetizations we investigate.In Fig. 21, we present the results from four simulationsthat differ in magnetization σ ∈ { , , , } . Herewe show only the results from 3D simulations, sincethose from 2D simulations are analogous. The rangesin γ are scaled with σ , which provides the typical en-ergy scale (e.g., the starting point of the high-energynonthermal tail, γ st , increases linearly with σ , as illus-trated by Eq. (8)). For γ ∼ ( σ / γ th (solid lines),we have a pitch-angle distribution peaked at cos α ∼ ± σ ). The butterfly distribution with min-ima at cos α = ± , γ ∼ σ / γ th (long-dashed lines). Finally, for γ (cid:29) σ / γ th , well into thenonthermal tail, the particle velocities become mostly perpendicular to the magnetic field and we can see thatall the distributions are peaked at cos α = 0 (see dashedlines). 6.2. Particle four-velocity distribution
The results on the anisotropy of the pitch-angle dis-tributions, computed with respect to the local magneticfield B = B + δ B , suggest that the four-velocity distri-bution function, with respect to the mean magnetic field B = B ˆ z , should also display significant anisotropy.Indeed, as the turbulence fluctuations decay, the localmagnetic field becomes progressively more aligned withthe direction of the mean magnetic field.We calculated the domain-averaged four-velocity dis-tributions in the xy plane, f (cid:0) γβ x , γβ y (cid:1) , and in the xz plane, f (cid:0) γβ x , γβ z (cid:1) , from the same samples of particlesused to analyze the local pitch-angle distributions. Theresults, for our reference 3D simulation (the 2D case isanalogous) are shown in Fig. 22 (and a zoom in Fig. 23).As for Fig. 20, we collected particle data in the timerange ct/l ∈ [3 ,
12] in order to increase statistics, but wehave also verified that decreasing this range (up to a sin-gle time snapshot taken at late times) does not modifythe results. As we expected, from Fig. 22 we find thatthe four-velocity distribution is isotropic in the planeperpendicular to B (top panel), while it develops morecomplex features with respect to planes that contain B , as for the case of f (cid:0) γβ x , γβ z (cid:1) (bottom panel). Theresults are analogous when considering f (cid:0) γβ y , γβ z (cid:1) or f (cid:0) ( γ β x + γ β y ) / , γβ z (cid:1) . The distribution f (cid:0) γβ x , γβ z (cid:1) displays a core region elongated in the γβ z direction, asparticle velocities are mostly aligned/antialigned withthe magnetic field at low energies. Furthermore, a closeinspection shows that there is an intermediate-energyregion with the majority of particles residing in a dou-ble cone whose axis is the direction of the mean mag-netic field B (see Fig. 23). This is the intermediate-energy range in which the peak of pitch-angle distribu-tion moves from cos α = ± α = 0. At evenhigher energies, Fig. 22 shows that f (cid:0) γβ y , γβ z (cid:1) becomeselongated in the direction perpendicular to B , consis-tently with the dominance of v ⊥ · E ⊥ energization athigher energies and the resulting anisotropy of the pitchangle cosine. 6.3. Particle mixing
We show that while the qualitative and quantitativefeatures of the pitch-angle distributions are similar inour 2D and 3D simulations, the turbulent mixing (inspace) of the energized particles is quite different. Par-ticle mixing in 2D is expected to be less efficient thanin 3D, since the translation-invariant symmetry along4 -100 200-200 0 100-200-10001002000.5 7.0 log [ f ( , )] γβ x y γβγβ x y γ β -100 200-200 0 100-200-10001002000.5 7.0 log [ f ( , )] γβ x z γβγβ x z γ β Figure 22.
Top frame: box-averaged four-velocity distribu-tion function f ( γβ x , γβ y ) for a 3D simulation with σ = 10, δB rms0 /B = 1, and L/d e = 820. Bottom frame: fromthe same simulation, box-averaged four-velocity distributionfunction f ( γβ x , γβ z ). The plots are obtained from data inthe time range ct/l ∈ [3 , log [ f ( , )] γβ x z γβ -20 40 -40 0 20 -40 -200 20 40 γβ x z γ β Figure 23.
Zoom around the intermediate-energy region ofthe four-velocity distribution function f ( γβ x , γβ z ) shown inFig. 22. ˆ z seriously constrains the 2D dynamics. As a conse-quence, regions of space devoid of high-energy particlescan be retained for a larger number of outer-scale eddyturnover times in 2D simulations.In Fig. 24, we show how the energized particles aredistributed in the spatial domain in 2D (left column)and 3D (right column). For both cases, we plot thecell-averaged kinetic energy per particle, (cid:104) γ − (cid:105) cell , atthree different times (from top to bottom). The meankinetic energy is normalized by mc . The time snap-shots are different for 2D ( ct/l =4.6,7.7,10.8) and 3D( ct/l =2.7,4.5,6.3), to account for the faster turbulencedecay in 3D. In both cases, the initial energization occursat current sheets, which display high values of (cid:104) γ − (cid:105) cell ,and then particles propagate outside current sheets inother regions of the domain (see top frames in Fig. 24).As time progresses, energized particles diffuse in thespatial domain, and the mean kinetic energy per par-ticle becomes more uniform (middle frames in Fig. 24).However, in 2D the cores of the large-scale flux tubes re-main essentially unaffected, as these overdense regionswith n (cid:29) n and with higher fluctuation magnetic en-ergy density δB / π (see Fig. 1) are mainly populatedby low energy particles that have not been processedby reconnecting current sheets. On the other hand, the3D domain does not present such isolated regions of low (cid:104) γ − (cid:105) cell . At quite early times, the mean kinetic en-5 y / l x/l 〈 γ - 〉 ce ll ≥ y / l x/l 〈 γ - 〉 ce ll ≥ y / l x/l 〈 γ - 〉 ce ll ≥ y / l x/l 〈 γ - 〉 ce ll ≥ y / l x/l 〈 γ - 〉 ce ll ≥ y / l x/l 〈 γ - 〉 ce ll ≥ Figure 24.
2D plots of the cell-averaged mean kinetic energy per particle normalized by mc , (cid:104) γ − (cid:105) cell , for 2D turbulence(left column) and 3D turbulence (right column). For 3D turbulence, the 2D plots refer to a slice of the domain at constant z/l = 0. The normalized times ct/l for the plots in the left column are (from top to bottom): ct/l = 4 . ct/l = 7 . ct/l = 10 .
8, while those for the plots in the right column are (from top to bottom): ct/l = 2 . ct/l = 4 . ct/l = 6 .
3. The2D simulation has a domain size
L/d e = 1640 (with l = L/ L/d e = 820 (with l = L/ σ = 10 and δB rms0 /B = 1. An animation showing (cid:104) γ − (cid:105) cell at ct/l = 2 . xy slices can be foundat https://doi.org/10.7916/d8-prt9-kn88. (cid:104) γ − (cid:105) cell for much longer (bottom framesin Fig. 24). This different behavior between 2D and 3Dturbulence is also reflected in the particle energy spec-trum, with 2D turbulence retaining more particles with γ (cid:46) γ th until late times (see Figs. 5 and 6). PARTICLE ENERGY DIFFUSION ANDSTOCHASTIC ACCELERATIONWe have seen that after the injection phase, the subse-quent energy gain is dominated by perpendicular electricfields via stochastic scattering off the turbulent fluctua-tions (Comisso & Sironi 2018). Here, in order to eluci-date the properties of the stochastic acceleration phase,we evaluate the energy diffusion coefficient directly fromthe self-consistent particle evolution of our PIC simu-lations. This allows us to determine the accelerationtimescale associated with stochastic acceleration. Thenwe show that the two-stage process that accelerates par-ticles is well modeled by an initial injection phase pow-ered by reconnection electric fields, followed by a secondacceleration phase modeled with the measured energydiffusion coefficient.7.1.
Particle energy diffusion
Particles that are stochastically scattered off the tur-bulent fluctuations experience a biased random walk inmomentum space, which can be modeled with a Fokker-Planck approach (e.g. Blandford & Eichler 1987), pro-vided that the fractional momentum change in a singlescattering is sufficiently small. In this case, one coulddescribe the process of stochastic acceleration from thepoint of view of a Fokker-Planck equation in energyspace (e.g. Ramaty 1979) ∂N∂t = − ∂∂γ ( A γ N ) + ∂ ∂γ ( D γ N ) . (34)Here, as usual, N is the particle spectrum differentialin energy, A γ is the energy convection coefficient, and D γ is the energy diffusion coefficient. Note that in gen-eral, the convection and diffusion coefficients are timedependent (this is indeed the case for the turbulencesimulations performed here). The convection coefficient A γ represents the mean energy gain due to stochasticacceleration, and is related to the diffusion coefficient inenergy space as A γ = d (cid:104) γ (cid:105) dt = 1 γ ∂∂γ (cid:0) γ D γ (cid:1) , (35)The diffusion coefficient in energy space D γ is also re-lated to the diffusion coefficient in momentum space D p , c ∆ t/l 〈 ( ∆ γ ) 〉 γ * γ / γ σ D γ σ = σ = σ = σ = ∝ γ Figure 25.
Diffusion in energy space from 2D simulationswith δB rms0 /B = 1 and different initial magnetizations σ . Top panel: mean square variation of the Lorentz fac-tor for particles binned in logarithmic intervals [ γ ∗ /ν, γ ∗ ν ]with ν = 1 . γ ∗ = 21 . →
84 (from blue to red) attime ct ∗ /l = 5 .
25 for the reference 2D simulation. Thedashed black lines indicate linear fits. Bottom panel: en-ergy diffusion coefficient D γ (in units of c/l ), as a func-tion of the Lorentz factor γ (divided by γ σ to align caseswith different magnetization), measured at the time interval c ∆ t/l = 1 .
875 from four simulations having initial magneti-zation σ ∈ { , , , } . with D γ (cid:39) D p for the ultra-relativistic particles con-sidered here. Given the fact that high-energy particlespreferentially lie in the plane perpendicular to the meanfield (Sec. 6), and that their energization is mostly con-tributed by perpendicular electric fields (Sec. 5), themomentum diffusion coefficient D p is essentially identi-cal to D p ⊥ , i.e., to the diffusion coefficient of momentaperpendicular to the mean field. The determination ofthis coefficient, or equivalently of the energy diffusioncoefficient, establishes the properties of the stochasticacceleration phase.7 c ∆ t/l 〈 ( ∆ γ ) 〉 γγ * γ / γ σ D γ σ = σ = σ = σ = ∝ γ Figure 26.
Diffusion in energy space from 3D simulationswith δB rms0 /B = 1 and different initial magnetizations σ . Top panel: mean square variation of the Lorentz fac-tor for particles binned in logarithmic intervals [ γ ∗ /ν, γ ∗ ν ]with ν = 1 . γ ∗ = 21 . →
84 (from blue to red) at time ct ∗ /l = 3 for the reference 3D simulation. The dashed blacklines indicate linear fits. Bottom panel: energy diffusion coef-ficient D γ (in units of c/l ), as a function of the Lorentz factor γ (divided by γ σ to align cases with different magnetization),measured at the time interval c ∆ t/l = 1 .
875 from four sim-ulations having initial magnetization σ ∈ { , , , } . We evaluate the energy diffusion coefficient directlyfrom PIC simulations (see also Wong et al. (2019)). Tothis aim, from each of the 2D and 3D simulations em-ployed for this analysis, we tracked in time the posi-tions, four-velocities, and electromagnetic field values ofabout 10 particles that were randomly selected at thebeginning of the simulation. From the time history ofthe particles evolution, we calculate the mean square γ -variation (cid:104) (∆ γ ) (cid:105) = 1 N p N p (cid:88) n =1 ( γ n ( t ) − γ n ( t ∗ )) (36) for particles grouped in such a way that at an initial time t ∗ , they belong to the same energy bin ( N p is the numberof particles in the selected bin). The energy bin at t ∗ ischosen accordingly to the particle energy calculated inthe frame comoving with the drift velocity v D = c E × B /B . For each particle, we perform a Lorentz boostfrom the observer/simulation frame to the local E × B frame, which results in the boosted Lorentz factor γ (cid:48) = γ D γ (cid:0) − v D · v /c (cid:1) , where γ D = 1 / (cid:113) − ( v D /c ) is theLorentz factor for the drift velocity. Then, we evaluateEq. (36) by selecting particles in a small energy bin with γ (cid:48) ∈ [ γ ∗ /ν, γ ∗ ν ], where γ ∗ is the characteristic Lorentzfactor of the energy bin, and ν is a constant factor thatshould be close to unity (we choose ν = 1 . D γ = (cid:104) (∆ γ ) (cid:105) t , (37)where ∆ t = t − t ∗ is a time interval that should be (i)long enough that the initial conditions become insignif-icant, and particles are in the diffusive regime; and (ii)short enough that the turbulence properties have notsignificantly changed. By using a large sample of par-ticles in each energy bin, non-secular variations of theparticle energy are averaged out and the mean energygain can be obtained.The results of our analysis of the particle energy dif-fusion are reported in Figs. 25 and 26, for 2D and 3Dsimulations, respectively. The top frames show the meansquare variation (cid:104) (∆ γ ) (cid:105) for particles binned accordingto their initial energy at ct ∗ /l = 5 .
25 for the reference2D simulation and ct ∗ /l = 3 for the reference 3D simu-lation. In both cases, at the selected time t ∗ , turbulenceis well developed and the time dependent magnetiza-tion calculated with the magnetic energy in turbulentfields is σ ( t ∗ ) ∼
1. The plots indicate that a diffusivebehavior in energy space, (cid:104) (∆ γ ) (cid:105) ∝ ∆ t (compare withdashed black lines), is achieved after c ∆ t/l ∼
1, in both2D and 3D reference simulations. For shorter time inter-vals, particles preserve memory of the initial conditionsand their motion is not diffusive. The slope at late times(dashed lines) depends on particle energy, and it allowsto quantify the energy dependence of the diffusion coef-ficient.The bottom frames on Figs. 25 and 26 show the parti-cle energy dependence of the energy diffusion coefficientfrom simulations with different initial magnetization σ (indicated with different colors in the figures). The dif-fusion coefficient is evaluated using Eq. (37) in the timeinterval c ∆ t/l = 1 . ct ∗ /l = 5 .
25 for 2Dsimulations and ct ∗ /l = 3 for 3D. We verified that theenergy dependence remains the same when taking dif-8ferent time intervals, or by fitting the slopes of (cid:104) (∆ γ ) (cid:105) as a function of time in the diffusive regime (as donewith the dashed lines in the top panels). In order toproperly compare different σ , we display the energydiffusion coefficient as a function of the Lorentz fac-tor normalized by γ σ (see Eq. (8)). The energy rangewhere stochastic acceleration occurs starts at the begin-ning of the power-law high-energy tail of the particlespectrum, i.e. for γ/γ σ (cid:38)
1. In the stochastic accel-eration range, the energy diffusion coefficient scales as D γ ∝ γ (compare with the dashed black lines in thebottom panels). A similar dependence on the particleenergy was also found in Lynn et al. (2014); Kimuraet al. (2016, 2019); Wong et al. (2019), and is consis-tent with particle acceleration by non-resonant and/orbroadened resonant interactions with the turbulent fluc-tuations (e.g. Skilling 1975; Blandford & Eichler 1987;Schlickeiser 1989; Chandran 2000; Cho & Lazarian 2006;Lemoine 2019). Then, at higher energies, near the high-energy cutoff of the power law, the energy dependence of D γ becomes weaker as the particle Larmor radius getscloser to the energy-containing scale of the turbulence.The energy diffusion coefficient depends also on theactual magnetization σ ( t ∗ ) = (cid:104) δB (cid:105) / πn wmc . In or-der to better understand this dependence, in Fig. 27 weplot the energy diffusion coefficient as a function themagnetization σ at four different times t ∗ (in the range ct ∗ /l ∈ [4 ,
6] for 2D, and ct ∗ /l ∈ [2 ,
4] for 3D) for the foursimulations having different initial magnetization. Both2D and 3D simulations show a clear trend of increasingdiffusion coefficient with increasing magnetization. The3D simulations are well fitted by a linear relation in σ (compare with dashed black line), D γ ∼ . σ (cid:16) cl (cid:17) γ . (38)This scaling can be understood by noting that for astochastic process akin to the original Fermi mechanism(e.g. Blandford & Eichler 1987; Lemoine 2019), the en-ergy diffusion coefficient is D γ = 13 (cid:104) γ V β V (cid:105) cλ mfp γ , (39)where (cid:104) γ V β V (cid:105) / is the typical four-velocity of the scat-terers, and λ mfp is the particle scattering mean-free-path. Therefore, if we estimate the scattering mean-free-path as λ mfp ∼ ( B /δB rms ) l and identify γ V β V with the dimensionless Alfv´enic four-velocity, (cid:104) γ V β V (cid:105) ∼(cid:104) B (cid:105) / πnwmc , from Eq. (39) we obtain D γ ∝ σ for (cid:104) B (cid:105) /B ∼
1, in agreement with Eq. (38). Then, fromthese results we can also estimate the stochastic accel- σ ( t ∗ ) σ = σ = σ = σ = σ D γ / γ σ ( t ∗ )10 -2 -1 D γ / γ σ = σ = σ = σ = σ Figure 27.
Diffusion coefficient in energy space as a func-tion of the actual magnetization σ ( t ∗ ) from 2D simulations(top) and 3D simulations (bottom) with same δB rms0 /B = 1but different initial magnetization σ ∈ { , , , } . Weemployed c ∆ t/l = 1 .
875 for all measurements of the energydiffusion coefficient D γ . Note that here D γ is in units of c/l .A linear fit is shown with a dashed black line. eration timescale t acc = (cid:12)(cid:12)(cid:12)(cid:12) γ d (cid:104) γ (cid:105) dt (cid:12)(cid:12)(cid:12)(cid:12) − ∼ σ lc . (40)In our simulations, the acceleration timescale increasesin time since σ decreases in time as a combined effectof the decaying turbulent fluctuations δB rms ( t ) and theincrease of the enthalpy per particle mc w ( t ).7.2. Injection and turbulence acceleration
As discussed in Sections 4 and 5, a large fraction ofparticles is preaccelerated by magnetic reconnection be-fore being accelerated by scattering off the turbulentfluctuations. This two-stages acceleration process isshown in Fig. 28 for both 2D and 3D simulations. Here,each colored curve represents the average Lorentz factorof particles having the same injection time t inj (within∆ t inj = 0 . c/l for 2D and ∆ t inj = 0 . c/l for 3D). The9 ct/l 〈 γ 〉 turbulence accelerationinjection phase ct/l 〈 γ 〉 turbulence accelerationinjection phase Figure 28.
Evolution of the mean Lorentz factor of differentgenerations of particles undergoing injection at early times( ct inj /l (cid:46)
2) for 2D turbulence (top) and 3D turbulence (bot-tom). Both simulations have σ = 10 and δB rms0 /B = 1.The initial energy gain, due to the reconnection electric field,can be modeled as in Eq. (41) with β R = 0 .
05 (dashed lines),while the subsequent evolution, governed by stochastic in-teractions with the turbulent fluctuations, follows Eq. (42)(dot-dashed line). linear growth from (cid:104) γ (cid:105) ∼ (cid:104) γ (cid:105) ∼
30 (i.e., the in-jection phase) is powered by field-aligned electric fields,whose magnitude is | E (cid:107) | (cid:39) β R δB rms , via d (cid:104) γ (cid:105) dt = emc β R δB rms . (41)The dashed black lines in Fig. 28 show Eq. (41) tak-ing a reconnection rate β R (cid:39) .
05, as appropriate forrelativistic reconnection with guide field comparable tothe alternating fields (Werner & Uzdensky 2017). Af-ter this first acceleration phase, stochastic accelerationtakes place, and, as discussed above, we can estimate d (cid:104) γ (cid:105) dt = 4 κ stoc σ (cid:16) cl (cid:17) γ . (42) with κ stoc ∼ .
03 from the 2D simulations and κ stoc ∼ . σ .A final remark concerns the acceleration timescalesassociated with magnetic reconnection and turbulencefluctuations. Fast magnetic reconnection leads to theacceleration timescale t acc = β − R ( ρ L /c ), where ρ L isthe particle Larmor radius. On the other hand, we haveseen that stochastic acceleration by turbulent fluctua-tion yields t acc = (3 /σ )( l/c ). Therefore, for the hypo-thetical case in which reconnection could drive particlesup to the highest energies ( ρ L ∼ l ), the accelerationtimescale of fast magnetic reconnection could actuallybe longer than the one associated with the turbulencefluctuations for σ (cid:38)
1. Indeed, in this magnetically-dominated regime, turbulence provides an exceptionallyfast acceleration mechanism that can potentially explainthe most extreme astrophysical accelerators. SUMMARYIn this article, we have presented the results of a se-ries of first-principles kinetic PIC simulations of decay-ing turbulence in magnetically-dominated plasmas, withthe goal of understanding how plasma turbulence, andits interplay with magnetic reconnection, can acceler-ate charged particles. We considered a pair (electron-positron) plasma, which is relevant for various astro-physical systems, such as jets from supermassive blackholes, pulsar and magnetar magnetospheres, winds, andwind nebulae. In this regime, our computational domain(2460 cells in 3D; from 16400 to 65600 cells in 2D)is large enough to capture the turbulence cascade fromlarge (MHD) scales to small (kinetic) scales.In the following, we itemize the main points of thispaper.1. The generation of a large population of nonthermalparticles is a self-consistent by-product of both 2D and3D magnetically-dominated turbulence. In particular,the late time particle energy spectrum displays a power-law high-energy range whose slope p , high-energy cutoff γ c , and fraction of particles in the power-law tail ζ nt aremarkedly similar in 2D and 3D, even though the timedevelopment of the particle energy spectrum is different.2. The power-law slope decreases (i.e., becomesharder) with increasing initial values of magnetizationand fractional strength of the turbulence fluctuations,0with slopes that can be as hard as p (cid:46)
2. In con-trast, the initial plasma temperature does not affect thepower-law slope, but only yields an overall energy shiftto larger energies for higher initial plasma temperatures.For power-law energy tails with p > π/k I d e , the larger the high-energy cutoff,which can extend up to γ c ∼ ( e/mc ) (cid:112) (cid:104) B (cid:105) π/k I , ifturbulence survives long enough to allow the particlesto reach this upper limit. The fact that the power-lawstarts close to the peak of the distribution yields a largefraction of particles in the nonthermal tail. For thephysical parameters explored in this work, we obtaina number fraction ζ nt ∼
15% - 31% of particles in thenonthermal tail.3. The majority of particles are injected into accel-eration at regions of high electric curent density. Morespecifically, a large fraction of particles is extracted fromthe thermal pool and injected into the acceleration pro-cess by reconnecting current sheets. These reconnectingcurrent sheets are strongly unstable to the formationof plasmoids, which allows fast magnetic reconnectionto occur. We observe the development of plasmoids incurrent sheets formed as a self-consistent result of mag-netized turbulence, both in 2D and in 3D. In 3D, theyappear as a chain of flux ropes elongated in the directionof the mean magnetic field.4. Reconnecting current sheets are efficient in inject-ing particles (i.e., they promote a large fraction of parti-cles in the nonthermal tail) in spite of their small fillingfraction, as they can process a large fraction of particleswithin the sheet lifetime. The efficiency remains highalso when increasing system size, as we have shown thatthe plasmoid instability (whose properties are obtainedfrom a tearing mode dispersion relation generalized forrelativistically hot plasmas) ensures the triggering of fastmagnetic reconnection within the lifetime of the large-scale current sheets, which are the ones that dominatethe particle injection census. As a consequence, mag-netic reconnection can process a large volume of plasmain few large (outer-scale) eddy turnover times (a volume V R ∼ β R L in one outer-scale eddy turnover time).5. Particle acceleration at reconnecting current sheetscan propel particles up to a typical Lorentz factor gain∆ γ inj = κσγ th , after which the acceleration is contin-ued by means of stochastic scattering off turbulent fluc-tuations. It is the stochastic acceleration process thatallows particles to reach the highest energies, up to aLarmor radius roughly equal to the energy-containingscale of the turbulence. The work done by the electricfield parallel to the magnetic field (which is expectedat reconnecting current sheets), W (cid:107) , is responsible for most of the early particle energy gain (injection). On theother hand, the second acceleration phase is powered byperpendicular electric fields. For high-energy particles,i.e., such that ∆ γ (cid:29) κσγ th , we find W ⊥ (cid:29) W (cid:107) , i.e., thework done by perpendicular electric fields dominates theoverall energy gain.6. An additional confirmation of the fact that theparallel electric field controls the injection physics butnot the subsequent acceleration process comes from anumerical simulation with extra (test) particles that donot feel parallel electric fields. This shows that the in-jection fraction is strongly suppressed. In fact, only asmall fraction of these test particles participate in theacceleration process ( ζ nt decreases by almost two ordersof magnitude). On the other hand, for those test parti-cles that can participate in the acceleration process, thepower-law slope p is very similar to that of the regularparticles. This indicates that acceleration by the per-pendicular electric field controls the slope of the power-law high-energy tail.7. The fact that different energization mechanismsdominate at different energy ranges affects the parti-cle pitch-angle distribution, f (cos α, γ ). We find thatthe pitch-angle distribution develops distinguishing fea-tures at low , intermediate , and high values of γ . Thesevalues depend on the initial mean Lorentz factor andmagnetization. For γ ∼ ( σ / γ th , particles velocitiesare strongly aligned/antialigned with the local magneticfield B , while at γ (cid:29) σ / γ th , particles velocitiesare mostly perpendicular to B . At intermediate energiessuch that γ ∼ σ / γ th , particles follow a distributionwhich has minima for both parallel and perpendiculardirections (i.e., at cos α = ± , α = ±
1, asthe low-energy population controls the number census.8. The different energization mechanisms are also re-sponsible for producing a gyrotropic four-velocity distri-bution with distinct features in the direction pertainingto the mean magnetic field B = B ˆ z . Specifically, thedomain-averaged four-velocity distribution is elongatedin the γβ z direction at low particle energies, due to the v · E (cid:107) energization, while it becomes elongated in thedirection perpendicular to the mean field at high parti-cle energies, due to the v · E ⊥ energization. At inter-mediate energies the distribution peaks at intermediateangles (i.e., at 45 degrees from the γβ z axis).9. After the injection phase, particles exhibit a dif-fusive energy behavior in both 2D and 3D turbulence.We measured the diffusion coefficient in energy space di-rectly from our PIC simulations, showing that D γ ∝ γ D γ ∝ σ , with σ being the time-dependent magnetiza-tion. The estimated energy diffusion coefficient D γ ∼ . σ ( c/l ) γ gives an acceleration timescale that can bevery fast, t acc ∼ (3 /σ )( l/c ), comparable to that of fastmagnetic reconnection or even higher, depending on theplasma magnetization.10. The mean energy gain of particles during thefirst acceleration phase (injection) is well described bylinear acceleration by the typical reconnection electricfield. Then, the subsequent mean energy gain due tostochastic scattering off the turbulent fluctuations fol-lows from the energy diffusion coefficient D γ . In oursimulations of decaying turbulence, as the plasma mag-netization decreases due to the magnetic field annihi-lation, the stochastic acceleration timescale gets longerover time and the stochastic acceleration process even-tually saturates.The aforementioned findings have implications for ourunderstanding of the generation of nonthermal particlesin high-energy astrophysical sources. The main astro-physical implications are: (i) the power-law slopes ofthe emitting particles, which are predicted to be harderfor larger plasma magnetizations and stronger turbu-lent fluctuations, can potentially explain the hard ra-dio spectrum of the Crab Nebula (e.g. Lyutikov et al.2019); (ii) the anisotropy of the particle pitch-angle dis-tribution, for which the synchrotron spectrum of theemitting particles is expected to be different than thecommonly-assumed case of isotropic particles, has con-sequences for our understanding of emission from AGNjets (e.g. Tavecchio & Sobacchi 2019); (iii) magnetically-dominated plasma turbulence leads to particle accelera-tion on rapid timescales, which can be even shorter thanthose associated with fast magnetic reconnection and arethen capable to explain particle acceleration in the mostextreme astrophysical accelerators (e.g. Takahashi et al.2009).We acknowledge fruitful discussions with MikhailMedvedev, Jonathan Zrake, Vah´e Petrosian, MartinLemoine, Aaron Tran, Chuanfei Dong, Yi-Min Huang,and Maxim Lyutikov. This research acknowledges sup-port from DoE DE-SC0016542, NSF ACI-1657507 andNASA ATP NNX17AG21G. The simulations were per-formed on Columbia University (Habanero and Terre-moto), NASA-HEC (Pleiades), NERSC (Cori and Edi-son), TACC (Stampede2) and ORNL (Titan) resources.APPENDIXIn the magnetically-dominated regime ( σ (cid:29)
1) stud-ied here, we have verified that the presented results are -2.0 -2.0 2.02.0 z J z J y /l y /l y /l x/l x/l Figure 29.
Formation of current sheets and plasmoids(in the central part of the zoomed domain) from two 2Dsimulations where the initial plasma skin depth d e is re-solved with 3 cells (left column) and 10 cells (right column).Top, middle, and bottom panels refer to frames taken at ct/l = 1 . , . , .
2. In both cases σ = 10, δB rms0 /B = 1, L/d e = 1640, and l = L/
8. In both cases we employ 16particles per cell. converged with the adopted grid resolution and numberof particles per cell.In this study, we presented results where the initialplasma skin depth d e is resolved with 10 cells in 2Dand 3 cells in 3D. However, 3 cells per initial skin depth d e (the skin depth increases during the simulation asthe mean Lorentz factor increases as a result of mag-netic field annihilation) are already sufficient to resolvecurrent sheets and plasmoids. In Fig. 8, we show theelectric current density J z taken at three different times( ct/l = 1 . , . , .
2) from two 2D simulations where d e is resolved with 3 cells (left column) and 10 cells (rightcolumn), which produce analogous structures.2 γ - -6 -5 -4 -3 -2 -1 d N / d l n ( γ - ) ppc=256ppc=128ppc=64ppc=32ppc=16ppc=4 p=2.7 Figure 30.
Particle spectra dN/d ln( γ −
1) at late times for2D simulations with σ = 10, δB rms0 /B = 1, L/d e = 820,and l = L/
8, using different values of computational particlesper cell, from ppc=4 to ppc=256.
We have also checked for convergence with respect tocomputational particles per cell. A comparison of thelate time spectra from simulations employing differentparticles per cell (up to 256) is shown in Fig. 30, forsimulations having domain size
L/d e = 820. We can seethat the particle spectra are converged with the adoptedparticle resolution. Indeed, noise-level fluctuations areon small scales and do not affect the acceleration processin the regime investigated here.REFERENCES Abdelhamid, H. M., Lingam, M., & Mahajan, S. M. 2016,ApJ, 829, 87, doi: 10.3847/0004-637X/829/2/87Ambrosiano, J., Matthaeus, W. H., Goldstein, M. L., &Plante, D. 1988, J. Geophys. Res., 93, 14383,doi: 10.1029/JA093iA12p14383Ara, G., Basu, B., Coppi, B., et al. 1978, Annals of Physics,112, 443, doi: 10.1016/S0003-4916(78)80007-4Armstrong, J. W., Rickett, B. J., & Spangler, S. R. 1995,ApJ, 443, 209, doi: 10.1086/175515Arzner, K., Knaepen, B., Carati, D., Denewet, N., &Vlahos, L. 2006, ApJ, 637, 322, doi: 10.1086/498341Arzner, K., & Vlahos, L. 2004, ApJL, 605, L69,doi: 10.1086/392506Baalrud, S. D., Bhattacharjee, A., & Daughton, W. 2018,Physics of Plasmas, 25, 022115, doi: 10.1063/1.5020777Balbus, S. A., & Hawley, J. F. 1998, Reviews of ModernPhysics, 70, 1, doi: 10.1103/RevModPhys.70.1Begelman, M. C., Blandford, R. D., & Rees, M. J. 1984,Reviews of Modern Physics, 56, 255,doi: 10.1103/RevModPhys.56.255Beresnyak, A., & Li, H. 2016, ApJ, 819, 90,doi: 10.3847/0004-637X/819/2/90Bessho, N., & Bhattacharjee, A. 2012, ApJ, 750, 129,doi: 10.1088/0004-637X/750/2/129Bhattacharjee, A., Huang, Y.-M., Yang, H., & Rogers, B.2009, Physics of Plasmas, 16, 112102,doi: 10.1063/1.3264103Birdsall, C. K., & Langdon, B. 1985, Plasma Physics viaComputer Simulation (McGraw-Hill) Biskamp, D. 2000, Magnetic Reconnection in Plasmas,Vol. 3 (Cambridge University Press)—. 2003, Magnetohydrodynamic Turbulence (CambridgeUniversity Press)Biskamp, D., & Schwarz, E. 2001, Physics of Plasmas, 8,3282, doi: 10.1063/1.1377611Biskamp, D., & Welter, H. 1989, Physics of Fluids B, 1,1964, doi: 10.1063/1.859060Blandford, R., & Eichler, D. 1987, PhR, 154, 1,doi: 10.1016/0370-1573(87)90134-7Boldyrev, S. 2006, PhRvL, 96, 115002,doi: 10.1103/PhysRevLett.96.115002Boldyrev, S., & Loureiro, N. F. 2017, ApJ, 844, 125,doi: 10.3847/1538-4357/aa7d02Brandenburg, A., Kahniashvili, T., & Tevzadze, A. e. G.2015, PhRvL, 114, 075001,doi: 10.1103/PhysRevLett.114.075001Brandenburg, A., & Subramanian, K. 2005, PhR, 417, 1,doi: 10.1016/j.physrep.2005.06.005B¨uhler, R., & Blandford, R. 2014, Reports on Progress inPhysics, 77, 066901, doi: 10.1088/0034-4885/77/6/066901Buneman, O. 1993, in “Computer Space Plasma Physics”,Terra Scientific, Tokyo, 67Carbone, V., Veltri, P., & Mangeney, A. 1990, Physics ofFluids A, 2, 1487, doi: 10.1063/1.857598Cerri, S. S., & Califano, F. 2017, New Journal of Physics,19, 025007, doi: 10.1088/1367-2630/aa5c4aCerri, S. S., Servidio, S., & Califano, F. 2017, ApJ, 846,L18, doi: 10.3847/2041-8213/aa87b0 Cerutti, B., Werner, G. R., Uzdensky, D. A., & Begelman,M. C. 2012, The Astrophysical Journal, 754, L33,doi: 10.1088/2041-8205/754/2/L33Chandran, B. D. G. 2000, PhRvL, 85, 4656,doi: 10.1103/PhysRevLett.85.4656Cho, J. 2005, ApJ, 621, 324, doi: 10.1086/427493Cho, J., & Lazarian, A. 2006, ApJ, 638, 811,doi: 10.1086/498967—. 2014, ApJ, 780, 30, doi: 10.1088/0004-637X/780/1/30Comisso, L., & Asenjo, F. A. 2014, PhRvL, 113, 045001,doi: 10.1103/PhysRevLett.113.045001Comisso, L., & Bhattacharjee, A. 2016, Journal of PlasmaPhysics, 82, 595820601, doi: 10.1017/S002237781600101XComisso, L., & Grasso, D. 2016, Physics of Plasmas, 23,032111, doi: 10.1063/1.4942940Comisso, L., Grasso, D., & Waelbroeck, F. L. 2015, Physicsof Plasmas, 22, 042109, doi: 10.1063/1.4918331Comisso, L., Huang, Y. M., Lingam, M., Hirvijoki, E., &Bhattacharjee, A. 2018, ApJ, 854, 103,doi: 10.3847/1538-4357/aaac83Comisso, L., Lingam, M., Huang, Y.-M., & Bhattacharjee,A. 2016, Physics of Plasmas, 23, 100702,doi: 10.1063/1.4964481Comisso, L., Lingam, M., Huang, Y. M., & Bhattacharjee,A. 2017, ApJ, 850, 142, doi: 10.3847/1538-4357/aa9789Comisso, L., & Sironi, L. 2018, PhRvL, 121, 255101,doi: 10.1103/PhysRevLett.121.255101Coppi, B., Galvao, R., Pellat, R., Rosenbluth, M., &Rutherford, P. 1976, Fizika Plazmy, 2, 961Cranmer, S. R., van Ballegooijen, A. A., & Edgar, R. J.2007, ApJS, 171, 520, doi: 10.1086/518001Dalena, S., Rappazzo, A. F., Dmitruk, P., Greco, A., &Matthaeus, W. H. 2014, ApJ, 783, 143,doi: 10.1088/0004-637X/783/2/143Daughton, W., & Karimabadi, H. 2007, Physics of Plasmas,14, 072303, doi: 10.1063/1.2749494Daughton, W., Roytershteyn, V., Albright, B. J., et al.2009, PhRvL, 103, 065004,doi: 10.1103/PhysRevLett.103.065004Daughton, W., Roytershteyn, V., Karimabadi, H., et al.2011, Nature Physics, 7, 539, doi: 10.1038/nphys1965Daughton, W., Scudder, J., & Karimabadi, H. 2006,Physics of Plasmas, 13, 072101, doi: 10.1063/1.2218817Dmitruk, P., & Matthaeus, W. H. 2006, Physics ofPlasmas, 13, 042307, doi: 10.1063/1.2192757Dmitruk, P., Matthaeus, W. H., & Seenu, N. 2004, ApJ,617, 667, doi: 10.1086/425301Dong, C., Wang, L., Huang, Y.-M., Comisso, L., &Bhattacharjee, A. 2018, PhRvL, 121, 165101,doi: 10.1103/PhysRevLett.121.165101 Ebrahimi, F. 2017, Physics of Plasmas, 24, 056119,doi: 10.1063/1.4983631Franci, L., Cerri, S. S., Califano, F., et al. 2017, ApJL, 850,L16, doi: 10.3847/2041-8213/aa93fbFraschetti, F., & Melia, F. 2008, MNRAS, 391, 1100,doi: 10.1111/j.1365-2966.2008.13987.xFurth, H. P., Killeen, J., & Rosenbluth, M. N. 1963,Physics of Fluids, 6, 459, doi: 10.1063/1.1706761Goldreich, P., & Sridhar, S. 1995, ApJ, 438, 763,doi: 10.1086/175121Gonz´alez, C. A., Dmitruk, P., Mininni, P. D., & Matthaeus,W. H. 2017, ApJ, 850, 19, doi: 10.3847/1538-4357/aa8c02Guo, F., Li, H., Daughton, W., & Liu, Y.-H. 2014, PhysicalReview Letters, 113, 155005,doi: 10.1103/PhysRevLett.113.155005Guo, F., Liu, Y.-H., Daughton, W., & Li, H. 2015, ApJ,806, 167, doi: 10.1088/0004-637X/806/2/167Haggerty, C. C., Parashar, T. N., Matthaeus, W. H., et al.2017, Physics of Plasmas, 24, 102308,doi: 10.1063/1.5001722Harris, E. G. 1962, Il Nuovo Cimento, 23, 115,doi: 10.1007/BF02733547Huang, Y.-M., & Bhattacharjee, A. 2010, Physics ofPlasmas, 17, 062104, doi: 10.1063/1.3420208—. 2016, ApJ, 818, 20, doi: 10.3847/0004-637X/818/1/20Huang, Y.-M., Comisso, L., & Bhattacharjee, A. 2017, ApJ,849, 75, doi: 10.3847/1538-4357/aa906dInoue, T., Asano, K., & Ioka, K. 2011, ApJ, 734, 77,doi: 10.1088/0004-637X/734/2/77Iroshnikov, P. S. 1963, AZh, 40, 742Isliker, H., Vlahos, L., & Constantinescu, D. 2017, PhRvL,119, 045101, doi: 10.1103/PhysRevLett.119.045101Ji, H., & Daughton, W. 2011, Physics of Plasmas, 18,111207, doi: 10.1063/1.3647505Kagan, D., Sironi, L., Cerutti, B., & Giannios, D. 2015,Space Science Reviews, doi: 10.1007/s11214-014-0132-9Kimura, S. S., Toma, K., Suzuki, T. K., & Inutsuka, S.-i.2016, ApJ, 822, 88, doi: 10.3847/0004-637X/822/2/88Kimura, S. S., Tomida, K., & Murase, K. 2019, MNRAS,485, 163, doi: 10.1093/mnras/stz329Koide, S. 2009, ApJ, 696, 2220,doi: 10.1088/0004-637X/696/2/2220Kowal, G., de Gouveia Dal Pino, E. M., & Lazarian, A.2012, PhRvL, 108, 241102,doi: 10.1103/PhysRevLett.108.241102Kraichnan, R. H. 1965, Physics of Fluids, 8, 1385,doi: 10.1063/1.1761412Kulsrud, R. M., & Ferrari, A. 1971, Ap&SS, 12, 302,doi: 10.1007/BF00651420 Kumar, P., & Narayan, R. 2009, MNRAS, 395, 472,doi: 10.1111/j.1365-2966.2009.14539.xKunz, M. W., Stone, J. M., & Quataert, E. 2016, PhysicalReview Letters, 117, 235101,doi: 10.1103/PhysRevLett.117.235101Lazarian, A., Vlahos, L., Kowal, G., et al. 2012, SSRv, 173,557, doi: 10.1007/s11214-012-9936-7Lemoine, M. 2019, PhRvD, 99, 083006,doi: 10.1103/PhysRevD.99.083006Lithwick, Y., & Goldreich, P. 2001, ApJ, 562, 279,doi: 10.1086/323470Liu, Y.-H., Guo, F., Daughton, W., Li, H., & Hesse, M.2015, Physical Review Letters, 114, 095002,doi: 10.1103/PhysRevLett.114.095002Liu, Y.-H., Hesse, M., Guo, F., et al. 2017, PhRvL, 118,085101, doi: 10.1103/PhysRevLett.118.085101Loureiro, N. F., & Boldyrev, S. 2017, PhRvL, 118, 245101,doi: 10.1103/PhysRevLett.118.245101Lynn, J. W., Quataert, E., Chandran, B. D. G., & Parrish,I. J. 2014, ApJ, 791, 71,doi: 10.1088/0004-637X/791/1/71Lyutikov, M., Sironi, L., Komissarov, S. S., & Porth, O.2017, Journal of Plasma Physics, 83, 635830602,doi: 10.1017/S002237781700071XLyutikov, M., Temim, T., Komissarov, S., et al. 2019,MNRAS, 2051, doi: 10.1093/mnras/stz2023MacDonald, N. R., & Marscher, A. P. 2018, ApJ, 862, 58,doi: 10.3847/1538-4357/aacc62Mallet, A., Schekochihin, A. A., & Chandran, B. D. G.2017, MNRAS, 468, 4862, doi: 10.1093/mnras/stx670Marscher, A. P., Jorstad, S. G., D’Arcangelo, F. D., et al.2008, Nature, 452, 966, doi: 10.1038/nature06895Matsumoto, Y., Amano, T., Kato, T. N., & Hoshino, M.2015, Science, 347, 974, doi: 10.1126/science.1260168Matthaeus, W. H., & Lamkin, S. L. 1986, Physics of Fluids,29, 2513, doi: 10.1063/1.866004Matthaeus, W. H., Zank, G. P., Oughton, S., Mullan, D. J.,& Dmitruk, P. 1999, ApJL, 523, L93, doi: 10.1086/312259Melrose, D. B. 1980, Plasma astrohysics. Nonthermalprocesses in diffuse magnetized plasmas - Vol.1: Theemission, absorption and transfer of waves in plasmas;Vol.2: Astrophysical applicationsMicha(cid:32)lek, G., & Ostrowsky, M. 1996, Nonlinear Processesin Geophysics, 3, 66N¨attil¨a, J. 2019, arXiv e-prints, arXiv:1906.06306.https://arxiv.org/abs/1906.06306Ni, L., Germaschewski, K., Huang, Y.-M., et al. 2010,Physics of Plasmas, 17, 052109, doi: 10.1063/1.3428553O’Sullivan, S., Reville, B., & Taylor, A. M. 2009, MNRAS,400, 248, doi: 10.1111/j.1365-2966.2009.15442.x Papini, E., Franci, L., Landi, S., et al. 2019, ApJ, 870, 52,doi: 10.3847/1538-4357/aaf003Passot, T., Sulem, P. L., & Tassi, E. 2017, Journal ofPlasma Physics, 83, 715830402,doi: 10.1017/S0022377817000514Pecora, F., Servidio, S., Greco, A., et al. 2018, Journal ofPlasma Physics, 84, 725840601,doi: 10.1017/S0022377818000995Petropoulou, M., & Sironi, L. 2018, MNRAS, 481, 5687,doi: 10.1093/mnras/sty2702Petrosian, V. 2012, SSRv, 173, 535,doi: 10.1007/s11214-012-9900-6Piran, T. 2004, Reviews of Modern Physics, 76, 1143,doi: 10.1103/RevModPhys.76.1143Politano, H., Pouquet, A., & Sulem, P. L. 1995, Physics ofPlasmas, 2, 2931, doi: 10.1063/1.871473Porcelli, F. 1991, PhRvL, 66, 425,doi: 10.1103/PhysRevLett.66.425Porth, O., Komissarov, S. S., & Keppens, R. 2014,MNRAS, 438, 278, doi: 10.1093/mnras/stt2176Ramaty, R. 1979, in American Institute of PhysicsConference Series, Vol. 56, Particle AccelerationMechanisms in Astrophysics, ed. J. Arons, C. McKee, &C. Max, 135–154Retin`o, A., Sundkvist, D., Vaivads, A., et al. 2007, NaturePhysics, 3, 236, doi: 10.1038/nphys574Roy, N., Bharadwaj, S., Dutta, P., & Chengalur, J. N. 2009,MNRAS, 393, L26, doi: 10.1111/j.1745-3933.2008.00591.xSchlickeiser, R. 1989, ApJ, 336, 243, doi: 10.1086/167009Servidio, S., Matthaeus, W. H., Shay, M. A., Cassak, P. A.,& Dmitruk, P. 2009, PhRvL, 102, 115003,doi: 10.1103/PhysRevLett.102.115003Servidio, S., Matthaeus, W. H., Shay, M. A., et al. 2010,Physics of Plasmas, 17, 032315, doi: 10.1063/1.3368798Servidio, S., Valentini, F., Califano, F., & Veltri, P. 2012,PhRvL, 108, 045001,doi: 10.1103/PhysRevLett.108.045001Shi, C., Velli, M., & Tenerani, A. 2018, ApJ, 859, 83,doi: 10.3847/1538-4357/aabd83Sironi, L., Giannios, D., & Petropoulou, M. 2016, MNRAS,462, 48, doi: 10.1093/mnras/stw1620Sironi, L., & Spitkovsky, A. 2014, ApJL, 783, L21,doi: 10.1088/2041-8205/783/1/L21Skilling, J. 1975, Monthly Notices of the RoyalAstronomical Society, 172, 557,doi: 10.1093/mnras/172.3.557Spitkovsky, A. 2005, in AIP Conf. Ser., Vol. 801,Astrophysical Sources of High Energy Particles andRadiation, ed. T. Bulik, B. Rudak, & G. Madejski, 3455