Ion versus electron heating in compressively driven astrophysical gyrokinetic turbulence

The partition of irreversible heating between ions and electrons in compressively driven (but subsonic) collisionless turbulence is investigated by means of nonlinear gyrokinetic simulations. We derive a prescription for the ion-to-electron heating ratio $Q_{\text{i}}/Q_{\text{e}}$ as a function of the compressive-to-Alfv\'enic driving power ratio $P_{\text{compr}}/P_{\text{AW}}$, of the ratio of ion thermal pressure to magnetic pressure $\beta_{\text{i}}$, and of the ratio of ion-to-electron background temperatures $T_{\text{i}}/T_{\text{e}}$. It is shown that $Q_{\text{i}}/Q_{\text{e}}$ is an increasing function of $P_{\text{compr}}/P_{\text{AW}}$. When the compressive driving is sufficiently large, $Q_{\text{i}}/Q_{\text{e}}$ approaches $\simeq P_{\text{compr}}/P_{\text{AW}}$. This indicates that, in turbulence with large compressive fluctuations, the partition of heating is decided at the injection scales, rather than at kinetic scales. Analysis of phase-space spectra shows that the energy transfer from inertial-range compressive fluctuations to sub-Larmor-scale kinetic Alfv\'en waves is absent for both low and high $\beta_{\text{i}}$, meaning that the compressive driving is directly connected to the ion entropy fluctuations, which are converted into ion thermal energy. This result suggests that preferential electron heating is a very special case requiring low $\beta_{\text{i}}$ and no, or weak, compressive driving. Our heating prescription has wide-ranging applications, including to the solar wind and to hot accretion disks such as M87 and Sgr A*.


I. INTRODUCTION
Most astrophysical systems, e.g., the solar wind, lowluminosity accretion disks, supernova remnants, and the intracluster medium, are in a collisionless turbulent state. The turbulent fluctuations are generally driven by a large-scale freeenergy source that is specific to each system. These fluctuations are cascaded to small scales via nonlinear interactions, and they are converted ultimately into thermal energy. This process is called turbulent heating. In a collisionless plasma, heat is generally deposited into ions and electrons unequally, resulting in a two-temperature state, e.g., in the solar wind [1], accretion disks around black holes [2,3], and the intracluster medium [4]. The partition of turbulent energy between ions and electrons is key to understanding many astrophysical phenomena. Particularly, in the context of accretion disks around black holes, determining the ion-to-electron heating ratio Q i /Q e is critical for interpreting radio images from the Event Horizon Telescope (EHT) [5]. While a recent EHT observation was reproduced numerically using generalrelativistic magnetohydrodynamic (GRMHD) simulations [6], the results strongly depend on the Q i /Q e prescription used (see [7][8][9][10] for the GRMHD simulations with different models of Q i /Q e ). Thus, a physical determination of Q i /Q e is required. * kawazura@tohoku.ac.jp Kinetic, rather than fluid, models must be used in order to calculate correctly the heating rates in a weakly collisional plasma. For the last few years, turbulent heating has been studied by means of particle-in-cell [11][12][13][14][15][16][17] and gyrokinetic (GK) [18][19][20][21] simulations. In these kinetic simulations, turbulence is excited by injection of artificially configured boxscale fluctuations. Such box-scale fluctuations are intended to mimic the fluctuations that cascade from larger scales. In most of the kinetic simulations referenced above [22], the boxscale fluctuations were Alfvénic, meaning that the inertialrange turbulence was assumed to be predominantly Alfvénic. Spacecraft measurements of the solar wind are qualitatively consistent with this assumption, with less than ten percent of the power contained in compressive (slow-mode-like) fluctuations in the inertial range [23,24]. However, there is no guarantee that inertial-range fluctuations of turbulence in other astrophysical systems are predominantly Alfvénic. For example, in our recent study of turbulence driven by the toroidal magnetorotational instability, we found that the Alfvénic and compressive fluctuations are nearly equipartitioned [25].
In this paper, we employ nonlinear GK [26,27] simulations to calculate Q i /Q e in collisionless, subsonic turbulence driven by a mixture of externally injected compressive and Alfvénic fluctuations (implicit in our use of GK is the assumption that the turbulent fluctuations of interest have sub-Larmor frequencies and small amplitudes; hence cyclotronresonance heating [28] and stochastic heating [29,30] are ignored). We use slow-mode-like fluctuations to drive the compressive component of the turbulent cascade. In our previ-ous, purely Alfvénically driven GK simulations [21], we determined the dependence of Q i /Q e on the ratio of the ion thermal pressure to the magnetic pressure, β i = 8πn i T i /B 2 0 , and on the ion-to-electron temperature ratio T i /T e . We found that Q i /Q e was an increasing function of β i , while the dependence on T i /T e was weak (similar to the result arising from linear analysis of Landau/Barnes damping [3,31]). In this work, we determine the dependence of Q i /Q e on the ratio of the compressive to Alfvénic injection power P compr /P AW . We also investigate the properties of the phase-space spectra to understand the heating mechanisms related to the compressive cascade.

