3.1. Change in perpendicular energy

We are interested in the change in $v_\perp ^2$ as $T\to \infty$ ,

(3.13) \begin{equation} \frac {v_\perp ^2}{v_{\perp 0}^2} = 1 + 2\varepsilon \left (\dot {X}_0\dot {X}_1 + \dot {Y}_0\dot {Y_1}\right )_{T\to \infty } + O(\varepsilon ^2). \end{equation}

Using (3.2) differentiated with respect to $T$ , (3.11) and (3.10),

(3.14) \begin{align} \varDelta = 2\left ( \dot {X}_0\dot {X}_1+\dot {Y}_0\dot {Y}_1\right )_{T\to \infty } = \frac {1}{\pi }\int _{-\infty }^\infty \int _{-\infty }^\infty \tilde {g}(K,0,\eta _K T') \sum _{n=-\infty }^\infty \frac {nJ_n(K)}{K} e^{inT'}\,\mathrm{d} K\,\mathrm{d} T'. \end{align}

Clearly, the contribution from $n=0$ in the sum vanishes. The integral is of the form

(3.15) \begin{equation} \varDelta = \frac {1}{\pi }\int _{-\infty }^\infty \int _{-\infty }^\infty \sum _{n=-\infty }^\infty A(n,K,\eta _KT')e^{inT'}\,\mathrm{d} K\,\mathrm{d} T', \end{equation}

which we can perform by closing the contour in the appropriate half-plane. The dominant contribution comes from the pole of $g(K,0,s)$ (say $s_*$ ) closest to the real axis, so that

(3.16) \begin{equation} \varDelta \sim 2i\int _{-\infty }^\infty \sum _{n=-\infty }^\infty \textrm {sgn}(n)\,\textrm {Res} [A(n, K, s), s_*]\exp \left (- \frac {|n \mathrm{Im}\left \{s_*\right \}|}{\eta _K}\right )\!\mathrm{d} K. \end{equation}

Since $A(n,K,s)=0$ for $n=0$ , if we have that $\eta _K\ll 1$ for all $K$ , this is exponentially small. At higher order in $\varepsilon$ , similar exponentially small expressions occur; this is a special case of the general conservation of adiabatic invariants to all orders (Kruskal Reference Kruskal1958, Reference Kruskal1962).

Moreover, if $\eta _K\ll 1$ for all $K$ , we need only keep the $n=\pm 1$ term, since it is obviously the largest. Therefore, we approximate $\varDelta$ as

(3.17) \begin{align} \varDelta &\approx \frac {2}{\pi } \int _{-\infty }^\infty \cos T' \int _{-\infty }^\infty \frac {J_1(K)}{K}\tilde {g}(K,0,\eta _K T')\mathrm{d} K\,\mathrm{d} T'\nonumber \\ &\sim c_1\int _{-\infty }^\infty \frac {J_1(K)}{K}\mathrm{Res}\left [\tilde {g}\right ]\exp (-c_2/\eta _K)\mathrm{d} K, \end{align}

where $c_1,c_2$ are (system-dependent) dimensionless constants of order unity and $Res[\tilde {g}]$ denotes the residue from the pole of $\tilde {g}$ closest to the real axis, and is a function of $K$ .

At this point, it is worth discussing when and how our solution breaks down. First, note that it is possible to have a situation where the exponential term arising from the pole is cancelled out by part of $g$ : for example, if $g(Y,Z,T) = \cos (T+\phi ) g'(Y,Z,T)$ . This is resonance and leads to the breakdown of the ordering if $g'$ is non-zero for a time $\delta T \sim O(1/\varepsilon )$ . We will assume this is not the case, but briefly discuss it in Appendix B.

Second, while we have shown that the first-order change in perpendicular kinetic energy (3.14) is exponentially small for $\eta _K\ll 1$ for all $K$ , the same is not true for our expressions for the perpendicular velocities (see (3.10)–(3.11)): the second integral in each case clearly has a non-zero $n=0$ term, meaning that if the fields are left on for a time $\varDelta T \gtrsim \varepsilon ^{-1}$ , the ordering of the solution will break down due to secularly growing terms in $\dot {Y}_1$ and $\dot {X}_1$ . This could be the case, for example, if the field varies so slowly that $\eta _K\sim \varepsilon$ ; hence, our formal restriction to larger $\eta _K$ in this section. These secular terms cancel out in the expression (3.14) for the change in kinetic energy. Unlike the case of a true resonance, this is simply a consequence of the naive perturbation method. In Appendix C, we use the Poincaré–Lindstedt method to extend the calculation to the case of arbitrarily small $\eta _K$ , showing that the expression (3.14) for the change in perpendicular kinetic energy does not change. This more involved calculation is therefore perhaps of more mathematical than physical interest, but is included for completeness.

3.2. Scale dependence

Because of the Bessel function, the contributions to $\varDelta$ from different perpendicular scale ${\sim} 1/K$ vary with $K$ . For small argument ( $K\ll 1$ ), $J_n(K) \approx (K/2)^n/n!$ , so that the only term that survives in (3.10) is $n=1$ , and we may replace $J_1(K)/K$ in (3.17) with a scale-independent factor $1/2$ . In the opposite limit of large argument $K\gg 1$ , the envelope of $|J_n(K)|\sim K^{-1/2}$ , and so all the terms in (3.8) become small. However, in turbulence, $\eta _K$ is typically an increasing function of $K$ , and the exponential suppression of the heating will be less effective for larger $K$ : the balance between these is system-specific, depending on $\eta _K$ and the $K$ -dependence of the Fourier amplitudes $\tilde {g}$ .

