3.1.1. Imbalanced turbulence

Although balanced turbulence simplifies the description of turbulence, many systems are imbalanced with $\varepsilon ^+\gg \varepsilon ^-$ . While there is currently no theory that describes imbalanced turbulence fully, phenomenological arguments can be made to explain qualitatively aspects seen in observations and numerical simulations. Here, we cover a simple version, as needed for later, and refer the reader to Schekochihin (Reference Schekochihin2022) for a complete description of the subtleties involved.

Using the constant-flux argument with (3.1) for $\tau _{\mathrm{nl}}^\pm$ , at each scale $\lambda$ , one assumes

(3.4) \begin{equation} \frac {(\delta z^\pm _\lambda )^2}{\tau _{\mathrm{nl}}^\pm }\sim \frac {(\delta z^\pm _\lambda )^2\delta z^\mp _ \lambda }{\lambda }\sim \varepsilon ^\pm , \end{equation}

where $\varepsilon ^\pm$ are the energy injection rates at the outer scales of the turbulence for the $\boldsymbol{z}^\pm$ fluctuations. The nonlinear time scales (3.1) also imply that

(3.5) \begin{equation} \frac {\tau _{\mathrm{nl}}^+}{\tau _{\mathrm{nl}}^-}\sim \frac {\delta z^+_\lambda }{\delta z^-_\lambda }\gg 1; \end{equation}

in other words, the cascade of the stronger field becomes less efficient with increasing imbalance as it is dependent on the weaker field to advect it. Although $\delta z^-_\lambda \ll \delta z^+_\lambda$ , the cascades can still be assumed strong, with $\tau _{\mathrm{nl}}^\pm \lesssim \tau _{\mathrm{A}}$ individually.

The constant-flux argument aligns with results from numerical simulations that the ratio of the energies of the individual fields scales with imbalance as $\alpha \sim \alpha _\varepsilon ^2$ (Schekochihin Reference Schekochihin2022), where

(3.6) \begin{equation} \alpha _\varepsilon \equiv \frac {\varepsilon ^+}{\varepsilon ^-} \end{equation}

is the injection ratio and

(3.7) \begin{equation} \alpha \equiv \frac {E^+}{E^-}=\frac {1+\sigma _{\mathrm{c}}}{1-\sigma _{\mathrm{c}}} \end{equation}

is the ratio of Elsasser energies (assumed to be independent of scale such that $\alpha \sim (\delta z^+_\lambda )^2/(\delta z^-_\lambda )^2$ ).

Despite this agreement with numerics, the fact that the weaker $\boldsymbol{z}^-$ field is able to advect and distort the stronger $\boldsymbol{z}^+$ field over the long time scale $\tau _{\mathrm{nl}}^+$ requires a separation of its correlation and spatial-distortion time scales (Lithwick, Goldreich & Sridhar Reference Lithwick, Goldreich and Sridhar2007). The stronger field $\boldsymbol{z}^+$ distorts the weaker $\boldsymbol{z}^-$ in space over a time scale $\tau _{\mathrm{nl}}^-$ . However, $\boldsymbol{z}^-$ is ‘swept up’ by the stronger $\boldsymbol{z}^+$ fluctuations: in the frame in which $\boldsymbol{z}^+$ fluctuations are stationary, the $\boldsymbol{z}^-$ fluctuations develop fine structures due to the spatial distortion, but remain approximately constant over $\tau _{\mathrm{nl}}^+$ to remain correlated with the $\boldsymbol{z}^+$ fluctuations.

The separation of these time scales means that the causality argument for critical balance in the balanced case (dubbed ‘causal CB’) no longer holds. Directly applying the assumptions of causal CB in the imbalanced case would imply that the parallel correlation lengths of the $\boldsymbol{z}^\pm$ fluctuations are disparate with $\ell _\|^+/\ell _\|^- \gg 1$ . However, Lithwick et al. (Reference Lithwick, Goldreich and Sridhar2007) argue that, because $\boldsymbol{z}^-$ perturbations are spatially distorted over a time $\tau _{\mathrm{nl}}^-$ (such that $\ell _\|^-\sim v_{\mathrm{A}}\tau _{\mathrm{nl}}^-$ ), $\boldsymbol{z}^+$ fluctuations separated by a distance $\ell _\|^-$ are advected by the now spatially decorrelated $\boldsymbol{z}^-$ and thus themselves decorrelate over a similar distance: $\ell _\|^+\sim \ell _\|^-\sim v_{\mathrm{A}}\tau _{\mathrm{nl}}^- \sim k_\perp \delta z^+$ . This effect, termed ‘propagation CB’ (Beresnyak & Lazarian Reference Beresnyak and Lazarian2008), allows a CB argument to be used to balance linear and nonlinear times using the fact that the $\boldsymbol{z}^-$ fluctuations are swept along by $\boldsymbol{z}^+$ fluctuations, rather than information travelling between points.

Although the theory outlined previously does well in explaining some results observed in simulations of imbalanced turbulence, it is far from a complete picture. The scaling $\alpha \sim \alpha _\varepsilon ^2$ appears to be robust, observed in multiple simulations of RMHD turbulence (Schekochihin Reference Schekochihin2022), but does not agree with alternative theories of imbalanced turbulence (Perez & Boldyrev Reference Perez and Boldyrev2009). The recent theory of Schekochihin (Reference Schekochihin2022), which assumes non-local nonlinear interactions, appears to account for some of these discrepancies, but further work is needed to test this thoroughly. Regardless, the phenomenology of propagation CB presented here is adequate for the model we use to investigate ion heating, with the simulations of Appendix A showing modest support for this theory.

3.1.2. Imbalance dependence of the wavevector–frequency spectrum

The form of the wavevector–frequency spectrum depends on the type of fluctuations present in the turbulence and the strength of the nonlinear interactions between them. Unbroadened fluctuations – waves following a dispersion relation $\omega(\boldsymbol{k})$ with a single frequency $\omega$ at a particular $\boldsymbol{k}$ – are sharply peaked along this dispersion relation in the sense of a delta function $\delta[\omega-\omega(\boldsymbol{k})]$ in the wavevector–frequency spectrum. In contrast, due to the nonlinear interactions between different wavenumbers in the energy cascade, fluctuations in a strongly turbulent system take on a spread of frequencies of width $\Delta \omega$ at a given $\boldsymbol{k}$ . This broadening is directly related to the nonlinear turnover time scale of the fluctuations $\tau _{\mathrm{nl}}$ , which itself is of the order of the lifetime of the eddies before they break apart: the nonlinear interactions generating the cascade are precisely what cause the fluctuations to decorrelate over this time. For general levels of imbalance, the broadening in fluctuations at a given $\boldsymbol{k}$ must then be $\Delta \omega \sim \omega _{\mathrm{nl}}=\tau _{\mathrm{nl}}^{-1}$ .

The dependence of the behaviour of the wavevector–frequency spectrum on the level of imbalance reveals a deep connection to the mechanism by which ions are heated. In imbalanced turbulence, the linear frequencies $\omega _{\mathrm{A}}\sim v_{\mathrm{A}} / \ell _\|^+ \sim k_\perp \delta z^+$ of $\boldsymbol{z}^+$ fluctuations far exceed their nonlinear rate $\omega _{\mathrm{nl}}\sim k_\perp \delta z^-$ (Lithwick et al. Reference Lithwick, Goldreich and Sridhar2007; Schekochihin Reference Schekochihin2022), leading the wavevector–frequency spectrum to become sharply peaked around the dispersion relation of the fluctuations (because $\Delta \omega \sim \omega _{\mathrm{nl}}\ll \omega _{\mathrm{A}}$ ). This is the cyclotron-resonant-heating limit, where ions are able to resonate with exactly one frequency at a given $\boldsymbol{k}$ . In contrast, using critical balance we have $\omega _{\mathrm{A}}\sim v_{\mathrm{A}}/\ell _\| \sim \omega _{\mathrm{nl}}\sim k_\perp \delta z^\pm \sim \Delta \omega$ in balanced turbulence (Goldreich & Sridhar Reference Goldreich and Sridhar1995). This is the stochastic-heating limit, where fluctuations are sufficiently broadened and decorrelated allowing for ions to be heated by $\rho _{\mathrm{i}}$ -scale fluctuations. A smooth transition between these limits is expected as the imbalance is adjusted and one of the Elsasser fields dominates over the other; this is confirmed in Appendix A, where the wavevector–frequency spectra of simulations of RMHD turbulence with varying imbalance are investigated.

3.2.1. Two-dimensional wavevector energy spectrum

The two-dimensional wavevector energy spectrum in RMHD is defined as

(3.8) \begin{equation} \mathcal{E}_{\mathrm{2D}}(k_\perp ,k_\|) = 2\pi k_\perp \langle |\boldsymbol{z}^\pm (k_\perp , k_\|)|^2 \rangle , \end{equation}

where $\boldsymbol{z}^\pm (k_\perp , k_\|)$ is the complex Fourier amplitude of the $\boldsymbol{z}^\pm$ fields (which are assumed to be statistically cylindrically symmetric about ${\boldsymbol{B}}_0$ ) at a given $\boldsymbol{k}$ ; the $2\pi k_\perp$ term comes from the Jacobian of cylindrical coordinates in $\boldsymbol{k}$ -space. By analogy to the one-dimensional spectrum, where $\mathcal{E}(k){\mathrm{d}} k$ measures the energy contained in modes with $k\lt k'\lt k+{\mathrm{d}} k$ , this spectrum measures the energy in modes around a given parallel and perpendicular wavenumber. Via the Wiener–Khinchin theorem, (3.8) can likewise be interpreted as a measure of how correlated are turbulent fluctuations of a given scale along ( $\ell _\|$ ) and perpendicular ( $\lambda$ ) to the background magnetic field (Thorne & Blandford Reference Thorne and Blandford2017).Footnote ³