II. HYBRID GYROKINETIC MODEL
We solve a hybrid-GK model in which ions are gyrokinetic, while electrons are treated as a massless, isothermal fluid [32]: The electromagnetic fields are determined via the quasineutrality condition and the (parallel and perpendicular) Ampère's law: where e is the elementary charge, Ze is the ion charge, c is the speed of light, B 0 is the ambient magnetic field, z is the coordinate along B 0 , (x, y) is the plane perpendicular to B 0 , v is the particle velocity, F i is the ion equilibrium distribution function, assumed to be Maxwellian, δ f i = h i − Zeφ/T i is the perturbed ion distribution function, n i and T i = m i v 2 thi /2 are the ion density and temperature associated with F i , χ = φ−v·A/c is the GK potential, C[. . . ] is the Coulomb collision operator, . . . R is the gyroaverage at fixed gyrocenter position R, . . . r is the gyroaverage at fixed particle position r, δn e is the electron density perturbation, u e is the parallel electron flow velocity, n e and T e are electron equilibrium density and temperature, φ is the perturbed electrostatic potential, A is the parallel component of the perturbed vector potential, d/dt = ∂ t + (c/B 0 ){φ, . . . }, ∇ = ∂ z − (1/B 0 ){A , . . . }, and { f, g} = (∂ x f )(∂ y g) − (∂ x g)(∂ y f ). The remaining symbols follow standard notation. The compressive fluctuations are driven by an external parallel acceleration a ext in the ion-GK equation (1) [33], while the Alfvénic fluctuations are driven by an external current j ext in the parallel Ampère's law (5) [18-21, 34, 35]. We consider an electron-proton plasma (Z = 1).
The energy budget of the hybrid-GK system is where δn e n e 2 + |δB| 2 8π is the free energy, is the Alfvénic injection power, is the compressive injection power, and is the ion heating rate [33]. The electron heating rate Q e is calculated via the hyperresistive and hyperviscous dissipation of the isothermal electron fluid, which are added to Eqs.
(2) and (3), respectively [36]. In a statistically steady state, P AW + P compr = Q i + Q e , where each term is time averaged. This hybrid model is valid at k ⊥ ρ −1 e . When k ⊥ ρ −1 i , the system follows the equations of kinetic reduced MHD (RMHD) wherein compressive fluctuations are passively advected by the Alfvénic ones ("Alfvén waves", AW), and the two types of fluctuations are energetically decoupled [32,37]. The free energy (8), therefore, can be split as W tot = W AW + W compr , where where δE ⊥ is the fluctuating perpendicular electric field, v A is the Alfvén speed, and g i = δ f i R . In the RMHD range, Alfvénic fluctuations follow fluid equations, whereas the compressive fluctuations are determined by the ion drift kinetic equation [32,37]. Therefore, only ion heating can occur through the phase mixing of compressive fluctuations. When ρ −1 i k ⊥ ρ −1 e , the system follows kinetic electron RMHD (ERMHD) [32], which includes two types of fluctuations, ion entropy fluctuations and kinetic AWs (KAWs) [38].
These fluctuations are again decoupled, and the former are passively advected by the latter. While the KAWs are ultimately dissipated into electron thermal energy, the ion entropy fluctuations lead to ion heating through phase mixing [32].
There are two types of phase mixing in the GK approximation that cause heating: linear Landau/Barnes damping [39,40] and nonlinear phase mixing [32,[41][42][43]. The former creates small-scale structure of the distribution function in the v direction of velocity space, which is thermalized via v derivatives in the collision operator C. The nonlinear phase mixing creates small-scale structure in v ⊥ , and the v ⊥ derivatives in C cause ion heating. Previous Alfvénic-turbulence simulations showed that ion heating occurs in the ERMHD range exclusively via nonlinear phase mixing for low to modest β i [19][20][21]37], while at high β i , there is finite ion heating at k ⊥ ρ −1 i via linear Landau damping [21]. We shall see shortly how this scenario is amended when there is compressive driving.

