Particle acceleration in MHD turbulence

MHD Turbulence describes dynamics of astrophysical plasmas on large scale. It is characterized by energy transfer between different scales and the exchange of energy with nonthermal population – typically cosmic rays. Recent progress in theory regarding almost all basic regimes of turbulence – from the dynamo to the decaying case and the asymptotic scaling laws, allowed us to proceed with more observationally motivated questions. One of them is why almost all strongly magnetized environments are indeed observable, e.g. why such environments are infused with high-energy particles, their distributions stretching to energies orders of magnitude higher than thermal. It turns out that there are generic mechanisms of acceleration in turbulence, both second and first order in v/c, the latter was somehow unnoticed. These generic mechanisms may explain the ubiquity of non-thermal tails in various magnetized astrophysical environments, e.g. solar chromosphere, pulsar magnetospheres, jets from supermassive black holes, gamma-ray bursts, etc.


Introduction
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 which dominates emission across most of the spectrum. Following the idea of [1] that collisionless particles can get energy by scattering in fluid motions, especially converging motions, the diffusive shock acceleration (DSA) mechanism has been proposed [2,3,4] and has become rather popular in 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 diffusion coefficient, with the latter being importantly bound by the so-called Bohm limit 1 which assumes that 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 suggested extremely fast acceleration timescales, incompatible with DSA. Few blazar jets powered by supermassive black holes exhibit ∼ 10 min variations in TeV emissions, e.g., [6,7]. Such fast time variabilities, along with other emission region constraints, have been suggested as evidence for mini-jets generated by reconnection [8]. Recent observations of gamma-ray flares from Crab [9] revealed that the impulsive nature of the energy release and the associated particle acceleration might need an alternative explanation as well [10]. It has also been suggested that reconnection plays a crucial role in producing high energy emissions from gamma-ray bursts [11]. In the environments which 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 energy source. The global energetics of the powerful X-ray flares confirms this assumption. Reconnection and associated phenomena has been a active field of study, see, e.g. [12] for a review.
Here we describe 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 [13] and ab-initio plasma simulations [14,15,16] of reconnection there is a regular acceleration of particles Particle acceleration is often classified into "first order Fermi" mechanism where particles are gaining energy regularly, e.g., by colliding with converging magnetic mirrors and "second order Fermi" where particles can both gain and lose energy [1]. 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 correspondingly. 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 particle 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 esc /r acc , see, e.g. [17]. Various environments, such as supernova shocks, were thought to satisfy this condition and produce power-law distributed cosmic rays, which become consistent with observations after being modified by diffusion. 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 [18,19]. Such acceleration lacks characteristic energy scale because there is only a single spacial scale -the scale of the system, which determines the maximum achievable energy.
In this paper, however, 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 statistical sense, such as turbulent systems. Normally, turbulent environments are expected to be regions of second-order acceleration, see, e.g., [20,21]. 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 magnetic field to kinetic motions and the sign of curvature drift acceleration are inherently related, so that in systems with the 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. A significant effort was put into understanding on non-ideal plasma effects that could both cause reconnection and create non-zero parallel component of the electric field [22,23]. In this paper we apply our model to reconnecting layers which are large-scale in a 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. 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 Three acceleration scenarious in magnetized environments. (upper) -reconnection with inflow, acceleration is due to gradient drift, the rate is determined by v z E z . (middle)acceleration driven by contracting field lines which drive fluid motion and at the same time cause the curvature drift of particles. The direction of drift is typically perpendicular to the fluid motion. The dominant term is associated with E x and the field curvature of B y component. The bottom panelsame mechanism as in the middle panel but in a homogeneous and isotropic setup of decaying MHD turbulence, which also has contracting field lines (all vector components contribute equally to acceleration).
as [14,15,16], 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.