Schekochihin (Reference Schekochihin2022) argues that $\mathcal{E}_{\mathrm{2D}}(k_\perp ,k_\|)$ should be a product of power laws of $k_\perp$ and $k_\|$ :

(3.9) \begin{equation} \mathcal{E}_{\mathrm{2D}}(k_\perp ,k_\|)\sim \begin{cases} |k_\|L_\||^{-s_1}(k_\perp \lambda _{\mathrm{CB}})^{s_2}, & |k_\|L_\|| \gt (k_\perp \lambda _{\mathrm{ CB}})^{s_{\mathrm{CB}}},\\ |k_\|L_\||^{{s_4}}(k_\perp \lambda _{\mathrm{CB}})^{-{s_3}}, & |k_\|L_\|| \leqslant (k_\perp \lambda _{\mathrm{CB}})^{s_{\mathrm{CB}}}. \end{cases} \end{equation}

Here, $L_\|$ is the outer scale parallel to ${\boldsymbol{B}}_0$ and $\lambda _{\mathrm{CB}}$ is the perpendicular scale at which weak turbulence becomes strong (Schekochihin Reference Schekochihin2022). Due to critical balance (CB), we have $k_\| \sim k_\perp ^{s_{\mathrm{CB}}}$ . For the standard Goldreich & Sridhar (Reference Goldreich and Sridhar1995) CB argument covered in § 3.1, $s_{\mathrm{CB}} = 2/3$ , whereas for dynamic alignment, $s_{\mathrm{CB}} = 1/2$ (Boldyrev Reference Boldyrev2006).

The exponents in (3.9) can be determined based on the phenomenology of CB (Schekochihin Reference Schekochihin2022). At scales where $k_\parallel \leqslant k_\perp ^{s_{\mathrm{CB}}}$ , the spectrum measures correlations between points along magnetic-field lines that are separated by a distance $\ell _\|\gtrsim v_{\mathrm{A}}\tau _{\mathrm{A}}$ . Due to CB, these points are causally disconnected and thus uncorrelated; the correlation thus measures white noise and the spectrum must be flat along $k_\|$ (Thorne & Blandford Reference Thorne and Blandford2017), so ${s_4} = 0$ . The perpendicular exponent $s_3$ is fixed by requiring that the two-dimensional spectrum reproduces the one-dimensional perpendicular scaling implied by the constant-flux condition and CB (which imply $\delta z_\lambda ^2 \sim k_\perp ^{-s_{\mathrm{CB}}}$ ), yielding ${s_3} = 2s_{\mathrm{CB}} + 1$ . Continuity across the CB line then requires $s_1 = (2s_{\mathrm{CB}} + 1 + s_2)/s_{\mathrm{CB}}$ . Finally, $s_2$ depends on whether a kinematic ( $s_2 = 3$ ) or thermodynamic (equal energy in all modes; $s_2 = 1$ ) regime is assumed (Schekochihin Reference Schekochihin2022); in this work, we adopt the thermodynamic viewpoint, which better matches our simulation results in Appendix A.

Figure 1.

Two-dimensional spectrum model (3.13) with the Goldreich & Sridhar (Reference Goldreich and Sridhar1995) critical balance scaling $s_{\mathrm{CB}} = 2/3$ , $C=1$ and $\xi _{\rho ,{\mathrm{th}}}=0.5$ (with critical balance line $k_\| d_{\mathrm{ i}}=0.5(k_\perp \rho _{\mathrm{i}})^{2/3}$ ). The grey boundaries represent the spectrum cutoff at $|k_\| d_{\mathrm{i}}| \gt 1$ or $k_\perp \rho _{\mathrm{i}} \gt 1$ assumed in (3.13). Side panels show slices through the spectra (normalised to the maximum value along the slice), showing the individual scalings of $k_\perp$ and $k_\|$ above and below the CB line.

With this, (3.9) can be written entirely in terms of the CB exponent $s_{\mathrm{CB}}$ :

(3.10) \begin{equation} \mathcal{E}_{\mathrm{2D}}(k_\perp ,k_\|)\sim \begin{cases} |k_\|L_\||^{-(2s_{\mathrm{CB}}+2)/s_{\mathrm{CB}}}(k_\perp \lambda _{\mathrm{CB}}), & |k_\|L_\|| \gt (k_\perp \lambda _{\mathrm{ CB}})^{s_{\mathrm{CB}}},\\ (k_\|L_\|)^0(k_\perp \lambda _{\mathrm{CB}})^{-2s_{\mathrm{CB}} - 1}, & |k_\|L_\|| \leqslant (k_\perp \lambda _{\mathrm{CB}})^{s_{\mathrm{CB}}}. \end{cases} \end{equation}

For the Goldreich & Sridhar (Reference Goldreich and Sridhar1995) CB argument with $s_{\mathrm{CB}}=2/3$ , the spectrum follows $|k_\|L_\||^{-5}k_\perp \lambda _{\mathrm{CB}}$ for $|k_\|L_\|| \gt (k_\perp \lambda _{\mathrm{CB}})^{2/3}$ and $(k_\|L_\|)^0(k_\perp \lambda _{\mathrm{CB}})^{-7/3}$ for $|k_\|L_\|| \leqslant (k_\perp \lambda _{\mathrm{CB}})^{2/3}$ (see figure 1).

The derivation of (3.10) assumes that the turbulence is balanced, and is identical for both the $\boldsymbol{z}^+$ and $\boldsymbol{z}^-$ fields. In Appendix A, we argue that this model can be extended to describe imbalanced turbulence by using RMHD turbulence simulations with different values of $\sigma _{\mathrm{c}}$ and comparing their spectra with (3.10). There, we show that the model of Schekochihin (Reference Schekochihin2022) works well for describing the two-dimensional energy spectra at all levels of imbalance of the $\boldsymbol{z}^+$ and $\boldsymbol{z}^-$ fluctuations (with discrepancies above the CB cone for the $\boldsymbol{z}^-$ spectra), with the measured scalings showing $s_2\approx 1$ and $1/2 \lesssim s_{\mathrm{CB}}\lesssim 2/3$ .

Normalisation – Assuming a collection of ions with thermal speed $v_{\mathrm{th}, \textrm{i}}$ , thermal gyroradius $\rho _{\mathrm{i}}$ , and $\beta _{\mathrm{i}} = v_{\mathrm{th}, \textrm{i}}^2/v_{\mathrm{A}}^2$ , we normalise (3.10) such that

(3.11) \begin{equation} \frac {1}{v_{\mathrm{th}, \textrm{i}}^2}\int _{e^{-1/2}/\rho _{\mathrm{i}}}^{e^{1/2}/\rho _{\mathrm{i}}} {\mathrm{d}} k_\perp \ \int _{-k_\|^{\mathrm{ CB}}(k_\perp )}^{k_\|^{\mathrm{CB}}(k_\perp )} {\mathrm{d}} k_\| \ \mathcal{E}_{\mathrm{2D}}(k_\perp ,k_\|)\approx \xi _{\rho ,{\mathrm{th}}}^2, \end{equation}

where $k_\|^{\mathrm{CB}}(k_\perp ) \propto k_\perp ^{s_{\mathrm{CB}}}$ is the CB cone boundary. This normalisation ensures that the energy in $\rho _{\mathrm{i}}$ -scale modes is given by the standard stochastic heating parameter (Chandran et al. Reference Chandran, Li, Rogers, Quataert and Germaschewski2010a )

(3.12) \begin{equation} \xi _{\rho ,{\mathrm{th}}} =\frac {\delta u_{\rho _{\mathrm{i}}}}{v_{\mathrm{th}, \textrm{i}}}. \end{equation}

Normalising the wavenumbers as $\tilde {k}_\perp \equiv k_\perp \rho _{\mathrm{i}}$ and $\tilde {k}_\|\equiv k_\|v_{\mathrm{A}}/\varOmega _{\mathrm{i}} = k_\|d_{\textrm {i}}$ (where $d_{\textrm {i}}\equiv v_{\mathrm{A}}/\varOmega _{\mathrm{i}}$ is the ion inertial length), the spectrum becomes

(3.13) \begin{align} \mathcal{E}_{\mathrm{2D}}(k_\perp ,k_\|)&=\frac {1}{2}Vv_{\mathrm{A}}^2\rho _{\mathrm{i}}^2\beta _{\mathrm{i}}^{1/2}\begin{cases} C^{-1}\xi _{\rho ,{\mathrm{th}}}\tilde {k}_\perp ^{-{s_3}}, & |\tilde {k}_\||\leqslant \tilde {k}_\|^{\mathrm{CB}},\nonumber \\ C^{s_1-1}\xi _{\rho ,{\mathrm{th}}}^{s_1+1}|\tilde {k}_\||^{-s_1}\tilde {k}_\perp ^{s_2}, & |\tilde {k}_\||\geqslant \tilde {k}_\|^{\mathrm{CB}},\\ 0, & \tilde {k}_\perp \gt 1 \text{ or } |\tilde {k}_\|| \gt 1, \end{cases}\nonumber \\ &\equiv \frac {1}{2}Vv_{\mathrm{A}}^2\rho _{\mathrm{i}}^2\beta _{\mathrm{i}}^{1/2}\tilde {\mathcal{E}}_{\mathrm{2D}}(k_\perp ,k_\|), \end{align}

where $\tilde {\mathcal{E}}_{\mathrm{2D}}(k_\perp ,k_\|)$ is the dimensionless part of the spectrum, $V$ is the volume of the plasma, $C$ is an arbitrary Kolmogorov-like constant to account for order-unity factors, and the exponents are ${s_3} = 2s_{\mathrm{CB}}+1$ and $s_1 = (2s_{\mathrm{CB}}+s_2+1)/s_{\mathrm{CB}}$ . With this normalisation, the critical balance line scales as

(3.14) \begin{equation} k_\|v_{\mathrm{A}}\sim k_\perp \delta z_\lambda \implies \tilde {k}_\|^{\mathrm{CB}} \equiv C \xi _{\rho ,{\mathrm{th}}} \tilde {k}_\perp ^{s_{\mathrm{CB}}}. \end{equation}