3.3. Scattering contours

Let us for the moment assume that the electromagnetic fluctuations are from a propagating wave or superposition of waves, with a parallel phase velocity $v_{\textrm{ph}}(K)$ , so that $\tilde {b}=\tilde {b}(K,Z-(v_{\textrm{ph}}(K)/v_{\perp 0})T)$ and $\tilde {g}=\tilde {g}(K,Z-(v_{\textrm{ph}}/v_{\perp 0})T)$ . The results derived in the previous sections do not require this, but it will allow us to make contact with the usual quasilinear theory of cyclotron heating. It may also be directly applicable to the turbulence in the solar wind, which can be highly imbalanced: dominated by outward-going Alfvén waves, for which $v_{\textrm{ph}} = v_{\textrm{A}}$ . More generally, this situation could potentially also apply to nonlinear solitary waves which have a single effective $v_{\textrm{ph}}$ due to the nonlinearity balancing dispersion (Kawahara Reference Kawahara1969; Kakutani & Ono Reference Kakutani and Ono1969; Hasegawa & Mima Reference Hasegawa and Mima1976; Mjølhus & Wyller Reference Mjølhus and Wyller1986; Mallet Reference Mallet2023). We note that we are ignoring the possibility of waves purely with phase velocity only in the perpendicular direction, and also electrostatic waves (which to have parallel phase velocities, must also have $E_z\neq 0$ so that Faraday’s law is satified).

In a frame moving at $v_{ \textrm{ref}}$ compared to the laboratory, to first order in $\varepsilon$ , as $T\to \infty$ , we have the energy change

(3.18) \begin{align} v_\perp ^2 + (v_\parallel - v_{\textrm {ref}})^2 &= v_{\perp 0}^2 + (v_{\parallel 0} - v_{ \textrm{ref}})^2 + 2 \varepsilon v_{\perp 0}^2\left [ \varDelta + 2 \dot {Z}_1 (v_{\parallel 0}-v_{\textrm {ref}})/v_{\perp 0}\right]\!. \end{align}

From Faraday’s law (2.9), we have

(3.19) \begin{equation} \tilde {g}= \frac {v_{\parallel 0} - v_{\textrm{ph}}(K)}{v_{\perp 0}} \tilde {b}, \end{equation}

so that (3.14) can be written as

(3.20) \begin{equation} \varDelta = \frac {1}{\pi } \int _{-\infty }^\infty \int _{-\infty }^\infty \frac {v_{\parallel 0} - v_{\textrm{ph}}(K)}{v_{\perp 0}}\tilde {b}(K,0,\eta _K T') \sum _{n=-\infty }^\infty \frac {nJ_n(K)}{K} e^{inT'}\,\mathrm{d} K\,\mathrm{d} T'. \end{equation}

Combining this expression with (3.12) as $T\to \infty$ , we find that

(3.21) \begin{align} &\varDelta + 2 Z_1(v_{\parallel 0} - v_{ \textrm{ref}})/v_{\perp 0}\nonumber\\ =& \frac {1}{\pi }\int _{-\infty }^\infty \int _{-\infty }^\infty \frac {v_{ \textrm{ref}} - v_{\textrm{ph}}(K)}{v_{\perp 0}}\tilde {b}(K,0,\eta _K T') \sum _{n=-\infty }^\infty \frac {nJ_n(K)}{K} e^{inT'}\,\mathrm{d} K\,\mathrm{d} T'. \end{align}

The integrand of this expression vanishes if $v_{\textrm{ph}}(K)=v_{ \textrm{ref}}$ . The left-hand side is the expression appearing in square brackets in (3.18). Therefore, for a propagating wave or coherent wavepacket (e.g. a soliton), diffusion occurs along the scattering contours $v_\perp ^2 + (v_\parallel - v_{\textrm{ph}})^2=\text{const}$ . This behaviour is lost if there is not a single $v_{\textrm{ph}}$ , as would be the case for a dispersive wavepacket where $v_{\textrm{ph}}(K)$ is not constant. This is also the case in strong, balanced Alfvénic turbulence, where while the linear and nonlinear frequencies are statistically in critical balance, so that $\partial _t \sim v_{\textrm{A}}\partial _z \sim u_x \partial _y$ , there is a broad range of effective frequencies; or equivalently, a distribution of effective phase velocities with mean zero and width ${\sim} v_{\textrm{A}}$ .

Obviously, for a particle moving at $v_{\parallel 0}=v_{\textrm{ph}}$ , the electric field is zero (again, provided that there is no electrostatic wave), the magnetic field is stationary in time and, thus, there is no change in the perpendicular or parallel energy of the particle. This is encoded in the fact that for a particle moving at $v_{\parallel 0}=v_{\textrm{ph}}$ , the phase of the wave is $z-v_{\textrm{ph}} t = z_0 + (v_{\parallel 0}-v_{\textrm{ph}})t=z_0$ , independent of $t$ , and so $\eta _K=0$ for all $K$ .

3.4. Example

Let us (for simplicity’s sake) assume that $\tilde {g}(K,0,\eta _KT)=f(\eta _K T) \tilde {h}(K)$ . As an example, we choose

(3.22) \begin{equation} f[\eta _K T] = \frac {1}{\pi }\left \{\arctan [\eta _K(T-a)]-\arctan [\eta _K(T-b)]\right \}\!. \end{equation}