III. NUMERICAL SETUP
We solve the hybrid-GK model using the AstroGK code [36,44] with two sizes of the simulation domain: the "fiducial" box 0.125 ≤ k x ρ i , k y ρ i ≤ 5.25 and the "doublesized" box 0.0625 ≤ k x ρ i , k y ρ i ≤ 5.25. The grid resolution of the phase-space is (n x , n y , n z , n λ , n ε ) = (128, 128, 32, 32, 16) for the fiducial box, where λ = v 2 ⊥ /v 2 is the pitch angle, and ε = v 2 /2 is the particle's kinetic energy. A recursive expansion procedure [45] is employed to reduce the numerical cost of achieving a statistically steady state.
An oscillating Langevin antenna [35] is employed to drive the Alfvénic and compressive fluctuations. We choose (k x /k x0 , k y /k y0 , k z /k z0 ) = (1, 0, ±1) and (0, 1, ±1) for the driving modes (k 0 is the box-size wave number), 0.9ω A0 for the driving frequency (ω A0 is the box-size Alfvén frequency), and 0.6ω A0 for the decorrelation rate. The amplitude of the Alfvén antenna and, therefore, the power of the Alfvénic driving, P AW , is tuned so that critical balance [46] holds at the box scale [35]. We set the same frequency for the compressive driving and Alfvénic driving because the compressive fluctuations are passively advected by AWs in the RMHD range.
The ion entropy fluctuations are dissipated by the ion collision operator C. In our code, we employ a fully conservative linearized collision operator [47,48] and set the collision frequency to 0.005ω A0 , meaning that ions are almost collisionless. Since the spatial resolution of our simulation is not sufficient to dissipate all of the ion entropy fluctuations via collisions, we add to C a hypercollisionality term proportional to k 8 ⊥ . Its contribution to ion heating is added to Eq. (11). For the dissipation of KAWs, we employ hyperresistivity and hyperviscosity terms proportional to k 8 ⊥ [36] in the isothermal electron fluid (2) and (3).
Given this setup, the free parameters are β i , T i /T e , and the relative amplitude of the compressive driving, which sets P compr /P AW . We investigate β i = (0.1, 1, 4) and T i /T e = (1, 10). For each case, we consider a range of values of P compr /P AW .
IV. ION VS. ELECTRON HEATING Figure 1 shows the dependence of Q i /Q e on P compr /P AW for various values of (β i , T i /T e ). When P compr /P AW = 0, we recover our previous Alfvénic results [21]. When compressive driving is present, Q i /Q e is an increasing function of P compr /P AW for all sets of (β i , T i /T e ) that we investigated. When β i = 0.1, Q i /Q e = P compr /P AW holds for all P compr /P AW , meaning that all of the compressive power is converted into ion heating, and all Alfvénic power is converted into electron heating. This result was theoretically predicted in [33], and is easy to understand physically: when β i 1, ions are too slow to resonate with AWs, and so the Alfvénic cascade goes from the RMHD to ERMHD (sub-ρ i ) regime without losing power and then gets dissipated on electrons. What is both new and surprising in our present numerical result is that, even for β i > 1, Q i /Q e approaches P compr /P AW when P compr /P AW is large. In other words, regardless of β i , almost all the compressive fluctuations in the inertial range are converted into ion heat, if the compressive fluctuations are sufficiently large compared to the Alfvénic fluctuations.
The comparison of the left and right panels in Fig. 1 suggests that Q i /Q e does not depend on T i /T e , which already has been seen for the purely Alfvénic case [21]; here we find that it appears to be true also for the compressively driven case. Admittedly, only two T i /T e cases (T i /T e = 1 and 10) have been investigated in our simulation campaign. To confirm the insensitivity of the heating ratio to T i /T e beyond reasonable doubt, a more extensive scan in T i /T e is needed. An investigation of the T i /T e 1 case is especially important because there is a theoretical expectation of Q i /Q e → 1 for any P compr /P AW when T i /T e 1 and β i 1, namely in the Hall limit [33]. Therefore, the weak dependence on T i /T e that is suggested by the present simulations covers only the values T i /T e 1.
Summarizing the parameter dependences that we have (a-1) FIG. 2. Top row: power spectra of the electric and magnetic fields, normalized by the total perpendicular magnetic energy. Second row: ratios of compressive-field spectra to the perpendicular-magnetic-field spectrum. The horizontal lines correspond to the theoretical predictions for KAWs [32] (dotted lines for (15) and dashed lines for (16)). Third row: spectrum of the ion heating rate normalized by the total heating rate. Bottom row: the ion-heating-rate spectrum integrated up to k ⊥ and normalized by Q i . Parameter values: (a-1)-(c-4) β i = 1, (d-1)-(e-4) β i = 4, and P compr /P AW is increased from left to right. The gray shaded region contains the corner modes in the (k x , k y ) plane. The results shown in (a-1)-(d-4) and (e-1)-(e-4) are from simulations done in the "fiducial" and "double-sized" boxes, respectively (see Sec. III).
found, we propose a simple fitting formula for Q i /Q e : The first term is our previous purely-Alfvénic formula [21]. One finds that Q i /Q e ≥ 1, when P compr /P AW ≥ 1 for any β i and T i /T e ; the implication is that preferential electron heating occurs only for Alfvénic-dominated turbulence at low β i .

