Ion versus electron heating in compressively driven astrophysical gyrokinetic turbulence

Most astrophysical systems, e.g., the solar wind, low-luminosity accretion disks, supernova remnants, and the intracluster medium, are in a collisionless turbulent state. The turbulent fluctuations are generally driven by a large-scale free-energy 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_{\mathrm{i}}/Q_{\mathrm{e}}$ is critical for interpreting radio images from the Event Horizon Telescope (EHT) [5]. While a recent EHT observation was reproduced numerically using general-relativistic magnetohydrodynamic (GRMHD) simulations [6], the results strongly depend on the $Q_{\mathrm{i}}/Q_{\mathrm{e}}$ prescription used (see [7, 8, 9, 10] for the GRMHD simulations with different models of $Q_{\mathrm{i}}/Q_{\mathrm{e}}$). Thus, a physical determination of $Q_{\mathrm{i}}/Q_{\mathrm{e}}$ is required.

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, 18] and gyrokinetic (GK) [19, 20, 21, 22] simulations. In these kinetic simulations, turbulence is excited by injection of artificially configured box-scale fluctuations. Such box-scale fluctuations are intended to mimic the fluctuations that cascade from larger scales. In most of the kinetic simulations referenced above ¹¹1With the exception of those of the freely decaying whistler turbulence [12], the box-scale fluctuations were Alfvénic, meaning that the inertial-range 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 [24, 25]. 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 [26].

In this paper, we employ nonlinear GK [27, 28] simulations to calculate $Q_{\mathrm{i}}/Q_{\mathrm{e}}$ in collisionless, subsonic turbulence driven by a mixture of externally injected compressive and Alfvénic fluctuations. The GK theory shows that the energy partition between ions and electrons is decided around the ion Larmor scale [29], meaning that the energy that is not destined for the ion heating at the ion Larmor scale will be channeled into the electron heating. Therefore, the amount of electron heating can be computed as $Q_{\mathrm{e}}=P_{\mathrm{AW}}+P_{\mathrm{compr}}-Q_{\mathrm{i}}$, where $P_{\mathrm{AW}}$ and $P_{\mathrm{compr}}$ are the Alfvénic and compressive energy injection. Since the electron kinetic effects are not necessary to obtain the electron heating, we use a hybrid-GK model in which electrons are treated as a massless, isothermal fluid [29].

To drive the compressive component of the turbulent cascade, we use slow-mode-like fluctuations. In our previous, purely Alfvénically driven GK simulations [22], we determined the dependence of $Q_{\mathrm{i}}/Q_{\mathrm{e}}$ on the ratio of the ion thermal pressure to the magnetic pressure, $\beta_{\mathrm{i}}=8\pi n_{\mathrm{i}}T_{\mathrm{i}}/B_{0}^{2}$, and on the ion-to-electron temperature ratio $T_{\mathrm{i}}/T_{\mathrm{e}}$. We found that $Q_{\mathrm{i}}/Q_{\mathrm{e}}$ was an increasing function of $\beta_{\mathrm{i}}$, while the dependence on $T_{\mathrm{i}}/T_{\mathrm{e}}$ was weak (similar to the result arising from linear analysis of Landau/Barnes damping [3, 30]). In this work, we determine the dependence of $Q_{\mathrm{i}}/Q_{\mathrm{e}}$ on the ratio of the compressive to Alfvénic injection power $P_{\mathrm{compr}}/P_{\mathrm{AW}}$. We also investigate the properties of the phase-space spectra to understand the heating mechanisms related to the compressive cascade.

