1 Introduction

# The Efficiency of Second-Order Fermi Acceleration by Weakly Compressible MHD Turbulence

## Abstract

We investigate the effects of pitch-angle scattering on the efficiency of particle heating and acceleration by MHD turbulence using phenomenological estimates and simulations of non-relativistic test particles interacting with strong, subsonic MHD turbulence. We include an imposed pitch-angle scattering rate, which is meant to approximate the effects of high frequency plasma waves and/or velocity space instabilities. We focus on plasma parameters similar to those found in the near-Earth solar wind, though most of our results are more broadly applicable. An important control parameter is the size of the particle mean free path relative to the scale of the turbulent fluctuations . For small scattering rates, particles interact quasi-resonantly with turbulent fluctuations in magnetic field strength. Scattering increases the long-term efficiency of this resonant heating by factors of a few-10, but the distribution function does not develop a significant non-thermal power-law tail. For higher scattering rates, the interaction between particles and turbulent fluctuations becomes non-resonant, governed by particles heating and cooling adiabatically as they encounter turbulent density fluctuations. Rapid pitch-angle scattering can produce a power-law tail in the proton distribution function but this requires fine-tuning of parameters. Moreover, in the near-Earth solar wind, a significant power-law tail cannot develop by this mechanism because the particle acceleration timescales are longer than the adiabatic cooling timescale set by the expansion of the solar wind. Our results thus imply that MHD-scale turbulent fluctuations are unlikely to be the origin of the tail in the proton distribution function observed in the solar wind.

plasmas – heating – acceleration of particles – (Sun:) solar wind
\DefineShortVerb

## 1. Introduction

Given the ubiquity of MHD turbulence in space and astrophysical plasmas, it is critical to understand the efficiency of particle heating and acceleration by MHD-scale turbulent fluctuations. In this paper we address this problem using numerical simulations of test particles interacting with strong MHD turbulence.

According to quasi-linear theory, high-gyrofrequency particles interacting with low-frequency, long-wavelength turbulence in collisionless plasmas do so primarily through resonant interactions with compressible modes (Achterberg, 1981). For low amplitude turbulent fluctuations, conservation of magnetic moment ensures that changes in perpendicular velocity are reversible, and so most acceleration is in the parallel direction, through the mirror force at and parallel electric fields at . Quasilinear theory predicts that only particles with precisely the wave velocity experience significant acceleration over long times (Kennel & Engelmann, 1966; Jokipii, 1966). However, because waves are not long-lived in realistic MHD turbulence, the resonance is significantly broadened, and many particles can experience acceleration (Bieber et al., 1994; Shalchi et al., 2004; Yan & Lazarian, 2008; Lynn et al., 2012).

As a result of this parallel diffusion, an initially isotropic distribution function becomes anisotropic over time, with more energy in the direction. However, as high-velocity particles are continually accelerated, the heating and acceleration of a distribution of particles becomes progressively less effective. This is because the heating efficiency is substantially lower for particles with velocities much larger than the resonant wave velocity. Moreover, the mirror force depends on but not .

In the solar wind and other astrophysical plasmas, small-scale microinstabilities such as the mirror, firehose, and ion cyclotron instabilities can become important if the distribution function becomes too anisotropic (Kasper et al., 2002; Sharma et al., 2007; Schekochihin et al., 2008; Bale et al., 2009). These instabilities will tend to quickly isotropize distributions functions upon their onset (Hellinger & Trávníček, 2008). A population of high frequency ion-cyclotron or whistler waves excited independent of velocity space instabilities can have the same effect for super-thermal particles.

If pitch-angle scattering is effective, particles which gain parallel velocity are efficiently isotropized, increasing the efficiency of particle heating and acceleration. If scattering rates are even higher, such that the mean free paths are smaller than the size of the turbulent eddies, then resonant wave-particle interactions no longer occur because scattering destroys the wave-particle phase coherence. Instead, particles are locked into the fluid motions and heat and cool adiabatically along with compressions and expansions of the fluid (Skilling, 1971). Particles then gain or lose energy irreversibly by spatially diffusing to neighboring uncorrelated eddies. This non-resonant process can contribute significantly to the overall heating rate. It has also been proposed as an efficient particle acceleration mechanism that can generate a population of non-thermal, high energy particles (Ptuskin, 1988). Recently, Fisk & Gloeckler (2008) argued that interactions of this kind with compressional turbulence could produce the tail in the proton distribution function observed at energies MeV in the solar wind and outer heliosphere.

In this paper, we determine the efficiency of particle heating and acceleration in MHD turbulence by implementing pitch-angle scattering of test particles evolving in numerical simulations of strong MHD turbulence; the scattering mimics the isotropization caused by plasma waves and/or instabilities. There is a substantial literature investigating the interaction of test particles with an ensemble of linear waves, intended to represent turbulence (e.g. Giacalone & Jokipii, 1994; Mace et al., 2000; Qin & Shalchi, 2012; Tautz et al., 2013). Our numerical results relax many of the simplifying assumptions of these simulations and of standard quasi-linear theory and thus provide a more robust estimate of particle energization by MHD-scale turbulent fluctuations in heliospheric and astrophysical plasmas. In our calculations, we restrict our analysis to non-relativistic particles and focus on low-frequency, weakly-compressible MHD turbulence that consists primarily of Alfvén and slow waves, with a small fraction of the energy in fast waves. This is qualitatively consistent with MHD turbulence in the solar wind (e.g., Yao et al. 2011; Howes et al. 2012).

In §2, we summarize phenomenological arguments from the literature on the heating and acceleration of test particles in the presence of turbulence and pitch-angle scattering. These analytic estimates provide valuable context for interpreting our test particle simulations. We describe our numerical methods in §3, and our test particle simulation results in §4. In §5 we discuss whether velocity diffusion in the presence of pitch-angle scattering can produce power-law tails in the distribution function. We then summarize our results and discuss their implications (§6).

## 2. Velocity diffusion: Analytic Estimates

The interaction of high-gyrofrequency test particles with low-frequency MHD turbulence leads to primarily parallel velocity diffusion. Perpendicular velocity diffusion is largely suppressed due to the approximate conservation of (see Chandran et al. 2010 for a detailed discussion of the conditions under which conservation is violated). Efficient pitch-angle scattering converts this parallel velocity diffusion into isotropic velocity diffusion, described by

 ∂f∂t=1v2∂∂v(v2D∂f∂v), (1)

where is the 3D distribution function, normalized by (Achterberg, 1981).