We will leave the spatial pattern of the fluctuation $\tilde {h}(K)$ arbitrary since we are mainly interested in the dependence of the energy change on $\eta _K$ . This is a model of a fluctuation that ‘turns on’ at a rate $\eta$ around $t=a$ and then ‘turns off’ at the same rate at $t=b$ , and is plotted in figure 1. We need to perform integrals of the form

(3.23) \begin{equation} \int _{-\infty }^\infty e^{inT} f(\eta _K T)\, {\rm d}T = \frac {i}{n}\int _{-\infty }^\infty e^{inT} \frac {{\rm d}}{\mathrm{d} T}f(\eta _KT)\, \mathrm{d} T, \end{equation}

for integer $n\neq 0$ , and we have integrated by parts to get the second expression.Footnote ² We have

(3.24) \begin{equation} \frac {{\rm d}}{\mathrm{d} T}f(\eta _KT) = \frac {\eta _K}{\pi }\left [\frac {1}{\eta _K^2(T-a)^2+1}-\frac {1}{\eta _K^2(T-b)^2+1}\right]\!. \end{equation}

Figure 1. The functional form (3.22) for the time dependence of the electric field, with $a=-2$ , $b=2$ .

The poles are at $T=a\pm i/\eta$ , $T=b\pm i/\eta$ . We close in the appropriate half-plane, and obtain for $n\neq 0$ ,

(3.25) \begin{equation} \int _{-\infty }^\infty e^{inT} f(\eta _K T)\, {\rm d}T = \frac {i}{n}\left [e^{ina} -e^{inb}\right ]e^{-|n|/\eta _K}. \end{equation}

Then, the perpendicular energy change (3.14) is

(3.26) \begin{align} \varDelta = \frac {1}{\pi }\int _{-\infty }^\infty \tilde {h}(K)\sum _{\substack {n=-\infty \\ n\neq 0}}^\infty \frac {iJ_n(K)}{K} \left [e^{ina}-e^{inb}\right ]e^{-|n|/\eta _K}\,\mathrm{d} K. \end{align}

Two simplified limits are of interest: first, if $\eta _K\ll 1$ and the fields vary slowly compared with the gyrofrequency, the $n=1$ terms dominate (and even they are exponentially small):

(3.27) \begin{align} \varDelta \approx \frac {2}{\pi }\left (\sin b - \sin a\right )\int _{-\infty }^\infty \tilde {h}(K)\frac {J_1(K)}{K} e^{-|n|/\eta _K}\,\mathrm{d} K, \quad \eta _K\ll 1 \,\,\forall \,\, K. \end{align}

Second, if $\tilde {h}(K)$ only has power at small $K$ , we can again drop all but the $n=1$ term, as discussed earlier in § 3.2, and additionally $J_1(K)/K\to 1/2$ :

(3.28) \begin{align} \varDelta \approx \frac {1}{\pi }\left (\sin b - \sin a\right )\int _{-\infty }^\infty \tilde {h}(K) e^{-|n|/\eta _K}\,\mathrm{d} K, \quad \tilde {h}(K)\ll 1 \,\,\forall \,\, K\gtrsim 1. \end{align}

We can understand the dependence on $a$ and $b$ as depending on the phase of the particle’s orbit at some reference time. Assuming that the ion velocity distribution is gyrotropic, each individual particle is as likely to gain or lose energy from the interaction: the average $\overline {\varDelta }=0$ .

3.5. Diffusion coefficient and heating rate

We have so far derived an expression for the change in perpendicular energy of a single particle, $\varDelta$ (see 3.14) and derived its explicit form for an example (3.26). Importantly, $\varDelta$ can be both positive and negative, and in fact, the average over the initial gyrophase of the particle $\overline {\varDelta }$ vanishes.

Repeated interactions with coherent fluctuations will cause diffusion in energy. To be more precise, let us suppose for the moment that there are a large number of identical coherent fluctuations present and the particle encounters one approximately every $\delta t$ . Each interaction with a fluctuation occurs with a random initial gyrophase and provides an uncorrelated kick in perpendicular kinetic energy of magnitude $\delta \kappa = m_i v_{\perp 0}^2\varepsilon \varDelta /2$ . This leads to an energy diffusion coefficient

(3.29) \begin{equation} D \sim \frac {\delta \kappa ^2}{\delta t}\sim \frac {1}{4}\frac {m_i^2 v_{\perp 0}^4 \varepsilon ^2\overline {\varDelta ^2}}{\delta t}, \end{equation}

where the overline denotes averaging over the (uniform) gyrophase distribution of the particles.

If all the fluctuations are characterised by a single perpendicular scale $\lambda \sim 1/k_\perp =\rho /K$ ,Footnote ³ then the normalised time scale for the interaction is $\tau _\lambda \sim 1/\varOmega _i\eta _K$ . Let us now suppose that these fluctuations are rare: in each time interval of length $\tau$ , an encounter with a fluctuation occurs with probability $P$ . Then, the $\delta t$ appearing in the formula for the diffusion coefficient is clearly just $\delta t = \tau /P=1/(P\varOmega _i\eta _K)$ , i.e.

(3.30) \begin{equation} D \sim \frac {1}{4}\varOmega _i m_i^2 v_{\perp 0}^4 P\varepsilon ^2\eta _K\overline {\varDelta ^2}. \end{equation}

Now, consider the more realistic case where at each scale $\lambda$ , there is an ensemble of different fluctuations, each with their own $\varepsilon$ , $\overline {\varDelta ^2}$ and $\delta t$ : we may characterise this ensemble by a joint probability distribution $P_\lambda (\varepsilon ,\overline {\varDelta ^2},\delta t)$ , noting that the arguments need not be independent. Then, we can generalise (3.30): denoting the average over the distribution of fluctuations $P_\lambda$ with angle brackets,

