Hybrid-kinetic Simulations of Ion Heating in Alfvénic Turbulence

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

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

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

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

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

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

We consider a nonrelativistic, quasi-neutral, collisionless, and initially homogeneous plasma of kinetic ions (mass m_i, charge e), and massless fluid electrons threaded by a uniform magnetic field ${{\boldsymbol{B}}}_{0}={B}_{0}\hat{{\boldsymbol{z}}}$ and subject to a stochastic driving force ${\boldsymbol{F}}(t,{\boldsymbol{r}})$. The model equations governing the evolution of the ion distribution function ${f}_{{\rm{i}}}(t,{\boldsymbol{r}},{\boldsymbol{v}})$ and the magnetic field ${\boldsymbol{B}}(t,{\boldsymbol{r}})$ are, respectively, the Vlasov equation,

$$\begin{eqnarray}&&\displaystyle \frac{\partial {f}_{{\rm{i}}}}{\partial t}+{\boldsymbol{v}}\,{\boldsymbol{\cdot }}\,{\boldsymbol{\nabla }}{f}_{{\rm{i}}}+\left[\displaystyle \frac{e}{{m}_{{\rm{i}}}}\left({\boldsymbol{E}}+\displaystyle \frac{{\boldsymbol{v}}}{c}{\boldsymbol{\times }}{\boldsymbol{B}}\right)+\displaystyle \frac{{\boldsymbol{F}}}{{m}_{{\rm{i}}}}\right]\ \,{\boldsymbol{\cdot }}\,\displaystyle \frac{\partial {f}_{{\rm{i}}}}{\partial {\boldsymbol{v}}}=0,\end{eqnarray} \tag{ 1 }$$

and Faraday's law of induction,

$$\begin{eqnarray}&&\displaystyle \frac{\partial {\boldsymbol{B}}}{\partial t}=-c{\boldsymbol{\nabla }}{\boldsymbol{\times }}{\boldsymbol{E}}.\end{eqnarray} \tag{ 2 }$$

The electric field,

$$\begin{eqnarray}&&{\boldsymbol{E}}=-\displaystyle \frac{{{\boldsymbol{u}}}_{{\rm{i}}}{\boldsymbol{\times }}{\boldsymbol{B}}}{c}+\displaystyle \frac{({\boldsymbol{\nabla }}{\boldsymbol{\times }}{\boldsymbol{B}})\,{\boldsymbol{\times }}{\boldsymbol{B}}}{4\pi {{en}}_{{\rm{i}}}}-\displaystyle \frac{{T}_{{\rm{e}}}{\boldsymbol{\nabla }}{n}_{{\rm{i}}}}{{{en}}_{{\rm{i}}}},\end{eqnarray} \tag{ 3 }$$

is obtained by expanding the electron momentum equation in ${({m}_{{\rm{e}}}/{m}_{{\rm{i}}})}^{1/2}$, enforcing quasi-neutrality,

$$\begin{eqnarray}&&{n}_{{\rm{e}}}={n}_{{\rm{i}}}\equiv \ \int {d}^{3}{\boldsymbol{v}}\,{f}_{{\rm{i}}},\end{eqnarray} \tag{ 4 }$$