In this paper we focus on cases in which the pitch-angle scattering rate is sufficiently high that the distribution function remains isotropic. The timescale over which the distribution function at velocity would become significantly anisotropic in the case of primarily parallel velocity diffusion can be approximated by . Thus, the distribution function will be roughly isotropic at velocity if . In what follows we present results without an explicit velocity dependence of the pitch-angle scattering rate . However, it is quite plausible that in reality the pitch-angle scattering rate does depend on velocity, and, in particular, that thermal and super-thermal particles have different scattering rates. Because our numerical results in §4 are based on test particle calculations, the case of a velocity-dependent scattering rate can be obtained by interpolating between our numerically determined velocity diffusion coefficients for different values of .

Analytic treatments of the heating and acceleration of particles by MHD turbulence have considered primarily two distinct limits. We briefly review these results as they provide important context for interpreting our numerical simulations.

In the low scattering rate limit, the underlying interactions are essentially still resonant, as in a collisionless plasma. In order for an interaction with fluctuations on a scale to be resonant, the mean free path of a particle due to scattering must satisfy , where the mean free path is given by . This nearly resonant case is discussed in §2.1.

The second limit is when the particles are effectively collisional due to the scattering, i.e., . Particles now heat and cool adiabatically along with density fluctuations in the fluid. This limit is discussed in more detail in §2.2.

### 2.1. Resonant velocity diffusion

For low scattering rates, the interactions between turbulent fluctuations and particles are similar to those in the collisionless limit, but the scattering maintains a nearly isotropic distribution function. The isotropic diffusion coefficient is given by an appropriate average of the parallel velocity diffusion coefficient over the pitch-angle cosine ,

 D=12∫1−1dαα2D∥, (2)

where is the parallel diffusion coefficient in the absence of scattering. We assume that and will only depend on the isotropic speed, .

In the low scattering rate limit, the key interactions that set (and hence the isotropic diffusion coefficient) are approximate resonances between magnetic compressions and particles (also sometimes known as transit time damping, TTD). In quasi-linear theory, these magnetic compressions are due to fast or slow magneto-sonic waves and the resonances show up as strict delta functions of the form where is the phase speed of the wave with the largest magnetic compressions. In the case of strong MHD turbulence, however, quasilinear theory is quantitatively inapplicable for describing the wave-particle interactions, due to the short lifetimes of waves in turbulence and to the nonlinearity of the interaction. This has motivated a number of phenomenological models for resonance-broadened interactions between particles and fluctuations in MHD turbulence. In Lynn et al. (2012) we calibrated analytic expressions for including resonance broadening against test particle simulations. Among other things, we showed that the exponential time decorrelation of waves often assumed in the literature, which leads to a Lorentzian broadening of the wave-particle resonance, is not consistent with test particle simulations; instead a much more rapid Gaussian time decorrelation is required.

We parameterize the resulting diffusion coefficient as

 Unknown environment '% (3)

where is the phase velocity of the magnetic compressions (fast or slow waves in quasilinear theory) and for a resonance-broadening model with a Gaussian time decorrelation for the waves comprising the turbulence. The functional form in equation 3 is quite different from the of TTD in linear theory; resonance broadening allows many more particles to interact with the turbulent fluctuations. Using this simplified form for the diffusion coefficient, we perform the integral in equation 2 to estimate the isotropic velocity diffusion coefficient in the low scattering rate limit:

 D∝⎧⎪⎨⎪⎩8v4/105if vvp, (4)

For high velocities, equation 4 predicts .

In addition to the quasi-resonant diffusion described by equation 4, there is a contribution from non-resonant interactions with large-scale eddies. These interactions consist of particles being flung along moving, curved magnetic field lines, like beads on a wire, and are most important for low velocity (sub-thermal) particles. We term these Fermi Type-B (FTB) interactions. Diffusion coefficient in this limit are calibrated against test particle simulations in Lynn et al. (2012).

### 2.2. Non-Resonant Velocity Diffusion

If the scattering rate is sufficiently high that the mean free path of a particle becomes smaller than the size of an eddy (i.e., ), the particle acceleration mechanism becomes entirely non-resonant and the results of the previous section do not apply. Instead, particles become approximately “locked” into the scattering frame and undergo adiabatic increases (decreases) in energy as the large scale eddies converge (diverge); spatial diffusion out of the compressible fluctuations makes this energy change diffusive rather than adiabatic. Note that in this case it is the turbulent density fluctuations rather than the magnetic field compressions that the particles tap into. In particular, the important density changes are those associated with compressions and rarefactions of the fluid, as in sound waves, not density jumps due to pressure-balanced entropy modes. Quantitatively, this energy change may be expressed as

 dydt=−∇⋅u3, (5)

where is the fluid velocity, and we have introduced the variable , where is the initial momentum of a particle (Skilling, 1975).

Ptuskin (1988) derived the velocity diffusion coefficient for these non-resonant interactions assuming isotropic spatial diffusion and that the scatterers are at rest relative to the ambient fluid (see also, e.g., Chandran 2003 and Jokipii & Lee 2010 for related studies). The latter assumption is plausibly appropriate if non-propagating instabilities such as the firehose and mirror instabilities dominate pitch-angle scattering but not necessarily if (say) Alfvén waves do. We consider moving scatterers in §2.2.1.

Although spatial diffusion in MHD is primarily along magnetic fields lines (and thus anisotropic), Ptuskin’s calculation with isotropic diffusion captures many of the key results found in the more general MHD calculation. The results are simplest in the limit that particles diffuse across the largest scale eddies before the crossing time of compressible waves, , i.e., , where is the phase velocity of compressible waves, is the outer-scale of the turbulent fluctuations, and . A simplified phenomenological version of Ptuskin’s derivation, as applied to the outer-scale eddies, is as follows. The correlation time of the particle-turbulence interaction is the time required to cross one eddy, because neighboring eddies are uncorrelated: , where is the spatial diffusion coefficient. The rms change in as a particle diffuses across an eddy is given from Equation 5 by , where is the rms velocity divergence associated with the outer-scale (driving-scale) eddies. Thus, the resulting velocity diffusion coefficient is

 D∼v2⟨δy2⟩2τdiff∼v2l2D(∇⋅u)2D36κ   (τdiff≪τw), (6)

which is very similar to equation 21 of Ptuskin (1988), up to dimensionless coefficients and the identification of as a proportionality constant times .7 Note that for a fixed pitch-angle scattering rate , so that equation 6 predicts a diffusion coefficient that is independent of velocity.

In addition to the result in equation 6 valid in the limit , Ptuskin also derived the diffusion coefficient for when i.e., for (still assuming ). In this regime, the magnitude of generated by the waves is , where is the wave period and is the amplitude of compressible component of the velocity on a given scale . The fluctuations in momentum with amplitude are reversible until particles diffuse out of a given wave into an uncorrelated fluctuation, which happens on the timescale . The resulting diffusion coefficient is thus

 D≃v22u2l,cκ9v2wl2   (τdiff≫τw). (7)