(3.31) \begin{equation} D \sim \frac {\varOmega _i m_i^2 v_{\perp 0}^4}{4}\left \langle \varepsilon ^2\eta _K\overline {\varDelta ^2}\right \rangle . \end{equation}

In the case of fluctuations that are propagating (linear or nonlinear) waves, so that $\omega = v_{\textrm{ph}} k_\parallel$ , the diffusion will be along the scattering contours (§ 3.3).

Finally, we can use our expression (3.17) for $\varDelta$ to estimate

(3.32) \begin{equation} D \sim \varOmega _i m_i^2v_{\perp 0}^4 \frac {J_1^2(k_\perp \rho )}{k_\perp ^2\rho ^2} \varepsilon ^2_{k_\perp }\eta _{k_\perp }\exp \left (-\frac {2}{\eta _{k_\perp }}\right )\!. \end{equation}

As $\eta _{k_\perp } \to 1$ from below, the diffusion becomes strong. To get an overall effective heating rate per unit mass, suppose $v_{\perp 0}\sim v_{\textrm{th}}$ ; then,

(3.33) \begin{equation} Q_\perp \sim \frac {D(v_{\textrm{th}})}{m_i^2v_{\textrm{th}}^2} \sim \varOmega _i v_{\textrm{th}}^2 \frac {J_1^2(k_\perp \rho _{\textrm{th}})}{k_\perp ^2\rho _{\textrm{th}}^2} \varepsilon ^2_{k_\perp }\eta _{k_\perp }\exp \left (-\frac {2}{\eta _{k_\perp }}\right )\!. \end{equation}

We have derived this diffusion coefficient for the energy of a single particle interacting with a distribution of fluctuations. If we consider the whole population of ions, with ion velocity distribution function $f(v_\perp )$ , interacting with a single coherent fluctuation, the behaviour is also diffusive. If there is a gradient in $f(v_\perp )$ , while individual particles are just as likely to gain or lose energy, the flux of particles from the region with larger $f(v_\perp )$ will be larger than the flux from the region with smaller $f(v_\perp )$ , smoothing the gradient. As an illustration, suppose that the initial distribution is uniform in gyrophase, but confined to a single $v_{\perp 0}$ : a ring distribution. Afterwards,

(3.34) \begin{equation} v_\perp = v_{\perp 0}\sqrt {1+\varepsilon \varDelta }\approx v_{\perp 0}(1+\varepsilon _{k_\perp }\varDelta /2). \end{equation}

The variance of this distribution is then

(3.35) \begin{equation} \overline { v_\perp ^2} - v_{\perp 0}^2 \approx \varepsilon _{k_\perp }^2v_{\perp 0}^2\overline {\varDelta ^2}/4, \end{equation}

i.e. the effective temperature has changed by an amount of order

(3.36) \begin{equation} \delta T_\perp \sim {m_i}\varepsilon _{k_\perp }^2v_{\perp 0}^2\overline {\varDelta ^2} \sim m_i v_{\perp 0}^2 \frac {J_1^2(k_\perp \rho )}{k_\perp ^2\rho ^2}\varepsilon _{k_\perp }^2 \exp (-2/\eta _{k_\perp }), \end{equation}

where we have used (3.17) to estimate $\varDelta$ . As $Q_\perp \sim \delta T_\perp /\delta t$ and $\delta t\sim \eta _{k_\perp }\varOmega _i$ , (3.33) and (3.36) agree with each other.

In the rest of the paper, we will use the theory already described to study ion heating in Alfvénic turbulence (using (3.33), since over a long time period, each particle will interact with many fluctuations) and reconnection (using (3.36), since the particles only interact with a reconnection exhaust once). To do so, it is necessary to have on hand estimates of $\varepsilon _{k_\perp }$ and $\eta _{k_\perp }$ , the accuracy of the latter being more critically important: due to its presence inside the exponential cutoff, our cavalier disregard of coefficients of order unity might lead to large inaccuracies in the estimated heating rates. For this reason, in much of the rest of the paper, we will (following the approach of Chandran et al. Reference Chandran, Li, Rogers, Quataert and Germaschewski2010) insert adjustable constants parametrising these unknown coefficients in our estimates for $\varepsilon _{k_\perp }$ and $\eta _{k_\perp }$ . Given a detailed enough knowledge of the system’s dynamics, it would in principle be possible to derive these coefficients from first principles; more practically, one can fit them numerically (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), although care must be taken to consider the unrealistically limited scale separation possible in simulations of turbulence and reconnection.

There are a few different approaches to estimating $\eta _{k_\perp }$ . In the first, we estimate $\eta _1=|\omega - k_\parallel v_{\parallel 0}|/\varOmega _i$ , where $\omega$ is some linear or nonlinear frequency of the system. If $\beta _i\ll 1$ and we have Alfvénic fluctuations with $\omega \sim k_\parallel v_{\textrm{A}}\gg k_\parallel v_{\parallel 0}$ , the $\omega$ term tends to dominate. In the second, we estimate $\eta _2\sim \varepsilon K$ as the (inverse of the) time it takes to $E\times B$ drift out of the structure, assuming that, in reality, the fields have structure in the $\hat {\boldsymbol{x}}$ direction too. This can only possibly be relevant once $k_\perp \rho _{\textrm{th}}\gtrsim 1$ , since for $k_\perp \rho _{\textrm{th}}\ll 1$ , the fields are frozen into the plasma flow.