assuming isothermal electrons (T_e = const), and using Ampe`re's law to solve for the mean electron velocity

$$\begin{eqnarray}&&{{\boldsymbol{u}}}_{{\rm{e}}}={{\boldsymbol{u}}}_{{\rm{i}}}-\displaystyle \frac{{\boldsymbol{j}}}{{{en}}_{{\rm{i}}}}\equiv \displaystyle \frac{1}{{n}_{{\rm{i}}}}\int {d}^{3}{\boldsymbol{v}}\,{\boldsymbol{v}}{f}_{{\rm{i}}}-\displaystyle \frac{c{\boldsymbol{\nabla }}{\boldsymbol{\times }}\,{\boldsymbol{B}}}{4\pi {{en}}_{{\rm{i}}}}\end{eqnarray} \tag{ 5 }$$

in terms of the mean ion velocity ${{\boldsymbol{u}}}_{{\rm{i}}}$ and the current density ${\boldsymbol{j}}$. Equations (1)–(5) constitute the hybrid-kinetic model of Vlasov ions and (massless) fluid electrons (Byers et al. 1978; Hewett & Nielson 1978).

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

Observations of turbulence in the solar wind show that it consists of a majority of nearly incompressible, Alfvénically polarized, spatially anisotropic fluctuations (typical density variations ≲10%; Celnikier et al. 1983; Roberts et al. 1987; Marsch & Tu 1990; Tu & Marsch 1994; Bieber et al. 1996). To excite such fluctuations and minimize the excitation of compressible motions, ${\boldsymbol{F}}$ is oriented in the x–y plane perpendicular ("⊥") to ${{\boldsymbol{B}}}_{0}$ and constrained to satisfy ${\boldsymbol{\nabla }}\,{\boldsymbol{\cdot }}\,{\boldsymbol{F}}=0$. At each simulation time step, the Fourier coefficients ${{\boldsymbol{F}}}_{k}$ are independently generated from a Gaussian-random field with power spectrum ${k}_{\perp }^{-5/3}$ in the wavenumber range ${k}_{z}\in [1,2]\times \left(2\pi /{L}_{z}\right)$ and k_(x,y) ∈ [1, 2] × (2π/L_(x,y)), where L_(x,y,z) is the size of the periodic computational domain in the (x, y, z) direction. The resulting force is then inverse-Fourier transformed, shifted to ensure no net momentum injection, and normalized to provide power per unit volume $\dot{\varepsilon }$. The force is time-correlated over t_corr using an Ornstein–Uhlenbeck process:

$$\begin{eqnarray}&&{\boldsymbol{F}}(t+{dt},{\boldsymbol{r}})=\theta {\boldsymbol{F}}(t,{\boldsymbol{r}})+\sqrt{1-{\theta }^{2}}{\boldsymbol{F}}^{\prime} ({\boldsymbol{r}}),\end{eqnarray} \tag{ 6 }$$

which has the autocorrelation (in the limit ${dt}\to 0$)

$$\begin{eqnarray}&&\langle {\boldsymbol{F}}({t}_{1},{\boldsymbol{r}})\,{\boldsymbol{\cdot }}\,{\boldsymbol{F}}({t}_{2},{\boldsymbol{r}})\rangle =\langle | {\boldsymbol{F}}(t,{\boldsymbol{r}}){| }^{2}\rangle \,{e}^{({t}_{1}-{t}_{2})/{t}_{\mathrm{corr}}},\end{eqnarray} \tag{ 7 }$$

where ${dt}$ is the timestep, $\theta \equiv \exp (-{dt}/{t}_{\mathrm{corr}})$, t_corr is the correlation time of the driving, and ${\boldsymbol{F}}^{\prime}$ is a new Gaussian-random field generated as detailed above. Time-correlated driving avoids spurious particle acceleration via resonances with high-frequency power in, e.g., δ-correlated driving (Lynn et al. 2012).

Most of the power in strong MHD turbulence resides in fluctuations satisfying ${k}_{\parallel }/{k}_{\perp }\lesssim {u}_{\perp }({k}_{\perp })/{v}_{{\rm{A}}}$ ("critical balance"; Goldreich & Sridhar 1995; Mallet et al. 2015), where k_∥ (k_⊥) is the wavenumber parallel (perpendicular) to the local magnetic field and ${v}_{{\rm{A}}}\equiv B/{(4\pi {m}_{{\rm{i}}}{n}_{{\rm{i}}})}^{1/2}$ is the Alfvén speed. We therefore choose $\dot{\varepsilon }$ and L_(x,y,z) such that, in saturation, the rms velocity fluctuation u_rms satisfies ${u}_{\mathrm{rms}}/{v}_{{\rm{A}}}\simeq {L}_{x}/{L}_{z}={L}_{y}/{L}_{z}\ll 1$ at the outer scale; likewise, ${t}_{\mathrm{corr}}={L}_{z}/2\pi {v}_{{\rm{A}}}\simeq {L}_{x}/2\pi {u}_{\mathrm{rms}}$. This is meant to mimic an energy cascade from larger scales that is present in the solar wind, whose inertial-range turbulence is consistent with critical balance (Horbury et al. 2008; Podesta 2009; Luo & Wu 2010; Wicks et al. 2010).

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

We present results from two simulations: ${\beta }_{{\rm{i}}0}\equiv 8\pi {n}_{{\rm{i}}0}{T}_{{\rm{i}}0}/{B}_{0}^{2}\;=0.3$ and 1, both with T_e = T_i0 (the subscript "0" denotes an initial value). For ${\beta }_{{\rm{i}}0}=1$, N_ppc = 512 particles per cell were drawn from a Maxwell distribution and placed on a 3D periodic grid of N_x × N_y × N_z = 200² × 1600 cells spanning ${L}_{x}\times {L}_{y}\times {L}_{z}\,={(20\pi {\rho }_{{\rm{i}}0})}^{2}\times 160\pi {\rho }_{{\rm{i}}0}$, where ${v}_{\mathrm{thi}0}\equiv {(2{T}_{{\rm{i}}0}/{m}_{{\rm{i}}})}^{1/2}$ is the ion thermal speed, ${\rho }_{{\rm{i}}0}\equiv {v}_{\mathrm{thi}0}/{{\rm{\Omega }}}_{{\rm{i}}0}$ is the ion Larmor radius, and Ω_i0 ≡ eB₀/m_ic is the ion gyrofrequency. These parameters provide reasonable scale separation between the grid scale, the ion-kinetic scales, and the driving scales: ${k}_{z}{d}_{{\rm{i}}0}\in [0.0125,10]$ and ${k}_{(x,y)}{d}_{{\rm{i}}0}\in [0.1,10]$, where ${d}_{{\rm{i}}0}={\rho }_{{\rm{i}}0}/{\beta }_{{\rm{i}}0}^{1/2}$ is the ion skin depth. For ${\beta }_{{\rm{i}}0}=0.3$, N_ppc = 216, N_x × N_y × N_z = 200² × 1200, and ${L}_{x}\times {L}_{y}\times {L}_{z}={(20\pi {\rho }_{{\rm{i}}0})}^{2}\,\times \,120\pi {\rho }_{{\rm{i}}0}$. These imply ${k}_{z}{d}_{{\rm{i}}0}\in [0.03,18.26]$ and ${k}_{(x,y)}{d}_{{\rm{i}}0}\in [0.18,18.26]$.

With these scales borne in mind, it is useful to compare with the observed spectral anisotropy in the near-Earth solar wind. There, spectral anisotropy of the turbulent cascade begins around kρ_i ∼ 10⁻³ and increases as ${k}_{\parallel }/{k}_{\perp }\propto {k}_{\perp }^{-1/3}$ toward smaller scales (Wicks et al. 2010). Adopting these scalings implies k_∥/k_⊥ ≃ 0.17 around k_⊥ρ_i ≃ 0.2, where our simulated inertial range begins. Thus, our chosen box aspect ratios of 1/8 (for ${\beta }_{{\rm{i}}0}=1$) and 1/6 (for ${\beta }_{{\rm{i}}0}=0.3$) accurately capture the spectral anisotropy of solar-wind turbulence arriving from larger scales to near ion-Larmor scales. (This point is revisited near the end of Section 3.1 in the context of the observed spectral anisotropy at ion-Larmor scales.)

Both simulations required at least $2{L}_{z}/{v}_{{\rm{A}}0}\equiv 2{t}_{\mathrm{cross}}$ to obtain a quasi-steady state, in which the properties of the turbulence exhibit only minimal secular evolution. An additional ∼10t_cross were run to procure statistically converged results (viz., 10t_cross for ${\beta }_{{\rm{i}}0}=1$ and 18t_cross for ${\beta }_{{\rm{i}}0}=0.3$). Note that the particle distribution function never reaches a true steady state, as there is no cooling and so the total particle energy grows monotonically in time. However, the changes in ion temperature at the end of both simulations are small enough to affect neither the properties of the turbulence nor the heating diagnostics measured in these simulations. In what follows, $\langle \,\cdot \,\rangle$ denotes a spatio-temporal average over all cells measured in this quasi-steady state.

An example of the turbulent quasi-steady state is given in Figure 1, which shows pseudo-color images of the x- and y-components (i.e., those perpendicular to the mean field) of the fluctuating magnetic field from the ${\beta }_{{\rm{i}}0}=1$ run. Spatial anisotropy is evident, with short perpendicular scales and long parallel scales consistent with a critically balanced cascade (i.e., δB_rms/B₀ ≈ 1/8, the box aspect ratio). This anisotropy plays an important role in shaping the energy spectra (Section 3.1) and the nature of ion heating (Section 3.2).

Zoom In Zoom Out Reset image size

3.1. Energy Spectra and Spectral Anisotropy

Figures 2(a), (b) present energy spectra of magnetic-field (${{ \mathcal E }}_{B}$) and electric-field (${{ \mathcal E }}_{E}$) fluctuations for ${\beta }_{{\rm{i}}0}=0.3$ and β_i0 = 1 versus the wavenumber k_⊥ perpendicular to ${{\boldsymbol{B}}}_{0}$. Their scale-dependent spectral indices α(k_⊥), computed about each k_⊥ using 21 neighboring points (except for the first 5 points for which we use 11 neighboring harmonics), are shown in Figures 2(c), (d). Accompanying these in Figures 2(e), (f) are the ion-flow-velocity (${{ \mathcal E }}_{u}$), density (${{ \mathcal E }}_{\widetilde{n}}$), and parallel-magnetic-field (${{ \mathcal E }}_{{B}_{\parallel }}$) energy spectra.

Zoom In Zoom Out Reset image size

In the inertial range (${k}_{\perp }{\rho }_{{\rm{i}}0}\lesssim 1$), the spectral slopes are close to −5/3, the spectral slope predicted for an anisotropic, critically balanced cascade of Alfvén waves (Goldreich & Sridhar 1995). A spectral break occurs at ${k}_{\perp }{\rho }_{{\rm{i}}0}\sim 1$, at which point ${{ \mathcal E }}_{E}$ flattens to take on a slope near −2/3, in agreement with measurements that show a flattening of the solar-wind electric-field power spectrum in the ion-kinetic range (Bale et al. 2005; Sahraoui et al. 2009; Salem et al. 2012) and predictions for intermittent KAW turbulence (Boldyrev & Perez 2012). In this range, ${{ \mathcal E }}_{B}$ also steepens to take on a slope comparable to the −2.8 commonly observed in the near-Earth solar wind (e.g., Alexandrova et al. 2009, 2012; Sahraoui et al. 2009) and within the range [−2.5, −3.1] found in Cluster spacecraft measurements of the β ∼ 1 solar wind (Sahraoui et al. 2013). However, this spectral steepening appears to be k_⊥-dependent, an attribute we revisit below when discussing ion heating. For now, we remark that such steepening is consistent with the sub-ion spectrum measured in Earth's magnetosheath (Chen et al. 2019). Despite −2.8 being steeper than predictions for standard (−7/3; Schekochihin et al. 2009) and intermittent (−8/3; Boldyrev & Perez 2012) KAW turbulence, the normalized perturbed density $\delta \widetilde{n}\equiv (\delta {n}_{{\rm{i}}}/{n}_{{\rm{i}}0})\sqrt{{\beta }_{{\rm{i}}0}(1+{\beta }_{{\rm{i}}0})}$ approximately follows the linear KAW eigenfunction $\delta \widetilde{n}=\delta {B}_{\perp }/{B}_{0}=(\delta {B}_{\parallel }/{B}_{0})\sqrt{1+1/{\beta }_{{\rm{i}}0}}$ (for T_e0 = T_i0; Schekochihin et al. 2009; Boldyrev & Perez 2013). This suggests that the sub-ion-Larmor-scale cascade is primarily composed of KAWs (at least up to ${k}_{\perp }{\rho }_{{\rm{i}}0}\approx 3$, where ${{ \mathcal E }}_{\widetilde{n}}$ starts to be affected by particle noise and digital filters), in agreement with combined analyses of magnetic-field fluctuations and of electric-field or density fluctuations in the ion-kinetic range of solar-wind turbulence (e.g., Salem et al. 2012; Chen et al. 2013). That being said, there is slightly more (about a factor ≲2) magnetic-fluctuation energy than anticipated for KAWs, suggesting the presence of additional wave modes (an excess of ≃1.33 was measured in the solar wind by Chen et al. 2013). At ${k}_{\perp }{\rho }_{{\rm{i}}0}\approx 6$, the spectra steepen further due to hyper-resistivity and spectral filters (ion heating is less important at these scales—see Section 3.2).

The predicted slopes of −5/3 in the inertial range and −7/3 (or −8/3) in the KAW range follow from assuming locality of interactions and constant energy flux in Fourier space, combined with a model for the spatial anisotropy of the turbulent fluctuations, k_∥(k_⊥). The latter is afforded by the critical balance assumption, which states that the scale-dependent nonlinear cascade time is comparable to the linear timescale of the dominant wave at that scale (Goldreich & Sridhar 1995)—essentially a causality argument (Boldyrev 2005). For the inertial range, in which the characteristic linear frequency is that of Alfvén waves, ω_AW = k_∥v_A, a perpendicular spectrum of −5/3 corresponds to ${k}_{\parallel }\propto {k}_{\perp }^{2/3}$. For the KAW range, in which the characteristic linear frequency is ${\omega }_{\mathrm{KAW}}={k}_{\parallel }{v}_{{\rm{A}}}{k}_{\perp }{\rho }_{{\rm{i}}}/\sqrt{1+{\beta }_{{\rm{i}}}}$, a perpendicular spectrum of −7/3 (−8/3) corresponds to ${k}_{\parallel }\propto {k}_{\perp }^{1/3}$ (${k}_{\parallel }\propto {k}_{\,\perp }^{2/3}$), the steeper parenthetical values reflecting the two-dimensionalization of intermittent KAW turbulence (Boldyrev & Perez 2012).

In Figures 2(g), (h), we show the spectral anisotropy computed from our simulations.⁵ In both runs, ${k}_{\parallel }\propto {k}_{\perp }^{2/3}$ in the inertial range. At the start of the KAW range, the wavevector anisotropy ${({k}_{\perp }/{k}_{\parallel })}_{\mathrm{KAW}}\approx 15$ and 12 for ${\beta }_{{\rm{i}}0}=0.3$ and 1, respectively, corresponding to wavevector obliquities ${\theta }_{(k,B)}\equiv {\tan }^{-1}({k}_{\perp }/{k}_{\parallel })\approx 86^\circ$ and 85°. These values are comparable to those measured in the solar wind (e.g., Sahraoui et al. 2010; Narita et al. 2011). Well into the sub-ion-Larmor range, however, we observe k_∥ ∝ k_⊥, steeper than current theoretical predictions and corresponding to a scale-independent anisotropy. Constant spectral anisotropy in the sub-ion-Larmor-scale range has also been seen in other hybrid-kinetic simulations of 3D Alfvénic turbulence (Franci et al. 2018).

Constant spectral anisotropy at ${k}_{\perp }{\rho }_{{\rm{i}}0}\gt 1$ implies that the turbulence there is more likely to attain ion-Larmor frequencies before reaching electron scales, since ${\omega }_{\mathrm{KAW}}\propto {k}_{\,\perp }^{2}$, rather than the standard ${\omega }_{\mathrm{KAW}}\propto {k}_{\,\perp }^{4/3}$ predicted by theories of low-frequency gyrokinetic turbulence (e.g., Schekochihin et al. 2009). Such high-frequency fluctuations facilitate additional energy-transfer channels that are not present in gyrokinetics (see Howes et al. 2008a), such as cyclotron resonances and high-frequency stochastic heating, to which we now turn.

3.2. Ion Heating

In this section we examine the turbulent heating of ions in our simulations. In doing so, we are careful to distinguish between particle energization, by which we mean the rate at which work is done on the particles by the electric field (viz., $Q\equiv {\boldsymbol{v}}\,{\boldsymbol{\cdot }}\,{\boldsymbol{E}}$) temporally averaged over long time intervals, and particle heating, by which we mean an increase in the Maxwellian temperature of the distribution function. This distinction is important. For example, while transit-time damping (TTD) is related to the work done by fluctuating perpendicular electric fields on Landau-resonant particles (i.e., perpendicular energization, ${Q}_{\perp }\equiv {{\boldsymbol{v}}}_{\perp }\,{\boldsymbol{\cdot }}\,{{\boldsymbol{E}}}_{\perp }$), it results in an increase in the parallel temperature T_∥. Likewise, changes in T_∥ need not be driven by particle energization; for example, such changes may occur via pitch-angle scattering of perpendicular particle energy into the parallel direction (as, indeed, is shown later in the paper).

Figure 3 shows time-averaged parallel (Q_∥) and perpendicular (Q_⊥) particle-energization rates as functions of ${k}_{\perp }{\rho }_{{\rm{i}}0}$. These are calculated by using Fourier transforms to spatially filter the electric field into 12 logarithmically spaced k_⊥-bins, giving ${\boldsymbol{E}}({k}_{\perp ,\mathrm{bin}})$, and then summing ${v}_{\parallel }{E}_{\parallel }({k}_{\perp ,\mathrm{bin}})$ and ${{\boldsymbol{v}}}_{\perp }\,{\boldsymbol{\cdot }}\,{{\boldsymbol{E}}}_{\perp }({k}_{\perp ,\mathrm{bin}})$ over all particles for each bin (the first and last bins are not shown). (The use of Fourier transforms results in "⊥" and "∥" being measured here with respect to the guide field.) Rates are normalized to the rate of energy injection by the large-scale forcing, Q_inj. In both runs, Q_⊥ ≫ Q_∥. For ${\beta }_{{\rm{i}}0}=1$, about 80% of injected energy goes into the ions; the remaining 20% is removed by hyper-resistivity. For β_i0 = 0.3, this ratio is roughly 75%:25%. In both cases, Q_⊥ peaks at ${k}_{\perp }{\rho }_{{\rm{i}}0}\approx 4$, approximately where the measured spectral anisotropies in the KAW range imply ${\omega }_{\mathrm{KAW}}\approx {{\rm{\Omega }}}_{{\rm{i}}0}$ (marked by the dotted lines in Figures 2(g), 2(h), and 3).⁶ This suggests that the majority of the energization is produced when ions move in the oscillating potential of high-frequency KAWs, a possibility discussed further in Section 4. For now, we note that the sub-ion-Larmor spectrum of magnetic-field fluctuations seen in Figures 2(a), (b) steepens beyond the anticipated ≈−2.8 at the values of ${k}_{\perp }{\rho }_{{\rm{i}}0}$ for which the perpendicular energization is largest—likely not a coincidence.

Zoom In Zoom Out Reset image size

In Figure 4, we present the velocity- and real-space tracks of two particles from the ${\beta }_{{\rm{i}}0}=1$ run, one ("1") that exhibits appreciable perpendicular energization and another ("2") that does not. Upper panels show snapshots of small-scale (${k}_{\perp }{\rho }_{{\rm{i}}0}\geqslant 1$) magnetic (left) and electric (right) field fluctuations in the x–y plane located at the z-coordinate of both tracked particles. The black lines trace projected particle trajectories over four gyro-periods. Particle 1 resides in a region of large-amplitude, sheet-like structures, whereas particle 2 samples only weak field fluctuations. (Note that the fluctuations experience the same ${\boldsymbol{E}}\,{\boldsymbol{\times }}\,{\boldsymbol{B}}$ drift motion as the particles, and so particle 1's guiding center is not actually sweeping across the small-scale fluctuations.) Lower panels show the evolution of the particles' peculiar ("thermal") velocity ${\boldsymbol{w}}\equiv {\boldsymbol{v}}-{{\boldsymbol{u}}}_{{\rm{i}}}$ over the final $1000{{\rm{\Omega }}}_{{\rm{i}}0}^{-1}$ of the run. The gray shaded region marks the time interval over which particle trajectories are plotted. At each point, thermal velocities are averaged over 4π/Ω_i0 to remove high-frequency oscillations in energy; these oscillations are shown for ${w}_{\mathrm{tot}}^{2}$ with the thin green line.

Zoom In Zoom Out Reset image size

Particle 1 sees an appreciable increase in its energy, with ${\rm{\Delta }}{w}_{\mathrm{tot}}^{2}\approx 0.3{v}_{\mathrm{thi}0}^{2}$ over only four gyro-periods. Most of this increase occurs perpendicularly to the local magnetic field. On the other hand, the energy of particle 2 stays nearly constant. Note further that the field-parallel and field-perpendicular energies do not grow monotonically. Often, they are subject to strong kicks during which the total energy is almost constant but the pitch angle of the particle changes dramatically, a feature we refer to as pitch-angle scattering. As a result, high-w_⊥ particles scatter and subsequently contribute to wings produced in the parallel distribution function, a feature we return to below in the context of Figures 6–9.

Figure 5 shows the time-averaged parallel (perpendicular) differential energization, ${{dQ}}_{\parallel (\perp )}/{dw}$, in the gyrotropic velocity space (w_∥, w_⊥). (Namely, ${dQ}/{dw}$ is the particle energization per interval of velocity space; these rates are normalized in the figure to the total energization rate Q_tot in each run. Here, "∥" and "⊥" are measured with respect to the local magnetic field at the position of the particle.) Figure 5 is to be read alongside Figure 6, which displays the parallel and perpendicular ion distribution functions measured at the ends of both runs, viz., $f({w}_{\parallel })\equiv \int {{dw}}_{\perp }^{2}\,f({w}_{\parallel },{w}_{\perp })$ and $f({w}_{\perp })\equiv \int {{dw}}_{\parallel }\,f({w}_{\parallel },{w}_{\perp })$, respectively (again, measured with respect to the local magnetic-field direction). The dashed lines in Figure 6 show best-fit Maxwellians to the core of f(w_∥) (with temperature ${T}_{\parallel }^{\mathrm{core}}$) and to the tail of f(w_⊥) (with temperature ${T}_{\perp }^{\mathrm{tail}}$). Both runs exhibit the following attributes: (i) resonant features in the particle-energization rates with fine-scale structure near ${w}_{\parallel }\sim {v}_{\mathrm{thi}0}$ and $\sim {v}_{{\rm{A}}0}$ (Figures 5(a), (c)); (ii) quasi-linear flattening of the parallel distribution at $| {w}_{\parallel }| /{v}_{{\rm{A}}0}\sim 1$ and nonthermal wings at $| {w}_{\parallel }| /{v}_{\mathrm{thi}0}\gtrsim 1$ (Figures 6(a), (c)); (iii) almost no change in the parallel temperature of the core ($| {w}_{\parallel }| /{v}_{\mathrm{thi}0}\lt 1;$ Figures 6(a), (c)), with very little parallel energization there (Figures 5(a), (c)); (iv) a broadened perpendicular distribution for $1\lesssim {w}_{\perp }/{v}_{\mathrm{thi}0}\lesssim 3$ (Figures 6(b), (d)), where the perpendicular energization peaks (Figures 5(b), (d)); and (v) a flattening of the core of the perpendicular distribution f(w_⊥) (Figures 6(b), (d)), with suppressed perpendicular energization at ${w}_{\perp }/{v}_{\mathrm{thi}0}\lesssim 1$ (Figures 5(b), (d)). Regarding this final point, we find evidence early in each run for perpendicular energization at ${w}_{\perp }/{v}_{\mathrm{thi}0}\lesssim 1$, which ultimately causes the observed flattening of the core of f(w_⊥) (see Figure 8 and the accompanying discussion below).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

All of these features can also be seen in the full 2D (gyrotropic) distribution function f(w_∥, w_⊥) shown in Figure 7, where the differences between the two runs are even more striking. In particular, non-Maxwellian features such as the flattened w_⊥ core and the parallel beams at ${w}_{\parallel }\sim {v}_{{\rm{A}}0}$ are readily apparent, with the latter driving $\partial f/\partial {w}_{\parallel }\approx 0$ for ${w}_{\perp }/{v}_{\mathrm{thi}0}\gtrsim 3$ and $| {w}_{\parallel }| \lesssim {v}_{{\rm{A}}0}$.

Zoom In Zoom Out Reset image size

We associate with these non-Maxwellian features a number of physical effects. The early-time flattening of f(w_⊥) for ${w}_{\perp }/{v}_{\mathrm{thi}0}\lesssim 1$ is most likely due to stochastic heating by low-frequency, Larmor-scale fluctuations, following the prediction by Klein & Chandran (2016) that such heating leads to a flattop distribution (see also Johnson & Cheng 2001). As stochastic heating is expected to be more important at lower β for fixed $\delta {B}_{\perp }/{B}_{0}$ (e.g., Chandran et al. 2010; Hoppock et al. 2018), it is noteworthy that $f({w}_{\perp })$ exhibits a larger amount of flattening for ${\beta }_{{\rm{i}}0}=0.3$ than for ${\beta }_{{\rm{i}}0}=1$. The broadening of the thermal tail of $f({w}_{\perp })$ is instead driven by the peak in ${{dQ}}_{\perp }/{{dw}}_{\perp }$ at ${w}_{\perp }/{v}_{\mathrm{thi}0}\approx 1.6$, which constitutes the bulk of the perpendicular heating and is related to the peak in ${dQ}/d\mathrm{log}{k}_{\perp }$ at ${k}_{\perp }{\rho }_{{\rm{i}}0}\approx 4\mbox{--}5$ (Figure 3).

This interpretation is further strengthened by examining the perpendicular-energy diffusion coefficient ${D}_{\perp \perp }^{{\rm{E}}}$. Assuming that the distribution function evolves according to a Fokker–Planck-like equation,

$$\begin{eqnarray}&&\displaystyle \frac{\partial f}{\partial t}=\displaystyle \frac{\partial }{\partial {e}_{\perp }}\left({D}_{\perp \perp }^{{\rm{E}}}\displaystyle \frac{\partial f}{\partial {e}_{\perp }}\right),\end{eqnarray} \tag{ 8 }$$

where ${e}_{\perp }\equiv {w}_{\perp }^{2}/2$ is the field-perpendicular kinetic energy per particle, the energization rate

$$\begin{eqnarray}\begin{array}{rcl}{Q}_{\perp } & \equiv & \displaystyle \frac{\partial }{\partial t}\displaystyle \int {{de}}_{\perp }\,{e}_{\perp }f({e}_{\perp })=\displaystyle \int {{de}}_{\perp }\,{e}_{\perp }\displaystyle \frac{\partial }{\partial {e}_{\perp }}\left({D}_{\perp \perp }^{{\rm{E}}}\displaystyle \frac{\partial f}{\partial {e}_{\perp }}\right)\\ & = & -\displaystyle \int {{de}}_{\perp }\,{D}_{\perp \perp }^{{\rm{E}}}\displaystyle \frac{\partial f}{\partial {e}_{\perp }},\end{array}\end{eqnarray} \tag{ 9 }$$

where the final equality is obtained after integrating by parts. The box-averaged perpendicular-energy diffusion coefficient can thus be computed from the derivatives of the perpendicular energization and the distribution function, as follows:

$$\begin{eqnarray}&&\langle {D}_{\perp \perp }^{{\rm{E}}}\rangle \equiv -\left\langle \displaystyle \frac{\partial {Q}_{\perp }}{\partial {e}_{\perp }}\right\rangle /\left\langle \displaystyle \frac{\partial f}{\partial {e}_{\perp }}\right\rangle .\end{eqnarray} \tag{ 10 }$$

This quantity is plotted as a function of ${w}_{\perp }/{v}_{\mathrm{thi}0}$ in Figure 8 both at early times in the ${\beta }_{{\rm{i}}0}=0.3$ run (panel (a)) and at late times in both runs (panel (b)). Accompanying these panels are (c), (d) the differential perpendicular particle energization (as in Figure 5) and (e), (f) the perpendicular distribution function (as in Figure 6) at those times. In the latter panels, the evolution of f(w_⊥) from its initial condition (black dotted line) to its profile at the beginning (solid line) and the end (dashed line) of the time-averaging window are shown. At early times, the diffusion coefficient and the accompanying ${{dQ}}_{\perp }/{{dw}}_{\perp }$ exhibit two distinct peaks. We associate the peak at ${w}_{\perp }/{v}_{\mathrm{thi}0}\lt 1$ with low-frequency stochastic heating, which flattens the core of f(w_⊥) by accelerating particles to larger w_⊥. The ${w}_{\perp }/{v}_{\mathrm{thi}0}\gt 1$ peak in ${{dQ}}_{\perp }/{{dw}}_{\perp }$ is associated with a diffusion coefficient that scales approximately as ${w}_{\,\perp }^{4}$. At late times in the ${\beta }_{{\rm{i}}0}=0.3$ run (red lines), after the core of the perpendicular distribution has been flattened, the diffusion coefficient is roughly constant with w_⊥ and very little energization happens for ${w}_{\perp }/{v}_{\mathrm{thi}0}\lt 1$ (note that $\langle {D}_{\perp \perp }^{{\rm{E}}}\rangle$ is not well-defined in this limit, as both $\partial {Q}_{\perp }/\partial {w}_{\perp }^{2}$ and $\partial f/\partial {w}_{\perp }^{2}$ are close to zero). At larger velocities, the perpendicular energization peaks (as in Figure 5) and the diffusion coefficient continues to scale approximately as ${w}_{\,\perp }^{4}$. In the ${\beta }_{{\rm{i}}0}=1$ run, which did not experience as great of a flattening in the core of f(w_⊥), $\langle {D}_{\perp \perp }^{{\rm{E}}}\rangle \propto {w}_{\perp }$ appears to be a good approximation for ${w}_{\perp }/{v}_{\mathrm{thi}0}\lt 1$.

Zoom In Zoom Out Reset image size

Most of the injected energy goes into the increase of perpendicular temperature and into the development of nonthermal tails in the parallel distribution function. (Overall, ≈60% of the total energization Q_tot ultimately finds its way into the nonthermal w_∥ wings, with the remaining ≈40% going into raising the perpendicular temperature.) Given that most of the energization is perpendicular, this strongly suggests that TTD and pitch-angle scattering of super-thermal ${w}_{\perp }^{2}$ into ${w}_{\parallel }^{2}$ are the mechanisms responsible for the non-Maxwellian features seen in $f({w}_{\parallel })$.

We tested these possibilities by examining the energization and pitch angles of 160,000 individually tracked particles in the ${\beta }_{{\rm{i}}0}=1$ run. We found that the mirror force $\mu \,{dB}({{\boldsymbol{x}}}_{{\rm{p}}})/{dt}$ (where $\mu \equiv {w}_{\perp }^{2}/2B$ is the magnetic moment) is responsible for ≲20% of Q_⊥, a value consistent with the quasi-linear flattening of the distribution function observed at ${w}_{\parallel }\sim {v}_{{\rm{A}}0}$ in Figures 6 and 7. The remaining (small) amount of increase in T_∥ can be accounted for by Landau damping (i.e., Q_∥ = v_∥E_∥ ≲ 0.1Q_tot). The other ≳80% of Q_⊥ is due to some other mechanism that leads to perpendicular heating at ${w}_{\perp }\gt {v}_{\mathrm{thi}0}$ (see Section 4 for a discussion of possibilities). To examine the effect of pitch-angle scattering on the distribution function, we divided each particle track into time intervals of 5 · 2π/Ω_i0 and computed ${w}_{\perp }^{2}/{w}_{\parallel }^{2}$ at the beginning and the end of each interval. We then filtered each interval based on whether or not ${w}_{\perp }^{2}/{w}_{\parallel }^{2}$ changes (up or down) by a certain minimal factor that we call "threshold." Large values of threshold represent large changes in a particle's pitch angle; small values of threshold include time intervals during which the pitch angle is almost constant. Figure 9 shows $\overline{{\rm{\Delta }}{w}_{(\parallel ,\perp ,\mathrm{tot})}^{2}/{\rm{\Delta }}t}$, the average⁷ rate of change of the parallel (blue), perpendicular (red), and total (green) thermal energies of tracked particles, segregated into super- and subthermal populations (${w}_{\perp }\gt {v}_{\mathrm{thi}0}$ and ${w}_{\perp }\lt {v}_{\mathrm{thi}0}$, respectively). Given a threshold value, Figure 9 provides the expected energization rates for an individual particle due to all events whose change in ${w}_{\perp }^{2}/{w}_{\parallel }^{2}$ are above that threshold. As the length of the time interval used to subdivide the particle tracks is somewhat arbitrary, we indicate with the shaded regions the variation of particle-energization rates for time intervals from 4 · 2π/Ω_i0 to 8 · 2π/Ω_i0. Note that $| \overline{{\rm{\Delta }}{w}_{\parallel ,\perp }^{2}/{\rm{\Delta }}t}|$ is a decreasing function of threshold, as it represents the cumulative energization for all events above the threshold. The energization per event (not shown) is an increasing function of threshold. For ${w}_{\perp }\lt {v}_{\mathrm{thi}0}$, there is a net conversion of thermal energy from parallel to perpendicular, consistent with the other diagnostics shown in Figures 5–8. For ${w}_{\perp }\gt {v}_{\mathrm{thi}0}$, however, the flow of thermal energy between parallel and perpendicular is reversed and, for thresholds ≳1.5, this flow takes place at constant total energy. This strongly suggests that pitch-angle scattering is responsible for this transfer of super-thermal perpendicular energy into super-thermal parallel energy (similar to recent work by Isenberg et al. 2019), ultimately producing the nonthermal wings seen in Figures 6 and 7.⁸

Zoom In Zoom Out Reset image size

We reemphasize that, while ${T}_{\perp }\gt {T}_{\parallel }^{\mathrm{core}}$ at the end of our simulations, more than half of the cascade energy ultimately goes into the development of the nonthermal wings in the parallel distribution function. To illustrate that, the blue line in Figure 10 traces the temporal evolution of the ratio of the total parallel and perpendicular energies of the particles, $\langle {v}_{\parallel }^{2}\rangle /\langle {v}_{\perp }^{2}/2\rangle$, from the ${\beta }_{{\rm{i}}0}=0.3$ run in which the nonthermal wings are most pronounced (the factor of 1/2 accounts for the number of degrees of freedom in the perpendicular direction, so that an isotropic Maxwellian has $\langle {v}_{\parallel }^{2}\rangle /\langle {v}_{\perp }^{2}/2\rangle =1$). The accompanying red line represents the ratio of best-fit Maxwellian temperatures to the core of the parallel distribution function (${T}_{\parallel }^{\mathrm{core}}$) and to the tail of the perpendicular distribution function (${T}_{\perp }^{\mathrm{tail}}\simeq {T}_{\perp }$). The initial drop in $\langle {v}_{\parallel }^{2}\rangle /\langle {v}_{\perp }^{2}/2\rangle$ is caused by the increase in perpendicular bulk motion in the initial stage of the simulation. (Recall that ${\boldsymbol{v}}$ denotes the total (bulk + thermal) velocity of the ion particles.) Once the turbulence obtains quasi-steady state (for t/t_cross ≳ 3), the perpendicular temperature steadily grows larger than the parallel temperature of the core. This is mirrored in the evolution of $\langle {v}_{\parallel }^{2}\rangle /\langle {v}_{\perp }^{2}/2\rangle$, at least until t/t_cross ≈ 8, after which the ratio of energies suddenly begins to increase, eventually becoming larger than 1 at t/t_cross ≈ 16. The distinction between the ratio of energies and the ratio of best-fit-Maxwellian temperatures is thus an important one. Indeed, while nonthermal v_∥-wings are often observed in ion velocity distribution functions in the solar wind, the "observed" T_∥ is often determined by a bi-Maxwellian fit to the core of the distribution (e.g., Bame et al. 1975; Marsch et al. 1982b).

Zoom In Zoom Out Reset image size

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

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

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

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

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

Given that complementary gyrokinetic theory and simulations show a predominance of electron heating over ion heating for β ≲ 1 (Howes 2010; Howes et al. 2011; Told et al. 2015; Navarro et al. 2016; Kawazura et al. 2019), it is worthwhile to contemplate which framework is better suited for describing near-Earth solar-wind turbulence (see Howes et al. 2008a, Section 3 for arguments in favor of a gyrokinetic description). Gyrokinetic theory is built on the assumption of asymptotically low-frequency, small-amplitude, and spatially anisotropic fluctuations; the accompanying computational savings affords a magnetokinetic description of a realistic hydrogenic plasma without prohibitive cost. However, the magnetic moment is an invariant in the gyrokinetic equations (in the absence of explicit collisions), and so perpendicular particle heating is not allowed. By contrast, the hybrid model makes no such simplifying assumptions, but saves computational expense by neglecting electron kinetics. While the latter precludes a truly rigorous study of ion versus electron heating, it is worth noting that there is empirical evidence for a majority fraction (≳60%) of ion versus electron heating between 0.3 and 5 au in the solar wind (Cranmer et al. 2009); recall that ≈75%–80% of our cascade energy goes into heating the ions (the rest is dissipated via hyper-resistivity). Moreover, the ion-Larmor-scale spectral anisotropy in the β_i ∼ 1 solar wind is not necessarily asymptotically small. Indeed, ${\theta }_{(k,B)}\approx 80^\circ \mbox{--}90^\circ$ for f_spacecraft ∼ 1 Hz fluctuations measured in the near-Earth solar wind (Sahraoui et al. 2010); our values are similar, ${\theta }_{(k,B)}\approx 85^\circ \mbox{--}86^\circ$. If k_∥ scales with k_⊥ to some power larger than the 1/3 adopted by Howes et al. (2008a) for k_⊥ρ_i0 > 1, as it does in our simulations and as is predicted in theories of intermittent KAW turbulence (Boldyrev & Perez 2012), it becomes all the more probable that the ion-cyclotron frequency will be attained in a KAW cascade. These considerations favor the hybrid-kinetic description over the gyrokinetic one.

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

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

L.A. and M.W.K. were supported by the National Aeronautics and Space Administration (NASA) under grant No. NNX16AK09G issued through the Heliophysics Supporting Research Program. Additional support for M.W.K. was provided by an Alfred P. Sloan Research Fellowship in Physics and, during the early stages of this project (c. 2014), by a Lyman Spitzer, Jr. Fellowship. B.D.G.C. was supported by NASA grants NNN06AA01C, NNX15AI80, NNX16AG81G, and NNX17AI18G, and by the National Science Foundation (NSF) under grant PHY-1500041. E.Q. was supported in part by a Simons Investigator award from the Simons Foundation. High-performance computing resources were provided by the Texas Advanced Computer Center at The University of Texas at Austin under grant numbers TG-AST140011 and TG-AST130058; the NASA High-End Computing (HEC) Program through the NASA Advanced Supercomputing (NAS) Division at Ames Research Center; and the PICSciE-OIT TIGRESS High Performance Computing Center and Visualization Laboratory at Princeton University. This work used the Extreme Science and Engineering Discovery Environment (XSEDE), which is supported by NSF grant OCI-1053575. This work benefited from useful conversations with Christopher Chen, Alfred Mallet, Robert Wicks, Vladimir Zhdankin, and especially Silvio Sergio Cerri, Gregory Howes, Philip Isenberg, Kristopher Klein, and Alexander Schekochihin.