# Hybrid-Kinetic Simulations of Ion Heating in Alfvénic Turbulence

###### Abstract

We present three-dimensional, hybrid-kinetic numerical simulations of driven Alfvén-wave turbulence of relevance to the collisionless near-Earth solar wind. Special attention is paid to the spectral transition that occurs near the ion-Larmor scale and to the origins of preferential perpendicular ion heating and of non-thermal wings in the parallel distribution function. Several novel diagnostics are used to show that the ion heating rate increases as the kinetic-Alfvén-wave fluctuations, which comprise the majority of the sub-ion-Larmor turbulent cascade, attain near-ion-cyclotron frequencies. We find that – of the cascade energy goes into heating the ions, broadly consistent with the near-Earth solar wind. This heating is accompanied by clear velocity-space signatures in the particle energization rates and the distribution functions, including a flattened core in the perpendicular-velocity distribution and non-Maxwellian wings in the parallel-velocity distribution. The latter are attributed to transit-time damping and the pitch-angle scattering of perpendicularly heated particles into the parallel direction. Accompanying these features is a steepening of the spectral index of sub-ion-Larmor magnetic-field fluctuations beyond the canonical , as field energy is transferred to thermal energy. These predictions may be tested by measurements in the near-Earth solar wind.

###### Subject headings:

…## 1. Introduction.

It has been 48 years since NASA’s Mariner 5 established definitively that the interplanetary medium plays host to a broadband spectrum of large-amplitude Alfvén waves propagating outwards from the Sun (Belcher & Davis 1971, following pioneering work using Mariner 2 data by Coleman 1968). We now know that the solar wind is turbulent (Tu & Marsch, 1995; Goldstein et al., 1995), exhibiting a power spectrum extending over several decades in scale (Bruno & Carbone, 2005; Alexandrova et al., 2013). Most of the energy at large scales is in the form of Alfvénic fluctuations, which have magnetic and velocity fields perpendicular to the mean magnetic-field direction. As this energy cascades down to smaller scales, an inertial range is set up, one whose defining characteristic is anisotropy with respect to the field direction (Bieber et al., 1996; Horbury et al., 2008; Podesta, 2009; Wicks et al., 2010; Chen et al., 2011). This anisotropy, central to modern theories of magnetohydrodynamic (MHD) turbulence (Goldreich & Sridhar, 1995; Boldyrev, 2006; Chandran et al., 2015; Mallet & Schekochihin, 2017) and manifest in the observed shapes of turbulent eddies and their spectral slopes (Horbury et al., 2012), extends all the way down to kinetic scales (Sahraoui et al., 2010; Chen et al., 2010), where the turbulence is ultimately dissipated as heat.

The amplitudes of turbulent fluctuations measured in the solar wind are positively correlated with solar-wind temperature (Grappin et al., 1990) and imply a turbulent heating rate comparable to the observationally inferred solar-wind heating rate (Smith et al., 2001; Breech et al., 2009; Cranmer et al., 2009; Stawarz et al., 2009). While these observations establish a close connection between the global evolution of the solar wind and the dissipation of turbulence within it, the precise nature of this dissipation is puzzling. Minor ions in coronal holes and protons in low-beta fast-solar-wind streams are heated in such a way that thermal motions perpendicular to the background magnetic field are more rapid than thermal motions along it (i.e., ; Kohl et al. 1998; Li et al. 1998; Antonucci et al. 2000; Marsch et al. 1982a, 2004; Hellinger et al. 2006). Moreover, the evolution of proton temperature anisotropy from to clearly indicates non-adiabatic particle heating preferentially in the field-perpendicular direction (see Fig. 1 of Matteini et al. 2007).

These observations present a challenge for models of solar-wind heating based upon theories of Alfvénic turbulence, which predict an anisotropic cascade of energy primarily to small perpendicular (rather than parallel) scales (i.e., ; Shebalin et al. 1983; Goldreich & Sridhar 1995; Ng & Bhattacharjee 1996; Matthaeus et al. 1998; Galtier et al. 2000; Cho et al. 2002; Maron & Goldreich 2001). Such an anisotropic cascade is inefficient at transporting energy to high frequencies traditionally considered necessary to explain the observed strong perpendicular heating (e.g., via ion-cyclotron resonant heating; Quataert 1998; Leamon et al. 1998; Isenberg 2001; Marsch & Tu 2001; Hollweg & Isenberg 2002; Kasper et al. 2013; Cranmer 2014).

An alternative explanation for the measured strong perpendicular heating of ions is that of low-frequency stochastic heating, which arises when particles interact with turbulent fluctuations whose characteristic frequencies are much smaller than the cyclotron frequency but whose amplitudes at the Larmor scale (i.e., for ions) are sufficiently large (McChesney et al., 1987). The particle’s motion in the field-perpendicular plane then becomes chaotic instead of quasi-periodic (Kruskal, 1962), leading to diffusion in the perpendicular energy space due to interactions with the time-varying electrostatic potential (Johnson & Cheng, 2001; Chen et al., 2001; White et al., 2002; Voitenko & Goossens, 2004; Bourouaine et al., 2008; Chandran et al., 2010). In the context of solar-wind turbulence, low-frequency stochastic heating has been studied numerically using test particles in randomly phased kinetic Alfvén waves (KAWs) (Chandran et al., 2010) and in reduced-MHD turbulence (Xia et al., 2013), and observationally in Helios-2 and Wind data (Chandran et al., 2013; Bourouaine & Chandran, 2013; Vech et al., 2017). It was also the focus of work by Vasquez (2015), who used hybrid-kinetic simulations of decaying turbulence in a low-beta, cold-electron plasma to find that the perpendicular heating rate scales with the cube of the turbulence amplitude evaluated at the ion-Larmor scale (as predicted by Chandran et al. 2010, but without their multiplicative factor that exponentially suppresses stochastic heating at small enough turbulence amplitude).

In this paper, we analyze ion heating in three-dimensional (3D) hybrid-kinetic simulations of driven, quasi-steady-state, magnetized turbulence, resolving scales above and below the ion Larmor radius. This analysis is done self-consistently, in that the evolution of the distribution function due to the ions’ interactions with the electromagnetic fields in turn modifies those fields in a reciprocal fashion (as opposed to test-particle calculations, in which the particles’ evolution does not feed back on the fluctuations driving it). Our goal is to understand the relationship between the properties of kinetic, Alfvénic turbulence and the character of the ion heating occurring within it. For this, it is important to note that our simulations do not adopt the oft-employed gyrokinetic approximation (see, e.g., Howes et al., 2006; Schekochihin et al., 2009), and so we allow for electromagnetic fluctuations having finite amplitudes on ion-Larmor scales and/or having ion-Larmor frequencies. This is crucial, as it is quantitatively unclear to what extent the observed anisotropy of solar-wind turbulence precludes an appreciably energetic component of high-frequency electromagnetic fluctuations (e.g., He et al., 2011).