We solve a hybrid-GK model in which ions are gyrokinetic, while electrons are treated as a massless, isothermal fluid [29]:

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, ${\mathbf{B}}_{0}$ is the ambient magnetic field, $z$ is the coordinate along ${\mathbf{B}}_{0}$, $(x,y)$ is the plane perpendicular to ${\mathbf{B}}_{0}$, ${\mathbf{v}}$ is the particle velocity, $F_{\mathrm{i}}$ is the ion equilibrium distribution function, assumed to be Maxwellian, $\delta f_{\mathrm{i}}=h_{\mathrm{i}}-Ze\phi/T_{\mathrm{i}}$ is the perturbed ion distribution function, $n_{\mathrm{i}}$ and $T_{\mathrm{i}}=m_{\mathrm{i}}v_{\mathrm{thi}}^{2}/2$ are the ion density and temperature associated with $F_{\mathrm{i}}$, $\chi=\phi-{\mathbf{v}}\cdot{\mathbf{A}}/c$ is the GK potential, $C[\dots]$ is the Coulomb collision operator, $\left<\dots\right>_{\mathbf{R}}$ is the gyroaverage at fixed gyrocenter position ${\mathbf{R}}$, $\left<\dots\right>_{\mathbf{r}}$ is the gyroaverage at fixed particle position ${\mathbf{r}}$, $\delta n_{\mathrm{e}}$ is the electron density perturbation, $u_{\|\mathrm{e}}$ is the parallel electron flow velocity, $n_{\mathrm{e}}$ and $T_{\mathrm{e}}$ are electron equilibrium density and temperature, $\phi$ is the perturbed electrostatic potential, $A_{\|}$ is the parallel component of the perturbed vector potential, $\mathrm{d}/\mathrm{d}t=\partial_{t}+(c/B_{0})\{\phi,\,\dots\}$, $\nabla_{\|}=\partial_{z}-(1/B_{0})\{A_{\|},\,\dots\}$, and $\{f,\,g\}=(\partial_{x}f)(\partial_{y}g)-(\partial_{x}g)(\partial_{y}f)$. The remaining symbols follow standard notation. The compressive fluctuations are driven by an external parallel acceleration $a_{\mathrm{ext}}$ in the ion-GK equation (1) [31], while the Alfvénic fluctuations are driven by an external current $j_{\mathrm{\|ext}}$ in the parallel Ampère’s law (5) [32, 19, 33, 20, 21, 22]. We consider an electron-proton plasma ($Z=1$).

The energy budget of the hybrid-GK system is

where

is the free energy,

is the Alfvénic injection power,

is the compressive injection power, and

is the ion heating rate [31]. The electron heating rate $Q_{\mathrm{e}}$ is calculated via the hyperresistive and hyperviscous dissipation of the isothermal electron fluid, which are added to Eqs. (2) and (3), respectively [34]. In a statistically steady state, $P_{\mathrm{AW}}+P_{\mathrm{compr}}=Q_{\mathrm{i}}+Q_{\mathrm{e}}$, where each term is time averaged.

This hybrid model is valid at $k_{\perp}\ll\rho_{\mathrm{e}}^{-1}$. When $k_{\perp}\ll\rho_{\mathrm{i}}^{-1}$, 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 [29, 35]. The free energy (8), therefore, can be split as $W_{\mathrm{tot}}=W_{\mathrm{AW}}+W_{\mathrm{compr}}$, where

where $\delta E_{\perp}$ is the fluctuating perpendicular electric field, $v_{\mathrm{A}}$ is the Alfvén speed, and $g_{\mathrm{i}}=\langle\delta f_{\mathrm{i}}\rangle_{\mathbf{R}}$. In the RMHD range, Alfvénic fluctuations follow fluid equations, whereas the compressive fluctuations are determined by the ion drift kinetic equation [29, 35]. Therefore, in the RMHD range, only ion heating can occur through the phase mixing of compressive fluctuations.