For the remainder of this paper, we will assume the Goldreich & Sridhar (Reference Goldreich and Sridhar1995) CB scaling $s_{\mathrm{CB}}=2/3$ and $s_2=1$ (which follows from assuming the spectrum above the CB cone has reached a thermodynamic limit, as discussed already).

This spectrum is shown in figure 1 with $C=1$ and $\xi _{\rho ,{\mathrm{th}}} = 0.5$ . A common misunderstanding of the CB argument is that energy must be concentrated at the CB boundary $k_\|\sim k_\perp ^{2/3}$ . This is not true, with energy contained in all modes below the CB cone, as is clearly seen in figure 1. The steep drop in the spectrum in $k_\|$ for fluctuations above the CB cone means that there is very little energy in these fluctuations.

The cutoff at $\tilde {k}_\perp \gt 1$ stems from the assumption that $k_\perp \rho _{\mathrm{i}}\gtrsim 1$ fluctuations do not contribute strongly to heating, because they are suppressed by particles sampling many small-scale eddies over an orbit (Chandran Reference Chandran2000). However, these fluctuations may still have some effect (Arzamasskiy et al. Reference Arzamasskiy, Kunz, Chandran and Quataert2019; Isenberg & Vasquez Reference Isenberg and Vasquez2019); we quantify this contribution in Appendix C, where we use a model spectrum of strong balanced CB turbulence that includes $k_\perp \rho _{\mathrm{i}}\geqslant 1$ kinetic-Alfvén wave fluctuations, arguing their effect is often small when $k_\|/k_\perp \ll 1$ at $k_\perp \rho _{\mathrm{i}} \sim 1$ . Additionally, by including a cutoff at $|\tilde {k}_\|| \gt 1$ , this model also neglects the transition of Alfvén waves to dispersive ion-cyclotron waves (ICWs) as $k_\|$ approaches $d^{-1}_{\mathrm{i}}$ . Hybrid-kinetic simulations show the details of ICWs near $k_\| d_{\mathrm{i}}\sim 1$ are important for some aspects of ion heating (Squire et al. Reference Squire, Meyrand, Kunz, Arzamasskiy, Schekochihin and Quataert2022; Squire, Meyrand & Kunz Reference Squire, Meyrand and Kunz2023; Zhang et al. Reference Zhang, Kunz, Squire and Klein2025); additionally, due to the difference in their dispersion relations, the parallel resonant wavenumber of the Alfvén waves and ICWs differ, which may bias the heating rate (details on this discrepancy are presented in Appendix E). Regardless, we neglect the effect of the Alfvén-wave dispersion (i.e. the change to ICWs at $k_\| d_{\mathrm{i}} \sim 1$ ) on ion heating in our model, and assume the heating is dominated by large-scale Alfvénic fluctuations.

3.2.2. Temporal correlation function

For a turbulent system with nonlinear frequency $\omega _{\mathrm{nl}}(\boldsymbol{k})$ , a general temporal correlation function $f[\tau \omega _{\mathrm{nl}}(\boldsymbol{k})]$ encodes how correlated fluctuations are expected to be after a normalised lag time $\tau \omega _{\mathrm{nl}}(\boldsymbol{k})$ . This function should satisfy the properties $f(0)=1$ , $\lim _{\tau \to \infty }f[\tau \omega _{\mathrm{nl}}(\boldsymbol{k})]=0$ and $\int _{0}^{\infty } {\mathrm{d}}\tau f[\tau \omega _{\mathrm{nl}}(\boldsymbol{k})] \approx \omega ^{-1}_{\mathrm{nl}}(\boldsymbol{k})$ . For our model of RMHD turbulence with variable imbalance, we take the nonlinear frequency to be that of the dominant $\boldsymbol{z}^+$ fluctuations:

(3.15) \begin{equation} \omega _{\mathrm{nl}}(k_\perp ) = k_\perp \delta z^-_\lambda = C\varOmega _{\mathrm{i}}\alpha ^{-1/2}\xi _{\rho ,{\mathrm{th}}}(k_\perp \rho _{\mathrm{i}})^{s_{\mathrm{CB}}}, \end{equation}

where we have assumed the CB scaling $\omega _{\mathrm{nl}}\sim \omega _{\mathrm{A}} = k_\|v_{\mathrm{A}}$ and used the normalisation for the above-mentioned spectrum. To capture the reduced efficiency of the $\boldsymbol{z}^+$ cascade and the corresponding decrease in their nonlinear frequency as imbalance increases, $\omega _{\mathrm{nl}}$ is scaled by $\alpha ^{-1/2}$ ; this scaling arises from a simple estimate using CB for the dominant $\boldsymbol{z}^+$ fluctuations:

(3.16) \begin{equation} \omega _{\mathrm{nl}} = k_\perp \delta z^-_\lambda \sim \frac {\delta z^-_\lambda }{\delta z^+_\lambda }k_\perp \delta z^+_\lambda \propto \alpha ^{-1/2}. \end{equation}

A commonly used choice for the temporal correlation function that satisfies the above-mentioned properties is (e.g. Schlickeiser & Achatz Reference Schlickeiser and Achatz1993; Chandran Reference Chandran2000)

(3.17) \begin{equation} f[\tau \omega _{\mathrm{nl}}(k_\perp )]=e^{-|\tau \omega _{\mathrm{nl}}(k_\perp )|}; \end{equation}

however, the discontinuity in slope of the function at $\tau =0$ causes an unphysically slow drop-off in the Fourier transform of $f$ at large frequencies, which can affect the calculated heating rates (see Appendix D for more details). In contrast, a formal renormalisation group calculation in hydrodynamic turbulence shows that $f(\tau \omega _{\mathrm{nl}})$ reduces to (3.17) for $\tau \omega _{\mathrm{nl}} \gg 1$ and to a Gaussian $e^{-(\tau \omega _{\mathrm{nl}})^2}$ for $\tau \omega _{\mathrm{nl}} \ll 1$ (Gorbunova et al. Reference Gorbunova, Balarac, Canet, Eyink and Rossetto2021), the latter of which is smooth at $\tau =0$ . Based on this, a convenient choice for our model spectrum is

(3.18) \begin{equation} f[\tau \omega _{\mathrm{nl}}(k_\perp )]={\mathrm{sech}}[\tau \omega _{\mathrm{nl}}(k_\perp )], \end{equation}

which has the asymptotic limits

(3.19) \begin{equation} {\mathrm{sech}}(x)\approx \begin{cases} e^{-x^2/2}, & |x|\ll 1,\\ 2 e^{-|x|}, & |x|\gg 1, \end{cases} \end{equation}

and satisfies $\int _{0}^{\infty } {\mathrm{d}} \tau \ f(\tau \omega _{\mathrm{nl}}) = ({\pi }/{2})\omega _{\mathrm{nl}}^{-1}$ . Taking the Fourier transform in time of (3.18) with the convention

(3.20) \begin{equation} F(\omega ) = \mathcal{F}\left [f(\tau )\right ]= \frac {1}{2\pi }\int _{-\infty }^{\infty }{\mathrm{d}}\tau e^{{\mathrm{i}}\omega \tau }f(\tau ) \end{equation}

gives

(3.21) \begin{equation} F(\omega ;k_\perp ,\sigma _{\mathrm{c}}) = \frac {1}{2\omega _{\mathrm{nl}}(k_\perp )}{\mathrm{sech}}\left (\frac {\pi \omega }{2\omega _{\mathrm{nl}}(k_\perp )}\right)\!. \end{equation}

The choice of function used for $f(\tau \omega _{\mathrm{nl}})$ can affect the results of the following sections; however, qualitatively similar results are obtained so long as the tails of the function’s temporal Fourier transform decay sufficiently fast for large frequencies. The effect of the choice on heating is studied in Appendix D.

3.2.3. Full model

With the mentioned choices for the energy spectrum and temporal correlation of the turbulence, the following functional form is set forward as a simple model to capture the general dependence of RMHD turbulence on imbalance. DefiningFootnote ⁴

(3.22) \begin{equation} \mathcal{E}^\pm \left (k_\perp ,k_\|,\omega ;\sigma _{\mathrm{c}}\right )\equiv \frac {\mathcal{E}_{\mathrm{ 2D}}(k_\perp ,k_\|)}{2\pi k_\perp v_{\mathrm{A}}^2}F\left (\omega \pm k_\|v_{\mathrm{A}};k_\perp ,\sigma _{\mathrm{c}}\right )\!, \end{equation}

where $F(\omega )$ is the temporal Fourier transform of the correlation function $f(\tau )$ , the total wavevector–frequency spectrum is

(3.23) \begin{align} \mathcal{E}_{\mathrm{tot}}\left (k_\perp , k_\|, \omega ;\sigma _{\mathrm{c}}\right ) &\equiv \frac {1}{2}\left [(1+\sigma _{\mathrm{c}})\mathcal{E}^+\left (k_\perp ,k_\|,\omega ;\sigma _{\mathrm{c}}\right ) +(1-\sigma _{\mathrm{c}})\mathcal{E}^-\left (k_\perp ,k_\|,\omega ;\sigma _{\mathrm{c}}\right )\right ]\nonumber \\ &= \frac {V\rho _{\mathrm{i}}^2\beta _{\mathrm{i}}^{1/2}}{4}\frac {\tilde {\mathcal{E}}_{\mathrm{2D}}(k_\perp ,k_\|)}{2\pi k_\perp } \sum _{\nu =\pm 1}(1+\nu \sigma _{\mathrm{c}})F(\omega +\nu k_\|v_{\mathrm{A}};k_\perp ,\sigma _{\mathrm{c}}), \end{align}

where (3.13) has been used for $\mathcal{E}_{\mathrm{2D}}(k_\perp ,k_\|)$ . The amplitudes of the $\boldsymbol{z}^+$ and $\boldsymbol{z}^-$ spectra are weighted such that their ratio is $\alpha$ .

Figure 2.

