The Efficiency of SecondOrder Fermi Acceleration by Weakly Compressible MHD Turbulence
Abstract
We investigate the effects of pitchangle scattering on the efficiency of particle heating and acceleration by MHD turbulence using phenomenological estimates and simulations of nonrelativistic test particles interacting with strong, subsonic MHD turbulence. We include an imposed pitchangle 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 nearEarth 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 quasiresonantly with turbulent fluctuations in magnetic field strength. Scattering increases the longterm efficiency of this resonant heating by factors of a few10, but the distribution function does not develop a significant nonthermal powerlaw tail. For higher scattering rates, the interaction between particles and turbulent fluctuations becomes nonresonant, governed by particles heating and cooling adiabatically as they encounter turbulent density fluctuations. Rapid pitchangle scattering can produce a powerlaw tail in the proton distribution function but this requires finetuning of parameters. Moreover, in the nearEarth solar wind, a significant powerlaw 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 MHDscale turbulent fluctuations are unlikely to be the origin of the tail in the proton distribution function observed in the solar wind.
Subject headings:
plasmas – heating – acceleration of particles – (Sun:) solar wind1. 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 MHDscale turbulent fluctuations. In this paper we address this problem using numerical simulations of test particles interacting with strong MHD turbulence.
According to quasilinear theory, highgyrofrequency particles interacting with lowfrequency, longwavelength 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 longlived 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 highvelocity 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, smallscale 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 ioncyclotron or whistler waves excited independent of velocity space instabilities can have the same effect for superthermal particles.
If pitchangle 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 waveparticle interactions no longer occur because scattering destroys the waveparticle 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 nonresonant 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 nonthermal, 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 pitchangle 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 quasilinear theory and thus provide a more robust estimate of particle energization by MHDscale turbulent fluctuations in heliospheric and astrophysical plasmas. In our calculations, we restrict our analysis to nonrelativistic particles and focus on lowfrequency, weaklycompressible 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 pitchangle 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 pitchangle scattering can produce powerlaw tails in the distribution function. We then summarize our results and discuss their implications (§6).
2. Velocity diffusion: Analytic Estimates
The interaction of highgyrofrequency test particles with lowfrequency 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 pitchangle scattering converts this parallel velocity diffusion into isotropic velocity diffusion, described by
(1) 
where is the 3D distribution function, normalized by (Achterberg, 1981).
In this paper we focus on cases in which the pitchangle 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 pitchangle scattering rate . However, it is quite plausible that in reality the pitchangle scattering rate does depend on velocity, and, in particular, that thermal and superthermal particles have different scattering rates. Because our numerical results in §4 are based on test particle calculations, the case of a velocitydependent 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 pitchangle cosine ,
(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 quasilinear theory, these magnetic compressions are due to fast or slow magnetosonic 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 waveparticle 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 resonancebroadened 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 waveparticle 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
(3) 
where is the phase velocity of the magnetic compressions (fast or slow waves in quasilinear theory) and for a resonancebroadening 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:
(4) 
For high velocities, equation 4 predicts .
In addition to the quasiresonant diffusion described by equation 4, there is a contribution from nonresonant interactions with largescale 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 (subthermal) particles. We term these Fermi TypeB (FTB) interactions. Diffusion coefficient in this limit are calibrated against test particle simulations in Lynn et al. (2012).
2.2. NonResonant 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 nonresonant 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 pressurebalanced entropy modes. Quantitatively, this energy change may be expressed as
(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 nonresonant 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 nonpropagating instabilities such as the firehose and mirror instabilities dominate pitchangle 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 outerscale of the turbulent fluctuations, and . A simplified phenomenological version of Ptuskin’s derivation, as applied to the outerscale eddies, is as follows. The correlation time of the particleturbulence 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 outerscale (drivingscale) eddies. Thus, the resulting velocity diffusion coefficient is
(6) 
which is very similar to equation 21 of Ptuskin (1988), up to
dimensionless coefficients and the identification of as a proportionality constant
times .
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
(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 pitchangle 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 largescale magnetic field, particles undergo a onedimensional 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 particleturbulence 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 largescale 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
(8) 
and the resulting velocity diffusion coefficient is
(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).
NonResonant 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., highfrequency Alfvén waves), and the particles scatter elastically in the scattering frame, nonrelativistic 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
(10) 
This is essentially Fermi scattering (Fermi, 1949) in the nonrelativistic 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 Fermilike 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 pitchangle 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 OrnsteinUhlenbeck 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 nearEarth 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.
Parameter  Value 

Resolution  
Volume ()  
() 
0.1 

1 


() 





() 


Measurement of turbulence properties
The predicted efficiency of nonresonant 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 powerlaw 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 realspace. This yields for our fiducial simulation summarized in Table 1.
To calculate the correlation time, we choose one point per cpu to output a timeseries 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 nearEarth 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 onedimensional 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 pitchangle 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 triangularshaped 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 restframe of the MHD fluid in the simulation. More details are provided in Lehe et al. (2009).
Pitchangle scattering
We implement pitchangle 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 pitchangle scattering rate is in fact given by our input parameter as long as , so that scattering is dominated by small changes in pitchangle. 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, pitchangle 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 highgyrofrequency 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
(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 pitchangle 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 pitchangle 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 leastsquares (as indicated in Fig. 2).
4. Testparticle Diffusion and Heating
4.1. Velocity Diffusion: Stationary Scatterers
Figures 3 and 4 show diffusion coefficients for a broad range of pitchangle 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 outerscale 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 pitchangle average of the parallel diffusion coefficient appropriate for resonance broadened TTD + Fermi TypeB diffusion described in §2.1 (cyan line in Fig. 3). In this limit, scattering isotropizes the velocityspace 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 largescale eddies produce the Fermi TypeB diffusion at low velocities – these interactions are thus fundamentally altered once the scattering rate is comparable to the frequency of the largescale 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 velocityspace 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 nonresonant 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 nonresonant regime of equation 9 requires simultaneously satisfying and , which restricts us to the following range of velocities:
(12) 
For the simulation with , the highvelocity nonresonant 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, nonresonant 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 nonresonant 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 Fermilike scattering predicted by equation 10, namely . At sufficiently high velocity, the nonresonant diffusion shown in Figures 3 and 4 would eventually dominate over the roughly constant diffusion coefficient shown in Figure 6.
4.3. Testparticle heating
Figure 7 shows how the pitchangle 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 pitchangle 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 pitchangle 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 outerscale 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. Powerlaw tails?
Figure 8 shows the latetime 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 pitchangle 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 (moreso with increasing ; see Fig. 7), but no significant nonthermal powerlaw 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 powerlaw component over almost a decade in velocity, particularly for . The nonthermal tail also becomes more energetically important at later times, after the distribution function has been evolved for . There is no evidence for a universal powerlaw slope in the velocity distribution function in this regime of rapid pitchangle scattering. Instead, the exact slope depends sensitively on time, the pitchangle 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 powerlaw at earlier times.
The significant nonthermal 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.
The absence of a nonthermal 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 nonthermal 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 powerlaw distribution function. Our simulations do not explicitly have escape and one might thus worry that the absence of a powerlaw 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 powerlaw 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 wellposed question is whether escape due to spatial diffusion, which is selfconsistently produced by pitchangle scattering, produces the requisite conditions for a powerlaw tail. Scattering implies so that . If is roughly independent of velocity, only a very steep diffusion coefficient of the form can produce a powerlaw tail in the presence of diffusive escape. This is the regime in which we already find a significant nonthermal tail in Figure 9. The velocity diffusion coefficient that results from most of our calculations is (roughly) , i.e., . This applies in both the quasiresonant 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 nonthermal 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 pitchangle scattering in modifying the efficiency of nonrelativistic particle heating and acceleration. To study this we have implemented explicit pitchangle 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 nonpropagating 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 pitchangle scattering may differ for thermal and superthermal particles. Results for a velocitydependent scattering rate can be obtained from the calculations presented here by appropriately interpolating between the diffusion coefficients for different scattering rates shown in Figures 35 (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 pitchangle 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 wavevector 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 pitchangle 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 MHDscale 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 pitchangle 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 quasiresonant 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 waveparticle interaction is often called transittime damping (TTD). Pitch angle scattering increases the efficiency of this quasiresonant 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 pitchangle scattering (; Fig. 7).
In the limit of rapid pitchangle 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, pitchangle scattering produces a Fermilike velocity kick that corresponds to a velocity diffusion coefficient that is independent of velocity for nonrelativistic particles (eq. 10 and Fig. 6).
One of our primary motivations for including the effects of pitchangle scattering was to assess whether pitchangle scattering increases the efficiency of particle acceleration by MHD turbulence, i.e., the generation of a nonthermal tail in the distribution function. We find that this is not the case at low scattering rates where quasiresonant 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 nonthermal tail. This conclusion from our test particle simulations is consistent with solutions of the FokkerPlanck 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 nonthermal tail, but only under certain relatively restrictive conditions. The powerlaw 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 powerlaw 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 pitchangle 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 wellstudied and wellconstrained astrophysical plasma to which our results can be applied? This depends to some extent on the effective scattering rate in the lowcollisionality 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 pitchangle 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 outerscale of turbulence in the nearEarth 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 superthermal particles. These estimates are probably only accurate to a factor of , but this is sufficient to draw several important conclusions.
First, pitchangle scattering likely enhances heating of the thermal particles by the compressible part of the cascade (which is primarily slow modes in the nearEarth solar wind; Yao et al. 2011; Howes et al. 2012) by factors of a few10 (Fig. 7). Secondly, Figure 8 shows that even in the presence of pitchangle scattering with , the particle distribution function does not develop a significant nonthermal tail due to MHDscale 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 pitchangle scattering would produce a powerlaw tail in the proton distribution function in the nearEarth 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 pitchangle 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 nearEarth solar wind, adiabatic cooling (not present in our calculations) dominates particle acceleration by MHD turbulence even for very rapid pitchangle scattering; this would preclude the powerlaw 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 kineticscale physics acting on small scales in the turbulent
cascade could produce such a nonthermal tail (e.g., via reconnection or kinetic Alfvén wave turbulence). In addition, our conclusions are based primarily on the wellstudied plasma conditions near AU. If the pitchangle 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 nonthermal tail.
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 “reacceleration” 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.
Footnotes
 affiliation: Physics Department, University of California, Berkeley, CA 94720; jacob.lynn@berkeley.edu
 affiliation: Astronomy Department and Theoretical Astrophysics Center, University of California, Berkeley, CA 94720
 affiliation: Astronomy Department and Theoretical Astrophysics Center, University of California, Berkeley, CA 94720
 affiliation: Space Science Center and Department of Physics, University of New Hampshire, Durham, NH 03824
 affiliation: Astronomy Department and Theoretical Astrophysics Center, University of California, Berkeley, CA 94720
 affiliation: Present address: Canadian Institute for Theoretical Astrophysics, 60 St. George Street, University of Toronto, Toronto, ON MS H, Canada
 Ptuskin’s notation for spatial and velocity diffusion coefficients is the opposite of ours (we label them as and , respectively).
 The turbulent energy input rate. This corresponds to a sonic Mach number of .
 Ratio of thermal to magnetic pressure.
 refers to the correlation time in the OrnsteinUhlenbeck 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.
 Outer (driving) scale of the turbulence.
 RMS density fluctuations in the saturated turbulence.
 RMS magnetic field fluctuations in the saturated turbulence.
 Test particle gyrofrequency.
 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 powerlaw tail would develop, based on analytic and numerical solutions of the FokkerPlanck equation with .
 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
 Achterberg, A. 1981, Astronomy & Astrophysics, 97, 259
 Bale, S., Kasper, J., Howes, G., Quataert, E., Salem, C., & Sundkvist, D. 2009, Physical Review Letters, 103, 211101
 Bellamy, B. R., Cairns, I. H., & Smith, C. W. 2005, Journal of Geophysical Research (Space Physics), 110, 10104
 Bieber, J. W., Matthaeus, W. H., Smith, C. W., Wanner, W., Kallenrode, M.B., & Wibberenz, G. 1994, The Astrophysical Journal, 420, 294
 Bochsler, P. 2007, A&A Rev., 14, 1
 Boris, J. 1970, in Proceedings of the Fourth Conference on Numerical Simulations of Plasmas, Naval Research Lab, 3–67
 Breech, B., Matthaeus, W. H., Cranmer, S. R., Kasper, J. C., & Oughton, S. 2009, Journal of Geophysical Research (Space Physics), 114, 9103
 Chandran, B. D. G. 2003, The Astrophysical Journal, 599, 1426
 Chandran, B. D. G., Li, B., Rogers, B., Quataert, E., & Germaschewski, K. 2010, The Astrophysical Journal, 720, 503
 Chandran, B. D. G., & Maron, J. L. 2004, The Astrophysical Journal, 603, 23
 Cho, J., & Lazarian, A. 2003, Monthly Notices of the Royal Astronomical Society, 345, 325
 Fermi, E. 1949, Physical Review, 75, 1169
 Fisk, L., & Gloeckler, G. 2006, The Astrophysical Journal Letters, 640, 79
 —. 2008, The Astrophysical Journal, 686, 1466
 Giacalone, J., & Jokipii, J. R. 1994, The Astrophysical Journal, 430, L137
 Gloeckler, G. 2003, AIP Conference Proceedings, 679, 583
 Gloeckler, G., Schwadron, N. A., Fisk, L. A., & Geiss, J. 1995, Geophys. Res. Lett., 22, 2665
 Goldreich, P., & Sridhar, S. 1995, The Astrophysical Journal, 438, 763
 Gruzinov, A., & Quataert, E. 1999, The Astrophysical Journal, 520, 849
 Hellinger, P., & Trávníček, P. M. 2008, Journal of Geophysical Research, 113, A10109
 Hockney, R., & Eastwood, J. 1981, Computer Simulation Using Particles (CRC Press)
 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
 Howes, G. G., Cowley, S. C., Dorland, W., Hammett, G. W., Quataert, E., & Schekochihin, A. A. 2008, Journal of Geophysical Research, 113, 1
 Jokipii, J. R. 1966, The Astrophysical Journal, 146, 480
 Jokipii, J. R., & Lee, M. A. 2010, The Astrophysical Journal, 713, 475
 Kasper, J. C., Lazarus, A. J., & Gary, S. P. 2002, Geophysical Research Letters, 29, 1839
 Kennel, C., & Engelmann, F. 1966, Physics of Fluids, 9, 2377
 Kowal, G., & Lazarian, a. 2010, The Astrophysical Journal, 720, 742
 Lehe, R., Parrish, I. J., & Quataert, E. 2009, The Astrophysical Journal, 707, 404
 Longair, M. S. 1992, High Energy Astrophysics, 2nd edn. (Cambridge University Press)
 Lynn, J. W., Parrish, I. J., Quataert, E., & Chandran, B. D. G. 2012, The Astrophysical Journal, 758, 78
 Mace, R. L., Matthaeus, W. H., & Bieber, J. W. 2000, The Astrophysical Journal, 538, 192
 Maechler, M., et al. 2012, —sfsmisc—: Utilities from Seminar fuer Statistik ETH Zurich, R Package Version 1.020
 Ng, C. S., Bhattacharjee, A., Munsi, D., Isenberg, P. A., & Smith, C. W. 2010, Journal of Geophysical Research (Space Physics), 115, 2101
 Ptuskin, V. 1988, Soviet Astronomy Letters, 14, 255
 Qin, G., & Shalchi, A. 2012, Advances in Space Research, 49, 1643
 R Core Team. 2012, R: A Language and Environment for Statistical Computing, R Foundation for Statistical Computing, Vienna, Austria, ISBN 3900051070
 Schekochihin, A., Cowley, S., Kulsrud, R., Rosin, M., & Heinemann, T. 2008, Physical Review Letters, 100, 081301
 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
 Shalchi, A., Bieber, J. W., Matthaeus, W. H., & Qin, G. 2004, The Astrophysical Journal, 616, 617
 Sharma, P., Quataert, E., Hammett, G. W., & Stone, J. M. 2007, The Astrophysical Journal, 667, 714
 Skilling, J. 1971, The Astrophysical Journal, 170, 265
 —. 1975, Monthly Notices of the Royal Astronomical Society
 Stone, J., Gardiner, T. A., Teuben, P., Hawley, J. F., & Simon, J. B. 2008, The Astrophysical Journal Supplement Series, 178, 137
 Tautz, R. C., Lerche, I., & Kruse, F. 2013, Astronomy & Astrophysics, 555, A101
 Tu, C.Y., & Marsch, E. 1995, Space Sci. Rev., 73, 1
 Yan, H., & Lazarian, A. 2008, The Astrophysical Journal, 673, 942
 Yao, S., He, J.S., Marsch, E., Tu, C.Y., Pedersen, A., Rème, H., & Trotignon, J. G. 2011, ApJ, 728, 146