This is consistent with the third inequality in eq. 28 of Ptuskin (1988) up to numerical factors of order unity. For a scattering rate that is independent of velocity, equation 7 implies that eddies of a given scale contribute a velocity diffusion coefficient . The net velocity diffusion coefficient in fully developed turbulence depends on the relative contribution of eddies on different scales. For density fluctuations associated with slow magnetosonic waves at , which roughly follow an anisotropic critical balanced cascade (Goldreich & Sridhar, 1995), this leads to a velocity diffusion coefficient (Chandran, 2003).

The results summarized here for the velocity diffusion coefficients produced by rapid pitch-angle scattering are consistent with those derived by Jokipii & Lee (2010) in the case of rapid spatial diffusion (our eq. 6; their eq. 20). We do not, however, agree with their results for slow spatial diffusion (our eq. 7; their eq. 19). Instead, our results in the limit of slow spatial diffusion are consistent with those originally derived by Ptuskin (1988).

For the magnetized plasmas of interest in this paper, the spatial diffusion is not isotropic but instead primarily along magnetic field lines. Chandran & Maron (2004) show that the resulting velocity diffusion coefficient differs somewhat from that of Ptuskin (1988); we briefly reproduce their argument here, focusing for simplicity on the case of rapid diffusion with . Because of the strong large-scale magnetic field, particles undergo a one-dimensional random walk along field lines and are likely to return to their original eddy several or more times. As a result, the correlation time for the particle-turbulence interaction is the decorrelation time of the turbulence, , rather than the diffusion time across an eddy. In one correlation time, a typical particle will traverse large-scale eddies, where is the spatial diffusion coefficient along magnetic field lines. The average velocity divergence felt by the particle during this longer correlation time is reduced by a factor of because the individual eddies are not correlated with one another. Thus the rms change in over a correlation time is given by

 δy∼13∫τcorrodt ∇⋅u∼13τcorr(∇⋅u)DN−1/2, (8)

and the resulting velocity diffusion coefficient is

 DNR∼v2lD18(∇⋅u)2D√τcorrκ∥, (9)

consistent with eq. 15 of Chandran & Maron (2004). For a scattering rate that is independent of velocity, equation 9 predicts , in contrast to the prediction of equation 6 for isotropic diffusion.

Equation 9 relies on the fact that for high gyrofrequency particles, spatial diffusion is primarily along field lines, and thus the particles undergo a 1D random walk, rather than a fully 3D random walk. This anisotropy, however, does not significantly affect the diffusion coefficient in the limit, given in equation 7, except that the eddy size must be interpreted as a parallel correlation length. For , even this is of little practical importance given that the compressible (fast mode) component of MHD turbulence cascades fairly isotropically (Kowal & Lazarian, 2010).

#### Non-Resonant Diffusion: Moving Scatterers

The Ptuskin (1988) and Chandran & Maron (2004) results summarized in the previous section apply if the scattering is isotropic in the local frame of the fluid. If, instead, the scatterers are randomly moving forward or backward in the fluid frame with a velocity (e.g., high-frequency Alfvén waves), and the particles scatter elastically in the scattering frame, non-relativistic particles will experience a change in velocity of due to the scattering. Given a scattering rate , the scattering will produce an additional contribution to the velocity diffusion of

 DFermi∼v2sν    (moving scatterers). (10)

This is essentially Fermi scattering (Fermi, 1949) in the non-relativistic limit. Note that because the velocity diffusion coefficient produced by interaction with turbulent density compressions scales as at high velocities and high scattering rates (eq. 9), the latter will dominate over the Fermi-like diffusion coefficient estimated in equation 10 at sufficiently high velocities.

An important distinction between scatterers that are stationary versus those moving relative to the fluid frame lies in the source of energy for the resulting velocity diffusion. Scattering that is elastic in the fluid frame does not transfer any energy to the particles and so the ultimate energy source for the velocity diffusion is the turbulence itself. By contrast, scattering that is elastic in a frame moving relative to the fluid frame leads to energy transfer from the scatterers to the test particles independent of the existence of turbulence. In our calculations with moving scatterers we specify a constant scattering rate and there is no back reaction on the scatterers associated with this energy transfer.

## 3. Numerical methods

Our simulations consist of charged test particles evolving in the macroscopic electric and magnetic fields of isothermal, subsonic MHD turbulence. Apart from our addition of explicit pitch-angle scattering, which we describe below, our computational approach is identical to that of Lynn et al. (2012). Dimensional quantities throughout the paper are expressed in units of the sound speed and the driving scale , when not explicitly stated (because we drive at multiple , the actual driving scale is slightly smaller than ; see Table 1).

### 3.1. Turbulence simulations

We simulate ideal MHD turbulence with the Athena code (Stone et al., 2008). The velocity field of the turbulence is driven solenoidally using an Ornstein-Uhlenbeck process, which has a characteristic autocorrelation time . Fiducial properties for the MHD simulations used in this work are summarized in Table 1, and any simulation with different parameters is explicitly noted. Our fiducial simulation parameters are broadly similar to those measured in the near-Earth solar wind.

In our calculations, we choose a correlation time similar to the outer scale eddy turnover time, to mimic a natural driving process. The simulation box is extended in the parallel direction along the mean magnetic field, because otherwise the particles (which obey periodic boundary conditions) would unphysically interact with the same eddies multiple times before the eddies decorrelate.

#### Measurement of turbulence properties

The predicted efficiency of non-resonant velocity diffusion, summarized in equation 9, depends on two properties of the turbulence, and . To order of magnitude, these can be approximated as and , respectively, for our turbulence. However, it is more accurate to measure these two parameters directly in the turbulence.

Recall that refers to the velocity divergence of the largest eddies. To measure this, we first Fourier transform the velocity field and then multiply by a lowpass window function of the form , where is chosen so that the window function is at . is defined by , where the average is taken over the power spectrum of the driving. Our driving power spectrum is a power-law in with a 1D power spectrum of defined between and , which gives . We then calculate the rms velocity divergence, which may be done easily in Fourier- or real-space. This yields for our fiducial simulation summarized in Table 1.

To calculate the correlation time, we choose one point per cpu to output a time-series of the velocity magnitude. The point chosen is fixed with respect to the origin of the local grid. Because the local grids periodically tile the entire domain, these points are arranged in a lattice. There are 512 points in our fiducial case, which should be approximately independent, as there is approximately one per physical eddy volume. For each of these lattice points, we calculate the autocorrelation of the time series in using the —sfsmisc— package in R (Maechler et al., 2012; R Core Team, 2012). These autocorrelation functions are then averaged over all lattice points, and fit with a functional form , which gives for our fiducial simulation (see Table 1). Quantitatively, this is roughly consistent with an estimate of .