An outline is as follows. In §2 we detail our numerical method and state the parameters used in the runs. §3 catalogues our results, whose interpretation is taken up in §4. We close in §5 with a summary of our results, placed in the context of particle heating in the turbulent solar wind and other hot, dilute astrophysical plasmas.

## 2. Hybrid-kinetic Simulations of Driven Alfvénic Turbulence

We consider a non-relativistic, quasi-neutral, collisionless, and initially homogeneous plasma of kinetic ions (mass , charge ) and massless fluid electrons threaded by a uniform magnetic field and subject to a stochastic driving force . The model equations governing the evolution of the ion distribution function and the magnetic field are, respectively, the Vlasov equation,

(1) |

and Faraday’s law of induction,

(2) |

The electric field,

(3) |

is obtained by expanding the electron momentum equation in , enforcing quasi-neutrality,

(4) |

assuming isothermal electrons (), and using Ampére’s law to solve for the mean electron velocity

(5) |

in terms of the mean ion velocity and the current density . Equations (1)–(5) constitute the hybrid-kinetic model of Vlasov ions and (massless) fluid electrons (Byers et al., 1978; Hewett & Nielson, 1978).

These equations are solved using the second-order–accurate, particle-in-cell code Pegasus (Kunz et al., 2014). A nonlinear method is used to reduce the impact of finite-particle-number noise on computed moments of . To remove small-scale power from the fluctuating fields, the zeroth () and first () moments of are low-pass filtered once per time step. A fourth-order hyper-resistivity is included to remove grid-scale magnetic energy, its value tuned to achieve steady state.

To excite Alfvénically polarized fluctuations and minimize the excitation of compressible motions, is oriented in the - plane perpendicular (“”) to and constrained to satisfy . At each simulation time step, the Fourier coefficients are independently generated from a Gaussian-random field with power spectrum in the wavenumber range and , where is the size of the periodic computational domain in the () direction. The resulting force is then inverse-Fourier transformed, shifted to ensure no net momentum injection, and normalized to provide power per unit volume . The force is time-correlated over using an Ornstein–Uhlenbeck process:

(6) |

which has the autocorrelation (in the limit )

(7) |

where is the timestep, , is the correlation time of the driving, and is a new Gaussian-random field generated as detailed above. Time-correlated driving avoids spurious particle acceleration via resonances with high-frequency power in, e.g., -correlated driving (Lynn et al., 2012).

Most of the power in strong MHD turbulence resides in fluctuations satisfying (Goldreich & Sridhar, 1995; Mallet et al., 2015), where () is the wavenumber parallel (perpendicular) to the local magnetic field and is the Alfvén speed. We therefore choose and such that, in saturation, the rms velocity fluctuation satisfies at the outer scale; likewise, . This is meant to mimic an energy cascade from larger scales that is present in the solar wind, whose inertial-range turbulence is consistent with critical balance (Horbury et al., 2008; Podesta, 2009; Luo & Wu, 2010; Wicks et al., 2010).

By replacing electron kinetics with an isothermal equation of state, the hybrid-kinetic model excludes electron Landau damping and its effect on the turbulent cascade (e.g., TenBarge et al., 2013; Told et al., 2016a; Grošelj et al., 2017). However, the hybrid-kinetic model affords a huge cost savings over the fully kinetic approach (see Vasquez et al. (2014), Cerri et al. (2017), and Franci et al. (2018) for recent examples of using the hybrid-kinetic approach to simulate 3D solar-wind-like turbulence, and Parashar et al. (2009) for an early hybrid-kinetic approach to studying particle heating in a 2D Orszag–Tang vortex). And it captures physics not described by the oft-employed gyrokinetic approach (e.g., Howes et al., 2006, 2008b, 2011; Schekochihin et al., 2009; Told et al., 2015; Kawazura et al., 2019), such as stochastic ion heating, ion-cyclotron resonances, and modes whose propagation angles are not asymptotically oblique. We refer the reader to Told et al. (2016b) and Camporeale & Burgess (2017) for comparison of linear modes in hybrid-kinetics, gyrokinetics, and full kinetics.

We present results from two simulations: and , both with (the subscript “0” denotes an initial value). For , particles per cell were drawn from a Maxwell distribution and placed on a 3D periodic grid of cells spanning , where is the ion thermal speed, is the ion Larmor radius, and is the ion gyrofrequency. These parameters provide reasonable scale separation between the grid scale, the ion-kinetic scales, and the driving scales: and , where is the ion skin depth. For , , , and . These imply and .

With these scales borne in mind, it is useful to compare with the observed spectral anisotropy in the near-Earth solar wind. There, spectral anisotropy of the turbulent cascade begins around and increases as towards smaller scales (Wicks et al., 2010). Adopting these scalings implies around , where our simulated inertial range begins. Thus, our chosen box aspect ratios of (for ) and (for ) accurately capture the spectral anisotropy of solar-wind turbulence arriving from larger scales to near ion-Larmor scales. (This point is revisited near the end of §3.1 in the context of the observed spectral anisotropy at ion-Larmor scales.)

Both simulations required at least to obtain steady-state turbulence. An additional were run to procure statistically converged results (viz., for and for ). In what follows, denotes a spatio-temporal average over all cells measured in steady state.

## 3. Results

An example of the turbulent steady state is given in Figure 1, which shows pseudo-color images of the - and -components (i.e., those perpendicular to the mean field) of the fluctuating magnetic field from the run. Spatial anisotropy is evident, with short perpendicular scales and long parallel scales consistent with a critically balanced cascade (i.e., , the box aspect ratio). This anisotropy plays an important role in shaping the energy spectra (§3.1) and the nature of ion heating (§3.2).

### 3.1. Energy spectra and spectral anisotropy

Figures 2(a,b) present energy spectra of magnetic-field () and electric-field () fluctuations for and versus the wavenumber perpendicular to . Their scale-dependent spectral indices , computed about each using 21 neighboring points (except for the first 5 points for which we use 11 neighboring harmonics), are shown in Figures 2(c,d). Accompanying these in Figures 2(e,f) are the ion-flow-velocity (), density (), and parallel-magnetic-field () energy spectra.