Slices at constant $k_\perp L_\perp = 10$ through the model wavevector–frequency spectrum of RMHD turbulence (3.23, top), showing it qualitatively reproduces features of the RMHD simulations presented in Appendix A (bottom). The turbulence has imbalance $\sigma _{\mathrm{c}} = 0,\ 0.59$ and $0.96$ (left, middle and right columns, respectively), and the model sets $\xi _{\rho ,{\mathrm{th}}} = 1$ . The solid black lines correspond to zero frequency and the Alfvén dispersion relation $\omega _{\mathrm{A}} = \pm k_zv_{\mathrm{A}}$ ; the red dashed line is the critical balance scaling $k_\|^{\mathrm{CB}}$ , such that everything with $k_\| \lt k_\|^{\mathrm{CB}}$ is within the CB cone. The frequencies are normalised to the outer-scale Alfvén frequency of the simulations, $\omega _{\mathrm{A0}}=v_{\mathrm{A}}/L_z$ . Note that the appearance of energy in high frequencies at small $k_z$ in the bottom left and right figures is an artefact of the method used to calculate the spectrum.

Figure 3.

Slices at constant $k_z L_z = 10$ through the model wavevector–frequency spectrum of RMHD turbulence (3.23, top), showing it qualitatively reproduces features of the RMHD simulations presented in Appendix A (bottom). The turbulence has imbalance $\sigma _{\mathrm{c}} = 0,\ 0.59$ and $0.96$ (left, middle and right columns, respectively), and the model sets $\xi _{\rho ,{\mathrm{th}}} = 1$ . Solid black lines represent the centre of the Alfvén dispersion relation $\omega _{\mathrm{A}}=\pm k_zv_{\mathrm{A}}$ at the given value of $k_z$ , and the red dashed line shows where $k_\perp = (k_\|^{\mathrm{CB}})^{3/2}$ . Frequencies are normalised to the outer-scale Alfvén frequency, $\omega _{\mathrm{A0}}=v_{\mathrm{A}}/L_z$ .

Figures 2 and 3 compare the model in (3.23) to the wavevector–frequency spectra of the simulations in Appendix A for $\sigma _{\mathrm{c}}=0$ , 0.59 and 0.96. Slices are taken through the spectra at $k_\perp L_\perp =10$ and $k_z L_z=10$ , respectively, with $s_{\mathrm{CB}}=2/3$ and $\xi _{\rho ,{\mathrm{th}}}=1$ . For $k_\| L_z \gt (k_\perp L_\perp )^{2/3}$ (to the right or left of the red dashed line in figures 2 or 3, respectively), the spectrum consists of two bands of fluctuations centred on the $\boldsymbol{z}^\pm$ dispersion relations $\omega _{\mathrm{A}}^\pm = \mp k_\|v_{\mathrm{A}}$ , corresponding to the wave-like nature of turbulence in this region. For modes below the CB cone with $k_\| L_z \lesssim (k_\perp L_\perp )^{2/3}$ (to the left or right of the red dashed line in figures 2 or 3, respectively) the turbulence becomes strong, with nonlinear interactions causing the spectrum to broaden; the width of this broadening increases with $k_\perp$ (figure 3) as the nonlinear frequency $\omega _{\mathrm{nl}}(k_\perp )$ is an increasing function of $k_\perp$ . Additionally, the width of the spectrum decreases as the imbalance increases, corresponding to the reduced efficiency of the cascade of $\boldsymbol{z}^+$ fluctuations as the $\boldsymbol{z}^-$ fluctuations become subdominant. In the limit of fully imbalanced turbulence ( $\sigma _{\mathrm{c}}\to 1$ ), the spectrum would become a delta function along the $\boldsymbol{z}^+$ dispersion relation $\omega _{\mathrm{A}} = -k_\|v_{\mathrm{A}}$ .

Figures 2 and 3 show that the model (3.23) qualitatively captures the features of the wavevector–frequency spectra of RMHD turbulence, as well as its dependence on the turbulence’s imbalance. The ability to express the spectrum in this form makes it useful for connecting theoretical models of RMHD turbulence to numerical or observational data across a range of plasma conditions.

4.1.1. Perpendicular heating rate of a Maxwellian distribution

The velocity-space diffusion of the distribution function results in particle heating. In this section, we derive the form of the perpendicular heating rate $Q_\perp$ arising from this diffusion. For simplicity, we assume the initial ion distribution is Maxwellian,

(4.11) \begin{equation} f_0(\tilde{v}_\perp , \tilde{v}_\|)=n_0\left (\pi \beta _{\mathrm{i}}\right )^{-3/2}\exp\! {\left (-\frac {\tilde{v}_\perp ^2+\tilde{v}_\|^2}{\beta _{\mathrm{i}}}\right )}, \end{equation}

where $\tilde{v}_{\perp ,\|} \equiv v_{\perp ,\|}/v_{\mathrm{A}}$ and $f_0$ is normalised such that $\int _{-\infty }^{\infty } {\mathrm{d}} \tilde{v}_\| \ \int _{0}^{\infty } {\mathrm{d}} \tilde{v}_\perp \ 2\pi \tilde{v}_\perp f_0 = v_{\mathrm{A}}^{3}n_0$ .

The perpendicular heating rate per unit mass is given by

(4.12) \begin{equation} Q_\perp \equiv \frac {{\mathrm{d}}\langle v_\perp ^2 /2 \rangle }{{\mathrm{d}} t}= \int _{-\infty }^{\infty } {\mathrm{d}} v_\| \ \int _{0}^{\infty } {\mathrm{d}} v_\perp \ 2\pi v_\perp \left (\frac {v_\perp ^2}{2}\right ) \frac {1}{n_0}\frac {\partial f_0}{\partial t}. \end{equation}

Assuming that ions undergo a diffusion in velocity space in the quasi-linear form (2.5), (4.12) can be written as

(4.13) \begin{align} Q_\perp &= \int _{-\infty }^{\infty } {\mathrm{d}} v_\| \ \int _{0}^{\infty } {\mathrm{d}} v_\perp \ 2\pi v_\perp \frac {v_\perp ^2}{2}\frac {1}{v_\perp }\frac {\partial }{\partial v_\perp }\left [v_\perp \frac {1}{n_0}\left (D_{\perp \perp }\frac {\partial f_0}{\partial v_\perp }+D_{\|\perp }\frac {\partial f_0}{\partial v_\|}\right )\right ]\nonumber \\[3pt] &= -2\pi \int _{-\infty }^{\infty } {\mathrm{d}} v_\| \ \int _{0}^{\infty } {\mathrm{d}} v_\perp \ v_\perp ^2\frac {1}{n_0}\left (D_{\perp \perp }\frac {\partial f_0}{\partial v_\perp }+D_{\|\perp }\frac {\partial f_0}{\partial v_\|}\right )\!, \end{align}

by assuming that the velocity gradients of a physically reasonable distribution function should go to zero as $|\boldsymbol{v}|\to \infty$ to remove $\partial / \partial v_\|$ terms via integration by parts; the second equality is obtained by further integration by parts over $v_\perp$ . Changing variables to $\tilde{v}_{\perp ,\|}$ and using the derivatives of the Maxwellian, (4.13) becomes

(4.14) \begin{equation} Q_\perp = \frac {4}{\sqrt {\pi \beta _{\mathrm{i}}^5}}\int _{-\infty }^{\infty } {\mathrm{d}} \tilde{v}_\| \ e^{-\tilde{v}_\|^2/\beta _{\mathrm{i}}}\int _{0}^{\infty } {\mathrm{d}} \tilde{v}_\perp \ \tilde{v}_\perp ^2 e^{-\tilde{v}_\perp ^2/\beta _{\mathrm{i}}} \left (\tilde{v}_\perp D_{\perp \perp }+\tilde{v}_\|D_{\|\perp }\right)\!. \end{equation}

Finally, inserting the diffusion coefficients (4.4) and noting that $\mathcal{D}$ in (4.10) is independent of $\tilde{v}_\perp$ , we obtain

(4.15) \begin{equation} \frac {Q_\perp }{\varOmega _{\mathrm{i}}v_{\mathrm{th}, \textrm{i}}^2} = \sqrt {\frac {4}{\pi \beta _{\mathrm{i}}}}\int _{-\infty }^{\infty } {\mathrm{d}} \tilde{v}_\| \ e^{-\tilde{v}_\|^2/\beta _{\mathrm{i}}}(1+\sigma _{\mathrm{c}}\tilde{v}_\|)\frac {\mathcal{D}(\tilde{v}_\|)}{\varOmega _{\mathrm{i}}v_{\mathrm{th}, \textrm{i}}^2}. \end{equation}

This form will be used to calculate $Q_\perp$ analytically and numerically in §§ 4.2 and 4.3. We note that the form of $Q_\perp$ in (4.15) only gives the instantaneous heating rate of a Maxwellian distribution and does not describe the heating once the diffusion causes $f_0$ to evolve away from a Maxwellian. This evolution is shown in more detail in § 4.4, where we numerically solve the quasi-linear diffusion equation in (2.5).

4.2.1. Fully imbalanced turbulence: $\sigma _{\mathrm{c}}=1$

The fully imbalanced limit may be viewed as the case $\sigma _{\mathrm{c}} \to 1$ , corresponding to a collection of $\boldsymbol{z}^\pm$ fluctuations undergoing a cascade that is infinitely slow compared with their propagation rate. In this limit, the quasi-linear framework reduces to studying the effect of strong interactions between ions and waves following a linear dispersion relation $\omega (\boldsymbol{k})$ , with the evolution of the distribution function controlled by the resonance condition (2.6). When this condition is satisfied, the distribution function diffuses in velocity space along constant energy contours in the wave frame (Kennel & Engelmann Reference Kennel and Engelmann1966; we prove this explicitly in Appendix B.2.1). For the Alfvénic $\boldsymbol{z}^+$ fluctuations, with $\omega (\boldsymbol{k}) = -k_\|v_{\mathrm{A}}$ , the resonance condition becomes