MHD flows and energy transfer
Well-conductive plasmas can be described on large scales as inviscid and perfectly conducting fluid (ideal MHD). The ideal MHD equations allows for exchange between kinetic and magnetic energies. The Lorentz force density, multiplied by the fluid velocity, u · [j × 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 which 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 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 [24], it is also true that kinetic and magnetic part 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 equipartition cascade sets in. It follows that the B to 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 u · [j × 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 −(u · ∇)B 2 /8π, advection of magnetic energy density by the fluid, and T bv = u · (B · ∇)B/4π, the actual energy transfer between B and v. For the purpose of future calculations we will separate the T bv in the following way: The term X contains cross helicity density u · B. We argue that in those systems where ⟨u · B⟩ = 0, which include many physically relevant cases, such as spontaneous reconnection, the whole X term could average out (see also Fig. 2). The second term D contains magnetic field curvature (B · ∇)b and will be important for subsequent calculation of curvature drift. to v accumulated down to scale l, for which reason it asymptotes to a constant at small l or large cutoff wavenumbers. The inset in the right panel shows kinetic and magnetic spectra in case A to demonstrate the range of scales within which magnetic energy is transferred to kinetic -down to k ≈ 30 with thick arrows depicting the energy transfer.

Curvature drift
To explore the implications of magnetic energy transfer in non-thermal particle acceleration, it is instructive to consider the particles motion in slowly-varying electric and magnetic field, which can be described in the so-called drift approximation. The leading drift terms are known as electric, gradient and curvature drifts. While electric drift, proportional to [E × B], can not produce acceleration, the other two drifts can. For example, imagine the configuration of the reconnection with the inflow , Fig 1 top panel. The gradient drift ∼ [B × ∇B] is along −z, and so is the electric field in the ideal case E = −[u × B]/c. Their product will be positive and will result in acceleration which is due to particles being compressed by the converging inflow. This mechanism does not work in the initial, most energetic stages of spontaneous reconnection which has negligible inflow, Fig. 1 middle panel [24]. It is this initial stage that has the highest volumetric dissipation rate, however. Fig 1 illustrates why curvature drift acceleration is important in this configuration. It also turns out that in any magneticallydriven turbulent environment, such as depicted on the bottom panel of Fig 1, curvature drift acceleration will accelerate particles on average. Let us look carefully at the term which is responsible for acceleration by curvature drift, see, e.g., [25], where E ∥ = v ∥ p ∥ /2 = γmv 2 ∥ /2 is a particle's parallel kinetic energy. With some manipulation, this expression could be equivalently transformed as 2(E ∥ /B)u · (B · ∇)b. It now becomes clear that this term is related to the transfer rate between magnetic and kinetic energies, in particular it is a fraction of D: The quantity D depends on the coordinate and can be negative or positive, however we can average the Equation 3 in a volume sufficient to obtain average energy transfer between v and In this volume, we can estimate E ∥ ≈ E/3, this implies averaged particle distribution to be isotropic, note that it can still be locally anisotropic. The acceleration rate, as the rate of growth of the total energy will be determined by 8πD/3B 2 , not including the contribution from the X term. The above quantity is a 1/3τ tr , where τ tr is the characteristic energy transfer timescale.

A case study
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 thin planar current sheet and small perturbations in u and B and was described in detail in [24], while the decaying MHD turbulence was similar to our previous incompressible driven simulations in [26], 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 magneticallydriven flows, from which we calculated the average T , D, and X terms and presented them in Fig. 2. The spontaneous reconnection case had fairly stable reconnection rate, this also corresponded to the approximately constant D integrated over the volume. The X term didn't seem to be signdefinite 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 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 [24], so the D term magnitude, pertaining to the reconnection region itself was much higher than that of an averaged D over the total volume. Given the reconnection layer thickness l(t), the local D could be estimated as (1/2)v r (B 2 /8π)/l(t), where v r is a reconnection rate, which was v r ≈ 0.015v A in [24] and the 1/2 comes from only half of magnetic energy being transferred into kinetic energy before physically dissipating on very small dissipation scale η ≪ l(t). This would correspond to 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 excessive amount of magnetic energy was converted into kinetic energy and the bounce back and partial inverse conversion afterwards; 2) the self-similar decay stage of MHD turbulence. The first stage had the strongest conversion term which was dominated by D. All terms integrated over time were mostly accumulated within the first 1-2 Alfvénic times.

Reconnection Simulations Setup
Earlier we simulated spontaneous reconnection in resistive MHD with Lundquist numbers up to 6 × 10 4 , in which which the turbulent current layer is expanding in time [24]. The setup of these simulations on Fig 3. Additionally to this, we obtained a range of statistical measures of current layer turbulence [24]. Most interestingly: a) the power spectra of turbulent perturbations resembled decaying magnetic turbulence with the magnetic contribution dominating over the kinetic contribution on large scales and with an approximate energy equipartition on small scales; b) the energy spectral slope was around −1.5 ÷ −1.7, roughly consistent with the [27] scaling, such spectral slopes are indicative of local-in-scale turbulence [28]; c) the turbulence was observed to have scale-dependent anisotropy, again, in accord with [27].
As we noted above, locality implied that, if the ion Larmor radius (r L ) and the ion skin depth (d i ) are both below the equilibrium layer width, 0.015L, the bulk quantities, such as the average Simulations were seeded with small initial perturbations. The time evolution of the current layer width ∆ (right) and the inferred reconnection speed v r = d∆/dt (inset). From [24].
reconnection rate or the average dissipation rate, are not affected by plasma processes. This has already been established for reconnection due to ambient turbulence [29], but we extended this result to the case where reconnection develops spontaneously without an external agent. Unlike the 2D case, in which large plasmoids couple large and small scales directly, the 3D turbulence that we observe appears to be local-in-scale. This implies that MHD reconnection rates are universal [24].