When $\rho_{\mathrm{i}}^{-1}\ll k_{\perp}\ll\rho_{\mathrm{e}}^{-1}$, the system follows kinetic electron RMHD (ERMHD) [29], which includes two types of fluctuations, ion entropy fluctuations and kinetic AWs (KAWs) ²²2 Here and in what follows, we refer, “colloquially”, to the RMHD- and ERMHD-range turbulent fluctuations as AWs or KAWs. This should not be read to imply that we believe turbulence in these regimes to be a collection of random phased weakly interacting waves. Turbulence in these regimes is critically balanced [60] and so always strong, and capable of producing intermittent structures, current sheets, etc. The reference to “waves” simply highlights the fact that even in this strong regime, the linear response relations between fluctuations of different fields, e.g., $\delta E_{\perp}$, $\delta B_{\perp}$, $\delta B_{||}$, $\delta u_{\|\mathrm{i}}$, $\delta n_{\mathrm{e}}$, etc., are of the same physical nature as in AWs or KAWs (and indeed follow those quite closely even quantitatively). This is because critical balance implies that linear and nonlinear physics are always of the same order. A detailed study of this topic can be found in [80]. . 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 [29].

There are two types of phase mixing in the GK approximation that cause heating: linear Landau/Barnes damping [37, 38] and nonlinear phase mixing [39, 29, 40, 41]. 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_{\perp}$, and the $v_{\perp}$ 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 $\beta_{\mathrm{i}}$ [20, 21, 22, 35], while at high $\beta_{\mathrm{i}}$, there is finite ion heating at $k_{\perp}\lesssim\rho_{\mathrm{i}}^{-1}$ via linear Landau damping [22]. We shall see shortly how this scenario is amended when there is compressive driving.

We solve the hybrid-GK model using the AstroGK code [56, 34] with two sizes of the simulation domain: the “fiducial” box $0.125\leq k_{x}\rho_{\mathrm{i}},k_{y}\rho_{\mathrm{i}}\leq 5.25$ and the “double-sized” box $0.0625\leq k_{x}\rho_{\mathrm{i}},k_{y}\rho_{\mathrm{i}}\leq 5.25$. The grid resolution of the phase-space is $(n_{x},n_{y},n_{z},n_{\lambda},n_{\varepsilon})=(128,128,32,32,16)$ for the fiducial box, where $\lambda=v_{\perp}^{2}/v^{2}$ is the pitch angle, and $\varepsilon=v^{2}/2$ is the particle’s kinetic energy. Although this grid resolution is lower than that used in some studies (e.g., [20, 57, 58]), this is the price for being able to carry out an adequately broad parameter scan. We will also simulate a single higher-velocity-space-resolution case $(n_{x},n_{y},n_{z},n_{\lambda},n_{\varepsilon})=(128,128,32,64,32)$ to check the numerical convergence. A recursive expansion procedure [59] is employed to reduce the numerical cost of achieving a statistically steady state.