In the inertial range (), the spectral slopes are close to , the spectral slope predicted for an anisotropic, critically balanced cascade of Alfvén waves (Goldreich & Sridhar, 1995). A spectral break occurs at , at which point flattens to take on a slope near , in agreement with measurements that show a flattening of the solar-wind electric-field power spectrum in ion-kinetic range (Bale et al., 2005; Sahraoui et al., 2009; Salem et al., 2012) and predictions for intermittent KAW turbulence (Boldyrev & Perez, 2012). In this range, also steepens to take on a slope comparable to the commonly observed in the near-Earth solar wind (e.g., Sahraoui et al., 2009; Alexandrova et al., 2009, 2012) and within the range found in Cluster spacecraft measurements of the solar wind (Sahraoui et al., 2013). (For , steepens further, an attribute we revisit below when discussing ion heating.) Despite being steeper than predictions for standard (; Schekochihin et al. 2009) and intermittent (; Boldyrev & Perez 2012) KAW turbulence, the normalized perturbed density approximately follows the linear KAW eigenfunction (for ; Schekochihin et al. 2009; Boldyrev & Perez 2013). This suggests that the sub-ion-Larmor-scale cascade is primarily composed of KAWs (at least up to where steepens), in agreement with combined analyses of magnetic-field fluctuations and of electric-field or density fluctuations in the ion-kinetic range of solar-wind turbulence (e.g., Salem et al., 2012; Chen et al., 2013). That being said, there is slightly more (about a factor ) magnetic-fluctuation energy than anticipated for KAWs, suggesting the presence of additional wave modes (an excess of was measured in the solar wind by Chen et al. (2013)). At , the spectra steepen further due to hyper-resistivity and spectral filters (ion heating is less important at these scales—see §3.2).

The predicted slopes of in the inertial range and (or ) in the KAW range follow from assuming locality of interactions and constant energy flux in Fourier space, combined with a model for the spatial anisotropy of the turbulent fluctuations, . The latter is afforded by the critical balance assumption, which states that the scale-dependent nonlinear cascade time is comparable to the linear timescale of the dominant wave at that scale (Goldreich & Sridhar, 1995)—essentially a causality argument (Boldyrev, 2005). For the inertial range, in which the characteristic linear frequency is that of Alfvén waves, , a perpendicular spectrum of corresponds to . For the KAW range, in which the characteristic linear frequency is , a perpendicular spectrum of () corresponds to (), the steeper parenthetical values reflecting the two-dimensionalization of intermittent KAW turbulence (Boldyrev & Perez, 2012).

In Figure 2(g,h), we show the spectral anisotropy computed from our simulations.^{1}^{1}1To compute the spectral anisotropy, we use equation (34) of Cho & Lazarian (2009) to compute the scale-dependent wavenumber parallel to the local mean magnetic field, ; those authors showed that this method works very well for steep spectra. We also computed two- and three-point second-order structure functions and measured the spectral anisotropy from their isocontours (similar to Cho & Vishniac 2000), finding similar results. In both runs, in the inertial range. At the start of the KAW range, the wavevector anisotropy and for and , respectively, corresponding to wavevector obliquities and . These values are comparable to those measured in the solar wind (e.g., Sahraoui et al., 2010; Narita et al., 2011). Well into the sub-ion-Larmor range, however, we observe , steeper than current theoretical predictions and corresponding to a scale-independent anisotropy. Constant spectral anisotropy in the sub-ion-Larmor-scale range has also been seen in other hybrid-kinetic simulations of 3D Alfvénic turbulence (Franci et al., 2018).

Constant spectral anisotropy at implies that the turbulence there is more likely to attain ion-Larmor frequencies before reaching electron scales, since , rather than the standard predicted by theories of low-frequency gyrokinetic turbulence (e.g., Schekochihin et al., 2009). Such high-frequency fluctuations facilitate additional energy-transfer channels that are not present in gyrokinetics (cf. Howes et al., 2008a), such as cyclotron resonances and high-frequency stochastic heating, to which we now turn.

### 3.2. Ion heating

Figure 3 shows parallel () and perpendicular () particle energization rates as functions of . These are calculated by spatially filtering the electric field into 12 logarithmically-spaced -bins, giving , and then summing and over all particles for each bin (the first and last bins are not shown). Rates are normalized to the rate of energy injection by the large-scale forcing, .^{2}^{2}2We are careful to distinguish between particle energization, by which we mean the rate at which work is done on the particles temporally averaged over long time intervals, and particle heating, by which we mean an increase in the Maxwellian temperature of the distribution function. For example, transit-time damping (TTD) is driven by but results in an increase in the parallel temperature , and changes in need not be driven by particle energization (e.g., pitch-angle scattering of perpendicular particle energy into the parallel direction). In both runs, . For , about 80% of injected energy goes into the ions; the remaining 20% is removed by hyper-resistivity. For , this ratio is roughly . In both cases, peaks at , approximately where the measured spectral anisotropies in the KAW range imply (marked by the dotted lines in Figures 2(g,h) and 3).^{3}^{3}3These values were calculated using the approximation and confirmed using the numerical linear hybrid-kinetic solver HYDROS (Told et al., 2016a, cf. their figure 8). This suggests that the majority of the energization is produced when ions move in the oscillating potential of high-frequency KAWs, a possibility discussed further in §4. For now, we note that the sub-ion-Larmor spectrum of magnetic-field fluctuations seen in Fig. 2(a,b) steepens beyond the anticipated at the values of for which the perpendicular energization is largest—likely not a coincidence.

In Figure 4, we present the velocity- and real-space tracks of two particles from the run, one (“1”) that exhibits appreciable perpendicular energization and another (“2”) that does not. Upper panels show snapshots of small-scale () magnetic (left) and electric (right) field fluctuations in the - plane located at the -coordinate of both tracked particles. The black lines trace projected particle trajectories over 4 gyro-periods. Particle 1 resides in a region of large-amplitude, sheet-like structures, whereas particle 2 samples only weak field fluctuations. (Note that the fluctuations experience the same drift motion as the particles, and so particle 1’s guiding center is not actually sweeping across the small-scale fluctuations.) Lower panels show the evolution of the particles’ peculiar (“thermal”) velocity over the final of the run. The gray shaded region marks the time interval over which particle trajectories are plotted. At each point, thermal velocities are averaged over to remove high-frequency oscillations in energy; these oscillations are shown for with the thin green line.

