We argue that two prominent theories of ion heating in lowbeta collisionless plasmas  stochastic and quasilinear heating  represent similar physical processes in turbulence with different normalized cross helicities. To capture both we propose a simple phenomenology based on the power in scales at which critically balanced fluctuations reach their smallest parallel scale. Simulations of test ions interacting with turbulence confirm our scalings across a wide range of different ion and turbulence properties including with a steep ionkinetic transition range as relevant to the solar wind.
['Astrophysics - Solar and Stellar Astrophysics', 'Physics - Plasma Physics', 'Physics - Space Physics']
A Unified Phenomenology of Ion Heating in Lowbeta Plasmas TestParticle Simulations
{'A Unified Phenomenology of Ion Heating in Lowβ Plasmas: Test-Particle Simulations': "Zade Johnston, ∗ Jonathan Squire, and Romain Meyrand Physics Department, University of Otago, Dunedin 9010, New Zealand  \n(Dated: September 12, 2024)  \nWe argue that two prominent theories of ion heating in lowβ collisionless plasmas-stochastic and quasi-linear heating-represent similar physical processes in turbulence with different normalized cross helicities. To capture both, we propose a simple phenomenology based on the power in scales at which critically balanced fluctuations reach their smallest parallel scale. Simulations of test ions interacting with turbulence confirm our scalings across a wide range of different ion and turbulence properties, including with a steep ion-kinetic transition range as relevant to the solar wind.  \nIons within the turbulent lowβ solar corona and wind are observed to be continously heated, with more energization perpendicular to the local magnetic field than parallel [1-3]. Additionally, heavy ions are heated more than protons [4-6] and protons more than electrons [7], despite the low-frequency gyrokinetic model predicting that electron heating dominates at low β [8, 9]. Many mechanisms have been proposed to explain these observations, including stochastic heating by uncorrelated turbulent fluctuations [10, 11], wave-particle interactions leading to quasi-linear heating [12-14], and others [1517]. Despite this previous work, aspects of the underlying processes remain unclear, representing an important deficiency for building understanding of coronal heating and solar-wind acceleration.  \nExtending previous studies [18-23], in this Letter we present a unified description of turbulent heating of ions in both balanced turbulence (with equal amplitudes of inwards and outwards Elsasser fluctuations z ± ) and imbalanced turbulence (where one dominates: | z + | ≫ | z -| ). Noting that stochastic heating (SH) and cyclotron-resonant quasi-linear heating (QLH) can be considered as two limits of similar physical processes and motivated by the principle of 'critical balance', we propose that the ion-heating rate is controlled by the power in scales at which turbulent fluctuations reach their maximum frequency, which may not be at ion-gyroradius scales in the presence of a steep ion-kinetic 'transitionrange' drop in the spectrum. We test this proposal with high-resolution simulations of test particles interacting with lowβ magnetized turbulence. Our proposed scaling accurately captures measured proton and minor-ion heating rates in turbulence across a wide range of conditions (balanced, imbalanced, with and without a 'helicity barrier'). These results clarify the physics of collisionless turbulent heating with clear application to both in-situ observational analyses and theoretical models of astrophysical plasmas [24].  \nBefore continuing we define important symbols, with species index s either ions or protons ( s = i or p): E and B are the electric and magnetic fields; k ∥ and k ⊥ are the components of the wavevector k parallel and perpendicular to B ; ρ s = v th ,s / Ω s is the thermal gyroradius with  \nv th ,s = √ 2 T s /m s and Ω s = | q s | B/m s c the ion thermal speed and gyrofrequency; q s = Ze and m s = Am p are the ion charge and mass in units of proton quantities; β s = v 2 th ,s /v 2 A is the ion beta; c is the speed of light; δu λ is the rms E × B velocity at scales k ⊥ λ ∼ 1; and v A is the Alfv'en speed.  \nTurbulence Theory. -Turbulence near the Sun is imbalanced, with | z + | ≫ | z -| . Recent work [25-27] argues that such turbulence forms a 'helicity barrier', which, due to its effect on turbulent fluctuations near ρ p -scales, is a key ingredient in understanding ion heating. The effect occurs because lowβ gyrokinetics conserves energy and a 'generalized helicity' [25, 28], which must be simultaneously cascaded with constant flux to be in a statistical steady-state. At scales larger than ρ p , helicity is the cross-helicity which cascades to small perpendicular scales. At scales smaller than ρ p , helicity becomes magnetic helicity, which would inverse cascade in imbalanced turbulence. For balanced turbulence, with zero helicity, energy is able to cascade to small perpendicular scales. However, in imbalanced turbulence only the balanced (zero-helicity) portion of the energy cascade is let through to small perpendicular scales; the remainder is blocked from cascading and trapped at scales above ρ p by the helicity barrier. This causes the amplitude of the turbulence to grow in time and, through critical balance [29], reach and dissipate at small parallel scales. Hybridkinetic imbalanced turbulence simulations [26] show that fluctuations at these scales become oblique ion-cyclotron waves, heating ions through quasi-linear interactions and allowing the turbulence to saturate. The helicity barrier is well supported observationally, explaining the steep 'transition-range' drop in fluctuation spectra [30, 31] and other features [32, 33].  \nIon Heating. -Two widely studied ion-heating mechanisms are stochastic heating and cyclotron-resonant quasi-linear diffusion. Stochastic heating involves ions being heated by uncorrelated kicks from fluctuations at scales k ⊥ ρ i ∼ 1. This leads to a diffusion perpendicular to B in velocity space, generally heating the plasma [10, 34, 35]. The theory explains ion heating well in lowβ balanced turbulence simulations [19, 21, 22, 36], and has also been applied to ion heating within the imbalanced near-Sun environment [37-41]. Quasi-linear theory pre-  \nthat ions and waves interact strongly if the waves' Doppler-shifted frequencies are resonant with the ion gyrofrequency [42-44]. This causes ions to diffuse ions in velocity space along contours of constant energy in the frame moving at the wave's phase velocity, heating the plasma (however, note that ions may also be heated by non-resonant interactions with waves [35, 45]).  \nWhile traditionally considered as separate mechanisms, we argue that SH and QLH should instead be considered as two limits of a continuum controlled by the nonlinear broadening of the fluctuations' frequency spectrum. Unbroadened fluctuations-waves with a single frequency ω at a particular k -cause QLH-like behavior; fluctuations that are broadened to a level comparable to ω itself, which thus have significant power in low-frequency fluctuations [46], cause SH-like behavior. This reveals a deep connection to the turbulence imbalance: in imbalanced turbulence z + -fluctuations have linear frequencies that exceed their nonlinear rate (the QLH-limit); in balanced turbulence the two are comparable [29] (the SH-limit). A smooth transition between these limits occurs as the imbalance is adjusted.  \nTo further illustrate the connection, we argue that SH and QLH rates scale similarly with turbulent amplitude, even though SH is traditionally discussed in the context of low-frequency turbulence, while QLH is generally considered to involve high-frequency fluctuations ( ω ∼ Ω i ). SH requires turbulent amplitudes δu ρ i comparable to v th,i (generally δu ρ i /v th,i ≳ 0 . 1 [10, 19]). To understand the amplitude at which QLH will be activated, we use the principle of 'propagation critical balance' (in which parallel scales are determined by the wandering of magnetic field lines by perpendicular perturbations [46-48]), giving ω ∼ k ∥ v A ∼ k ⊥ δu λ throughout the inertial range. Thus, in order to reach ω ∼ Ω i and activate QLH requires Ω i ∼ δu ρ i /ρ i = ⇒ δu ρ i ∼ v th,i , where we used the fact that the ion-eddy interaction is strongest at ρ i scales because ω increases with k ⊥ but interactions with k ⊥ ρ i > 1 fluctuations are suppressed [49] (there may still be contributions from such eddies, however [14, 21]). This means QLH requires similar turbulent amplitudes as SH to occur; likewise, for critically balanced turbulence, the amplitudes needed to cause SH imply turbulent frequencies that approach Ω i , although due to eddies decorrelating, there remains significant power in ω < Ω i fluctuations in balanced turbulence [46]. The argument also implies that heating is strongest at the smallest parallel scale l ∥ ∼ 1 /k ∥ ∼ v A / Ω i with λ ≲ ρ i .  \nBased on the argument above, we propose to extend the usual SH formula for the ion heating rate per unit mass [10],  \nQ ⊥ = Ω i v 2 th,i ˆ c 1 ˆ ξ 3 i e -ˆ c 2 / ˆ ξ i , (1)  \nto also describe QLH in imbalanced turbulence. Here, ˆ c 1 and ˆ c 2 are empirical constants and the exponential  \nsuppression factor accounts phenomenologically for the conservation of magnetic moment at small amplitudes and/or the exponentially small number of particles resonant with lower-frequency fluctuations (in reality, this suppression factor presumably depends on imbalance and other parameters). A key change, however, is needed to correctly describe the helicity barrier, with the definition of ξ i extended from its usual one ( ξ i = δu ρ i /v th,i ) to instead capture the scale at which the fluctuations reach their highest frequency. Thus ˆ ξ i ≡ δu λ /v th,i , where δu λ is the turbulent amplitude at perpendicular scales λ ≳ ρ i corresponding to the smallest parallel scale l ∥ , estimated from the maximum of δu λ /λ via propagation critical balance: v A /l ∥ ∼ δu λ /λ ∼ k ⊥ ˜ z + k ⊥ . Here, ˜ z + k ⊥ is the rms value at scales k ⊥ λ ∼ 1 of ˜ z + , the generalized Elsasser eigenmode of the FLR-MHD model (introduced below), which reduces to its MHD counterpart z + at k ⊥ ρ i ≪ 1.  \nFor turbulence without a helicity barrier (such as balanced turbulence or imbalanced RMHD), the energy cascade is able to dissipate at perpendicular scales smaller than ρ p , forming a smooth spectrum with continuously increasing δu λ /λ . Based on the argument above, we expect ions to be heated by fluctuations at scales λ ∼ ρ p , and Eq. (1) reduces to the normal SH formula with ˆ ξ i = δu ρ i /v th,i . However, with a helicity barrier, the resulting steep spectral break causes k ⊥ ˜ z + k ⊥ to peak at an intermediate scale λ > ρ p before dropping off steeply at smaller scales. We thus use λ ∼ max( ρ i , 1 /k nl ⊥ ) where  \nk nl ⊥ ≡ ∫ ρ -1 i 0 d k ⊥ k ⊥ ( k ⊥ ˜ z + k ⊥ ) ∫ ρ -1 i 0 d k ⊥ k ⊥ ˜ z + k ⊥ , (2)  \nthe average value of k ⊥ weighted by k ⊥ ˜ z + k ⊥ over scales k ⊥ ρ i ≤ 1, is used as a proxy for the scale at which the heating is most efficient. This also takes into account that minor ions interact with fluctuations at scales larger than the spectral break, due to their larger gyroradii.  \nTo verify this proposal, we now present numerical simulations of test protons and minor ions interacting with both balanced and imbalanced turbulence with or without a helicity barrier.  \nNumerical Setup. -We use the 'finite Larmor radius MHD' (FLR-MHD) model [25, 28, 50], which can be formally derived from gyrokinetics in the limit β ≪ 1 and k ⊥ d e ≫ 1 (where d e is the electron inertial length). At k ⊥ ρ i ≪ 1, FLR-MHD reduces to the well-known RMHD model [51]; for k ⊥ ρ i ≫ 1 it reduces to 'electron RMHD' [16, 22, 52]. The FLR-MHD equations are advanced in time using a pseudospectral method [25, 53] (see Suppl. Mat.). The simulation domain is a three-dimensional periodic grid of size L ⊥ = L z = 2 π and a resolution of N 2 ⊥ × N z , where the background magnetic field is along the z -direction and ⊥ represents the x - and y -directions. To test our proposed heating model across a wide range of regimes, we simulate three cases: balanced and imbalanced FLR-MHD with ρ p,gk /L ⊥ = 0 . 02, and imbalanced  \nRMHD (we assume a majority hydrogen plasma, with ρ p,gk the ρ i in the gyrokinetic model). These simulations have a resolution of N 2 ⊥ × N z = 1024 2 × 1024 for the balanced FLR-MHD case and 1024 2 × 512 for the other two cases, refined from an initial resolution of 256 2 × 512 to obtain spectra that extend to scales smaller than ρ p,gk (see Suppl. Mat.). In the imbalanced FLR-MHD simulations we run ions at three different stages of the helicity barrier's evolution as it grows, at times tv A /L z ≈ 3 , 6 , and 10.  \nOnce the refined turbulence simulations have reached a quasi-steady-state, we introduce test ions (for details on implementation, see Suppl. Mat.). For a given collection of ions, we can choose the initial ρ i /L ⊥ , β i , and δu ρ i /v th,i [19]. For FLR-MHD, the choice of ρ i /L ⊥ is fixed to match ρ p,gk /L ⊥ = 0 . 02; we also choose this value for ions in RMHD simulations for comparison. This scale determines δu ρ i via δu 2 ρ i = 2 ∫ k + k -d k ⊥ E E ( k ⊥ ), where E E ( k ⊥ ) is the perpendicular electric-field spectrum (equivalent to the E × B velocity spectrum in FLR-MHD) and using k ± ≡ e ± 1 / 2 /ρ i in line with previous work [10, 19]. This value of δu ρ i with the choice of δu ρ i /v th,i determines v th,i . Due to the invariance of the quantity v A ∇ ∥ in both RMHD and FLR-MHD, we are free to choose v A based on the values of v th,i and β i so long as length-scales along the background field are scaled by the same factor [46]; this rescaling is done within the ion integrator. Once these quantities are determined, the ions are uniformly distributed throughout the simulation domain and given an isotropic Maxwellian velocity distribution. For a given choice of ρ i /L ⊥ and β i , we initialize separate cohorts of N = 10 6 particles with 0 . 01 ≲ δu ρ i /v th,i ≲ 0 . 3.  \nWe introduce minor ions with mass m i = Am p and charge q i = Ze by taking v th,i to be the same of that of protons within the same simulation, so that they have similar initial δu ρ i /v th,i . This choice, along with Ω i = ( Z/A )Ω p , determines the required ion gyroradius  \nTABLE I. f Simulation parameters studied in this Letter. F and R correspond to the FLR-MHD and RMHD models. Coefficients c 1 , c 2 and ˆ c 1 , ˆ c 2 are fits of simulation results to Eq. (1) using δu ρ i /v th,i and ˆ ξ i , respectively; the best fit of all data using ˆ ξ i is ˆ c 1 , ˆ c 2 = 1 . 11 , 0 . 10.  \nFIG. 1. Left: Space-time Fourier transforms E tot ( k z , ω ) from balanced (top) and imbalanced (bottom) FLR-MHD turbulence. Each is normalized to the maximum value at each k z , and the dotted lines are the Alfv'en wave dispersion relation ω = ± k z v A . Right: Comparison of electric-field spectra (top) and k ⊥ ˜ z + k ⊥ (bottom) as used to estimate l ∥ from the different turbulence simulations in this Letter.  \n<!-- image -->  \n⊥  \nρ i = ( A/Z ) ρ p corresponding to ρ p /L ⊥ . Once this is determined, the process of initializing minor ions is identical to protons. Table I summarizes the turbulence and particle parameters studied in this Letter.  \nWe note that because FLR-MHD is derived from gyrokinetics it can only simulate low-frequency dynamics, and does not contain the transition to ion-cyclotron waves. Although this still allows for QLH of protons by Alfv'enic fluctuations, it does not capture their dispersion or polarization changes at k ∥ d i ∼ 1 (where d i is the ion inertial length; see Ref. [26] for the effects of ion-cyclotron waves). This approximation is less impactful for minor ions, because their smaller gyrofrequency means they interact with lower-frequency fluctuations.  \nResults. -To highlight how imbalance affects the frequency spectrum of fluctuations, which was argued above to be important for controlling the heating, the left-hand plot of Fig. 1 shows the space-time Fourier transform E tot ( k z , ω ) = E + ( k z , ω ) + E -( k z , ω ), with E ± ( k z , ω ) = 1 2 〈 | ˜ z ± ( k z , k ⊥ , ω ) | 2 〉 ⊥ ; here, ⟨·⟩ ⊥ denotes an average over all k ⊥ . In balanced turbulence, E tot ( k z , ω ) has significant power at small ω , a signature that the nonlinear and linear times are comparable [46]. In contrast, imbalanced turbulence is dominated by fluctuations at a particular frequency, as seen by the band peaked around ω ≈ -k z v A . These show that increasing imbalance reduces the nonlinear broadening of fluctuations relative to their linear frequency, causing the ion heating to transition from a SH-like to QLH-like mechanism.  \nFIG. 2. Evolution of the ensemble-averaged perpendicular energy per unit mass ⟨ w 2 ⊥ ⟩ of all particle cases with 0 . 02 ≤ δu ρ p /v th,p ≤ 0 . 2 (dark to light) from F-6 β p 0 . 05, along with the linear fit for Q ⊥ in dashed lines.  \n<!-- image -->  \nThe electric-field spectra E E ( k ⊥ ) and k ⊥ ˜ z + k ⊥ for different turbulence simulations are compared on the right of Fig. 1, where we compute k ⊥ ˜ z + k ⊥ from its spectrum: k ⊥ ˜ z + k ⊥ ≈ k ⊥ √ k ⊥ E + ( k ⊥ ) [54]. In contrast to the wide inertial range of the imbalanced RMHD simulation, FLRMHD captures the transition from an Alfv'en-wave to kinetic-Alfv'en-wave cascade [16], with the electric-field specturm approaching a -1 / 3 scaling at k ⊥ ρ i ≳ 1. The presence of a helicity barrier is clearly seen in the imbalanced FLR-MHD simulations via the steep drop in spectra at k ⊥ ρ p ≲ 1, as observed in the solar wind [31, 55] and previous simulations [25-27]. The steep spectral slope of the helicity barrier also causes k ⊥ ˜ z + k ⊥ to reach a maximum at k ⊥ ρ p ≲ 1 scales. This maximum, and the corresponding spectral break in E E ( k ⊥ ), shifts to larger scales as the energy grows in time.  \nFigure 2 illustrates how the ion heating rate Q ⊥ is calculated from the evolution of the ensemble-averaged perpendicular energy per unit mass ⟨ w 2 ⊥ ⟩ for protons from F-6 β p 0 . 05 with 0 . 02 ≤ δu ρ p /v th,p ≤ 0 . 2. Here w ⊥ is the component of w ≡ v -u E × B , the thermal velocity of the protons after subtracting the local E × B velocity, perpendicular to B . Protons experience an initial heating transient lasting around one gyroperiod as they 'pick up' the local E × B velocity (also seen in Ref. [19]). Because of this, we calculate Q ⊥ by a linear fit to ⟨ w 2 ⊥ ⟩ via Q ⊥ = 0 . 5( ⟨ w 2 ⊥ ,f ⟩ - ⟨ w 2 ⊥ ,0 ⟩ ) / ( t f -t 0 ). Here t f is either the end of the simulation or when ⟨ w 2 ⊥ ,f ⟩ = 1 . 2 ⟨ w 2 ⊥ ,0 ⟩ ; we choose this as Q ⊥ decreases as ⟨ w 2 ⊥ ⟩ increases [10, 19]. For δu ρ p /v th,p ≲ 0 . 05, the initial heating transient exhibits oscillations, so to remove the effect of these from the measurement, we choose the start time t 0 Ω i / 2 π = 5 . 25 for δu ρ i /v th,i < 0 . 02, 3 . 25 for 0 . 02 ≤ δu ρ i /v th,i ≤ 0 . 05, and 0 . 75 for δu ρ i /v th,i > 0 . 05 (see Suppl. Mat.).  \nFigure 3 compares the dependence of Q ⊥ on δu ρ i /v th,i , the measure commonly used in SH theory, to our proposed definition ˆ ξ i = δu λ /v th,i . We calculate k nl ⊥ in Eq. (2) using k ⊥ ˜ z + k ⊥ from Fig. 1. It is clear that the  \nstandard normalization is inadequate for describing Q ⊥ , particularly for cases with a helicity barrier. These show greater heating than expected from the power in fluctuations at these scales, indicating the true ˆ ξ i is being underestimated. Minor ions show similar heating rates in both imbalanced FLR-MHD and RMHD, due to their gyroradii lying at scales larger than ρ p and close to or above 1 /k nl ⊥ . In contrast, when plotted against our proposed ˆ ξ i , the data are well described by Eq. (1) with best fit parameters ˆ c 1 = 1 . 11 , ˆ c 2 = 0 . 10 (fits of individual simulations to Eq. (1) comparing δu ρ i /v th,i and ˆ ξ i are listed in Table I); our value of ˆ c 2 is smaller than previous test particle simulations [10, 19] but similar to recent hybridkinetic simulations [22]. That Eq. (1) holds despite the difference in simulation conditions, from turbulence that is balanced, imbalanced, with or without a helicity barrier, and for both protons and minor ions, helps verify our unified model of ion heating. We note that the worst argreement is seen for imbalanced RMHD, which is also unphysical, neglecting the effect of ρ p on the turbulence.  \nDespite the success of the unified model in reproducing Q ⊥ , there remains important differences between balanced- and imbalanced-turbulent heating, as demonstrated in Fig. 4 which compares the proton velocity distribution functions (VDFs) f ( w ⊥ , w ∥ ) in these two regimes. In balanced turbulence, the VDF undergoes diffusion in perpendicular velocity, as seen in the 1D VDF f ( w ⊥ ) = ∫ d w ∥ f ( w ⊥ , w ∥ ) which has a characteristic flattopped shape [22, 38]. In contrast, in imbalanced turbulence the VDF diffuses along constant-energy contours within the frame of the waves, a hallmark of QLH [42]. For Alfv'enic fluctuations, these contours are semicircles centered on the phase speed v ph = ω/k ∥ = -v A . The 1D VDF exhibits a flattopped shape similar to the balanced case (dropping off at slightly earlier upon close inspection), presumably because the contours are nearly vertical at small w ⊥ . This picture also allows an understanding of the β p -dependent dropoff in Q ⊥ seen at large ˆ ξ p in Fig. 3 (representing the largest deviation from the phenomenological model): the constant-energy contours become increasingly curved as β p increases, where v A approaches v th,p , causing ions to gain parallel energy at the expense of perpendicular energy (see the Suppl. Mat.).  \nConclusion. -In this Letter we provide a unified picture for describing ion heating in collisionless plasma turbulence across a wide range of regimes. The heating mechanism changes character depending on the frequency spectrum of turbulent fluctuations, transitioning from SH to QLH as the turbulence imbalance increases and reduces the relative nonlinear broadening of fluctuations. We propose a phenomenological model, controlled by the amplitude of fluctuations at the scale where k ⊥ ˜ z + k ⊥ ∼ v A /l ∥ peaks for k ⊥ ρ i ≲ 1, which accurately captures the ion-heating rate across both regimes. In the absence of a transition-range break or for minor ions,  \n‖  \nFIG. 3. Comparison of the dependence of the measured ion heating rate Q ⊥ against the commonly used stochastic heating parameter δu ρ i /v th,i (left) and the new scale-dependent version ˆ ξ i (right). Filled and hollow markers represent simulations with and without a helicity barrier, respectively. The inset shows the ratio Q ⊥ / (ˆ c 1 ˆ ξ 3 i ), with ˆ c 1 = 1 . 11 the best fit to Eq. (1) for all data, highlighting the suppression of heating for small ξ i . The solid black lines show the best fit of Eq. (1) to F-Balβ p 0 . 05 (left) and all data (right and inset).  \n<!-- image -->  \nFIG. 4. Comparison of two-dimensional velocity distributions f ( w ⊥ , w ∥ ) of ˆ ξ p ≈ 0 . 3 protons from F-Balβ p 0 . 05 (top) and ˆ ξ p ≈ 0 . 7 protons from F-6 β p 0 . 1 (bottom). Dashed contours represent the initial Maxwellian distribution, and solid contours represent constant-energy contours in the frame of the wave, centered on the phase speed v ph = -v A (dotted line). Insets show the 1D perpendicular distribution function f ( w ⊥ ) in blue ( f ( w ⊥ ) at t = 0 is in gray), normalized to the maximum of the initial distribution.  \n<!-- image -->  \nthis scale is usually at k ⊥ ρ i ∼ 1 and our model reduces to the standard SH formula; however, with a transition range, caused by the helicity barrier, this scale can lie above at k ⊥ ρ i < 1. Using high-resolution simulations of test particles interacting with turbulence, we show that our proposal works well in describing measured heating  \nrates for both protons and minor ions in turbulence that is balanced, imbalanced, with or without a helicity barrier. The fact that this holds, despite the very different 'arced' VDFs that result in imbalanced turbulence, hints at the universality of the heating mechanism.  \nWe note that the FLR-MHD model, on which our results are based, is highly simplified, missing out the transition to ion-cyclotron waves when ω ∼ Ω i . Additionally, we also only measure the instantaneous ion Q ⊥ from an initial Maxwellian, rather than tracking large changes of the velocity distribution function (as occurs for strong heating). Both of these issues have been addressed in a recent hybrid-kinetic study [56], which shows similar features. Nonetheless, the results of this Letter help to disentangle the signatures of ion heating and provide a base for more complex theories. The ideas presented here can also be used to understand helicity-barrier saturation on ion heating [26, 27], which is not captured by FLR-MHD and may be important for ab-initio models of coronal and solar-wind heating.  \nThe authors would like to thank T. Adkins, B. D. G. Chandran, M. W. Kunz, N. F. Loureiro, M. Zhang, and M. Zhou for helpful discussions over the course of this work. Research is supported by the University of Otago, through a University of Otago Doctoral Scholarship (ZJ), and the Royal Society Te Ap¯arangi, through Marsden-Fund grant MFP-UOO2221 (JS) and MFP-U0020 (RM), as well as through the Rutherford Discovery Fellowship RDF-U001004 (JS). High-performance computing resources were provided by the New Zealand eScience Infrastructure (NeSI) under project grant uoo02637.  \n- ∗ johza721@student.otago.ac.nz\n- [1] E. Marsch, K.-H. Muhlhauser, R. Schwenn, H. Rosenbauer, W. Pilipp, and F. M. Neubauer, Solar wind protons: Three-dimensional velocity distributions and derived plasma parameters measured between 0.3 and 1 AU, J. Geophys. Res. 87 , 52 (1982).\n- [2] E. Marsch, On the temperature anisotropy of the core part of the proton velocity distribution function in the solar wind, J. Geophys. Res. 109 , 10.1029/2003ja010330 (2004).\n- [3] P. Hellinger, P. Tr'avn'ıˇcek, J. C. Kasper, and A. J. Lazarus, Solar wind proton temperature anisotropy: Linear theory and wind/swe observations, Geophys. Res. Lett. 33 , 10.1029/2006gl025925 (2006).\n- [4] J. L. Kohl, G. Noci, E. Antonucci, G. Tondello, M. C. E. Huber, S. R. Cranmer, L. Strachan, A. V. Panasyuk, L. D. Gardner, M. Romoli, S. Fineschi, D. Dobrzycka, J. C. Raymond, P. Nicolosi, O. H. W. Siegmund, D. Spadaro, C. Benna, A. Ciaravella, S. Giordano, S. R. Habbal, M. Karovska, X. Li, R. Martin, J. G. Michels, A. Modigliani, G. Naletto, R. H. O'Neal, C. Pernechele, G. Poletto, P. L. Smith, and R. M. Suleiman, UVCS/ SOHO empirical determinations of anisotropic velocity distributions in the solar corona, Astrophys. J. 501 , L127 (1998).\n- [5] R. Esser, S. Fineschi, D. Dobrzycka, S. R. Habbal, R. J. Edgar, J. C. Raymond, J. L. Kohl, and M. Guhathakurta, Plasma properties in coronal holes derived from measurements of minor ion spectral lines and polarized white light intensity, Astrophys. J. 510 , L63 (1999).\n- [6] E. Antonucci, Fast solar wind velocity in a polar coronal hole during solar minimum, Solar Phys. 197 , 115 (2000).\n- [7] S. R. Cranmer, W. H. Matthaeus, B. A. Breech, and J. C. Kasper, Empirical constraints on proton and electron heating in the fast solar wind, Astrophys. J. 702 , 1604 (2009).\n- [8] G. Howes, A prescription for the turbulent heating of astrophysical plasmas, Monthly Notices of the Royal Astronomical Society: Letters 409 , L104 (2010), 1009.4212.\n- [9] Y. Kawazura, M. Barnes, and A. A. Schekochihin, Thermal disequilibration of ions and electrons by collisionless plasma turbulence, Proc. Natl. Acad. Sci. U. S. A. 116 , 771 (2019).\n- [10] B. D. G. Chandran, B. Li, B. N. Rogers, E. Quataert, and K. Germaschewski, Perpendicular ion heating by low-frequency Alfv'en-wave turbulence in the solar wind, Astrophys. J. 720 , 503 (2010).\n- [11] B. D. G. Chandran, D. Verscharen, E. Quataert, J. C. Kasper, P. A. Isenberg, and S. Bourouaine, Stochastic heating, differential flow, and the alpha-to-proton temperature ratio in the solar wind, ApJ 776 , 45 (2013).\n- [12] J. V. Hollweg and P. A. Isenberg, Generation of the fast solar wind: A review with emphasis on the resonant cyclotron interaction, Journal of Geophysical Research: Space Physics 107 , SSH 12 (2002).\n- [13] P. A. Isenberg and B. J. Vasquez, A kinetic model of solar wind generation by oblique ion-cyclotron waves, ApJ 731 , 88 (2011).\n- [14] P. A. Isenberg and B. J. Vasquez, Perpendicular ion heating by cyclotron resonant dissipation of turbulently generated kinetic Alfv'en waves in the solar wind, Astrophys.  \n- J. 887 , 63 (2019).\n- [15] B. D. G. Chandran, Weak compressible magnetohydrodynamic turbulence in the solar corona, Phys. Rev. Lett. 95 , 265004 (2005).\n- [16] A. A. Schekochihin, S. C. Cowley, W. Dorland, G. W. Hammett, G. G. Howes, E. Quataert, and T. Tatsuno, Astrophysical gyrokinetics: kinetic and fluid turbulent cascades in magnetized weakly collisional plasmas, Astrophys. J. Supp. 182 , 310 (2009).\n- [17] J. M. TenBarge and G. G. Howes, Current sheets and collisionless damping in kinetic plasma turbulence, ApJL 771 , L27 (2013).\n- [18] R. Lehe, I. J. Parrish, and E. Quataert, The heating of test particles in numerical simulations of Alfv'enic turbulence, Astrophys. J. 707 , 404 (2009).\n- [19] Q. Xia, J. C. Perez, B. D. G. Chandran, and E. Quataert, Perpendicular ion heating by reduced magnetohydrodynamic turbulence, Astrophys. J. 776 , 90 (2013).\n- [20] B. Teaca, M. S. Weidl, F. Jenko, and R. Schlickeiser, Acceleration of particles in imbalanced magnetohydrodynamic turbulence, Phys. Rev. E Stat. Nonlin. Soft Matter Phys. 90 , 021101 (2014).\n- [21] L. Arzamasskiy, M. W. Kunz, B. D. G. Chandran, and E. Quataert, Hybrid-kinetic simulations of ion heating in Alfv'enic turbulence, Astrophys. J. 879 , 53 (2019).\n- [22] S. S. Cerri, L. Arzamasskiy, and M. W. Kunz, On stochastic heating and its phase-space signatures in lowbeta kinetic turbulence, Astrophys. J. 916 , 120 (2021).\n- [23] M. S. Weidl, F. Jenko, B. Teaca, and R. Schlickeiser, Cosmic-ray pitch-angle scattering in imbalanced MHD turbulence simulations, ApJ 811 , 8 (2015).\n- [24] G. G. Howes, The fundamental parameters of astrophysical plasma turbulence and its dissipation: Nonrelativistic limit, arXiv [astro-ph.SR] (2024), arXiv:2402.12829 [astro-ph.SR].\n- [25] R. Meyrand, J. Squire, A. A. Schekochihin, and W. Dorland, On the violation of the zeroth law of turbulence in space plasmas, J. Plasma Phys. 87 , 10.1017/S0022377821000489 (2021).\n- [26] J. Squire, R. Meyrand, M. W. Kunz, L. Arzamasskiy, A. A. Schekochihin, and E. Quataert, High-frequency heating of the solar wind triggered by low-frequency turbulence, Nature Astronomy 6 , 715 (2022).\n- [27] J. Squire, R. Meyrand, and M. W. Kunz, Electron-ion heating partition in imbalanced solar-wind turbulence, Astrophys. J. Lett. 957 , L30 (2023).\n- [28] A. A. Schekochihin, Y. Kawazura, and M. A. Barnes, Constraints on ion versus electron heating by plasma turbulence at low beta, J. Plasma Phys. 85 , 10.1017/S0022377819000345 (2019).\n- [29] P. Goldreich and S. Sridhar, Toward a theory of interstellar turbulence. 2: Strong Alfv'enic turbulence, Astrophys. J. 438 , 763 (1995).\n- [30] R. J. Leamon, C. W. Smith, N. F. Ness, W. H. Matthaeus, and H. K. Wong, Observational constraints on the dynamics of the interplanetary magnetic field dissipation range, Journal of Geophysical Research: Space Physics 103 , 4775 (1998).\n- [31] T. A. Bowen, A. Mallet, S. D. Bale, J. W. Bonnell, A. W. Case, B. D. G. Chandran, A. Chasapis, C. H. K. Chen, D. Duan, T. Dudok de Wit, K. Goetz, J. S. Halekas, P. R. Harvey, J. C. Kasper, K. E. Korreck, D. Larson, R. Livi, R. J. MacDowall, D. M. Malaspina, M. D. McManus, M. Pulupa, M. Stevens, and P. Whittlesey, Constraining\n- ion-scale heating and spectral energy transfer in observations of plasma turbulence, Phys. Rev. Lett. 125 , 025102 (2020).\n- [32] T. A. Bowen, B. D. G. Chandran, J. Squire, S. D. Bale, D. Duan, K. G. Klein, D. Larson, A. Mallet, M. D. McManus, R. Meyrand, J. L. Verniero, and L. D. Woodham, In situ signature of cyclotron resonant heating in the solar wind, Phys. Rev. Lett. 129 , 165101 (2022).\n- [33] T. A. Bowen, S. D. Bale, B. D. G. Chandran, A. Chasapis, C. H. K. Chen, T. Dudok de Wit, A. Mallet, R. Meyrand, and J. Squire, Mediation of collisionless turbulent dissipation through cyclotron resonance, Nat. Astron. , 1 (2024).\n- [34] J. M. McChesney, R. A. Stern, and P. M. Bellan, Observation of fast stochastic ion heating by drift waves, Phys. Rev. Lett. 59 , 1436 (1987).\n- [35] J. R. Johnson and C. Z. Cheng, Stochastic ion heating at the magnetopause due to kinetic Alfv'en waves, Geophys. Res. Lett. 28 , 4421 (2001).\n- [36] I. W. Hoppock, B. D. G. Chandran, K. G. Klein, A. Mallet, and D. Verscharen, Stochastic proton heating by kinetic-Alfv'en-wave turbulence in moderately highβ plasmas, J. Plasma Phys. 84 , 10.1017/S0022377818001277 (2018).\n- [37] S. Bourouaine and B. D. G. Chandran, Observational test of stochastic heating in lowβ fast-solar-wind streams, ApJ 774 , 96 (2013).\n- [38] K. G. Klein and B. D. G. Chandran, Evolution of the proton velocity distribution due to stochastic heating in the near-sun solar wind, Astrophys. J. 820 , 47 (2016).\n- [39] D. Vech, K. G. Klein, and J. C. Kasper, Nature of stochastic ion heating in the solar wind: Testing the dependence on plasma beta and turbulence amplitude, Astrophys. J. Lett. 850 , L11 (2017).\n- [40] M. M. Martinovi'c, K. G. Klein, and S. Bourouaine, Radial evolution of stochastic heating in lowβ solar wind, ApJ 879 , 43 (2019).\n- [41] M. M. Martinovi'c, K. G. Klein, J. C. Kasper, A. W. Case, K. E. Korreck, D. Larson, R. Livi, M. Stevens, P. Whittlesey, B. D. G. Chandran, B. L. Alterman, J. Huang, C. H. K. Chen, S. D. Bale, M. Pulupa, D. M. Malaspina, J. W. Bonnell, P. R. Harvey, K. Goetz, T. Dudok de Wit, and R. J. MacDowall, The enhancement of proton stochastic heating in the near-Sun solar wind, Astrophys. J. Suppl. Ser. 246 , 30 (2020).\n- [42] C. F. Kennel and F. Engelmann, Velocity space diffusion from weak plasma turbulence in a magnetic field, Phys. Fluids 9 , 2377 (1966).\n- [43] T. H. Stix, Waves in Plasmas , 1992nd ed. (American Institute of Physics, New York, NY, 1992).\n- [44] R. Schlickeiser and U. Achatz, Cosmic-ray particle transport in weakly turbulent plasmas. part 1. theory, J. Plasma Phys. 49 , 63 (1993).\n- [45] L. Chen, Z. Lin, and R. White, On resonant heating below the cyclotron frequency, Phys. Plasmas 8 , 4713 (2001).\n- [46] A. A. Schekochihin, MHD turbulence: a biased review, J. Plasma Phys. 88 , 155880501 (2022).\n- [47] Y. Lithwick, P. Goldreich, and S. Sridhar, Imbalanced strong MHD turbulence, Astrophys. J. 655 , 269 (2007).\n- [48] A. Beresnyak and A. Lazarian, Strong imbalanced turbulence, Astrophys. J. 682 , 1070 (2008).\n- [49] B. D. Chandran, Scattering of energetic particles by anisotropic magnetohydrodynamic turbulence with a  \n- Goldreich-Sridhar power spectrum, Phys. Rev. Lett. 85 , 4656 (2000).\n- [50] T. Passot, P. L. Sulem, and E. Tassi, Gyrofluid modeling and phenomenology of lowβ e Alfv'en wave turbulence, Phys. Plasmas 25 , 042107 (2018).\n- [51] H. R. Strauss, Nonlinear, three-dimensional magnetohydrodynamics of noncircular tokamaks, Phys. Fluids 19 , 134 (1976).\n- [52] S. Boldyrev, K. Horaites, Q. Xia, and J. C. Perez, Toward a theory of astrophysical plasma turbulence at subproton scales, Astrophys. J. 777 , 41 (2013).\n- [53] B. Teaca, M. K. Verma, B. Knaepen, and D. Carati, Energy transfer in anisotropic magnetohydrodynamic turbulence, Phys. Rev. E Stat. Nonlin. Soft Matter Phys. 79 , 046312 (2009).\n- [54] A. F. A. Bott, L. Arzamasskiy, M. W. Kunz, E. Quataert, and J. Squire, Adaptive critical balance and firehose instability in an expanding, turbulent, collisionless plasma, Astrophys. J. Lett. 922 , L35 (2021).\n- [55] K. H. Kiyani, K. T. Osman, and S. C. Chapman, Dissipation and heating in solar wind turbulence: from the macro to the micro and back again, Philos. Trans. A Math. Phys. Eng. Sci. 373 , 10.1098/rsta.2014.0155 (2015).\n- [56] M. F. Zhang, M. W. Kunz, J. Squire, and K. G. Klein, Extreme heating of minor ions in imbalanced solar-wind turbulence, arXiv [astro-ph.SR] (2024), arXiv:2408.04703 [astro-ph.SR].\n- [57] G. G. Howes, S. C. Cowley, W. Dorland, G. W. Hammett, E. Quataert, and A. A. Schekochihin, Astrophysical gyrokinetics: Basic equations and linear theory, Astrophys. J. 651 , 590 (2006).\n- [58] A. Zocco and A. A. Schekochihin, Reduced fluid-kinetic equations for low-frequency dynamics, magnetic reconnection, and electron heating in low-beta plasmas, Phys. Plasmas 18 , 102309 (2011).\n- [59] R. Meyrand, A. Kanekar, W. Dorland, and A. A. Schekochihin, Fluidization of collisionless plasma turbulence, Proc. Natl. Acad. Sci. U. S. A. 116 , 1185 (2019).\n- [60] V. David, S. Galtier, and R. Meyrand, Monofractality in the solar wind at electron scales: Insights from kinetic Alfv'en waves turbulence, Phys. Rev. Lett. 132 , 085201 (2024).\n- [61] J. H. Williamson, Low-storage Runge-Kutta schemes, J. Comput. Phys. 35 , 48 (1980).\n- [62] T. Adkins, R. Meyrand, and J. Squire, The effects of finite electron inertia on helicity-barrier-mediated turbulence, arXiv [physics.plasm-ph] (2024), arXiv:2404.09380 [physics.plasm-ph].\n- [63] A. Mallet, A. A. Schekochihin, and B. D. G. Chandran, Disruption of Alfv'enic turbulence by magnetic reconnection in a collisionless plasma, J. Plasma Phys. 83 , 905830609 (2017).\n- [64] N. F. Loureiro and S. Boldyrev, Collisionless reconnection in magnetohydrodynamic and kinetic turbulence, Astrophys. J. 850 , 182 (2017).\n- [65] M. Zhou, Z. Liu, and N. F. Loureiro, Electron heating in kinetic-Alfv'en-wave turbulence, Proc. Natl. Acad. Sci. U. S. A. 120 , e2220927120 (2023).\n- [66] J. P. Boris, Relativistic plasma simulation-optimization, in Proc. Fourth Conf. Num. Sim. Plasmas (1970) pp. 367.\n- [67] C. K. Birdsall and A. B. Langdon, Plasma physics via computer simulation (CRC Press, 2018).", 'The FLR-MHD Model': "We use the finite Larmor radius MHD (FLR-MHD) model [25, 28, 50], a reduced fluid model that includes effects at scales both larger and smaller than ρ i with-  \nut needing to solve kinetic equations. This theory can be derived from the Vlasov-Maxwell equations by assuming all quantities vary on timescales slower than the ioncyclotron frequency, a strong background magnetic field, and that perturbations are elongated along magnetic field lines [16, 57]. Assuming a background magnetic field B 0 = B 0 ˆ z and taking the lowβ limit [28, 58] gives the equations of FLR-MHD:  \n( ∂ ∂t + u ⊥ · ∇ ⊥ ) δn e n 0 e = -c 4 πen 0 e ( ∂ ∂z + b ⊥ · ∇ ⊥ ) ∇ 2 ⊥ A ∥ + D 6 ν δn e n 0 e , (3)  \n( ∂ ∂t + u ⊥ · ∇ ⊥ ) A ∥ = -c ∂φ ∂z + cT 0 e e ( ∂ ∂z + b ⊥ · ∇ ⊥ ) δn e n 0 e + D 6 η A ∥ , (4)  \nδn e n 0 e = -Z τ ( 1 -ˆ Γ 0 ) eφ T 0 e , (5)  \nwhere τ = T 0i /T 0e is the ratio of the ion and electron temperatures, δn e /n 0 e = δn i /n 0i is the perturbed electron density (equal to the ion density by quasi-neutrality), A ∥ is the ˆ z component of the magnetic vector potential, φ is the electrostatic potential, u ⊥ = ( c/B 0 )ˆ z × ∇ ⊥ φ is the perpendicular E × B flow, and b ⊥ = -ˆ z × ∇ ⊥ A ∥ /B 0 is the magnetic field perturbation. The gyrokinetic Poisson operator 1 -ˆ Γ 0 = 1 -I 0 ( α ) e -α , where I 0 is the modified Bessel function of the first kind and α = -ρ 2 i ∇ 2 ⊥ / 2, models the transition from an Alfv'enic to a kinetic-Alfv'enwave cascade that occurs at k ⊥ ρ i ∼ 1. In our simulations a hyperdiffusion operator D 6 ν = ν 6 ⊥ ∇ 6 ⊥ + ν 6 z ∇ 6 z (and similarly for D 6 η ) is used to dissipate energy above grid scales, but does not model any specific physical process; we set η ⊥ ,z = ν ⊥ ,z in all cases except balanced FLRMHD, discussed below.  \nThis model describes the evolution of the Alfv'enic component of turbulence both above and below ρ i ; the assumption of small β leads to minimal coupling between this and the compressive component of the turbulence [16]. The FLR-MHD model is valid as long as the electron inertial length d e ≪ ρ i . This is because d e sets the length scale below which the magnetic field is no longer frozen into the electron flow; in this regime the equations are coupled to the electron distribution function and are no longer closed. Further details on the properties of FLR-MHD are given in Ref. [25].", 'Eigenmodes and Conserved Invarients': "For a given wavenumber k , the FLR-MHD model supports eigenmodes that travel parallel and anti-parallel to  \nB 0 given by  \nΘ ± k = -Ω i v ph ( k ⊥ ) k 2 ⊥ δn e n 0 e ∓ A ∥ √ 4 πm i n 0i , (6)  \nwith frequency ω = ± k z v ph ( k ⊥ ) v A . The normalized phase speed v ph ( k ⊥ ) = k ⊥ ρ i [(1 / [1 -ˆ Γ 0 ] + Z/τ ) / 2] 1 / 2 reduces to 1 at k ⊥ ρ i ≪ 1 scales and is proportional to k ⊥ ρ i at k ⊥ ρ i ≫ 1 scales, reflecting the transition from Alfv'en to kinetic-Alfv'en waves. These eigenmodes allow the definition of generalized Elsasser variables ˜ z ± = ˆ z × ∇ ⊥ Θ ± that reduce to the standard Elsasser variables of MHD ( z ± = u ⊥ ± B ⊥ / √ 4 πm i n 0i ) in the limit k ⊥ ρ i ≪ 1.  \nThis model also has two conserved invariants: the free energy E = 1 4 ∑ k ( | k ⊥ Θ + | 2 + | k ⊥ Θ -| 2 ), and the generalized helicity H = 1 4 ∑ k ( | k ⊥ Θ + | 2 -| k ⊥ Θ -| 2 ) /v ph ( k ⊥ ). The helicity barrier results from the inability of the system to support a forward cascade of both E and H simultaneously at k ⊥ ρ i ≫ 1, due to the scaling v ph ∝ k ⊥ in this range [25].", 'Additional Simulation Details': 'We use asterix [25, 59, 60], a modified verison of the pseudospectral code turbo [53], to solve the equations of RMHD and FLR-MHD, which are advanced in time with a third-order modified Williamson algorithm (a four-step, low-order Runge-Kutta method [61]). We force the imbalanced FLR-MHD and RMHD simulations with energy injection ε and σ ε = 0 . 88, where σ ε = | ε H | /ε is the ratio of the helicity and energy injection rates. In contrast, the balanced FLR-MHD simulation has no helicity injected into the system ( ε H = 0).  \nFigure 5 shows the magnetic-field spectra of the highresolution simulations after refinement (discussed below),  \nFigure 6 compares the time-evolution of dissipation in all balanced and imbalanced RMHD and FLR-MHD cases during the refinement process. The imbalanced cases take longer to saturate compared to the balanced case for the reason discussed above. As we increase the perpendicular resolution while refining the cascade can reach smaller corresponding parallel scales due to critical balance, leading to an increase in parallel dissipation. This effect is clearly seen in the balanced FLR-MHD case:  \n<!-- image -->  \n⊥  \nFIG. 5. Comparison of magnetic-field spectra from the different turbulence simulations in this Letter.  \ncorresponding to the electric-field spectra in Figure 1 of the Letter. All simulations exhibit an approximate -3 / 2 scaling for k ⊥ ρ i ≲ 1. With the addition of a helicity barrier, the magnetic-field spectra show similar scalings to electric-field spectra, steeping to a -4 scaling around k ⊥ ρ i ∼ 1.', 'Choice of Parallel Viscosity': "̸  \nWe use different choices of the parallel viscosity ν 6 z in the different turbulent regimes to reflect how they arrive at a steady state. For imbalanced RMHD turbulence, the majority of the energy cascade is able to reach small perpendicular scales and dissipate there. However, due to critical balance this energy also reaches small parallel scales. To ensure correct dissipation at small scales, we use a non-zero parallel viscosity. For the helicity barrier, however, we are interested in how it affects the heating of test particles as it grows in time during the 'pseudostationary' phase, before it saturates unphysically on the parallel hyperdissipation. We thus set ν 6 z = 0 for the helicity barrier simulation, studying it well before it would saturate with ν 6 z = 0. Various tests (not shown) found that this produces results similar to using a non-zero ν 6 z .", 'Simulation Refinement': "A highly imbalanced forcing causes the turbulence to take longer to reach a steady-state, as the energy cascade is less efficient for the stronger Elsasser field [46, 47]. The addition of a helicity barrier complicates this further as the system evolves through a 'pseudo-stationary' state where the balanced portion of the cascade dissipates at small perpendicular scales, whereas the imbalanced portion grows at scales k ⊥ ρ i ≲ 1 before saturating unphysically on parallel hyperdissipation [25, 26]. To run a highresolution simulation, such as the ones used in this work, to this saturation point would be time-consuming and  \nFIG. 6. Time evolution of turbulent dissipation for the imbalanced RMHD (top), and the balanced and imbalanced FLR-MHD cases (middle, bottom respectively). Contributions to the parallel (orange) and perpendicular (blue) dissipation and their total (black) are normalized to the energy input ε . The red dotted lines show when the simulations were refined, doubling the perpendicular resolution from 256 2 × 512 up to 1024 2 × 512 for the imbalanced cases or 1024 3 for the balanced case. Note the difference in scale on both axes; see text for discussion.  \n<!-- image -->  \ncomputationally expensive.  \nInstead, we start from simulations with a resolution of N 2 ⊥ × N z = 256 2 × 512 and let these evolve to a steadystate (or pseudo-stationary state for simulations with a helicity barrier). Once this point is reached, we refine the simulation by doubling the perpendicular resolution and decreasing the perpendicular viscosity ν 6 ⊥ , allowing the energy cascade to reach smaller scales before dissipating. Once the cascade has reached a new steady-state, we repeat this refinement procedure until the desired resolution is reached, which is then used as the starting point for our test particle simulations. The final resolutions are N ⊥ = 1024 , N z = 512 for the imbalanced RMHD and FLR-MHD cases, and N ⊥ = N z = 1024 for the balanced FLR-MHD case (the reason for this is discussed below).  \nFigure 7 compares two-dimensional dissipation spectra D 2D ( k ⊥ , k z ) from the high-resolution imbalanced FLRMHD and RMHD simulations. In the absence of a helicity barrier, all energy is dissipated at grid scales in the RMHD case. In contrast, the presence of a helicity barrier causes a flat dissipation spectrum around k ⊥ ρ i ∼ 1, due to the steep energy spectrum slope in this range.  \n<!-- image -->  \n⊥  \n⊥  \n⊥  \n⊥  \nFIG. 7. Two-dimensional dissipation spectra D 2D ( k ⊥ , k z ) from the 1024 2 × 512 imbalanced FLR-MHD (left) and RMHD (right) simulations, normalized to their maximum values. The vertical dotted lines represent the test particle gyroradius scale k ⊥ ρ i = 1.  \nthe initial refinement in N ⊥ leads to an increase in ε z , and the subsequent refinement in both N ⊥ and N z resulting in a decrease (due to dissipation scales occuring at larger k z ). In the imbalanced RMHD and balanced FLR-MHD cases, all energy is dissipated near the grid at predominately small perpendicular scales with ε diss ≈ ε . In contrast, the helicity barrier in the imbalanced FLRMHDcase blocks the majority of the cascade, only letting the balanced portion of the flux ( ∼ 1 -σ ε ) through to be dissipated at grid scales; the remainder is 'trapped' at large perpendicular scales. The lower-resolution simulations allow more of the cascade to be dissipated than predicted as their dissipation scales are close to ρ i -scales; this effect decreases as the range between these scales increases with resolution, as seen by the lower ε diss /ε at each refinement level, converging to ε diss /ε = 1 -σ ε .", 'Viscosity in the Balanced FLR-MHD Case': "A notable exception to the methods described above was the balanced FLR-MHD turbulence case. Due to the transition from an Alfv'enic to a kinetic-Alfv'en-wave cascade at scales k ⊥ ρ i ≳ 1 we expect the slope of the magnetic-field spectrum to steepen from -5 / 3 to -7 / 3, and the electric-field specturm to flatten from -5 / 3 to -1 / 3 [62]. However, when refining the balanced FLRMHD simulation to a resolution of 1024 2 × 512 with the method above we found that both the magnetic- and electric-field spectra were steeper than expected (dashed  \n̸  \nblue and red lines in Fig. 8, respectively), with the electric field reaching a slope of around -7 / 3. These spectra (solid lines in Fig. 8) approach their theoretical predictions after refining the simulation to a resolution of 1024 3 and switching off only the viscosity ( ν ⊥ ,z = 0 and η = 0), leaving the electric field unaffected by dissipation (a method used in hybrid-kinetic simulations [21, 22, 26]). This steepening appears to be an unexpected nonlinear effect previously unseen in the FLR-MHD model. Previous work has shown that the reconnection of current sheets by subρ i effects can lead to a steepening of the energy spectrum similar to what is seen in the transition range [63-65]; however, these require the electron inertial length d e which is ordered out of FLR-MHD. Whether what we observe here is related to this physics remains unclear, and more work is required to clarify the details of this mechanism. Because of this, we integrate particles using the ν = 0 balanced FLR-MHD simulation in the main Letter to remove any influence on heating due to this effect.", 'Additional Test Particle Details': 'In order to study particle heating, we introduced the ability to integrate test particles simultaneously with the fields solved by asterix . As in previous test particle simulations (e.g. [18, 19]), the equations of motion of the particles are integrated using the Boris push [66] with electric and magnetic fields interpolated from the grid to their positions via the triangular-shaped-cloud (TSC) method [18, 67]. The test particles are run simultaneously with the fluid code so that particles can use the current fields without needing to save the full history of the simulation. The perpendicular fields are obtained via B ⊥ = B 0 b ⊥ and E ⊥ = -u ⊥ × B /c (with u ⊥ and B ⊥ as defined above via φ and A ∥ in  \nFIG. 8. Magnetic- (blue) and electric-field (red) spectra and spectral slope α (measured by fitting over a range ± log 10 ( k ⊥ ) ≈ 0 . 1 around each point), highlighting the dependence on dissipation by viscosity. Dashed lines show the spectra from a N 2 ⊥ × N z = 1024 2 × 512 balanced FLR-MHD simulation with viscosity dissipation activated, and solid lines from the 1024 3 balanced simulation with ν = 0 (used in the main Letter).  \n<!-- image -->  \n⊥  \nFLR-MHD); any artificial parallel electric field component due to interpolation is then removed using the method of Ref. [18]. We calculate the required timestep via ∆ t particle = min(2 π/ ( N 1 Ω i ) , ∆ t fluid /N 2 ) [18], where ∆ t fluid is the fluid timestep calculated by asterix . Setting N 1 = 50 and N 2 = 5 was found to give accurate results; these results converge so long as N 1 ≳ 10 (also seen by Ref. [18]).  \nThere are some important differences between RMHD and FLR-MHD that need to be taken into account when implementing test particles. Firstly, the initial ρ i /L ⊥ of the particle distribution is no longer a free parameter due to the inclusion of kinetic scales in the FLR-MHD model (in comparison to RMHD, which assumes ρ i → 0). We fix this parameter during the initialization of particles in an FLR-MHD simulation. Secondly, the inclusion of k ⊥ ρ i ≫ 1 scales introduces a electric field component parallel to B = B 0 + B 0 b ⊥ that arises from the conservation of magnetic flux in gyrokinetics [16, 22], given by  \nE ∥ = E · ˆ b = -( ∂ ∂z + b ⊥ · ∇ )( T 0 e e δn e n 0 e )  \nwhere ˆ b is the unit vector along B . We include this term when interpolating fields to the particles, although it is generally small in comparison to E ⊥ .  \nFIG. 9. Measured heating rate Q ⊥ of protons in 1024 3 balanced FLR-MHD turbulence (top), demonstrating the greater-than-expected heating at low turbulent amplitudes. The solid line shows the best fit to the empirical Q ⊥ (Eq. 1 in Letter) for F-Balβ p 0 . 05 with c 1 , c 2 = 1 . 83 , 0 . 18. These heating rates are measured using a linear fit (dashed lines; bottom) to the ensemble-averaged perpendicular energy per unit mass ⟨ w 2 ⊥ ⟩ for the protons with 5 × 10 -3 ≤ δu ρ p /v th,p ≤ 2 × 10 -2 (dark to light). The oscillations are clearly seen and have mostly decayed after approximately 5 gyroperiods; Q ⊥ is measured from 5 . 25 gyroperiods onward for these smallδu ρ i /v th,i particles. For clarity, ⟨ w 2 ⊥ ⟩ is normalized to its initial values; the values used for the heating rate calculations are normalized to v 2 A .  \n<!-- image -->', 'Testing of Implementation': 'After its implementation, the particle code was extensively tested to verify its correctness. Simple cases with known analytical properties (such as motion in constant E ⊥ and B 0 fields as well as conservation of energy in E = -∇ φ fields) are reproduced with a relative error of less than 10 -3 . The implementation of the TSC interpolation scheme was found to exhibit second-order convergence up to a resolution of 128 3 , before additional errors dominate (e.g., from timestepping). The accuracy of this scheme was tested by comparing the use of gridinterpolated fields to those at the particle position via more complicated analytical fields (such as a linear combination of sinusoidal electromagnetic modes); for a grid resolution of 128 3 or greater, the maximum relative error in energy evolution is less than 10 -2 .  \nThe convergence of statistics with particle number was tested by running cohorts of N particles (with N ranging from 100 to 10 6 in powers of 10) in the same 128 3 balanced RMHD turbulence simulation. The relative error in the statistics of an N = 10 5 and N = 10 6 run is ≲ 1%; we use N = 10 6 particles to allow for better statistics when calculating the velocity distribution functions.  \nFIG. 10. Comparison of two-dimensional velocity distributions f ( w ⊥ , w ∥ ) of protons with ˆ ξ p ≈ 0 . 7 from F-6 β p 0 . 01 (top), F-6 β p 0 . 05 (middle), and F-6 β p 0 . 1 (bottom; same as Fig. 4 of main paper). Solid contours represent constant-energy contours (Eq. 7) in the frame of the wave, centered on the phase speed v ph = -v A (dashed line).  \n<!-- image -->  \n‖  \nParticle Heating at Small Turbulent Amplitudes  \nIons experience a greater-than-expected heating rate than the empirical exponential suppression in Q ⊥ = Ω i v 2 th,i ˆ c 1 ˆ ξ 3 i e -ˆ c 2 / ˆ ξ i (Eq. 1 of Letter) when initialized with very low δu ρ i /v th,i . In Figure 9 we compare the measured heating rates of protons with 5 × 10 -3 ≤ δu ρ p /v th,p ≤ 1 × 10 -2 interacting with the same 1024 3 balanced FLR-MHD turbulence simulation as those with  \n0 . 02 ≤ δu ρ p /v th,p ≤ 0 . 3 (F-Balβ p 0 . 05 from the main Letter). Due to the low-amplitude turbulence at ρ p -scales, the evolution of ⟨ w 2 ⊥ ⟩ for these protons exhibits oscillations from the large-scale E × B flow that introduces some ambiguity in measuring the heating rate; these oscillations are reduced with increasing δu ρ p /v th,p as they experience greater heating (e.g., Figure 2 of the main Letter). We find that the oscillations decay after ∼ 5 gyroperiods for ions initialized with δu ρ i /v th,i < 0 . 02, ∼ 3 for those with 0 . 02 ≤ δu ρ i /v th,i ≤ 0 . 05, and ∼ 0 . 5 for those with δu ρ i /v th,i > 0 . 05. To remove any contributions from the oscillations when measuring the linear growth of ⟨ w 2 ⊥ ⟩ , we choose the start time of the linear fit for Q ⊥ to be 5 . 25, 3 . 25, and 0 . 75 gyroperiods respectively for the above cases (with the quarter-gyroperiod offset ensuring we are in the middle of any further oscillation cycle). Despite this, the measured heating rate is greater than that expected by exponential suppression ( Q ⊥ / Ω p v 2 th,p ≲ 10 -12 for δu ρ p /v th,p ≲ 10 -2 ).', 'Beta-dependent Perpendicular Heating in Imbalanced FLR-MHD Turbulence': "In Fig. 10, we show how the diffusion of the proton velocity distribution function in imbalanced FLR-MHD turbulence depends on β p = v 2 th,p /v 2 A . Quasi-linear theory predicts that the evolution of the distribution function f ( w ⊥ , w ∥ ) is a diffusion along constant-energy contours  \nK ' = m p 2 [ w 2 ⊥ +( w ∥ -v ph ) 2 ] (7)  \nwithin the frame of the dominant waves that the ions interact with [42, 43]. In Fig. 10, these are circles centred on the phase speed of Alfv'enic fluctuations parallel to the magnetic field. At smaller β p , the contours of K ' are nearly vertical for ions with small w ⊥ , with diffusion then leading to an increase in perpendicular energy and a greater measured heating rate (as seen in Figure 3 of the Letter)."}