<figure><img src="image_1.png" /></figure>

# On Stochastic Heating and Its Phase-space Signatures in Low-beta Kinetic Turbulence

S. S. Cerri¹, L. Arzamasskiy¹,², and M. W. Kunz¹,³

¹Department of Astrophysical Sciences, Princeton University, 4 Ivy Lane, Princeton, NJ 08544, USA; scerri@astro.princeton.edu

²Institute for Advanced Study, 1 Einstein Drive, Princeton, NJ 08540, USA

³Princeton Plasma Physics Laboratory, PO Box 451, Princeton, NJ 08543, USA

Received 2021 February 18; revised 2021 April 8; accepted 2021 April 25; published 2021 August 5

## Abstract

We revisit the theory of stochastic heating of ions and investigate its phase-space signatures in kinetic turbulence of relevance to low-$\beta$ portions of the solar wind. In particular, we retain a full scale-dependent approach in our treatment, and we explicitly consider the case in which electric-field fluctuations can be described by a generalized Ohm's law that includes Hall and thermoelectric effects. These two electric-field terms provide the dominant contributions to stochastic ion heating when the ion-Larmor scale is much smaller than the ion skin depth, $\rho_i \ll d_i$, which is the case at $\beta \ll 1$. Employing well-known spectral scaling laws for Alfvén-wave and kinetic-Alfvén-wave turbulent fluctuations, we obtain scaling relations characterizing the field-perpendicular particle-energization rate and energy diffusion coefficient associated with stochastic heating in these two regimes. Phase-space signatures of ion heating are then investigated using three-dimensional hybrid-kinetic simulations of continuously driven Alfvénic turbulence at low $\beta$ (namely, $\beta_i = \beta_e = 0.3$ and $\beta_i = \beta_e = 1/9$). In these simulations, energization of ions parallel to the magnetic field is subdominant compared to its perpendicular counterpart ($Q_{\parallel,i} \ll Q_{\perp,i}$), and the fraction of turbulent energy that goes into ion heating is $\approx 75\%$ at $\beta_i = 0.3$ and $\approx 40\%$ at $\beta_i \simeq 0.1$. The phase-space signatures of ion energization are consistent with Landau-resonant collisionless damping and a ($\beta$-dependent) combination of ion-cyclotron and stochastic heating. We demonstrate good agreement between our scale-dependent theory and various signatures associated with the stochastic portion of the heating. We discuss briefly the effect of intermittency on stochastic heating and the implications of our work for the interpretation of stochastic heating in solar-wind spacecraft data.

*Unified Astronomy Thesaurus concepts: Plasma astrophysics (1261); Plasma physics (2089); Space plasmas (1544); Solar wind (1534)*

## 1. Introduction

The solar wind is arguably the most well-diagnosed weakly collisional magnetized plasma, both in terms of the electromagnetic fluctuations it hosts and the thermodynamics of its constituent particles. It therefore serves as an excellent (and, with some effort, directly accessible) laboratory with which one may discriminate between different theories of magnetized turbulence and the various ways in which such turbulence energizes plasma particles. Indeed, a persistent puzzle in solar-wind research is why the temperature of the solar wind evolves nonadiabatically as it expands, and why this heating occurs preferentially in the direction perpendicular to the local magnetic field (e.g., Marsch et al. 1982; Matteini et al. 2007; Hellinger et al. 2011; Maruca et al. 2011). While the solution to this puzzle is known to be connected to the pervasive Alfvénic turbulence that is now routinely measured by in situ spacecraft (e.g., Goldstein et al. 1995; Alexandrova et al. 2013; Bruno & Carbone 2013; Chen et al. 2020; Sahraoui et al. 2020), the relative contributions to this turbulent heating from different wave-particle interactions are debated.

Much of this debate has been centered on the nature of the turbulent fluctuations and their relative energetic importance at various stages during their nonlinear cascade to increasingly finer scales in both configuration and velocity space (e.g., Leamon et al. 1999; Howes et al. 2008; Schekochihin et al. 2009; Chandran et al. 2011; Cranmer 2014). Namely, how spatially anisotropic are typical fluctuations at a given scale? What fraction of those fluctuations ultimately attain cyclotron frequencies? Are the fluctuations at Larmor scales of sufficient amplitude to disrupt the particles' otherwise smooth gyromotion and heat the plasma

appreciably? How do the answers to these questions depend on the plasma properties, such as the ratio of thermal and magnetic pressures, $\beta \doteq 8\pi p/B^2$? This is an indirect way of understanding particle energization in the solar wind; guided by observational constraints (e.g., Horbury et al. 2012; Chen 2016), one postulates the characteristics of the fluctuations in the turbulent cascade, models the various particle-energization channels available to those fluctuations, and then infers whether those channels are thermodynamically important by comparing the implied heating and any unique features with the data. Such an approach has been used to find evidence for ion-cyclotron-resonant heating in the solar wind via measured correlations between plasma heating, differential flow between ion species, and magnetic-field-biased temperature anisotropy (Kasper et al. 2013). Similarly, correlations between the amplitudes of ion-Larmor-scale magnetic fluctuations and enhanced proton and minor-ion temperatures measured in coronal holes and the bulk solar wind have been taken as evidence for the stochastic heating of ions by low-frequency Alfvén-wave (AW) and kinetic-Alfvén-wave (KAW) fluctuations (Chandran 2010; Bourouaine & Chandran 2013; Chandran et al. 2013; Vech et al. 2017; Martinović et al. 2019, 2020).

A more direct but more technically challenging way of distinguishing between different particle-energization mechanisms is through their imprint on the velocity-space structure of the plasma (e.g., Klein & Howes 2016; Howes et al. 2017; Howes 2017; Klein et al. 2017; Adkins & Schekochihin 2018; Servidio et al. 2017; Cerri et al. 2018; Pezzi et al. 2018; Kawazura et al. 2019; Li et al. 2019). For example, it is well known that collisionless Landau damping flattens the particle distribution
function in the vicinity of “Landau resonances,” at which a particle’s velocity (in a magnetized plasma, the velocity component parallel to the local magnetic-field direction) matches the phase speed of a wave. This flattening is a consequence of the secular transfer of free energy from the electromagnetic waves to the particles, whether it be via parallel electric fields (Landau 1946) or parallel gradients in magnetic-field strength (Barnes 1966). Recently, a clear signature of this transfer (in this case, to the electron population) has been found in data taken in the Earth’s turbulent magnetosheath (Chen et al. 2019). This follows on pioneering work by Marsch & Tu (2001; see also Heuer & Marsch 2007 and He et al. 2015) showing plateaus in solar-wind particle distribution functions near the Alfvén speed, suggesting velocity-space diffusion due to Alfvén/ion-cyclotron fluctuations (e.g., Isenberg 2001; Isenberg & Vasquez 2019). Similar velocity-space signatures of ion-cyclotron damping, revealed by applying field-particle correlation techniques to hybrid-kinetic simulations, have been discussed by Klein et al. (2020).

Nonresonant energization mechanisms, such as stochastic heating, also make an imprint on the velocity space. Adopting the theory of Chandran et al. (2010), Klein & Chandran (2016) showed that the stochastic heating of ions by moderate-amplitude, Larmor-scale, electric-field fluctuations ultimately flattens the core of their velocity distribution function along the field-perpendicular direction. Such a flat-top distribution has been observed recently by Martinović et al. (2020) using data from the Parker Solar Probe. Formulating and testing such velocity-space diagnostics is particularly important in the case of stochastic heating, since it provides an attractive alternative to other (namely, resonant) mechanisms of particle energization whose phase-space signatures have long drawn the attention of the heliophysics community. This becomes particularly true for situations in which the turbulent cascade exhibits strong spatial anisotropy that inhibits the production of high-frequency waves and/or for values of $\beta \ll 1$ at which ions are unable to obtain the Landau resonance (Quataert 1998; Hollweg 1999).

Accordingly, the purpose of this paper is to further elucidate the consequences of stochastic ion heating for the organization of phase space and sharpen certain aspects of how the theory of stochastic heating can be tested using solar-wind data. The paper is written in two parts. First, we extend the work of Chandran et al. (2010) and Klein & Chandran (2016) to make further predictions for the phase-space signatures of stochastic heating and their dependence on the properties of the plasma ($\beta$, ion-to-electron temperature ratio) and turbulence (Section 2). Second, we present results from a new hybrid-kinetic simulation of driven Alfvénic turbulence, which we use to test these predictions (Section 3). We also demonstrate that intermittency, as revealed in the statistics of the electrostatic potential, enhances stochastic heating, with some particles acquiring large amounts of energy in spatially and temporally localized events. A corollary of our analysis is that an oft-employed conversion of measured ion-Larmor-scale magnetic-field fluctuation amplitudes to bulk ion-velocity fluctuations, which are then used in a formula to determine the expected amount of stochastic heating, becomes increasingly inaccurate at low values of $\beta$, precisely where stochastic heating is expected to be most important (Section 4). For $\beta \ll 1$, noninductive components of the electric field—namely, the Hall effect and the thermoelectric field—contribute appreciably to the total electrostatic potential with which the particles interact.

Our work follows on that of Arzamasskiy et al. (2019). Those authors presented results from hybrid-kinetic simulations of driven Alfvénic turbulence and employed several novel diagnostics to quantify the roles of Landau and Barnes damping, stochastic heating, and cyclotron heating—all of which appeared to be in play—in the energization and differential heating of plasma particles at $\beta \lesssim 1$. Taken together, this set of simulations and their analyses suggest that stochastic heating plays an important role in modifying both the velocity distribution function of the ions and the cascade of turbulent energy to sub-ion-Larmor scales in low-$\beta$, collisionless plasmas.

## 2. Theory of Stochastic Ion Heating in AW/KAW Turbulence

Chandran et al. (2010) presented a theory for perpendicular ion heating in the solar wind caused by finite-amplitude, low-frequency, AW/KAW fluctuations occurring on scales comparable to the ion-Larmor scale (following on work by Chen et al. 2001; Johnson & Cheng 2001; White et al. 2002; Voitenko & Goossens 2004; Bourouaine et al. 2008). In this theory, the ions interact stochastically with a time-varying electrostatic potential, break their magnetic moments, and execute a random walk in perpendicular energy. Here we generalize this theory to account for a spectrum of critically balanced fluctuations whose electrostatic potential satisfies a generalized Ohm’s law. We compute the perpendicular heating rate and energy diffusion coefficient as functions of the perpendicular plasma beta parameter of the ions, $\beta_{\perp i} \doteq 8\pi p_{\perp i}/B^2$, which is the ratio of thermal pressure of the ions perpendicular to the magnetic-field direction, $p_{\perp i} \doteq nT_{\perp i}$, where $n$ is the ion number density, and the magnetic pressure, $B^2/8\pi$; the electron-to-ion temperature ratio, $\tau_{\perp} \doteq Z_i T_c / T_{\perp i}$, where $Z_i$ is the ion charge in units of $e$; and the energy cascade rate, $\epsilon$. (We take the electron temperature $T_c$ to be isotropic, for reasons that will be explained in Section 2.2.) Before doing so, we briefly recapitulate the theory presented in Chandran et al. (2010) in a way that establishes the notation used in the remainder of the paper.

### 2.1. Stochastic Heating Revisited

Consider an ion with mass $m_i$ and charge $q_i = Z_i e$ that is interacting with electric-field fluctuations $\delta E_{\perp,\lambda}$ having perpendicular wavelength $\lambda$ of the order of the ion’s gyroradius $\rho_i \doteq w_{\perp}/\Omega_i$, i.e., $k_{\perp} \rho_i \sim 1$. Here $w_{\perp}$ is the component of the ion’s random velocity perpendicular to a background magnetic field $B_0$, $\Omega_i \doteq q_i B_0 / m_i c$ is the ion-cyclotron frequency, and $k_{\perp} = 2\pi/\lambda$ is the field-perpendicular wavenumber associated with $\lambda$. If the amplitude of these fluctuations is sufficiently large (just how large is quantified in Section 2.1.2), the ion’s gyromotion about $B_0$ becomes chaotic, its magnetic moment $\mu \doteq m_i w_{\perp}^2 / 2B$ is no longer conserved, and the ion is stochastically heated in the field-perpendicular direction. Such stochasticity is the result of a sequence of “random kicks” that the ion experiences due to the fluctuating field within a turbulent eddy of size $\lambda \sim \rho_i$.

In what follows, we assume that the main contribution to this heating is from the potential part of the fluctuating electric field, so that $\delta E_{\perp,\lambda} \sim k_{\perp} \delta \Phi_{\lambda}$. This is justified (and verified a posteriori using our simulations) if $\beta_{\perp i}$ is not much larger than unity and/or if the fluctuations’ frequency $\omega$ remains smaller than $\sim \Omega_i / \beta_{\perp i}$ (Hoppock et al. 2018). Such electrostatic fluctuations on the scale of an ion’s gyroradius induce a change
in an ion's perpendicular kinetic energy, $\Delta K_{\perp}$, that is directly related to the average change of the potential over the time $\tau_{\lambda}$ that the particle spends within the turbulent eddy of size $\lambda$, viz., $\Delta K_{\perp} \sim q_i(\partial\delta\Phi_{\lambda}/\partial t)\tau_{\lambda}$. We estimate $\tau_{\lambda}$ as the time required for the ion's guiding center to drift in the direction perpendicular to $B_0$ by a distance of order $\lambda$. Taking this drift to be of the $B \times \nabla\Phi$ type, so that $u_{\text{dr},\lambda} \sim (c/B_0)(|\delta\Phi_{\lambda}|/\lambda)$, we find that

$$ \tau_{\lambda} \sim \Omega_{i}^{-1} \left( \frac{\lambda}{\rho_{\text{th},i}} \right)^{2} \left( \frac{m_{i} v_{\text{th},i}^{2}}{q_{i} |\delta\Phi_{\lambda}|} \right), \quad (1) $$

where $v_{\text{th},i}^2 \doteq 2T_{\perp,i}/m_i$ is square of the (perpendicular) ion-thermal speed and $\rho_{\text{th},i} \doteq v_{\text{th},i}/\Omega_i$ is the thermal ion-Larmor radius. For the change in perpendicular kinetic energy to be effective, the turbulent fluctuations must be as coherent as possible over this timescale. Denoting the typical frequency of the turbulent fluctuations at scale $\lambda$ by $\omega_{\lambda}$, this requirement may be written as $\omega_{\lambda}\tau_{\lambda} \sim 1$. In this case, $\Delta K_{\perp} \sim q_i\omega_{\lambda}\delta\Phi_{\lambda}\tau_{\lambda} \sim q_i\delta\Phi_{\lambda}$. (For a lengthier discussion of these arguments, see Section 2 and Equations (12)–(16) and (24), in particular, of Chandran et al. (2010).)

Using this information and assuming that the stochastic gain of the perpendicular kinetic energy of a single ion during the time $\tau_{\lambda}$ can be seen as a random walk in perpendicular-energy space, we determine the perpendicular-energy diffusion coefficient and heating rate as follows.

### 2.1.1. Perpendicular Diffusion Coefficient and Heating Rate

We quantify the stochastic gain in an ion's perpendicular kinetic energy using the diffusion coefficient $D_{\perp}^E \sim \Delta K_{\perp}^2 / \tau_{\lambda}$. With $\Delta K_{\perp} \sim q_i \delta\Phi_{\lambda}$ and $\tau_{\lambda}$ being given by Equation (1), we find

$$ D_{\perp}^{E}(\lambda) \sim \Omega_{i} \left( \frac{\rho_{\text{th},i}}{\lambda} \right)^{2} \frac{q_{i}^{3} |\delta\Phi_{\lambda}|^{3}}{m_{i} v_{\text{th},i}^{2}}. \quad (2a) $$

Alternatively, $D_{\perp}^E$ may be expressed in velocity space by using the condition $k_{\perp}\rho_i = k_{\perp}w_{\perp}/\Omega_i \sim 1$ to replace $\lambda$ with $(w_{\perp}/v_{\text{th},i})\rho_{\text{th},i}$. Then, denoting the resulting velocity-space potential $\delta\Phi_{\lambda}|_{\lambda \sim w_{\perp}/\Omega_i}$ as $\delta\Phi_w$, Equation 2(a) may be reinterpreted as

$$ D_{\perp}^{E}(w_{\perp}) \sim \Omega_{i} \frac{q_{i}^{3} |\delta\Phi_{w}|^{3}}{m_{i} w_{\perp}^{2}}. \quad (2b) $$

This equation states that particles drawn from different regions of the perpendicular distribution function experience different perpendicular energization, depending on the part of the spectrum of the fluctuations that they sample during their orbits and off of which they stochastically diffuse.

To obtain an equation for how this diffusion affects the evolution of the perpendicular-energy distribution function, $f^E$, we insert Equation 2(b) into the Fokker–Planck-like equation

$$ \frac{\partial f^E}{\partial t} = \frac{\partial}{\partial e_{\perp}} \left( D_{\perp}^E \frac{\partial f^E}{\partial e_{\perp}} \right), \quad (3) $$

where $e_{\perp} \doteq w_{\perp}^2/2$ is the ion's perpendicular kinetic energy per unit mass. Then, using Equation (3), we may write the total perpendicular heating as

$$ Q_{\perp} = - \int d e_{\perp} D_{\perp}^{E} \frac{\partial f^{E}}{\partial e_{\perp}}. \quad (4) $$

Alternatively, one may introduce a differential heating rate in $w_{\perp}$ via<fnref n="4" />

$$ \frac{\partial Q_{\perp}}{\partial w_{\perp}} \doteq - D_{\perp}^{E}(w_{\perp}) \frac{\partial f^{E}(w_{\perp})}{\partial w_{\perp}}, \quad (5) $$

with $D_{\perp}^E(w_{\perp})$ given by Equation 2(b). Equation (5) will be used in Section 3 to compute $D_{\perp}^E(w_{\perp})$ using the functions $f^E(w_{\perp})$ and $\partial Q_{\perp}/\partial w_{\perp}$ obtained directly from our numerical simulations.

It is helpful at this stage to work through a simple estimate for how $D_{\perp}^E$ and $\partial Q_{\perp}/\partial w_{\perp}$ would scale with $w_{\perp}$ for a particular scaling law of the fluctuating potential. Let us assume that the dominant contribution to the electric field is due to $\delta u_i \times B/c$ induction from a fluctuating ion-velocity field $\delta u_i$, such that $\delta\Phi_{\lambda} \sim \lambda\delta u_{\perp,i}(B_0/c)$. Adopting the Kolmogorov-like scaling $\delta u_{\perp,i} \propto \lambda^{1/3}$ for these fluctuations, we find that $\delta\Phi_{\lambda} \propto \lambda^{4/3}$. Enacting the transformation to velocity space described above, $\delta\Phi_w \propto (w_{\perp}/v_{\text{th},i})^{4/3}$. Equation 2(b) then gives $D_{\perp}^E \propto (w_{\perp}/v_{\text{th},i})^2$, which is a scaling that matches the one of Klein & Chandran (2016) when the induction term, $\sim u_{\perp,i}\lambda B_0/c$, is the dominant contribution to the electrostatic potential.<fnref n="5" /> Further assuming a Maxwellian distribution in $w_{\perp}$ yields a differential heating rate, $\partial Q_{\perp}/\partial w_{\perp} \propto (w_{\perp}/v_{\text{th},i})^3 \exp(-w_{\perp}^2/v_{\text{th},i}^2)$. In this case, ion particles whose perpendicular velocities satisfy $v_{\perp}^2 = (3/2)v_{\text{th},i}^2$ would experience the largest differential heating rate.