Finally, we can decompose the turbulence into components corresponding to the linear MHD Alfvén, slow, and fast modes, following the Fourier space method of Cho & Lazarian (2003). This method makes the simplifying assumption that the magnetic field points along an axis of the simulation domain, which will be asymptotically appropriate as . For our fiducial simulation summarized in Table 1, we find that 50%, 45%, and 5% of the turbulent energy is in the Alfvén, slow, and fast modes respectively. This is somewhat more slow and fast mode energy than measured in the near-Earth solar wind (Howes et al., 2012), likely because of the strong plasma heating (and thus turbulent dissipation) produced by the compressible component of the turbulence (see §4). An alternative decomposition of the turbulence is into solenoidal and compressive () fluctuations: we find that of the energy is in compressive fluctuations. This contains both slow and fast mode contributions for our turbulence. The corresponding rms density fluctuation is , somewhat larger than the in situ value of measured in the solar wind (e.g., Tu & Marsch 1995).

#### Effect of limited dynamic range

Given our periodic boundary conditions and a box of finite size, a test particle will cross the box multiple times. If a particle returns to nearly the same location following a magnetic field line, before the turbulence has had time to randomize, it will experience artificial correlations in the turbulence and therefore artificial acceleration along its trajectory. To minimize this effect, we extend the box in the parallel direction, so that particles interact with many uncorrelated eddies as they cross the box. However, for fixed computational time, this effectively means that we are reducing the dynamic range of the simulation in terms of the number of simulation elements that comprise each turbulent eddy. Figure 1 shows the one-dimensional kinetic energy power spectrum of the turbulence in our fiducial simulation. A clear consequence of extending the box along the mean magnetic field is that the inertial range of the turbulence is modest, less than a decade in k. We have found, however, that not extending box along the mean field significantly changes the acceleration of high energy particles, so this compromise of limited dynamic range is necessary.

Our analytical estimates in §2 are given in terms of arbitrary power spectra for the turbulent fluctuations. In many cases, the diffusion and acceleration of particles in the presence of pitch-angle scattering is dominated by interactions with a particular scale in the turbulence, typically either the smallest or largest scale in the inertial range. In this case, simulating a large inertial range may not be critical, though it would obviously be preferable. Later in the paper we highlight the specific results that we suspect may be most influenced by our limited dynamic range.

### 3.2. Test particle integration

Our particle integration methods are described in detail in Lehe et al. (2009). Particles are initialized in the fully saturated turbulence, and then evolved according to the Lorentz force. The and -fields are those on the MHD grid, interpolated to the particle’s location using the triangular-shaped cloud (Hockney & Eastwood, 1981) method in space and time. Particles are integrated using the Boris (1970) particle pusher, which is sympletic and symmetric in time, and conserves energy and the magnetic moment adiabatic invariant to machine precision in tests with constant fields. We initialize particles randomly across the simulation domain with the initial values of and defined relative to the local rest-frame of the MHD fluid in the simulation. More details are provided in Lehe et al. (2009).

#### Pitch-angle scattering

We implement pitch-angle scattering by specifying a scattering rate, . We consider two cases to bracket the likely range of physical scattering mechanisms. The first is when the scatterers are at rest in the fluid frame. This seems likely to be appropriate for scattering by mirror or firehose instabilities, but not by, e.g., Alfvén waves, which move relative to the fluid frame at the Alfvén speed. We consider scatterers moving along the local magnetic field direction to model the latter.

At each scattering timestep (defined below), and for each test particle, we subtract off the local fluid velocity, to perform the scattering in the fluid frame. For the simulations with moving scatterers, we move into the frame of the scatterer by subtracting off a further component of the velocity , where is the local magnetic field direction and the scatterer is randomly chosen to be moving forward or backward in the fluid frame. We choose an isotropic random unit vector . We then calculate , i.e., a unit vector perpendicular to the particle velocity . We rotate the test particle’s velocity vector around through a random angle uniformly distributed between 0 and . is chosen to be perpendicular to so that the rotation of occurs along a great circle on the unit sphere. We choose the timestep for the scattering procedure according to . We implement this rotation using Rodrigues’ rotation formula. Finally, the local fluid velocity (in addition to any scatterer velocity as appropriate) is added back to the particle’s velocity to return us to the lab frame.

This procedure ensures that the pitch-angle scattering rate is in fact given by our input parameter as long as , so that scattering is dominated by small changes in pitch-angle. We typically choose , which is sufficiently small that diffusion coefficients differ by less than 1% from their saturated values at , while minimizing computational effort associated with frequent scatterings.

In addition to diffusion in velocity space, pitch-angle scattering also leads to spatial diffusion with a mean free path and a diffusion coefficient , where the parallel subscript indicates that the spatial diffusion is still primarily along magnetic field lines in our calculations. In particular, perpendicular spatial diffusion is , where is the particle gyroradius. We simulate high-gyrofrequency particles with small and so perpendicular diffusion is negligible even at the highest we consider.

### 3.3. Calculating Velocity Diffusion Coefficients

The velocity diffusion coefficients are calculated according to the formal definition

 D≡⟨δv2⟩2δt, (11)

where the average is over many particles with the same initial velocity. However, one subtlety is that a particle which experiences no acceleration can nevertheless appear to change velocity when measured with respect to the local fluid frame. This is because the fluid’s parallel velocity can change even if a given test particle feels no electromagnetic acceleration. Thus, absent pitch-angle scattering, after a time of order the correlation time of the fluid, an ensemble of particles initially at rest relative to the fluid will develop a finite dispersion in parallel velocity given by of the turbulence. (This effect becomes less significant at high pitch-angle scattering rates because the particles are ‘pinned’ near the location in the fluid where they are initialized.) Figure 2 plots two representative curves and shows a rapid increase in at early times for the low scattering rate calculation. To not include this in our diffusion coefficient calculations, we measure diffusion coefficients by instead fitting a straight line after the initial jump using least-squares (as indicated in Fig. 2).

## 4. Test-particle Diffusion and Heating

### 4.1. Velocity Diffusion: Stationary Scatterers

Figures 3 and 4 show diffusion coefficients for a broad range of pitch-angle scattering rates for our fiducial simulation (Table 1); the scattering is isotropic in the fluid frame. The chosen values of include scattering times significantly longer than, and significantly shorter than, the outer-scale turnover time of the turbulence. At the lowest scattering rates (, corresponding to scattering times somewhat smaller than the nonlinear timescale of the turbulent fluctuations), the diffusion coefficient is well approximated by the pitch-angle average of the parallel diffusion coefficient appropriate for resonance broadened TTD + Fermi Type-B diffusion described in §2.1 (cyan line in Fig. 3). In this limit, scattering isotropizes the velocity-space diffusion but does not fundamentally alter the physics of how particles interact with the turbulent fluctuations.