(4.16) \begin{equation} k_\|(v_\|+v_{\mathrm{A}})+n\varOmega _{\mathrm{i}} = 0. \end{equation}

Taking the limit $\sigma _{\mathrm{c}}\to 1$ in the model of the wavevector–frequency spectrum in § 3.2.3, the nonlinear frequency $\tilde {\omega }_{\mathrm{nl}}(\tilde{k}_\perp )$ goes to zero ( $\alpha = (1+\sigma _{\mathrm{c}})/(1-\sigma _{\mathrm{c}}) \to \infty$ in 3.15) and the temporal correlation function approaches a delta function (Wheeler Reference Wheeler2002),

(4.17) \begin{equation} \frac {1}{2\tilde {\omega }_{\mathrm{nl}}(\tilde{k}_\perp )}{\mathrm{sech}}\left (\frac {\pi }{2\tilde {\omega }_{\mathrm{nl}}(\tilde{k}_\perp )}\left [\tilde{k}_\|\left (\tilde{v}_\|+\nu \right ) + 1\right ]\right )\to \delta \big(\tilde{k}_\|(\tilde{v}_\|+\nu) + 1\big), \end{equation}

which is exactly the resonance condition (4.16) with $n=1$ . With this, in the limit $\sigma _{\mathrm{c}}\to 1$ , the generalised diffusion coefficient in (4.10) becomes

(4.18) \begin{equation} \frac {\mathcal{D}}{\varOmega _{\mathrm{i}}v_{\mathrm{th}, \textrm{i}}^2} = \frac {\pi }{16}\int _{0}^{1} {\mathrm{d}} \tilde{k}_\perp \ \int _{-1}^{1} {\mathrm{d}} \tilde{k}_\| \ \tilde {\mathcal{E}}_{\mathrm{2D}}(\tilde{k}_\perp ,\tilde{k}_\|)\delta \big[\tilde{k}_\|(\tilde{v}_\| +1)+1\big]. \end{equation}

The delta function enforces that ions with parallel speed $v_\|$ interact with fluctuations with resonant wavenumber

(4.19) \begin{equation} \tilde{k}_\|^{(1)}\equiv -(1+\tilde{v}_\|)^{-1}. \end{equation}

Note that, because $|\tilde{k}_\||\leqslant 1$ in this model (as for $|\tilde{k}_\|| = k_\|d_{\textrm {i}}\gtrsim 1$ , the fluctuations are no longer Alfvénic, transitioning to ICWs), resonance only occurs for ions with $\tilde{v}_\|\gt 0$ , counter-propagating along magnetic-field lines with respect to the waves.

The resonance condition (4.19) shows that the speed of the ions controls the part of the turbulent spectrum with which they are able to interact. If ions are too slow,

(4.20) \begin{equation} \xi _{\rho ,{\mathrm{th}}} \leqslant |\tilde{k}_\|^{(1)}|\leqslant 1 \implies 0 \leqslant \tilde{v}_\|\leqslant \xi _{\rho ,{\mathrm{th}}}^{-1}-1, \end{equation}

they can only interact with the portion of the spectrum above the CB cone; however, if $|\tilde{k}_\|^{(1)}| \lt \xi _{\rho ,{\mathrm{th}}}$ , they are also able to interact with that below the cone. Performing the integrals, setting $\tilde{k}_\| = \tilde{k}_\|^{(1)}$ using the delta function and allowing for $\xi _{\rho ,{\mathrm{th}}}\gtrsim 1$ , we obtain (setting $C=1$ in $\tilde {\mathcal{E}}_{\mathrm{2D}}$ , 3.13)

(4.21) \begin{equation} \frac {\mathcal{D}}{\varOmega _{\mathrm{i}}v_{\mathrm{th}, \textrm{i}}^2}=\frac {\pi }{16(s_2+1)}\begin{cases} 0, & \tilde{v}_\| \lt 0,\\ \xi _{\rho ,{\mathrm{th}}}^{s_1+1}(1+\tilde{v}_\|)^{s_1}, & 0\leqslant \tilde{v}_\|\leqslant \max (0,\xi _{\rho ,{\mathrm{th}}}^{-1}-1),\\[3pt] \dfrac {1}{2}\xi _{\rho ,{\mathrm{th}}}^3\left [s_1(1+\tilde{v}_\|)^2-\dfrac {s_2+1}{s_{\mathrm{CB}}}\xi _{\rho ,{\mathrm{th}}}^{-2}\right ]\!, & \tilde{v}_\| \gt \max (0,\xi _{\rho ,{\mathrm{th}}}^{-1}-1); \end{cases} \end{equation}

using the scalings for the Goldreich & Sridhar (Reference Goldreich and Sridhar1995) spectrum ( $s_{\mathrm{CB}} = 2/3$ and $s_2=1$ ) gives

(4.22) \begin{equation} \frac {\mathcal{D}}{\varOmega _{\mathrm{i}}v_{\mathrm{th}, \textrm{i}}^2}=\frac {\pi }{32}\begin{cases} 0, & \tilde{v}_\| \lt 0,\\ \xi _{\rho ,{\mathrm{th}}}^6(1+\tilde{v}_\|)^5, & 0\leqslant \tilde{v}_\|\leqslant \max (0,\xi _{\rho ,{\mathrm{th}}}^{-1}-1),\\[3pt] \dfrac {1}{2}\xi _{\rho ,{\mathrm{th}}}^3\left [5(1+\tilde{v}_\|)^2-3\xi _{\rho ,{\mathrm{th}}}^{-2}\right ]\!, & \tilde{v}_\| \gt \max (0,\xi _{\rho ,{\mathrm{th}}}^{-1}-1).\\ \end{cases} \end{equation}

This form of $\mathcal{D}$ transitions smoothly from a $\xi _{\rho ,{\mathrm{th}}}^3$ to a $\xi _{\rho ,{\mathrm{th}}}^6$ scaling with increasing $\tilde{v}_\|$ .Footnote ⁶

Figure 4.

Perpendicular heating rate $Q_\perp$ in the fully imbalanced limit, where the general quasi-linear theory of velocity-space diffusion reduces to diffusion along contours of constant energy in the wave frame (Kennel & Engelmann Reference Kennel and Engelmann1966). $Q_\perp$ is calculated numerically from (4.23) using the imbalanced limit of the generalised diffusion coefficient $\mathcal{D}$ (4.22) for different values of $\beta _{\mathrm{i}}$ (with $C=1$ in $\tilde {\mathcal{E}}_{\mathrm{2D}}$ , 3.13). Dotted lines show only the contribution of fluctuations below the CB cone to $Q_\perp$ .

In the limit $\sigma _{\mathrm{c}}\to 1$ , the general heating rate in (4.15) reduces to

(4.23) \begin{equation} \frac {Q_\perp }{\varOmega _{\mathrm{i}}v_{\mathrm{th}, \textrm{i}}^2} = \sqrt {\frac {4}{\pi \beta _{\mathrm{i}}}}\int _{-\infty }^{\infty } {\mathrm{d}} \tilde{v}_\| \ e^{-\tilde{v}_\|^2/\beta _{\mathrm{i}}}(1+\tilde{v}_\|)\frac {\mathcal{D}(\tilde{v}_\|)}{\varOmega _{\mathrm{i}}v_{\mathrm{th}, \textrm{i}}^2}. \end{equation}

Figure 4 shows $Q_\perp$ as a function of $\xi _{\rho ,{\mathrm{th}}}$ and $\beta _{\mathrm{i}}$ , calculated numerically using the integral in (4.23). The solid lines are calculated using the form of $\mathcal{D}$ in (4.22), which shows the contribution of modes in $\boldsymbol{k}$ -space both above and below the CB cone to the heating rate. These inherit the scaling with $\xi _{\rho ,{\mathrm{th}}}$ present in $\mathcal{D}$ , transitioning from a $\xi _{\rho ,{\mathrm{th}}}^6$ scaling to a $\xi _{\rho ,{\mathrm{th}}}^3$ scaling with increasing turbulent amplitude. The presence of the $\xi _{\rho ,{\mathrm{th}}}^6$ scaling at small $\xi _{\rho ,{\mathrm{th}}}$ (a feature also seen in the balanced limit and general calculation later) suggests that ions can be heated by weak fluctuations above the CB cone (although the assumptions used in deriving this model may mean this is not quite correct; see § 5 for a discussion of the caveats of this model). The heating rates in figure 4 show little to no variation with $\beta _{\mathrm{i}}$ for $\beta _{\mathrm{i}}\ll 1$ . This is because only ions with $\tilde{v}_\|\ll 1$ contribute to the integral; ions with large $\tilde{v}_\|$ will contribute less to the total heating rate $Q_\perp$ as there are comparatively fewer of them.

If we consider ions only interacting with strong fluctuations below the CB cone (dotted lines in figure 4), we instead see that the heating rate is strongly suppressed at small $\xi _{\rho ,{\mathrm{th}}}$ . This follows from an analogous calculation for $\mathcal{D}$ in (4.22) where only the contribution of modes below the CB cone ( $|\tilde{k}_\|| \leqslant \tilde{k}_\|^{\mathrm{CB}}$ ) are taken into account, which gives

(4.24) \begin{equation} \frac {\mathcal{D}}{\varOmega _{\mathrm{i}}v_{\mathrm{th}, \textrm{i}}^2} = \frac {3\pi }{64}\xi _{\rho ,{\mathrm{th}}}^3\left [(1+\tilde{v}_\|)^2 - \xi _{\rho ,{\mathrm{th}}}^{-2}\right ] \end{equation}

