Hybrid-Kinetic Simulations of Ion Heating in Alfv\'{e}nic Turbulence

TL;DR

This paper uses hybrid-kinetic simulations to study ion heating mechanisms in Alfvénic turbulence relevant to the near-Earth solar wind, revealing how energy cascades and particle distributions evolve at sub-ion scales.

Contribution

It introduces novel diagnostics and detailed analysis of ion heating processes, especially near the ion-Larmor scale, highlighting the role of kinetic-Alfvén waves and non-thermal distribution features.

Findings

75-80% of turbulent energy heats ions

Spectral steepening beyond -2.8 in sub-ion-Larmor range

Non-Maxwellian velocity distribution wings observed

Abstract

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

Equations

$$\frac{\partial f_{\rm i}}{\partial t}+\mbox{\boldmath{$v$}}\,\mbox{\boldmath{$\cdot$}}\,\mbox{\boldmath{$\nabla$}}f_{\rm i}+\left[\frac{e}{m_{\rm i}}\left(\mbox{\boldmath{$E$}}+\frac{\mbox{\boldmath{$v$}}}{c}\,\mbox{\boldmath{$\times$}}\,\mbox{\boldmath{$B$}}\right)+\frac{\mbox{\boldmath{$F$}}}{m_{\rm i}}\right]\!\,\mbox{\boldmath{$\cdot$}}\,\frac{\partial f_{\rm i}}{\partial\mbox{\boldmath{$v$}}}=0,$$

$$\frac{\partial\mbox{\boldmath{$B$}}}{\partial t}=-c\mbox{\boldmath{$\nabla$}}\,\mbox{\boldmath{$\times$}}\,\mbox{\boldmath{$E$}}.$$

$$\mbox{\boldmath{$E$}}=-\frac{\mbox{\boldmath{$u$}}_{\rm i}\,\mbox{\boldmath{$\times$}}\,\mbox{\boldmath{$B$}}}{c}+\frac{(\mbox{\boldmath{$\nabla$}}\,\mbox{\boldmath{$\times$}}\,\mbox{\boldmath{$B$}})\,\mbox{\boldmath{$\times$}}\,\mbox{\boldmath{$B$}}}{4\pi en_{\rm i}}-\frac{T_{\rm e}\mbox{\boldmath{$\nabla$}}n_{\rm i}}{en_{\rm i}},$$

$$n_{\rm e}=n_{\rm i}\equiv\!\int{\rm d}^{3}\mbox{\boldmath{$v$}}\,f_{\rm i},$$

$$\mbox{\boldmath{$u$}}_{\rm e}=\mbox{\boldmath{$u$}}_{\rm i}-\frac{\mbox{\boldmath{$j$}}}{en_{\rm i}}\equiv\frac{1}{n_{\rm i}}\int{\rm d}^{3}\mbox{\boldmath{$v$}}\,\mbox{\boldmath{$v$}}f_{\rm i}-\frac{c\mbox{\boldmath{$\nabla$}}\,\mbox{\boldmath{$\times$}}\,\mbox{\boldmath{$B$}}}{4\pi en_{\rm i}}$$

$$\mbox{\boldmath{$F$}}(t+{\rm d}t,\mbox{\boldmath{$r$}})=\theta\mbox{\boldmath{$F$}}(t,\mbox{\boldmath{$r$}})+\sqrt{1-\theta^{2}}\mbox{\boldmath{$F$}}^{\prime}(\mbox{\boldmath{$r$}}),$$

$$\langle\mbox{\boldmath{$F$}}(t_{1},\mbox{\boldmath{$r$}})\,\mbox{\boldmath{$\cdot$}}\,\mbox{\boldmath{$F$}}(t_{2},\mbox{\boldmath{$r$}})\rangle=\langle|\mbox{\boldmath{$F$}}(t,\mbox{\boldmath{$r$}})|^{2}\rangle\,e^{(t_{1}-t_{2})/t_{\rm corr}},$$

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

$$Q_{\perp }$$