### 2.1.2. Exponential Suppression of Stochastic Heating

In order to take into account the reduction of stochastic heating due to the near-conservation of the particles' magnetic moments when the fluctuations' amplitudes at the $\rho_{\text{th},i}$ scale are “sufficiently small,” Chandran et al. (2010) proposed a multiplicative exponential suppression term of the type $\exp(-c_2/\xi_{\text{th}})$ in Equation 2(a), where $c_2$ is a (small, scale-independent) constant. This quasi-conservation condition is quantified by a so-called stochasticity parameter $\xi$, which in our theory would read as a scale-dependent parameter defined by<fnref n="6" />

$$ \xi_w \doteq \frac{q_i |\delta\Phi_w|}{m_i w_{\perp}^2}. \quad (6) $$

The parameter $\xi_{\text{th}}$, which is $\xi_w$ evaluated at the ion-thermal speed $w_{\perp} \sim v_{\text{th},i}$ (or, equivalently, at the ion-thermal gyroradius, $\lambda \sim \rho_{\text{th},i}$), provides an estimate of the amount of energy in

<footnote n="4">⁴ This definition is consistent with the diagnostics implemented in our simulations (see Section 3). Vasquez et al. (2020) argued for an alternative definition of $\partial Q_{\perp}/\partial w_{\perp}$, which nevertheless results in the same total heating rate given by Equation (4). Further discussion of this alternative definition and its use in analyzing our simulation results is provided in Appendix A.</footnote>

<footnote n="5">⁵ Note that Klein & Chandran (2016) adopt $\delta u_{\perp,i} \propto \lambda^{1/4}$, consistent with the dynamic-alignment argument of Boldyrev (2006). Then $\delta\Phi_{\lambda} \propto \lambda\delta u_{\perp,i} \propto \lambda^{5/4}$, and Equation 2(b) gives $D_{\perp}^E \propto w_{\perp}^{7/4}$, consistent with Equation (17) of Klein & Chandran (2016).</footnote>

<footnote n="6">⁶ When the induction term provides the dominant contribution to the electrostatic fluctuations, and using the condition $\lambda \sim (w_{\perp}/v_{\text{th},i})\rho_{\text{th},i}$ to obtain $\delta\Phi_w$, our definition of $\xi$ reduces to (a scale-dependent version of) the definition $\xi = \delta u_{\perp}/w_{\perp}$ of Chandran et al. (2010). In that work, this parameter (evaluated at the ion-thermal Larmor scale) is called $\epsilon$. However, in order to avoid confusion with the symbol typically used for the cascade rate, as well as to differentiate the generalized stochasticity parameter based on potential fluctuations from that based on ion flow velocity fluctuations, we use $\xi$ instead. When the need arises to refer specifically to Chandran et al.'s stochasticity parameter (namely, in Section 4), we adopt the notation $\epsilon_i$.</footnote>
the electrostatic potential fluctuations that goes into stochastic heating, weighted by the particles' thermal energy, viz., $\xi_{th} \sim q_i|\delta\Phi_{th}|/m_i v_{th,i}^2$, where $\delta\Phi_{th}$ is the velocity-space potential $\delta\Phi_w$ evaluated at $w_{\perp} \sim v_{th,i}$. An exponential suppression factor would be justified if $\xi_{th} \ll c_2$. One may then obtain a rough estimate for when the amplitude of the potential fluctuations is “sufficiently large” for stochastic heating to be important, that is, when the (thermal-)Larmor-scale potential satisfies $q_i|\delta\Phi_{th}|/m_i v_{th,i}^2 \gtrsim c_2$. An assortment of test-particle calculations (Chandran et al. 2010; Xia et al. 2013) has suggested values for $c_2$ in the range of $\approx 0.1-0.3$. Analyses of solar-wind data in the context of stochastic heating have adopted similar values of $c_2$ (Chandran 2010; Bourouaine & Chandran 2013; Martinović et al. 2019, 2020).<fnref n="7" />

In our theory, we allow for an analogous, scale-dependent exponential suppression term, so that Equation 2(b) becomes (after using Equation (6) to replace $q_i|\delta\Phi_w|$ with $m_i w_{\perp}^2 \xi_w$)

$$ \frac{D_{\perp}^{E}(w_{\perp})}{\Omega_{i} m_{i}^{2} v_{\text{th},i}^{4}} \sim \left(\frac{w_{\perp}}{v_{\text{th},i}}\right)^{4} \xi_{w}^{3} \exp\left(-\frac{c_{*}}{\xi_{w}}\right), \quad (7) $$

where $c_*$ is a constant to be determined. The notation $c_*$ differs from the notation $c_2$ used by Chandran et al. (2010) to emphasize that the exponential correction is being applied within the scale-dependent formulation of $D_{\perp}^{E}(w_{\perp})$, rather than within the scale-independent formulation with $w_{\perp} \approx v_{th,i}$ (or, equivalently, $\lambda \approx \rho_{th,i}$; cf. Equations (20)–(25) of Chandran et al. 2010). For this reason, the value of $c_*$ does not necessarily match that of $c_2$ found in previous work.<fnref n="8" /> Within this scale-dependent formulation, a potential fluctuation is “sufficiently large” to perpendicularly heat an ion with velocity $w_{\perp}$ effectively when its amplitude satisfies $(c/B_0)|\delta\Phi_w|$ $\gtrsim c_* w_{\perp}^2/\Omega_i$. In terms of perpendicular scales $\lambda$, this corresponds to the range for which $(c/B_0)|\delta\Phi_{\lambda}|/\lambda^2 \gtrsim c_* \Omega_i$. We further caution that this “constant” may be dependent upon $\beta$ and/or the level of intermittency in the ion-Larmor-scale fluctuations, the two possibly being related to each other as $\beta$ decreases (e.g., Cerri et al. 2017b; Grošelj et al. 2017). Such intermittency could indeed partially compensate for the simultaneous decrease of $\xi_{th}$ that would be associated with the enhanced separation between injection and $\rho_{th,i}$ scales in the $\beta \ll 1$ regime, which is precisely the regime in which stochastic heating is likely to be most relevant. This possibility seems to be supported by our simulation results (see Section 3.3); future kinetic simulations with yet larger scale separations, and thus statistically smaller values of $\xi$, than those performed here are needed to further investigate the behavior of this exponential correction.

It is worth noting that, while the exponential suppression factor was originally introduced to account for the reduction in perpendicular heating when ion-Larmor-scale fluctuations are small, this factor also serves to suppress stochastic heating by

larger-scale fluctuations (despite their larger relative amplitudes). Qualitatively, the lower frequencies of these fluctuations allow the ions to drift smoothly in a quasi-static potential, precluding chaotic motion and preserving approximate adiabatic invariance. Quantitatively, we may rewrite the argument of the exponential term in Equation (7) as $-c_* \Omega_i \tau_w$, where $\tau_w$ is given by Equation (1) with $\lambda/\rho_{th,i} \sim w_{\perp}/v_{th,i}$. Then the requirement for strong suppression of stochastic heating becomes $\omega_w/\Omega_i \ll c_* \omega_w \tau_w \lesssim c_*$, where $\omega_w$ is the frequency of gyroscale fluctuations as seen by particles with a gyroradius $\rho_i = w_{\perp}/\Omega_i \sim \lambda$. Conversely, fluctuations whose frequencies satisfy $\omega_w/\Omega_i \gtrsim c_*$ are the most effective at stochastically heating the ions.

## 2.2. Generalized Ohm's Law and Contributions to Stochastic Ion Heating

While the example given at the end of Section 2.1.1 is illustrative, the $\delta u_i \times B/c$ inductive electric field contributes just one piece to a more general Ohm's law. In particular, because the mechanism of stochastic ion heating occurs primarily at ion-kinetic scales (which are much smaller than the injection scales), contributions to the electric field from, e.g., the Hall effect may be important, particularly at low values of $\beta_i$ at which the ion skin depth $d_i \gg \rho_{th,i}$. To quantify these contributions, we adopt the following generalized Ohm's law for the electric field $E$ in which electron-inertia effects have been neglected but contributions from the Hall and thermo-electric fields are retained:

$$ E = - \frac{u_i \times B}{c} + \frac{J \times B}{enc} - \frac{\nabla p_e}{en}. \quad (8) $$

Here we have used quasi-neutrality to replace the electron number density $n_e$ with the ion number density $n$. Equation (8) is valid at scales $\lambda$ much larger than the electron kinetic scales, viz., $\lambda \gg d_e$, $\rho_{th,e}$, where $d_e$ and $\rho_{th,e}$ are the electron skin depth and thermal Larmor radius, respectively.<fnref n="9" /> To simplify matters further, we adopt an isothermal equation of state for the electrons, so that the electron pressure $p_e = nT_e$ with $T_e = \text{const}$. This is a good approximation for KAW fluctuations at perpendicular scales satisfying $\rho_{th,e} \ll \lambda \ll \rho_{th,i}$, for which the electron response is Boltzmann and therefore isothermal (see, e.g., Section 7.2 of Schekochihin et al. 2009).

To obtain the potential contribution to the electric field (Equation (8)), we consider AW/KAW turbulence in which the fluctuations are anisotropic with respect to the magnetic-field direction, with $k_{\parallel} \ll k_{\perp}$. As in Section 2.1, we therefore assume that the electric field is dominated by its potential contribution and write $|E| \approx |\delta E_{\perp,\lambda}| \sim \delta\Phi_{\lambda}/\lambda$. The other terms on the right-hand side of Equation (8) are then ordered as follows:

$$ \frac{|u_i \times B|}{c} \sim \delta u_{\perp,i,\lambda} \frac{B_0}{c}, \quad (9) $$

<footnote n="7">⁷ In contrast, perpendicular ion heating measured in low-resolution hybrid-kinetic simulations of decaying Alfvén-wave turbulence by Vasquez (2015) suggests that $c_2 \lesssim 0.03$, if $\xi_{th}$ is calculated using the $E \times B_0$ drift evaluated on scales in the vicinity of $\rho_{th,i}$.</footnote>