for $\tilde{v}_\| \gt \max (0,\xi _{\rho ,{\mathrm{th}}}^{-1}-1)$ and $0$ otherwise. For ions that are fast enough to interact with modes below the CB cone, $\mathcal{D}$ scales as $\xi _{\rho ,{\mathrm{th}}}^3$ and drops off rapidly as $\xi _{\rho ,{\mathrm{th}}}$ approaches $1/(1+\tilde{v}_\|)$ . This is quite restrictive compared with the $\xi _{\rho ,{\mathrm{th}}}\gtrsim 0.1$ requirement of stochastic heating for $\tilde{v}_\|\ll 1$ , although the exact threshold depends on the choice of $C$ . It is a result of the resonance condition (4.19) allowing ions to interact with a single frequency.

4.2.2. Balanced turbulence: $\sigma _{\mathrm{c}}=0$

In contrast to the imbalanced limit, balanced turbulence consists of a collection of $\boldsymbol{z}^+$ and $\boldsymbol{z}^-$ fluctuations with equal energies. The nonlinear interactions between the two populations lead to broadening of the wavevector–frequency spectrum, which is captured by the temporal correlation function (3.18). Using (4.10) for the generalised diffusion coefficient with $\sigma _{\mathrm{c}}=0$ (where $\alpha =1$ in $\tilde {\omega }_{\mathrm{nl}}$ ), we have

(4.25) \begin{align} \frac {\mathcal{D}}{\varOmega _{\mathrm{i}}v_{\mathrm{th}, \textrm{i}}^2} &= \frac {1}{16}\sum _{\nu =\pm 1} \nonumber \\ &\quad \times\int _{0}^{1} {\mathrm{d}} \tilde{k}_\perp \ \frac {\pi }{2\tilde {\omega }_{\mathrm{nl}}(\tilde{k}_\perp )}\int _{-1}^{1} {\mathrm{d}} \tilde{k}_\| \ \tilde {\mathcal{E}}_{\mathrm{2D}}(\tilde{k}_\perp ,\tilde{k}_\|){\mathrm{sech}}\!\left (\frac {\pi }{2\tilde {\omega }_{\mathrm{nl}}(\tilde{k}_\perp )}\big[\tilde{k}_\|(\tilde{v}_\|+\nu)+1\big]\right)\!. \end{align}

The integral (4.25) is now symmetric in $\tilde{v}_\|$ because there are fluctuations travelling in both directions along the magnetic field. The general heating rate (4.15) in this limit reduces to

(4.26) \begin{equation} \frac {Q_\perp }{\varOmega _{\mathrm{i}}v_{\mathrm{th}, \textrm{i}}^2} = \sqrt {\frac {4}{\pi \beta _{\mathrm{i}}}}\int _{-\infty }^{\infty } {\mathrm{d}} \tilde{v}_\| \ e^{-\tilde{v}_\|^2/\beta _{\mathrm{i}}}\frac {\mathcal{D}(\tilde{v}_\|)}{\varOmega _{\mathrm{i}}v_{\mathrm{th}, \textrm{i}}^2}. \end{equation}

The full integral of (4.25) is not able to be completed analytically over all $\boldsymbol{k}$ -space and can only be calculated numerically. To investigate its limiting behaviour, we first use the assumption of low- $\beta _{\mathrm{i}}$ ions such that $\tilde{v}_\| = v_\|/v_{\mathrm{A}}\ll 1$ , removing the dependence of $\mathcal{D}$ on $\tilde{v}_\|$ in (4.25). This approximation essentially ignores the equivalent resonant condition in the imbalanced limit, where ions are only able to interact with a single frequency at a given $\boldsymbol{k}$ (see 4.19). However, because the turbulent fluctuations in the balanced case have a greater spread in frequency due to their nonlinear broadening, ions are still able to diffuse in velocity space as they interact with a variety of fluctuations with differing frequencies at a given $\boldsymbol{k}$ .

With the $\tilde{v}_\|\ll 1$ limit, the contributions of the integral (4.25) to ion heating from below the CB cone can be investigated by looking at fluctuations such that the ratio

(4.27) \begin{equation} \frac {\tilde{k}_\|}{\tilde {\omega }_{\mathrm{nl}}(\tilde{k}_\perp )} = \frac {\tilde{k}_\|}{C\xi _{\rho ,{\mathrm{th}}} \tilde{k}_\perp ^{s_{\mathrm{CB}}}}\lesssim 1. \end{equation}

In this region, $\tilde {\mathcal{E}}_{\mathrm{2D}} = C^{-1}\xi _{\rho ,{\mathrm{th}}}\tilde{k}_\perp ^{-2s_{\mathrm{CB}}-1}$ and, after taking the $x\gg 1$ limit of ${\mathrm{sech}}(x)$ in (3.19) as $1/ \tilde {\omega }_{\mathrm{nl}} = \varOmega _{\mathrm{i}}/ \omega _{\mathrm{nl}} \gtrsim 1$ , (4.25) becomesFootnote ⁷

(4.28) \begin{align} \frac {\mathcal{D}}{\varOmega _{\mathrm{i}}v_{\mathrm{th}, \textrm{i}}^2} &\approx \frac {1}{4}\int _{0}^{1} {\mathrm{d}} \tilde{k}_\perp \ C^{-1}\xi _{\rho ,{\mathrm{th}}}\tilde{k}_\perp ^{-2s_{\mathrm{CB}}-1}\frac {\pi }{2\tilde {\omega }_{\mathrm{nl}}(\tilde{k}_\perp )}\exp\!\left (-\frac {\pi }{2\tilde {\omega }_{\mathrm{nl}}(\tilde{k}_\perp )}\right )\int _{-C\xi _{\rho ,{\mathrm{th}}}\tilde{k}_\perp ^{s_{\mathrm{CB}}}}^{C\xi _{\rho ,{\mathrm{th}}}\tilde{k}_\perp ^{s_{\mathrm{CB}}}} {\mathrm{d}} \tilde{k}_\| \ \nonumber \\ &= \frac {C}{\pi s_{\mathrm{CB}}}\xi _{\rho ,{\mathrm{th}}}^3\int _{\pi /(2C\xi _{\rho ,{\mathrm{th}}})}^{\infty } {\mathrm{d}} u \ u e^{-u}\nonumber \\ &= \frac {C}{\pi s_{\mathrm{CB}}}\xi _{\rho ,{\mathrm{th}}}^3 \left (1 + \frac {\pi }{2C\xi _{\rho ,{\mathrm{th}}}}\right )e^{-\pi / (2C\xi _{\rho ,{\mathrm{th}}})}, \end{align}

where a change of variables to $u = \pi / (2\tilde {\omega }_{\mathrm{nl}}(\tilde{k}_\perp )) = (\pi /2)(C\xi _{\rho ,{\mathrm{th}}})^{-1}\tilde{k}_\perp ^{-s_{\mathrm{CB}}}$ was made in the second equality. Finally, as $\mathcal{D}$ is independent of $\tilde{v}_\|$ for low- $\beta _{\mathrm{i}}$ ions, the heating rate (4.26) becomes

(4.29) \begin{equation} \frac {Q_\perp }{\varOmega _{\mathrm{i}}v_{\mathrm{th}, \textrm{i}}^2} = \hat {c}_1\xi _{\rho ,{\mathrm{th}}}^3 \left (1 + \frac {\hat {c}_2}{\xi _{\rho ,{\mathrm{th}}}}\right )e^{-\hat {c}_2/ \xi _{\rho ,{\mathrm{th}}}}, \end{equation}

where $\hat {c}_1 \equiv 2C/(\pi s_{\mathrm{CB}}) \approx C$ (for the Goldreich & Sridhar (Reference Goldreich and Sridhar1995) CB scaling with $s_{\mathrm{CB}}=2/3$ ) and $\hat {c}_2 \equiv \pi / (2C) \approx 1.5C^{-1}$ .

Figure 5.

Perpendicular heating rate $Q_\perp$ in the balanced limit for different values of $\beta _{\mathrm{i}}$ , calculated numerically using (4.26) with $C=5$ in $\tilde {\mathcal{E}}_{\mathrm{2D}}$ (3.13). To better highlight the suppression, only the contribution from modes below the CB cone are considered. The black dotted line shows the $\beta _{\mathrm{i}}$ -independent analytic expression for $Q_\perp$ (4.29) (obtained in the limit $|\tilde{v}_\|| \ll 1$ and $\tilde{k}_\|\ll \tilde{k}_\|^{\mathrm{CB}}$ ) with $\hat {c}_1\approx 5$ and $\hat {c}_2 \approx 0.3$ . These show qualitatively similar behaviour to the empirical stochastic heating formula (1.3) (Chandran et al. Reference Chandran, Li, Rogers, Quataert and Germaschewski2010a ), shown by the blue dashed line with the same values for $c_1$ and $c_2$ .

Remarkably, (4.29) has an almost identical form to the empirical stochastic heating formula (1.3) (Chandran et al. Reference Chandran, Li, Rogers, Quataert and Germaschewski2010a ). The coefficient $c_2$ in (1.3) is generally found to be $\lesssim 0.5$ in numerical simulations of stochastic heating (Chandran et al. Reference Chandran, Li, Rogers, Quataert and Germaschewski2010a ; Xia et al. Reference Xia, Perez, Chandran and Quataert2013; Cerri et al. Reference Cerri, Arzamasskiy and Kunz2021; Johnston et al. Reference Johnston, Squire and Meyrand2025); using this to constrain $C$ , we find $C\gtrsim 3$ and $\hat {c}_1\gtrsim 3$ . This may be larger than results reported in past simulations, although there is no agreement on exact values for $c_1$ ; for example, Xia et al. (Reference Xia, Perez, Chandran and Quataert2013) find $0.5\lesssim c_1 \lesssim 1.5$ in low- $\beta _{\mathrm{i}}$ balanced RMHD turbulence and $c_1\approx 4$ at $\beta _{\mathrm{i}}=1$ , but Johnston et al. (Reference Johnston, Squire and Meyrand2025) find $c_1\approx 4$ in balanced turbulence at $\beta _{\mathrm{i}} = 0.05$ . Given the drastic assumptions of quasi-linear theory and the likely dependence on other properties such as the choice of forcing or inertial-range length (for instance, as a consequence of intermittency; Mallet et al. Reference Mallet, Klein, Chandran, Grošelj, Hoppock, Bowen, Salem and Bale2019; Bowen et al. Reference Bowen, Ervin, Mallet, Chandran, Sioulas, Isenberg, Bale, Squire, Klein and Pezzi2025), the general agreement of the functional form and order-unity agreement of the coefficients is encouraging.