The particle-MHD approach
Recently we developed a novel kinetic particle solver with a broad range of physical capabilities and coupled it fully dynamically with a Hall-MHD code. This code allows us to study environments which have significant mutual interaction between particles and the MHD fluid, such as reconnection with efficient particle acceleration. The code was developed from the ground up with a focus on scalability. Furthermore, it implements the most general prescription for particle motion. In particular, apart from the Lorentz force there is the possibility of imposing various stochastic diffusion terms intended to reproduce the subgrid diffusion, i.e., we solve for trajectories of pseudoparticles whose distribution satisfies a Fokker-Plank type of equation. Due to this, particle splitting and merging comes naturally in this code and it allows for better coverage of the particle energy space by splitting highly energetic particles. In addition, for all particles we use a flexible precision-controlled solver, which significantly improves accuracy and stability in situations when particle energies vary by orders of magnitude. The MHD part of the code solves the incompressible MHD equations with source terms from pseudoparticles. The pseudoparticle part solves stochastic differential equations (SDE), which are essentially equivalent to the Fokker-Planck equation [30]. The MHD equations are modified by adding −(J cr × B)/c + en cr (U × B)/c to the RHS of the momentum equation, where J cr , n cr are the particle current and density, obtained from volume-sampled pseudoparticle distribution first and zeroth moments. The first term reflects the fact that only the fluid part of the current should be included in the Lorentz force, while the second part is the Lorentz force due to extra electrons in the fluid to maintain quasi-neutrality. The term (cη∇ × J cr )/4π, where η is the magnetic diffusivity, could be added to the RHS of the induction equation and is due to the fact that only the fluid current is subject to resistivity. In most non-relativistic flows only the −J cr × B term is relevant and other terms could be omitted. The pseudoparticles are evolved using equations mcdu/dt = e(E + [u × B]/ 1 + u 2 ) + dp w /dt, dx/dt = uc + dx w /dt, where E and B are the MHD fields and u is a space part of the particle 4-velocity. The stochastic terms dp w /dt and dx w /dt are numerically calculated by increments in the frame associated with the local mean magnetic field, where the diffusion tensor has a diagonal √ ∆t, and are then transformed to the laboratory frame 2 . Here g * are Gaussian random numbers and∆ t is a time step.∆p w have similar form with momentum diffusion tensor instead of spacial diffusion tensor. Each pseudoparticle also carries "weight", which is accounted for, when calculating the RHS side of the MHD equations. Due to the weight, pseudoparticles could be "split", e.g., two pseudoparticles are created in place of one with the weight split between them. Such split particles follow different trajectories due to the stochastic part. The splitting is usually required in the acceleration scenario when higher energy particles are represented by fewer numbers, and splitting is necessary to maintain good statistics for high energies. Physically, stochastic terms represent evolution due to the unresolved, subgrid-scale electromagnetic perturbations. The so-called Monte Carlo simulations do not solve for the fluid flow, but use the prescribed velocity structure. In this case stochastic terms represent all magnetic perturbations, which are conjectured to exist but are not calculated explicitly by the model. Ordinarily Monte Carlo calculations use a simplistic prescription for the diffusion terms, such as Bohm scattering, implicitly assuming that scattering is very efficient. Our approach does not imply such strong assumptions and solves a coupled MHD-pseudoparticle system, which naturally includes particle-fluid interactions. The diffusion terms could be either set ad hoc as functions of energy or obtained from earlier local box simulations. In the reconnection experiments we describe below the particle diffusion terms has been set to zero.