<footnote n="8">⁸ Klein & Chandran (2016) also allowed for a velocity-dependent exponential suppression in their formulation of $D_{\perp}^{E}(w_{\perp})$ (see their Equations (8) and (17), associating $c_*$ with $c_2 = 0.2$.</footnote>

<footnote n="9">⁹ Here we are considering scales relevant to stochastic ion heating, i.e., $k_{\perp} \rho_{th,i} \sim 1$. In our treatment, electron-inertia terms and electron finite-Larmor radius corrections can be neglected in Equation (8), if $k_{\perp} d_e \ll 1$ and $k_{\perp} \rho_{th,e} \ll 1$ hold at ion scales. This means that we are considering a range of $\beta_i$ that is still larger than the (small) electron-to-ion mass ratio, i.e., $m_e/m_i \ll \beta_i \lesssim 1$, as well as a range of temperature ratio, $\tau_{\perp}$, that is smaller than the (large) inverse of such a mass ratio, i.e., $0 \ll \tau_{\perp} \ll m_i/m_e$.</footnote>
$$ \frac{|\mathbf{J} \times \mathbf{B}|}{enc} \sim v_{\text{A0}} \frac{d_i}{\lambda} \left( \frac{\delta B_{\parallel, \lambda}}{B_0} + \frac{\lambda}{\ell_{\parallel, \lambda}} \frac{\delta B_{\perp, \lambda}}{B_0} \right) \frac{B_0}{c}, \quad (10) $$

$$ \frac{|\nabla p_e|}{en} \sim c_s \frac{\rho_s}{\lambda} \frac{\delta n_\lambda}{n} \frac{B_0}{c}, \quad (11) $$

where $v_{A0} \doteq B_0/\sqrt{4\pi m_i n}$ is the Alfvén speed, $c_s \doteq \sqrt{T_c/m_i}$ is the sound speed, and $\rho_s \doteq c_s/\Omega_i$ is the sound radius. In the Hall term (Equation (10)), $\ell_{\parallel, \lambda}$ is the characteristic length scale along the magnetic-field direction of a fluctuation with perpendicular extent $\lambda$; the ratio $\Theta_\lambda \doteq \lambda/\ell_{\parallel, \lambda}$ is related to the (possibly scale-dependent) anisotropy of the turbulent cascade.

Finally, we assume that the sub-ion-scale fluctuations are composed primarily of KAWs, an assumption supported by measurements in the solar wind (e.g., Chen 2016, and references therein). Such a cascade satisfies an approximate perpendicular pressure balance (Schekochihin et al. 2009; Kunz et al. 2018): $\delta n_\lambda/n \approx -(2/\beta_\perp)\delta B_{\parallel, \lambda}/B_0$, where $\beta_\perp \doteq (1+\tau_\perp)/\beta_{\perp,i}$. This allows one to combine the thermoelectric potential with the $\delta B_{\parallel, \lambda}$ term in the electrostatic piece of the Hall field to obtain

$$ \frac{\delta \Phi_\lambda}{\lambda} \sim \delta u_{\perp,i,\lambda} \frac{B_0}{c} + v_{A0} \frac{d_i}{\lambda} \left( \frac{1}{1+\tau_\perp} \frac{\delta B_{\parallel, \lambda}}{B_0} + \Theta_\lambda \frac{\delta B_{\perp, \lambda}}{B_0} \right) \frac{B_0}{c}. \quad (12) $$

For a critically balanced Alfvénic cascade with enough separation between the outer scale $L$ and $\rho_{\text{th,i}}$, the spectral anisotropy $\Theta_\lambda$ becomes $\ll 1$ as the ion-kinetic scales are approached. As a result, the contribution from the $\Theta_\lambda \delta B_{\perp, \lambda}$ term in Equation (12) at a given perpendicular scale $\lambda \ll L$ may be small enough when compared to that of the field-parallel fluctuations, $\delta B_{\parallel, \lambda}$, to be neglected. (Note that $\delta B_\perp/\delta B_\parallel \approx \sqrt{1+2/\beta_\perp}$ for KAW-like fluctuations, e.g., see Section 3.6.2 of Kunz et al. 2018.) We make this assumption in the remainder of the paper and drop the term $\propto \delta B_{\perp, \lambda}$ in Equation (12).<fnref n="10" />

Converting Equation (12) without the $\Theta_\lambda$ term into the velocity-space potential $\delta\Phi_w$ and inserting it into Equation 2(b) (i.e., neglecting the multiplicative exponential suppression factor in Equation (7) for the moment), one obtains an analytic formula for the perpendicular-energy diffusion coefficient,

$$ \frac{D_{\perp \parallel}^{E}}{\Omega_i m_i^2 v_{\text{th,i}}^4} \sim \left( \frac{w_{\perp}}{v_{\text{th,i}}} \right) \left( \frac{\delta u_{\perp,i,w}}{v_{\text{th,i}}} + \frac{1}{\beta_{\perp}} \frac{v_{\text{th,i}}}{w_{\perp}} \frac{\delta B_{\parallel,w}}{B_0} \right)^3. \quad (13) $$

Equation (13) implies that, depending on the spectral slopes of the fluctuation spectra at the ion gyroradii, ions with different perpendicular energies will diffuse differently in velocity space. This dependence is computed in Section 2.3, where we assign various spectral scaling laws to $\delta u_{\perp,i,\lambda}$ and $\delta B_{\parallel, \lambda}$ that correspond to different regimes of AW/KAW turbulence. These are then substituted into Equation (12) with $\lambda \sim (w_\perp/v_{\text{th,i}})\rho_{\text{th,i}}$, thereby yielding the velocity-scale dependence of $\delta\Phi_w$ and, through

Equations 2(b) and (4), $D_{\perp \parallel}^E$ and $Q_\perp$. In preparation for this exercise, we first advance arguments for which of the terms in Equation (12) provides the dominant contribution to the potential as seen by a particle with Larmor radius $\rho_i$ (when compared to the thermal gyroradius, $\rho_{\text{th,i}}$, and the ion skin depth, $d_i$, and thus depending upon $\beta_\perp$ as well).

Stochastic heating of an ion with perpendicular random velocity $w_\perp$ involves fluctuations that occur on scales comparable to that ion's gyroradius, $\lambda \sim \rho_i = w_\perp/\Omega_i$. This scale must be compared with the ion-kinetic scales of the background plasma, namely $\rho_{\text{th,i}}$ and $d_i$, which determine the nature of the turbulent fluctuations at scale $\lambda$ and thus the corresponding ordering of the different terms in Equation (12). These background spatial scales also have a corresponding scale in perpendicular velocity, namely the ion-thermal and Alfvén speeds, $v_{\text{th,i}} = \Omega_i \rho_{\text{th,i}}$ and $v_A = \Omega_i d_i = v_{\text{th,i}}/\sqrt{\beta_{\perp,i}}$, respectively. Just as the spatial scales determine the type of fluctuations that are responsible for the stochastic heating, these background velocity scales—and how they compare with the ion's velocity $w_\perp$—determine the corresponding ordering of the different terms in Equation (13). Moreover, as discussed in Section 2.1, for a quasi-Maxwellian distribution, we expect that the largest contribution to the total stochastic heating is provided by ions with $w_\perp \sim v_{\text{th,i}}$. The contribution from those ions whose perpendicular velocity exceeds a few times the ion-thermal speed, $w_\perp \gg v_{\text{th,i}}$, is exponentially suppressed. Similarly, the contribution from low-$w_\perp$ ions (i.e., those with $w_\perp \ll v_{\text{th,i}}$) to the overall heating would be progressively less important due to the strong dependence of $D_{\perp \parallel}^E \propto \delta\Phi_w^E$ on the fluctuations' amplitudes (viz., the lower the $w_\perp$, the smaller the spatial scale $\lambda \sim \rho_i \propto w_\perp$ at which the potential is sampled). Therefore, based on these arguments and what we know about the cascade of Alfvénic fluctuations, we may anticipate the following features of stochastic heating in the different $\beta_\perp$ regimes.

We first consider Equation (13) at $w_\perp \approx v_{\text{th,i}}$. When $\beta_{\perp,i} \gtrsim 1$, we have $\rho_{\text{th,i}} \gtrsim d_i$, so the ion-thermal gyroradius is encountered sooner by the cascading fluctuations than is the ion skin depth. At such a scale, the incompressive AW-like $\delta u_{\perp,i}$ fluctuations are still dominant over their compressive KAW-like $\delta B_\parallel$ counterparts (e.g., Cerri et al. 2017a, 2017b; which are also suppressed by an additional factor $\beta_{\perp,i}^{-1}$ in Equation (13) when $\beta_\perp > 1$). As a result, for $\beta_{\perp,i} \gtrsim 1$, we expect that the main contribution to the overall stochastic heating of ions is provided by the potential associated with the inductive term in Equation (8).

On the other hand, if $\beta_{\perp,i} \ll 1$, then the ion-thermal Larmor radius is much smaller than the ion skin depth, $\rho_{\text{th,i}} \ll d_i$, and turbulent fluctuations encounter $d_i$ as the first ion-kinetic scale in their cascade. Because the ions decouple from the dynamics of the magnetic field at sub-$d_i$ scales, the spectrum of ion flow velocity fluctuations becomes much steeper than its magnetic counterpart (an effect captured by the $J \times B/enc$ Hall term in Equation (8)). Accordingly, $\delta u_{\perp,i}$ fluctuations are negligibly small at $\lambda \ll d_i$ relative to magnetic-field fluctuations, a feature that has been seen in both in situ measurements of solar-wind turbulence (e.g., Šafránková et al. 2016; Chen & Boldyrev 2017) and kinetic numerical simulations of Alfvénic turbulence (e.g., Cerri et al. 2017a; Franci et al. 2018; Arzamasskiy et al. 2019). Moreover, at $\beta_\perp < 1$, the compressive KAW-like $\delta B_\parallel$ contribution to Equation (13) is now further enhanced by the factor $\beta_{\perp,i}^{-1}$. As a result, in the low-$\beta$ regime,

<footnote n="10">¹⁰ In our simulations (see Section 3), $\Theta \approx 0.05$ at $k_\perp \rho_{\text{th,i}} \approx 1$. This corresponds to an angle between the fluctuations' wavevector, $\mathbf{k}$, and the local background magnetic-field direction (i.e., using a scale-dependent definition of the background magnetic field, $B_{\text{loc}}(\mathbf{r}, t)$, computed via 5-point increments; Cerri et al. 2019) of $\vartheta_{(k_B)} = \arctan(\Theta^{-1}) \approx 87^\circ$. We note that $\vartheta_{(k_B)} \approx 80^\circ-90^\circ$ for fluctuations measured in the near-Earth solar wind with spacecraft-frame frequencies $f_{\text{spacecraft}} \sim 1$ Hz (Sahraoui et al. 2010).</footnote>
we anticipate the main contribution to the overall stochastic heating of ions to be provided by the potential associated with the nonideal terms in Equation (8).

## 2.3. Explicit Scalings for Stochastic Ion Heating from a Critically Balanced Alfvénic Cascade

In this section, we utilize well-known spectral scaling relations for $\delta u_{\perp,i,\lambda}$ and $\delta B_{\parallel,\lambda}$ in AW and KAW turbulence to evaluate Equation (13) and the associated perpendicular heating rate, Equation (4). To keep our expressions compact, we neglect for the time being the exponential suppression factor. A brief comment on how this factor modifies the results is then provided in Section 2.3.3; the full calculation with the factor included is reported in Appendix B. Strictly speaking, the contents of this section (Section 2.3) are not fully self-consistent, in that the transfer of turbulent energy to the thermal energy of the particles via stochastic heating is not accounted for in the adopted spectral scalings (which are power-law in form). However, it does allow us to gain some intuition for how $D_{\perp}$ might scale with $w_{\perp}$ and how the perpendicular heating rate per unit mass $Q_{\perp}$ depends on the plasma parameters. In doing so, we are most closely following Klein & Chandran (2016), who noted that their approach neglects the back reaction of the heating process on the turbulent power spectrum. A self-consistent determination of $D_{\perp}$ and $Q_{\perp}$ follows in Section 3, where we obtain spectral scalings for $\delta\Phi_{\lambda}$, $\delta u_{\perp,i,\lambda}$, and $\delta B_{\parallel,\lambda}$ from self-consistent numerical simulations and use them in Equations (5), (7), and (13) to determine $D_{\perp}$ and $\partial Q_{\perp}/\partial w_{\perp}$.

Consider an inertial-range cascade of large-scale (MHD) Alfvénic fluctuations characterized by a constant energy cascade rate per unit mass $\epsilon_{AW}$ and $\delta u_{\perp,i,\lambda} \sim (\epsilon_{AW}\lambda)^{1/3}$. This cascade is taken to exhibit a scale-dependent spectral anisotropy governed by critical balance (Goldreich & Sridhar 1995; Horbury et al. 2008), such that the characteristic field-parallel length scale of a fluctuation of perpendicular size $\lambda$ satisfies $\ell_{\parallel,\lambda} \sim L^{1/3}\lambda^{2/3}$, where $L \sim v_A^3/\epsilon_{AW}$ is the outer scale. As the ion-kinetic scales are approached, the AWs mutate into KAWs, with a fraction $\epsilon_{KAW}/\epsilon_{AW}$ of the inertial range cascade energy penetrating down into the dispersive range.<fnref n="11" /> For the sub-ion-scale KAW cascade, we do not adhere to any particular prescription for the associated wavevector anisotropy, using instead a generalized version of Equation (4.47) of Kunz et al. (2018),

$$ \ell_{\parallel, \lambda} \sim \left( \frac{\epsilon_{AW}}{\epsilon_{KAW}} \right)^{1/3} \frac{(1 + \tau_{\perp})^{1/6}}{(2 + \beta_{\perp})^{1/6}} L^{1/3} \rho_{\text{th,i}}^{2/3} \left( \frac{\lambda}{\rho_{\text{th,i}}} \right)^{\alpha/3}, \quad (14) $$

in which the anisotropy is parameterized by the exponent $\alpha$ (Cerri et al. 2018). Different values of $\alpha$ may result from assuming different nonlinear energy transfer timescales that govern the critically balanced cascade. For example, $\alpha = 1$ corresponds to a conservative KAW cascade with spectral

slope $-7/3$, as predicted by the gyrokinetic theory (e.g., Schekochihin et al. 2009). Accounting for a scale-dependent volume-filling factor of the KAW fluctuations instead yields $\alpha = 2$, with an associated KAW spectrum having a slope of $-8/3$ (Boldyrev & Perez 2012). Finally, $\alpha = 3$ corresponds to a scale-independent anisotropy, a feature sometimes seen in hybrid-kinetic simulations of AW/KAW turbulence (e.g., Franci et al. 2018; Arzamasskiy et al. 2019) and predicted by theories of reconnection-mediated Alfvénic turbulence (Loureiro & Boldyrev 2017; Mallet et al. 2017).

### 2.3.1. Stochastic Heating in $\beta \gtrsim 1$ AW Turbulence

When $\beta_{\perp} \gtrsim 1$, the nonlinear fluctuations approaching the ion-Larmor scale are composed primarily of AWs. Therefore, the main contribution to the electrostatic potential in Equation (12) is from the $\delta u_{\perp}$ fluctuations, and the diffusion coefficient can be approximated by

$$ \frac{D_{\perp}^{E}}{\Omega_{i} m_{i}^{2} v_{\text{th},i}} \approx \left( \frac{w_{\perp}}{v_{\text{th},i}} \right) \left( \frac{\delta u_{\perp, i, w}}{v_{\text{th},i}} \right)^{3}, \quad (15) $$

with the Alfvénic fluctuations satisfying

$$ \frac{\delta u_{\perp,i,\lambda}}{v_{\text{th},i}} \sim \left( \frac{\epsilon_{\text{AW}}}{\Omega_i v_{\text{A0}}^2} \right)^{1/3} \beta_{\perp,i}^{-1/3} \left( \frac{\lambda}{\rho_{\text{th},i}} \right)^{1/3}. \quad (16) $$

Substituting this expression into Equation (15) with $\lambda/\rho_{\text{th,i}} \sim w_{\perp}/v_{\text{th,i}}$ yields

$$ D_{\perp}^{(\text{AW})} \sim \epsilon_{\text{AW}} m_i^2 v_{\text{th,i}}^2 \left( \frac{w_{\perp}}{v_{\text{th,i}}} \right)^2. \quad (17) $$

By using Equation (4) and adopting for simplicity a Maxwellian distribution function in $w_{\perp}$, $f^E(w_{\perp}) = \exp(-w_{\perp}^2/v_{\text{th}}^2)/(m_i v_{\text{th,i}}^2)$, we find that the perpendicular heating rate per unit mass is given by

$$ \frac{Q_{\perp}^{(\text{AW})}}{\epsilon_{\text{AW}}} = \Lambda_{\text{AW}}, \quad (18) $$

where $\Lambda_{\text{AW}}$ is a constant independent of $\beta_{\perp}$ and $\tau_{\perp}$ that takes into account the various coefficients neglected in our scaling arguments. Therefore, at any $\beta_{\perp} \gtrsim 1$, the stochastic-heating rate (associated with AW-like fluctuations only) obtains an approximately constant fraction of the energy cascade rate. This result is consistent with the one in Chandran et al. (2010) for the case in which the dominant contribution to the electric-field fluctuations is due to the $\delta u_i \times B_0/c$ induction (and the exponential suppression factor is neglected; cf. their Equation (31)).

### 2.3.2. Stochastic Heating in low-$\beta$ KAW Turbulence

When $\beta_{\perp} \ll 1$, the ion-Larmor radius is smaller than the ion skin depth, $\rho_i \sim \rho_{\text{th,i}} \ll d_i$. As a result, the fluctuating potential (Equation (12)) evaluated at ion-Larmor scales is dominated by the contribution from the $\delta B_{\parallel,\lambda}$ fluctuations, and the diffusion coefficient can be approximated by

$$ \frac{D_{\perp}^{E}}{\Omega_{i} m_{i}^{2} v_{\text{th},i}} \approx \beta_{\perp}^{-3} \left( \frac{w_{\perp}}{v_{\text{th},i}} \right)^{-2} \left( \frac{\delta B_{\parallel, w}}{B_0} \right)^{3}, \quad (19) $$

<footnote n="11">¹¹ In gyrokinetic turbulence, the AW energy that does not make its way into the KAW cascade channel while going through the ion-kinetic scales is transferred into ion-thermal energy through Landau damping and/or a perpendicular phase-space cascade of ion-entropy fluctuations (Schekochihin et al. 2009). Here we allow for a portion of the cascading energy to also go into perpendicular stochastic heating of the ions. According to the discussion that follows Equation (13) in Section 2.2, when $\beta_{\perp} \gtrsim 1$, this heating mechanism drains a portion of the energy carried by the AW cascade ($\epsilon_{\text{AW}}$), while it is a portion of the KAW cascade ($\epsilon_{\text{KAW}}$) that is going into such an ion-energy channel at $\beta \ll 1$.</footnote>
with compressive KAW-like fluctuations satisfying

$$ \frac{\delta B_{\parallel, \lambda}}{B_0} \sim \left( \frac{\varepsilon_{\text{KAW}}}{\Omega_i v_{\text{A0}}^2} \right)^{1/3} \frac{\beta_{\perp i}^{1/3}}{(1 + 2/\beta_{\perp})^{1/3}} \left( \frac{\lambda}{\rho_{\text{th,i}}} \right)^{(3+\alpha)/6} . \quad (20) $$

Substituting this expression into Equation (19) with $\lambda/\rho_{\text{th,i}} \sim w_{\perp}/v_{\text{th,i}}$ yields

$$ D_{\perp}^{(\text{KAW})} \sim \varepsilon_{\text{KAW}} m_i^2 v_{\text{th,i}}^2 \frac{(1 + \tau_{\perp})^{-2}}{(2 + \beta_{\perp})} \left( \frac{w_{\perp}}{v_{\text{th,i}}} \right)^{(\alpha-1)/2} . \quad (21) $$

For $\alpha = 1$, $D_{\perp}^{(\text{KAW})}$ is independent of $w_{\perp}$; for $\alpha = 2$, $D_{\perp}^{(\text{KAW})} \propto w_{\perp}^{1/2}$; and for $\alpha = 3$, $D_{\perp}^{(\text{KAW})} \propto w_{\perp}$. Again adopting a Maxwellian distribution function in $w_{\perp}$, we may estimate the perpendicular heating rate per unit mass in low-$\beta$ KAW turbulence as

$$ \frac{Q_{\perp}^{(\text{KAW})}}{\varepsilon_{\text{KAW}}} = \Lambda_{\text{KAW}} (1 + \tau_{\perp})^{-2} (2 + \beta_{\perp})^{-1}, \quad (22) $$

where $\Lambda_{\text{KAW}}$ is a constant independent of $\beta_{\perp}$ and $\tau_{\perp}$.

If we further make the assumption that the transition from the AW cascade to the KAW cascade occurs at and is continuous across $k_{\perp}\rho_{\text{th,i}} \sim 1$, then we may estimate $\varepsilon_{\text{KAW}}/\varepsilon_{\text{AW}} \sim (\tau_{\text{AW}}/\tau_{\text{KAW}})k_{\perp}\rho_{\text{th,i}} \sim 1 \sim (2 + \beta_{\perp})^{-1/2}(1 + \tau_{\perp})^{1/2}$, in which case

$$ \frac{Q_{\perp}^{(\text{KAW})}}{\varepsilon_{\text{AW}}} \propto (1 + \tau_{\perp})^{-3/2} (2 + \beta_{\perp})^{-3/2}. \quad (23) $$

If instead the transition were to occur at $k_{\perp}d_i \sim 1$ (e.g., Chen et al. 2014), then

$$ \frac{Q_{\perp}^{(\text{KAW})}}{\varepsilon_{\text{AW}}} \propto \beta_{\perp}^{1/2} (1 + \tau_{\perp})^{-2} (2 + \beta_{\perp})^{-3/2}. \quad (24) $$

### 2.3.3. Exponential Attenuation

As forewarned at the start of Section 2.3, we have been omitting the exponential suppression factor introduced in Equation (7) to keep the limiting expressions for $D_{\perp}^E$ and $Q_{\perp}$ in different $\beta$ regimes compact. When this correction is included, the diffusion coefficient $D_{\perp}^E(w_{\perp})$ acquires a peak at a certain velocity, $w_{\perp}^{(\text{peak})}$, corresponding to the “most affected” (or “quasi-resonant”) ion population. For example, an exponentially corrected diffusion coefficient of the form $D_{\perp}^E(w_{\perp}) \propto (w_{\perp}/v_{\text{th,i}})^a \exp[-c(w_{\perp}/v_{\text{th,i}})^b]$, with constants $a$, $b$, and $c \ge 0$, displays a peak at $w_{\perp} = (a/bc)^{1/b} v_{\text{th,i}}$ (except in the case of standard KAW anisotropy, for which $\alpha = 1$, $a = 0$, and $D_{\perp}^E$ is just an exponentially decreasing function of $w_{\perp}$). If the exponential suppression were important, then the differential perpendicular-heating rate, $\partial Q_{\perp}/\partial w_{\perp}$, would also peak, at $w_{\perp}^{(\text{peak})} = ((a+1)/bc)^{1/b} v_{\text{th,i}}$. This would result in stochastic heating occurring most strongly on length scales $\chi^{(\text{peak})} \approx (w_{\perp}^{(\text{peak})}/v_{\text{th,i}})\rho_{\text{th,i}}$. On the other hand, if the fluctuations are in a regime in which the exponential correction is not important, then, to the lowest order, we recover the cases discussed in Sections 2.3.1 and 2.3.2, viz., a power-law diffusion coefficient of the form $D_{\perp}^E(w_{\perp}) \propto (w_{\perp}/v_{\text{th,i}})^a$ and a differential heating $\partial Q_{\perp}/\partial w_{\perp}$ peaking at $w_{\perp}^{(\text{peak})} = \sqrt{(a+1)/2} v_{\text{th,i}}$

because of the $\partial f^E/\partial w_{\perp}$ factor. We refer the reader to Appendix B for details.

## 3. Numerical Verification

We test the theory presented in Section 2 using hybrid-kinetic simulations with the particle-in-cell code PEGASUS++ (Kunz et al. 2014; Arzamasskiy et al. 2021, in preparation). Our hybrid model consists of fully kinetic ions coupled to a massless, charge-neutralizing, isothermal electron fluid via the generalized Ohm’s law (Equation (8); see Arzamasskiy et al. 2019, for the model equations). While hybrid kinetics excludes electron kinetic effects such as electron Landau damping (e.g., TenBarge & Howes 2013; Told et al. 2016; Grošelj et al. 2017), it retains certain ion-energization mechanisms (such as stochastic heating and ion-cyclotron resonances) that are not included in other models often used to study turbulent dissipation in collisionless plasmas (e.g., gyrokinetics; Howes et al. 2008; Told et al. 2015; Kawazura et al. 2019). We refer the interested reader to Told et al. (2016) and Camporeale & Burgess (2017) for a comparison of linear modes in hybrid kinetics, gyrokinetics, and full kinetics. Similarly, a comparative study of fluctuations’ properties in three-dimensional hybrid- and full-kinetic turbulence at sub-ion scales can be found in Cerri et al. (2019).

### 3.1. Simulation Setup

We consider an initially uniform plasma with ion density $n_0$, threaded by a uniform background magnetic field $B_0 = B_0e_z$ and placed within a three-dimensional, periodic computational domain of size $L_{\perp}^2 \times L_z$ with $L_x = L_y = L_z$. Turbulence is driven continuously in this plasma via a random, incompressible external force $F_{\text{ext}}$, which excites ion momentum fluctuations in the $x-y$ plane perpendicular to $B_0$. The forcing is time decorrelated over the interval $\tau_{\text{corr}}$ using an Ornstein-Uhlenbeck process (see Section 2 of Arzamasskiy et al. 2019). Only the largest-scale modes with $k_{\perp}^F = 2\pi/L_{\parallel}$ and $k_{\perp}^F = [1, 2] \times 2\pi/L_{\perp}$ are driven. A critical balance of the largest-scale fluctuations is assured by choosing a forcing amplitude such that the rms mean velocity fluctuation, $u_{\text{rms}}$, satisfies $u_{\text{rms}}/v_{\text{A0}} \approx L_{\perp}/L_{\parallel}$ in the quasi-steady turbulent state. Accordingly, $\tau_{\text{corr}} = L_{\perp}/2\pi u_{\text{rms}} \approx L_{\parallel}/2\pi v_{\text{A0}}$ is proportional to the Alfvén crossing time $\tau_A = L_{\parallel}/v_{\text{A0}}$. At the smallest scales, dissipation of turbulent energy is achieved by means of a fourth-order hyperresistivity on the magnetic field and low-pass filters on the first two moments of $f_i$ (viz., $n_i$ and $n_i u_i$).

In this paper, we combine results from two simulations of low-$\beta$ turbulence: a simulation with $\beta_{i0} = 0.3$ presented by Arzamasskiy et al. (2019) and a new simulation with $\beta_{i0} = 1/9$. This new simulation employs an elongated box with $L_{\parallel} = 6L_{\perp} = 48\pi d_{i0} = 144\pi\rho_{i0}$, discretized into $N_x = N_y = 288$ and $N_z = 1728$ cells, achieving an isotropic resolution $\Delta x \approx 0.087d_{i0}$ ($\approx 0.26\rho_{i0}$). The simulated wavenumber space is then $0.25 \le k_{\perp}d_{i0} \le 36$ and $0.04 \le k_{\parallel}d_{i0} \le 36$ (corresponding to $0.083 \le k_{\perp}\rho_{i0} \le 12$ and $0.014 \le k_{\parallel}\rho_{i0} \le 12$). In each cell, the initial ion distribution function is represented with 512 particles (giving $\approx 73$ billion particles in total).<fnref n="12" /> We run this simulation for $\approx 7.6 \tau_A$, with the quasi-steady state developing around

<footnote n="12">¹² The $\beta_{i0} = 0.3$ simulation of Arzamasskiy et al. (2019) utilized a $\delta f$ method to reduce the impact of particle noise on the fluctuations. This new $\beta_{i0} = 1/9$ simulation adopts a full-$f$ scheme in order to better handle potentially strong local density variations that arise in this low-beta regime.</footnote>
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Compensated energy spectra (top panels) and local spectral slopes (bottom panels) for the $\beta_{i0} = 1/9$ simulation. (a) Wavenumber spectra (compensated by $(k_{\perp}\rho_{i0})^{3/2}$) of the magnetic field **B** (blue solid), the electric field **E** (red solid), the MHD component of the electric field $E_{mhd} = -u_i \times \mathbf{B}/c$ (green dashed), and the kinetic component of the electric field $E_{kin} = (\mathbf{J} \times \mathbf{B})/c - T_e \nabla n$ (purple dashed). (b) Local spectral slopes vs. $k_{\perp}\rho_{i0}$. (c) Frequency spectra (compensated by $(\omega/\Omega_{i0})^2$) of **B** and **E** (blue and red, respectively). (d) Local spectral slopes vs. $\omega/\Omega_{i0}$. The light red (light blue) shaded regions highlight the wavenumber/frequency ranges where stochastic (cyclotron) heating is considered to be important.</figcaption></figure>

$\approx 4.3 \tau_A$. Our results are time-averaged over the remaining $\approx 3.3 \tau_A$ (corresponding to $\approx 500 \Omega_{i0}^{-1}$). For the $\beta_{i0} = 0.3$ run, we define the quasi-steady state as starting from $t/\tau_A \approx 4.4$ and continuing to the end of the simulation at $t/\tau_A \approx 19.9$ (corresponding to $\approx 4110 \Omega_{i0}^{-1}$).

## 3.2. Fluctuation Spectra for $\beta_{i0} = 1/9$

Figure 1 presents energy spectra and scale-dependent spectral indices ("local slopes") for the $\beta_{i0} = 1/9$ run versus (a) and (b) the wavenumber $k_{\perp}$ perpendicular to $\mathbf{B}_0$ and (c) and (d) the frequency $\omega$ measured in the plasma frame. These fluctuations exhibit significantly different spectra than in the corresponding $\beta \sim 1$ case (e.g., see Cerri et al. 2019, and references therein). First, the MHD-range spectra of electric and magnetic fluctuations both show a slope shallower than the usual anisotropic-MHD $-5/3$ scaling (e.g., Goldreich & Sridhar 1995) and closer to $-3/2$. (This may be due to the limited scale separation between the driving scales and the ion skin depth.) Second, while the spectral slope of the electric-field energy in the kinetic range is extremely close to $-2/3$, the corresponding magnetic-field spectrum steepens continuously beyond the $-8/3$ predicted to accompany the $-2/3$ electric spectrum.

We interpret this sub-ion-Larmor steepening as a signature of energy dissipation due to ion-heating mechanisms. This interpretation is supported by the frequency spectra in Figure 1(c), which exhibit slopes close to the $-2$ corresponding to a conservative energy cascade at frequencies $\omega/\Omega_{i0} \lesssim 0.2$ but steepen progressively through the sub-ion-Larmor range.

As we will show in Section 3.3, there are two ion-heating mechanisms operating simultaneously in this range, namely stochastic and cyclotron heating. The corresponding approximate wavenumber ranges in which one of these mechanisms is measured to be dominant over the other one are indicated in Figures 1(a) and 1(b) as light red (light blue) shaded regions for stochastic (cyclotron) heating. These ranges have been determined via direct measurement of the ions' perpendicular heating versus $k_{\perp}$, which shows a first peak around $k_{\perp}\rho_{i0} \sim 1$ that we associate with stochastic heating and a second peak around $k_{\perp}\rho_{i0} \sim 3$ that we associate with cyclotron heating (see Figure 3 and accompanying discussion in Section 3.3). Although there would likely be an overlap between the actual ranges over which these mechanisms operate at sub-ion scales, for the sake of clarity, the extent of these regions in Figures 1(a) and 1(b) is taken to be between $k_0/\sqrt{e}$ and $k_0/\sqrt{e}$ (k_0 being the peak wavenumber of each mechanism), a range previously used to estimate the total amount of stochastic ion heating (see, e.g., Xia et al. 2013; Martinović et al. 2020). The highlighted wavenumber ranges also have corresponding frequency ranges, highlighted in panels (c) and (d). These frequency ranges are obtained using an approximate AW/KAW dispersion relation for the stochastic-heating range<fnref n="13" /> and, for cyclotron heating associated with the n = 1 resonance, considering a resonance broadening of roughly $\Delta\omega/\omega_0 \sim$

<footnote n="13">¹³ Namely, $\omega^2 = k_{\parallel}^2 v_A^2 [1 + (1 + \tau_{\perp}) k_{\perp}^2 / (2 + \beta_{\perp})]$ (this formula smoothly interpolates between the AW and KAW limits; cf. Equations (4)–(5) in Howes et al. 2008). Different approximations for the KAW limit (see, e.g., Lysak & Lotko 1996) provide similar qualitative results, viz., that $\omega \simeq \Omega_i$ at $k_{\perp} \rho_i \approx 3$.</footnote>
$1/k_0\rho_{i0}$ (light blue region in panels (c) and (d)). We mention that there are also higher-$n$ resonances (shown as vertical dotted lines), likely contributing to the overall cyclotron heating.<fnref n="14" /> A detailed analysis of the fluctuations' spectral features, structure functions, and turbulence-related dynamics (e.g., magnetic reconnection) will be reported on elsewhere.

Before providing diagnostic evidence supporting this claim —that the ion- and sub-ion-Larmor-scale spectral steepening we observe is attributable to particle energization via stochastic and cyclotron heating—we note that such an association between changes in spectral slopes and energy dissipation is a relatively old idea in the solar-wind context (Coleman 1968), one that continues to be employed today (e.g., Woodham et al. 2018). Indeed, the steepness of the magnetic spectrum has been shown to correlate with both the energy cascade rate and power level in the inertial range (Smith et al. 2006; Bruno & Trenchi 2014) and the thermal proton temperature (Leamon et al. 1998). A more recent example may be found in Figure 5 of Chen et al. (2019), which shows a gradual steepening of the magnetic-field power spectrum in the Earth's magnetosheath throughout the sub-ion-Larmor range. While this kind of steepening has been attributed in some theoretical models to electron Landau damping (Sahraoui et al. 2009; Howes et al. 2011; TenBarge et al. 2013; Passot & Sulem 2015), the resemblance between our Figure 1(b) and Figure 5(b) of Chen et al. (2019) is notable given that our simulations do not include electron kinetics.

### 3.3. Ion Heating in Low-$\beta$ Turbulence

In Section 3.2, we attributed the steepening of the magnetic spectrum in the sub-ion-Larmor range to the energization of ion particles through stochastic and cyclotron heating. Here we provide evidence for this interpretation using data taken from both the $\beta_{i0} = 1/9$ and 0.3 simulations. In particular, we examine the (gyrotropized) ion distribution function $f(w_\|, w_\perp)$ alongside direct measures of $dQ_\perp/dw_\perp$ and $dQ_\perp/d\log k_\perp$ from these simulations, which in turn enable the evaluation of $D_{\perp}^E$ via Equation (5). These quantities are then compared to the theoretical predictions presented in Section 2. Namely, the actual $\delta\Phi_\lambda$ fluctuation spectrum obtained from 80 (50) snapshots of the $\beta_{i0} = 1/9$ (0.3) simulation during its quasi-steady state is employed in the expression for the diffusion coefficient (Equation (7)) and the associated differential heating (Equation (5)), including the exponential correction; these quantities are then time-averaged. At the same time, we employ an analogous procedure that considers only the $\delta u_{\perp,\lambda}$ or $\delta B_{\perp,\lambda}$ fluctuations' spectrum in the approximate expression for $D_{\perp}^E$ (Equation (13), including the exponential suppression term); this allows us to separate out the MHD and kinetic (non-MHD) contributions to the diffusion coefficient and the associated differential heating (Equation (5)).

#### 3.3.1. Ion-heating Diagnostics

To obtain the differential heating rate in the simulations, the following procedures have been implemented in the PEGASUS+ code (see also Arzamasskiy et al. 2019). At a given time, the differential rate of perpendicular heating in velocity space is computed as the sum of the instantaneous rate of work done by

the electric field on each particle $p$. Namely, we compute

$$ \tilde{Q}_{\perp} = \frac{\partial^2 Q_{\perp}}{\partial w_{\|} \partial w_{\perp}} \doteq \sum_{p} E_{\perp p} \cdot w_{\perp p} \quad (25) $$

and

$$ \tilde{Q}_{\|} = \frac{\partial^2 Q_{\|}}{\partial w_{\|} \partial w_{\perp}} \doteq \sum_{p} E_{\|p} w_{\|p}, \quad (26) $$

where $E_p \doteq E(x_p)$ is the electric field at the position $x_p$ of the particle $p$ with peculiar velocity $w_p \doteq v_p - u_p$, where $u_p \doteq u(x_p)$ is the mean-flow velocity at the particle's position. Here $\perp$ and $\|$ are defined with respect to the actual magnetic-field direction at location $x_p$: $w_p = w_{\|p}b_p + w_{\perp p}$ and $E_p = E_{\|p}b_p + E_{\perp p}$, with $b_p \doteq B(x_p)/|B(x_p)|$ being the local magnetic-field unit vector. Each of the above quantities are then binned in a two-dimensional $(w_\|, w_\perp)$ space so that they are a function of the gyrotropic (peculiar) velocity space: $\tilde{Q}_\perp(w_\|, w_\perp)$ and $\tilde{Q}_\|$($w_\|$), $w_\perp$). The total perpendicular or parallel heating rate is obtained as their integrals over the whole $(w_\|, w_\perp)$ space. (Thus, for instance, the one-dimensional $w_\|$ integral of $\partial^2 Q_\perp / \partial w_\| \partial w_\perp$ provides $dQ_\perp/dw_\perp$.) To obtain the differential rate of heating in wavenumber space, e.g., $dQ_\perp/d\log k_\perp$, the electric field is Fourier-transformed and then evaluated in different log-spaced $k_\perp \doteq (k_x^2 + k_y^2)^{1/2}$ bins, $E_\perp(k_{\perp,\text{bin}})$, which are then used to compute the associated rate of work on all of the simulation particles. (In this case, the rate of work is integrated over the whole $w$ space during runtime, so that the simulation output is a function of the $k_\perp$ bins only; an updated version of this diagnostic that outputs the heating rate in the whole three-dimensional $(w_\|, w_\perp, k_\perp)$ space is currently under development.) In the following analysis, all of the above quantities are time-averaged over the quasi-steady state (hereafter denoted by $\langle \cdot \rangle$).

#### 3.3.2. Free Parameters in Theoretical Predictions

When the theoretical predictions presented in Section 2 are computed from the actual fluctuation spectra obtained from the simulations, the theory has essentially three free parameters: (i) a normalization constant in Equation (7), (ii) an order-unity constant $\kappa_0$ that specifies the “resonance-like condition” $k_\perp w_\perp / \Omega_{i0} = \kappa_0$ that is used to transform the fluctuations' spectra from wavenumber to perpendicular-velocity space, viz. $\delta\Phi(w_\perp) \leftrightarrow \delta\Phi(k_\perp)|_{k_\perp = \kappa_0 \Omega_{i0} / w_\perp}$, and (iii) the constant $c_*$ in the exponential suppression factor. The constant in the first item is determined by normalizing the perpendicular-energy diffusion coefficient obtained from the $\delta\Phi_{\text{tot}}$ fluctuations' spectra (Equation (7)) to the $D_{\perp}^E$ directly obtained from the simulation at a single velocity point in the $w_\perp \le v_{\text{th},i0}$ range (the exact point used in the following being $w_\perp / v_{\text{th},i0} = 0.8$, but we verified that using any value in the range $0.5 \lesssim w_\perp / v_{\text{th},i0} \lesssim 1$ did not qualitatively change the results). This very same normalization constant is then used consistently for all of the theoretical curves, i.e., $dQ_\perp/dw_\perp$ and $dQ_\perp/d\log k_\perp$, as well as for the theoretical predictions obtained via the different contributions to the total potential (viz., $\delta\Phi_{\text{mhd}}$ and $\delta\Phi_{\text{kin}}$). Concerning the value of $\kappa_0$ and $c_*$, we show the plots when $(\kappa_0, c_*) = (1.1, 0.09)$ are adopted for the $\beta_{i0} = 0.3$ simulation, and $(\kappa_0, c_*) = (1.25, 0.05)$ are used in the $\beta_{i0} = 1/9$ case.

<footnote n="14">¹⁴ The $n > 1$ resonances are not formally associated with KAW-like fluctuations but rather with other types of fluctuations being relevant at low $\beta$ (see, e.g., Cerri et al. 2016, 2017b; Grošelj et al. 2017).</footnote>
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Left column: comparison between the stochastic-heating theory presented in Section 2 and $\beta_{i0} = 1/9$ simulation results vs. $w_{\perp}/v_{th,i0}$. (a) Perpendicular distribution function averaged over the quasi-steady turbulent state, $\langle f(w_{\perp}) \rangle$ (solid line; the dotted line shows the initial Maxwellian distribution for reference). (b) Averaged perpendicular-energy diffusion coefficient, $\langle D_{\perp}^E \rangle$, from simulation (solid black line) and theory (using Equation (13) with the exponential suppression factor) when the full potential ($\delta\Phi_{\text{tot}}$, solid orange line) or only its ideal ($\delta\Phi_{\text{mhd}}$; dashed-dotted green line) or nonideal ($\delta\Phi_{\text{kin}}$; dashed purple line) contribution is used. (c) Averaged differential perpendicular heating, $\langle dQ_{\perp}/dw_{\perp} \rangle$. Right column: same as left column but using results from the $\beta_{i0} = 0.3$ simulation.</figcaption></figure>

These values seem to best fit the simulations' results. The difference in the two values of $\kappa_0$ accounts somewhat for the different duration of the quasi-steady turbulent stage in the two simulations and thus of the consequent total absolute heating of the ions during the runs (i.e., how $\rho_{\text{th,i}}$ changes in the longer $\beta_{i0} = 0.3$ simulation). Nevertheless, we have verified that as long as it is in the range $0.9 \lesssim \kappa_0 \lesssim 1.4$, the results do not change qualitatively. For what concerns the difference in the two values of $c_*$, we interpret it as the result of a different level of intermittency within the two runs (being larger at lower $\beta$). We have verified that, when varying $\kappa_0$, $c_*$ can also be slightly adjusted without qualitatively changing the results; values in the range $0.045 \lesssim c_* \lesssim 0.055$ are allowed at $\beta_{i0} = 1/9$, while the same holds for a range of values $0.08 \lesssim c_* \lesssim 0.11$ at $\beta_{i0} = 0.3$ (this case being less well constrained due to the

higher errors associated with the $dQ_{\perp}/d\log k_{\perp}$ diagnostics around $k_{\perp}\rho_{\text{th,i}} \lesssim 1$; see Figure 3 and footnote 18).