Particle 1 sees an appreciable increase in its energy, with over only 4 gyro-periods. Most of this increase is perpendicular. On the other hand, the energy of particle 2 stays nearly constant. Note further that the parallel and perpendicular energies do not grow monotonically. Often, they are subject to strong kicks during which the total energy is almost constant but the pitch angle of the particle changes dramatically, a feature we refer to as pitch-angle scattering. As a result, high- particles scatter and subsequently contribute to wings produced in the parallel distribution function, a feature we return to below in the context of Figures 6–9.

Figure 5 shows the parallel (perpendicular) differential energization, , in the gyrotropic velocity space . (Namely, is the particle energization per interval of velocity space; these rates are normalized in the figure to the total energization rate in each run.) Figure 5 is to be read alongside Figure 6, which displays the parallel and perpendicular ion distribution functions measured at the ends of both runs, viz., and , respectively. The dashed lines in Figure 6 show best-fit Maxwellians to the core of (with temperature ) and to the tail of (with temperature ). Both runs exhibit the following attributes: (i) resonant features in the particle-energization rates with fine-scale structure near and (Figs 5(a),(c)); (ii) quasi-linear flattening of the parallel distribution at and non-thermal wings at (Figs 6(a),(c)); (iii) almost no change in the parallel temperature of the core (; Figs 6(a),(c)), with very little parallel energization there (Figs 5(a),(c)); (iv) a broadened perpendicular distribution for (Figs 6(b),(d)), where the perpendicular energization peaks (Figs 5(b),(d)); and (v) a flattening of the core of the perpendicular distribution (Figs 6(b),(d)), with suppressed perpendicular energization at (Figs 5(b),(d)). Regarding this final point, we find evidence early in each run for perpendicular energization at , which ultimately causes the observed flattening the core of (see Fig. 8 and accompanying discussion below).

All of these features can also be seen in the full 2D (gyrotropic) distribution function shown in Figure 7, where the differences between the two runs are even more striking. In particular, non-Maxwellian features such as the flattened core and the parallel beams at are readily apparent, with the latter driving for and .

We associate with these non-Maxwellian features a number of physical effects. The early-time flattening of for is most likely due to stochastic heating by low-frequency, Larmor-scale fluctuations, following the prediction by Klein & Chandran (2016) that such heating leads to a flattop distribution (see also Johnson & Cheng 2001). As stochastic heating is expected to be more important at lower for fixed (e.g., Chandran et al., 2010; Hoppock et al., 2018), it is noteworthy that exhibits a larger amount of flattening for than for . The broadening of the thermal tail of is instead driven by the peak in at , which constitutes the bulk of the perpendicular heating and is related to the peak in at – (Fig. 3).

This interpretation is further strengthened by examining the box-averaged perpendicular-energy diffusion coefficient

(8) |

as a function of , shown in Figure 8 both at early times in the run (panel a) and at late times in both runs (panel b). Accompanying these panels are (c,d) the differential perpendicular particle energization (as in Figure 5) and (e,f) the perpendicular distribution function (as in Figure 6) at those times. In the latter panels, the evolution of from its initial condition (black dotted line) to its profile at the beginning (solid line) and the end (dashed line) of the time-averaging window are shown. At early times, the diffusion coefficient and the accompanying exhibit two distinct peaks. We associate the peak at with low-frequency stochastic heating, which flattens the core of by accelerating particles to larger . The peak in is associated with a diffusion coefficient that scales approximately as . At late times in the run (red lines), after the core of the perpendicular distribution has been flattened, the diffusion coefficient is roughly constant with and very little energization happens for (note that is not well-defined in this limit, as both and are close to zero). At larger velocities, the perpendicular energization peaks (as in Figure 5) and the diffusion coefficient continues to scale approximately as . In the run, which did not experience as great of a flattening in the core of , appears to be a good approximation for .

Most of the injected energy goes into the increase of perpendicular temperature and into the development of non-thermal tails in the parallel distribution function. (Overall, of the total energization ultimately finds its way into the non-thermal wings, with the remaining going into raising the perpendicular temperature.) Given that most of the energization is perpendicular, this strongly suggests that transit-time damping (TTD) and pitch-angle scattering of super-thermal into are the mechanisms responsible for the non-Maxwellian features seen in .

We tested these possibilities by examining the energization and pitch angles of 160,000 individually tracked particles in the run. We found that the mirror force (where is the magnetic moment) is responsible for of , a value consistent with the quasi-linear flattening of the distribution function observed at in Figures 6 and 7. The remaining (small) amount of increase in can be accounted for by Landau damping (i.e., ). The other of is due to some other mechanism that leads to perpendicular heating at (see §4 for a discussion of possibilities). To examine the effect of pitch-angle scattering on the distribution function, we divided each particle track into time intervals of and computed at the beginning and the end of each interval. We then filtered each interval based on whether or not changes (up or down) by a certain minimal factor that we call “threshold”. Large values of threshold represent large changes in a particle’s pitch angle; small values of threshold include time intervals during which the pitch angle is almost constant. Figure 9 shows , the average^{4}^{4}4To compute as a function of threshold, we sum and for all events that change by a factor larger than threshold, and then divide this energy increment by the total time over which we examine the particle tracks () and by the total number of tracked particles. rate of change of the parallel (blue), perpendicular (red), and total (green) thermal energies of tracked particles, segregated into super- and sub-thermal populations ( and , respectively). Given a threshold value, Figure 9 provides the expected energization rates for an individual particle due to all events whose change in are above that threshold. As the length of the time interval used to sub-divide the particle tracks is somewhat arbitrary, we indicate with the shaded regions the variation of particle-energization rates for time intervals from to . Note that is a decreasing function of threshold, as it represents the cumulative energization for all events above the threshold. The energization per event (not shown) is an increasing function of threshold. For , there is a net conversion of thermal energy from parallel to perpendicular, consistent with the other diagnostics shown in Figures 5–8. For , however, the flow of thermal energy between parallel and perpendicular is reversed and, for thresholds , this flow takes place at constant total energy. This strongly suggests that pitch-angle scattering is responsible for this transfer of super-thermal perpendicular energy into super-thermal parallel energy, ultimately producing the non-thermal wings seen in Figures 6 and 7.^{5}^{5}5Evidence of pitch-angle scattering can also be seen in the evolution of particle 1 at – in Figure 4; and there is the additional circumstantial evidence that the more pronounced non-thermal wings in seen in the run (as compared to the run) coincides with a larger perpendicular temperature measured in the tail of .