Test results
Due to the flexibility of our code it can be used in a variety of circumstances, from a simple Monte Carlo solver to calculate particle dynamics in a prescribed flows to a fully self-consistent, dynamically relevant interplay between particles and the fluid. This flexibility also makes validation easier, as the separate parts of the code can be independently tested. The particle Note small-scale perturbations created in the lower half of the datacube by particles. Plot on the right shows evolution of the maximum energy (red) and average energy (green) of the particles. The particle distribution function also confirms that in this scenario all particles in the current layer are accelerated with the same rate, consistent with the theory presented in this paper.
code can be tested in a problem with a prescribed flow, the MHD code can be tested by introducing zero particles, and the test particle case can be studied by setting the pseudo-particle "weights" to zero. We already completed a variety of such tests, e.g. confirming the momentum conservation and energy conservation for the coupled system. Figure 4 show two tests. The first is the Monte Carlo test verifying that our particle part indeed solves Fokker-Plank equations and gives standard diffusive shock acceleration solutions for the test particle case. Note that this case feature spectra with 10 orders of magnitude range of energies. This is not unusual for the Monte Carlo with efficient particle splitting, however our code actually solves for trajectories. The second test features fully dynamic interplay between the fluid turbulence and particles, where energy is provided by driving the fluid, while particles are accelerated by the 2nd order mechanism and provide back-reaction to the fluid in the form of small-scale micro-turbulence, which in turn greatly enhances particle scattering. Figure 5 shows results of our simulations of spontaneous reconnection with particles. Note that coupled PIC-MHD system have four extra dimensionless numbers, characterizing the system, in addition to the MHD numbers. These are: m p , mass load of particles, the ratio of particle mass to the fluid mass, e p , energy load, the ratio of particle energy density to characteristic kinetic plus magnetic energy density of the fluid, r p , the ratio of Larmor radius of the particles to the dissipation/grid scale of the fluid. v p , the ratio of characteristic particle speed and characteristic fluid speed. Of these four numbers, r p is naturally set to values somewhat below one, to ensure that particles are accelerated from "thermal pool" with subgrid Larmor radii. Also v p is set to a large number (unless it limits numerical efficiency), as we are dealing with particles much faster that the fluid, e.g. relativistic particles. We typically inject particles with γ ∼ 1 that are further accelerated to relativistic energies. Finally, m p and e p are free parameters dependent on the problem. In simulations described on Figure 5 we chose a very small mass load of 10 −4 and a sizable initial energy load of 1.3 × 10 −2 , which increased by an order of magnitude during the

Discussion
In this paper we showed 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., in the beginning of spontaneous reconnection. This mechanism of acceleration of collisionless non-thermal particles by 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. Eq. (2) RHS will be zero. In the bulk fluid acceleration the energy gained by each particle does not depend on its initial energy, while for drift acceleration scenario it is proportional to the particle energy. 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 which produce local E∥B. Also, our mechanism is not tied to a special X-point, but instead volumetric.
The regular acceleration due to converging field lines have been suggested in [18], later it was pointed by [31] that an outflow cooling should also be included. In this paper we do not rely on simple transport equation, such as Parker's, therefore we relax the above approach's requirement that particle distribution need to be quasi-isotropic. In fact, in cases with significant acceleration, such as described in Section 8 we confirmed directly that particle distribution have significant deviations from isotropy locally.
The acceleration in turbulent reconnection has been further numerically studied in [13]. Plasma PIC simulations has 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 [32,22,23,33,34,35,36]. The change of energy due to curvature drift in a single collision of a particle with magnetic island was estimated analytically in [33]. 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 cutoffi s not related to the island size, instead it is related to the system size.
PIC simulations are limited in the range of scales and energies they cover. Recent simulations in [14,15] demonstrated acceleration up to 100 MeV in electron energies, which is below maximum energy in most astrophysical sources. Theory, therefore, is necessary to supplement conjectures based on the observed PIC distribution tails, explaining the underlying physical mechanism and making predictions for astrophysical systems which often feature gigantic scale separation between plasma scales and the size of the system. The feedback from simulations, nevertheless, was very useful, in particular the recent simulations [15] that reached MHD scales and confirmed the prediction that the curvature drift acceleration will dominate compared to the non-ideal electric field acceleration.