### 3.3.3. Velocity-space Dependence of Ion Heating

We begin by examining how the ion perpendicular distribution function $\langle f(w_{\perp}) \rangle$, the perpendicular-energy diffusion coefficient $\langle D_{\perp}^E \rangle$, and the associated differential perpendicular heating $\langle dQ_{\perp}/dw_{\perp} \rangle$ behave in $w_{\perp}$ space. These quantities are traced by the solid black lines in Figure 2; the results from $\beta_{i0} = 1/9$ (0.3) are in the left (right) column. These are to be compared with the theoretical predictions derived in Section 2 for the diffusion and heating coefficients obtained using the spectra of the total electrostatic potential (solid orange line), the
MHD part of the potential (dashed–dotted green line), and the kinetic (i.e., non-MHD) part of the potential (dashed purple line).

In both simulations, we observe an evolution of the
perpendicular distribution function, $f(w_{\perp})$, from its initial
Maxwellian (dotted black lines) toward a broader shape with a
flat-top core (solid black lines). This evolution is the
consequence of the heating mechanisms operating in the
turbulence. In particular, we attribute the development of a
flattened core to stochastic heating, following Klein &
Chandran (2016). This interpretation is supported by the two
lower panels of this figure, in which both the diffusion
coefficient $D_{\perp}^{E}$ and the differential heating $dQ_{\perp}/dw_{\perp}$ are fit
reasonably well by the theoretical curve for $w_{\perp} \lesssim v_{\text{th,i0}}$, i.e.,
where the flat-top core develops.<fnref n="15" /> From these curves, it is also
evident how the relative importance of the contribution to the
total stochastic ion heating from different fluctuations changes
with the plasma beta; as $\beta_{i0}$ decreases, the nonideal contrib-
ution to the electrostatic potential responsible for the stochastic
heating of the ions, $\delta\Phi_{\lambda,\text{kin}} \propto (1 + \tau_{\perp})^{-1}\delta B_{\|,\lambda}$, becomes
progressively more important than its ideal counterpart,
$\delta\Phi_{\lambda,\text{mhd}} \propto \lambda \delta u_{\perp,\lambda}$ (cf. Equations (12) and (13) and the
accompanying discussion). This is highlighted by explicitly
plotting the theoretical perpendicular diffusion coefficient (and
the associated differential perpendicular heating) when only the
ideal ($\delta\Phi_{\text{mhd}}$; dashed–dotted green line) or nonideal ($\delta\Phi_{\text{kin}}$;
dashed purple line) contributions to the total electrostatic
potential ($\delta\Phi_{\text{tot}}$; solid orange line) are used<fnref n="16" /> However,
Figure 2 also shows that theoretical curves fit neither the
diffusion coefficient $D_{\perp}^{E}$ nor the differential heating $dQ_{\perp}/dw_{\perp}$
over the full range of $w_{\perp}$. This can be understood by
considering the fact that (i) stochastic heating is not the only
mechanism involved in the heating of ions in our simulation,
and (ii) the differential heating in Figure 2 is the result of an
integration over $w_{\|}$ of a more structured $\tilde{Q}_{\perp}(w_{\|}, w_{\perp})$. A
discussion of heating signatures within the two-dimensional
$(w_{\|}, w_{\perp})$ space is provided in Section 3.5.