We re-emphasize that, while at the end of our simulations, more than half of the cascade energy ultimately goes into the development of the non-thermal wings in the parallel distribution function. To illustrate that, the blue line in Figure 10 traces the temporal evolution of the ratio of the total parallel and perpendicular energies of the particles, , from the run in which the non-thermal wings are most pronounced (the factor of accounts for the number of degrees of freedom in the perpendicular direction, so that an isotropic Maxwellian has ). The accompanying red line represents the ratio of best-fit Maxwellian temperatures to the core of the parallel distribution function () and to the tail of the perpendicular distribution function (). The initial drop in is caused by the increase in perpendicular bulk motion in the initial stage of the simulation. (Recall that denotes the total (bulk + thermal) velocity of the ion particles.) Once the turbulence obtains steady state (for ), the perpendicular temperature steadily grows larger than the parallel temperature of the core. This is mirrored in the evolution of , at least until , after which the ratio of energies suddenly begins to increase, eventually becoming larger than 1 at . The distinction between the ratio of energies and the ratio of best-fit-Maxwellian temperatures is thus an important one. Indeed, while non-thermal -wings are often observed in ion velocity distribution functions in the solar wind, the “observed” is often determined by a bi-Maxwellian fit to the core of the distribution (e.g., Bame et al., 1975; Marsch et al., 1982b).

## 4. Interpretation of Sub-Ion-Larmor-Scale Perpendicular Heating

The perpendicular energization that occurs at early times for appears to be consistent with the predictions of stochastic ion heating by low-frequency fluctuations. Namely, the core of the perpendicular distribution function becomes flattened as particles there are promoted to higher perpendicular energies via energization, with lower and larger amplitudes leading to more flattening. This finding is notable, in that it was achieved in a self-consistent setting in which the evolution of the distribution function is allowed to feed back on the electromagnetic fields that drive that evolution. This is in contrast to other studies of stochastic heating that have employed a test-particle approach.

By contrast, the origin of the perpendicular energization that occurs for is less clear. Despite the highly suggestive correlation between the at which (i) the perpendicular ion heating peaks, (ii) the sub-ion-Larmor magnetic-energy spectrum steepens beyond the canonical , and (iii) the linear KAW frequency attains the ion-Larmor frequency, it is nevertheless difficult for us to definitively conclude that the majority of the ion heating is due to cyclotron heating by high-frequency KAW turbulence. The reason for this reluctance is twofold. First, while the scale separation and dynamic range afforded by our simulations bears resemblance to solar-wind data on the side, on the side there is a less-than-optimal scale separation between the grid scale (near which hyper-resistivity and low-pass filtering are important) and the location of maximal perpendicular energization and the concomitant spectral steepening. But more important than this computational concern is a physical one: we currently have no theory explaining the dependencies of the perpendicular particle-energization rate and the perpendicular energy-diffusion coefficient on the particles’ perpendicular energy for . Both of these profiles appear to be inconsistent with previous studies of stochastic heating (Klein & Chandran, 2016) and cyclotron damping (e.g., Isenberg, 2004; Isenberg & Vasquez, 2007).

That being said, it is notable that Stereo measurements in high-speed solar-wind streams show a rapid decrease in the power anisotropy (i.e., the difference between the energy stored in the perpendicular and parallel fluctuations at a given scale) near , which Podesta (2009) associated with the strong linear dissipation of KAWs occurring at when (parameters not unlike ours). It is in this wavenumber range that KAWs are known to couple to ion-Bernstein waves (IBWs; Bernstein 1958), which Podesta (2012) showed provide a channel for mode coupling and thus energy transfer from KAWs. Once excited, IBWs are strongly damped through a combination of ion-cyclotron and electron-Landau resonances (the latter of which are of course absent in our simulations). This is rather suggestive, but more work is needed to predict the velocity-space signatures of ion energization by high-frequency KAWs/IBWs.

## 5. Summary: Implications for the Solar Wind and for Simulations of Kinetic Turbulence

We have presented 3D numerical simulations of driven, quasi-steady-state, hybrid-kinetic turbulence in a magnetized plasma of relevance to the solar wind. Despite the more general hybrid-kinetic framework employed, many aspects of the spectral scalings are in rough agreement with those obtained using gyrokinetics (Howes et al., 2008b, 2011; Told et al., 2015). These include a critically balanced, spatially anisotropic, inertial-range cascade of Alfvénic fluctuations; a spectral break occurring near the ion-Larmor scale, beyond which the magnetic spectrum steepens and the electric spectrum flattens; signatures of linear phase mixing in the ion distribution function (e.g., flattening of the parallel distribution function near ); and what appears to be a sub-ion-Larmor-scale cascade composed primarily of KAWs. However, there are differences, most notably in the efficiency and mechanism of ion heating.

Although these simulations were originally designed to test the theory of stochastic ion heating occurring at —which we do find—our results also make the case for perpendicular ion heating at sub-ion-Larmor scales, as a cascade of KAWs approaches the ion-cyclotron frequency. This heating, alongside contributions from TTD, Landau damping, and pitch-angle scattering, simultaneously heat the particles perpendicularly, produce non-Maxwellian wings in the ion parallel distribution function, and steepen the sub-ion-Larmor-range magnetic-energy spectra beyond the typically observed power-law scaling by transferring electromagnetic energy into thermal energy. It is perhaps no coincidence, then, that the ion-kinetic-range spectral index as measured in the solar wind correlates with the amount of inferred energy dissipation, with more dissipation going hand-in-hand with steeper spectra (e.g., Smith et al., 2006). We predict that such spectra will be found to be accompanied by non-Maxwellian wings in , a preponderance of perpendicular heating (over parallel heating), and a flattened core in .

Given that complementary gyrokinetic theory and simulations show a predominance of electron heating over ion heating for (Howes, 2010; Howes et al., 2011; Told et al., 2015; Navarro et al., 2016; Kawazura et al., 2019), it is worthwhile to contemplate which framework is better suited for describing near-Earth solar-wind turbulence (see Howes et al. (2008a, §3) for arguments in favor of a gyrokinetic description). Gyrokinetic theory is built on the assumption of asymptotically low-frequency, small-amplitude, spatially anisotropic fluctuations; adiabatic invariance is enforced by the resulting equations. The accompanying computational savings affords a magneto-kinetic description of a realistic hydrogenic plasma without prohibitive cost. By contrast, the hybrid model makes no such simplifying assumptions, but saves computational expense by neglecting electron kinetics. While the latter precludes a truly rigorous study of ion versus electron heating, it is worth noting that there is empirical evidence for a majority fraction (60%) of ion versus electron heating between and in the solar wind (Cranmer et al., 2009); recall that – of our cascade energy goes into heating the ions (the rest is dissipated via hyper-resistivity). Moreover, the ion-Larmor-scale spectral anisotropy in the solar wind is not necessarily asymptotically small. Indeed, – for fluctuations measured in the near-Earth solar wind (Sahraoui et al., 2010); our values are similar, –. If scales with to some power larger than the adopted by Howes et al. (2008a) for , as it does in our simulations and as is predicted in theories of intermittent KAW turbulence (Boldyrev & Perez, 2012), it becomes all the more probable that the ion-cyclotron frequency will be attained in a KAW cascade. These considerations favor the hybrid-kinetic description over the gyrokinetic one.