For scattering rates , the diffusion coefficient increases at low velocity but not significantly for . The reason is that large-scale eddies produce the Fermi Type-B diffusion at low velocities – these interactions are thus fundamentally altered once the scattering rate is comparable to the frequency of the large-scale eddies. By contrast, eddies of all scales (including, in particular, small scales) contribute significantly to the high velocity TTD diffusion. The small scale eddies have higher nonlinear frequencies and thus a higher scattering rate is needed to modify the velocity-space diffusion of high speed particles.

For , there is a qualitative change in the properties of the velocity diffusion coefficients in Figures 3 and 4, particularly for particles with . This corresponds to the onset of non-resonant diffusion governed by interactions between particles and turbulent density fluctuations. The properties of the high scattering rate diffusion depend on the ratio of the diffusion time to the wave period . Lower velocity particles are in the regime where and equation 7 predicts that the velocity diffusion coefficient is a strong function of velocity, with . This is consistent with the numerical results in Figures 3 and 4 for high scattering rates and intermediate velocities . Note that we do not numerically recover Chandran (2003)’s prediction that compressible slow mode turbulence at should produce a scaling due to the contribution of eddies of different scale in equation 7. This may be because the finite dynamic range in our simulations means that we do not capture the full inertial range of the turbulence (this is exacerbated by our anisotropic box [see Table 1], which is required in order to avoid artifacts due to periodic boundary conditions along the mean field). In the limit that eddies of a single spatial scale dominate the velocity diffusion, equation 7 implies , consistent with our numerical results.

For yet higher velocity particles, and equation 9 predicts in the rapid scattering limit. This is consistent with the numerical results in Figures 3 and 4. Quantitatively, the non-resonant regime of equation 9 requires simultaneously satisfying and , which restricts us to the following range of velocities:

 νl2D/τw≪v2≪ν2l2D. (12)

For the simulation with , the high-velocity non-resonant diffusion is appropriate in our test particle calculations for , while for , particles with should satisfy equation 9. In Figure 4, we include the prediction of equation 9 for the simulation as a dashed curve (measuring and as in §3.1.1 to determine the predicted normalization of ). The velocity diffusion is reasonably well approximated by equation 9 in the appropriate velocity interval.

Figure 5 shows how the diffusion coefficient in the high scattering rate, non-resonant regime depends on (for fixed ). The turbulence has the same driving as in the fiducial simulation in Table 1. The key difference in the resulting turbulence properties is that the fraction of the turbulent energy in compressional fluctuations (those ; see §3.1.1) is only for versus and for and , respectively. This is because fast modes dominate the compressional fluctuations at high and are not efficiently excited in subsonic turbulence with solenoidal driving.

The results for the diffusion coefficient in Figure 5 are similar for and . This is consistent with the fact that at low slow magnetosonic modes dominate the compressional fluctuations and are energetically important in MHD turbulence. By contrast, the diffusion coefficient changes significantly at , where fast modes dominate the compressional fluctuations. In particular, the results with in Figure 5 are similar to the results at a much lower (see Fig. 3). This is because of the much lower compressional energy in the simulation, which decreases the efficiency of the non-resonant diffusion (see §2.2).

### 4.2. Velocity Diffusion: Moving Scatterers

Figure 6 shows test particle results from simulations with moving scatterers. The velocity diffusion coefficient is nearly independent of particle velocity. These results are consistent with Fermi-like scattering predicted by equation 10, namely . At sufficiently high velocity, the non-resonant diffusion shown in Figures 3 and 4 would eventually dominate over the roughly constant diffusion coefficient shown in Figure 6.

### 4.3. Test-particle heating

Figure 7 shows how the pitch-angle scattering rate affects the heating of test particles that have an initially Maxwellian distribution with thermal speed . The results shown are nearly identical to the heating rates in the absence of pitch-angle scattering. For the calculations in Figure 7, we assume that the scattering is at rest in the local fluid frame. Note that for our fiducial simulation parameters (Table 1), protons correspond to , electrons to and minor ions (with temperatures of order the proton temperature) to in Figure 7.

Not surprisingly, the enhanced velocity diffusion shown in Figures 3 and 4 leads to a corresponding increase in the net particle heating rate. At the scattering rates we focus on here, this increase in the heating rate is a consequence of pitch-angle scattering enforcing isotropy, which leads to a substantially enhanced TTD heating rate at and an enhanced FTB heating rate at low velocities. Quantitatively, Figure 7 shows that even a modest scattering rate, with comparable to the outer-scale turbulent correlation time (; see Table 1), can increase the net heating rate significantly so that a large fraction of the compressible turbulent energy can go into heating protons. Moreover, when for protons, our test particle assumption breaks down because the protons will absorb the majority of the compressible turbulent power on large scales in the turbulent cascade. Figure 7 also shows that minor ion heating is particularly effective if the ions and protons have comparable temperatures, so that the minor ions have lower thermal velocities. This is true, e.g., in the slow speed solar wind but not in the fast solar wind, where the ion and proton thermal velocities are comparable (Bochsler, 2007); in this case, our results imply comparable proton and minor ion heating rates per particle.

## 5. Power-law tails?

Figure 8 shows the late-time distribution function that results when we evolve test particles which initially have a Maxwellian with in our fiducial turbulence simulation (Table 1), for several different scattering rates. The evolution is for , which is many eddy turnover times. These results are very similar to the analogous results without pitch-angle scattering in our earlier work (Lynn et al., 2012). In particular, Figure 8 demonstrates that at low to moderate scattering rates the particles are efficiently heated by the turbulence (more-so with increasing ; see Fig. 7), but no significant non-thermal power-law tail develops. Recall that this scattering rate regime corresponds to the particles gaining energy from turbulent magnetic compressions (as in linear TTD of fast and slow magnetosonic waves).

Figure 9 shows the evolution of the distribution function for for two initial thermal velocities, and (these two different thermal velocities are included to highlight the physics of the distribution function evolution at high ; below we also summarize results for calculations with , more relevant for protons). In the high scattering rate calculations in Figure 9, the resulting distribution function is qualitatively different: there is a nearly power-law component over almost a decade in velocity, particularly for . The non-thermal tail also becomes more energetically important at later times, after the distribution function has been evolved for . There is no evidence for a universal power-law slope in the velocity distribution function in this regime of rapid pitch-angle scattering. Instead, the exact slope depends sensitively on time, the pitch-angle scattering rate, and the thermal velocity of the test particle distribution: for , at and at while for , at late times and is not convincingly a power-law at earlier times.