Finally, one might think to estimate the time it takes for the polarisation drift ( ${\sim} \varepsilon _{k_\perp }\eta _{k_\perp }$ ) to cause the particle to leave the structure in the $\hat {\boldsymbol{y}}$ direction (McChesney et al. Reference McChesney, Stern and Bellan1987; Chen, Lin & White Reference Chen, Lin and White2001; White, Chen & Lin Reference White, Chen and Lin2002), $\eta _3 \sim \varepsilon \eta _{k_\perp } K$ , where probably $\eta _{k_\perp }\sim \eta _1,\eta _2$ . As can be seen, this is only comparable to $\eta _1$ or $\eta _2$ if $\varepsilon _{k_\perp } K \sim 1$ , i.e. only at very large amplitude.

To preview the approach of the next two sections, in Alfvénic turbulence (§ 4), we will find that $\eta _1\sim \eta _2$ , while in reconnection (§ 5), we will use $\eta _2$ exclusively, assuming little structure in the parallel direction.

4.1. Balanced turbulence

We will first examine the case of ‘balanced’ turbulence, where the fluxes of Alfvénic fluctuations propagating parallel and antiparallel to the magnetic field are comparable: the imbalanced case is rather different and is discussed in § 4.3. Typical scalings for the fluctuation amplitudes $\delta b_k$ and $c\delta E_k/B_0$ in balanced turbulence are plotted in figure 2. At small scales, $k_\perp \rho _p\gg 1$ , we have $\delta b_k \sim k_\perp ^{-2/3}$ or steeper (Schekochihin et al. Reference Schekochihin, Cowley, Dorland, Hammett, Howes, Quataert and Tatsuno2009; Boldyrev & Perez Reference Boldyrev and Perez2012; Zhou, Liu & Loureiro Reference Zhou, Liu and Loureiro2023); say $\delta b_k \sim \delta b_* (k_\perp \rho _p)^{-a}$ , where $\delta b_*$ is then the amplitude at $k_\perp \rho _p = 1$ . Then,

(4.11) \begin{equation} Q_{\perp p} \sim (k_\perp \rho _p)^{-3a+1}\frac {\delta b_*^3}{\rho _{\textrm{th}}} \exp \left (-\frac {c_2 \varOmega _p}{k_\perp ^{2-a}\rho _p^{1-a}\delta b_*}\right )\!, \quad k_\perp \rho _p\gg 1. \end{equation}

This reaches a maximum when

(4.12) \begin{equation} \frac {c_2\varOmega _p}{k_\perp ^{2-a}\rho _p^{1-a}\delta b_*} = \frac {3a-1}{2-a}, \end{equation}

or at

(4.13) \begin{equation} k_{max}\rho _p = \left (\frac {2-a}{3a-1}\frac {c_2\varOmega _p}{\omega _*}\right )^{1/(2-a)}, \end{equation}

where $\omega _* = \delta b_*/\rho _p$ is the frequency of gyroscale fluctuations. The maximum proton heating rate therefore occurs when $\eta _{k_\perp } = \omega /\varOmega _p \sim 1$ (in an order-of-magnitude sense). Here, we remind the reader that in reality, this estimate will become inaccurate since our estimates will fail around $k_\perp d_e \sim 1$ , where the spectrum steepens again (Stawarz et al. Reference Stawarz2019) and our estimates for $\alpha _k$ also break down (Adkins, Meyrand & Squire Reference Adkins, Meyrand and Squire2024).

Figure 3. The proton heating rate $Q_{\perp p}$ normalised to the turbulent energy flux through scales $\epsilon =\delta b_L^3/L$ for different values of the normalised outer-scale amplitude $A$ , with $L/\rho _p = 10^{4}$ and $v_{\textrm{th p}}/v_{\textrm{A}}=0.1$ . The horizontal black line denotes $Q_{\perp p}/\epsilon =1$ , complete damping of the turbulent cascade: in reality, if the heating approaches this line, the power-law behaviour of the spectra and the constancy of $\epsilon$ will no longer be accurate. The vertical solid black line denotes $k_\perp \rho _p=1$ and the vertical dashed line denotes $k_\perp \rho _e=1$ : the small heating rates at or beyond the electron scales in our model (which neglects electron-scale physics) are an overestimate due to the much steeper spectrum in this range.

At large scales, $k_\perp \rho _p\ll 1$ , the magnetic fluctuations scale as a power law, between $\delta b_k \propto k_\perp ^{-1/3}$ (Goldreich & Sridhar Reference Goldreich and Sridhar1995; Chen et al. Reference Chen, Mallet, Yousef, Schekochihin and Horbury2011) and $\delta b_k\propto k_\perp ^{-1/4}$ (Boldyrev Reference Boldyrev2006; Mallet et al. Reference Mallet, Schekochihin, Chandran, Chen, Horbury, Wicks and Greenan2016; Chen et al. Reference Chen2020), where the exact scaling likely depends on the regime of turbulence. Then, $Q_{\perp p}$ is an increasing function of $k_\perp$ . In balanced turbulence, there is typically a smooth join between the $k_\perp \rho _p\ll 1$ and $k_\perp \rho _p\gg 1$ spectra, and so we expect the heating rate as a function of $k_\perp$ to reach a maximum at relatively small scales, when $\omega \sim \varOmega _p$ . This agrees (qualitatively at least) with the peak of the heating rate observed in the hybrid-kinetic numerical simulations of Arzamasskiy et al. (Reference Arzamasskiy, Kunz, Chandran and Quataert2019) and Cerri et al. (Reference Cerri, Arzamasskiy and Kunz2021). We have plotted the proton heating rate (4.8) as a function of (inverse) scale $k$ in figure 3 for several different outer scale amplitudes $A=\delta b_L/v_{\rm {A}}$ , assuming that $L/\rho _p = 10^{4}$ and $v_{\textrm{th p}}/v_{\textrm{A}}=0.1$ . The energy flux into the turbulent cascade is defined as $\epsilon = \delta b_L^3/L$ . The trend agrees with our previous analysis: if the amplitude is high enough that $Q_{\perp p}/\epsilon \sim 1$ at $k_\perp \rho _p\sim 1$ , the peak heating is at the proton gyroradius scale (e.g. the red line in figure 3). At lower amplitudes, the heating is less efficient and peaks at a smaller scale (e.g. the blue line in figure 3). The oscillations in $Q_{\perp p}$ when $k_\perp \rho _p \gt 1$ are due to the Bessel function; in reality, structures will have contributions from a range of $k_\perp \sim 1/\lambda$ and this behaviour will be smoothed out.