More broadly, our study of ion heating in kinetic, Alfvénic turbulence may be applicable to a number of collisionless astrophysical plasmas in which the collisional mean free path is comparable to or even larger than the system size, the canonical example of which being the low-luminosity accretion flow onto the supermassive black hole at the Galactic center, Sgr A. The observed low luminosity of this system can be explained if the gravitational energy released during accretion is stored primarily in the poorly radiating ions rather than the electrons (Ichimaru, 1977; Narayan & Yi, 1994; Narayan et al., 1998). As angular-momentum transport in myriad accretion disks is thought to be driven by magnetorotational turbulence (Balbus & Hawley, 1998), the question of ion versus electron heating in such turbulent flows is thus important for models of low-luminosity accretion (Quataert & Gruzinov, 1999; Sharma et al., 2007). Recently, there have been several studies implementing various particle-heating prescriptions in general relativistic MHD simulations of black-hole accretion flows (e.g., Ressler et al., 2015; Sa̧dowski et al., 2017; Chael et al., 2018). These “sub-grid” prescriptions are based either on gyrokinetic theory and simulation (Howes, 2010; Told et al., 2015), which predict preferential electron (ion) heating at low (high) plasma beta driven by Landau-resonant damping, or on models of collisionless reconnection (Rowan et al., 2017; Werner et al., 2018; Rowan et al., 2019), in which the amount of ion versus electron heating is sensitive to the presence of a strong guide field. It would be interesting to study the effect of the non-Landau-resonant processes presented in this paper, which predominantly heat the ions, on the imaging and evolution of collisionless accretion flows. However, the application of our results to these systems is not straightforward. While the wavevector anisotropies in our simulations are consistent with the scale separation observed in the solar wind (viz., a factor of between the outer scale and the ion-Larmor scale), the scale separation in black-hole accretion flows is expected to be even larger (e.g., a factor of between the outer scale and the ion-Larmor scale for Sgr A; see Quataert 1998). This implies that Alfvén-wave/KAW frequencies near the ion-Larmor scale are a factor smaller than in the solar wind, possibly inhibiting particle heating via cyclotron resonances. We defer a study of ion heating in this regime to future work.

In a subsequent publication, a wider parameter study will be conducted alongside further analysis of field-particle correlations (following Klein & Howes 2016, Howes et al. 2017, and Klein et al. 2017). In the meantime, we hope that the various particle energization diagnostics employed herein will spur their application to both existing and future simulation data of solar-wind-like kinetic turbulence.

## References

