FIRST-ORDER PARTICLE ACCELERATION IN MAGNETICALLY DRIVEN FLOWS

Magnetically dominated environments are fairly common in astrophysics and a sizable fraction of astrophysical objects are made of rarefied plasma. These objects are only visible because they contain a non-thermal particle population that dominates emission across most of the spectrum. The basic idea that non-thermal emission requires particle acceleration has been around for some time. One of the main elements of the explanation of how energy is being transferred to particles is the mechanism of energy dissipation. Indeed, if astrophysical fluids may be considered almost inviscid and perfectly conducting, how can large-scale energy be lost to thermal and non-thermal particles? A few mechanisms have been suggested: (a) discontinuities in the fluid motion, e.g., shocks, (b) turbulence, and (c) discontinuities in the magnetic field, i.e., current sheets. Following the ideas of Fermi (1949), who posited that collisionless particles can be energized by scattering in fluid motions, especially converging motions, the diffusive shock acceleration (DSA) mechanism has been proposed (Krymskii 1977; Bell 1978; Malkov & O'C Drury 2001) and has become rather popular for explaining non-thermal electromagnetic emission, as well as cosmic rays—energetic charged particles, detected at Earth. In the DSA, however, the acceleration rate is related to the scattering rate, with the latter being importantly bound by the so-called Bohm limit, which assumes that the scattering rate must be lower than Larmor gyration frequency. This makes acceleration progressively slower at higher energies.

On the other hand, the observations of many variable astrophysical objects suggest extremely fast acceleration timescales which are incompatible with DSA. Few blazar jets powered by supermassive black holes exhibit ∼10 minute variations in TeV emissions, e.g., Aharonian et al. (2007), Aleksić et al. (2011). Such fast time variabilities, along with other emission region constraints, have been suggested as evidence for mini-jets generated by reconnection (Giannios 2013). Recent observations of gamma-ray flares from Crab (Abdo et al. 2011) revealed that the impulsive nature of the energy release and the associated particle acceleration might need an alternative explanation as well (Clausen-Brown & Lyutikov 2012). It has also been suggested that reconnection plays a crucial role in producing high-energy emissions from gamma-ray bursts (Zhang & Yan 2011).

In the environments that are known to be magnetically dominated, e.g., the solar corona or the pulsar wind nebula, it is natural to expect the magnetic energy to be the main source. The global energetics of the powerful X-ray flares confirm this natural assumption. Reconnection and associated phenomena have been an active field of study, (see, e.g., Uzdensky (2015 for a review). It would be interesting to find out if there is a generic mechanism to transfer magnetic energy into particles and to conceptually understand the nature of recent numerical results that demonstrated that in both MHD fluid simulations (Kowal et al. 2012) and ab-initio plasma simulations of reconnection (Guo et al. 2014; Sironi & Spitkovsky 2014; Uzdensky & Rightley 2014) there is a regular acceleration of particles.

Particle acceleration is often classified into a "first-order Fermi" mechanism where particles are gaining energy regularly, e.g., by colliding with converging magnetic mirrors and a "second-order Fermi" mechanism where particles can both gain and lose energy (Fermi 1949). These two are not mutually exclusive and represent two different terms in the equation for the evolution of the distribution function: the terms describing advection and diffusion in energy space, respectively. In practical terms, the first-order mechanism usually dominates, if present, as the acceleration rate is proportional to the first order of u/v, where u is the typical velocity of the scatterer and v is the particle velocity. The outcome of the first-order acceleration can be described in terms of the rate at which particles gain energy, acceleration rate, r_acc, and the rate of particle escaping from the system, r_esc. If escape is negligible and r_acc is constant with energy, the energy of each particle grows exponentially. Also, if r_esc/r_acc does not depend on energy, the stationary solution for the particle distribution is a power law, with the power-law index determined by $-1-{r}_{\mathrm{esc}}/{r}_{\mathrm{acc}}$; see, e.g., Drury et al. (1999). Various environments, such as supernova shocks, were once thought to satisfy this condition and produce power-law distributed cosmic rays, which become consistent with observations after being modified by diffusion from the Galaxy. Acceleration within many orders of magnitude in energy was regarded as a result of a large-scale physical layout of the acceleration site, e.g., the planar shock can be thought of as a set of large-scale converging mirrors. The very same picture could also be applied to the large-scale reconnection site, where the two sides of the inflow effectively work as converging mirrors (de Gouveia dal Pino & Lazarian 2005). In this paper we deviate from this mindset in which the problem of achieving scale-free acceleration is easily resolved just because there is only a single scale—the scale of the system. Instead we will try to find regular acceleration over large energy ranges in systems that do not necessarily possess global regular structure—however, they could still be scale-free in a statistical sense, such as turbulent systems. Normally, turbulent environments are expected to be regions of second-order acceleration; see, e.g., Schlickeiser (2002), Cho & Lazarian (2006). In this paper we point to the mechanism of regular or first-order acceleration that was overlooked in the literature. This mechanism is inherently related to a certain statistical measure of energy transfer in turbulence, and therefore does not rely on a particular geometry and is very robust. As we will show below, the direction of energy transfer from a magnetic field to kinetic motions and the sign of curvature drift acceleration are inherently related, so in systems with an average positive energy transfer from magnetic energy to kinetic motions there is an average positive curvature drift acceleration, while in the opposite case, there is an average curvature drift cooling.

One of the commonly considered cases of magnetically driven flows is magnetic reconnection. Significant effort was put into understanding non-ideal plasma effects that could both cause reconnection and create non-zero parallel components of the electric field (Fu et al. 2006; Pritchett 2006). In this paper we apply our model to reconnecting layers that are large-scale in the sense that the current layer is many orders of magnitude thicker than the ion skin depth d_i. In this case the reconnecting layer will have multiple X-points and while non-ideal effects are indeed required for the individual field lines to break and reconnect, their influence is limited to fairly small scales, typically around the ion skin depth d_i. In this paper, we instead focus on larger scales and ignore non-ideal effects for the following two reasons. First, it has been argued that understanding the global energetics of large-scale reconnection, such as the amount of magnetic energy dissipated per unit time, does not necessarily require detailed knowledge on how individual field lines break and reconnect (Lazarian & Vishniac 1999; Loureiro et al. 2007; Eyink et al. 2011; Beresnyak 2013). These global energetic parameters, as we will show below, could be more important for acceleration to high energies than local non-ideal effects.  Second, in order to understand high-energy particle acceleration, one normally has to consider plasma dynamics on scales much larger than d_i, that is, on MHD scales. While modern simulations such as Sironi & Spitkovsky (2014), Guo et al. (2014), and Uzdensky & Rightley (2014), are able to reach box sizes of several hundred d_i, some theory work is needed to disentangle acceleration from MHD and the non-ideal effects.

Conductive plasmas can be described on large scales as inviscid and perfectly conducting fluids (ideal MHD). Ideal MHD equations allow for exchange between kinetic and magnetic energies. The Lorentz force density, multiplied by the fluid velocity, ${\boldsymbol{u}}\cdot [{\boldsymbol{j}}\times {\boldsymbol{B}}]/c$, is the amount of energy transferred from magnetic to kinetic energy. While macroscopic (i.e., kinetic plus magnetic) energy is expected to be conserved in the ideal MHD, it is not the case in real systems that have non-zero dissipation coefficients. This is qualitatively explained by the nonlinear turbulent cascade that brings macroscopic energy to smaller and smaller scales until it dissipates into thermal energy. One of the most important examples of this is the spontaneous reconnection where the thin current layer becomes turbulent and starts dissipating magnetic energy at a constant rate. The small scales of these turbulent flows resemble "normal" MHD turbulence, which has equipartition between magnetic and kinetic energies (Beresnyak 2013); it is also true that the kinetic and magnetic parts of the cascade each contribute around half of the total cascade rate. Therefore, if we assume that the turbulent cascade is being fed with magnetic energy, approximately half of the magnetic energy has to be transferred into kinetic energy before the equipartition cascade sets in. It follows that the ${\boldsymbol{B}}$ to ${\boldsymbol{v}}$ energy transfer must be positive on average and could be approximated by one half of the volumetric energy dissipation rate , the main parameter of turbulence. The term ${\boldsymbol{u}}\cdot [{\boldsymbol{j}}\times {\boldsymbol{B}}]/c$ is the Eulerian expression for the work done by magnetic tension upon the fluid element. This term can be rewritten as the sum of $-({\boldsymbol{u}}\cdot {\boldsymbol{\nabla }}){B}^{2}/8\pi$, advection of magnetic energy density by the fluid, and

$$\begin{eqnarray}&&{{ \mathcal T }}_{{bv}}={\boldsymbol{u}}\cdot ({\boldsymbol{B}}\cdot {\boldsymbol{\nabla }}){\boldsymbol{B}}/4\pi ,\end{eqnarray} \tag{ 1 }$$

the actual energy transfer between ${\boldsymbol{B}}$ and ${\boldsymbol{v}}$. For the purpose of future calculations we will separate the ${{ \mathcal T }}_{{bv}}$ in the following way:

$$\begin{eqnarray}{ \mathcal T } & = & \displaystyle \frac{1}{4\pi }{\boldsymbol{u}}\cdot ({\boldsymbol{B}}\cdot {\boldsymbol{\nabla }}){\boldsymbol{B}}\\ & = & \displaystyle \frac{1}{4\pi }({\boldsymbol{u}}\cdot {\boldsymbol{B}})({\boldsymbol{b}}\cdot {\boldsymbol{\nabla }})B+\displaystyle \frac{B}{4\pi }{\boldsymbol{u}}\cdot ({\boldsymbol{B}}\cdot {\boldsymbol{\nabla }}){\boldsymbol{b}}={ \mathcal X }+{ \mathcal D },\end{eqnarray} \tag{ 2 }$$

where we have designated ${\boldsymbol{b}}$ = ${\boldsymbol{B}}$/B, a unit magnetic vector. The term ${ \mathcal X }$ contains cross helicity density ${\boldsymbol{u}}\cdot {\boldsymbol{B}}$. We argue that in those systems where $\langle {\boldsymbol{u}}\cdot {\boldsymbol{B}}\rangle =0$, which include many physically relevant cases, such as spontaneous reconnection, the whole ${ \mathcal X }$ term could average out (see also Figure 2). The second term ${ \mathcal D }$ contains magnetic field curvature $({\boldsymbol{B}}\cdot {\boldsymbol{\nabla }}){\boldsymbol{b}}$ and will be important for subsequent calculations of curvature drift.

To explore the implications of magnetic energy transfer in non-thermal particle acceleration, it is instructive to consider the particles' motions in a slowly varying electric and magnetic field, which can be described by the so-called drift approximation. The leading drift terms are known as electric, gradient, and curvature drifts. While electric drift, proportional to $[{\boldsymbol{E}}\times {\boldsymbol{B}}]$, cannot produce acceleration, the other two drifts can. For example, imagine the configuration of the reconnection with the inflow (Figure 1, top panel). The gradient drift $\sim [{\boldsymbol{B}}\times {\boldsymbol{\nabla }}{\boldsymbol{B}}]$ is along $-z$, and so is the electric field in the ideal case $E=-[{\boldsymbol{u}}\times {\boldsymbol{B}}]/c$. Their product will be positive and will result in acceleration that is due to particles being compressed by the converging inflow. This mechanism does not work in the initial, most energetic stage of spontaneous reconnection, which has negligible inflow (Figure 1, middle panel; Beresnyak 2013). It is this initial stage that has the highest volumetric dissipation rate, however. Figure 1 illustrates why curvature drift acceleration is important in this configuration. It also turns out that in any magnetically driven turbulent environment, such as that depicted in the bottom panel of Figure 1, curvature drift acceleration will accelerate particles on average.

Zoom In Zoom Out Reset image size

Let us look carefully at the term that is responsible for acceleration by curvature drift,

$$\begin{eqnarray}&&\displaystyle \frac{d{ \mathcal E }}{{dt}}=-2\displaystyle \frac{{{ \mathcal E }}_{\parallel }}{B}[{\boldsymbol{u}}\times {\boldsymbol{B}}]\cdot [{\boldsymbol{b}}\times ({\boldsymbol{b}}\cdot {\boldsymbol{\nabla }}){\boldsymbol{b}}]\end{eqnarray} \tag{ 3 }$$

(see, e.g., Sivukhin 1965), where ${{ \mathcal E }}_{\parallel }={v}_{\parallel }{p}_{\parallel }/2=\gamma {{mv}}_{\parallel }^{2}/2$ is a particle's parallel kinetic energy. With some manipulation, this expression could be equivalently transformed to $2({{ \mathcal E }}_{\parallel }/B){\boldsymbol{u}}\cdot ({\boldsymbol{B}}\cdot {\boldsymbol{\nabla }}){\boldsymbol{b}}$. It now becomes clear that this term is related to the transfer rate between magnetic and kinetic energies; particularly, it is a fraction of ${ \mathcal D }$:

$$\begin{eqnarray}&&\displaystyle \frac{d{ \mathcal E }}{{dt}}={{ \mathcal E }}_{\parallel }\displaystyle \frac{8\pi }{{B}^{2}}{ \mathcal D }.\end{eqnarray} \tag{ 4 }$$

The physical meaning of this equation is that, given efficient particle scattering, so that ${{ \mathcal E }}_{\parallel }={ \mathcal E }/2$, the acceleration rate is determined by the half of the local energy transfer rate $8\pi { \mathcal D }/{B}^{2}$, not including the ${ \mathcal X }$ term.

We can test the general ideas outlined above in two physical cases that feature turbulent energy transfer from magnetic to kinetic energies. Using spontaneous reconnection and the decaying MHD turbulence simulated numerically, we can directly calculate the discussed terms and compare them. The spontaneous reconnection numerical experiment was started with a thin planar current sheet and small perturbations in ${\boldsymbol{u}}$ and ${\boldsymbol{B}}$ and was described in detail in Beresnyak (2013), while the decaying MHD turbulence was similar to our previous incompressible driven simulations in Beresnyak & Lazarian (2009), except that there was no driving and the initial conditions were set as a random magnetic field with wavenumbers 1 < k < 5 and zero velocity. Both simulations developed magnetically driven flows, from which we calculated the average ${ \mathcal T }$, ${ \mathcal D }$, and ${ \mathcal X }$ terms and presented them in Figure 2.

Zoom In Zoom Out Reset image size

The spontaneous reconnection case had a fairly stable reconnection rate; this also corresponded to the approximately constant ${ \mathcal D }$ integrated over the volume. The ${ \mathcal X }$ term did not seem to be sign-definite and contributed relatively little. Gradient drift acceleration was also negligible, possibly due to the absence of global compression. Keeping in mind that all the energy had to come from magnetic energy, and given that the dissipation rate was approximately constant, it was no surprise that the average ${ \mathcal D }$ term evolved relatively little. It should be noted, however, that in the spontaneous reconnection experiment the width of the reconnection region was growing approximately linear with time (Beresnyak 2013), so the ${ \mathcal D }$ term magnitude pertaining to the reconnection region itself was much higher than that of an averaged ${ \mathcal D }$ over the total volume. Given the reconnection layer thickness $l(t)$, the local ${ \mathcal D }$ could be estimated as $(1/2){v}_{r}({B}^{2}/8\pi )/l(t)$, where v_r is a reconnection rate, which was ${v}_{r}\approx 0.015{v}_{A}$ in Beresnyak (2013), and the 1/2 comes from only half of the magnetic energy being transferred into kinetic energy before physically dissipating on a very small dissipation scale $\eta \ll l(t)$. This would correspond to an acceleration rate of $(1/4){v}_{r}/l(t)$ and can be very high, because the $l(t)$ could be as small as the Sweet-Parker current sheet width or the ion skin depth.

The decaying MHD turbulence experienced two regimes: (1) the initial oscillation when an excessive amount of magnetic energy was converted into kinetic energy and the bounce back and partial inverse conversion afterward; (2) the self-similar decay stage of MHD turbulence. The first stage had the strongest conversion term, which was dominated by ${ \mathcal D }$. All terms integrated over time were mostly accumulated within the first 1–2 Alfvénic times.

It is known from turbulence theory that the energy transfer rate ${ \mathcal T }$ can be demonstrated to be local in scale, under relatively weak assumptions (Aluie & Eyink 2010; Beresnyak 2012). The scale-locality means that each scale contributes to the transfer independently. We also know empirically that most of the transfer between magnetic and kinetic energies happens on relatively large scales that are comparable to the outer scale of the system, while below the outer scale there is little average transfer due to an approximate equipartition between kinetic and magnetic energies. For example, the reconnecting turbulent current layer has most of its ${ \mathcal T }$ transfer within a factor of a few of the scale of the layer thickness, while in the decaying turbulence problem it is within a factor of a few from the outer scale of the system. Let us designate ${{ \mathcal T }}_{l}$ as a transfer calculated from quantities that were filtered by a low-pass Gaussian filter with a cutoff wavenumber of 1/l. Keeping locality in mind, we will conclude that only scales larger than l will contribute to ${{ \mathcal T }}_{l}$. This means that ${{ \mathcal T }}_{l}$ will be constant for small l and will start decreasing when l approaches the outer scale of the problem L. In general, we cannot deduce the same for ${{ \mathcal D }}_{l}$ and ${{ \mathcal X }}_{l}$, but keeping in mind that ${{ \mathcal X }}_{l}$ contributes relatively little in the two cases that we considered (${{ \mathcal D }}_{l}\approx {{ \mathcal T }}_{l}$ in those cases). Figure 2 demonstrates this behavior in the bottom panel, where energy transfer is operating between wavenumbers 2 and 20 in the reconnection case (1/20 is approximately the layer width at this point), and the ${{ \mathcal D }}_{l}$ mostly changes within the same range of scales. The decaying turbulence case has the transfer more localized around the outer scales.

In terms of drift, the particle with Larmor radius r_L will "feel" magnetic and electric fields on scales larger than r_L, while the scales smaller or equal to r_L will contribute to particle scattering. It follows that the "effective" ${ \mathcal D }$ will be ${{ \mathcal D }}_{{r}_{L}}$, an implicit function of energy. Combining this with the result obtained above that ${{ \mathcal D }}_{l}$ goes to a constant for small l we conclude that the acceleration rate will also go to a constant for particles with an r_L smaller than the system size. A similar, more hand-waving argument, is that the term ${{ \mathcal D }}_{l}$ could be roughly approximated as ${B}_{l}^{2}{v}_{l}/l$, resembling a turbulent energy transfer rate, which is scale-independent. Interestingly, starting with scale-independent energy transfer in turbulence, we arrived at the energy-independent acceleration rates. Given the generality of the arguments presented above it is not surprising that energy-independent rates were indeed observed in simulations (Guo et al. 2014).

The development of the thin current sheet instability results in turbulence and reconnection in the sense of dissipated magnetic energy. This process will come through two distinct regimes, the regime without significant outflow for times smaller than $L/{v}_{A}$ (Beresnyak 2013) and the stationary reconnection with an outflow for larger times (Lazarian & Vishniac 1999). Let us consider the first regime, which has a larger dissipation rate per unit volume, because the current layer thickness $l(t)={v}_{r}t$ is relatively small. We use v_r for the reconnection rate and t as the time since the beginning of spontaneous reconnection. Let us assume that the current layer thickness is limited from below by the ion skin depth. We will have an acceleration rate of $1/(4t)$ for all times larger than d_i/v_r but smaller than L/v_A. The solution for energy, therefore, will be ${ \mathcal E }={{ \mathcal E }}_{0}{({{tv}}_{r}/{d}_{i})}^{1/4}$, where ${{ \mathcal E }}_{0}$ is the initial energy, e.g., the thermal energy. The particle's energy will be ${{ \mathcal E }}_{\mathrm{max}}=0.35{{ \mathcal E }}_{0}{(L/{d}_{i})}^{1/4}$ given the reconnection rate ${v}_{r}=0.015{v}_{A}$. This, however, is only the maximum energy of accelerated particles, as only a tiny fraction of particles were contained in the original thin current sheet and started off as being accelerated from initial time ${d}_{i}/{v}_{r}$ (Guo et al. 2014). With proper stochastic treatment and assuming that the escape rate r_esc is negligible compared with acceleration rate in this no-outflow problem, we expect to see the power-law particle spectrum with the −1 slope and a cutoff at ${{ \mathcal E }}_{\mathrm{max}}.$

The subsequent development of an outflow will do three things: first, it will enable the inflow and therefore the extra acceleration term associated with gradient drift or converging magnetic mirrors (de Gouveia dal Pino & Lazarian 2005); second, it will stabilize the acceleration rate for the curvature drift acceleration at ${v}_{A}/4L$, as the current layer is no longer expanding; and third, it will enable particle escape through the outflow. In this regime the spectrum will extend from ${{ \mathcal E }}_{\mathrm{max}}$ to higher energies, up to Larmor radii of the large scale of the current sheet L. The spectral slope of this extension will be determined by $-1-{r}_{\mathrm{esc}}/{r}_{\mathrm{acc}}$, where r_acc should account for all acceleration and cooling mechanisms, such as gradient drift acceleration and outflow cooling. Detailed analysis of this stage will be the subject of a future work.

The electron spectra observed in solar X-ray flares are fitted with the thermal component with a temperature of several keV and the steep power-law component. This is consistent with our picture, as the rather shallow −1 slope, containing most of its energy near ${{ \mathcal E }}_{\mathrm{max}}$, is likely to thermalize. Also, the outflow phase will extend this distribution as a power law to higher energies.

We demonstrated that magnetic configurations that relax to the lower states of magnetic energy will also regularly accelerate particles, on timescales which are, typically, Alfvénic, but can be much shorter, e.g., at the beginning of a spontaneous reconnection. This mechanism of acceleration of collisionless non-thermal particles by an MHD electric field should not be confused with the acceleration of the bulk of the plasma by magnetic tension. Indeed, for particles with low energies, the drift terms could be neglected, i.e., the right side of Equation (3) will be trivially zero. In the bulk fluid acceleration the energy gained by each particle does not depend on its initial energy, while for the drift acceleration scenario it is proportional to the particle energy. Another way to understand the difference between plasma heating and our acceleration mechanism is to consider a turbulent dynamo where the kinetic energy of fluid motions is partially converted into magnetic energy and partially dissipated by the turbulent cascade. In this case plasma will be heated due to the turbulent cascade, while we expect energetic particles to be cooled or decelerated by the curvature drift, since the ${ \mathcal T }$ term will be negative. This means that the relation between energy dissipation and acceleration of particles is not trivial.

Some recent observations, e.g., Aharonian et al. 2007, Abdo et al. 2011, and Aleksić et al. 2011, suggested that high-energy emission variability could be as fast as variability at lower energies, which is at odds with DSA, which predicts an acceleration timescale proportional to the diffusion coefficient, which is typically a positive power of energy. This has been pointed out as a motivation for the reconnection scenario (Zhang & Yan 2011; Clausen-Brown & Lyutikov 2012; Giannios 2013).

Particle acceleration during reconnection is a topic under intense study, but the mechanism discussed in this paper is distinctly different from the direct acceleration by the reconnecting electric field at the X-line. In fact, we completely ignore non-ideal effects that produce local ${{\boldsymbol{E}}}_{\parallel }$. Also, our mechanism is not tied to a special X-point, but is instead volumetric. An interesting first-order acceleration mechanism in ideal MHD turbulence related to imbalanced turbulence and convergent field lines has recently been suggested by Schlickeiser (2009).

Regular acceleration due to converging field lines was suggested in de Gouveia dal Pino & Lazarian (2005); later it was pointed out by Drury (2012) that an outflow cooling should also be included. In this paper we do not rely on a simple transport equation, such as Parker's; therefore we relax the above approach's requirement that particles need to be almost isotropic. The acceleration in turbulent reconnection has been further numerically studied in Kowal et al. (2012). Plasma PIC simulations have also been increasingly used to understand particle acceleration. The emphasis was mostly on the non-ideal effects near X-line regions and interaction with magnetic islands (Zenitani & Hoshino 2001; Drake et al. 2006; Fu et al. 2006; Pritchett 2006; Huang et al. 2010; Oka et al. 2010; Dahlin et al. 2014). The change of energy due to curvature drift in a single collision of a particle with a magnetic island was estimated analytically in Drake et al. (2006). An important question that was left out in that study was whether this term will result in acceleration or cooling, on average. Without understanding this, it was not clear whether this process results in acceleration or deceleration or the second-order effect. In this paper we unambiguously decide this by relating the answer to a certain well-known statistical measure in turbulence–the direction of transfer between magnetic and kinetic energies. We also showed that the curvature drift acceleration does not require particles to be trapped in contracting magnetic islands, so their energy is not limited by this requirement, meaning that the energy cutoff is not related to the island size; instead it is related to the system size (see also Uzdensky & Rightley 2014).

PIC simulations are limited in the range of scales and energies they cover. Recent simulations in Sironi & Spitkovsky (2014) and Guo et al. (2014) demonstrated acceleration up to 100 MeV in electron energies, which is below the maximum energy in most astrophysical sources. Theory, therefore, is necessary to supplement conjectures based on the observed PIC distribution tails, potentially explaining the underlying physical mechanism and making predictions for astrophysical systems that often feature gigantic scale separation between plasma scales and the size of the system. The feedback from simulations, nevertheless, was very useful, particularly the recent simulations (Guo et al. 2014) that reached MHD scales and confirmed the prediction that the curvature drift acceleration will dominate compared to the non-ideal electric field acceleration.

We are grateful to Fan Guo for sharing his results during a preliminary stage. The work is supported by the LANL/LDRD program and the DoE/Office of Fusion Energy Sciences through CMSO. Computing resources at LANL were provided through the Institutional Computing Program. We also acknowledge support from XSEDE computational grant TG-AST110057.