4.2. Proton heating fraction in balanced turbulence

The maximum proton heating rate is

(4.14) \begin{equation} Q_{\perp max,p}\sim \frac {\delta b_*^3}{\rho _p} \left (\frac {\omega _*}{\varOmega _p}\right )^{s}, \end{equation}

ignoring prefactors of order unity, and with $s=(3a-1)/(2-a)$ . Now, suppose the turbulent cascade (in the absence of dissipation) had a constant energy flux through scale $\epsilon$ ; by a Kolmogorov-style argument, dimensionally (as a reminder, we are neglecting intermittency in the distribution of fluctuation amplitudes),

(4.15) \begin{equation} \epsilon \sim \delta b_*^3/\rho _p. \end{equation}

If $\omega _*/\varOmega _p\ll 1$ , $Q_{\perp max}/\epsilon \ll 1$ and there is no significant ion heating; the energy must be dissipated at smaller scales onto the electrons. Writing $\omega _* \sim \delta b_*/\rho _p$ ,

(4.16) \begin{equation} \frac {Q_{\perp max,p}}{\epsilon } \sim \left (\frac {\delta b_*}{v_{\textrm{th p}}}\right )^{s}. \end{equation}

With $a=3/4$ (a not unreasonable value given the observed spectrum, e.g. Chen et al. Reference Chen, Horbury, Schekochihin, Wicks, Alexandrova and Mitchell2010), $s=1$ and this agrees with the expression for the ion heating fraction given by Matthaeus et al. (Reference Matthaeus, Parashar, Wan and Wu2016). It may be more useful to write this in terms of the amplitude at the outer scale $L$ of the turbulence, parametrised as

(4.17) \begin{equation} \delta b_L = A v_{\textrm{A}}, \end{equation}

where $v_{\textrm{A}}$ is the Alfvén velocity. Assuming $\delta b_k \propto k_\perp ^{-b}$ for $k_\perp \rho _i\ll 1$ (with $b=1/3$ corresponding to Goldreich & Sridhar Reference Goldreich and Sridhar1995 and $b=1/4$ corresponding to Boldyrev Reference Boldyrev2006 including dynamic alignment), we have

(4.18) \begin{equation} \delta b_* \sim A v_{\textrm{A}} \left (\frac {\rho _p}{L}\right )^{b} \end{equation}

and

(4.19) \begin{equation} \frac {Q_{\perp max,p}}{\epsilon } \sim A^{s}\beta _p^{-s/2}\left (\frac {\rho _p}{L}\right )^{sb}, \end{equation}

where the proton plasma beta $\beta _p = v_{\textrm{th p}}^2/v_{\textrm{A}}^2$ . Inserting $a=3/4$ ( $s=1$ ) for simplicity,

(4.20) \begin{equation} \frac {Q_{\perp max,p}}{\epsilon } \sim A\beta ^{-1/2}\left (\frac {\rho _p}{L}\right )^{b}. \end{equation}

This estimate can be compared with other mechanisms, for example, with the approach outlined by Howes (Reference Howes2024): for a dissipation mechanism to be important, we require the turbulent system to pass a threshold in parameter space beyond which $Q_{\perp max}/\epsilon \sim 1$ . Here, this threshold is described in terms of $\beta$ and $\rho _p/L$ . For perpendicular ion heating to be important, we need (taking $b=1/4$ )

(4.21) \begin{equation} \frac {\rho _p}{L} \gtrsim A^{-4} \beta ^{2}, \end{equation}

which is the same dependence on $\beta$ as found by Howes (Reference Howes2024) for stochastic heating (Chandran et al. Reference Chandran, Li, Rogers, Quataert and Germaschewski2010). Typically, $\rho _p\ll L$ : for example, in the solar wind, $\rho _p/L \sim 10^{-4}$ , while in the interstellar medium (ISM), $\rho _p/L \sim 10^{-11}$ . Thus, this simple estimate suggests that ion heating should be negligible in the ISM, and accounts for only $5\,\%{-}10\,\%$ of the energy budget in the $\beta \sim 1$ solar wind, in stark contrast with the available evidence (Cranmer et al. Reference Cranmer, Matthaeus, Breech and Kasper2009). This suggests that we need to incorporate additional physics into our model: in § 4.4, we show that intermittency (Mallet et al. Reference Mallet, Klein, Chandran, Grošelj, Hoppock, Bowen, Salem and Bale2019) results in much higher ion heating fractions, because fluctuations attain larger amplitudes such that $\omega _*/\varOmega \sim 1$ , leading to $Q_{\perp max,p} \sim \epsilon$ .