An oscillating Langevin antenna [33] 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,\pm 1)$ and $(0,1,\pm 1)$ for the driving modes (${\mathbf{k}}_{0}$ is the box-size wave number), $0.9\omega_{\mathrm{A0}}$ for the driving frequency ($\omega_{\mathrm{A0}}$ is the box-size Alfvén frequency), and $0.6\omega_{\mathrm{A0}}$ for the decorrelation rate. The amplitude of the Alfvén antenna and, therefore, the power of the Alfvénic driving, $P_{\mathrm{AW}}$, is tuned so that critical balance [60] holds at the box scale [33]. 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 [61, 62] and set the collision frequency to $0.005\omega_{\mathrm{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_{\perp}^{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_{\perp}^{8}$ [34] in the isothermal electron fluid (2) and (3).

Given this setup, the free parameters are $\beta_{\mathrm{i}}$, $T_{\mathrm{i}}/T_{\mathrm{e}}$, and the relative amplitude of the compressive driving, which sets $P_{\mathrm{compr}}/P_{\mathrm{AW}}$. We investigate $\beta_{\mathrm{i}}=(0.1,1,4)$ and $T_{\mathrm{i}}/T_{\mathrm{e}}=(1,10)$. For each case, we consider a range of values of $P_{\mathrm{compr}}/P_{\mathrm{AW}}$.

Figure 1 shows the dependence of $Q_{\mathrm{i}}/Q_{\mathrm{e}}$ on $P_{\mathrm{compr}}/P_{\mathrm{AW}}$ for various values of $(\beta_{\mathrm{i}},\,T_{\mathrm{i}}/T_{\mathrm{e}})$. When $P_{\mathrm{compr}}/P_{\mathrm{AW}}=0$, we recover our previous Alfvénic results [22]. When compressive driving is present, $Q_{\mathrm{i}}/Q_{\mathrm{e}}$ is an increasing function of $P_{\mathrm{compr}}/P_{\mathrm{AW}}$ for all sets of $(\beta_{\mathrm{i}},\,T_{\mathrm{i}}/T_{\mathrm{e}})$ that we investigated. When $\beta_{\mathrm{i}}=0.1$, $Q_{\mathrm{i}}/Q_{\mathrm{e}}=P_{\mathrm{compr}}/P_{\mathrm{AW}}$ holds for all $P_{\mathrm{compr}}/P_{\mathrm{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 [31], and is easy to understand physically: when $\beta_{\mathrm{i}}\ll 1$, ions are too slow to resonate with AWs, and so the Alfvénic cascade goes from the RMHD to ERMHD 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 $\beta_{\mathrm{i}}>1$, $Q_{\mathrm{i}}/Q_{\mathrm{e}}$ approaches $P_{\mathrm{compr}}/P_{\mathrm{AW}}$ when $P_{\mathrm{compr}}/P_{\mathrm{AW}}$ is large. In other words, regardless of $\beta_{\mathrm{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_{\mathrm{i}}/Q_{\mathrm{e}}$ does not depend on $T_{\mathrm{i}}/T_{\mathrm{e}}$, which already has been seen for the purely Alfvénic case [22]; here we find that it appears to be true also for the compressively driven case. Admittedly, only two $T_{\mathrm{i}}/T_{\mathrm{e}}$ cases ($T_{\mathrm{i}}/T_{\mathrm{e}}=$ 1 and 10) have been investigated in our simulation campaign. Therefore, the weak dependence on $T_{\mathrm{i}}/T_{\mathrm{e}}$ that is suggested by the present simulations may cover only the values $T_{\mathrm{i}}/T_{\mathrm{e}}\gtrsim 1$. In contrast, when $T_{\mathrm{i}}/T_{\mathrm{e}}\ll 1$ and $\beta_{\mathrm{i}}\ll 1$, namely in the Hall limit [31], there is a theoretical expectation of $Q_{\mathrm{i}}/Q_{\mathrm{e}}\to 1$ for any $P_{\mathrm{compr}}/P_{\mathrm{AW}}$. Since this theoretical expectation is inconsistent with Eq. (14), we need to examine $T_{\mathrm{i}}/T_{\mathrm{e}}<1$ to see whether and when the insensitivity of $Q_{\mathrm{i}}/Q_{\mathrm{e}}$ to $T_{\mathrm{i}}/T_{\mathrm{e}}$ breaks.

Summarizing the parameter dependences that we have found, we propose a simple fitting formula for $Q_{\mathrm{i}}/Q_{\mathrm{e}}$:

The first term is our previous purely-Alfvénic formula [22]. One finds that $Q_{\mathrm{i}}/Q_{\mathrm{e}}\geq 1$, when $P_{\mathrm{compr}}/P_{\mathrm{AW}}\geq 1$ for any $\beta_{\mathrm{i}}$ and $T_{\mathrm{i}}/T_{\mathrm{e}}$; the implication is that preferential electron heating occurs only for Alfvénic-dominated turbulence at low $\beta_{\mathrm{i}}$.

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_{\mathrm{compr}}/P_{\mathrm{AW}}=0$). As expected, the compressive field, $\delta B_{\|}$, is negligible compared to the Alfvénic fields, $\delta B_{\perp}$ and $\delta E_{\perp}$, in the RMHD range. Alfvénic and compressive fluctuations merge at $k_{\perp}\rho_{\mathrm{i}}\sim 1$ and are reorganized into KAWs and ion entropy fluctuations in the sub-$\rho_{\mathrm{i}}$ range. In the RMHD range, the spectra of AW turbulence are $E_{\delta B_{\perp}}\sim E_{\delta E_{\perp}}\sim k_{\perp}^{-5/3}$, while in the sub-$\rho_{\mathrm{i}}$ range, the spectra are $E_{\delta B_{\perp}}\sim k_{\perp}^{-7/3}$ and $E_{\delta E_{\perp}}\sim k_{\perp}^{-1/3}$, which match the standard predictions for KAW turbulence [29]. 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 [63] or the -7/3 vs. -8/3 [64] debates.

The panels in the second row of Fig. 2 show the spectral ratios: $E_{\delta B_{\|}}$, $E_{\delta n_{\mathrm{e}}}$, $E_{g_{\mathrm{i}}}$, and $E_{u_{\|\mathrm{i}}}$ divided by $E_{\delta B_{\perp}}$ ($E_{u_{\|\mathrm{i}}}$ is the power spectrum of $m_{\mathrm{i}}n_{\mathrm{i}}u_{\|\mathrm{i}}^{2}/2$). In the ERMHD range, the theoretical predictions based on the linear response for KAW [29],

are quite accurately satisfied. Furthermore, $u_{\|\mathrm{i}}$ rapidly drops in the sub-$\rho_{\mathrm{i}}$ range, which is also consistent with the KAW turbulence theory, where $u_{\|\mathrm{i}}=0$ [29]. While the transition from AW to KAW turbulence is transparent at $k_{\perp}\rho_{\mathrm{i}}\simeq 1$ for $\beta_{\mathrm{i}}=1$, the AW scaling starts to break at $k_{\perp}\rho_{\mathrm{i}}\simeq 0.5$, and then KAW scaling starts at $k_{\perp}\rho_{\mathrm{i}}\simeq 2$ for $\beta_{\mathrm{i}}=4$. This “intermediate” range at high $\beta_{\mathrm{i}}$ was discovered in our previous purely Alfvénic $\beta_{\mathrm{i}}=100$ simulation [22].

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_{\|\mathrm{i}}$ increases more rapidly than those of $\delta B_{\|}$ and $\delta n_{\mathrm{e}}$, and dominates $E_{g_{\mathrm{i}}}$. This is because we drive the compressive fluctuations through an external parallel acceleration of ions, $a_{\mathrm{ext}}$ [see Eq. (1)]. The amplitude of $g_{\mathrm{i}}$ is much greater than those of $\delta B_{\|}$ and $\delta n_{\mathrm{e}}$ when the compressive driving is large, meaning that the compressive driving primarily goes to $g_{\mathrm{i}}$ as it includes the contribution from $u_{\|\mathrm{i}}$. On the other hand, examining the top panels of Fig. 2, one finds that the Alfvénic fields do not change as $P_{\mathrm{compr}}/P_{\mathrm{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 [29], the theoretical prediction is based on an asymptotic expansion in $k_{\perp}\rho_{\mathrm{i}}\ll 1$ and relies on a number of assumptions — most importantly, locality of nonlinear interactions, which is not a completely uncontroversial approach (e.g., [65]). Thus, the numerical confirmation of the decoupling shown in Fig. 2 is a nontrivial result, and a confirmation that a certain way of thinking about plasma turbulence problems is a reasonable one, and that reassuringly, asymptotic theory gives one a decent grasp of the problem even when the small parameter — $k_{\perp}\rho_{\mathrm{i}}$ in this case — is only moderately small.

Let us now examine the effect of compressive driving on the sub-$\rho_{\mathrm{i}}$-range cascade. Even with sufficiently large $P_{\mathrm{compr}}/P_{\mathrm{AW}}$, the spectra of the KAW fields, $E_{\delta B_{\perp}}$ and $E_{\delta E_{\perp}}$, 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_{\mathrm{i}}}$ increases at all scales as $P_{\mathrm{compr}}/P_{\mathrm{AW}}$ increases. This result means that the compressive fluctuations in the RMHD range are directly connected to the ion entropy fluctuations in the sub-$\rho_{\mathrm{i}}$ range, while the connection with KAWs appears to be absent. If there were an energy-transfer path from the inertial-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_{\delta B_{\perp}}$ and $E_{\delta E_{\perp}}$ are proportional to $\varepsilon_{\mathrm{KAW}}^{2/3}$, where $\varepsilon_{\mathrm{KAW}}$ is the energy flux of the KAW cascade [29]. Nonetheless, the comparison of Fig. 2 (d-1) and (e-1) shows that $E_{\delta B_{\perp}}$ and $E_{\delta E_{\perp}}$ in the compressively driven case are less than double the purely Alfvénic ones even for $P_{\mathrm{compr}}/P_{\mathrm{AW}}$ larger than 30. In the low-$\beta_{\mathrm{i}}$ regime, the absence of a path between the inertial-range compressive fluctuations and KAWs was analytically proven in [31]. Here, even at $\beta_{\mathrm{i}}=4$, we find that compressive driving affects only the ion-entropy fluctuations. This is the reason why $Q_{\mathrm{i}}/Q_{\mathrm{e}}\simeq P_{\mathrm{compr}}/P_{\mathrm{AW}}$ is satisfied for $P_{\mathrm{compr}}/P_{\mathrm{AW}}\gg 1$ even at $\beta_{\mathrm{i}}\gtrsim 1$.

The panels in the third row of Fig. 2 show the spectrum of the ion heating rate. For $\beta_{\mathrm{i}}=1$ and $P_{\mathrm{compr}}/P_{\mathrm{AW}}=0$, most of the ion heating occurs at sub-$\rho_{\mathrm{i}}$ scales. This heating-rate spectrum is consistent with the full GK simulation at the same parameters, spanning both the ion and electron kinetic scales [19, 20, 21]. As $P_{\mathrm{compr}}/P_{\mathrm{AW}}$ increases, the heating rate both in the RMHD range and at sub-$\rho_{\mathrm{i}}$ scales increases. For $\beta_{\mathrm{i}}=4$ and $P_{\mathrm{compr}}/P_{\mathrm{AW}}=0$, there is ion heating in the RMHD range with comparable amplitude to the sub-$\rho_{\mathrm{i}}$ heating. The ion heating in $k_{\perp}\rho_{\mathrm{i}}\lesssim 0.3$ is due to the Landau damping of AWs since there are no compressive fluctuations at $k_{\perp}\rho_{\mathrm{i}}\lesssim 0.3$. This indicates that the box scale of $\beta_{\mathrm{i}}=4$ simulation is not precisely asymptotically in the RMHD range even though the electromagnetic spectra at $k_{\perp}\rho_{\mathrm{i}}\lesssim 0.3$ look like those of RMHD turbulence, viz., $E_{\delta B_{\perp}}\sim E_{\delta E_{\perp}}\sim k_{\perp}^{-5/3}$. Similar to the $\beta_{\mathrm{i}}=1$ case, the heating rate both in the RMHD range and at sub-$\rho_{\mathrm{i}}$ scales increases as $P_{\mathrm{compr}}/P_{\mathrm{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 [35] 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 ion-heating-rate spectrum integrated up to $k_{\perp}$ and normalized by $Q_{\mathrm{i}}$, viz., $\int_{k_{{\perp}0}}^{k_{\perp}}\mathrm{d}k_{\perp}^{\prime}\,E_{Q_{\mathrm{i}}}(k_{\perp}^{\prime})/Q_{\mathrm{i}}$, where $k_{{\perp}0}^{2}=k_{x0}^{2}+k_{y0}^{2}$. This is the fraction of ion heating rate contained at the scales larger than $k_{\perp}^{-1}$. We find for all cases, most of the ion heating ($\sim$80%) occurs at sub-$\rho_{\mathrm{i}}$ scales. While the compressive driving increases the heating rate both in the RMHD and sub-$\rho_{\mathrm{i}}$ ranges (the third row of Fig. 2), the contribution to the total ion heating is predominantly from the sub-$\rho_{\mathrm{i}}$ range. It is also evident that the (possibly artificial) box-scale heating in the presence of the compressive driving is negligible, being only $\simeq$5% of the total.

In order to investigate the heating process, we show the velocity-space structure of $g_{\mathrm{i}}$. We are particularly interested in the small-scale structures of $g_{\mathrm{i}}$ in velocity space as they are the route to heating (i.e., to activating the collision operator) in weakly collisional plasmas [66, 39]. Figure 3 shows snapshots of $g_{\mathrm{i}}$ and $g_{\mathrm{i}}/F_{\mathrm{i}}$ in the $z=0$ plane for zero and large compressive driving when $\beta_{\mathrm{i}}=$ 0.1 and 4. The normalization by $F_{\mathrm{i}}$ helps accentuate the structure at large $|{\mathbf{v}}|$ [67]. In all panels, the top half is taken at $k_{x}=k_{y}=0.375\rho_{\mathrm{i}}^{-1}$, and the bottom half is taken at $k_{x}=k_{y}=5.25\rho_{\mathrm{i}}^{-1}$. A rough trend is common for both low and high $\beta_{\mathrm{i}}$, with and without the compressive driving: in the RMHD range, $g_{\mathrm{i}}$ has small-scale structure in the $v_{\|}$ direction and little structure in the $v_{\perp}$ direction; in contrast, in the ERMHD range, there is small-scale structure both in $v_{\|}$ and $v_{\perp}$. The small-scale structure in $v_{\|}$ is due to linear Landau damping [37, 68, 69]; the small-scale structure in $v_{\perp}$ is created by nonlinear phase mixing [39, 29, 40, 41].

In order to investigate quantitatively the heating mechanism, we examine the Hermite and Laguerre spectra [70, 69, 71, 72, 73, 58, 74, 35, 22] of $g_{\mathrm{i}}$, viz., $|\hat{g}_{m,\ell}|^{2}$, defined by

where $H_{m}(x)$ and $L_{\ell}(x)$ are Hermite and Laguerre polynomials, respectively. The top panels of Fig. 4 show the Hermite spectra in the RMHD and ERMHD (sub-$\rho_{\mathrm{i}}$) ranges when the compressive driving is on or off for $\beta_{\mathrm{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}$ [69]; a steeper spectrum, which in our simulations is measured to be $m^{-1}$ (cf. [71, 74]), may be an indication that Landau damping (phase mixing) is suppressed by the stochastic echo effect [71, 74, 35, 22]. We find that, at both high and low $\beta_{\mathrm{i}}$, the compressive driving does not change the Hermite spectral slope, viz., $m^{-1}$ for $\beta_{\mathrm{i}}=0.1$ and $m^{-1/2}$ for $\beta_{\mathrm{i}}=4$ in the RMHD range and $m^{-1/2}$ both for $\beta_{\mathrm{i}}=0.1$ and $\beta_{\mathrm{i}}=4$ in the ERMHD range. Therefore, regardless of whether the compressive driving exists or not, ion Landau damping is suppressed for $\beta_{\mathrm{i}}=0.1$ but is active for $\beta_{\mathrm{i}}=4$ in the RMHD regime ³³3 In our previous paper [22], we noted that the correspondence between the ion heating and $m^{-1/2}$ spectrum must be viewed cautiously at high $\beta_{\mathrm{i}}$ because the stochastic echo effect may not be computed correctly when the effective collisional cutoff $m_{\mathrm{c}}$ is smaller than $\beta_{\mathrm{i}}$. In this study, the cutoff is $m_{\mathrm{c}}\sim 10$ while the maximum $\beta_{\mathrm{i}}$ is 4; hence we consider our simulation to be just about safe from this concern ⁴⁴4 We also note that the hypercollisions in our simulations are effective only around the grid scale and do not play any role in the Hermite spectrum at the RMHD range . 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_{\perp}$ directions in ERMHD, and we use $(\lambda,\,\varepsilon)$ grid rather than $(v_{\|},\,v_{\perp})$ grid, the small scale structure in $v_{\perp}$ 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_{\perp}$ 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 $\beta_{\mathrm{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_{\perp}$, presumably triggered by nonlinear phase mixing. Whereas nonlinear phase mixing has been considered to start at $k_{\perp}\rho_{\mathrm{i}}\sim 1$ in Alfvénic turbulence, we find that RMHD-range compressive fluctuations triggers nonlinear phase mixing at $k_{\perp}\rho_{\mathrm{i}}\ll 1$. We believe that this is due to the effect of $\nabla\delta B_{\|}$ drifts [28] 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 $\beta_{\mathrm{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 [29].

In this paper, we have obtained the ion-to-electron irreversible-heating ratio $Q_{\mathrm{i}}/Q_{\mathrm{e}}$ in compressively driven (but subsonic) gyrokinetic turbulence. Summarizing the dependence on the free parameters, $Q_{\mathrm{i}}/Q_{\mathrm{e}}$ is (i) an increasing function of $P_{\mathrm{compr}}/P_{\mathrm{AW}}$, (ii) an increasing function of $\beta_{\mathrm{i}}$, and (iii) almost independent of $T_{\mathrm{i}}/T_{\mathrm{e}}$. With regard to (i), $Q_{\mathrm{i}}/Q_{\mathrm{e}}\simeq P_{\mathrm{compr}}/P_{\mathrm{AW}}$ for any $\beta_{\mathrm{i}}$ when the compressive driving is sufficiently large. This result suggests that preferential electron heating, $Q_{\mathrm{i}}/Q_{\mathrm{e}}\ll 1$, occurs only when $\beta_{\mathrm{i}}\ll 1$ and $P_{\mathrm{compr}}/P_{\mathrm{AW}}\ll 1$, a fairly special case. A very simple fitting 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, such as the solar wind, AGN jets [77, 78], and accretion disks around black holes. Especially for accretion disks, $Q_{\mathrm{i}}/Q_{\mathrm{e}}$ is important for interpreting observations by the EHT. We offer this prescription to the modelers with a word of caution that our results may not be precisely, quantitatively applicable beyond the limitations discussed in Sec. II.1. We also note that the parameter sets used for determining our $Q_{\mathrm{i}}/Q_{\mathrm{e}}$ function are limited, i.e., $\beta_{\mathrm{i}}=(0.1,1,4)$ and $T_{\mathrm{i}}/T_{\mathrm{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_{\mathrm{i}}/T_{\mathrm{e}}\ll 1$ and $\beta_{\mathrm{i}}\ll 1$, which may be a special case [31].

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 compressive fluctuations in the RMHD range and the ion entropy fluctuations in the sub-$\rho_{\mathrm{i}}$ range, while AWs in the RMHD range and KAWs in the sub-$\rho_{\mathrm{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-$\rho_{\mathrm{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) $\beta_{\mathrm{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 [29, 71, 31], ours appears to be the first study in which some of them have received their numerical corroboration.

In order for results like those reported here to be useful in large-scale modelling, the modeler 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 a 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 [79] with near-azimuthal mean magnetic field in a forthcoming publication [26].