The significant non-thermal tail to the distribution function in Figure 9 is caused by the steep scaling of the velocity diffusion coefficient with for intermediate velocities seen in Figure 4.15 The ‘strength’ of this power-law tail depends, however, on the assumed pitch-angle scattering rate and the thermal speed of the test particles since this determines how well the diffusion coefficient regime in Figure 4 overlaps with the initial Maxwellian of the test particles. Higher scattering rates lead to setting in at higher velocity. This implies significantly less efficient acceleration of particles out of the thermal population. Thus a strong power-law tail requires fine-tuning between the scattering rate and the thermal velocity of the species of interest. Quantitatively, the power-law tail in Figure 9 is much less significant for relative to . In addition, in calculations with (more appropriate for protons), which are not explicitly shown here, we find that after interacting with our fiducial turbulence (Table 1), the distribution function roughly satisfies for , for , and for over the velocity range . This highlights the strong dependence of the resulting distribution function on all of the parameters of the problem.

The absence of a non-thermal tail at low scattering rates in Figure 8 is consistent with the general form of the diffusion coefficient in Figure 3, in particular the fact that the diffusion coefficient scales at most . Green’s function solutions of the diffusion equation show that only if the diffusion coefficient scales as or steeper does velocity diffusion alone generically produce a strong non-thermal tail to the distribution function (e.g., Gruzinov & Quataert 1999; Jokipii & Lee 2010). This is consistent with our numerical results in Figure 8.

In his original analysis, Fermi invoked a competition between velocity diffusion and particle escape in order to obtain a power-law distribution function. Our simulations do not explicitly have escape and one might thus worry that the absence of a power-law tail is an artifact of the periodic boundary conditions in our calculations. This is not the case.

In the presence of particle escape, the combination of diffusion and escape produces a power-law tail only if the ratio of the acceleration time to the escape time is independent of particle energy (e.g. Longair, 1992). In our problem, a well-posed question is whether escape due to spatial diffusion, which is self-consistently produced by pitch-angle scattering, produces the requisite conditions for a power-law tail. Scattering implies so that . If is roughly independent of velocity, only a very steep diffusion coefficient of the form can produce a power-law tail in the presence of diffusive escape. This is the regime in which we already find a significant non-thermal tail in Figure 9. The velocity diffusion coefficient that results from most of our calculations is (roughly) , i.e., . This applies in both the quasi-resonant limit and in the high scattering rate limit at high velocities. Thus, unless for some reason , so that , weakly compressible MHD turbulence will not generically produce a significant non-thermal tail in the regime even in the presence of particle escape.

## 6. Conclusions

We have studied the heating and acceleration of particles by subsonic MHD turbulence by evolving an ensemble of test particles in real time in MHD turbulence simulations. Our focus in this paper has been on the role of pitch-angle scattering in modifying the efficiency of non-relativistic particle heating and acceleration. To study this we have implemented explicit pitch-angle scattering (at a fixed rate ), with the scatterers either at rest in the fluid frame or moving relative to the fluid frame along the local magnetic field. We suspect that these two cases bracket reality in problems of interest, representing the role of mirror/firehose instabilities (which are non-propagating instabilities and thus roughly at rest in the fluid frame) and high frequency ion cyclotron, Alfvén, or whistler waves, respectively. In our calculations the scattering rate is independent of velocity, but this may not be the case in real systems. In particular, the mechanisms that dominate pitch-angle scattering may differ for thermal and super-thermal particles. Results for a velocity-dependent scattering rate can be obtained from the calculations presented here by appropriately interpolating between the diffusion coefficients for different scattering rates shown in Figures 3-5 (this is formally valid for our test particle calculations).

In all of our calculations, we have focused on the evolution of particles with gyrofrequencies much larger than the resolved frequencies of the turbulent fluctuations. In the absence of pitch-angle scattering, the magnetic moment of our test particles is thus conserved: particles only undergo diffusion in . The reason for restricting our analysis to these high gyrofrequency particles is that the outer scale turbulent fluctuations in astrophysical and heliospheric systems have frequencies much smaller than the particle gyrofrequencies. Thus the high gyrofrequency limit is the correct limit for studying particle heating and acceleration using MHD simulations. In addition, much of the energy in MHD turbulence cascades anisotropically in wave-vector space, and there is little power at frequencies comparable to the ion cyclotron frequency, even when the turbulent fluctuations have wavelengths comparable to the ion Larmor radii (e.g., Schekochihin et al. 2009). A corollary of this is that it is difficult to predict a priori what the pitch-angle scattering rate should be in a given astrophysical or heliospheric context, since it depends on kinetic processes that excite very high frequency fluctuations (relative to the MHD-scale turbulence that is better understood). We return to this point below when discussing the implications of our results for the solar wind.

We find that the physics of particle energization by MHD turbulence depends on the ratio of the mean free path set by pitch-angle scattering () to the scale of the turbulent fluctuations . For low scattering rates, i.e., , the particle energization is analogous to that in the collisionless limit and occurs primarily through quasi-resonant interactions with turbulent fluctuations in the local strength of the magnetic field. In linear theory, these field compressions would be associated with slow () or fast () magnetosonic modes and the resulting wave-particle interaction is often called transit-time damping (TTD). Pitch angle scattering increases the efficiency of this quasi-resonant heating mechanism by reducing the parallel acceleration of particles out of resonance and converting heating into heating, thus increasing the magnitude of the force associated with magnetic field compressions (Achterberg, 1981). The net heating rate for thermal particles can be a factor of larger in the presence of modest pitch-angle scattering (; Fig. 7).

In the limit of rapid pitch-angle scattering, the physical process leading to particle energization is fundamentally different. Our numerical results in this limit are consistent with the analytic calculations of Ptuskin (1988) and Chandran & Maron (2004). Physically, at high scattering rates, particles are nearly pinned to the frame of the scatterers. To first approximation, particles thus adiabatically heat and cool as they move through density fluctuations in the plasma. Spatial diffusion of particles relative to the density fluctuations ultimately makes this energy change diffusive rather than reversible. Note that the energy source for particle energization in the rapid scattering limit is turbulent density compressions and rarefactions, rather than the turbulent magnetic field compressions that are important at low scattering rates. In linear theory, these density fluctuations are associated with fast modes at and both slow and fast modes at . In addition to coupling to the turbulent density compressions, if the scatterers are moving relative to the fluid frame, pitch-angle scattering produces a Fermi-like velocity kick that corresponds to a velocity diffusion coefficient that is independent of velocity for non-relativistic particles (eq. 10 and Fig. 6).

One of our primary motivations for including the effects of pitch-angle scattering was to assess whether pitch-angle scattering increases the efficiency of particle acceleration by MHD turbulence, i.e., the generation of a non-thermal tail in the distribution function. We find that this is not the case at low scattering rates where quasi-resonant interactions with magnetic field compressions dominate the velocity diffusion (see Fig. 8). Although particle heating is quite effective, with heating rates for thermal protons comparable to the energy cascade rate in the turbulence, the distribution function does not develop a significant non-thermal tail. This conclusion from our test particle simulations is consistent with solutions of the Fokker-Planck equation given the velocity diffusion coefficient we find at low scattering rates, namely for (Fig. 3).