V. POWER SPECTRA
In order to investigate the nature of our simulated turbulence, we plot its free energy spectra in the top row of panels of Fig. 2. The energy spectrum of each integrand in Eqs. (12) and (13) is denoted by E with a corresponding subscript. We start by looking at the case of purely Alfvénic driving (P compr /P AW = 0). As expected, the compressive field, δB , is negligible compared to the Alfvénic fields, δB ⊥ and δE ⊥ , in the RMHD range. Alfvénic and compressive fluctuations merge at k ⊥ ρ i ∼ 1 and are reorganized into KAWs and ion entropy fluctuations in the sub-ρ i range. In the RMHD range, the spectra of AW turbulence are E δB ⊥ ∼ E δE ⊥ ∼ k −5/3 ⊥ , while in the sub-ρ i range, the spectra are E δB ⊥ ∼ k −7/3 ⊥ and E δE ⊥ ∼ k −1/3 ⊥ , which match the standard predictions for KAW turbulence [32]. We are not primarily interested in the accuracy of the spectral slopes because the dynamic ranges of either AW or KAW cascades in our simulations are not wide, so these results are not to be viewed as a contribution to the -5/3 vs. -3/2 [49] or the -7/3 vs. -8/3 [50] debates.
The panels in the second row of Fig. 2 show the spectral ratios: E δB , E δn e , E g i , and E u i divided by E δB ⊥ (E u i is the power spectrum of m i n i u 2 i /2). In the ERMHD range, the theoretical predictions based on the linear response for KAW [32], are quite accurately satisfied. Furthermore, u i rapidly drops in the sub-ρ i range, which is also consistent with the KAW turbulence theory, where u i = 0 [32]. While the transition from AW to KAW turbulence is transparent at k ⊥ ρ i 1 for β i = 1, the AW scaling starts to break at k ⊥ ρ i 0.5, and then KAW scaling starts at k ⊥ ρ i 2 for β i = 4. This "intermediate" range at high β i was discovered in our previous purely Alfvénic β i = 100 simulation [21].
Next, we examine how the spectra change when the compressive driving is present. We start by focusing on the RMHD range. As the compressive driving increases, the amplitudes of the compressive fields increase. One finds that the amplitude of u i increases more rapidly than those of δB and δn e , and dominates E g i in the RMHD range. This is because we drive the compressive fluctuations through an external parallel acceleration of ions, a ext [see Eq.
(1)]. The amplitude of g i is much greater than those of δB and δn e when the compressive driving is large, meaning that the compressive driving primarily goes to g i as it includes the contribution from u i . On the other hand, examining the top panels of Fig. 2, one finds that the Alfvénic fields do not change as P compr /P AW increases, indicating that the compressive driving does not contaminate the Alfvénic fields and confirming that the compressive and Alfvénic fields are indeed decoupled in the RMHD range. While this is a theoretical result that has been accepted for some time [32], Fig. 2 appears to be the first confirmation of it based on a hybrid-GK simulation.
Let us now examine the effect of compressive driving on the sub-ρ i -range cascade. Even with sufficiently large P compr /P AW , the spectra of the KAW fields, E δB ⊥ and E δE ⊥ , do not change. The absolute values of spectral amplitude are also preserved. Therefore, the effect of the compressive driving on KAWs is minor. In contrast, E g i increases at all scales as P compr /P AW increases. This result means that the compressive fluctuations in the RMHD range are directly connected to the ion entropy fluctuations in the sub-ρ i range, while the connection with KAWs appears to be absent. If there were an energytransfer path from the RMHD-range compressive fluctuations to KAWs, the amplitudes of KAWs in the compressively driven case would be larger than those in the purely Alfvénic case because E δB ⊥ and E δE ⊥ are proportional to ε 2/3 KAW , where ε KAW is the energy flux of the KAW cascade [32]. Nonetheless, the comparison of Fig. 2 (d-1) and (e-1) shows that E δB ⊥ and E δE ⊥ in the compressively driven case are less than double the purely Alfvénic ones even for P compr /P AW larger than 30. In the low-β i regime, the absence of a path between the inertial-range compressive fluctuations and KAWs was analytically proven in [33]. Here, even at β i = 4, we find that compressive driving affects only the ion-entropy fluctuations. This is the reason why Q i /Q e P compr /P AW is satisfied for P compr /P AW 1 even at β i 1. The panels in the third row of Fig. 2 show the spectrum of the ion heating rate. For β i = 1 and P compr /P AW = 0, most of the ion heating occurs at sub-ρ i scales. This heatingrate spectrum is consistent with the full GK simulation at the same parameters, spanning both the ion and electron kinetic scales [18][19][20]. As P compr /P AW increases, the heating rate both in the RMHD range and at sub-ρ i scales increases. For β i = 4 and P compr /P AW = 0, there is ion heating in the RMHD range with comparable amplitude to the sub-ρ i heating. Similar to the β i = 1 case, the heating rate both in the RMHD range and at sub-ρ i scales increases as P compr /P AW increases.
We note that ion heating near the injection scale may be an artifact when the compressive driving is present: recent drift-kinetic simulations [37] showed that compressive driving directly heated the ions at the injection scale because the turbulent cascade was not yet well developed at that scale. However, in our simulations, the contribution of the heating at the injection scale to the total heating rate is negligible. To show this, we plot, in the bottom panels of Fig. 2, the ionheating-rate spectrum integrated up to k ⊥ and normalized by Q i , viz., . This is the fraction of ion heating rate contained at the scales larger than k −1 ⊥ . We find for all cases, most of the ion heating (∼80%) occurs at sub-ρ i scales. While the compressive driving increases the heating rate both in the RMHD and sub-ρ i ranges (the third row of Fig. 2), the contribution to the total ion heating is predominantly from the sub-ρ i range. It is also evident that the (possibly artificial) box-scale heating in the presence of the compressive driving is negligible, being only 5% of the total.