Figure 5 shows the integral (4.26) used to calculate $Q_\perp$ for different values of $\beta _{\mathrm{i}}$ . To compare with previous simulation results discussed in the previous paragraph, we determine $C$ from a representative value of the coefficient $\hat {c}_2\approx 0.3$ using the relation $\hat {c}_2\approx 1.5C^{-1}$ from (4.29), giving $C=5$ ; the value $\hat {c}_1\approx 5$ then follows from the same relations. Additionally, we only take into account the contribution from modes below the CB cone to better show the form of the heating-rate suppression. The calculated curves show little to no variation with $\beta _{\mathrm{i}}$ when $\beta _{\mathrm{i}}\ll 1$ , with small changes at $\xi _{\rho ,{\mathrm{th}}} \ll 1$ when $\beta _{\mathrm{i}}=1$ . This behaviour is exactly what is expected from stochastic heating, where $Q_\perp /(\varOmega _{\mathrm{i}}v_{\mathrm{th}, \textrm{i}}^2)$ is independent of $\beta _{\mathrm{i}}$ for $\beta _{\mathrm{i}}\ll 1$ and heating in small-amplitude turbulence is suppressed due to the conservation of the magnetic moment of the ions (Chandran et al. Reference Chandran, Li, Rogers, Quataert and Germaschewski2010a ) (however, due to the term $1+\hat {c}_2/\xi _{\rho ,{\mathrm{th}}}$ in (4.29), the suppression of $Q_\perp$ calculated by this model decreases more slowly than that of the empirical form for the choice of parameters used). If the contribution of modes above the CB cone are also taken into account, the model predicts the heating of ions by weak fluctuations as in the imbalanced case; this is explicitly shown by the $\xi _{\rho ,{\mathrm{th}}}^6$ scaling for small $\xi _{\rho ,{\mathrm{th}}}$ in figure 6.

Finally, we note that the suppression in (4.28) and (4.29) arises in part from the form of the temporal correlation function, rather than just from the exponential dependence of the Maxwellian distribution function assumed. This indicates that the effect is fundamentally determined by the dynamical properties of the turbulence and its coupling to the ions, as noted previously (Cerri et al. Reference Cerri, Arzamasskiy and Kunz2021; Mallet et al. Reference Mallet, Klein, Chandran, Ervin and Bowen2026). This perspective contrasts with the original interpretation in (1.3), where the suppression is attributed to the conservation of the magnetic moment of ions in the presence of small-amplitude fluctuations (Chandran et al. Reference Chandran, Li, Rogers, Quataert and Germaschewski2010a ).

4.3.1. Physical origin of the heating-rate scalings

The reason for this scaling behaviour can be understood by studying the properties of the integrand of the generalised diffusion coefficient $\mathcal{D}$ (4.10). This is shown in figure 7 for the same values of $\sigma _{\mathrm{c}}$ as in figure 6 with $\xi _{\rho ,{\mathrm{th}}}=0.5$ and $\tilde{v}_\|=0.1$ ; we show only $k_\| \lt 0$ modes, which correspond to the sign of the dispersion relation of the $\boldsymbol{z}^+$ fluctuations ( $-\omega _{\mathrm{A}}=-k_\|v_{\mathrm{A}}$ ). The integrand of $\mathcal{D}$ represents the contribution that fluctuations with a given $(k_\perp ,k_\|)$ provide to the heating of ions with parallel velocity $\tilde{v}_\|$ . The integrand consists of a peak at $|\tilde{k}_\|^{(1)}| = (1+\tilde{v}_\|)^{-1}$ (approximately $0.9$ for $\tilde{v}_\|=0.1$ , red dashed line in figure 7), the resonant parallel wavenumber in the imbalanced limit (4.19), as well as a region of $\boldsymbol{k}$ -space both above and below the CB cone that is dependent on the broadening of the wavevector–frequency spectrum arising from nonlinear interactions. Because these interactions become less dominant as $\sigma _{\mathrm{c}}$ is increased, the width of the broadening decreases and the integrand becomes increasingly sharply peaked around $\tilde{k}_\|^{(1)}$ .

Figure 7.

Integrand of the generalised diffusion coefficient $\mathcal{D}$ , (4.10), which therefore shows the contribution of fluctuations at a given $(\tilde{k}_\perp ,\tilde{k}_\|)$ to the overall heating. Each integrand is normalised to its maximum value and is plotted with $\xi _{\rho ,{\mathrm{th}}} = 0.5$ , $\tilde{v}_\|=0.1$ , and $\sigma _{\mathrm{c}} = (a)\, 0.0$ , (b) $0.75$ and (c) $0.999$ . The black line corresponds to the CB cone $\tilde{k}_\|^{\mathrm{CB}} = \xi _{\rho ,{\mathrm{th}}}\tilde{k}_\perp ^{2/3}$ . The horizontal red dashed line is $\tilde{k}_\|^{(1)}\approx 0.9$ , the resonant $k_\|$ with which ions interact in the imbalanced limit defined in (4.19). Note that we are looking at the $\tilde{k}_\| \lt 0$ modes.

The scaling of $Q_\perp$ and $\mathcal{D}$ with $\xi _{\rho ,{\mathrm{th}}}$ in figure 6 can be intuitively understood by noting that varying $\xi _{\rho ,{\mathrm{th}}}$ changes the position of the CB cone, $\tilde{k}_\|^{\mathrm{CB}} = \xi _{\rho ,{\mathrm{th}}}\tilde{k}_\perp ^{2/3}$ , in figure 7. As shown previously in the balanced limit (§ 4.2), the region below the CB cone contributes a term to $\mathcal{D}$ proportional to $\xi _{\rho ,{\mathrm{th}}}^3$ ; similarly, the region above the cone contributes to the $\xi _{\rho ,{\mathrm{th}}}^6$ scaling in both the $\sigma _{\mathrm{c}}=0$ and $1$ limits. As $\xi _{\rho ,{\mathrm{th}}}$ is increased, the cone shifts upwards in figure 7, causing fluctuations with $\tilde{k}_\|\lt \tilde{k}_\|^{\mathrm{CB}}$ to have a greater contribution to the integral, and leading to the $\xi _{\rho ,{\mathrm{th}}}^3$ scaling in (4.28) and (4.22). Similarly, the region above the cone becomes the dominant contribution as $\xi _{\rho ,{\mathrm{th}}}$ is decreased and the cone shifts downwards, leading to the $\xi _{\rho ,{\mathrm{th}}}^6$ scaling. Because the width of the broadening function decreases with increasing $\sigma _{\mathrm{c}}$ , a larger value of $\xi _{\rho ,{\mathrm{th}}}$ is required so that the region below the cone can contribute to the integral. This leads to the transition to a $\xi _{\rho ,{\mathrm{th}}}^3$ scaling at larger values of $\xi _{\rho ,{\mathrm{th}}}$ for highly imbalanced turbulence, as seen in figure 6.

Figure 7 also illustrates the suppression of $Q_\perp$ that arises when only modes below the CB cone are included. This suppression depends on the contribution from the tails of the Fourier transform of the temporal correlation function, which – at fixed imbalance – becomes increasingly small as $\xi _{\rho ,{\mathrm{th}}}$ decreases. The imbalance-dependent width of these tails also plays an important role in determining the strength of the suppression: as $\sigma _{\mathrm{c}}$ increases, the integrand below the cone contributes progressively less to the total heating, since the spectrum becomes more sharply peaked around $\tilde{k}_\|^{(1)}$ . This leads to the increased suppression seen in both $\mathcal{D}$ and $Q_\perp$ in figure 6. It should be emphasised that, although the integral for $Q_\perp$ in (4.15) includes a Gaussian weighting $\exp (-\tilde{v}_\|^2/\beta _{\mathrm{i}})$ arising from the assumed Maxwellian distribution, the observed suppression is primarily a consequence of the broadening of the wavevector–frequency spectrum discussed previously, rather than the form of the distribution (see § 4.2.2).

Figure 8.

Differential heating rate ${\mathrm{d}} Q_\perp /{\mathrm{d}} v_{\|}$ (4.31), which quantifies the heat given to ions with parallel velocities between $v_{\|}$ and $v_{\|}+{\mathrm{d}}v_{\|}$ , calculated using $\beta _{\mathrm{i}}=0.1$ , and $\sigma _{\mathrm{c}}= (a)\,0.0$ , (b) $0.75$ and (c) $0.999$ . The red dashed lines correspond to where $|\tilde{v}_{\|}|=v_{\mathrm{th}, \textrm{i}}/v_{\mathrm{A}}=\sqrt {\beta _{\mathrm{i}}}$ . These plots are normalised to the maximum value in the balanced case in panel (a).

4.3.2. Differential heating rate ${\mathrm{d}}^{\!}Q_\perp /{\mathrm{d}} v_\|$

An additional metric of ion heating is the differential heating rate ${\mathrm{d}} Q_\perp /{\mathrm{d}} v_\|$ , which quantifies the heat given to ions with parallel velocities between $v_\|$ and $v_\|+{\mathrm{d}}v_\|$ . This can be calculated from (4.15) where, after setting the upper integral limit to $\tilde{v}_\|$ , we have

(4.31) \begin{equation} \frac {1}{\varOmega _{\mathrm{i}}v_{\mathrm{th}, \textrm{i}}^3}\frac {{\mathrm{d}} Q_\perp }{{\mathrm{d}} v_\|}(\tilde{v}_\|) = \sqrt {\frac {4}{\pi }}e^{-\tilde{v}_\|^2/\beta _{\mathrm{i}}}(1+\sigma _{\mathrm{c}}\tilde{v}_\|)\mathcal{D}(\tilde{v}_\|). \end{equation}