For high scattering rates, , we find that the distribution function can develop a non-thermal tail, but only under certain relatively restrictive conditions. The power-law tail arises due to the strong scaling that occurs at intermediate velocities in the rapid scattering limit (eq. 7 and Fig. 4). At high particle velocities, the diffusion coefficient varies much more weakly with velocity () and so does not produce efficient acceleration. Physically, the transition between and occurs when the spatial diffusion time () across an eddy becomes shorter than the eddy turnover time, as predicted in Ptuskin (1988)’s original analysis. This produces a significant power-law tail only if the transition between the two diffusion regimes occurs at a velocity comparable to the thermal velocity of the particles of interest (Fig. 9). This requires fine tuning between the pitch-angle scattering rate and the thermal velocity of the particles.

What are the implications of our results for particle heating and acceleration in the solar wind, the most well-studied and well-constrained astrophysical plasma to which our results can be applied? This depends to some extent on the effective scattering rate in the low-collisionality solar wind. Given the ubiquitous in situ evidence for significant temperature anisotropies in the thermal plasma in the solar wind (e.g., Kasper et al. 2002) it is likely that the effective scattering rate for the thermal plasma cannot be much larger than the expansion rate of the solar wind (where is the local radius, i.e., AU near Earth). On the other hand, the measured temperature anisotropies in the solar wind also appear to be constrained by the thresholds for mirror, firehose, and (to a lesser extent) ion cyclotron instabilities, suggesting that there is some minimum level of pitch-angle scattering even in nearly collisionless epochs. Expansion of the solar wind inevitably drives solar wind plasma towards these instability thresholds (e.g., Hellinger & Trávníček 2008). Together, these arguments suggest that the effective scattering rate for thermal plasma is , which is in turn where is the outer-scale of turbulence in the near-Earth solar wind (e.g., Howes et al. 2008) and . Measurements of anisotropies in the distribution function at (at distances of a few AU) suggest that the mean free path is of order an AU (Gloeckler et al., 1995), implying a comparable effective scattering rate for super-thermal particles. These estimates are probably only accurate to a factor of , but this is sufficient to draw several important conclusions.

First, pitch-angle scattering likely enhances heating of the thermal particles by the compressible part of the cascade (which is primarily slow modes in the near-Earth solar wind; Yao et al. 2011; Howes et al. 2012) by factors of a few-10 (Fig. 7). Secondly, Figure 8 shows that even in the presence of pitch-angle scattering with , the particle distribution function does not develop a significant non-thermal tail due to MHD-scale turbulent fluctuations. Moreover, even if we optimistically assume that (in spite of the evidence for anisotropic distribution functions in the solar wind), our results do not imply that rapid pitch-angle scattering would produce a power-law tail in the proton distribution function in the near-Earth solar wind. The reason is that the particle acceleration time implied by Figure 4 () is longer than the expansion time of the solar wind over the entire velocity regime in which . This is even more true when we take into account that the compressional energy in the solar wind – which is the important energy source in the presence of pitch-angle scattering – is somewhat smaller than in our simulations: e.g., in the solar wind vs. in our fiducial simulation (see §3.1.1). Thus in the context of the near-Earth solar wind, adiabatic cooling (not present in our calculations) dominates particle acceleration by MHD turbulence even for very rapid pitch-angle scattering; this would preclude the power-law tail seen in Figure 9 from developing.

Taken together, these considerations strongly argue against the MHD scale fluctuations as the source of the tail in the proton distribution function observed throughout the heliosphere (Gloeckler, 2003; Fisk & Gloeckler, 2006). It does not, of course, preclude that kinetic-scale physics acting on small scales in the turbulent cascade could produce such a non-thermal tail (e.g., via reconnection or kinetic Alfvén wave turbulence). In addition, our conclusions are based primarily on the well-studied plasma conditions near AU. If the pitch-angle scattering rate and/or the compressibility of the turbulence are significantly larger at larger radii in the solar wind, the compressible part of the turbulence may be important for producing a significant non-thermal tail.16

Although we have focused the discussion of our results on the solar wind application because of the wealth of in situ measurements, the calculations presented here are much more broadly applicable. For example, particle heating and acceleration by MHD turbulence has been widely used to model solar flares and the “re-acceleration” of cosmic rays in the Milky Way and clusters of galaxies. Some of these applications would require generalizing our techniques to relativistic particles, but this is straightforward.

We thank Phil Isenberg, Marty Lee, and Anatoly Spitkovsky for useful conversations. This material is based on work supported by the National Science Foundation Graduate Research Fellowship under Grant No. DGE-1106400. Additionally, this work was supported in part by NSF-DOE grant PHY-0812811, NASA HTP grant NNX11AJ37G, and NSF grant ATM-0752503. EQ was supported in part by a Simons Investigator award from the Simons Foundation, the David and Lucile Packard Foundation, and the Thomas Alison Schneider Chair in Physics. Computing time was provided by the National Science Foundation TeraGrid/XSEDE resource on the Frost, Trestles, and Kraken supercomputers.

### Footnotes

1. affiliation: Physics Department, University of California, Berkeley, CA 94720; jacob.lynn@berkeley.edu
2. affiliation: Astronomy Department and Theoretical Astrophysics Center, University of California, Berkeley, CA 94720
3. affiliation: Astronomy Department and Theoretical Astrophysics Center, University of California, Berkeley, CA 94720
4. affiliation: Space Science Center and Department of Physics, University of New Hampshire, Durham, NH 03824
5. affiliation: Astronomy Department and Theoretical Astrophysics Center, University of California, Berkeley, CA 94720
6. affiliation: Present address: Canadian Institute for Theoretical Astrophysics, 60 St. George Street, University of Toronto, Toronto, ON MS H, Canada
7. Ptuskin’s notation for spatial and velocity diffusion coefficients is the opposite of ours (we label them as and , respectively).
8. The turbulent energy input rate. This corresponds to a sonic Mach number of .
9. Ratio of thermal to magnetic pressure.
10. refers to the correlation time in the Ornstein-Uhlenbeck turbulence forcing, while is the measured correlation time in the saturated state of the turbulence. We choose so that the driving does not artificially reduce the correlation time.
11. Outer (driving) scale of the turbulence.
12. RMS density fluctuations in the saturated turbulence.
13. RMS magnetic field fluctuations in the saturated turbulence.
14. Test particle gyrofrequency.
15. As noted in §4.1, Chandran (2003) argued that for a large turbulent inertial range, compressible slow mode fluctuations at would lead to a diffusion coefficient , rather than the somewhat steeper result we find numerically (Fig. 4). Even if this were to be the case, we believe that a significant power-law tail would develop, based on analytic and numerical solutions of the Fokker-Planck equation with .
16. There is some ambiguity in the literature on the radial evolution of the turbulence in the solar wind. Based on Ulysses data, Breech et al. (2009) quote a factor of decrease in from AU, which suggests that the turbulence is as or more compressible at large radii as at 1 AU. However, using Voyager 2 data Ng et al. (2010) find that decreases by a factor of over the same radial range, suggesting significantly lower compressibility at large radii. By contrast, direct measurements of the density fluctuations using Voyager 2 data suggest that increases with increasing heliospheric radius, albeit with somewhat large error bars (Bellamy et al., 2005).