$$=-\int \mathrm{d}{e}_{\perp }{D}_{\perp \perp }^{\mathrm{E}}\frac{\partial f}{\partial e_{\perp }},$$

$$\langle D^{\rm E}_{\perp\perp}\rangle\equiv-\left\langle\frac{\partial Q_{\perp}}{\partial e_{\perp}}\right\rangle\bigg{/}\left\langle\frac{\partial f}{\partial e_{\perp}}\right\rangle.$$

Full text

1. Introduction.

It has been 48 years since NASA’s Mariner 5 established definitively that the interplanetary medium plays host to a broadband spectrum of large-amplitude Alfvén waves propagating outwards from the Sun (Belcher & Davis 1971, following pioneering work using Mariner 2 data by Coleman 1968). We now know that the solar wind is turbulent (Tu & Marsch, 1995; Goldstein et al., 1995), exhibiting a power spectrum extending over several decades in scale (Bruno & Carbone, 2005; Alexandrova et al., 2013). Most of the energy at large scales is in the form of Alfvénic fluctuations, which have magnetic and velocity fields perpendicular to the mean magnetic-field direction. As this energy cascades down to smaller scales, an inertial range is set up, one whose defining characteristic is anisotropy with respect to the field direction (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 (Sahraoui et al., 2010; Chen et al., 2010), where the turbulence is ultimately dissipated as heat.

The amplitudes of turbulent fluctuations measured in the solar wind are positively correlated with solar-wind temperature (Grappin et al., 1990) and imply a turbulent heating rate comparable to the observationally inferred solar-wind heating rate (Smith et al., 2001; Breech et al., 2009; Cranmer et al., 2009; Stawarz et al., 2009). While these observations establish a close connection between the global evolution of the solar wind and the dissipation of turbulence within it, the precise nature of this dissipation is puzzling. Minor ions in coronal holes and protons in low-beta fast-solar-wind streams are heated in such a way that thermal motions perpendicular to the background magnetic field are more rapid than thermal motions along it (i.e., $T_{\perp}>T_{\parallel}$; 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~{}{\rm au}$ to $0.9~{}{\rm au}$ clearly indicates non-adiabatic particle heating preferentially in the field-perpendicular direction (see Fig. 1 of Matteini et al. 2007).

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

An alternative explanation for the measured strong perpendicular heating of ions is that of low-frequency stochastic heating, which arises when particles interact with turbulent fluctuations whose characteristic frequencies are much smaller than the cyclotron frequency but whose amplitudes at the Larmor scale (i.e., $k_{\perp}\rho_{\rm i}\sim 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 (Johnson & Cheng, 2001; Chen et al., 2001; White et al., 2002; Voitenko & Goossens, 2004; Bourouaine et al., 2008; Chandran et al., 2010). In the context of solar-wind turbulence, low-frequency stochastic heating has been studied numerically using test particles in randomly phased kinetic Alfvén waves (KAWs) (Chandran et al., 2010) and in reduced-MHD turbulence (Xia et al., 2013), and observationally in Helios-2 and Wind data (Chandran et al., 2013; Bourouaine & Chandran, 2013; Vech et al., 2017). It was also the focus of work by Vasquez (2015), who used hybrid-kinetic simulations of decaying turbulence in a low-beta, cold-electron plasma to find that the perpendicular heating rate scales with the cube of the turbulence amplitude evaluated at the ion-Larmor scale (as predicted by Chandran et al. 2010, but without their multiplicative factor that exponentially suppresses stochastic heating at small enough turbulence amplitude).

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

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

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

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

[TABLE]

and Faraday’s law of induction,

[TABLE]

The electric field,

[TABLE]

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

[TABLE]

assuming isothermal electrons ($T_{\rm e}={\rm const}$), and using Ampére’s law to solve for the mean electron velocity

[TABLE]

in terms of the mean ion velocity \mbox{\boldmath{u}}_{\rm i} and the current density $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 $\delta f$ method is used to reduce the impact of finite-particle-number noise on computed moments of $f_{\rm i}$. To remove small-scale power from the fluctuating fields, the zeroth ($n_{\rm i}$) and first (n_{\rm i}\mbox{\boldmath{u}}_{\rm i}) moments of $f_{\rm 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 $\lesssim 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, $F$ is oriented in the $x$-$y$ plane perpendicular (“$\perp$”) to \mbox{\boldmath{B}}_{0} and constrained to satisfy \mbox{\boldmath{\nabla}}\,\mbox{\boldmath{\cdot}}\,\mbox{\boldmath{F}}=0. At each simulation time step, the Fourier coefficients \mbox{\boldmath{F}}_{\boldsymbol{k}} are independently generated from a Gaussian-random field with power spectrum $k^{-5/3}_{\perp}$ in the wavenumber range $k_{z}\in[1,2]\times\left(2\pi/L_{z}\right)$ and $k_{(x,y)}\in[1,2]\times(2\pi/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_{\rm corr}$ using an Ornstein–Uhlenbeck process:

[TABLE]

which has the autocorrelation (in the limit ${\rm d}t\rightarrow 0$)

[TABLE]

where ${\rm d}t$ is the timestep, $\theta\equiv\exp(-{\rm d}t/t_{\rm corr})$, $t_{\rm corr}$ is the correlation time of the driving, and \mbox{\boldmath{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., $\delta$-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_{\parallel}$ ($k_{\perp}$) 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_{\rm rms}$ satisfies $u_{\rm rms}/v_{\rm A}\simeq L_{x}/L_{z}=L_{y}/L_{z}\ll 1$ at the outer scale; likewise, $t_{\rm corr}=L_{z}/2\pi v_{\rm A}\simeq L_{x}/2\pi u_{\rm 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 i0}\equiv 8\pi n_{{\rm i}0}T_{{\rm i}0}/B^{2}_{0}=0.3$ and $1$, both with $T_{\rm e}=T_{{\rm i}0}$ (the subscript “0” denotes an initial value). For $\beta_{\rm i0}=1$, $N_{\rm ppc}=512$ particles per cell were drawn from a Maxwell distribution and placed on a 3D periodic grid of $N_{x}\times N_{y}\times N_{z}=200^{2}\times 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_{{\rm thi}0}\equiv(2T_{{\rm i}0}/m_{\rm i})^{1/2}$ is the ion thermal speed, $\rho_{{\rm i}0}\equiv v_{{\rm thi}0}/\Omega_{{\rm i}0}$ is the ion Larmor radius, and $\Omega_{{\rm i}0}\equiv eB_{0}/m_{\rm i}c$ 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 i0}^{1/2}$ is the ion skin depth. For $\beta_{\rm i0}=0.3$, $N_{\rm ppc}=216$, $N_{x}\times N_{y}\times N_{z}=200^{2}\times 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\rho_{\rm i}\sim 10^{-3}$ and increases as $k_{\parallel}/k_{\perp}\propto k^{-1/3}_{\perp}$ towards smaller scales (Wicks et al., 2010). Adopting these scalings implies $k_{\parallel}/k_{\perp}\simeq 0.17$ around $k_{\perp}\rho_{\rm i}\simeq 0.2$, where our simulated inertial range begins. Thus, our chosen box aspect ratios of $1/8$ (for $\beta_{\rm i0}=1$) and $1/6$ (for $\beta_{\rm i0}=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 §3.1 in the context of the observed spectral anisotropy at ion-Larmor scales.)

Both simulations required at least $2L_{z}/v_{{\rm A}0}\equiv 2t_{\rm cross}$ to obtain a quasi-steady state, in which the properties of the turbulence exhibit only minimal secular evolution. An additional ${\sim}10t_{\rm cross}$ were run to procure statistically converged results (viz., $10t_{\rm cross}$ for $\beta_{\rm i0}=1$ and $18t_{\rm cross}$ for $\beta_{\rm i0}=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.

3. Results

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 i0}=1$ run. Spatial anisotropy is evident, with short perpendicular scales and long parallel scales consistent with a critically balanced cascade (i.e., $\delta B_{\rm rms}/B_{0}\approx 1/8$, the box aspect ratio). This anisotropy plays an important role in shaping the energy spectra (§3.1) and the nature of ion heating (§3.2).

3.1. Energy spectra and spectral anisotropy

Figures 2(a,b) present energy spectra of magnetic-field ($\mathcal{E}_{B}$) and electric-field ($\mathcal{E}_{E}$) fluctuations for $\beta_{\rm i0}=0.3$ and $\beta_{\rm i0}=1$ versus the wavenumber $k_{\perp}$ perpendicular to \mbox{\boldmath{B}}_{0}. Their scale-dependent spectral indices $\alpha(k_{\perp})$, computed about each $k_{\perp}$ 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.

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 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., Sahraoui et al., 2009; Alexandrova et al., 2009, 2012) and within the range $[-2.5,-3.1]$ found in Cluster spacecraft measurements of the $\beta\sim 1$ solar wind (Sahraoui et al., 2013). However, this spectral steepening appears to be $k_{\perp}$-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 i0})\sqrt{\beta_{\rm i0}(1+\beta_{\rm i0})}$ approximately follows the linear KAW eigenfunction $\delta\widetilde{n}=\delta B_{\perp}/B_{0}=(\delta B_{\parallel}/B_{0})\sqrt{1+1/\beta_{\rm i0}}$ (for $T_{\rm e0}=T_{\rm 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 ${\lesssim}2$) magnetic-fluctuation energy than anticipated for KAWs, suggesting the presence of additional wave modes (an excess of ${\simeq}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 §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_{\parallel}(k_{\perp})$. 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, $\omega_{\rm AW}=k_{\parallel}v_{\rm A}$, a perpendicular spectrum of $-5/3$ corresponds to $k_{\parallel}\propto k^{2/3}_{\perp}$. For the KAW range, in which the characteristic linear frequency is $\omega_{\rm 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^{1/3}_{\perp}$ ($k_{\parallel}\propto k^{2/3}_{\perp}$), the steeper parenthetical values reflecting the two-dimensionalization of intermittent KAW turbulence (Boldyrev & Perez, 2012).

In Figure 2(g,h), we show the spectral anisotropy computed from our simulations.111To compute the spectral anisotropy, we use equation (34) of Cho & Lazarian (2009) to compute the scale-dependent wavenumber parallel to the local mean magnetic field, $k_{\parallel}(k_{\perp})$; those authors showed that this method works very well for steep spectra. We also computed two- and three-point second-order structure functions and measured the spectral anisotropy from their isocontours (similar to Cho & Vishniac 2000), finding similar results. In both runs, $k_{\parallel}\propto k^{2/3}_{\perp}$ in the inertial range. At the start of the KAW range, the wavevector anisotropy $(k_{\perp}/k_{\parallel})_{\rm KAW}\approx 15$ and $12$ for $\beta_{\rm i0}=0.3$ and $1$, respectively, corresponding to wavevector obliquities $\theta_{(\boldsymbol{k},\boldsymbol{B})}\equiv\tan^{-1}(k_{\perp}/k_{\parallel})\approx 86^{\circ}$ and $85^{\circ}$. 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_{\parallel}\propto k_{\perp}$, 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}>1$ implies that the turbulence there is more likely to attain ion-Larmor frequencies before reaching electron scales, since $\omega_{\rm KAW}\propto k^{2}_{\perp}$, rather than the standard $\omega_{\rm KAW}\propto k^{4/3}_{\perp}$ predicted by theories of low-frequency gyrokinetic turbulence (e.g., Schekochihin et al., 2009). Such high-frequency fluctuations facilitate additional energy-transfer channels that are not present in gyrokinetics (cf. Howes et al., 2008a), such as cyclotron resonances and high-frequency stochastic heating, to which we now turn.

3.2. Ion heating

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\mbox{\boldmath{v}}\,\mbox{\boldmath{\cdot}}\,\mbox{\boldmath{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\mbox{\boldmath{v}}_{\perp}\,\mbox{\boldmath{\cdot}}\,\mbox{\boldmath{E}}_{\perp}), it results in an increase in the parallel temperature $T_{\parallel}$. Likewise, changes in $T_{\parallel}$ 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_{\parallel}$) and perpendicular ($Q_{\perp}$) 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_{\perp}$-bins, giving \mbox{\boldmath{E}}(k_{\perp,{\rm bin}}), and then summing $v_{\parallel}E_{\parallel}(k_{\perp,{\rm bin}})$ and \mbox{\boldmath{v}}_{\perp}\,\mbox{\boldmath{\cdot}}\,\mbox{\boldmath{E}}_{\perp}(k_{\perp,{\rm bin}}) over all particles for each bin (the first and last bins are not shown). (The use of Fourier transforms results in “$\perp$” and “$\parallel$” being measured here with respect to the guide field.) Rates are normalized to the rate of energy injection by the large-scale forcing, $Q_{\rm inj}$. In both runs, $Q_{\perp}\gg Q_{\parallel}$. For $\beta_{\rm i0}=1$, about 80% of injected energy goes into the ions; the remaining 20% is removed by hyper-resistivity. For $\beta_{\rm i0}=0.3$, this ratio is roughly $75\%:25\%$. In both cases, $Q_{\perp}$ peaks at $k_{\perp}\rho_{{\rm i}0}\approx 4$, approximately where the measured spectral anisotropies in the KAW range imply $\omega_{\rm KAW}\approx\Omega_{\rm i0}$ (marked by the dotted lines in Figures 2(g,h) and 3).222These values were calculated using the approximation $\omega_{\rm KAW}=k_{\parallel}v_{\rm A}k_{\perp}\rho_{\rm i}/\sqrt{1+\beta_{\rm i}}$ and confirmed using the numerical linear hybrid-kinetic solver HYDROS (Told et al., 2016a, cf. their figure 8). This suggests that the majority of the energization is produced when ions move in the oscillating potential of high-frequency KAWs, a possibility discussed further in §4. For now, we note that the sub-ion-Larmor spectrum of magnetic-field fluctuations seen in Fig. 2(a,b) steepens beyond the anticipated ${\approx}-2.8$ at the values of $k_{\perp}\rho_{{\rm i}0}$ for which the perpendicular energization is largest—likely not a coincidence.

In Figure 4, we present the velocity- and real-space tracks of two particles from the $\beta_{\rm i0}=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}\geq 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 4 gyro-periods. Particle 1 resides in a region of large-amplitude, sheet-like structures, whereas particle 2 samples only weak field fluctuations. (Note that the fluctuations experience the same $E$ $\times$ $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 \mbox{\boldmath{w}}\equiv\mbox{\boldmath{v}}-\mbox{\boldmath{u}}_{\rm i} over the final $1000\Omega_{\rm i0}^{-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\pi/\Omega_{\rm i0}$ to remove high-frequency oscillations in energy; these oscillations are shown for $w_{\rm tot}^{2}$ with the thin green line.

Particle 1 sees an appreciable increase in its energy, with $\Delta w_{\rm tot}^{2}\approx 0.3v_{{\rm thi}0}^{2}$ over only 4 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_{\perp}$ 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, ${\rm d}Q_{\parallel(\perp)}/{\rm d}w$, in the gyrotropic velocity space $(w_{\parallel},w_{\perp})$. (Namely, ${\rm d}Q/{\rm d}w$ is the particle energization per interval of velocity space; these rates are normalized in the figure to the total energization rate $Q_{\rm tot}$ in each run. Here, “$\parallel$” and “$\perp$” 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{\rm d}w^{2}_{\perp}\,f(w_{\parallel},w_{\perp})$ and $f(w_{\perp})\equiv\int{\rm d}w_{\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_{\parallel})$ (with temperature $T_{\parallel}^{\rm core}$) and to the tail of $f(w_{\perp})$ (with temperature $T_{\perp}^{\rm tail}$). Both runs exhibit the following attributes: (i) resonant features in the particle-energization rates with fine-scale structure near $w_{\parallel}\sim v_{{\rm thi}0}$ and ${\sim}v_{{\rm A}0}$ (Figs 5(a),(c)); (ii) quasi-linear flattening of the parallel distribution at $|w_{\parallel}|/v_{{\rm A}0}\sim 1$ and non-thermal wings at $|w_{\parallel}|/v_{{\rm thi}0}\gtrsim 1$ (Figs 6(a),(c)); (iii) almost no change in the parallel temperature of the core ($|w_{\parallel}|/v_{{\rm thi}0}<1$; Figs 6(a),(c)), with very little parallel energization there (Figs 5(a),(c)); (iv) a broadened perpendicular distribution for $1\lesssim w_{\perp}/v_{{\rm thi}0}\lesssim 3$ (Figs 6(b),(d)), where the perpendicular energization peaks (Figs 5(b),(d)); and (v) a flattening of the core of the perpendicular distribution $f(w_{\perp})$ (Figs 6(b),(d)), with suppressed perpendicular energization at $w_{\perp}/v_{{\rm thi}0}\lesssim 1$ (Figs 5(b),(d)). Regarding this final point, we find evidence early in each run for perpendicular energization at $w_{\perp}/v_{{\rm thi}0}\lesssim 1$, which ultimately causes the observed flattening the core of $f(w_{\perp})$ (see Fig. 8 and accompanying discussion below).

All of these features can also be seen in the full 2D (gyrotropic) distribution function $f(w_{\parallel},w_{\perp})$ 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_{\perp}$ 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_{{\rm thi}0}\gtrsim 3$ and $|w_{\parallel}|\lesssim v_{{\rm A}0}$.

We associate with these non-Maxwellian features a number of physical effects. The early-time flattening of $f(w_{\perp})$ for $w_{\perp}/v_{{\rm 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 $\beta$ 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 i0}=0.3$ than for $\beta_{\rm i0}=1$. The broadening of the thermal tail of $f(w_{\perp})$ is instead driven by the peak in ${\rm d}Q_{\perp}/{\rm d}w_{\perp}$ at $w_{\perp}/v_{{\rm thi}0}\approx 1.6$, which constitutes the bulk of the perpendicular heating and is related to the peak in ${\rm d}Q/{\rm d}\log k_{\perp}$ at $k_{\perp}\rho_{{\rm i}0}\approx 4$–$5$ (Fig. 3).

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

[TABLE]

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

[TABLE]

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:

[TABLE]

This quantity is plotted as a function of $w_{\perp}/v_{{\rm thi}0}$ in Figure 8 both at early times in the $\beta_{\rm i0}=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_{\perp})$ 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 ${\rm d}Q_{\perp}/{\rm d}w_{\perp}$ exhibit two distinct peaks. We associate the peak at $w_{\perp}/v_{{\rm thi}0}<1$ with low-frequency stochastic heating, which flattens the core of $f(w_{\perp})$ by accelerating particles to larger $w_{\perp}$. The $w_{\perp}/v_{{\rm thi}0}>1$ peak in ${\rm d}Q_{\perp}/{\rm d}w_{\perp}$ is associated with a diffusion coefficient that scales approximately as $w^{4}_{\perp}$. At late times in the $\beta_{\rm i0}=0.3$ run (red lines), after the core of the perpendicular distribution has been flattened, the diffusion coefficient is roughly constant with $w_{\perp}$ and very little energization happens for $w_{\perp}/v_{{\rm thi}0}<1$ (note that $\langle D^{\rm E}_{\perp\perp}\rangle$ is not well-defined in this limit, as both $\partial Q_{\perp}/\partial w^{2}_{\perp}$ and $\partial f/\partial w^{2}_{\perp}$ 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^{4}_{\perp}$. In the $\beta_{\rm i0}=1$ run, which did not experience as great of a flattening in the core of $f(w_{\perp})$, $\langle D^{\rm E}_{\perp\perp}\rangle\propto w_{\perp}$ appears to be a good approximation for $w_{\perp}/v_{{\rm thi}0}<1$.

Most of the injected energy goes into the increase of perpendicular temperature and into the development of non-thermal tails in the parallel distribution function. (Overall, ${\approx}60\%$ of the total energization $Q_{\rm tot}$ ultimately finds its way into the non-thermal $v_{\parallel}$ wings, with the remaining ${\approx}40\%$ going into raising the perpendicular temperature.) Given that most of the energization is perpendicular, this strongly suggests that transit-time damping (TTD) and pitch-angle scattering of super-thermal $w^{2}_{\perp}$ into $w^{2}_{\parallel}$ 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 i0}=1$ run. We found that the mirror force \mu\,{\rm d}B(\mbox{\boldmath{x}}_{\rm p})/{\rm d}t (where $\mu\equiv w^{2}_{\perp}/2B$ is the magnetic moment) is responsible for ${\lesssim}20\%$ of $Q_{\perp}$, 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_{\parallel}$ can be accounted for by Landau damping (i.e., $Q_{\parallel}=v_{\parallel}E_{\parallel}\lesssim 0.1Q_{\rm tot}$). The other ${\gtrsim}80\%$ of $Q_{\perp}$ is due to some other mechanism that leads to perpendicular heating at $w_{\perp}>v_{{\rm thi}0}$ (see §4 for a discussion of possibilities). To examine the effect of pitch-angle scattering on the distribution function, we divided each particle track into time intervals of $5\cdot 2\pi/\Omega_{\rm i0}$ and computed $w^{2}_{\perp}/w^{2}_{\parallel}$ at the beginning and the end of each interval. We then filtered each interval based on whether or not $w^{2}_{\perp}/w^{2}_{\parallel}$ 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{\Delta w^{2}_{(\parallel,\perp,{\rm tot})}/\Delta t}$, the average333To compute $\overline{\Delta w_{(\parallel,\perp)}^{2}/\Delta t}$ as a function of threshold, we sum $\Delta w_{\parallel}^{2}$ and $\Delta w_{\perp}^{2}$ for all events that change $w^{2}_{\perp}/w^{2}_{\parallel}$ by a factor larger than threshold, and then divide this energy increment by the total time over which we examine the particle tracks ($\Delta t=2000\Omega_{\rm i0}^{-1}$) and by the total number of tracked particles. rate of change of the parallel (blue), perpendicular (red), and total (green) thermal energies of tracked particles, segregated into super- and sub-thermal populations ($w_{\perp}>v_{{\rm thi}0}$ and $w_{\perp}<v_{{\rm 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^{2}_{\perp}/w^{2}_{\parallel}$ are above that threshold. As the length of the time interval used to sub-divide the particle tracks is somewhat arbitrary, we indicate with the shaded regions the variation of particle-energization rates for time intervals from $4\cdot 2\pi/\Omega_{\rm i0}$ to $8\cdot 2\pi/\Omega_{\rm i0}$. Note that $|\overline{\Delta w^{2}_{\parallel,\perp}/\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}<v_{{\rm 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}>v_{{\rm thi}0}$, however, the flow of thermal energy between parallel and perpendicular is reversed and, for thresholds ${\gtrsim}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 non-thermal wings seen in Figures 6 and 7.444Evidence of pitch-angle scattering can also be seen in the evolution of particle 1 at $\Omega_{\rm i0}t\approx 5400$–$5600$ in Figure 4; and there is the additional circumstantial evidence that the more pronounced non-thermal wings in $f(w_{\parallel})$ seen in the $\beta_{\rm i0}=0.3$ run (as compared to the $\beta_{\rm i0}=1$ run) coincides with a larger perpendicular temperature measured in the tail of $f(w_{\perp})$.

We re-emphasize that, while $T_{\perp}>T_{\parallel}^{\rm core}$ at the end of our simulations, more than half of the cascade energy ultimately goes into the development of the non-thermal wings in the parallel distribution function. To illustrate that, the blue line in Figure 10 traces the temporal evolution of the ratio of the total parallel and perpendicular energies of the particles, $\langle v_{\parallel}^{2}\rangle/\langle v_{\perp}^{2}/2\rangle$, from the $\beta_{\rm i0}=0.3$ run in which the non-thermal 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}^{\rm core}$) and to the tail of the perpendicular distribution function ($T_{\perp}^{\rm 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 $v$ denotes the total (bulk + thermal) velocity of the ion particles.) Once the turbulence obtains quasi-steady state (for $t/t_{\rm cross}\gtrsim 3$), the perpendicular temperature steadily grows larger than the parallel temperature of the core. This is mirrored in the evolution of $\langle v^{2}_{\parallel}\rangle/\langle v^{2}_{\perp}/2\rangle$, at least until $t/t_{\rm cross}\approx 8$, after which the ratio of energies suddenly begins to increase, eventually becoming larger than 1 at $t/t_{\rm cross}\approx 16$. The distinction between the ratio of energies and the ratio of best-fit-Maxwellian temperatures is thus an important one. Indeed, while non-thermal $v_{\parallel}$-wings are often observed in ion velocity distribution functions in the solar wind, the “observed” $T_{\parallel}$ is often determined by a bi-Maxwellian fit to the core of the distribution (e.g., Bame et al., 1975; Marsch et al., 1982b).

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

The perpendicular energization that occurs at early times for $w_{\perp}\lesssim v_{{\rm 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 \mbox{\boldmath{v}}_{\perp}\,\mbox{\boldmath{\cdot}}\,\mbox{\boldmath{E}}_{\perp} energization, with lower $\beta_{\rm i0}$ 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_{{\rm 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_{{\rm 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^{\rm E}_{\perp\perp}\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~{}{\rm Hz}$, which Podesta (2009) associated with the strong linear dissipation of KAWs occurring at $k_{\perp}\rho_{{\rm i}0}\approx 4$ when $k_{\perp}/k_{\parallel}=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.

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

We have presented 3D numerical simulations of driven, quasi-steady-state, hybrid-kinetic turbulence in a magnetized plasma of relevance to the $\beta_{\rm i}\lesssim 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_{\perp}\rho_{\rm i}\sim 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_{\parallel})$, a preponderance of perpendicular heating (over parallel heating), and a flattened core in $f(w_{\perp})$.

Given that complementary gyrokinetic theory and simulations show a predominance of electron heating over ion heating for $\beta\lesssim 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, §3) for arguments in favor of a gyrokinetic description). Gyrokinetic theory is built on the assumption of asymptotically low-frequency, small-amplitude, spatially anisotropic fluctuations; the accompanying computational savings affords a magneto-kinetic 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 ($\gtrsim$60%) of ion versus electron heating between $0.3$ and $5~{}{\rm au}$ in the solar wind (Cranmer et al., 2009); recall that ${\approx}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 $\beta_{\rm i}\sim 1$ solar wind is not necessarily asymptotically small. Indeed, $\theta_{(\boldsymbol{k},\boldsymbol{B})}\approx 80^{\circ}$–$90^{\circ}$ for $f_{\rm spacecraft}\sim 1~{}{\rm Hz}$ fluctuations measured in the near-Earth solar wind (Sahraoui et al., 2010); our values are similar, $\theta_{(\boldsymbol{k},\boldsymbol{B})}\approx 85^{\circ}$–$86^{\circ}$. If $k_{\parallel}$ scales with $k_{\perp}$ to some power larger than the $1/3$ adopted by Howes et al. (2008a) for $k_{\perp}\rho_{\rm 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 “sub-grid” prescriptions are based either on gyrokinetic theory and simulation (Howes, 2010; Told et al., 2015), which predict preferential electron (ion) heating at low (high) plasma beta driven by Landau-resonant damping, or on models of collisionless reconnection (Rowan et al., 2017; Werner et al., 2018; Rowan et al., 2019), in which the amount of ion versus electron heating is sensitive to the presence of a strong guide field. It would be interesting to study the effect of the non-Landau-resonant processes presented in this paper, which predominantly heat the ions, on the imaging and evolution of collisionless accretion flows. However, the application of our results to these systems is not straightforward. While the wavevector anisotropies in our simulations are consistent with the scale separation observed in the solar wind (viz., a factor of ${\sim}10^{4}$ 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 ${\sim}10^{7}$ 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 ${\sim}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.