We have neglected all heating apart from that at $k_{max}$ : while amending this might increase the heating rates slightly, in practice, because of the exponential suppression of heating from low-frequency fluctuations, $Q_{\perp p}$ is quite sharply peaked at $k_{max}$ .

4.3. Forced imbalanced Alfvénic turbulence and the helicity barrier

Meyrand et al. (Reference Meyrand, Squire, Schekochihin and Dorland2021) found that in numerical simulations of forced imbalanced turbulence, the helicity barrier causes energy to build up at scales larger than $\rho _p$ with a spectral break to a very steep ( $\delta b_k\propto k_\perp ^{-3/2}$ or steeper) spectrum beyond $k_h\sim 2 (1-\sigma _c)^{1/4}/\rho _p\lesssim 1/\rho _p$ in the ‘transition range’, before returning to a shallower spectrum at $k_\perp \rho _p\sim 1$ . We can study this situation with our heating rate expression (4.9) also. Writing $\delta b_k \sim \delta b_h (k_\perp /k_h)^{-3/2}$ and taking $k_\perp \rho _p\lesssim 1$ , the frequency $\omega \sim k_\perp \delta b_k\propto k_\perp ^{-1/2}$ in the transition range, as a decreasing function of $k$ . Inserting into our expression for $Q_{\perp p}$ for $k_\perp \rho _p\ll 1$ , we find

(4.22) \begin{equation} Q_{\perp p} \sim k_h \delta b_h^3 \left (\frac {k_\perp }{k_h}\right )^{-7/2} \exp \left ( - \frac {c_2 \varOmega _p}{k_h \delta b_h (k_\perp /k_h)^{-1/2}}\right )\!, \end{equation}

which decreases with increasing $k_\perp$ in the transition range. The assumption $\delta b_k\propto k_\perp ^{-3/2}$ could be weakened; $Q_\perp$ always decreases with $k$ so long as $\delta b_k \propto k_\perp ^{-1}$ or steeper. Thus, the maximum heating occurs at $k_h$ and is given by

(4.23) \begin{equation} Q_{\perp max} \sim k_h \delta b_h^3 \exp \left (-\frac {c_2\varOmega }{k_h\delta b_h}\right )\!. \end{equation}

If $k_h\delta b_h\ll \varOmega$ , then $Q_{\perp max}/\epsilon \sim \exp (-c_2\varOmega / k_h\delta b_h)\ll 1$ . If energy is being injected into the turbulent cascade at large scales, this situation cannot be in steady state due to the presence of the helicity barrier, and thus $\delta b_h$ (and the whole large-scale spectrum) must be growing in time. A steady state can be achieved once $Q_{\perp max}\sim \epsilon \sim k_h\delta b_h$ , which happens when $c_2\varOmega / k_h\delta b_h \sim 1$ , or roughly $k_h\delta b_h \sim \varOmega$ . Thus, in forced imbalanced turbulence, essentially all the energy flux from the turbulent cascade goes into ion heating. A more refined analysis, taking into account the electron inertia effects entering at $d_e$ , shows that in fact there is a critical level of imbalance below which the helicity barrier will not form (Adkins et al. Reference Adkins, Meyrand and Squire2024): we assume that in the systems in which we are interested (for example, the solar corona), $1\gg \beta _e\gg m_e/m_i$ and the turbulence is extremely imbalanced, such that we can ignore this effect.

This picture of ion heating ‘switching on’ due to the helicity barrier is, as we mentioned previously, not new, and already predicted by Meyrand et al. (Reference Meyrand, Squire, Schekochihin and Dorland2021) and Squire et al. (Reference Squire, Meyrand, Kunz, Arzamasskiy, Schekochihin and Quataert2021). Even more similar to this work, Johnston et al. (Reference Johnston, Squire and Meyrand2025) found using test particle simulations that ion heating in imbalanced turbulence depends on a phenomenological exponential suppression factor, controlled by the fluctuation amplitude at the scale at which the maximum frequency is reached, i.e. the above-mentioned transition-range break scale $k_h$ .

Physical systems of imbalanced turbulence, for example, the solar wind, are not in fact homogeneous nor forced in the usual way, and so the helicity barrier may not have sufficient time to reach steady state and cause as sharp a transition to ion heating as suggested based on the forced numerical simulations (Meyrand et al. Reference Meyrand, Squire, Schekochihin and Dorland2021; Squire et al. Reference Squire, Meyrand, Kunz, Arzamasskiy, Schekochihin and Quataert2021). The scaling for the maximum heating rate (4.23) still applies, since the helicity barrier does still cause a steep transition range at scales larger than $\rho _p$ .

4.4. Intermittency

So far, we have assumed that the turbulence is characterised by a single amplitude at each scale, $\delta b_k$ . In reality, $\delta b_k$ is a random variable at each scale, typically with a heavy large-amplitude tail; both in solar wind observations (Salem et al. Reference Salem, Mangeney, Bale and Veltri2009; Zhdankin, Boldyrev & Mason Reference Zhdankin, Boldyrev and Mason2012; Sioulas et al. Reference Sioulas2022) and in numerical simulations (Mallet et al. Reference Mallet, Schekochihin and Chandran2015, Reference Mallet, Schekochihin, Chandran, Chen, Horbury, Wicks and Greenan2016; Zhdankin et al. Reference Zhdankin, Boldyrev and Chen2016a , Reference Zhdankin, Boldyrev and Uzdenskyb ). This can dramatically increase the ion heating fraction (Mallet et al. Reference Mallet, Klein, Chandran, Grošelj, Hoppock, Bowen, Salem and Bale2019; Cerri et al. Reference Cerri, Arzamasskiy and Kunz2021).