### References

1. Achterberg, A. 1981, Astronomy & Astrophysics, 97, 259
2. Bale, S., Kasper, J., Howes, G., Quataert, E., Salem, C., & Sundkvist, D. 2009, Physical Review Letters, 103, 211101
3. Bellamy, B. R., Cairns, I. H., & Smith, C. W. 2005, Journal of Geophysical Research (Space Physics), 110, 10104
4. Bieber, J. W., Matthaeus, W. H., Smith, C. W., Wanner, W., Kallenrode, M.-B., & Wibberenz, G. 1994, The Astrophysical Journal, 420, 294
5. Bochsler, P. 2007, A&A Rev., 14, 1
6. Boris, J. 1970, in Proceedings of the Fourth Conference on Numerical Simulations of Plasmas, Naval Research Lab, 3–67
7. Breech, B., Matthaeus, W. H., Cranmer, S. R., Kasper, J. C., & Oughton, S. 2009, Journal of Geophysical Research (Space Physics), 114, 9103
8. Chandran, B. D. G. 2003, The Astrophysical Journal, 599, 1426
9. Chandran, B. D. G., Li, B., Rogers, B., Quataert, E., & Germaschewski, K. 2010, The Astrophysical Journal, 720, 503
10. Chandran, B. D. G., & Maron, J. L. 2004, The Astrophysical Journal, 603, 23
11. Cho, J., & Lazarian, A. 2003, Monthly Notices of the Royal Astronomical Society, 345, 325
12. Fermi, E. 1949, Physical Review, 75, 1169
13. Fisk, L., & Gloeckler, G. 2006, The Astrophysical Journal Letters, 640, 79
14. —. 2008, The Astrophysical Journal, 686, 1466
15. Giacalone, J., & Jokipii, J. R. 1994, The Astrophysical Journal, 430, L137
16. Gloeckler, G. 2003, AIP Conference Proceedings, 679, 583
17. Gloeckler, G., Schwadron, N. A., Fisk, L. A., & Geiss, J. 1995, Geophys. Res. Lett., 22, 2665
18. Goldreich, P., & Sridhar, S. 1995, The Astrophysical Journal, 438, 763
19. Gruzinov, A., & Quataert, E. 1999, The Astrophysical Journal, 520, 849
20. Hellinger, P., & Trávníček, P. M. 2008, Journal of Geophysical Research, 113, A10109
21. Hockney, R., & Eastwood, J. 1981, Computer Simulation Using Particles (CRC Press)
22. Howes, G. G., Bale, S. D., Klein, K. G., Chen, C. H. K., Salem, C. S., & TenBarge, J. M. 2012, The Astrophysical Journal, 753, L19
23. Howes, G. G., Cowley, S. C., Dorland, W., Hammett, G. W., Quataert, E., & Schekochihin, A. A. 2008, Journal of Geophysical Research, 113, 1
24. Jokipii, J. R. 1966, The Astrophysical Journal, 146, 480
25. Jokipii, J. R., & Lee, M. A. 2010, The Astrophysical Journal, 713, 475
26. Kasper, J. C., Lazarus, A. J., & Gary, S. P. 2002, Geophysical Research Letters, 29, 1839
27. Kennel, C., & Engelmann, F. 1966, Physics of Fluids, 9, 2377
28. Kowal, G., & Lazarian, a. 2010, The Astrophysical Journal, 720, 742
29. Lehe, R., Parrish, I. J., & Quataert, E. 2009, The Astrophysical Journal, 707, 404
30. Longair, M. S. 1992, High Energy Astrophysics, 2nd edn. (Cambridge University Press)
31. Lynn, J. W., Parrish, I. J., Quataert, E., & Chandran, B. D. G. 2012, The Astrophysical Journal, 758, 78
32. Mace, R. L., Matthaeus, W. H., & Bieber, J. W. 2000, The Astrophysical Journal, 538, 192
33. Maechler, M., et al. 2012, —sfsmisc—: Utilities from Seminar fuer Statistik ETH Zurich, R Package Version 1.0-20
34. Ng, C. S., Bhattacharjee, A., Munsi, D., Isenberg, P. A., & Smith, C. W. 2010, Journal of Geophysical Research (Space Physics), 115, 2101
35. Ptuskin, V. 1988, Soviet Astronomy Letters, 14, 255
36. Qin, G., & Shalchi, A. 2012, Advances in Space Research, 49, 1643
37. R Core Team. 2012, R: A Language and Environment for Statistical Computing, R Foundation for Statistical Computing, Vienna, Austria, ISBN 3-900051-07-0
38. Schekochihin, A., Cowley, S., Kulsrud, R., Rosin, M., & Heinemann, T. 2008, Physical Review Letters, 100, 081301
39. Schekochihin, A. A., Cowley, S. C., Dorland, W., Hammett, G. W., Howes, G. G., Quataert, E., & Tatsuno, T. 2009, The Astrophysical Journal Supplement Series, 182, 310
40. Shalchi, A., Bieber, J. W., Matthaeus, W. H., & Qin, G. 2004, The Astrophysical Journal, 616, 617
41. Sharma, P., Quataert, E., Hammett, G. W., & Stone, J. M. 2007, The Astrophysical Journal, 667, 714
42. Skilling, J. 1971, The Astrophysical Journal, 170, 265
43. —. 1975, Monthly Notices of the Royal Astronomical Society
44. Stone, J., Gardiner, T. A., Teuben, P., Hawley, J. F., & Simon, J. B. 2008, The Astrophysical Journal Supplement Series, 178, 137
45. Tautz, R. C., Lerche, I., & Kruse, F. 2013, Astronomy & Astrophysics, 555, A101
46. Tu, C.-Y., & Marsch, E. 1995, Space Sci. Rev., 73, 1
47. Yan, H., & Lazarian, A. 2008, The Astrophysical Journal, 673, 942
48. Yao, S., He, J.-S., Marsch, E., Tu, C.-Y., Pedersen, A., Rème, H., & Trotignon, J. G. 2011, ApJ, 728, 146
104212