3.3.4. Fourier-space Dependence of Ion Heating

Figure 3 displays the complementary diagnostic, the (aver-
aged) differential heating in wavenumber space $\langle dQ_{\perp}/dlogk_{\perp} \rangle$,
measured in the $\beta_{i0} = 1/9$ (top panel; black solid line) and the
$\beta_{i0} = 0.3$ (bottom panel; black solid line) runs. Overlaid are the
theoretical curves corresponding to Equation (5) using the total
fluctuating potential (orange solid line), the MHD part of the

<figure><img src="image_4.png" /></figure>

**Figure 3.** Top panel: differential perpendicular energization averaged over the quasi-steady turbulent state, $\langle dQ_{\perp}/dlogk_{\perp} \rangle$, vs. $k_{\perp}\rho_{i0}$ in the $\beta_{i0} = 1/9$ simulation. The numerical result (black solid line) is compared with the theoretical prediction using the spectrum of the total electrostatic potential fluctuations in Equation (12), $\delta\Phi_{\lambda,\text{tot}}$ (orange solid line), and the approximations considering only the $\delta\Phi_{\lambda,\text{mhd}}$ (green dotted-dashed line) or $\delta\Phi_{\lambda,\text{kin}}$ (purple dashed line) spectrum; the exponential correction is included, with $c_{*} = 0.05$. The plots are obtained using the relation $k_{\perp}w_{\perp}/\Omega_{i0} = \kappa_{0}$ with $\kappa_{0} = 1.25$ to best fit the simulation's results in both the velocity and wavenumber spaces. The light red (light blue) shaded regions show the $k_{\perp}$ range where stochastic (ion-cyclotron) heating is considered to be important. Bottom panel: same as the top panel but for the $\beta_{i0} = 0.3$ simulation. Here $c_{*} = 0.09$ and $\kappa_{0} = 1.1$ have been adopted.

potential (green dotted–dashed line), and the kinetic part of the
potential (purple dashed line).<sup>17</sup>

At $\beta_{i0} = 1/9$, the differential heating clearly exhibits two distinct peaks in the perpendicular-wavenumber space: one at $k_{\perp}\rho_{i0} \approx 1$ and a second at $k_{\perp}\rho_{i0} \approx 3$. We interpret the first peak as the result of stochastic ion heating, consistent with the theoretical curves obtained when the actual fluctuations' spectra are employed in the expressions derived in Section 2. The second peak at $k_{\perp}\rho_{i0} \approx 3$ is interpreted as being due to ion-cyclotron heating associated with the $n=1$ cyclotron resonance, consistent with the fact that the frequency of the fluctuations reaches $\omega \approx \Omega_{i0}$ at such a value of $k_{\perp}\rho_{i0}$ (see Figure 1 and accompanying discussion). An additional (minor) contribution to the total ion heating can be seen at $k_{\perp}\rho_{i0} \gtrsim 6$, likely associated with the $n > 1$ cyclotron resonances discussed in Section 3.2. These two mechanisms, stochastic and ion-

<footnote n="15">¹⁵ The differential perpendicular energization $dQ_\perp/dw_\perp$, as measured in our simulations, exhibits some (subdominant) cooling effects at $w_\perp/v_{th,i0} \gtrsim 2$. Because these cooling features are also present at very early times (including the initial time, $t=0$), they are likely due to errors associated with numerical noise and interpolation of the fields to the particle positions. We have modeled this cooling feature using the first few snapshots of a simulation and removed it from $\langle dQ_\perp/dw_\perp \rangle$ in the quasi-steady state. While we have verified that this cooling correction behaves sensibly when applied at late times (see Figure 9 in Appendix A), one should consider the simulation curves in Figure 2 to be most reliable for $w_\perp/v_{th,i0} \lesssim 2$.</footnote>

<footnote n="16">¹⁶ Note that, while $\delta\Phi_{\text{tot},\lambda}$ is obtained as the potential part of the actual $\delta E_{\perp,\lambda}$ fluctuations, the two components $\delta\Phi_{\text{MHD},\lambda}$ and $\delta\Phi_{\text{kin},\lambda}$ are obtained via the approximate formulas using $\delta u_{\perp,i,\lambda}$ and $\delta B_{\|,\lambda}$, respectively (i.e., where an approximate perpendicular pressure balance has been used to rewrite $\delta n_\lambda$ fluctuations in terms of $\delta B_{\|,\lambda}$ and neglecting the anisotropy correction $\Theta_\lambda \delta B_{\perp,\lambda}$; see Equation (12)). For this reason, the curves obtained via the approximate formulae do not exactly overlap with the one obtained using the actual $\delta\Phi_{\text{tot}}$, especially at small $w_\perp$ (corresponding to small-scale wavelengths $\lambda$) where different fields (namely, $\delta n$ and $\delta B_\|$) are affected differently by numerical filters in the code.</footnote>

<footnote n="17">¹⁷ To obtain the theoretical predictions plotted in Figure 3, the theoretical lines of $dQ_{\perp}/dw_{\perp}$ corresponding to Equation (5), which are plotted in Figure 2, have been interpolated into log $k_{\perp}$ space. This procedure also takes into account the logarithmic spacing of the volume in passing from $dw_{\perp}$ to $dlogk_{\perp}$, i.e., that $dQ_{\perp}/dlogk_{\perp} = (\kappa_0 \Omega_{i0}/k_{\perp})[dQ_{\perp}/dw_{\perp}]_{w_{\perp}=\kappa_0 \Omega_{i0}/k_{\perp}}$.</footnote>
cyclotron heating, contribute roughly equally to the overall perpendicular heating of the ions at $\beta_{i0} = 1/9$: $Q_{\perp}^{\text{stoch}}/Q_{\perp}^{\text{tot}} \approx Q_{\perp}^{\text{cyccl}}/Q_{\perp}^{\text{tot}} \approx 50\%$.

The overall perpendicular ion heating at $\beta_{i0} = 0.3$ (Figure 3, bottom) is dominated by scales at which we expect ion-cyclotron heating to be important; stochastic heating accounts for at most a quarter of the total heating: $Q_{\perp}^{\text{cyccl}}/Q_{\perp}^{\text{tot}} \gtrsim 75\%$ and $Q_{\perp}^{\text{stoch}}/Q_{\perp}^{\text{tot}} \lesssim 25\%$.<fnref n="18" />

An important trend that arises from the above analysis is that (i) stochastic ion heating should become progressively more important than ion-cyclotron heating as the plasma $\beta$ decreases, and (ii) this result is mainly due to contributions from the nonideal electric field (and associated potential, $\delta\Phi_{\text{kin}}$) arising from the Hall and thermoelectric effects in Equation (8). In fact, while the ideal contribution to the stochastic heating from $\delta\Phi_{\text{mhd}}$ is nearly constant when passing from $\beta_{i0} = 0.3$ to $1/9$, the heating associated with $\delta\Phi_{\text{kin}}$ nearly doubles in its contribution. This in turn lowers the amount of the fluctuations' energy that is available when the ion-cyclotron frequency is reached in the cascade, consequently diminishing the contribution of the ion-cyclotron mechanism to the overall ions' perpendicular heating.

### 3.4. Intermittency Contributions to Stochastic Heating

To explore the degree of intermittency of the potential fluctuations (and its effect on the stochastic heating) in the $\beta_{i0} = 1/9$ simulation, in Figure 4, we report the probability density function (PDF) of the normalized total potential fluctuations, $q_i \delta\Phi_{\text{tot}}/m_i v_{\text{th},i0}^2$ (orange solid line), and its ideal and nonideal parts, $q_i \delta\Phi_{\text{mhd}}/m_i v_{\text{th},i0}^2$ (green dotted-dashed line) and $q_i \delta\Phi_{\text{kin}}/m_i v_{\text{th},i0}^2$ (purple dashed line), respectively. Equivalent Gaussian distributions are also drawn as gray lines (with the same line style as the potential contribution to which they correspond). These PDFs are computed on two different ranges of scales: (i) $1/\sqrt{e} \le k_{\perp} \rho_{i0} \le \sqrt{e}$ (upper panel), corresponding to the range where stochastic heating is considered to be the dominant ion-heating mechanism, and (ii) $k_{\perp} \rho_{i0} \ge 1$ (lower panel), corresponding to the entire sub-ion-gyroradius ("kinetic") range.

From a statistical point of view, Figure 4 clearly shows that, while the width of the overall fluctuation-amplitude distribution decreases toward smaller scales, the degree of intermittency of these fluctuations simultaneously increases. Both aspects are relevant for the enhancement of stochastic ion heating. Let us consider the range of scales reported in the upper panel of Figure 4 (viz. $1/\sqrt{e} \le k_{\perp} \rho_{i0} \le \sqrt{e}$). In this range around $k_{\perp} \rho_{i0} \sim 1$, the quantity $q_i \delta\Phi_{\text{tot}}/m_i v_{\text{th},i0}^2$ corresponds to (a generalized version of) the stochasticity parameter that has been previously used to estimate the efficiency of stochastic heating (e.g., Xia et al. 2013; Vasquez 2015; Martinović et al. 2020). First, one notices that the distribution of fluctuations' amplitudes is relatively broad in this simulation, even for an equivalent Gaussian distribution; this implies that, even without taking into account intermittency, gyroscale fluctuations are not negligibly small. This is further quantified by computing both the rms stochastic-heating parameter, $\xi^{\text{(rms)}}$, and an effective

<figure><img src="image_5.png" /><figcaption>**Figure 4.** The PDF of (normalized) potential fluctuations, $q_i \delta\Phi/m_i v_{\text{th},i0}^2$, from our $\beta_{i0} = 1/9$ simulation in the range of scales in which stochastic heating is considered to be the dominant ion-heating mechanism, $1/\sqrt{e} \le k_{\perp} \rho_{i0} \le \sqrt{e}$ (upper panel), and at all sub-ion-gyroradius scales, $k_{\perp} \rho_{i0} \ge 1$ (lower panel). Statistics of both the total potential, $\delta\Phi_{\text{tot}}$ (orange solid), and its ideal and nonideal parts, $\delta\Phi_{\text{mhd}}$ (green dotted-dashed) and $\delta\Phi_{\text{kin}}$ (purple dashed), respectively, are reported. Equivalent Gaussian statistics are also drawn as gray lines (with corresponding line styles). Both the "effective" and rms value of the stochasticity parameter (computed using the actual PDF of the fluctuations) is reported in each plot as $\xi^{\text{(eff)}}$ and $\xi^{\text{(rms)}}$, respectively. In the range $1/\sqrt{e} \le k_{\perp} \rho_{i0} \le \sqrt{e}$ (upper panel), $q_i \delta\Phi/m_i v_{\text{th},i0}^2$ corresponds roughly to the thermal stochasticity parameter, $\xi_{\text{th}}$, which estimates the overall efficiency of stochastic heating (see Section 2.1).</figcaption></figure>

value, $\xi^{\text{(eff)}}$, that takes into account the non-Gaussian nature of the actual fluctuations' PDFs.<fnref n="19" /> These values are reported in each panel for the different scale ranges considered. Even in its rms version, within both scale ranges, the stochasticity parameter is large enough ($\xi \gtrsim 0.1$) that the overall effect of an exponential suppression term in Equation (7) should be small if $c_* \approx 0.01-0.1$. Second, intermittency does enhance the effective stochasticity parameter (and the associated heating). In fact, in the range of scales around $k_{\perp} \rho_{i0} \sim 1$ (upper panel of Figure 4), intermittency increases $\xi^{\text{(rms)}}$ by ≈19%. This effect is more important when the whole sub-ion range of scales is considered, $k_{\perp} \rho_{i0} \ge 1$ (lower panel of Figure 4); over this range of scales, $\xi^{\text{(eff)}}$ is increased beyond its equivalent-rms value $\xi^{\text{(rms)}}$ by ≈22% (although the absolute values of $\xi$ in this range are indeed smaller than the corresponding values in the range around $k_{\perp} \rho_{i0} \sim 1$).

The degree of intermittency also appears to depend on $\beta_i$. In the top panel of Figure 5, we report a comparison between the

<footnote n="18">¹⁸ The older $\beta_{i0} = 0.3$ simulation employed a heating diagnostic that used the total particle velocity $v_p$ in Equations (25) and (26) rather than its peculiar velocity $w_p$ (as in the version of the diagnostic employed in the new $\beta_{i0} = 1/9$ run). Also, the $k_{\perp}$ resolution used to compute this diagnostic was lower in the $\beta_{i0} = 0.3$ run (12 bins) than for $\beta_{i0} = 1/9$ (40 bins). As a result, the error bars on the heating at $k_{\perp} \rho_{i0} \lesssim 1$ are much larger in the $\beta_{i0} = 0.3$ run.</footnote>

<footnote n="19">¹⁹ Because the heating is proportional to $|q_i \delta\Phi/m_i v_{\text{th},i0}^2|^3 \approx |\xi|^3$, we define this effective parameter by $\xi^{\text{(eff)}} = [\int d\xi \, |\xi|^3 \mathcal{P}(\xi)]^{1/3}$, where $\mathcal{P}$ is the actual PDF of $q_i \delta\Phi/m_i v_{\text{th},i0}^2$.</footnote>
<figure><img src="image_6.png" /><figcaption>**Figure 5.** Top: comparison between the PDF of (normalized) total potential fluctuations, $\xi \doteq q_i \delta\Phi_{\text{tot}}/m_i v_{\text{th},i0}^2$, from the $\beta_{i0} = 1/9$ (orange solid line) and 0.3 (purple dashed line) simulations in the range of scales $1/\sqrt{e} \le k_{\perp} \rho_{i0} \le \sqrt{e}$. (The corresponding value of $\xi_{\text{tot}}^{(\text{eff})}$ is also reported below each simulation label). Gray lines (with corresponding line styles) represent the equivalent Gaussian distribution characterized by the same standard deviation $\sigma$ as the actual PDF. Bottom: comparison of the deviation from Gaussian statistics for the potential fluctuations (still in the range $1/\sqrt{e} \le k_{\perp} \rho_{i0} \le \sqrt{e}$) in the $\beta_{i0} = 1/9$ (orange solid) and 0.3 (purple dashed) simulations. This deviation is quantified both by the ratio of the actual PDF and the equivalent Gaussian vs. $\xi/\sigma$ (colored lines) and by the “excess kurtosis,” $K-3$ (reported in the plot below each simulation label; $K \doteq (\xi^4)/(\xi^2)^2 = 3$ for a zero-mean Gaussian distribution).</figcaption></figure>

PDFs of the normalized total potential fluctuations, $\xi \doteq q_i \delta\Phi_{\text{tot}}/m_i v_{\text{th},i0}^2$, around $k_{\perp} \rho_{\text{th},i} \sim 1$ in the $\beta_{i0} = 1/9$ (orange solid line) and 0.3 (purple dashed line) simulations. It is evident that the fluctuations' distribution broadens significantly at lower $\beta_{i0}$, passing from $\xi_{\text{tot}}^{(\text{eff})} \approx 0.1$ at $\beta_{i0} = 0.3$ to $\xi_{\text{tot}}^{(\text{eff})} \approx 0.17$ at $\beta_{i0} = 1/9$. This demonstrates that stochastic heating is enhanced as the plasma $\beta$ decreases, as expected. But we also find that the level of intermittency increases at lower $\beta$. In the bottom panel of Figure 5, we report the ratio between the actual PDF of $\xi$ and an equivalent-width Gaussian distribution characterized by the same standard deviation $\sigma$ of the actual PDF (because $\sigma$ depends on $\beta_{i0}$, the ratio is plotted versus $\xi/\sigma$ for the comparison to be meaningful). This PDF-to-Gaussian ratio exhibits larger deviations from unity at $\beta_{i0} = 1/9$ (orange solid line) than at $\beta_{i0} = 0.3$, a feature we further quantify by calculating the so-called “excess kurtosis,” $K-3$ (with the kurtosis defined by $K \doteq (\xi^4)/(\xi^2)^2$; $K=3$ for a Gaussian distribution with zero mean). This quantity doubles passing from $\beta_{i0} = 0.3$ (for which $K-3 \approx 0.28$) to $\beta_{i0} = 1/9$ (being

$K-3 \approx 0.56$). We interpret this enhanced intermittency as being responsible for decreasing the effective value of $c_*$ needed to fit our simulation results at different $\beta_{i0}$. Further numerical and observational studies are needed to determine the exact dependence of $c_*$ on the plasma parameters.

To further illustrate the partially intermittent nature of the stochastic ion heating in the $\beta_{i0} = 1/9$ run, we show the evolution of one of the simulation particles in Figure 6. This particle was specifically chosen because it increased its energy significantly over a short period of time by interacting with an intense spatially and temporally localized potential fluctuation.

Figure 6(a) shows the perpendicular ion energization $\langle q_i E_{\perp} \cdot w_{\perp} \rangle$ in the perpendicular (to the guide field) plane through which the tracked particle passed at that moment. The energization is averaged over multiple cells in the simulation (a volume of $18^3$ cells) to reduce the noise; it is normalized to $\delta u_{\text{rms}}^2/t_{\text{cross}}$, which serves as a proxy for the cascade rate. It is clear that the majority of the perpendicular energization happens in a spatially localized region. In Figure 6(b), we show the (normalized) potential fluctuations, $q_i \delta\Phi/m_i v_{\text{th},i0}^2$, in the same plane represented in panel (a); in this case, $\delta\Phi$ has been filtered to select only those modes satisfying $k_{\perp} \rho_{i0} > 1$. Comparing this contour plot with the one in panel (a), one can see a clear correlation between the region in which the amplitude of the potential fluctuations is larger and where most of the ion energization occurs. The majority of the energization happens in the region in which the Larmor-scale potential fluctuations are comparable to the thermal kinetic energy of a typical particle, $q_i \delta\Phi \sim m_i v_{\text{th},i0}^2$ (i.e., $\xi \sim 1$). As discussed earlier, the reason why such potential fluctuations can be so large, even though $q_i \delta\Phi_{\text{rms}}|_{k_{\perp} \rho_{i0} > 1} \sim 0.1 m_i v_{\text{th},i0}^2$, is because the turbulence is intermittent (cf. the bottom panel of Figure 4). This picture is supported by solar-wind measurements, which show a clear correlation between coherent magnetic structures generated intermittently by strong turbulence and plasma (anisotropic) heating (e.g., Osman et al. 2012; Greco et al. 2018; Qudsi et al. 2020).<fnref n="20" />