VI. VELOCITY-SPACE STRUCTURE
In order to investigate the heating process, we show the velocity-space structure of g i . We are particularly interested in the small-scale structures of g i in velocity space as they are the route to heating (i.e., to activating the collision operator) in weakly collisional plasmas [41,51]. Figure 3 shows snapshots of g i and g i /F i in the z = 0 plane for zero and large compressive driving when β i = 0.1 and 4. The normalization by F i helps accentuate the structure at large |v| [52]. In all panels, the top half is taken at k x = k y = 0.375ρ −1 i (RMHD range), and the bottom half is taken at k x = k y = 5.25ρ −1 i (ERMHD range). A rough trend is common for both low and high β i , with and without the compressive driving: in the RMHD range, g i has small-scale structure in the v direction and little structure in the v ⊥ direction; in contrast, in the ERMHD range, there is small-scale structure both in v and v ⊥ . The small-scale structure in v is due to linear Landau damping [39,53,54]; the small-scale structure in v ⊥ is created by nonlinear phase mixing [32,[41][42][43].
In order to investigate quantitatively the heating mechanism, we examine the Hermite and Laguerre spectra [21,37,[54][55][56][57][58][59][60] of g i , viz., |ĝ m, | 2 , defined bŷ where H m (x) and L (x) are Hermite and Laguerre polynomials, respectively. The top panels of Fig. 4 show the Hermite spectra in the RMHD and ERMHD (sub-ρ i ) ranges when the compressive driving is on or off for β i = 0.1 and 4. The Hermite spectrum quantifies the filamentation in v and indicates whether Landau damping is significant or not: the signature of Landau damping is m −1/2 [54]; a steeper spectrum, which in our simulations is measured to be m −1 (cf. [56,60]), may be an indication that Landau damping (phase mixing) is suppressed by the stochastic echo effect [21,37,56,60]. We find that, at both high and low β i , the compressive driving does not change the Hermite spectral slope, viz., m −1 for β i = 0.1 and m −1/2 for β i = 4 in the RMHD range and m −1/2 both for β i = 0.1 and β i = 4 in the ERMHD range. Therefore, regardless of whether the compressive driving exists or not, ion Landau damping is suppressed for β i = 0.1 but is active for β i = 4 in the RMHD. Note, however, that the m −1/2 spectrum in the ERMHD range should be viewed subject to the following caveat. Since there is small-scale structure both in v and v ⊥ directions in ERMHD, and we use (λ, ε) grid rather than (v , v ⊥ ) grid, the small scale structure in v ⊥ may contaminate the Hermite spectrum, and thus the m −1/2 spectrum may turn out to be a numerical artifact. Higher velocity-space resolution (currently too expensive) is necessary to determine if this is the case.
The bottom panels of Fig. 4 show the Laguerre spectrum, which quantifies the filamentation in v ⊥ and is, thus, a diagnostic of nonlinear phase mixing. In contrast to the Hermite spectrum, the Laguerre spectrum in the RMHD range is noticeably modified by compressive driving; for both β i = 0.1 and 4, the Laguerre spectrum becomes shallower when the compressive driving is present. This result indicates that the additional heating in the RMHD range due to compressive driving [ Fig. 2 (b-3), (c-3), and (e-3)] is caused by the emergence of small-scale structures in v ⊥ , presumably triggered by nonlinear phase mixing. Whereas nonlinear phase mixing has been considered to start at k ⊥ ρ i ∼ 1 in Alfvénic turbulence, we find that RMHD-range compressive fluctuations triggers nonlinear phase mixing at k ⊥ ρ i 1. We believe that this is due to the effect of ∇δB drifts [27] but leave further investigation of this detail to future work. In the ERMHD range, on the other hand, compressive driving does not change the Laguerre spectrum. For both β i = 0.1 and 4, the Laguerre spectrum is shallower in the ERMHD range than that in the RMHD range, indicating that the ion heating in the ERMHD range is mediated by the nonlinear phase mixing, as indeed expected theoretically [32].

VII. CONCLUSIONS
In this paper, we have obtained the ion-to-electron irreversible-heating ratio Q i /Q e in compressively driven (but subsonic) gyrokinetic turbulence. Summarizing the dependence on the free parameters, Q i /Q e is (i) an increasing function of P compr /P AW , (ii) an increasing function of β i , and (iii) almost independent of T i /T e . With regard to (i), Q i /Q e P compr /P AW for any β i when the compressive driving is sufficiently large. This result suggests that preferential electron heating, Q i /Q e 1, occurs only when β i 1 and P compr /P AW 1, a fairly special case. A very simple fit- The crosses (circles) correspond to the cases without (with) compressive driving. m and stand for the Hermite and Laguerre moments, respectively. The spectra are integrated over z and v ⊥ for the Hermite spectrum and over z and v for the Laguerre spectrum. For the Hermite spectra, the auxiliary lines m −1 (suggesting suppressed Landau damping/phase mixing [37,56,60]) and m −1/2 (suggesting strong Landau damping [56]) are are shown for reference.
ting formula for the heating ratio is presented in Eq. (14) and is shown to work remarkably well by Fig. 1. This function can be useful in modeling a variety of astrophysical systems because it is applicable to any collisionless turbulent system where the GK ordering holds. Examples of possible applications are the solar wind, AGN jets [61,62], and accretion disks around black holes. Especially for accretion disks, Q i /Q e is important for interpreting observations by the EHT. We note that the parameter sets used for determining our Q i /Q e function are limited, i.e., β i = (0.1, 1, 4) and T i /T e = (1, 10). A wider parameter scan is necessary to extend our prescription Eq. (14) beyond this range, e.g., to the Hall limit, T i /T e 1 and β i 1, which may be a special case [33].
We have also analyzed the phase-space spectra of our turbulence to quantify the distribution, and flows, of free energy. The spectra show that compressive driving affects the com-pressive fluctuations in the RMHD range and the ion entropy fluctuations in the sub-ρ i range, while AWs in the RMHD range and KAWs in the sub-ρ i range are unaffected. This result indicates that compressively injected energy is predominantly converted to ion heating. The spectra of the ion heating rate (Fig. 2) show that most heating happens in the sub-ρ i range, regardless of whether compressive driving is applied or not. The analysis of the ion distribution function and its velocity-space spectra quantifies various phase mixing processes, which are routes to free energy thermalization. We have found that compressive driving does not change the linear phase mixing in the RMHD range, viz., the presence (absence) of phase mixing at high (low) β i ; however a new channel of heating through the enhanced nonlinear phase mixing in the RMHD range emerges when compressive driving is present. While most of these results conform to theoretical expectations [32,33,56], ours appears to be the first study in which some of them have received their numerical corroboration.
In order for result like those reported here to be useful in large-scale modelling, the modeller must know how the turbulent energy injected into their plasma system at large (system-size) scales is partitioned into Alfvénic and compressive (slow-wave-like) cascades in the inertial range. This is an unsolved problem in the majority of astrophysical contexts, but it is solvable one: such a partition is decided at fluid (MHD) rather than kinetic scales. We hope to present a solution to this problem for turbulence driven by the magnetorotational instability [63] with near-azimuthal mean magnetic field in a forthcoming publication [25].