- Alexandrova et al. (2013) Alexandrova, O., Chen, C. H. K., Sorriso-Valvo, L., Horbury, T. S., & Bale, S. D. 2013, SSRv, 178, 101
- Alexandrova et al. (2012) Alexandrova, O., Lacombe, C., Mangeney, A., Grappin, R., & Maksimovic, M. 2012, ApJ, 760, 121
- Alexandrova et al. (2009) Alexandrova, O., Saur, J., Lacombe, C., et al. 2009, PhRvL, 103, 165003
- Antonucci et al. (2000) Antonucci, E., Dodero, M. A., & Giordano, S. 2000, SoPh, 197, 115
- Balbus & Hawley (1998) Balbus, S. A., & Hawley, J. F. 1998, Reviews of Modern Physics, 70, 1
- Bale et al. (2005) Bale, S. D., Kellogg, P. J., Mozer, F. S., Horbury, T. S., & Reme, H. 2005, PhRvL, 94, 215002
- Bame et al. (1975) Bame, S. J., Asbridge, J. R., Feldman, W. C., Montgomery, M. D., & Gary, S. P. 1975, GeoRL, 2, 373
- Belcher & Davis (1971) Belcher, J. W., & Davis, Jr., L. 1971, JGR, 76, 3534
- Bernstein (1958) Bernstein, I. B. 1958, PhRv, 109, 10
- Bieber et al. (1996) Bieber, J. W., Wanner, W., & Matthaeus, W. H. 1996, JGR, 101, 2511
- Boldyrev (2005) Boldyrev, S. 2005, ApJ, 626, L37
- Boldyrev (2006) —. 2006, PhRvL, 96, 115002
- Boldyrev & Perez (2012) Boldyrev, S., & Perez, J. C. 2012, ApJ, 758, L44
- Boldyrev & Perez (2013) Boldyrev, S., & Perez, J. C. 2013, in Astronomical Society of the Pacific Conference Series, Vol. 474, Numerical Modeling of Space Plasma Flows (ASTRONUM2012), ed. N. V. Pogorelov, E. Audit, & G. P. Zank, 3
- Bourouaine & Chandran (2013) Bourouaine, S., & Chandran, B. D. G. 2013, ApJ, 774, 96
- Bourouaine et al. (2008) Bourouaine, S., Marsch, E., & Vocks, C. 2008, ApJ, 684, L119
- Breech et al. (2009) Breech, B., Matthaeus, W. H., Cranmer, S. R., Kasper, J. C., & Oughton, S. 2009, JGR, 114, 9103
- Bruno & Carbone (2005) Bruno, R., & Carbone, V. 2005, LRSP, 2, 4
- Byers et al. (1978) Byers, J. A., Cohen, B. I., Condit, W. C., & Hanson, J. D. 1978, JCoPh, 27, 363
- Camporeale & Burgess (2017) Camporeale, E., & Burgess, D. 2017, JPlPh, 83, 535830201
- Cerri et al. (2017) Cerri, S. S., Servidio, S., & Califano, F. 2017, ApJ, 846, L18
- Chael et al. (2018) Chael, A., Rowan, M., Narayan, R., Johnson, M., & Sironi, L. 2018, MNRAS, 478, 5209
- Chandran et al. (2010) Chandran, B. D. G., Li, B., Rogers, B. N., Quataert, E., & Germaschewski, K. 2010, ApJ, 720, 503
- Chandran et al. (2015) Chandran, B. D. G., Schekochihin, A. A., & Mallet, A. 2015, ApJ, 807, 39
- Chandran et al. (2013) Chandran, B. D. G., Verscharen, D., Quataert, E., et al. 2013, ApJ, 776, 45
- Chen et al. (2013) Chen, C. H. K., Boldyrev, S., Xia, Q., & Perez, J. C. 2013, PhRvL, 110, 225002
- Chen et al. (2010) Chen, C. H. K., Horbury, T. S., Schekochihin, A. A., et al. 2010, PhRvL, 104, 255002
- Chen et al. (2011) Chen, C. H. K., Mallet, A., Yousef, T. A., Schekochihin, A. A., & Horbury, T. S. 2011, MNRAS, 415, 3219
- Chen et al. (2001) Chen, L., Lin, Z., & White, R. 2001, PhPl, 8, 4713
- Cho & Lazarian (2009) Cho, J., & Lazarian, A. 2009, ApJ, 701, 236
- Cho et al. (2002) Cho, J., Lazarian, A., & Vishniac, E. T. 2002, ApJ, 564, 291
- Cho & Vishniac (2000) Cho, J., & Vishniac, E. T. 2000, ApJ, 539, 273
- Coleman (1968) Coleman, Jr., P. J. 1968, ApJ, 153, 371
- Cranmer (2014) Cranmer, S. R. 2014, ApJS, 213, 16
- Cranmer et al. (2009) Cranmer, S. R., Matthaeus, W. H., Breech, B. A., & Kasper, J. C. 2009, ApJ, 702, 1604
- Franci et al. (2018) Franci, L., Landi, S., Verdini, A., Matteini, L., & Hellinger, P. 2018, ApJ, 853, 26
- Galtier et al. (2000) Galtier, S., Nazarenko, S. V., Newell, A. C., & Pouquet, A. 2000, JPlPh, 63, 447
- Goldreich & Sridhar (1995) Goldreich, P., & Sridhar, S. 1995, ApJ, 438, 763
- Goldstein et al. (1995) Goldstein, M. L., Roberts, D. A., & Matthaeus, W. H. 1995, ARA&A, 33, 283
- Grappin et al. (1990) Grappin, R., Mangeney, A., & Marsch, E. 1990, JGR, 95, 8197
- Grošelj et al. (2017) Grošelj, D., Cerri, S. S., Bañón Navarro, A., et al. 2017, ApJ, 847, 28
- He et al. (2011) He, J., Marsch, E., Tu, C., Yao, S., & Tian, H. 2011, ApJ, 731, 85
- Hellinger et al. (2006) Hellinger, P., Trávníček, P., Kasper, J. C., & Lazarus, A. J. 2006, GeoRL, 33, 9101
- Hewett & Nielson (1978) Hewett, D. W., & Nielson, C. W. 1978, JCoPh, 29, 219
- Hollweg & Isenberg (2002) Hollweg, J. V., & Isenberg, P. A. 2002, JGR, 107, 1147
- Hoppock et al. (2018) Hoppock, I. W., Chandran, B. D. G., Klein, K. G., Mallet, A., & Verscharen, D. 2018, Journal of Plasma Physics, 84, 905840615
- Horbury et al. (2008) Horbury, T. S., Forman, M., & Oughton, S. 2008, PhRvL, 101, 175005
- Horbury et al. (2012) Horbury, T. S., Wicks, R. T., & Chen, C. H. K. 2012, SSRv, 172, 325
- Howes (2010) Howes, G. G. 2010, MNRAS, 409, L104
- Howes et al. (2006) Howes, G. G., Cowley, S. C., Dorland, W., et al. 2006, ApJ, 651, 590
- Howes et al. (2008a) —. 2008a, JGRA, 113, A05103
- Howes et al. (2008b) Howes, G. G., Dorland, W., Cowley, S. C., et al. 2008b, PhRvL, 100, 065004
- Howes et al. (2017) Howes, G. G., Klein, K. G., & Li, T. C. 2017, JPlPh, 83, 705830102
- Howes et al. (2011) Howes, G. G., Tenbarge, J. M., Dorland, W., et al. 2011, PhRvL, 107, 035004
- Ichimaru (1977) Ichimaru, S. 1977, ApJ, 214, 840
- Isenberg (2001) Isenberg, P. A. 2001, SSRv, 95, 119
- Isenberg (2004) —. 2004, JGRA, 109, A03101
- Isenberg & Vasquez (2007) Isenberg, P. A., & Vasquez, B. J. 2007, ApJ, 668, 546
- Johnson & Cheng (2001) Johnson, J. R., & Cheng, C. Z. 2001, GeoRL, 28, 4421
- Kasper et al. (2013) Kasper, J. C., Maruca, B. A., Stevens, M. L., & Zaslavsky, A. 2013, PhRvL, 110, 091102
- Kawazura et al. (2019) Kawazura, Y., Barnes, M., & Schekochihin, A. A. 2019, PNAS, 116, arXiv:1807.07702
- Klein & Chandran (2016) Klein, K. G., & Chandran, B. D. G. 2016, ApJ, 820, 47
- Klein & Howes (2016) Klein, K. G., & Howes, G. G. 2016, ApJ, 826, L30
- Klein et al. (2017) Klein, K. G., Howes, G. G., & Tenbarge, J. M. 2017, JPlPh, 83, 535830401
- Kohl et al. (1998) Kohl, J. L., Noci, G., Antonucci, E., et al. 1998, ApJ, 501, L127
- Kruskal (1962) Kruskal, M. 1962, JMP, 3, 806
- Kunz et al. (2014) Kunz, M. W., Stone, J. M., & Bai, X.-N. 2014, JCoPh, 259, 154
- Leamon et al. (1998) Leamon, R. J., Smith, C. W., Ness, N. F., Matthaeus, W. H., & Wong, H. K. 1998, JGR, 103, 4775
- Li et al. (1998) Li, X., Habbal, S. R., Kohl, J. L., & Noci, G. 1998, ApJ, 501, L133
- Luo & Wu (2010) Luo, Q. Y., & Wu, D. J. 2010, ApJ, 714, L138
- Lynn et al. (2012) Lynn, J. W., Parrish, I. J., Quataert, E., & Chandran, B. D. G. 2012, ApJ, 758, 78
- Mallet & Schekochihin (2017) Mallet, A., & Schekochihin, A. A. 2017, MNRAS, 466, 3918
- Mallet et al. (2015) Mallet, A., Schekochihin, A. A., & Chandran, B. D. G. 2015, MNRAS, 449, L77
- Maron & Goldreich (2001) Maron, J., & Goldreich, P. 2001, ApJ, 554, 1175
- Marsch et al. (2004) Marsch, E., Ao, X.-Z., & Tu, C.-Y. 2004, JGR, 109, 4102
- Marsch et al. (1982a) Marsch, E., Rosenbauer, H., Schwenn, R., Muehlhaeuser, K.-H., & Neubauer, F. M. 1982a, JGR, 87, 35
- Marsch et al. (1982b) Marsch, E., Schwenn, R., Rosenbauer, H., et al. 1982b, JGR, 87, 52
- Marsch & Tu (2001) Marsch, E., & Tu, C.-Y. 2001, JGR, 106, 8357
- Matteini et al. (2007) Matteini, L., Landi, S., Hellinger, P., et al. 2007, GeoRL, 34, 20105
- Matthaeus et al. (1998) Matthaeus, W. H., Oughton, S., Ghosh, S., & Hossain, M. 1998, PhRvL, 81, 2056
- McChesney et al. (1987) McChesney, J. M., Stern, R. A., & Bellan, P. M. 1987, PhRvL, 59, 1436
- Narayan et al. (1998) Narayan, R., Mahadevan, R., Grindlay, J. E., Popham, R. G., & Gammie, C. 1998, ApJ, 492, 554
- Narayan & Yi (1994) Narayan, R., & Yi, I. 1994, ApJ, 428, L13
- Narita et al. (2011) Narita, Y., Gary, S. P., Saito, S., Glassmeier, K.-H., & Motschmann, U. 2011, GeoRL, 38, L05101
- Navarro et al. (2016) Navarro, A. B., Teaca, B., Told, D., et al. 2016, PhRvL, 117, 245101
- Ng & Bhattacharjee (1996) Ng, C. S., & Bhattacharjee, A. 1996, ApJ, 465, 845
- Parashar et al. (2009) Parashar, T. N., Shay, M. A., Cassak, P. A., & Matthaeus, W. H. 2009, PhPl, 16, 032310
- Podesta (2009) Podesta, J. J. 2009, ApJ, 698, 986
- Podesta (2012) —. 2012, JGRA, 117, A07101
- Quataert (1998) Quataert, E. 1998, ApJ, 500, 978
- Quataert & Gruzinov (1999) Quataert, E., & Gruzinov, A. 1999, ApJ, 520, 248
- Ressler et al. (2015) Ressler, S. M., Tchekhovskoy, A., Quataert, E., Chandra, M., & Gammie, C. F. 2015, MNRAS, 454, 1848
- Rowan et al. (2017) Rowan, M. E., Sironi, L., & Narayan, R. 2017, ApJ, 850, 29
- Rowan et al. (2019) —. 2019, arXiv e-prints, arXiv:1901.05438
- Sahraoui et al. (2010) Sahraoui, F., Goldstein, M. L., Belmont, G., Canu, P., & Rezeau, L. 2010, PhRvL, 105, 131101
- Sahraoui et al. (2009) Sahraoui, F., Goldstein, M. L., Robert, P., & Khotyaintsev, Y. V. 2009, PhRvL, 102, 231102
- Sahraoui et al. (2013) Sahraoui, F., Huang, S. Y., Belmont, G., et al. 2013, ApJ, 777, 15
- Salem et al. (2012) Salem, C. S., Howes, G. G., Sundkvist, D., et al. 2012, ApJ, 745, L9
- Sa̧dowski et al. (2017) Sa̧dowski, A., Wielgus, M., Narayan, R., et al. 2017, MNRAS, 466, 705
- Schekochihin et al. (2009) Schekochihin, A. A., Cowley, S. C., Dorland, W., et al. 2009, ApJS, 182, 310
- Sharma et al. (2007) Sharma, P., Quataert, E., Hammett, G. W., & Stone, J. M. 2007, ApJ, 667, 714
- Shebalin et al. (1983) Shebalin, J. V., Matthaeus, W. H., & Montgomery, D. 1983, JPlPh, 29, 525
- Smith et al. (2006) Smith, C. W., Hamilton, K., Vasquez, B. J., & Leamon, R. J. 2006, ApJ, 645, L85
- Smith et al. (2001) Smith, C. W., Matthaeus, W. H., Zank, G. P., et al. 2001, JGR, 106, 8253
- Stawarz et al. (2009) Stawarz, J. E., Smith, C. W., Vasquez, B. J., Forman, M. A., & MacBride, B. T. 2009, ApJ, 697, 1119
- TenBarge et al. (2013) TenBarge, J. M., Howes, G. G., & Dorland, W. 2013, ApJ, 774, 139
- Told et al. (2016a) Told, D., Cookmeyer, J., Astfalk, P., & Jenko, F. 2016a, NJPh, 18, 075001
- Told et al. (2016b) Told, D., Cookmeyer, J., Muller, F., Astfalk, P., & Jenko, F. 2016b, NJPh, 18, 065011
- Told et al. (2015) Told, D., Jenko, F., TenBarge, J. M., Howes, G. G., & Hammett, G. W. 2015, PhRvL, 115, 025003
- Tu & Marsch (1995) Tu, C.-Y., & Marsch, E. 1995, SSRv, 73, 1
- Vasquez (2015) Vasquez, B. J. 2015, ApJ, 806, 33
- Vasquez et al. (2014) Vasquez, B. J., Markovskii, S. A., & Chandran, B. D. G. 2014, ApJ, 788, 178
- Vech et al. (2017) Vech, D., Klein, K. G., & Kasper, J. C. 2017, ApJ, 850, L11
- Voitenko & Goossens (2004) Voitenko, Y., & Goossens, M. 2004, ApJ, 605, L149
- Werner et al. (2018) Werner, G. R., Uzdensky, D. A., Begelman, M. C., Cerutti, B., & Nalewajko, K. 2018, MNRAS, 473, 4840
- White et al. (2002) White, R., Chen, L., & Lin, Z. 2002, PhPl, 9, 1890
- Wicks et al. (2010) Wicks, R. T., Horbury, T. S., Chen, C. H. K., & Schekochihin, A. A. 2010, MNRAS, 407, L31
- Xia et al. (2013) Xia, Q., Perez, J. C., Chandran, B. D. G., & Quataert, E. 2013, ApJ, 776, 90