Panels (c) and (d) of this figure show this tracked particle's magnetic moment $\mu$ (green line), normalized to its initial value $\mu_{\text{th},i0} \doteq m_i v_{\text{th},i0}^2/2B_0$, and the particle's parallel and perpendicular thermal energies (blue and red lines, respectively), normalized to $v_{\text{th},i0}^2$, versus time. All of these quantities are approximately constant during particle gyration.<fnref n="21" /> However, once the particle enters the region with strong potential fluctuations (the gray shaded region in these panels), its perpendicular energy and magnetic moment oscillate with large amplitude. After $\approx 6$ gyrations, the particle's perpendicular energy and magnetic moment change by a factor of $\approx 4.2$.

Figures 4–6 further highlight the importance of intermittency in reducing the effectiveness of the exponential suppression factor introduced by Chandran et al. (2010), at least under the conditions realized in our simulations (see Section 2). This is

<footnote n="20">²⁰ From Chandran et al. (2010): “...in strong AW/KAW turbulence (as opposed to randomly phased waves), a significant fraction of the cascade power may be dissipated in coherent structures in which the fluctuating fields are larger than their rms values. Proton orbits in the vicinity of such structures are more stochastic than in average regions, and thus $c_2$ may be smaller in AW/KAW turbulence than in our test-particle simulations, indicating stronger heating.”</footnote>

<footnote n="21">²¹ These quantities ($\mu$, $w_\parallel$, $w_\perp$) are calculated using the magnetic field interpolated to the particle position, rather than to the particle's guiding center. The difference is responsible for the small periodic variations seen in these quantities on timescales of $\sim 2\pi/\Omega_{i0}$.</footnote>
<figure><img src="image_7.png" /><figcaption>**Figure 6.** Example of a simulation particle undergoing stochastic heating in the $\beta_{i0} = 1/9$ simulation. (Top) Snapshots of ion energization averaged over $18^3$ cells. (left) and small-scale ($k_{\perp} \rho_{i0} > 1$) potential fluctuations (right) in a plane perpendicular to the guide field. The black line shows a trajectory of the particle located in this plane. This particle starts in a region with small potential fluctuations and moves through a localized region with large $\delta\Phi$. (Bottom) Evolution of the particle's magnetic moment (left), along with parallel and perpendicular energies (right). The time over which the particle trajectory is plotted in the top panels is indicated by the gray shaded region. As the particle moves through strong potential fluctuations, it undergoes nonadiabatic perpendicular heating, which changes the particle's energy by a factor of a few over a timescale of several orbits.</figcaption></figure>

because of the relatively large rms amplitude of gyroscale potential fluctuations ($\xi_{th} = q_i \delta\Phi_{\rho_i} / m_i v_{th,i0}^2 \sim 0.1$) and the intermittent nature of those fluctuations, the latter of which causes a nonnegligible fraction of heating to occur in localized regions exhibiting large potential fluctuations. As a result, a particle's energy often changes considerably during just a few gyrations, and their orbits become stochastic, so that exponential conservation of the magnetic moment no longer holds.

As a final remark, we speculate that intermittency may allow stochastic heating to remain an important energization mechanism for low-$\beta$ turbulent systems even at scale separations much larger than what was achieved in our simulations (see, e.g., Mallet et al. 2019). As the scale separation increases, $\delta\Phi_{rms}|_{k_{\perp} \rho_{i0} \sim 1}$ decreases, but $\delta\Phi$ becomes localized within a smaller volume, creating larger potential drops within this volume. In other words, the trend outlined in Figure 4 for our $\beta_{i0} = 1/9$ run suggests that, while the PDF of the fluctuations' amplitude at ion/sub-ion scales may become progressively narrower as the scale separation $L/\rho_{i0}$ increases, the

intermittency effects will become simultaneously more and more important in enhancing $\xi^{(\text{eff})}$ with respect to $\xi^{(\text{rms})}$.

### 3.5. Other Signatures of Wave–Particle Interaction

The parallel and perpendicular ion-energization rates in the two-dimensional velocity space, $\tilde{Q}_{\|}(w_{\|}, w_{\perp})$ and $\tilde{Q}_{\perp}(w_{\|}, w_{\perp})$, respectively (see Equations (25) and (26)), can also be used to uncover the phase-space signatures of different wave–particle interactions. Their time-averaged values in the quasi-steady state, $\langle \tilde{Q}_{\|}(w_{\|}, w_{\perp}) \rangle$ and $\langle \tilde{Q}_{\perp}(w_{\|}, w_{\perp}) \rangle$, are reported in Figure 7 for both the $\beta_{i0} = 1/9$ (left column) and $\beta_{i0} = 0.3$ (right column) simulations. Figure 8 additionally provides this information for $\beta_{i0} = 0.3$ during its “early phase,” which refers to times $t/\tau_A \approx 3.8–4.4$ before the core of the perpendicular distribution function becomes appreciably flattened and stochastic heating is consequently reduced (see Figure 8 of Arzamasskiy et al. 2019).

The velocity-space patterns of $\langle \tilde{Q}_{\|}(w_{\|}, w_{\perp}) \rangle$ seen in the quasi-steady state of both simulations (Figures 7(c) and (d))
<figure><img src="image_8.png" /><figcaption>**Figure 7.** Ion-energization rate as a function of parallel ($w_\|\|$) and perpendicular ($w_\perp$) velocities. Panels (a) and (c) show parallel and perpendicular energization in the $\beta_{i0} = 1/9$ simulation. Panels (b) and (d) show the same quantities for the $\beta_{i0} = 0.3$ simulation from Arzamasskiy et al. (2019). The ion-thermal speed is marked by dashed lines; the Alfvén speed in the $\beta_{i0} = 0.3$ run is marked by dotted lines.</figcaption></figure>

display the signature of collisionless damping at the Landau resonances, $w_\| \approx \pm v_{\text{th,i0}}$ (cf. Howes et al. 2017). We interpret this structure as being due to the collisionless damping of slow-mode fluctuations. In the $\beta_{i0} = 1/9$ simulation, the amount of parallel energization associated with this Landau-resonant damping is extremely subdominant, contributing only $\lesssim 2\%$ of the total ion-heating rate. In the $\beta_{i0} = 0.3$ simulation, this percentage is $\lesssim 10\%$. During the early phase of the $\beta_{i0} = 0.3$ run (Figure 8(b)), there is an additional signature of wave-particle interaction in the vicinity of $w_\| \approx \pm v_{\text{A0}}$. We instead attribute the majority of the measured increase in parallel temperature to a combination of transit-time damping, which is driven by $Q_\perp$ (note the vertical resonant-like red and blue “stripes” in Figures 7(a) and (b) and 8(a)), and pitch-angle scattering of perpendicularly energized particles (as in Arzamasskiy et al. 2019, Section 3.2).<fnref n="22" />

In contrast with the $\beta_{i0} = 0.3$ case, the parallel ion distribution, $f(w_\|)$, does not develop significant nonthermal tails at $\beta_{i0} = 1/9$ (not shown). This can be explained by the inefficient Landau damping of Alfvénic fluctuations at very low values of $\beta$; at $\beta \ll 1$, the Alfvén speed is much larger than the ion-thermal velocity, $v_A \gg v_{\text{th,i}}$, and only the very tail of the ion distribution can effectively resonate with the phase velocity $v_{\text{ph}} \sim v_A$ of Alfvénic fluctuations. Since this population is energetically unimportant for the overall thermal budget of the plasma, we do not expect to find significant (parallel) heating from this process at very low $\beta$.<fnref n="23" />

Finally, both runs display signatures that may be interpreted as the superposition of (i) stochastic heating and (ii) ion-cyclotron heating. Stochastic heating presents in both runs as a horizontal feature close to $w_\| \sim v_{\text{th,i0}}$. For $\beta_{i0} = 0.3$, this signature is much more pronounced during its “early phase” (Figure 8(a)) than in its quasi-steady state, in which the core of the perpendicular distribution function is substantially flattened and stochastic heating is reduced. Ion-cyclotron heating, on the other hand, presents as a (fuzzy) circular halo centered around $w_\| \sim 0$ and $w_\perp \gtrsim v_{\text{th,i0}}$ (cf. Klein et al. 2020). However, the

<footnote n="22">²² Isenberg et al. (2019) suggested that the perpendicularly heated ion distribution functions with $T_\perp > T_\|\|$ that are naturally generated by ion stochastic heating would be unstable to the ion-cyclotron anisotropy instability, which would then generate quasi-parallel-propagating ion-cyclotron waves and thereby scatter ions in the parallel direction. The connection between this suggestion and the pitch-angle scattering of perpendicularly energized particles measured by Arzamasskiy et al. (2019) and also seen here is not clear for two main reasons. First, the temperature anisotropies measured in our simulations never become as large as those found in the model devised by Isenberg et al. (2019); for example, $T_\perp/T_\|\| \lesssim 1.12$ in our $\beta_{i0} = 1/9$ simulation. Second, our steady-state perpendicular distribution functions retain flattened cores similar to those predicted by Klein & Chandran (2016); Isenberg et al. (2019) predicted that pitch-angle scattering from unstable ion-cyclotron waves would erase this distinctive feature.</footnote>

<footnote n="23">²³ The same argument can also explain why, within gyrokinetic theory and simulations, the ion-to-electron heating dramatically drops at low $\beta$ (e.g., Howes 2010; Kawazura et al. 2019): because species’ heating in gyrokinetics relies only on the Landau damping of the fluctuations (which can thus provide only parallel heating), Alfvénic turbulence will be damped inefficiently by ions as the plasma $\beta$ decreases. (The large-scale injection of compressive fluctuations, which may be collisionlessly damped even at low $\beta$, at energy levels comparable to those of the Alfvénic fluctuations modifies this expectation; Kawazura et al. 2020.)</footnote>
<figure><img src="image_9.png" /><figcaption>Figure 8. Parallel and perpendicular energization at early times in the βᵢ₀ = 0.3 simulation from Arzamasskiy et al. (2019), before flattening of the perpendicular-velocity core of the distribution function suppresses stochastic heating. The ion-thermal (Alfvén) speed is marked by the dashed (dotted) lines.</figcaption></figure>

position and extension of this halo in $w_⊥$ seem to vary between βᵢ₀ = 0.3 and 1/9; this feature is not well understood and should be investigated in future work.

## 4. A Comment on the Interpretation of Stochastic Heating in Spacecraft Data

Before summarizing our main findings, we pause here to offer a comment on how spacecraft data might be best interpreted when looking for evidence of stochastic ion heating in the low-β solar wind. We begin by summarizing the method adopted by Bourouaine & Chandran (2013), Vech et al. (2017), and Martinović et al. (2019, 2020). Those authors used spacecraft-measured amplitudes of magnetic-field fluctuations near the proton gyroscale, $δB_{ρi}$, as a proxy for the gyroscale velocity fluctuations, $δu_{i,ρi}$. The latter was then divided by the field-perpendicular proton thermal speed, $v_{th,i} = \sqrt{2T_{⊥,i}/m_i}$, to obtain estimates for the stochasticity parameter $ϵ_i$ originally introduced by Chandran et al. (2010; recall footnote 6). Specifically, they set

$$δu_{i,ρi} = \sigma \frac{δB_{ρi}}{\sqrt{4πm_i n}}, \qquad (27)$$

where $σ$ is an order-unity constant (typically 1.19), so that

$$ϵ_i = \sigma β_{⊥,i}^{-1/2} \frac{δB_{ρi}}{B_0}, \qquad (28)$$

where $B_0$ is the mean magnetic-field strength. The amplitudes of the gyroscale magnetic-field fluctuations were defined using

$$δB_{ρi} \doteq \left[ \int_{\rho_i^{-1}/\sqrt{\epsilon}}^{\sqrt{\epsilon}\rho_i^{-1}} dk_{\perp} E_B(k_{\perp}) \right]^{1/2}, \qquad (29)$$