Figure 8 plots (4.31) as a function of $\tilde{v}_\|$ and $\xi _{\rho ,{\mathrm{th}}}$ for $\beta _{\mathrm{i}}=0.1$ , and $\sigma _{\mathrm{c}}=0.0$ , $0.75$ and $0.999$ . All cases show increased heating at larger $\xi _{\rho ,{\mathrm{th}}}$ for ions with $|v_\||\lesssim v_{\mathrm{th}, \textrm{i}}$ (corresponding to a $\xi _{\rho ,{\mathrm{th}}}^3$ scaling). In the balanced case (figure 8 a), ${\mathrm{d}} Q_\perp /{\mathrm{d}} v_\|$ is symmetric around $v_\|=0$ due to the presence of equal populations of $\boldsymbol{z}^\pm$ fluctuations, with a moderate peak at $|v_\||\sim v_{\mathrm{th}, \textrm{i}}$ for small $\xi _{\rho ,{\mathrm{th}}}$ and at $|v_\||\sim 0$ for large $\xi _{\rho ,{\mathrm{th}}}$ . This symmetry is broken as $\sigma _{\mathrm{c}}$ increases, with the peak at large $\xi _{\rho ,{\mathrm{th}}}$ shifting slightly to positive $v_\|$ ; as the amplitude of $\boldsymbol{z}^-$ fluctuations decrease, ions with $v_\| \lt 0$ are heated less. This transition occurs slowly with increasing imbalance and only becomes starkly apparent when $\sigma _{\mathrm{c}}\gtrsim 0.9$ , as seen in figure 8(c). The tail at $\xi _{\rho ,{\mathrm{th}}}\ll 1$ is sharply cut off for particles with $\tilde{v}_\| \lt 0$ ; the remaining heating for $\tilde{v}_\| \lt 0$ particles at $\xi _{\rho ,{\mathrm{th}}}\sim 1$ arises from the fact that, although the peak of the correlation function at $|\tilde{k}_\|^{(1)}| \gt 1$ is outside the spectrum domain in (3.13), the finite width of the temporal correlation function allows for a contribution from modes with $|\tilde{k}_\|| \lt 1$ . This contribution disappears completely in the limit $\sigma _{\mathrm{c}}\to 1$ as the temporal correlation function approaches a delta function (cf. (4.22)).

5.1.1. Imbalanced turbulence and the helicity barrier

RMHD is well-suited for describing large-amplitude Alfvénic turbulent fluctuations relevant to ion heating as discussed in this work. However, it excludes key sub- $\rho _{\mathrm{i}}$ physics and, in particular, the helicity barrier: a mechanism that plays a key role in imbalanced, low- $\beta _{\mathrm{i}}$ turbulence (Meyrand et al. Reference Meyrand, Squire, Schekochihin and Dorland2021; Squire et al. Reference Squire, Meyrand, Kunz, Arzamasskiy, Schekochihin and Quataert2022, Reference Squire, Meyrand and Kunz2023; Adkins, Meyrand & Squire et al. Reference Adkins, Meyrand and Squire2024, Reference Adkins, Meyrand and Squire2025; Johnston et al. Reference Johnston, Squire and Meyrand2025). The helicity barrier arises due to the failure of the conservation of a gyrokinetic ‘generalised helicity’ at perpendicular scales above and below $\rho _{\mathrm{i}}$ in imbalanced turbulence (Schekochihin, Kawazura & Barnes Reference Schekochihin, Kawazura and Barnes2019; Meyrand et al. Reference Meyrand, Squire, Schekochihin and Dorland2021), which prevents turbulent energy from cascading to small perpendicular scales through $\rho _{\mathrm{i}}$ , trapping it at scales such that $k_\perp \rho _{\mathrm{i}}\lesssim 1$ , resulting in a decrease in the energy of $\rho _{\mathrm{i}}$ -scale fluctuations. The helicity barrier is well-supported observationally within the solar wind, manifesting as a steep ‘transition range’ drop in measured energy spectra (Leamon et al. Reference Leamon, Smith, Ness, Matthaeus and Wong1998; Bowen et al. Reference Bowen2020a ) as well as in other features (Bowen et al. Reference Bowen2022, Reference Bowen2024; McIntyre et al. Reference McIntyre, Chen, Squire, Meyrand and Simon2025).

Using simulations of a reduced gyrokinetic model able to capture the effects of the helicity barrier, Johnston et al. (Reference Johnston, Squire and Meyrand2025) showed that the heating of ions is unaffected by the barrier, despite the drop in energy of turbulent fluctuations at $k_\perp \rho _{\mathrm{i}}\sim 1$ scales. This is because ions are most efficiently heated by fluctuations whose frequencies approach $\varOmega _{\mathrm{i}}$ , which occurs at scales $\lambda \sim \rho _{\mathrm{i}}$ in balanced turbulence and shifts to $\lambda \gt \rho _{\mathrm{i}}$ in the presence of a helicity barrier (Johnston et al. Reference Johnston, Squire and Meyrand2025). Since the blocked cascade accumulates more power at these larger scales than at $\rho _{\mathrm{i}}$ -scales, the ion heating rate is maintained.

The effect of the helicity barrier on ion heating can be captured in the quasi-linear framework used in this paper by shifting the cutoff in the 2-D wavevector energy spectrum $\mathcal{E}_{\mathrm{2D}}$ to perpendicular scales $\lambda \gt \rho _{\mathrm{i}}$ (i.e. by changing the cutoff from $\tilde{k}_\perp \gt 1$ to $\tilde{k}_\perp \gt \rho _{\mathrm{i}}/\lambda$ in (3.13)), where the value of $\lambda$ can be measured from simulations of the helicity barrier. Because the fluctuations that ions interact with at this scale contain more power than those at $\rho _{\mathrm{i}}$ , the heating rate will be increased compared with the predictions in § 4.3 of imbalanced RMHD without a helicity barrier. Although this cutoff provides a good approximation to the helicity barrier for $k_\perp \rho _{\mathrm{i}} \lesssim 1$ fluctuations, the corresponding modification for the spectrum at $k_\perp \rho _{\mathrm{i}}\gtrsim 1$ (as would be needed in the calculations of Appendix C) is not clear.

5.1.2. Future work

Recent work has shown that intermittency – the heavy-tailed distribution of fluctuation amplitudes intrinsic to turbulence – can significantly affect the efficiency of perpendicular ion heating (Mallet et al. Reference Mallet, Klein, Chandran, Grošelj, Hoppock, Bowen, Salem and Bale2019; Cerri et al. Reference Cerri, Arzamasskiy and Kunz2021; Bowen et al. Reference Bowen, Ervin, Mallet, Chandran, Sioulas, Isenberg, Bale, Squire, Klein and Pezzi2025). Since our formalism relies on the wavevector–frequency spectrum of the turbulence, it does not directly capture such effects; extending it to account for intermittency remains an important direction for future work.

We note that the quasi-linear framework employed in this work assumes that cyclotron-resonant heating dominates in imbalanced turbulence, with ions directly resonating with the sharply peaked wavevector–frequency spectrum of the fluctuations. However, non-resonant heating of ions from wave-like fluctuations is also possible, given sufficiently large wave amplitude (Chen et al. Reference Chen, Lin and White2001; Johnson & Cheng Reference Johnson and Cheng2001). This mechanism cannot be captured using quasi-linear theory as it involves higher-order nonlinear interactions beyond the quadratic products retained in quasi-linear theory (see § 2 and Appendix B). Detailed comparisons between quasi-linear theory predictions and idealised single-wave simulations could help quantify these effects.

An interesting direction for future work is to investigate how ion heating is modified in the presence of purely weak turbulence, rather than in addition to strong turbulence as in this work. In weak turbulence, the wavevector–frequency spectrum is sharply peaked around the linear dispersion relation, even in the balanced limit (in contrast to the broadened spectrum of strong turbulence), reflecting the wave-like nature of its fluctuations (Meyrand, Galtier & Kiyani Reference Meyrand, Galtier and Kiyani2016). Our results predict cyclotron-resonant heating would dominate in this case for all $\sigma _{\mathrm{c}}$ . Exploring this regime by, for example, comparing heating in weak and strong turbulence with similar $\xi _{\rho ,{\mathrm{th}}}$ could therefore help to further disentangle the difference between stochastic heating and cyclotron-resonant heating.

While this work has focused primarily on the theoretical aspects of ion heating in turbulent plasmas, its predictions can be connected to space and astrophysical observations. Observations indicate that turbulent ion heating in the solar corona and wind is likely driven by a combination of stochastic heating (Bourouaine & Chandran Reference Bourouaine and Chandran2013; Klein & Chandran Reference Klein and Chandran2016; Vech, Klein & Kasper Reference Vech, Klein and Kasper2017; Martinović et al. Reference Martinović, Klein and Bourouaine2019, Reference Martinović, Klein and Bourouaine2020) and cyclotron-resonant heating by ICWs (Kasper et al. Reference Kasper, Maruca, Stevens and Zaslavsky2013; Bale et al. Reference Bale2019; Bowen et al. Reference Bowen, Ervin, Mallet, Chandran, Sioulas, Isenberg, Bale, Squire, Klein and Pezzi2020b , Reference Bowen2022, Reference Bowen2024). The heating rates associated with these mechanisms can be directly inferred from spacecraft measurements – for example, Bowen et al. (Reference Bowen, Ervin, Mallet, Chandran, Sioulas, Isenberg, Bale, Squire, Klein and Pezzi2025) explicitly quantify both the ICW- and stochastic-heating rates of ions in Parker Solar Probe data. Extending such analyses across a wider range of plasma conditions (varying imbalance $\sigma _{\mathrm{c}}$ , plasma beta $\beta _{\mathrm{i}}$ , turbulent amplitude $\xi _{\rho ,{\mathrm{th}}}$ and minor-ion composition) could allow direct tests of the general phenomenological ion-heating formula (5.1).