As an (extreme) toy example, suppose that the turbulence were characterised by fluctuations of a fixed amplitude $\delta b$ , filling a certain scale-dependent fraction $f_k=(k_\perp L)^{-d}$ of the volume at scale $1/k_\perp$ . In other words, the distribution of fluctuation amplitudes is

(4.24) \begin{equation} \delta b_k = \begin{cases} \delta b &\text{with probability } (k_\perp L)^{-d},\\ 0 &\text{with probability } 1-(k_\perp L)^{-d}. \end{cases} \end{equation}

Requiring the energy flux $\epsilon = \langle k_\perp \delta b_k^3\rangle = k_\perp f_k\delta b^3$ to be independent of $k$ , we have $d=1$ . Then, the root-mean-square value of the fluctuation amplitude measured at scale $1/k_\perp$ is

(4.25) \begin{equation} \delta b_{rms, k} = \sqrt {\langle \delta b_k^2 \rangle } = \delta b (k_\perp L)^{-1/2}. \end{equation}

Meanwhile, the overall heating rate at large scales, given by (4.9), is

(4.26) \begin{equation} \langle Q_{\perp p}\rangle /\epsilon \sim c_1 \exp \left (-\frac {c_2\varOmega _p}{k_\perp \delta b}\right )\!, \end{equation}

where the average is over the distribution of fluctuation amplitudes. Writing this solely in terms of $\delta b_{rms,k}$ ,

(4.27) \begin{equation} \langle Q_{\perp p} \rangle /\epsilon \sim c_1 \exp \left (-\frac {c_2\varOmega _p (k_\perp L)^{-1/2}}{k_\perp \delta b_{rms,k}}\right )\!. \end{equation}

Since $k_\perp L \gg 1$ , this dramatically increases the overall heating rate for a given observed $\delta b_{rms,k}$ , compared with the estimate without taking account of the intermittency of $\delta b_k$ . This is not, in fact, a realistic model of intermittent Alfvénic turbulence – we have included it here to show, in a transparent way, the difference that intermittency makes to the efficiency of ion heating.

It is worth mentioning that intermittency models including dynamic alignment (e.g. Chandran et al. Reference Chandran, Schekochihin and Mallet2015; Mallet & Schekochihin Reference Mallet and Schekochihin2017) do not necessarily increase the ion heating fraction, because the dynamic alignment between $\delta \boldsymbol{z}^\pm$ increases the nonlinear time scale of the turbulent structures, decreasing the effectiveness of the heating due to the exponential suppression factor. This decrease in ion heating fraction due to dynamic alignment is the opposite to what was found earlier by Mallet et al. (Reference Mallet, Klein, Chandran, Grošelj, Hoppock, Bowen, Salem and Bale2019), who used an exponential suppression factor based solely on the amplitude (Chandran et al. Reference Chandran, Li, Rogers, Quataert and Germaschewski2010), rather than the rate of change of the fluctuations, thus finding that dynamic alignment acted to increase the stochastic heating rate.

A more promising intermittency model in this regard (which nevertheless obtains the same $-3/2$ perpendicular spectrum) is the reflection-driven turbulence model of Chandran et al. (Reference Chandran, Sioulas, Bale, Bowen, David, Meyrand and Yerger2025), in which dynamic alignment does not play a role. There, larger-amplitude fluctuations naturally have higher frequencies, as required to make the ion heating more effective. Another approach, independent of any particular intermittency model, would be to use in situ measurements of the distribution of turbulent fluctuations to calculate the heating rate $\langle Q_\perp \rangle$ directly. The Parker Solar Probe has recently started to explore the plasma environment very close to the sun (Kasper et al. Reference Kasper2021), where ion heating is thought to be particularly important (Kasper & Klein Reference Kasper and Klein2019), and it will be interesting to assess the ion heating in this newly explored regime.

4.5. Minor ions

Observationally, as mentioned in § 1, minor ions appear to be heated even more strongly than the protons. Returning to (4.7) and writing the ion gyrofrequency in terms of the proton gyrofrequency, $\varOmega _i = Z_im_p/m_i \varOmega _p$ ,

(4.28) \begin{equation} Q_\perp \sim c_1 \frac {J_1^2(k_\perp \rho _{\textrm{th}})}{k_\perp ^2\rho _{\textrm{th}}^2} \frac {k_\perp ^4\rho _p^4}{\alpha _k[1-\varGamma _0(k_\perp ^2\rho _p^2/2)]^2} k_\perp \delta b_k^3 \exp \left (-\frac {c_2\varOmega _p}{k_\perp \alpha _k\delta b_k}\frac {Z_im_p}{m_i}\right )\!. \end{equation}

Since $Z_i m_p/m_i \lt 1$ , the exponential suppression is less effective for minor ions, exponentially increasing the heating rate relative to the protons. This is similar to what happens in the stochastic heating model applied to minor ions (Chandran Reference Chandran2010). Note that $\rho _{th} = \rho _p (v_{th i}/v_{\textrm{th p}})(Z_i m_p/m_i)$ , so that the gyroradius of the ions is typically larger than that of the protons: however, in the case of imbalanced turbulence where the cascade is already cut off due to the helicity barrier (§ 4.3), this may have less of an effect. An exception is in forced homogeneous numerical simulations, where the saturated state of the turbulence has such a high amplitude that fluctuations attain the minor ion’s gyrofrequency at large scales, causing extreme heating of minor species (Zhang et al. Reference Zhang, Kunz, Squire and Klein2025).