where $E_B(k_⊥)$ is the (appropriately normalized) one-dimensional magnetic energy spectrum in the plasma rest frame (obtained by applying Taylor's hypothesis to the frequency spectrum measured by the spacecraft). The amount of stochastic heating associated with these fluctuations was then inferred using

$$Q_{\perp} = \frac{v_{th,i}^3}{\rho_i} \left[ c_1 \epsilon_i^3 \exp\left(-\frac{c_2}{\epsilon_i}\right) \right], \qquad (30)$$

with $c_1 \sim 1$ (typically 0.75) and $c_2 \approx 0.1-0.3$ (typically 0.34 or $\simeq 0.2$). (Recall that the value of $c_*$ that best fits our simulation results is $\approx 0.05-0.1$.) Average values of $ϵ_i$ inferred between $\sim 0.3$ and $\sim 1$ au from the Sun were in the range of $\approx 0.03-0.05$.

The results of our paper suggest that the following refinements to this procedure may improve its accuracy. First, it is not necessarily the case that the fluctuations on ion gyroscales are accurately described by the Alfvénic relation (Equation (27)). Indeed, the argument in Section 2.2 is that the gyroscale potential fluctuations may be better inferred at low beta using $q_i \delta\Phi_{ρi}/m_i v_{th,i}^2 \sim β_{\perp}^{-1}(δB_{ρi}/B_0)$, rather than $\sim δu_{i,ρi}/v_{th,⊥,i} = σβ_{⊥,i}^{-1/2}(δB_{ρi}/B_0)$. (Recall the definition $β_{\perp} = (1 + τ_{\perp})(β_{\perp,i})$. While it is true that there are combinations of $τ_{\perp}$ and $β_{\perp,i}$ for which these two formulae return comparable inferred potential fluctuations, the interpretative difference is notable; at very low values of $β_{\perp}$, the electrostatic potential with which particles interact on their gyroscale has little to do with fluctuations in the ion flow velocity. When in doubt, a generalized Ohm's law that accounts for sub-$d_i$ contributions to the electrostatic potential, such as Equation (8), should be used.

To give concrete numbers, the rms fluctuation levels centered about the ion-thermal Larmor scale in our $β_{i0} = 1/9$ simulation (calculated as in Equation (29)) are $δB_{ρi}/B_0 \simeq 0.042$ and $δu_{i,ρi}/v_{A0} \simeq 0.024$; in our $β_{i0} = 0.3$ simulation, they are $δB_{ρi}/B_0 \simeq 0.043$ and $δu_{i,ρi}/v_{A0} \simeq 0.021$. Neither of these sets of values satisfy Equation (27) when $σ = 1.19$, and both suggest $σ < 1$. In this context, it is worth noting that these ion-Larmor-scale magnetic-field fluctuation amplitudes are typical of (if just slightly larger than) those in the low-beta solar wind; Bourouaine & Chandran (2013) used Helios data to report $δB_{ρi}/B_0 \approx 0.03$ at $\approx 0.3$ au, while Martinović et al. (2020) used Parker Solar Probe data to find strong evidence in the ion distribution function for stochastic heating at $\approx 0.2$ au when $δB_{ρi}/B_0 \simeq 0.049$ (see their Figure 5(a)). Both authors used Equation (28) with $σ = 1.19$ to compute $ϵ_i$, reporting values in the range $\approx 0.04-0.08$ when $β_i \approx 0.3-0.5$. The stochasticity parameter in our $β_{i0} = 1/9$ run, based on rms potential fluctuations centered about $ρ_i$, is notably larger at $ξ^{(rms)} \simeq 0.146$; accounting for intermittency raises its value to $ξ^{(eff)} \simeq 0.173$. In our $β_{i0} = 0.3$ run, we measured $ξ^{(rms)} \simeq 0.085$ and $ξ^{(eff)} \simeq 0.10$. Whether the difference between the observationally inferred $ϵ_i$ and the values of $ξ$ we obtained from the potential fluctuations in our simulations is primarily because Equation (27) is an inaccurate proxy for electrostatic potential fluctuations at low $β$, because intermittency effects must be
taken into account, or perhaps because our simulations could benefit from slightly larger scale separation awaits more data (both actual and numerical) and further scrutiny. Given the exponential sensitivity of $Q_{\perp}$ in Equation (30) to $\epsilon_i$, obtaining an accurate value of $c_2$ relies on an accurate definition of the stochasticity parameter.

## 5. Conclusions

We have derived a generalization of the theory of stochastic ion heating originally presented in Chandran et al. (2010), adapted to the case in which electric-field fluctuations can be described by a generalized Ohm's law that includes Hall and thermoelectric effects. We argued that these nonideal terms provide the dominant contribution to the stochastic heating of ions at sub-$d_i$ scales, which are the relevant scales at which stochastic heating operates in low-$\beta$ turbulence (i.e., when $\rho_i \ll d_i$). By keeping a fully scale-dependent approach, both in configuration space and in velocity space, we have derived the perpendicular-heating rate $Q_{\perp}$ and perpendicular-energy diffusion coefficient $D_{\perp}^E$ as functions of the perpendicular ion velocity $w_{\perp}$ and the perpendicular plasma beta $\beta_{\perp}$, adopting certain well-established properties of inertial- and dispersion-range turbulent fluctuations.

The predictions of this theory were then tested using three-dimensional hybrid-kinetic particle-in-cell simulations of continuously driven Alfvénic turbulence at low $\beta$, namely, the $\beta_{i0} = 0.3$ simulation presented by Arzamasskiy et al. (2019) and a newly performed $\beta_{i0} = 1/9$ simulation. In these simulations, parallel heating of ions is primarily associated with Landau/Barnes damping of turbulent fluctuations and is always subdominant with respect to its perpendicular counterpart, $Q_{\parallel,i} \ll Q_{\perp,i}$. Two perpendicular-heating mechanisms are shown to operate simultaneously on ions and provide most of their heating: ion-cyclotron and stochastic heating. While ion-cyclotron dominates over stochastic heating at $\beta_{i0} = 0.3$ ($Q_{\perp,i}^{\text{cycl}}/Q_{\perp,i}^{\text{tot}} \gtrsim 75\%$ and $Q_{\perp,i}^{\text{stoch}}/Q_{\perp,i}^{\text{tot}} \lesssim 25\%$), in the $\beta_{i0} = 1/9$ simulation, these two mechanisms contribute roughly equally to the perpendicular heating of ions ($Q_{\perp,i}^{\text{stoch}}/Q_{\perp,i}^{\text{tot}} \approx Q_{\perp,i}^{\text{cycl}}/Q_{\perp,i}^{\text{tot}} \approx 50\%$). As far as stochastic ion heating is concerned, the theoretical predictions derived in this work describe reasonably well the associated features emerging from the simulations and characterized by various heating diagnostics, both in perpendicular-velocity and in perpendicular-wavevector spaces. These diagnostics also emphasize the important role of non-MHD contributions to the electrostatic potential in stochastically heating the ions at low $\beta$ and demonstrate that intermittency in the turbulence enhances this heating. Finally, the fraction of injected energy that is channeled into total ion heating strongly depends on the plasma $\beta$, passing from $Q_i^{\text{tot}}/\epsilon_{AW} \approx 75\%$ at $\beta_i = 0.3$ to $Q_i^{\text{tot}}/\epsilon_{AW} \approx 40\%$ at $\beta_i \approx 0.1$.

Our work has three main implications for the interpretation of spacecraft data in the context of stochastic heating. First, we have provided a number of phase-space diagnostics that one may use to supplement the presently employed technique of inferring stochastic heating in the solar wind via correlations between the amplitudes of ion-Larmor-scale magnetic fluctuations and plasma heating. These diagnostics supplement concurrent work on field-particle correlations by Klein & Howes (2016), Howes et al. (2017), and others, which show great promise in their ability to distinguish between various particle-energization mechanisms and their contributions to the heating of the solar wind. Second, the precise way in which

spacecraft-measured, ion-Larmor-scale magnetic-field fluctuations are translated into electric potential fluctuations to calculate stochastic heating deserves careful reexamination, especially at $\beta$ values small enough that $d_i \gg \rho_{th,i}$. In particular, we advocate for the use of a generalized Ohm's law that accounts for the (sometimes dominant) contributions from the Hall and thermoelectric effects to the electric potential. We find that the implied stochasticity parameter $\xi_{th} = q_i \delta\Phi_{\rho_i}/m_i v_{th,i}^2$ obtained from the full potential fluctuations is generally larger than that implied by Equation (28), particularly when intermittency effects are taken into account. Third, our simulation results suggest a link between preferential perpendicular heating, magnetic spectra that exhibit sub-ion-Larmor steepening, and perpendicular distribution functions with flattened cores—a link that, if due to stochastic heating, should be pronounced when the amplitude of ion-Larmor-scale magnetic fluctuations is relatively large (viz., $\xi \gtrsim 0.1$).

With the gradual decrease in the perihelion of the Parker Solar Probe (Fox et al. 2016) and the increasing level of turbulent activity toward the Alfvén point (Tu & Marsch 1995; Chandran et al. 2011; Bruno & Carbone 2013; Chen et al. 2020), the importance of understanding the phase-space signatures of stochastic heating will only become greater. It is our hope that the predictions and diagnostics presented here will help to sharpen this understanding and facilitate a more robust analysis of current and future spacecraft data.

It is a pleasure to thank Kristopher Klein, Mihailo Martinović, and Benjamin Chandran for extremely useful discussions on stochastic ion heating and the interpretation of spacecraft measurements and Phil Isenberg, Bernie Vasquez, and the referee for constructive comments on the manuscript. We further acknowledge PRACE for awarding us access to the supercomputer Marconi, CINECA, Italy, where our $\beta_{i0} = 1/9$ simulation was performed, under grant No. 2017174107. This research was supported by NASA grant No. NNX16AK09G issued through the Heliophysics Supporting Research Program, the Max-Planck/Princeton Center for Plasma Physics (NSF grant PHY-1804048), and an Alfred P. Sloan Research Fellowship in Physics to M.W.K. L.A. gratefully acknowledges support from the Institute for Advanced Study.

# Appendix A

## A. Alternative Heating Diagnostics

In this Appendix, we summarize the implementation of the heating diagnostic in PEGASUS++ and discuss its limitations. We begin by reviewing our definition of the perpendicular-energy diffusion coefficient. If the particle heating occurs through a diffusion-like process, the distribution function evolves according to

$$ \frac{\partial f^E}{\partial t} = \frac{\partial}{\partial e_{\perp}} \left( D_{\perp}^E \frac{\partial f^E}{\partial e_{\perp}} \right), \quad (A1) $$

where $f^E$ is the perpendicular-energy distribution function. Here we assume that energization is only perpendicular to the magnetic field, $e_{\perp} \doteq w_{\perp}^2/2$. The total energy of the distribution
<figure><img src="image_10.png" /></figure>

<figure><img src="image_11.png" /><figcaption>**Figure 9.** Comparison between two methods for computing the perpendicular-energy diffusion coefficient. The blue line is obtained from the evolution of the energy distribution function following the method of Vasquez et al. (2020). The black line is computed using the $E \cdot w$ diagnostic used throughout this paper and in Arzamasskiy et al. (2019). The orange line represents the theoretical prediction for stochastic heating based on the electrostatic potential fluctuations (Equation (2b)). (Left) Energy diffusion coefficient. (Right) Velocity-space dependence of ion energization. The shaded regions represent the time variability of plotted quantities (computed as a standard deviation).</figcaption></figure>

is $E \equiv \int de_{\perp} e_{\perp} f^E$. The total energization rate is thus

$$Q_{\perp} \doteq \frac{dE}{dt} = \frac{d}{dt} \int de_{\perp} e_{\perp} f^E = \int de_{\perp} e_{\perp} \frac{\partial f^E}{\partial t} \\ = \int de_{\perp} e_{\perp} \frac{\partial}{\partial e_{\perp}} \left( D_{\perp}^{E} \frac{\partial f^E}{\partial e_{\perp}} \right). \qquad (A2)$$

In a collisionless plasma, this energization is provided by the work performed by the electric fields on the particle distribution,

$$q_{\perp}(e_{\perp}) \doteq q_i v_{\perp} \cdot E_{\perp}. \qquad (A3)$$

The perpendicular-energy distribution function of ions evolves according to a Vlasov equation of the form

$$\frac{\partial f^E}{\partial t} + \frac{\partial}{\partial e_{\perp}}(q_{\perp}(e_{\perp})f^E) = 0; \qquad (A4)$$

this equation neglects terms related to advection and exchange between parallel and perpendicular energies (e.g., pitch-angle scattering) but retains some basic properties of the kinetic equation, such as the conservation of the particle number. The equation for the diffusion coefficient is then simply

$$\frac{\partial}{\partial e_{\perp}} \left( D_{\perp}^{E} \frac{\partial f^E}{\partial e_{\perp}} \right) + \frac{\partial}{\partial e_{\perp}} (q_{\perp}(e_{\perp})f^E) = 0 \\ \Rightarrow \qquad D_{\perp}^{E} = - \frac{q_{\perp}(e_{\perp})f^E}{\partial f^E / \partial e_{\perp}}, \qquad (A5)$$

which could also be written as

$$D_{\perp}^{E} = - \frac{\partial Q_{\perp}}{\partial f^{E}} \frac{\partial}{\partial e_{\perp}}, \qquad (A6)$$

where we define the differential heating rate as

$$\frac{\partial Q}{\partial v} \doteq q_i v \cdot E f. \qquad (A7)$$

We use the definition (Equation (A6)) of the energy diffusion coefficient throughout this paper (as well as in Arzamasskiy et al. 2019) to describe the velocity-space dependence of ion heating. At large scales, energization (Equation (A7)) is

dominated by the conversion between bulk-kinetic and magnetic energies related to Alfvénic motions. In order to remove this nondissipative process, we use $w = v - u$ instead of $v$ for the heating diagnostics and perform a long-time average.

Recently, Vasquez et al. (2020) argued that the energy diffusion coefficient should be defined differently. They argued that differential energization should be equal to

$$\frac{\partial Q_{\perp}}{\partial e_{\perp}} = \frac{\partial e_{\perp} f^E}{\partial t}. \qquad (A8)$$

Using this definition, they arrived at a more complicated equation for the diffusion coefficient,

$$\frac{\partial Q_{\perp}}{\partial e_{\perp}} = -D_{\perp}^{E} \frac{\partial f^{E}}{\partial e_{\perp}} + \frac{\partial}{\partial e_{\perp}} \left( e_{\perp} D_{\perp}^{E} \frac{\partial f^{E}}{\partial e_{\perp}} \right), \qquad (A9)$$

which has an additional term relative to Equation (A6). The difference comes from the definition of $\partial Q_{\perp}/\partial e_{\perp}$. This quantity is not well defined; if one adds any derivative of the form $\partial F/\partial e_{\perp}$ to $\partial Q_{\perp}/\partial e_{\perp}$ with $F|_{0}^{\infty} = 0$, the total heating rate $Q_{\perp}$ remains unchanged. Indeed, Equation (A9) differs from Equation (A6) by such a term. If one were to use the definition (Equation (A7)), then the appropriate definition of the diffusion coefficient is Equation (A6).

These two methods for calculating $D_{\perp}^{E}$ require very different numerical implementations. In order to use the method of Vasquez et al. (2020), one only needs to measure the distribution function at different moments in time and then solve Equation (A8). In contrast, to use Equation (A6), one needs both $\partial Q_{\perp}/\partial e_{\perp}$ and $f^E$, but the equation for $D_{\perp}^{E}$ becomes much easier to solve.

Figure 9 shows the comparison of energy diffusion coefficients (left) and energization (right) as functions of velocity space computed using the evolution of the energy distribution function (blue) and our $E \cdot v$ diagnostic (black). The blue curve is normalized to the total heating rate, while the black curve has slightly different normalization so that the diffusion coefficient has the same magnitude in the $w_{\perp} \ll v_{th,i0}$ part of the plot. We conclude that both methods produce very similar results in the $w_{\perp} \lesssim v_{th,i0}$ part of velocity space, where our stochastic-heating theory is expected to work best.
## Appendix B
### B. Exact Calculation of $Q_{\perp}$ with Exponential Correction and Its Limits

In this Appendix, we use Equation (7) for the diffusion coefficient including the exponential correction to derive formulae for the implied perpendicular heating. We begin with Equation (7) written in terms of the potential fluctuations,

$$D_{\perp}^{E}(w_{\perp}) \sim \Omega_{i} m_{i}^{2} v_{\text{th},i}^{4} \left[ \left( \frac{w_{\perp}}{v_{\text{th},i}} \right)^{-2} \frac{q_{i}^{3} |\delta\Phi_{w}|^{3}}{m_{i}^{3} v_{\text{th},i}^{6}} \right] \exp \left[ -c_{*} \left( \frac{w_{\perp}}{v_{\text{th},i}} \right)^{2} \frac{m_{i} v_{\text{th},i}^{2}}{q_{i} \delta\Phi_{w}} \right], \quad (B1)$$

which is then substituted into the perpendicular-heating integral,

$$Q_{\perp} = - \int_{0}^{\infty} dw_{\perp} D_{\perp}^{E} \frac{\partial f^{E}}{\partial w_{\perp}}. \quad (B2)$$

We then evaluate the result in the two limits considered in Section 2.3, namely, $\beta_{\perp} \gtrsim 1$, for which the inductive electric field dominates the ion-gyroscale electrostatic potential, and $\beta_{\perp} \ll 1$, for which the ion-gyroscale fluctuations are predominantly sub-$d_i$ KAWs.

#### B.1. Stochastic Heating with Exponential Correction in $\beta \gtrsim 1$ AW Turbulence

In this limit, the electrostatic potential evaluated at perpendicular velocity $w_{\perp} \sim v_{\text{th},i}(\lambda/\rho_{\text{th},i})$ is given by

$$\delta\Phi_w \sim \rho_{\text{th},i} \left( \frac{w_{\perp}}{v_{\text{th},i}} \right) \frac{\delta u_{\perp,i,w}}{c} B_0 \sim \rho_{\text{th},i} B_0 \frac{v_{\text{th},i}}{c} \left( \frac{\varepsilon_{\text{AW}}}{\Omega_i v_A^2} \right)^{1/3} \beta_{\perp i}^{-1/3} \left( \frac{w_{\perp}}{v_{\text{th},i}} \right)^{4/3}, \quad (B3)$$

where we have used Equation (16) to rewrite $\delta\Phi_{\lambda} \rightarrow \delta\Phi_{w}$. The corresponding diffusion coefficient is then

$$D_{\perp}^{(\text{AW})}(w_{\perp}) \sim \varepsilon_{\text{AW}} m_i^2 v_{\text{th},i}^2 \left(\frac{w_{\perp}}{v_{\text{th},i}}\right)^2 \exp\left[-\mu_*^{(\text{AW})} \left(\frac{w_{\perp}}{v_{\text{th},i}}\right)^{2/3}\right], \quad (B4)$$

where

$$\mu_*^{(\text{AW})} \doteq c_* \beta_{\perp i}^{1/3} \left( \frac{\Omega_i v_A^2}{\varepsilon_{\text{AW}}} \right)^{1/3} = c_* \beta_{\perp i}^{1/2} \left( \frac{L}{\rho_{\text{th},i}} \right)^{1/3}. \quad (B5)$$

In the final step above, we have used $\varepsilon_{\text{AW}} = v_A^3/L$ to relate $\mu_*$ to the separation of scales in the system. Using $f^E(w_{\perp}) = \exp(-w_{\perp}^2/v_{\text{th},i}^2)/(m_i v_{\text{th},i}^2)$ and rewriting $\exp[-\mu_*^{(\text{AW})}(w_{\perp}/v_{\text{th},i})^{2/3}]$ using the definitions of the exponential and gamma functions, viz., $\exp(x) = \sum_{n=0}^{\infty} x^n/n!$ and $n! = \Gamma(n+1)$, respectively, we can perform the integral in Equation (B2) to determine the heating rate per unit mass of stochastic heating off of AW fluctuations, $Q_{\perp}^{(\text{AW})}$:

$$\begin{aligned} \frac{Q_{\perp}^{(\text{AW})}}{\varepsilon_{\text{AW}}} &\sim 2 \int_{0}^{\infty} \frac{dw_{\perp}}{v_{\text{th},i}} \left(\frac{w_{\perp}}{v_{\text{th},i}}\right)^3 \exp\left[-\left(\frac{w_{\perp}}{v_{\text{th},i}}\right)^2 - \mu_*^{(\text{AW})} \left(\frac{w_{\perp}}{v_{\text{th},i}}\right)^{2/3}\right] \\ &= \sum_{n=0}^{\infty} \frac{(-\mu_*^{(\text{AW})})^n}{\Gamma(n+1)} 2 \int_{0}^{\infty} dx \, x^{3+2n/3} e^{-x^2} \quad (\text{with } x \doteq w_{\perp}/v_{\text{th},i}) \\ &= \sum_{n=0}^{\infty} \frac{\Gamma(2+n/3)}{\Gamma(n+1)} (-\mu_*^{(\text{AW})})^n \Rightarrow \frac{Q_{\perp}^{(\text{AW})}}{\varepsilon_{\text{AW}}} = \Lambda_{\text{AW}} \sum_{n=0}^{\infty} \frac{\Gamma(2+n/3)}{\Gamma(n+1)} (-\mu_*^{(\text{AW})})^n. \end{aligned} \quad (B6)$$

As in Equation (18), $\Lambda_{\text{AW}}$ is a constant independent of $\beta_{\perp}$ and $\tau_{\perp}$ that takes into account the various coefficients neglected in our scaling arguments. Note that $\Gamma(2+n/3)/\Gamma(n+1)$ is a function that quickly decreases for $n > 2$; for $n = 0, 1$, and 2, its values are 1, $\simeq 1.19$, and $\simeq 0.75$. Although the result of the integral in Equation (B6) is exact, it is worth specifying its approximations in two regimes, $\mu_*^{(\text{AW})} \lesssim 1$ and $\mu_*^{(\text{AW})} \gg 1$.

$\mu_*^{(\text{AW})} \lesssim 1$ regime.—This is the regime in which stochastic heating is most efficient. The condition $\mu_*^{(\text{AW})} \lesssim 1$ is met as long as

$$\varepsilon_{\text{AW}} \gtrsim \varepsilon_{\text{crit}}^{(\text{AW})} \doteq c_*^3 \beta_{\perp i} \Omega_i v_A^2 \quad \text{or} \quad \frac{\rho_{\text{th},i}}{L} \gtrsim \chi_{\text{crit}}^{(\text{AW})} \doteq c_*^3 \beta_{\perp i}^{3/2}. \quad (B7)$$

Therefore, if the system is such that the energy injected in the Alfvénic cascade exceeds a certain critical value $\varepsilon_{\text{crit}}^{(\text{AW})}$ (or, equivalently, if the scale separation $\rho_{\text{th},i}/L$ remains above a critical value $\chi_{\text{crit}}^{(\text{AW})}$), then the dominant contributions to Equation (B6) are the $n = 0, 1$ terms. As a result,

$$\frac{Q_{\perp}^{(\text{AW})}}{\varepsilon_{\text{AW}}} \approx \Lambda_{\text{AW}} \left[ 1 - \Gamma(7/3) c_* \beta_{\perp i}^{1/2} \left( \frac{L}{\rho_{\text{th},i}} \right)^{1/3} \right]. \quad (B8)$$
The second term in brackets is a small correction to the expression (Equation (18)) obtained by neglecting the exponential suppression.

$\mu_{*}^{(AW)} \gg 1$ regime.—This is the regime in which stochastic heating is strongly suppressed for most of the ion population by the quasi-conservation of their magnetic moment. This regime holds, for instance, when the separation between the injection scale and the ion-Larmor scale in a system is significantly larger than the critical value derived above, i.e., when

$$ \frac{\rho_{\text{th,i}}}{L} \ll \chi_{\text{crit}}^{(\text{AW})} \doteq c_{*}^{3} \beta_{\perp i}^{3/2}. \qquad (\text{B9}) $$

To obtain $Q_{\perp}^{(AW)}$ in this limit, it is easier to make some approximations before performing the integral. Namely, when the exponential suppression factor is important, we may safely neglect the $\exp(-w_{\perp}^2/v_{\text{th,i}}^2)$ term in the integral (Equation (B6)). In this case, the resulting heating is

$$ \frac{Q_{\perp}^{(AW)}}{\varepsilon_{\text{AW}}} \approx \frac{\Lambda_{\text{AW}}}{(\mu_{*}^{(\text{AW})})^6} = \Lambda_{\text{AW}} c_{*}^{-6} \beta_{\perp i}^{-3} \left(\frac{\rho_{\text{th,i}}}{L}\right)^2. \qquad (\text{B10}) $$

## B.2. Stochastic Heating with Exponential Correction in Low-$\beta$ KAW Turbulence

In this limit, the electrostatic potential fluctuations may be approximated by

$$ \delta\Phi_w \approx d_i(1 + \tau_{\perp})^{-1} \frac{\delta B_{\parallel,w}}{B_0} \frac{v_A}{c} B_0 \sim \rho_{\text{th,i}} B_0 \frac{v_{\text{th,i}}}{c} \left( \frac{\varepsilon_{\text{KAW}}}{\Omega_i v_A^2} \right)^{1/3} \frac{(1 + \tau_{\perp})^{-2/3} \beta_{\perp i}^{-1/3}}{(2 + \beta_{\perp})^{1/3}} \left( \frac{w_{\perp}}{v_{\text{th,i}}} \right)^{(3+\alpha)/6}, \qquad (\text{B11}) $$

and thus

$$ D_{\perp\perp}^{(\text{KAW})}(w_{\perp}) \sim \varepsilon_{\text{KAW}} m_i^2 v_{\text{th,i}}^2 \frac{(1 + \tau_{\perp})^{-2}}{(2 + \beta_{\perp})} \left(\frac{w_{\perp}}{v_{\text{th,i}}}\right)^{(\alpha-1)/2} \exp\left[-\mu_{*}^{(\text{KAW})} \left(\frac{w_{\perp}}{v_{\text{th,i}}}\right)^{(9-\alpha)/6}\right], \qquad (\text{B12}) $$

where the parameter $\mu_{*}$ is now defined by

$$ \mu_{*}^{(\text{KAW})} \doteq c_{*} \beta_{\perp i}^{1/3} (1 + \tau_{\perp})^{2/3} (2 + \beta_{\perp})^{1/3} \left(\frac{\Omega_i v_A^2}{\varepsilon_{\text{KAW}}}\right)^{1/3} = c_{*} \beta_{\perp i}^{1/3} (1 + \tau_{\perp})^{2/3} (2 + \beta_{\perp})^{1/3} \left(\frac{\varepsilon_{\text{AW}}}{\varepsilon_{\text{KAW}}}\right)^{1/3} \left(\frac{L}{\rho_{\text{th,i}}}\right)^{1/3}. \qquad (\text{B13}) $$

We remind the reader that $1 \le \alpha \le 3$ is the parameter taking into account different models for the spectral anisotropy of the cascading KAW fluctuations (see Equation (14) in Section 2.3). Performing the $w_{\perp}$ integral of $-D_{\perp\perp}^{(\text{KAW})}(\partial f^E / \partial w_{\perp})$ and proceeding as in the AW case, we find that the heating rate per unit mass of stochastic heating off of KAW fluctuations satisfies

$$ \begin{aligned} \frac{Q_{\perp}^{(\text{KAW})}}{\varepsilon_{\text{KAW}}} &\sim 2 \int_{0}^{\infty} \frac{dw_{\perp}}{v_{\text{th,i}}} \left(\frac{w_{\perp}}{v_{\text{th,i}}}\right)^{(\alpha+1)/2} \exp\left[-\left(\frac{w_{\perp}}{v_{\text{th,i}}}\right)^2 - \mu_{*}^{(\text{KAW})} \left(\frac{w_{\perp}}{v_{\text{th,i}}}\right)^{(9-\alpha)/6}\right] \\ &= \frac{(1 + \tau_{\perp})^{-2}}{(2 + \beta_{\perp})} \sum_{n=0}^{\infty} \frac{(-\mu_{*}^{(\text{KAW})})^n}{\Gamma(n+1)} 2 \int_{0}^{\infty} dx \, x^{(\alpha+1)/2+n(9-\alpha)/6} e^{-x^2} \qquad (\text{B14}) \\ &= \frac{(1 + \tau_{\perp})^{-2}}{(2 + \beta_{\perp})} \sum_{n=0}^{\infty} \frac{\Gamma\left(\frac{3+\alpha}{4} + n\frac{9-\alpha}{12}\right)}{\Gamma(n+1)} (-\mu_{*}^{(\text{KAW})})^n \Rightarrow \frac{Q_{\perp}^{(\text{KAW})}}{\varepsilon_{\text{KAW}}} = \Lambda_{\text{KAW}} \frac{(1 + \tau_{\perp})^{-2}}{(2 + \beta_{\perp})} \sum_{n=0}^{\infty} \frac{\Gamma\left(\frac{3+\alpha}{4} + n\frac{9-\alpha}{12}\right)}{\Gamma\left(\frac{3+\alpha}{4}\right) \Gamma(n+1)} (-\mu_{*}^{(\text{KAW})})^n. \end{aligned} $$

As in Equation (22), $\Lambda_{\text{KAW}}$ is a constant independent of $\beta_{\perp i}$ and $\tau_{\perp}$ that takes into account the various coefficients neglected in our scaling arguments (note that a factor $\Gamma(\frac{3+\alpha}{4})$ has been introduced in the denominator within the sum, so that the $n=0$ term exactly matches the expression in Equation (22); this is also absorbed within the constant $\Lambda_{\text{KAW}}$). Once again, Equation (B14) is exact, but it is instructive to derive explicit analytical expressions for $Q_{\perp}^{(\text{KAW})}$ in the two interesting limits, $\mu_{*}^{(\text{KAW})} \lesssim 1$ and $\mu_{*}^{(\text{KAW})} \gg 1$.

$\mu_{*}^{(AW)} \lesssim 1$ regime.—This is the case for which the quasi-conservation of the magnetic moment does not effectively hold, making the stochastic heating of ions more effective. Such a regime occurs if the energy cascading as KAW fluctuations exceeds a critical energy cascade rate $\varepsilon_{\text{crit}}^{(\text{KAW})}$ given by

$$ \varepsilon_{\text{KAW}} \gtrsim \varepsilon_{\text{crit}}^{(\text{KAW})} \doteq c_{*}^{3} \beta_{\perp i}^{3} (1 + \tau_{\perp})^{2} (2 + \beta_{\perp}) \Omega_{i} v_{A}^{2}, \qquad (\text{B15}) $$

or, in other words, if the scale separation $\rho_{\text{th,i}}/L$ in the system remains above a critical value $\chi_{\text{crit}}^{(\text{KAW})}$ given by

$$ \frac{\rho_{\text{th,i}}}{L} \gtrsim \chi_{\text{crit}}^{(\text{KAW})} \doteq \left( \frac{\varepsilon_{\text{AW}}}{\varepsilon_{\text{KAW}}} \right) c_{*}^{3} \beta_{\perp i}^{3/2} (1 + \tau_{\perp})^{2} (2 + \beta_{\perp}). \qquad (\text{B16}) $$
Retaining only the $n=0, 1$ terms in Equation (B14), we may approximate the heating in this limit as

$$ \frac{Q_{\perp}^{(\text{KAW})}}{\varepsilon_{\text{AW}}} \approx \Lambda_{\text{KAW}} \left( \frac{\varepsilon_{\text{KAW}}}{\varepsilon_{\text{AW}}} \right) \frac{(1 + \tau_{\perp})^{-2}}{(2 + \beta_{\perp})} \left[ 1 - c_{*} \frac{\Gamma\left(\frac{9+\alpha}{6}\right)}{\Gamma\left(\frac{3+\alpha}{4}\right)} \beta_{\perp i}^{1/2} (1 + \tau_{\perp})^{2/3} (2 + \beta_{\perp})^{1/3} \left( \frac{\varepsilon_{\text{AW}}}{\varepsilon_{\text{KAW}}} \right)^{1/3} \left( \frac{L}{\rho_{\text{th},i}} \right)^{1/3} \right] \quad (\text{B17}) $$

The second term in brackets is a small correction to the expression (Equation (22)) obtained by neglecting the exponential suppression.

$\mu_{*}^{(\text{AW})} \gg 1$ regime.—Here we consider once more the regime in which the ions' magnetic moments are quasi-conserved, i.e., the regime of asymptotically weak stochastic heating from KAW fluctuations. Proceeding as in the AW case, we neglect the $\exp(-w_{\perp}^2/v_{\text{th},i}^2)$ term with respect to the suppression $\exp[-\mu_{*}^{(\text{KAW})}(w_{\perp}/v_{\text{th},i})^{(9-\alpha)/6}]$ in the integral leading to Equation (B14) and obtain the following approximate expression:

$$ \begin{aligned} \frac{Q_{\perp}^{(\text{KAW})}}{\varepsilon_{\text{AW}}} &\approx \Lambda_{\text{KAW}} \left( \frac{\varepsilon_{\text{AW}}}{\varepsilon_{\text{KAW}}} \right) \frac{(1 + \tau_{\perp})^{-2}}{(2 + \beta_{\perp})} (\mu_{*}^{(\text{KAW})})^{-3(3+\alpha)/(9-\alpha)} \\ &= \Lambda_{\text{KAW}} \left[ \left( \frac{\varepsilon_{\text{AW}}}{\varepsilon_{\text{KAW}}} \right) (1 + \tau_{\perp})^{-2} (2 + \beta_{\perp})^{-1} \right]^{12/(9-\alpha)} \left[ c_{*}^3 \beta_{\perp i}^{-3/2} \left( \frac{\rho_{\text{th},i}}{L} \right)^{(3+\alpha)/(9-\alpha)} \right]. \end{aligned} \quad (\text{B18}) $$

ORCID iDs

S. S. Cerri https://orcid.org/0000-0003-0562-6574

L. Arzamasskiy https://orcid.org/0000-0002-5263-9274

M. W. Kunz https://orcid.org/0000-0003-1676-6126

References

Adkins, T., & Schekochihin, A. A. 2018, JPIPh, 84, 905840107

Alexandrova, O., Chen, C. H. K., Sorriso-Valvo, L., Horbury, T. S., & Bale, S. D. 2013, SSRv, 178, 101

Arzamasskiy, L., Kunz, M. W., Chandran, B. D. G., & Quataert, E. 2019, ApJ, 879, 53

Barnes, A. 1966, PhFl, 9, 1483

Boldyrev, S. 2006, PhRvL, 96, 115002

Boldyrev, S., & Perez, J. C. 2012, ApJL, 758, L44

Bourouaine, S., & Chandran, B. D. G. 2013, ApJ, 774, 96

Bourouaine, S., Marsch, E., & Vocks, C. 2008, ApJL, 684, L119

Bruno, R., & Carbone, V. 2013, LRSP, 10, 2

Bruno, R., & Trenchi, L. 2014, ApJL, 787, L24

Camporeale, E., & Burgess, D. 2017, JPIPh, 83, 535830201

Cerri, S. S., Califano, F., Jenko, F., Told, D., & Rincon, F. 2016, ApJL, 822, L12

Cerri, S. S., Franci, L., Califano, F., Landi, S., & Hellinger, P. 2017a, JPIPh, 83, 705830202

Cerri, S. S., Grošelj, D., & Franci, L. 2019, FrASS, 6, 64

Cerri, S. S., Kunz, M. W., & Califano, F. 2018, ApJL, 856, L13

Cerri, S. S., Servidio, S., & Califano, F. 2017b, ApJL, 846, L18

Chandran, B. D. G. 2010, ApJ, 720, 548

Chandran, B. D. G., Dennis, T. J., Quataert, E., & Bale, S. D. 2011, ApJ, 743, 197

Chandran, B. D. G., Li, B., Rogers, B. N., Quataert, E., & Germaschewski, K. 2010, ApJ, 720, 503

Chandran, B. D. G., Verscharen, D., Quataert, E., et al. 2013, ApJ, 776, 45

Chen, C. H. K. 2016, JPIPh, 82, 535820602

Chen, C. H. K., Bale, S. D., Bonnell, J. W., et al. 2020, ApJS, 246, 53

Chen, C. H. K., & Boldyrev, S. 2017, ApJ, 842, 122

Chen, C. H. K., Klein, K. G., & Howes, G. G. 2019, NatCo, 10, 740

Chen, C. H. K., Leung, L., Boldyrev, S., Maruca, B. A., & Bale, S. D. 2014, GeoRL, 41, 8081

Chen, L., Lin, Z., & White, R. 2001, PhPl, 8, 4713

Coleman, Paul J. J. 1968, ApJ, 153, 371

Cranmer, S. R. 2014, ApJS, 213, 16

Fox, N. J., Velli, M. C., Bale, S. D., et al. 2016, SSRv, 204, 7

Franci, L., Landi, S., Verdini, A., Matteini, L., & Hellinger, P. 2018, ApJ, 853, 26

Goldreich, P., & Sridhar, S. 1995, ApJ, 438, 763

Goldstein, M. L., Roberts, D. A., & Matthaeus, W. H. 1995, ARA&A, 33, 283

Greco, A., Matthaeus, W. H., Perri, S., et al. 2018, SSRv, 214, 1

Grošelj, D., Cerri, S. S., Bañón Navarro, A., et al. 2017, ApJ, 847, 28

He, J., Wang, L., Tu, C., Marsch, E., & Zong, Q. 2015, ApJL, 800, L31

Hellinger, P., Matteini, L., Štverák, Š., Trávníček, P. M., & Marsch, E. 2011, JGRA, 116, A09105

Heuer, M., & Marsch, E. 2007, JGRA, 112, A03102

Hollweg, J. V. 1999, JGR, 104, 14811

Hoppock, I. W., Chandran, B. D. G., Klein, K. G., Mallet, A., & Verscharen, D. 2018, JPIPh, 84, 905840615

Horbury, T. S., Forman, M., & Oughton, S. 2008, PhRvL, 101, 175005

Horbury, T. S., Wicks, R. T., & Chen, C. H. K. 2012, SSRv, 172, 325

Howes, G. G. 2010, MNRAS, 409, L104

Howes, G. G. 2017, PhPl, 24, 055907

Howes, G. G., Cowley, S. C., Dorland, W., et al. 2008, JGRA, 113, A05103

Howes, G. G., Klein, K. G., & Li, T. C. 2017, JPIPh, 83, 705830102

Howes, G. G., Tenbarge, J. M., & Dorland, W. 2011, PhPl, 18, 102305

Isenberg, P. A. 2001, SSRv, 95, 119

Isenberg, P. A., & Vasquez, B. J. 2019, ApJ, 887, 63, (Erratum: 2020, ApJ, 901, 172)

Isenberg, P. A., Vasquez, B. J., & Hollweg, J. V. 2019, ApJ, 870, 119

Johnson, J. R., & Cheng, C. Z. 2001, GeoRL, 28, 4421

Kasper, J. C., Maruca, B. A., Stevens, M. L., & Zaslavsky, A. 2013, PhRvL, 110, 091102

Kawazura, Y., Barnes, M., & Schekochihin, A. A. 2019, PNAS, 116, 771

Kawazura, Y., Schekochihin, A. A., Barnes, M., et al. 2020, PhRvX, 10, 041050

Klein, K. G., & Chandran, B. D. G. 2016, ApJ, 820, 47

Klein, K. G., & Howes, G. G. 2016, ApJL, 826, L30

Klein, K. G., Howes, G. G., & Tenbarge, J. M. 2017, JPIPh, 83, 535830401

Klein, K. G., Howes, G. G., Tenbarge, J. M., & Valentini, F. 2020, JPIPh, 86, 905860402

Kunz, M. W., Abel, I. G., Klein, K. G., & Schekochihin, A. A. 2018, JPIPh, 84, 715840201

Kunz, M. W., Stone, J. M., & Bai, X.-N. 2014, JCoPh, 259, 154

Landau, L. D. 1946, ZhETF, 16, 574

Leamon, R. J., Smith, C. W., Ness, N. F., Matthaeus, W. H., & Wong, H. K. 1998, JGR, 103, 4775

Leamon, R. J., Smith, C. W., Ness, N. F., & Wong, H. K. 1999, JGR, 104, 22331

Li, T. C., Howes, G. G., Klein, K. G., Liu, Y.-H., & Tenbarge, J. M. 2019, JPIPh, 85, 905850406

Loureiro, N. L., & Boldyrev, S. 2017, ApJ, 850, 182

Lysak, R. L., & Lotko, W. 1996, JGR, 101, 5085

Mallet, A., Klein, K. G., Chandran, B. D. G., et al. 2019, JPIPh, 85, 175850302

Mallet, A., Schekochihin, A. A., & Chandran, B. D. G. 2017, MNRAS, 468, 4862

Marsch, E., Schwenn, R., Rosenbauer, H., et al. 1982, JGR, 87, 52
Marsch, E., & Tu, C. Y. 2001, JGR, 106, 227

Martinović, M. M., Klein, K. G., & Bourouaine, S. 2019, ApJ, 879, 43

Martinović, M. M., Klein, K. G., Kasper, J. C., et al. 2020, ApJS, 246, 30

Maruca, B. A., Kasper, J. C., & Bale, S. D. 2011, PhRvL, 107, 201101

Matteini, L., Landi, S., Hellinger, P., et al. 2007, GeoRL, 34, L20105

Osman, K. T., Matthaeus, W. H., Hnat, B., & Chapman, S. C. 2012, PhRvL, 108, 261103

Passot, T., & Sulem, P. L. 2015, ApJL, 812, L37

Pezzi, O., Servidio, S., Perrone, D., et al. 2018, PhPl, 25, 060704

Quataert, E. 1998, ApJ, 500, 978

Qudsi, R. A., Maruca, B. A., Matthaeus, W. H., et al. 2020, ApJS, 246, 46

Šafránková, J., Němeček, Z., Němec, F., et al. 2016, ApJ, 825, 121

Sahraoui, F., Goldstein, M. L., Belmont, G., Canu, P., & Rezeau, L. 2010, PhRvL, 105, 131101

Sahraoui, F., Goldstein, M. L., Robert, P., & Khotyaintsev, Y. V. 2009, PhRvL, 102, 231102

Sahraoui, F., Hadid, L., & Huang, S. 2020, RvMPP, 4, 4

Schekochihin, A. A., Cowley, S. C., Dorland, W., et al. 2009, ApJS, 182, 310

Servidio, S., Chasapis, A., Matthaeus, W. H., et al. 2017, PhRvL, 119, 205101

Smith, C. W., Hamilton, K., Vasquez, B. J., & Leamon, R. J. 2006, ApJL, 645, L85

TenBarge, J. M., & Howes, G. G. 2013, ApJL, 771, L27

TenBarge, J. M., Howes, G. G., & Dorland, W. 2013, ApJ, 774, 139

Told, D., Cookmeyer, J., Muller, F., Astfalk, P., & Jenko, F. 2016, NJP, 18, 065011

Told, D., Jenko, F., TenBarge, J. M., Howes, G. G., & Hammett, G. W. 2015, PhRvL, 115, 025003

Tu, C. Y., & Marsch, E. 1995, SSRv, 73, 1

Vasquez, B. J. 2015, ApJ, 806, 33

Vasquez, B. J., Isenberg, P. A., & Markovskii, S. A. 2020, ApJ, 893, 71

Vech, D., Klein, K. G., & Kasper, J. C. 2017, ApJL, 850, L11

Voitenko, Y., & Goossens, M. 2004, ApJL, 605, L149

White, R., Chen, L., & Lin, Z. 2002, PhPl, 9, 1890

Woodham, L. D., Wicks, R. T., Verscharen, D., & Owen, C. J. 2018, ApJ, 856, 49

Xia, Q., Perez, J. C., Chandran, B. D. G., & Quataert, E. 2013, ApJ, 776, 90