The effects of finite electron inertia on helicity-barrier-mediated turbulence

Abstract

Understanding the partitioning of turbulent energy between ions and electrons in weakly collisional plasmas is crucial for the accurate interpretation of observations and modelling of various astrophysical phenomena. Many such plasmas are ‘imbalanced’, wherein the large-scale energy input is dominated by Alfvénic fluctuations propagating in a single direction. In this paper, we demonstrate that when strongly-magnetised plasma turbulence is imbalanced, nonlinear conservation laws imply the existence of a critical value of the electron plasma beta (the ratio of the thermal to magnetic pressures) that separates two dramatically different types of turbulence in parameter space. For betas below the critical value, the free energy injected on the largest scales is able to undergo a familiar Kolmogorov-type cascade to small scales where it is dissipated, heating electrons. For betas above the critical value, the system forms a ‘helicity barrier’ that prevents the cascade from proceeding past the ion Larmor radius, causing the majority of the injected free energy to be deposited into ion heating. Physically, the helicity barrier results from the inability of the system to adjust to the disparity between the perpendicular-wavenumber scalings of the free energy and generalised helicity below the ion Larmor radius; restoring finite electron inertia can annul, or even reverse, this disparity, giving rise to the aforementioned critical beta. We relate this physics to the ‘dynamic phase alignment’ mechanism (that operates under yet lower beta conditions and in pair plasmas), and characterise various other important features of the helicity barrier, including the nature of the nonlinear wavenumber-space fluxes, dissipation rates, and energy spectra. The existence of such a critical beta has important implications for heating, as it suggests that the dominant recipient of the turbulent energy, ions or electrons, can depend sensitively on the characteristics of the plasma at large scales.

1 Introduction

Many astrophysical plasma systems are weakly collisional, with their characteristic dynamical timescales approaching those associated with inter-particle collisions. A key question in the context of such plasmas is what determines the partitioning of turbulent free energy between ions and electrons, given that they lack an obvious means of thermal equilibration. Indeed, two-temperature states are expected or observed in a variety of contexts, e.g. accretion discs around black holes (Ichimaru 1977; Quataert & Gruzinov 1999), the intracluster medium (Takizawa 1999; Kunz, Jones & Zhuravleva 2022), and the solar wind (Cranmer 2009). In the latter context, the Alfvénic turbulence launched by low-frequency motions in the corona (De Pontieu et al. 2007; Tomczyk et al. 2007) is observed to preferentially heat protons over electrons (Hansteen & Leer 1995; Cranmer et al. 2009; Bandyopadhyay et al. 2023), while heavier minor ions (e.g. helium or oxygen) are heated even more efficiently (Kohl et al. 1997). This is somewhat puzzling, however, as theories of Alfvénic turbulence at low plasma beta (the ratio of the thermal to magnetic pressures) predict that all of the heating from a Kolmogorov (1941) style cascade of free energy occurs on electrons (Howes et al. 2008a; Schekochihin et al. 2009; Howes 2010; Kawazura, Barnes & Schekochihin 2019; Schekochihin, Kawazura & Barnes 2019). Furthermore, such a cascade is unable to easily transfer energy to the small parallel scales required to excite ion-cyclotron waves (ICWs) that can cause efficient perpendicular ion heating (Kennel & Engelmann 1966; Isenberg & Vasquez 2011). Although there are other possible mechanisms that will preferentially heat ions, they either impose constraints on the turbulence that are contradicted by observations or the understanding thereof remains limited. For example, compressive fluctuations are able to cause parallel ion heating (Schekochihin et al. 2009), but are unlikely to have sufficient power to explain the temperature difference between ions and electrons (Howes et al. 2012). In a similar vein, mechanisms such as random-walk scattering from ion-Larmor-radius-scale electric-field fluctuations (so-called ‘stochastic heating’; Chandran et al. 2010) and sub-ion-Larmor-radius kinetic-Alfvén-wave (KAW) turbulence (Arzamasskiy et al. 2019; Isenberg & Vasquez 2019) are capable of dissipating a significant fraction of the turbulent energy so long as fluctuations remain of sufficiently large amplitude. Whether solar-wind fluctuations are capable of doing this robustly enough to explain the observed ion heating remains unclear (Howes et al. 2008b; Chandran et al. 2011).

A possible explanation for the preferential heating of ions that also explains a number of other solar-wind observations is the so-called ‘helicity barrier’ (Meyrand et al. 2021): when the turbulence is imbalanced (i.e. when there is a significant disparity in the energies of the forwards- and backwards-propagating fluctuations), as it is in the solar wind, free energy is prevented from cascading past the ion Larmor radius $\rho _{\rm i}$ and reaching smaller perpendicular scales, where it would, presumably, dissipate on electrons. Instead, the turbulence grows to large amplitudes, creating fine parallel structure that excites ICW fluctuations and heats the ions, which then absorb the majority of the injected power (Squire et al. 2022; Squire, Meyrand & Kunz 2023). This helicity-barrier-mediated turbulence has many features that agree with measurements of the low-beta solar wind, including those of the ion velocity distribution function (see, e.g., Marsch 2006; He et al. 2015; Bowen et al. 2022), helicity (Huang et al. 2021; Zhao et al. 2021), and properties of the steep spectral slopes of the electromagnetic fields in the ‘transition range’ on scales comparable to the ion Larmor radius. These transition range spectra have been observed for decades (Leamon et al. 1998; Alexandrova et al. 2009; Sahraoui et al. 2009) and, more recently, by Parker Solar Probe (PSP) (Bowen et al. 2020a; Duan et al. 2021; Bowen et al. 2024). The helicity barrier may also have an important impact on plasmas in other astrophysical contexts, such as in the interpretation of images from the Event Horizon Telescope (Wong & Arzamasskiy 2024).

This paradigm remains incomplete, however, as research exploring the impact of the helicity barrier in imbalanced solar-wind turbulence has thus far been conducted without accounting for the effects of electron kinetics. In particular, the assumption of a vanishing electron inertial scale $d_{\rm e} = \rho _{\rm e}/\sqrt {\beta _{\rm e}}$ (where $\rho _{\rm e}$ is the electron Larmor radius and $\beta _{\rm e}$ is the electron plasma beta) has resulted in the neglect of electron Landau damping, which is particularly significant on scales comparable to $d_{\rm e}$ (Schekochihin et al. 2009; Zocco & Schekochihin 2011; Zhou, Liu & Loureiro 2023a). Given that these effects will undoubtedly play a role in determining the partitioning of turbulent energy between ions and electrons, accounting for them is a necessary extension of the helicity-barrier paradigm to the low-beta regime most relevant to the lower corona. In this paper, we consider the effects of finite electron inertia in imbalanced helicity-barrier-mediated turbulence. Using equations derived in a low-beta asymptotic limit of gyrokinetics, we demonstrate the existence of a critical value of the electron beta $\beta _{\rm e}$ below which, for a given value of the energy imbalance in the outer-scale fluctuations, the helicity barrier will not form, allowing free energy to cascade to small perpendicular scales. This effect is shown to arise from the constraints placed on the turbulent dynamics by the simultaneous conservation of both free energy and (generalised) helicity across the ion Larmor scale $\rho _{\rm i}$ and the electron inertial scale $d_{\rm e}$. At perpendicular scales significantly larger than $\rho _{\rm i}$ or smaller than $d_{\rm e}$, both the free energy and the helicity display the same perpendicular-wavenumber scaling, meaning that they can, in principle, cascade simultaneously. This is not the case on scales between $\rho _{\rm i}$ and $d_{\rm e}$: below the aforementioned critical beta, the helicity exhibits a shallower scaling with perpendicular wavenumber than the free energy, whereas above it, it exhibits a steeper scaling. In the former case, fluctuations are able to compensate for this disparity in scaling by becoming increasingly misaligned at small scales through ‘dynamic phase alignment’ (Loureiro & Boldyrev 2018; Milanese et al. 2020), allowing the cascade of energy to proceed to small scales uninterrupted. In the latter case, however, fluctuations cannot account for this disparity because they are unable to become more than maximally aligned, eventually leading to the breakdown of the constant-flux solution and the formation of a helicity barrier that prevents energy from cascading past $\rho _{\rm i}$. The existence of this critical beta thus has important consequences for turbulent heating, acting as a ‘switch’ that determines whether the majority of turbulent fluctuations are dissipated on ions (above the critical beta) or electrons (below it) as a function of the equilibrium parameters of the system.

The remainder of this paper is organised as follows. Section 2 motivates and outlines the model equations used for theoretical arguments and simulations throughout this work (§ 2.1), before considering their linear phase velocity (§ 2.2) and nonlinearly conserved invariants (§ 2.3). Details of the numerical implementation are briefly discussed in § 2.4. Section 3 considers the dynamics of imbalanced turbulence within these model equations. We begin by outlining a theory of balanced turbulence assuming a constant-flux cascade of energy (§ 3.1), which exhibits good agreement with numerical simulations. The effect of helicity conservation is then considered in § 3.2, from which the critical value of $\beta _{\rm e}$ discussed previously is shown to follow. The principal characteristics of the helicity barrier are outlined in § 3.3 in order to motivate our testing of this critical beta that occurs in § 3.4, with the prediction being robustly supported by numerical simulations. We then briefly discuss the relationship between the helicity barrier and dynamic phase alignment, before summarising and discussing the implications of our results in § 4.

2 Isothermal kinetic reduced electron heating model

Standard magnetised plasma turbulence phenomenologies (see, e.g., Goldreich & Sridhar 1995; Boldyrev 2006; Schekochihin 2022) imply that fluctuations at scales much smaller than those associated with the variation of the local magnetic field $\boldsymbol {B}_0 = B_0 \boldsymbol {b}_0$ are highly anisotropic in space, viz., they satisfy $k_{\parallel } \ll k_{\perp }$ for characteristic wavenumbers $k_{\parallel }$ and $k_{\perp }$ parallel and perpendicular to $\boldsymbol {b}_0$. This anisotropy, which appears to be satisfied in the solar wind (Chen et al. 2013; Chen 2016), allows even small-scale fluctuations to have frequencies well below the frequency $\varOmega _s$ of the Larmor motion of the particles. Averaging, then, the Vlasov–Maxwell system over this fast Larmor timescale in the presence of such anisotropy leads to the gyrokinetic system of equations (Howes et al. 2006; Schekochihin et al. 2009; Abel et al. 2013), which has seen recent application to the study of kinetic plasma turbulence in astrophysical plasmas (see, e.g., Howes et al. 2008b, 2011; Kawazura et al. 2019). A further simplification can be made by expanding in the limit of low plasma beta, wherein there is minimal coupling between Alfvénic and ion-compressive fluctuations because the ion-thermal speed is much lower than the Alfvén speed (Schekochihin et al. 2019). This allows ion kinetics to be neglected even at the ion-Larmor scale. The resultant system of equations is known as the ‘kinetic reduced electron heating model’ (KREHM, Zocco & Schekochihin 2011; Loureiro et al. 2016), which couples the equations of reduced magnetohydrodynamics (RMHD, Kadomtsev & Pogutse 1974; Strauss 1976) and electron reduced magnetohydrodynamics (ERMHD, Schekochihin et al. 2009; Boldyrev et al. 2013) to the electron kinetic physics. Although the KREHM equations are formally derived with the electron beta ordered comparable to the electron–ion mass ratio $\beta _{\rm e} = 8{\rm \pi} n_{0{\rm e}} T_{0{\rm e}}/B_0^2 \sim m_{\rm e}/m_{\rm i}$ (Zocco & Schekochihin 2011; Adkins et al. 2022), they remain valid for all $\beta _{\rm e} \ll 1$ assuming an order-unity equilibrium-temperature-ratio between ions and electrons $\tau \equiv T_{0{\rm i}}/T_{0{\rm e}} \sim 1$ ($n_{0{\rm e}}$ is the equilibrium density of electrons).

In pursuit of further simplicity, we will assume the electrons to be isothermal along exact (equilibrium plus perturbed) field lines. Formally, this amounts to neglecting the parallel gradient of the parallel electron-temperature perturbation that would otherwise appear in the electron momentum equation (see (2.2)). The effect of this is to decouple the lowest-order fluid moments from the remainder of the kinetic hierarchy (higher-order velocity moments of the kinetic distribution function). Such an approximation is easily justified when investigating dynamics on timescales much slower than the electron parallel-streaming rate, which is typically the case at scales comparable to, or larger than, the ion-Larmor radius $k_{\perp } \rho _{\rm i} \lesssim 1$ (see, e.g., Schekochihin et al. 2009, 2019; Abel & Cowley 2013; Zielinski et al. 2017). We note, however, that the typical dynamical timescale associated with the resultant isothermal KREHM equations becomes comparable to the electron parallel-streaming rate on scales around the electron-inertial scale $k_{\perp } d_{\rm e} \sim 1$ (see § 2.2). This means that the isothermal approximation breaks down at smaller perpendicular scales, where the effects of electron Landau damping are significant (Zhou et al. 2023a). Nevertheless, even though the isothermal approximation cannot be formally derived, we view the isothermal KREHM equations (which are three-dimensional, having neglected the kinetic physics) as a useful, perhaps even necessary, intermediate step between the simpler RMHD or ERMHD systems and the more complicated KREHM system (which is four-dimensional). In particular, it allows us to investigate the dynamics of the ‘fluid’ moments of the system without Landau damping obscuring important features of the dynamics. We will discuss the limitations of the isothermal approximation in detail in § 4.1.

2.1 Model equations

The equations of isothermal KREHM are

(2.1)\begin{align} &\frac{\text{d} \delta n_{\rm e}}{\text{d} t} + n_{0{\rm e}} \boldsymbol{\nabla}_{{\parallel}} u_{{\parallel} {\rm e}} =0, \end{align}

(2.2)\begin{align}& m_{\rm e} n_{0{\rm e}} \frac{\text{d} u_{{\parallel} {\rm e}}}{\text{d} t} + T_{0{\rm e}} \boldsymbol{\nabla}_{{\parallel}} \delta n_{\rm e} ={-} e n_{0{\rm e}}\left(\frac{1}{c} \frac{\partial A_\parallel}{\partial t} + \boldsymbol{\nabla}_{{\parallel}} \phi \right), \end{align}

where $v_{\text {the}} = \sqrt {2T_{0{\rm e}}/m_{\rm e}}$ is the thermal speed of electrons, $m_{\rm e}$ their mass, and $-e$ their charge. Equation (2.1) is the electron continuity equation, describing the advection of the perturbed electron density $\delta n_{\rm e}$ by the $\boldsymbol {E} \times \boldsymbol {B}$ flow due to the perturbed electrostatic potential $\phi$:

(2.3)\begin{equation} \frac{\text{d}}{\text{d} t} = \frac{\partial}{\partial t} + \boldsymbol{u}_{{\perp}} \boldsymbol{\cdot} \boldsymbol{\nabla}_{{\perp}}, \quad \boldsymbol{u}_{{\perp}} = \frac{c}{B_0} \boldsymbol{b}_0 \times \boldsymbol{\nabla}_{{\perp}} \phi, \end{equation}

and their compression or rarefaction due to the perturbed parallel electron flow $u_{\parallel {\rm e}}$ along the exact magnetic field. The latter includes the perturbation of the magnetic-field direction arising from the parallel component of the magnetic-vector potential ${A_{\parallel }}$:

(2.4)\begin{equation} \boldsymbol{\nabla}_{{\parallel}} = \boldsymbol{b}\boldsymbol{\cdot} \boldsymbol{\nabla} = \frac{\partial}{\partial z} + \frac{{\delta \! \boldsymbol{B}_{\!\perp}}}{B_0} \boldsymbol{\cdot} \boldsymbol{\nabla}_{{\perp}}, \quad {\delta \! \boldsymbol{B}_{\!\perp}} ={-} \boldsymbol{b}_0 \times \boldsymbol{\nabla}_{{\perp}} {A_{{\parallel}}}. \end{equation}

Equation (2.2) is the electron parallel momentum equation, consisting of a balance between electron inertia, the parallel pressure gradient (which, here, is simply the parallel density gradient, due to the assumption of isothermality), and the parallel electric field appearing on the right-hand side. Since an electron flow uncompensated by an ion flow is a current (the ion thermal speed is formally small), $u_{\parallel {\rm e}}$ is related to ${A_{\parallel }}$ by Ampére's law:

(2.5)\begin{equation} -e n_{0{\rm e}} u_{{\parallel} {\rm e}} = j_\parallel= \frac{c}{4{\rm \pi}} \boldsymbol{b}_0 \boldsymbol{\cdot} \left( \boldsymbol{\nabla}_{{\perp}} \times {\delta \! \boldsymbol{B}_{\!\perp}}\right) \quad \Rightarrow \quad u_{{\parallel} {\rm e}} = \frac{c}{4{\rm \pi} e n_{0{\rm e}}} \boldsymbol{\nabla}_{{\perp}}^2 {A_{{\parallel}}}. \end{equation}

Finally, the electron-density perturbation is related to $\phi$ by quasineutrality:

(2.6)\begin{equation} \frac{\delta n_{\rm e}}{n_{0{\rm e}}} = \frac{\delta n_{\rm i}}{n_{0{\rm i}}} ={-} \bar{\tau}^{{-}1} \frac{e \phi}{T_{0{\rm e}}} \equiv{-} \frac{Z}{\tau}(1 - \hat{\varGamma}_0) \frac{e \phi}{T_{0{\rm e}}}, \end{equation}

where the operator $\hat {\varGamma }_0$ can be expressed, in Fourier space, in terms of the modified Bessel function of the first kind: $\varGamma _0 = {\rm I}_0(\alpha _{\rm i}) e^{-\alpha _{\rm i}}$, where $\alpha _{\rm i} = (k_{\perp } \rho _{\rm i})^2/2$. This becomes $1 - \hat {\varGamma }_0 \approx -\rho _{\rm i}^2 \boldsymbol {\nabla }_{\perp }^2/2$ at large scales $k_{\perp } \rho _{\rm i} \ll 1$, and $1 - \hat {\varGamma }_0 \approx 1$ at small scales $k_{\perp } \rho _{\rm i} \gg 1$; the former limit is why (2.6) is sometimes referred to as the gyrokinetic Poisson equation.

Using (2.5) and (2.6), we can write (2.1) and (2.2) as

(2.7)\begin{align} &\frac{\text{d}}{\text{d} t} \bar{\tau}^{{-}1} \frac{e\phi}{T_{0{\rm e}}} - \frac{c}{4 {\rm \pi}e n_{0{\rm e}}} \boldsymbol{\nabla}_{{\parallel}} \boldsymbol{\nabla}_{{\perp}}^2 {A_{{\parallel}}} = 0, \end{align}

(2.8)\begin{align}&\frac{\text{d}}{\text{d} t}\left({A_{{\parallel}}} - d_{\rm e}^2 \boldsymbol{\nabla}_{{\perp}}^2 {A_{{\parallel}}}\right) ={-} c \left(\frac{\partial \phi}{\partial z} + \boldsymbol{\nabla}_{{\parallel}} \bar{\tau}^{{-}1} \phi\right). \end{align}

Together, (2.7) and (2.8) form a closed pair of equations describing the evolution of the electromagnetic potentials $\phi$ and ${A_{\parallel }}$ of a strongly-magnetised, low-beta plasma in the absence of electron Landau damping. At scales $k_{\perp } d_{\rm e} \ll 1$, we recover from (2.7) and (2.8) the finite-Larmor-radius-MHD (FLR-MHD) system studied in Meyrand et al. (2021), which itself reduces to RMHD and ERMHD at large ($k_{\perp } \rho _{\rm i} \ll 1$) and small ($k_{\perp } \rho _{\rm i} \gg 1$) scales, respectively. In the ultra-low-beta limit $\beta _{\rm e} \ll m_{\rm e}/m_{\rm i}$, in which the electron-inertial length becomes larger than the ion-Larmor radius $d_{\rm e} \gg \rho _{\rm i}$, we also recover equations describing inertial-Aflvèn-wave turbulence considered by Loureiro & Boldyrev (2018) and Milanese et al. (2020). To reiterate, (2.7) and (2.8) can be simply obtained from the KREHM system (which itself can be derived from gyrokinetics in the low-beta limit) by neglecting higher-order moments of the electron kinetic distribution function. Thus, our model equations share all of KREHM's physical characteristics apart from electron Landau damping.

2.2 Phase velocity

Linearising and Fourier-transforming (2.7) and (2.8), we find forwards- and backwards-propagating modes of frequency $\omega = \pm k_{\parallel } v_{\text {ph}} (k_{\perp })$, where the perpendicular-wavenumber-dependent phase velocity is given by

(2.9)\begin{equation} v_{\text{ph}}(k_{{\perp}}) = k_{{\perp}} \rho_{s} \left(\frac{1 + \bar{\tau}}{1 + k_{{\perp}}^2 d_{\rm e}^2}\right)^{1/2}v_{\rm A} . \end{equation}

Here, $\rho _{s} = \sqrt {Z/2\tau } \rho _{\rm i}$ is the ion-sound radius, related to the (thermal) sound speed by ${c_{\rm s} = \rho _{s}/\varOmega _{\rm i}}$ and $v_{\rm A} = B_0/\sqrt {4{\rm \pi} n_{0{\rm i}} m_{\rm i}}$ is the Alfvén speed. Equation (2.9) has the asymptotic behaviour

(2.10)\begin{equation} \lim_{k_{{\perp}} \rightarrow 0} v_{\text{ph}}(k_{{\perp}}) = v_{\rm A}, \quad \lim_{k_{{\perp}} \rightarrow \infty} v_{\text{ph}}(k_{{\perp}}) =\left(\frac{1+\tau/Z}{2}\right)^{1/2} v_{\text{the}}, \end{equation}

meaning that we recover Alfvén waves and the electron parallel-streaming rate at the largest and smallest scales, respectively. Its behaviour in the intermediate region, however, depends on the relative sizes of the ion-sound radius $\rho _{s}$ and electron-inertial scale $d_{\rm e}$, a competition that is controlled (ignoring any $\tau$ or $Z$ dependence) by the ratio of the electron beta to the electron–ion mass ratio $\beta _{\rm e} (m_{\rm e}/m_{\rm i})^{-1}$. In particular, (2.9) is an increasing function of (perpendicular) wavenumber for $\beta _{\rm e} \gg m_{\rm e}/m_{\rm i}$, and a decreasing one for $\beta _{\rm e} \ll m_{\rm e}/m_{\rm i}$, viz.,

(2.11)\begin{equation} v_{\text{ph}}(k_{{\perp}}) \propto \left\{ \begin{array}{@{}ll@{}} \displaystyle k_{{\perp}} \rho_{s} , & \beta_{\rm e} \gg m_{\rm e}/m_{\rm i}, \\ \displaystyle (k_{{\perp}} d_{\rm e})^{{-}1}, & \beta_{\rm e} \ll m_{\rm e}/m_{\rm i}, \end{array} \right. \end{equation}

for wavenumbers in the intermediate region. The associated waves are known as kinetic Alfvén waves (KAWs), for $\beta _{\rm e} \gg m_{\rm e}/m_{\rm i}$, and inertial Alfvén waves, for $\beta _{\rm e} \ll m_{\rm e}/m_{\rm i}$. This behaviour of the phase velocity is manifest in figure 1, where we plot (2.9) for different values of $\beta _{\rm e}$.¹ Note that for $\beta _{\rm e} \sim m_{\rm e}/m_{\rm i}$, the curve is non-monotonic; see the inset of figure 1. Crucially, the increase, or otherwise, of (2.9) with $k_{\perp }$ determines when the helicity barrier must form, and so whether or not a significant fraction of energy is able to cascade towards small scales; this is discussed in § 3.2.

Figure 1. The phase velocity (2.9), normalised to the Alfvén speed, plotted as a function of perpendicular wavenumber $k_{\perp } \rho _{\rm i}$, and for $\tau = Z =1$. The colours indicate the value of $\beta _{\rm e} (m_{\rm e}/m_{\rm i})^{-1}$ for a given line, with the solid black line showing the FLR-MHD case ($d_{\rm e} \rightarrow 0$). The dotted lines are the scalings (2.11). The vertical shaded region indicates wavenumbers $k_{\perp } \rho _{\rm e} > 1$ for which the model ceases to apply. The inset panel shows the case of $\beta _{\rm e} = m_{\rm e}/m_{\rm i}$, with the horizontal dashed line indicating $v_{\text {ph}}/v_{\rm A} = 1$.

The eigenfunctions associated with the forwards- and backwards-propagating modes can be expressed, in Fourier space, as

(2.12)\begin{equation} {\varTheta_{\boldsymbol{k}}^\pm} \equiv \sqrt{1 + k_{{\perp}}^2 d_{\rm e}^2 }\left[\frac{c}{B_0} \frac{v_{\text{ph}}(k_{{\perp}})/v_{\rm A}}{(k_{{\perp}} \rho_{s})^2} \bar{\tau}^{{-}1} {\phi_{\boldsymbol{k}}} \mp \frac{{A_{{\parallel} \boldsymbol{k}}}}{\sqrt{4 {\rm \pi}n_{0{\rm i}} m_{\rm i}}}\right]. \end{equation}

Apart from the prefactor multiplying the square brackets and the difference in $v_{\text {ph}}(k_{\perp })$, this definition is identical to that adopted in Meyrand et al. (2021). These generalised Elsässer potentials have the property that on the largest scales $k_{\perp } \ll \rho _{\rm i}^{-1}, d_{\rm e}^{-1}$, they reduce to the standard RMHD Elsässer potentials (Elsässer 1950), viz.,

(2.13)\begin{equation} \lim_{k_{{\perp}} \rightarrow 0} \boldsymbol{b}_0 \times \boldsymbol{\nabla}_{{\perp}} \varTheta^\pm= \boldsymbol{z}^\pm \equiv \boldsymbol{u}_{{\perp}} \pm \frac{{\delta \! \boldsymbol{B}_{\!\perp}}}{\sqrt{4{\rm \pi} n_{0{\rm i}} m_{\rm i}}}. \end{equation}

These potentials (2.12) provide a natural basis for our investigation of imbalanced turbulence in isothermal KREHM.

2.3 Nonlinear invariants

Most turbulent systems possess at least one nonlinear invariant: a quantity that is conserved by nonlinear interactions but may have localised sources (e.g. forcing or equilibrium gradients) and sinks (e.g. viscosity or particle collisions). Gyrokinetics conserves the so-called free energy, which is the sum of quadratic norms of the magnetic perturbations, and the perturbations of the distribution functions of both ions and electrons. In the isothermal KREHM system (2.7) and (2.8), the free energy takes the form (Zocco & Schekochihin 2011; Loureiro et al. 2016; Adkins et al. 2022):

(2.14)\begin{equation} W = \int \frac{\text{d}^3 \boldsymbol{r}}{V} \left[\frac{e^2 n_{0{\rm e}}}{2 T_{0{\rm e}}} \left(\phi \bar{\tau}^{{-}1} \phi\right) + \frac{e^2 n_{0{\rm e}}}{2 T_{0{\rm e}}} \left(\bar{\tau}^{{-}1} \phi \right)^2 + \frac{ \left| \boldsymbol{\nabla}_{{\perp}} {A_{{\parallel}}}\right|^2 + d_{\rm e}^2 \left(\boldsymbol{\nabla}_{{\perp}}^2 {A_{{\parallel}}} \right)^2}{8 {\rm \pi}} \right], \end{equation}

where $V$ is the plasma volume. The contributions to (2.14) are, from left to right, the energies associated with perturbations of the electrostatic potential, electron density, (perpendicular) magnetic field and electron parallel velocity. At large scales $k_{\perp } \ll \rho _{\rm i}^{-1}, d_{\rm e}^{-1}$, this becomes

(2.15)\begin{equation} W \approx \frac{n_{0{\rm i}} m_{\rm i}}{2} \int \frac{\text{d}^3 \boldsymbol{r}}{V} \left(\left| \boldsymbol{u}_{{\perp}} \right|^2 + \frac{\left|{\delta \! \boldsymbol{B}_{\!\perp}}\right|^2}{4 {\rm \pi}n_{0{\rm i}} m_{\rm i}}\right) = \frac{n_{0{\rm i}} m_{\rm i}}{4} \int \frac{\text{d}^3 \boldsymbol{r}}{V} \left( \left|\boldsymbol{z}^+\right|^2 + \left|\boldsymbol{z}^-\right|^2 \right), \end{equation}

recovering the usual expression for the free energy in RMHD (see, e.g., Schekochihin et al. 2009).

Free energy is normally the quantity whose cascade from large (injection) to small (dissipation) scales determines the properties of the plasma's turbulent state (Schekochihin et al. 2008, 2009), as in hydrodynamic turbulence (see, e.g., Alexakis & Biferale 2018, and references therein). Isothermal KREHM, however, possesses another invariant that also plays an important role: the generalised helicity,

(2.16)\begin{equation} H ={-}\frac{e^2 n_{0{\rm e}} v_{A}}{cT_{0{\rm e}}} \int \frac{\text{d}^3 \boldsymbol{r}}{V} \bar{\tau}^{{-}1}\phi \left({A_{{\parallel}}}- d_{\rm e}^2 \boldsymbol{\nabla}_{{\perp}}^2 {A_{{\parallel}}} \right). \end{equation}

This reduces to the MHD cross-helicity at $k_{\perp } \ll \rho _{\rm i}^{-1}, d_{\rm e}^{-1}$, viz.,

(2.17)\begin{equation} H \approx n_{0{\rm i}} m_{\rm i} \int \frac{\text{d}^3 \boldsymbol{r}}{V} \frac{\boldsymbol{u}_{{\perp}} \boldsymbol{\cdot} {\delta \! \boldsymbol{B}_{\!\perp}}}{\sqrt{4 {\rm \pi}n_{0{\rm i}} m_{\rm i}}} = \frac{n_{0{\rm i}} m_{\rm i}}{4} \int \frac{\text{d}^3 \boldsymbol{r}}{V} \left( \left|\boldsymbol{z}^+\right|^2 - \left|\boldsymbol{z}^-\right|^2 \right) , \end{equation}

and is proportional to the magnetic helicity at $\rho _{\rm i}^{-1} \ll k_{\perp } \ll d_{\rm e}^{-1}$ (due to perpendicular pressure balance; see Schekochihin et al. 2009). The presence of this generalised helicity places an additional constraint on the dynamical states accessible by the system, as the turbulence must now evolve in such a way as to conserve both (2.14) and (2.16) simultaneously. The remainder of this paper is devoted to demonstrating how and when these constraints give rise to different turbulent states, and the consequences that this may have for plasma heating.

2.4 Numerical setup

In what follows, the isothermal KREHM system (2.7)–(2.8) is solved using a modified version of the pseudospectral code TURBO (Teaca et al. 2009) in a triply-periodic box of size $L_x = L_y = L_z = L$ with $n_\perp ^2 \times n_z$ Fourier modes. A third-order Williamson (1980) algorithm is used for the time-stepping. Time is measured in units of the parallel Alfvén time $t_{\rm A} = L_z/v_{\rm A}$. Hyperdissipation in both the perpendicular and parallel directions is introduced by replacing the time derivative of the left-hand sides of (2.7) and (2.8) by

(2.18)\begin{equation} \frac{\text{d}}{\text{d} t} + \nu_\perp \boldsymbol{\nabla}_{{\perp}}^6 + \nu_z \frac{\partial^6}{\partial z^6}. \end{equation}

The coefficients $\nu _\perp$ and $\nu _z$ are adaptive, viz., they are re-evaluated at each timestep to ensure that dissipation occurs near the grid scale, maximising the inertial range (details of the numerical implementation can be found in Meyrand et al. 2024). Fluctuations are forced at large scales at $k_{\perp } = 4{\rm \pi} /L$, $|k_z| = 2{\rm \pi} /L$ through the form of negative damping (Meyrand et al. 2021); this method allows the rates of free-energy and helicity injection to be controlled exactly while producing sufficiently random motions to generate turbulence. All of the simulations listed in table 1 have $\tau = Z = 1$, though we will retain dependencies on these parameters in analytical expressions for the sake of completeness.

Table 1. The parameters used for the isothermal KREHM simulations considered in this paper. All simulations have $\tau = Z = 1$. Values in parentheses indicate the minimum and maximum values for the corresponding column, with the final column (‘Sims’) indicating the number of simulations in a given set. A dash in an entry indicates that the physical simulation being considered does not contain that physical parameter.

3 Imbalanced Alfvénic turbulence

Observations show that solar-wind turbulence is imbalanced, meaning it is energetically dominated by outward-propagating Alfvénic structures associated with $\boldsymbol {z}^+$ (inward-propagating structures are associated with $\boldsymbol {z}^-$). It must, therefore, possess a non-zero cross-helicity (2.17). Writing the free energy (2.14) and (generalised) helicity (2.16) in terms of the generalised Elsässer potentials (2.12) as

(3.1)\begin{align} W &= \frac{n_{0{\rm i}} m_{\rm i}}{4} \sum_{\boldsymbol{k}} \left( \left|k_{{\perp}}{\varTheta_{\boldsymbol{k}}^+}\right|^2 + \left|k_{{\perp}} {\varTheta_{\boldsymbol{k}}^-}\right|^2 \right), \end{align}

(3.2)\begin{align}H &= \frac{n_{0{\rm i}} m_{\rm i}}{4} \sum_{\boldsymbol{k}} \frac{ \left|k_{{\perp}}{\varTheta_{\boldsymbol{k}}^+}\right|^2 - \left|k_{{\perp}} {\varTheta_{\boldsymbol{k}}^-}\right|^2 }{v_{\text{ph}}(k_{{\perp}})/v_{\rm A}}, \end{align}

it is clear that the same will be true of the isothermal KREHM system given a difference in the ${\varTheta ^\pm }$ energies. We quantify this energy imbalance by the ratio of the free energy to the helicity

(3.3)\begin{equation} \tilde{\sigma}_c = \frac{H}{W}, \end{equation}

which reduces to the normalised cross-helicity (or RMHD imbalance) $\sigma _c$ at large scales:

(3.4)\begin{equation} \lim_{k_{{\perp}} \rightarrow 0} \tilde{\sigma}_c = \sigma_c = \frac{\int \,\text{d}^3 \boldsymbol{r} \left( \left|\boldsymbol{z}^+\right|^2 - \left|\boldsymbol{z}^-\right|^2 \right)}{\int \,\text{d}^3 \boldsymbol{r} \left( \left|\boldsymbol{z}^+\right|^2 + \left|\boldsymbol{z}^-\right|^2 \right)} \leqslant 1. \end{equation}

Measured values of (3.4) often exceed $|\sigma _c| \gtrsim 0.8$ in the solar-wind, particularly in near-Sun regions (McManus et al. 2020). Despite this, a comprehensive theory of imbalanced turbulence remains elusive, even in the (simpler) context of RMHD (Lithwick, Goldreich & Sridhar 2007; Chandran 2008; Beresnyak & Lazarian 2009b; Perez & Boldyrev 2009; Chandran & Perez 2019; Schekochihin 2022). As such, the phenomenological theory presented in the following sections lays no claim to being comprehensive; instead, it should be viewed as a useful framework through which the effect of finite electron inertia can be explored in imbalanced turbulence.

3.1 Constant-flux cascade

Consider the case where (free) energy (2.14) and (generalised) helicity (2.8) are injected into our isothermal KREHM system at constant rates $\varepsilon _{W}$ and $\varepsilon _{H}$, respectively, by some large-scale stirring of turbulent fluctuations (due to, e.g., reflection of outwards-propagating fluctuations; Velli, Grappin & Mangeney 1989). We denote the resultant injection imbalance, the ratio of the injected flux of helicity to that of the energy, as $\sigma _\varepsilon = |\varepsilon _{H}|/\varepsilon _{W}$. Given that both the energy and helicity are nonlinear invariants, the only available route to their dissipation is through some nonlinear transfer to small scales. Motivated by this, we assume, for the moment, that there is a local, Kolmogorov (1941) style cascade that carries a constant flux of injected energy and helicity from the outer (injection) scale, through some putative inertial range, to the dissipation scale.

It follows immediately from this that the rates of energy injection into the forward- and backward-propagating fluctuations

(3.5)\begin{equation} \varepsilon^\pm= \frac{\varepsilon_{W} \pm \varepsilon_{H}}{2} =\frac{1 \pm \sigma_{\varepsilon}}{2}\varepsilon_{W}, \end{equation}

will be equal to the associated flux of ${\varTheta ^\pm }$ energy through the inertial range. We estimate these energy fluxes from (3.1) as

(3.6)\begin{equation} \frac{1}{n_{0{\rm e}} T_{0{\rm e}}} \frac{\text{d} W_\pm}{\text{d} t} \sim \left(t_{\text{nl}}^\pm \right)^{{-}1} \frac{\left(k_{{\perp}} {\varTheta_{k_{{\perp}}}^\pm} \right)^2}{c_{\rm s}^2} \sim \varepsilon^\pm= \text{const.}, \end{equation}

where here, and in what follows, ${\varTheta _{k_{\perp }}^\pm }$ refers to the characteristic amplitude of the Elsässer potentials at the scale $k_{\perp }^{-1}$, rather than to the Fourier transform of the field (2.12). Formally, ${\varTheta _{k_{\perp }}^\pm }$ can be defined by

(3.7)\begin{align} \frac{\left(k_{{\perp}} {\varTheta_{k_{{\perp}}}^\pm} \right)^2}{c_{\rm s}^2} = \frac{1}{n_{0{\rm e}}T_{0{\rm e}}}\int_{k_{{\perp}}}^\infty \,\text{d} k_{{\perp}}' E_\perp^\pm (k_{{\perp}}') , \quad E_\perp^\pm(k_{{\perp}}) = 2{\rm \pi} k_{{\perp}} \int_{-\infty}^\infty \,\text{d} k_{{\parallel}}\frac{n_{0{\rm i}}m_{\rm i}}{4}\left\langle \left| k_{{\perp}} {\varTheta_{\boldsymbol{k}}^\pm} \right|^2 \right\rangle, \end{align}

where $E_\perp ^\pm (k_{\perp })$ is the one-dimensional perpendicular energy spectrum of ${\varTheta ^\pm }$ (cf. (3.1)), and in which the brackets denote an ensemble average. An alternative definition would be via a second-order structure function (see, e.g., Davidson 2013). Perturbations of other quantities, such as potential, velocity, and magnetic field, will similarly be taken to refer to their characteristic amplitude at a given perpendicular scale.

In order to proceed, we need an expression for the nonlinear times appearing in (3.6). However, as we will discuss shortly, determining exactly what these nonlinear times are is not a straightforward task, remaining an open research question even in the RMHD regime (Schekochihin 2022). As such, let us henceforth consider the balanced regime, assuming that the rates of energy injection into the backwards- and forwards-propagating fluctuations are comparable, and that they have the same scaling with perpendicular wavenumber, such that their ratio is constant at all scales, viz.,

(3.8)\begin{equation} {\varepsilon^{+}} \sim {\varepsilon^{-}} \sim \varepsilon_{W}, \quad \frac{{\varTheta_{k_{{\perp}}}^+}}{{\varTheta_{k_{{\perp}}}^-}} = \text{const}. \end{equation}

Then, the nonlinear times in (3.6) can straightforwardly be taken to be the nonlinear $\boldsymbol {E} \times \boldsymbol {B}$ advection rate associated with either field, which, comparing (2.3) and (2.12), and neglecting any possible anisotropy in the perpendicular plane, can be written as

(3.9)\begin{equation} \left(t_{\text{nl}}^\pm \right)^{{-}1} \sim k_{{\perp}} u_\perp{\sim} \varOmega_{\rm i} \left(\frac{\bar{\tau}^2}{1+\bar{\tau}}\right)^{1/2} \left(k_{{\perp}} \rho_{s}\right)^3 \left(\frac{{\varTheta_{k_{{\perp}}}^\mp}}{\rho_{s} c_{\rm s}}\right). \end{equation}

Combining (3.6), (3.9) and (3.8), we obtain an estimate for the fluctuations of the Elsässer potentials

(3.10)\begin{equation} \frac{{\varTheta_{k_{{\perp}}}^\pm}}{\rho_{s} c_{\rm s}} \sim \left(\frac{\varepsilon_{W}}{\varOmega_{\rm i}}\right)^{1/3}\left(\frac{1+\bar{\tau}}{\bar{\tau}^2}\right)^{1/6} \left(k_{{\perp}} \rho_{s}\right)^{{-}5/3}, \end{equation}

and the associated spectra

(3.11)\begin{equation} E_\perp^\pm(k_{{\perp}}) \propto \left(\frac{1+ \bar{\tau}}{\bar{\tau}^2}\right)^{1/3} \left(k_{{\perp}} \rho_{s}\right)^{{-}7/3} \propto \left\{ \begin{array}{@{}ll@{}} \displaystyle k_{{\perp}}^{{-}5/3} , & k_{{\perp}} \rho_{\rm i} \ll 1, \\ \displaystyle k_{{\perp}}^{{-}7/3}, & k_{{\perp}} \rho_{\rm i} \gg 1. \end{array} \right. \end{equation}

These are the standard Kolmogorov-style scalings for both RMHD and KAW turbulence, respectively (see, e.g., Goldreich & Sridhar 1995; Cho & Lazarian 2004; Schekochihin et al. 2009), which hold even in the presence of finite electron inertia (only the magnetic energy exhibits a transition at $k_{\perp } d_{\rm e} \sim 1$; see (3.16)). However, simulations of forced RMHD turbulence have consistently shown scaling exponents closer to $-3/2$, whereas those of KAW turbulence appear to be closer to $-8/3$. The departure of both exponents from those derived here has been shown to arise due to intermittency effects (see, e.g., Boldyrev 2006; Chandran, Schekochihin & Mallet 2015; Mallet & Schekochihin 2017 in the former case, or Boldyrev & Perez 2012; Meyrand & Galtier 2013; Zhou, Liu & Loureiro 2023b, in the latter), which we have not attempted to account for here. Given that our aim is not a comprehensive theory of turbulence, but instead to highlight the effect of helicity conservation in such an environment, we consider agreement with either the derived spectra (3.11) or those corrected for intermittency to be sufficient evidence that the system is undergoing a constant-flux cascade. We will compare with a scaling exponent of $-3/2$ where necessary as it is well-motivated in balanced turbulence and supported by observations of imbalanced turbulence (Chen et al. 2020).

Perhaps the most important feature not captured by the scalings (3.10) and (3.11) is the imbalance, viz., the difference in amplitudes between each of the fields: we have assumed that ${\varTheta _{k_{\perp }}^+} \sim {\varTheta _{k_{\perp }}^-}$, which clearly will not be the case if ${\varepsilon ^{+}} \gg {\varepsilon ^{-}}$. One way of rectifying this would be to relax the first assumption in (3.8) and assume that the relevant nonlinear time for each field is the $\boldsymbol {E} \times \boldsymbol {B}$ advection rate associated with the counter-propagating field, which yields, via a straightforward generalisation of the argument given previously,

(3.12)\begin{equation} {\varTheta_{k_{{\perp}}}^\pm} \propto \left[\frac{\left(\varepsilon^\pm\right)^2}{\varepsilon^\mp}\right]^{1/3} \quad \Rightarrow \quad \frac{{\varTheta_{k_{{\perp}}}^+}}{{\varTheta_{k_{{\perp}}}^-}} \sim \frac{{\varepsilon^{+}}}{{\varepsilon^{-}}} = \text{const}. \end{equation}

This is the conclusion at which Lithwick et al. (2007) arrived in the context of strongly-imbalanced RMHD turbulence. Though there is numerical evidence to suggest that this approximately holds for the ratio of their associated energies (Beresnyak & Lazarian 2008, 2009a, 2010; Mallet & Schekochihin 2011; Schekochihin 2022), evidence for it holding throughout the inertial range is less clear. Previous work (see, e.g., Perez & Boldyrev 2010; Mallet & Schekochihin 2011) suggests that the stronger field typically has a steeper spectrum than the weaker one, although this difference in their slopes tends to decrease as numerical resolution is increased. Another potential issue following from (3.12) is that the ratio of the nonlinear times for each field must scale as $t_{\text {nl}}^+/t_{\text {nl}}^- \sim {\varepsilon ^{+}}/{\varepsilon ^{-}}$ (this follows from combining (3.9) and (3.12)). As pointed out by Lithwick et al. (2007), this has the counterintuitive implication that the weaker ${\varTheta _{k_{\perp }}^-}$ perturbation, which is advected by ${\varTheta _{k_{\perp }}^+}$ at a faster rate $(t_{\text {nl}}^-)^{-1}$, can nevertheless coherently advect ${\varTheta _{k_{\perp }}^+}$ at the slower rate $(t_{\text {nl}}^+)^{-1}$. Though they propose a disparity between the spatial and temporal coherences of each field as a possible explanation for this, such a state is hard to justify in general; for a possible alternative explanation, see § 9.6 of Schekochihin (2022). Finally, while the assumption that the counter-propagating field is the only source of nonlinear advection is guaranteed to be satisfied in the RMHD regime ($k_{\perp } \rho _{\rm i} \ll 1$), this is not obviously true at sub-ion scales ($k_{\perp } \rho _{\rm i} \gg 1$). The dispersive nature of KAWs makes nonlinear interactions between co-propagating perturbations (i.e. ${\varTheta ^\pm }$ with ${\varTheta ^\pm }$) possible at these scales, meaning that it is, in principle, possible to support a turbulent cascade with a single component ${\varTheta ^\pm }$ (Cho 2011; Kim & Cho 2015; Voitenko & De Keyser 2016). Given, however, that these subtleties take us beyond the main focus of this work, we will not engage further with them here, having highlighted our reasons for presenting only a balanced turbulence phenomenology when deriving the scaling predictions (3.11).

For completeness, we include the scalings of the two contributions to the free energy at RMHD scales (see (2.15)): the kinetic energy associated with the $\boldsymbol {E} \times \boldsymbol {B}$ flow $\propto |\boldsymbol {u}_{\perp }|^2$, and the energy contained in perpendicular magnetic field fluctuations $\propto |{\delta \! \boldsymbol {B}_{\!\perp }}|^2$. Defining their spectra, analogously to (3.7), as

(3.13)\begin{align} E_\perp^K(k_{{\perp}}) &= 2{\rm \pi} k_{{\perp}}\int_{-\infty}^\infty \,\text{d} k_{{\parallel}} \frac{e^2 n_{0{\rm e}}}{2 T_{0{\rm e}}} \rho_{s}^2 \left\langle |k_{{\perp}} {\phi_{\boldsymbol{k}}}|^2 \right\rangle, \end{align}

(3.14)\begin{align}E_\perp^B(k_{{\perp}}) &= 2{\rm \pi} k_{{\perp}}\int_{-\infty}^\infty \,\text{d} k_{{\parallel}} \frac{\left\langle |k_{{\perp}} {A_{{\parallel} \boldsymbol{k}}}|^2 \right\rangle}{8{\rm \pi}} , \end{align}

respectively, and comparing the definitions of the free-energy (2.14) and (3.1), we find the following scalings:

(3.15)\begin{equation} E_\perp^K(k_{{\perp}}) \sim \frac{\bar{\tau}^2}{1+\bar{\tau}} (k_{{\perp}} \rho_{s})^2 E_\perp^\pm(k_{{\perp}}) \propto \left\{ \begin{array}{@{}ll@{}} \displaystyle k_{{\perp}}^{{-}5/3} , & k_{{\perp}} \rho_{\rm i} \ll 1, \\ \displaystyle k_{{\perp}}^{{-}1/3}, & k_{{\perp}} \rho_{\rm i} \gg 1, \end{array} \right. \end{equation}

and

(3.16)\begin{equation} E_\perp^B(k_{{\perp}}) \sim \frac{1}{1 + k_{{\perp}}^2 d_{\rm e}^2} E_\perp^\pm(k_{{\perp}}) \propto \left\{ \begin{array}{@{}ll@{}} \displaystyle k_{{\perp}}^{{-}5/3} , & k_{{\perp}} \ll \rho_{\rm i}^{{-}1}, d_{\rm e}^{{-}1}, \\ \displaystyle k_{{\perp}}^{{-}7/3} , & \rho_{\rm i}^{{-}1} \ll k_{{\perp}} \ll d_{\rm e}^{{-}1}, \\ \displaystyle k_{{\perp}}^{{-}11/3} , & d_{\rm e}^{{-}1} \ll k_{{\perp}} \ll \rho_{\rm i}^{{-}1}, \\ \displaystyle k_{{\perp}}^{{-}13/3}, & k_{{\perp}} \gg \rho_{\rm i}^{{-}1}, d_{\rm e}^{{-}1}. \end{array} \right. \end{equation}

To test these predictions, we consider the simulations in table 1 labelled ‘high resolution’. The first two simulations, CF-res (‘constant flux’) and HB-res (‘helicity barrier’), differ only in their values of the electron inertial length (or, equivalently, the electron beta), having $d_{\rm e} = \rho _{\rm i}$ and $d_{\rm e} =\rho _{\rm i}/2$, respectively. The third is a simulation of RMHD, included for comparison. As we explain in the following, both simulations CF-res and RMHD are expected to saturate via a constant-flux cascade: this is what we indeed find, with their spectra, plotted in figure 2, showing good agreement with the theoretical predictions (3.10), (3.15), and (3.16), up to the aforementioned corrections due to intermittency. This serves as numerical confirmation of the turbulence phenomenology presented above across the RMHD, ERMHD, ultra-low-beta and sub-$d_{\rm e}$ regimes (see, e.g., Schekochihin et al. 2009; Meyrand & Galtier 2010; Loureiro & Boldyrev 2018), at least without detailed consideration of the difficulties relating to imbalance discussed previously. We note, however, that figure 2(a) appears to support the idea that both ${\varTheta ^\pm }$ fields have the same scaling despite their difference in amplitude, as in (3.8), although with the caveat that they do not have the same dissipation scale.

Figure 2. One-dimensional perpendicular energy spectra for the ‘high-resolution’ simulations in table 1: CF (blue), HB (red) and RMHD (green). Solid and dashed lines correspond to (a) $E_\perp ^+(k_{\perp })$ and $E_\perp ^-(k_{\perp })$ or (b) $E_\perp ^K(k_{\perp })$ and $E_\perp ^B(k_{\perp })$, respectively. The insets panels show the local scaling exponents $\alpha _{(\dots )} = \text {d} \log E^{(\dots )}_\perp /\text {d} \log k_{\perp }$ for each spectrum. The inertial-range scalings (3.10), (3.15) and (3.16) are shown by black dotted lines (with $-5/3$ scalings replaced with $-3/2$, as discussed following (3.11)). Note that the axis of the (scale-invariant) RMHD simulation has been rescaled for comparison with the other two cases. Both simulations CF (constant flux) and RMHD are expected to saturate via a constant-flux cascade, and show good agreement with the predicted scalings, whereas simulation HB (helicity barrier) exhibits entirely different scalings.

3.2 Effect of helicity conservation

As we noted in § 2.3, the fact that the generalised helicity (2.16) must be conserved by the system places an additional constraint on the dynamics that we did not account for in the theory presented in the preceding section. Let us rectify this now. It will be instructive to consider, for the moment, a theory of constant-flux turbulence in $\phi$ and ${A_{\parallel }}$ variables, rather than the ${\varTheta ^\pm }$ ones more appropriate for imbalanced turbulence that were used in § 3.1.

To begin, we can, from the definition (2.14), estimate the energy flux as

(3.17)\begin{equation} \frac{1}{n_{0{\rm e}} T_{0{\rm e}}} \frac{\text{d} W}{\text{d} t} \sim t_{\text{nl}}^{{-}1} \left(\frac{1+ \bar{\tau}}{\bar{\tau}^2}\right) \left(\frac{e{\phi}_{k_{{\perp}}}}{T_{0{\rm e}}}\right)^2 \sim \varepsilon_{W} = \text{const.}, \end{equation}

where $t_{\text {nl}}$ is some nonlinear time associated with the cascade; we will remain agnostic to exactly what this is. Note that in writing (3.17), we have assumed that the energy can be adequately represented by the first two terms in (2.14), viz., those associated with the electrostatic potential fluctuations. This is justified by the fact that the contributions to the free energy from $\phi$ and ${A_{\parallel }}$ must at least be in equipartition in order for the turbulence to be imbalanced (recall (2.12) and (3.1)). Analogously to (3.17), we can estimate the helicity flux from (2.16) as

(3.18)\begin{equation} \frac{1}{n_{0{\rm e}} T_{0{\rm e}}} \frac{\text{d} H}{\text{d} t} \sim t_{\text{nl}}^{{-}1} \frac{v_{\rm A}}{c_{\rm s}} \left(\frac{1 + k_{{\perp}}^2 d_{\rm e}^2}{\bar{\tau}}\right) \left(\frac{e\phi_{k_{{\perp}}}}{T_{0{\rm e}}}\right) \left(\frac{A_{{\parallel} k_{{\perp}}}}{\rho_{s} B_0}\right) \cos\alpha_{k_{{\perp}}} \sim \varepsilon_{H} = \text{const.}, \end{equation}

where $\cos \alpha _{k_{\perp }}$ is (the cosine of) some perpendicular-wavenumber-dependent phase angle between the $\phi$ and ${A_{\parallel }}$ fluctuations. This could be formally defined in terms of the Fourier components ${\phi _{\boldsymbol {k}}}$ and ${A_{\parallel \boldsymbol {k}}}$ as

(3.19)\begin{equation} \cos \alpha_{k_{{\perp}}} = \frac{\displaystyle\int_{-\infty}^\infty \,\text{d} k_{{\parallel}} \mathrm{Re}\left\langle {\phi_{\boldsymbol{k}}} ({A_{{\parallel} \boldsymbol{k}}})^*\right\rangle}{\displaystyle \left(\int_{-\infty}^\infty \,\text{d} k_{{\parallel}} \left\langle | {\phi_{\boldsymbol{k}}}|^2\right\rangle\right)^{1/2} \left(\int_{-\infty}^\infty \,\text{d} k_{{\parallel}} \left\langle|{A_{{\parallel} \boldsymbol{k}}}|^2 \right\rangle\right)^{1/2}}, \end{equation}

wherein ‘$*$’ denotes the complex conjugate, and the brackets denote an ensemble average, as previously. Taking the ratio of (3.18) to (3.17), using equipartition between the energies to relate the amplitudes of $\phi$ and ${A_{\parallel }}$, and recalling the definition of the injection imbalance $\sigma _\varepsilon = |\varepsilon _{H}|/\varepsilon _{W}$, it is straightforward to show that

(3.20)\begin{equation} \varepsilon_{H} \sim \varepsilon_{W} \left( \frac{v_{\text{ph}}}{v_{\rm A}} \right)^{{-}1} \cos \alpha_{k_{{\perp}}} \quad \Rightarrow \quad \cos\alpha_{k_{{\perp}}} \sim \sigma_{\varepsilon} \frac{v_{\text{ph}}}{v_{\rm A}}. \end{equation}

Finally, given that $0 < | \cos \alpha _{k_{\perp }} | \leqslant 1$, it follows directly from (3.20) that

(3.21)\begin{equation} \sigma_{\varepsilon} \frac{v_{\text{ph}}(k_{{\perp}})}{v_{\rm A}}\lesssim 1. \end{equation}

This inequality must be satisfied everywhere in perpendicular-wavenumber space in order for the system to support a constant flux of helicity from large to small scales; if it is anywhere violated, then the assumption of constant flux breaks down. It is important to clarify that this latter statement does not only apply to the helicity flux, but also that of the free energy: if the system is not able to simultaneously cascade both invariants via constant flux, it is unable to support a constant flux of either free energy or helicity individually. This means that the nature of the turbulence is fundamentally different depending on whether or not the inequality (3.21) is satisfied; if it is, we obtain the constant-flux-type turbulence discussed in § 3.1; if it is not, then the system will inevitably form a helicity barrier (Meyrand et al. 2021; Squire et al. 2022), the exact dynamics of which we shall discuss in detail in § 3.3. For now, the phrase ‘helicity barrier’ can simply serve as a placeholder term for the turbulent state that occurs in the absence of a constant-flux cascade within isothermal KREHM.

Due to the perpendicular-wavenumber dependence of the phase velocity (2.9), there are a number of regimes of (3.21) to consider.

1. (i) FLR-MHD ($\beta _{\rm e} \gg m_{\rm e}/m_{\rm i}$): In this regime, the phase velocity is a strictly increasing function of perpendicular wavenumber (this is the solid black in figure 1; see also (2.11)), meaning that, no matter the injection imbalance, (3.21) will always be violated at sufficiently small scales. The physical reason for this is clear from the first expression in (3.20). Given a constant energy flux $\varepsilon _{W}$, the increase in the phase velocity causes the helicity flux $\varepsilon _{H}$ to decrease with increasing $k_{\perp }$, which cannot be compensated for by, e.g., increasing $\cos \alpha _{k_{\perp }}$ (i.e. further aligning the fluctuations), since the latter is bounded from above. This means that the assumption of constant flux cannot be satisfied, which, in this system, manifests itself through the formation of a helicity barrier. In a finite simulation domain, however, the resolution may not be sufficient to allow the phase velocity to increase to a value at which it would violate (3.21). As such, if the dissipation scale of the turbulence $k_{\perp }^\text {diss}$ is less than some critical perpendicular wavenumber

   (3.22)\begin{equation} k_{{\perp}}^\text{diss} \lesssim k_{{\perp}}^\text{crit}, \quad k_{{\perp}}^\text{crit} \rho_{\rm i} = \frac{1}{\sigma_{\varepsilon}} \left(\frac{2}{1+Z/\tau}\right)^{1/2}, \end{equation}
   then the system is able to support a constant-flux cascade. This prediction is tested numerically in § 3.4.

2. (ii) Isothermal KREHM ($\beta _{\rm e} \gtrsim m_{\rm e}/m_{\rm i}$): In this intermediate regime, the phase velocity increases with perpendicular wavenumber, but eventually reaches a constant value at small scales (see figure 1). Using the fact that the maximum value of the phase velocity is given by the second expression in (2.10), we can rewrite (3.21) as

   (3.23)\begin{equation} \beta_{\rm e} \lesssim \beta_{\rm e}^{\text{crit}}, \quad \beta_{\rm e}^{\text{crit}} = \frac{2Z}{1+\tau/Z} \frac{m_{\rm e}}{m_{\rm i}} \frac{1}{\sigma_{\varepsilon}^2}. \end{equation}
   There is thus a critical value of the plasma beta below which a constant-flux cascade can occur, allowing the free energy and helicity to reach the smallest scales. In particular, (3.23) predicts a constant-flux solution is always possible for $\beta _{\rm e} \lesssim m_{\rm e}/m_{\rm i}$: it is below this value of $\beta _{\rm e}$ that the phase velocity begins to decrease with perpendicular wavenumber, and we find ourselves in the ultra-low beta regime.

3. (iii) Ultra-low beta ($\beta _{\rm e} \ll m_{\rm e}/m_{\rm i}$): This is the opposite of regime (i), in that the phase velocity is a strictly decreasing function of perpendicular wavenumber (see (2.11)). This means that, since $\sigma _{\varepsilon } \leqslant 1$, the inequality (3.21) will always be satisfied, regardless of the injection imbalance. This once again follows from the first expression in (3.20): the decrease in the phase velocity would cause $\varepsilon _{H}$ to increase with increasing $k_{\perp }$, but, unlike in regime (i), this always can be compensated for by decreasing $\cos \alpha _{k_{\perp }}$, and so a constant-flux cascade is always allowed. This is the phenomenon of ‘dynamic phase alignment’ (Loureiro & Boldyrev 2018; Milanese et al. 2020), wherein $\phi$ and ${A_{\parallel }}$ fluctuations become increasingly misaligned at small scales in order to maintain a constant flux of both free energy and helicity. The relationship between the helicity barrier and dynamic phase alignment is discussed further in § 3.5.

It is worth noting that while the inequality (3.21) was derived here using the above heuristic scaling arguments, it can be put on a rigorous footing. In regimes (i) and (ii), where the phase velocity (2.9) is an increasing function of perpendicular wavenumber, it is possible to prove, following an argument identical to that used in Meyrand et al. (2021) (see also Alexakis & Biferale 2018), that the perpendicular fluxes of free energy $\varPi _W(k_{\perp })$ and helicity $\varPi _H(k_{\perp })$ must satisfy the inequality $|\varPi _H(k_{\perp })|/|\varPi _W(k_{\perp })| \leqslant v_{\rm A}/v_{\text {ph}}(k_{\perp })$ in order for a constant-flux solution to exist (the presence of the phase velocity arises from the perpendicular-wavenumber dependence of the scale-by-scale helicity; see (3.2)). This inequality implies (3.21), since $|\varPi _H(k_{\perp })|/|\varPi _W(k_{\perp })| = \sigma _{\varepsilon }$ in a constant-flux cascade, by definition. One cannot apply the same argument in regime (iii), where the phase velocity is a decreasing function of perpendicular wavenumber: this is inconsequential, however, because (3.21) will always be satisfied in this regime, and so a constant-flux solution is always possible. Given, then, that the inequality (3.21) can only be violated in the regimes where rigorous proof exists, we regard it as a stronger condition than might be implied by the heuristic arguments that were initially used to derive it. The resultant prediction, and subsequent numerical confirmation, of the existence of the critical beta (3.23) is the central result of this paper. However, in order to be able to test this prediction, one must have an understanding of the dynamics associated with the helicity barrier, to which we now turn our focus. The impatient reader, or one already familiar with the helicity barrier, may wish to skip ahead to § 3.4, working backwards where further clarification is required.

3.3 The helicity barrier

To explore the dynamics of the helicity barrier within isothermal KREHM, we compare simulations CF and HB from the set ‘comparison’ in table 1. Importantly, these simulations lie on either side of the critical beta (3.23), with values of $d_{\rm e}$ that differ by a factor of two, but have parameters that are otherwise identical. The injection imbalance of $\sigma _{\varepsilon } = 0.8$ used for both simulations corresponds to a critical beta $\beta _{\rm e}^{\text {crit}} (m_{\rm e}/m_{\rm i})^{-1} = 1.56$: simulation CF, which saturates via constant flux, has $\beta _{\rm e} / \beta _{\rm e}^{\text {crit}} = 0.64$, whereas simulation HB, which forms a helicity barrier, has $\beta _{\rm e}/\beta _{\rm e}^{\text {crit}} = 2.56$. The choice of these specific values, however, is inconsequential: simulations on either side of $\beta _{\rm e}^{\text {crit}}$ display almost identical behaviour regardless of the choice of $\sigma _{\varepsilon }$. We have chosen these particular simulations because they provide illustrative examples of the behaviour of the system with and without the helicity barrier, despite their very similar parameters.

3.3.1 Energy fluxes

In figure 3, we plot the spectral energy fluxes $\varPi ^\pm (k_{\perp })$ associated with the potentials ${\varTheta ^\pm }$ from these two simulations, calculated by summing the contributions of the nonlinear transfers above and below the particular $k_{\perp }$ of interest. Simulation CF, with $\beta _{\rm e}/\beta _{\rm e}^{\text {crit}} < 1$, behaves as expected: the fluxes of both ${\varTheta ^+}$ and ${\varTheta ^-}$ are approximately stationary (in time), constant (as a function of perpendicular wavenumber), and equal to their injected values (3.5). Note, from figure 3(a), that the entirety of the injected free energy $\varepsilon _{W} = {\varepsilon ^{+}} + {\varepsilon ^{-}}$ is carried by the turbulence to small scales, where it is then dissipated by the perpendicular hyperdissipation (2.18). This provides a posteriori justification of the constant-flux assumption made in § 3.1.

Figure 3. Time evolution of the spectral energy fluxes $\varPi ^\pm (k_{\perp })$ computed directly from the nonlinear terms in (2.7) and (2.8), normalised to the total energy flux $\varepsilon _{W}$. Simulations CF and HB from table 1 are shown in panels (a) and (b), respectively. The colours indicate the value of time corresponding to a given line, whereas the solid black lines correspond to the average value of each flux over the last $20\,\%$ of the simulation time. The horizontal dashed lines indicate the values of the flux (3.5) expected if the system is able to maintain a constant flux; we have included a line corresponding to ${\varepsilon ^{-}}$ in the upper panel of (b) for ease of comparison. It is clear that the total flux reaching small scales in the presence of the helicity barrier is significantly smaller than in the constant flux case (note that, in both cases, the decrease in the flux at small scales is due to the presence of perpendicular hyperdissipation).

The simulation HB, with $\beta _{\rm e} /\beta _{\rm e}^{\text {crit}} > 1$, however, displays very different behaviour. While the flux of ${\varTheta ^-}$ remains approximately stationary and constant at ${\varepsilon ^{-}}$ (more so, in fact, than simulation CF), the flux of ${\varTheta ^+}$ in the inertial range is a decreasing function of both time and perpendicular wavenumber, whereas the largest scales display rapid fluctuations. Of particular note, readily apparent from figure 3(b), is the fact that only a small fraction of the injected energy is cascaded to small scales, remaining limited to $\approx 2{\varepsilon ^{-}} = (1-\sigma _{\varepsilon })\varepsilon _{W}$, irrespective of the turbulent amplitudes, which are much larger than in simulation CF (see figure 5). This is the helicity barrier: the breakdown of the constant-flux assumption, due to a violation of the inequality (3.21), causes the system to form a ‘barrier’ that prevents all but the balanced portion of the injected energy from cascading to small perpendicular scales.

3.3.2 Dissipation

The remaining energy must, therefore, find another route to thermalisation, which it does so by accessing small parallel scales. In figure 4, we plot the free energy $W$ and its parallel and perpendicular (hyper)dissipation rates, denoted $D_{\parallel }$ and $D_{\perp }$, respectively, as a function of time for both simulations. It is clear that $D_{\parallel } \ll D_{\perp }$ at all times in simulation CF, allowing the free energy to quickly saturate on perpendicular dissipation, as expected from a system undergoing a constant-flux cascade. Conversely, in simulation HB, the ratio of the parallel dissipation rate to the total dissipation rate, viz.,

(3.24)\begin{equation} R_{\text{diss}} = \frac{D_{{\parallel}}}{D_{{\parallel}} + D_{{\perp}}}, \end{equation}

is an increasing function of time; we plot this explicitly in the lower panel of figure 4. Both this ratio (henceforth termed the dissipation ratio) and the free energy eventually saturate at late times when the large-scale turbulent amplitudes reach sufficiently high levels that the energy removed by parallel dissipation (at small parallel scales) can balance the fraction of the injected energy that is unable to cascade (Meyrand et al. 2021). This saturation is, of course, unphysical, since it breaks the assumption of anisotropy ($k_{\parallel } \ll k_{\perp }$) used to derive the isothermal KREHM system, and depends on the details of the parallel dissipation (e.g. the specific value of $\nu _z$; see Meyrand et al. 2021). As such, the only dynamics relevant to real physical systems are those that occur before this saturation, a period of time that we shall henceforth refer to as the pseudostationary phase.

Figure 4. The free energy and its dissipation rates as a function of time for simulations CF (blue) and HB (red). The top panel shows the free energy, whereas middle panel shows the parallel and perpendicular dissipation rates, $D_{\parallel }$ and $D_{\perp }$, in dashed and solid lines, respectively. The bottom panel shows the dissipation ratio (3.24). It is clear that $R_{\text {diss}} \ll 1$ for the constant-flux case, while $R_{\text {diss}}$ grows slowly to $\approx 1$ for the helicity-barrier one.

Figure 5. Time evolution of the perpendicular spectra $E_\perp ^\pm (k_{\perp })$ for simulations CF and HB in table 1, shown in panels (a) and (b), respectively. The colours indicate the value of time corresponding to a given line, whereas the dotted black lines indicate approximate scalings of the spectra. The inset panels show the scaling exponents $\alpha ^\pm = \text {d} \log E_\perp ^\pm /\text {d} \log k_{\perp }$ for each spectrum. In the helicity barrier case, the $E_\perp ^+(k_{\perp })$ spectrum clearly forms a spectral break that moves towards large scales in time.

3.3.3 Perpendicular energy spectra

What is the impact of these dynamics on the measured perpendicular energy spectra? In figure 5, we plot the time evolution of the $E_\perp ^\pm (k_{\perp })$ spectra (3.7) for both simulations, with their high-resolution counterparts (see table 1) plotted in figure 2(a). As expected, simulation CF agrees well with the constant-flux cascade predictions of § 3.1: both $E_\perp ^\pm (k_{\perp })$ have approximately the same $\sim k_{\perp }^{-3/2}$ slope in the inertial range, differing by only their outer-scale amplitudes. In simulation HB, the weaker spectrum $E_\perp ^-(k_{\perp })$ behaves similarly, quickly saturating with a $\sim k_{\perp }^{-3/2}$ slope. However, given that the ${\varTheta ^+}$ flux is a decreasing function of perpendicular wavenumber, the associated $E_\perp ^+(k_{\perp })$ spectrum cannot reach such a stationary state. Instead, during the start of the pseudostationary phase, it forms a spectral break around $k_{\perp } \rho _{\rm i} \sim 1$, with the location of this spectral break then migrating to larger scales. The slope below the break is consistently measured to be $\sim k_{\perp }^{-4}$, across all simulations in table 1 that exhibited a helicity barrier. The higher-resolution simulation HB-res (see table 1) plotted in figure 2 also shows a spectral flattening at smaller scales $k_{\perp } \rho _{\rm i} \gtrsim 1$, consistent with observations of the spectral ‘transition range’ in near-Sun solar-wind plasmas (Bowen et al. 2020a; Duan et al. 2021; Bowen et al. 2024). In addition, the reduction of the ${\varTheta ^+}$ flux to small scales (in comparison with simulation CF; see figure 3) causes the outer-scale energy, and thus the normalised cross-helicity $\sigma _c$, to increase in time (while both the energy injection rate $\varepsilon _{W}$ and the injection imbalance $\sigma _{\varepsilon }$ remain constant). As discussed previously, this will continue until the large-scale amplitudes reach sufficient levels that saturation can occur on parallel dissipation. Both Meyrand et al. (2021) and Squire et al. (2023) found that the position of the break is correlated with $\sigma _c$ finding that, approximately, its location in perpendicular-wavenumber space evolved as $k_{\perp } \rho _{\rm i} \sim (1-\sigma _c)^{1/4}$. We have not attempted to verify such a scaling here due to the relatively small inertial range at this resolution, and the fact that said scaling may be complicated by the finite $d_{\rm e}$ effects present in these simulations. Nevertheless, although a dynamical theory that explains these features remains the focus of ongoing work, the existence of the break, as well as the steep spectral scaling below it, are both persistent features of helicity-barrier-mediated turbulence.

It should be clear from the previous discussion that the helicity barrier state is dramatically different from that associated with a constant-flux cascade, viz., the nature of the turbulence is fundamentally changed depending on whether or not the inequality (3.21) is satisfied. The surprising element of this is the fact that these two states can lie so close to one another in parameter space: we recall that the electron inertial scales for the simulations that we have been considering in this section only differ by a factor of two, being $d_{\rm e} = \rho _{\rm i}$ for the constant-flux case and $d_{\rm e} = \rho _{\rm i}/2$ for the helicity-barrier case. Indeed, the real-space snapshots of the turbulence shown in figure 6 are completely different, despite this small difference in physical parameters.

Figure 6. Real-space snapshots of the $\boldsymbol {E} \times \boldsymbol {B}$ flow $\boldsymbol {u}_{\perp }$ (see (2.3)) for the simulations CF-res (left) and HB-res (right). The colours indicate the magnitude of $\boldsymbol {u}_{\perp }$ relative to its (spatial) root-mean-square value, whereas the coordinate directions are as shown. Although the structure of the turbulence in simulation CF-res is typical of a constant-flux cascade (though note the significant small-scale plasmoid activity due to finite $d_{\rm e}$ effects; see Zhou et al. 2023b), this is not the case for simulation HB-res: the majority of the energy resides in large-scale structures because it is prevented from cascading to small scales by the helicity barrier. The dramatic difference between these two cases is made even more surprising by the fact that the simulations differ only in their value of the electron inertial length, having $d_{\rm e} = \rho _{\rm i}$ and ${d_{\rm e} = \rho _{\rm i}/2}$, respectively.

3.4 Breaking the barrier

To briefly summarise the findings of the previous section, the helicity barrier has three key features: (i) it only allows the balanced fraction ($\approx 2{\varepsilon ^{-}}$) of the free energy to cascade to small scales (§ 3.3.1); (ii) the remainder of the free energy remains at large perpendicular scales where it dissipates on small parallel ones, meaning that the ratio of the parallel dissipation to the total dissipation $R_{\text {diss}}$ (see (3.24)) is an increasing function of time during the pseudostationary phase (§ 3.3.2); and (iii) the spectrum of ${\varTheta ^+}$ displays a sharp spectral break, with an approximate $\sim k_{\perp }^{-4}$ scaling below it, which widens (moves towards larger scales) over time (§ 3.3.3).

We now wish to test the predictions of § 3.2 across a wide range of $\beta _{\rm e}$ and $\sigma _{\varepsilon }$. This requires the ability to efficiently determine whether or not a helicity barrier has formed in a given simulation. While feature (i) is the clearest measure of the helicity barrier, calculating the spectral energy fluxes $\varPi ^\pm (k_{\perp })$ is computationally expensive, making it unfeasible to use across a large simulation set. Similarly, feature (iii) is not a reliable measure of helicity-barrier formation at early times because it requires the spectral slopes both above and below the break to be sufficiently well resolved, which is not always possible at lower resolutions. As such, we choose to exploit the fact that the dissipation ratio (3.24) is an increasing function of time during the pseudostationary phase, viz., simulations that, on average, have $\text {d}R_{\text {diss}}/\text {d} t \approx 0$ will be undergoing a constant-flux cascade, while those with $\text {d} R_{\text {diss}}/\text {d} t \gtrsim 0$ will have formed a helicity barrier. This is by no means a unique measure of helicity-barrier formation, but will prove sufficient and appropriate for our purposes here. We emphasise that while the dissipation ratio provides a useful measure when applied to our simulations, it should not be viewed as a diagnostic to be measured with spacecraft data, wherein other features of the helicity barrier (such as energy fluxes, spectral slopes, or implied heating rates) would be more appropriate. Indeed, it is only a useful measure within subsidiary limits of gyrokinetics such as isothermal KREHM, since artificial parallel dissipation must be added to allow saturation in the absence of other dissipative mechanisms on small parallel scales at $k_{\perp } \rho _{\rm i} \lesssim 1$ (e.g. ICW heating of ions around $k_{\parallel } d_{\rm e} \sim 1$, which lies outside the gyrokinetic approximation).

To explicitly demonstrate the utility of using $R_{\text {diss}}$ as a diagnostic, we consider the series of simulations with $\sigma _{\varepsilon } = 0.8$ from the set labelled ‘beta scan’ in table 1, in which $\beta _{\rm e}/\beta _{\rm e}^{\text {crit}}$ is varied between $0.53$ and $4.00$. In figure 7(a), we plot $R_{\text {diss}}$ as a function of time for each simulation. It is clear the simulations are split between those that are approximately steady in time (constant flux), and those that increase with time (helicity barrier). This is illustrated in more detail in figure 7(b), in which we plot the two-dimensional spectrum of the dissipation $D_\text {2D}(k_{\perp }, k_z)$ for the same set of simulations. Those whose dissipation ratio is approximately constant in time dissipate almost all of their energy at small perpendicular scales (right-hand side of the plot), whereas those whose dissipation ratio is growing have more dissipation at small parallel scales and $k_{\perp } \rho _{\rm i} < 1$ (top left of the plot), a clear signature of the helicity barrier.² The behaviour of the simulation with $\beta /\beta _{\rm e}^{\text {crit}} = 1$, somewhat unsurprisingly, lies intermediate between these two states: the dissipation ratio initially grows in time during a short pseudostationary phase (figure 7a), before achieving saturation with similar levels of parallel and perpendicular dissipation (figure 7b).

Figure 7. Measures of the dissipation associated with the set of eight simulations with $\sigma _{\varepsilon } = 0.8$ from the set labelled ‘beta scan’ in table 1. (a) The dissipation ratio (3.24) plotted as a function of time. The colours indicate the value of $\beta /\beta _{\rm e}^{\text {crit}}$ for each simulation. (b) The two-dimensional spectrum of the dissipation in the $(k_{\perp }, k_z)$ plane, averaged over the last $20\,\%$ of the simulation time, and normalised to the energy injection rate $\varepsilon _{W}$, with amplitude as indicated by the colourbar.

It is clear that the rate of change of the dissipation ratio gives a clear measure of helicity barrier formation within this set of simulations. However, we require a criterion for helicity-barrier formation that will apply across the full range of imbalances considered in table 1. For this, we choose to consider the non-dimensionalised rate of change of the dissipation ratio given by $(\text {d} R_{\text {diss}}/\text {d} t) (\varepsilon _{H}/H)^{-1}$. The normalising factor $\varepsilon _{H}/H$, the ratio of the helicity injection rate to the helicity itself, is effectively a measure of the timescale over which the imbalanced portion of the energy (compare (3.1) and (3.2)) builds up on scales $k_{\perp } \rho _{\rm i} \lesssim 1$ (recall, from § 3.3.1, that only the balanced portion $\approx 2{\varepsilon ^{-}}$ is allowed to cascade to small scales). Given that $\varepsilon _{H} = \sigma _{\varepsilon } \varepsilon _{W}$, this normalisation accounts for both variations in the overall energy injection rate $\varepsilon _{W}$, as well as the fact that systems with lower imbalances will typically have smaller $\text {d} R_{\text {diss}}/\text {d} t$ even when a helicity barrier has formed (the lower imbalance means that it takes longer for sufficient energy to build up on the largest perpendicular scales and begin dissipating on small parallel scales). Using the above set of simulations, we find that a time-averaged value of

(3.25)\begin{equation} \left[\frac{\text{d} R_{\text{diss}}}{\text{d} t } \left(\frac{\varepsilon_{H}}{H}\right)^{{-}1}\right]^\text{crit} = 0.25 \end{equation}

is reasonable to distinguish between the constant-flux and helicity-barrier regimes. The time average is performed over the last 80 % of the simulation time in order to exclude the initial transient phase that occurs in the constant-flux simulations.³ We note that although simulations with $\beta /\beta _{\rm e}^{\text {crit}} \approx 1$ could be classified differently depending on the exact value chosen in (3.25), the classification of simulations with $\beta /\beta _{\rm e}^{\text {crit}}$ very different from unity is robust to such choices. As such, we will henceforth use the numerical value (3.25) as our criterion for determining helicity-barrier formation, as applied to the simulations considered in table 1. Let us now test the two predictions made in § 3.2: that the formation of the helicity barrier depends on: (i) having a plasma beta above the critical value (3.23); and (ii) resolving the critical perpendicular wavenumber (3.22) the dissipation range.

For the former, we consider the set of simulations labelled ‘beta scan’ in table 1, in which the injection imbalance $\sigma _{\varepsilon }$ is varied between $0.1$ and $0.99$, whereas $\beta _{\rm e} (m_{\rm e}/m_{\rm i})^{-1}$ is varied between $0.83$ and $625$ (corresponding to a variation in $\beta _{\rm e}/\beta _{\rm e}^{\text {crit}}$ between $0.03$ and $25$). Note that as $\sigma _{\varepsilon }$ is varied, the forcing is also modified in such a way as to keep the amplitude of the stronger field ${\varTheta ^+}$ approximately constant according to the first expression in (3.12), i.e. the forcing is scaled to keep $({\varepsilon ^{+}})^2/{\varepsilon ^{-}} \propto (1-\sigma _{\varepsilon })/(1+\sigma _{\varepsilon })^2 = \text {const}$. This was done in order to minimise the possible effect of the outer-scale forcing on helicity barrier formation. The results of this scan are plotted in figure 8, which displays excellent agreement with (3.23) over multiple orders of magnitude in $\beta _{\rm e} (m_{\rm e}/m_{\rm i})^{-1}$. The inset panel demonstrates that the criterion that we have applied to determine helicity barrier formation is very well satisfied: there is a rapid increase in $(\text {d} R_{\text {diss}}/\text {d} t) (\varepsilon _{H}/H)^{-1}$ around $\beta _{\rm e}/\beta _{\rm e}^{\text {crit}} = 1$ across all sets of simulations, consistent with this being the location in parameter space where the helicity barrier forms. We thus confirm the prediction (3.23) of the critical line in $(\sigma _{\varepsilon },\beta _{\rm e})$ parameter space above which a helicity barrier will always form; the implications of this result are discussed in § 4.

Figure 8. Data from a two-dimensional parameter scan of injection imbalance $\sigma _{\varepsilon }$ versus the (normalised) electron plasma beta $\beta _{\rm e} (m_{\rm e}/m_{\rm i})^{-1}$. Each circle corresponds to a simulation from the set ‘beta scan’ in table 1, with filled/open circles indicating the presence/absence of the helicity barrier. The solid black line corresponds to the value of $\sigma _{\varepsilon }$ below which, for a given $\beta _{\rm e}$, the helicity barrier should not form (cf. (3.23)). The shaded region below the line thus corresponds to saturation via constant flux, whereas above it a helicity barrier forms. The inset plot shows the time-averaged average value of the rate-of-change of the dissipation ratio (3.24) normalised to $\varepsilon _{H}/H$ for each set of simulations, as indicated by the colour, with the horizontal axis rescaled by the critical beta (3.23). The horizontal dotted line therein corresponds to the value (3.25) above which the helicity barrier is determined to have formed.

For prediction (ii), that a helicity barrier forms only when the critical perpendicular wavenumber (3.22) is well-resolved, we consider the set of simulations labelled ‘resolution scan’ in table 1. These all lie in the FLR-MHD limit (finite $\rho _{\rm i}$, $d_{\rm e} \rightarrow 0$). In these simulations, as well as all others considered in this paper, the adaptive dissipation (see § 2.4) is implemented in such a way as to ensure that the dissipation scale $k_{\perp }^\text {diss}$ is approximately half that of the maximum wavenumber in the simulation $k_{\perp }^\text {max}$, irrespective of forcing amplitude and injection imbalance. This means that (3.22) can be rewritten as a condition on $k_{\perp }^\text {max}$, viz.,

(3.26)\begin{equation} k_{{\perp}}^\text{max} \rho_{\rm i} \lesssim 2 k_{{\perp}}^\text{crit} \rho_{\rm i} = \frac{2}{\sigma_{\varepsilon}} \left(\frac{2}{1+Z/\tau}\right)^{1/2}, \end{equation}

where $k_{\perp }^\text {crit}$ is as defined in (3.22). We expect no helicity barrier to be present if the inequality in (3.26) is satisfied. To confirm this, we varied the injection imbalance $\sigma _{\varepsilon }$ across a different set of resolutions from $64^3$ to $256^3$, as in table 1. The results of this scan are shown in figure 9, which show good agreement with (3.26): deviations from this prediction are due to the difficulty of ensuring adequate separation between regions of parallel and perpendicular dissipation at lower resolutions. That being said, given that real physical systems are not limited by resolution, we do not consider exact agreement with (3.26) necessary. Rather, these results provide clear evidence supporting the general principles of helicity barrier formation outlined previously.

Figure 9. Data from a two-dimensional parameter scan of injection imbalance $\sigma _{\varepsilon }$ versus the maximum wavenumber $k_{\perp }^\text {max} \rho _{\rm i}$ set by the numerical resolution. Each circle corresponds to a simulation from the set ‘resolution scan’ in table 1, with filled/open circles indicating the presence/absence of the helicity barrier. The solid black line corresponds to the value of $\sigma _{\varepsilon }$ below which, for a given $k_{\perp }^\text {max} \rho _{\rm i}$, the helicity barrier should not form (see (3.26)). As in figure 8, the shaded region below the line thus corresponds to saturation via constant flux, whereas above it a helicity barrier forms. The inset plot shows the time-averaged average value of the rate-of-change of the dissipation ratio (3.24) normalised to $\varepsilon _{H}/H$ for each set of simulations, as indicated by the colour, with the horizontal axis rescaled by $2k_{\perp }^\text {crit} \rho _{\rm i}$. The horizontal dotted line therein corresponds to the value (3.25) above which the helicity barrier is determined to have formed.

3.5 A comment on dynamic phase alignment

Before discussing the implications of these findings, let us briefly comment on the relationship between the helicity barrier and the concept of dynamic phase alignment (Loureiro & Boldyrev 2018; Milanese et al. 2020). As discussed in § 3.2, dynamic phase alignment refers to the phenomenon whereby, in the ultra-low-beta regime $\beta _{\rm e} \ll m_{\rm e}/m_{\rm i}$, the fluctuations of the electrostatic potential $\phi$ and parallel magnetic vector potential ${A_{\parallel }}$ become increasingly misaligned at small scales in order to maintain a constant flux of both free energy and helicity. The principal finding of Milanese et al. (2020) was that this alignment is directly manifest in the phase angle (3.19), in that it becomes a decreasing function of perpendicular wavenumber, as per the theoretical prediction (3.20). As we discussed in § 3.2, the helicity barrier arises as a consequence of such an alignment becoming impossible, breaking the constant-flux solution.

To illustrate this, we plot, in figure 10, the phase angle (3.19) for the set of three simulations labelled ‘comparison’ in table 1: the first two are those examined in detail in § 3.3, whereas the third is in the ultra-low-beta regime considered by Loureiro & Boldyrev (2018) and Milanese et al. (2020). It is clear that in the latter case, $\cos \alpha _{k_{\perp }}$ decreases with perpendicular wavenumber, as demonstrated by Milanese et al. (2020), whereas the former cases both have $\cos \alpha _{k_{\perp }} \approx 1$: simulation CF is sufficiently close to the critical boundary predicted by (3.21) (or, alternatively, (3.20)) that it must be as highly aligned in order to maintain a constant flux, whereas simulation HB has already passed beyond this threshold and formed a helicity barrier, with the fluctuations remaining maximally aligned on the largest scales. We have also plotted (dashed lines) the theoretical scalings (3.20) for the phase angle in figure 10; the fact that this curve exceeds unity for simulation HB implies the formation of a helicity barrier, i.e. it violates the criterion (3.20), which is what we indeed find dynamically. In some sense, the helicity barrier can be viewed as the ‘opposite’ of dynamic phase alignment, in that it is the state that occurs when the system cannot sufficiently align the $\phi$ and ${A_{\parallel }}$ fluctuations. Indeed, heuristically, these cases can easily be distinguished by recalling the definitions of the free energy (3.1) and helicity (3.2). Although it always remains possible to decrease the difference between the ${\varTheta ^\pm }$ energies at small scales to compensate for the decrease in the phase velocity at $\beta /\beta _{\rm e}^{\text {crit}} \lesssim 1$, the opposite is not always true at $\beta /\beta _{\rm e}^{\text {crit}} \gtrsim 1$: at a certain perpendicular wavenumber, the difference in energies would need to be greater than the sum, and so the constant-flux solution must break down. There is, however, a crucial difference between these two cases: whereas dynamic phase alignment modifies the dynamics at small perpendicular scales in order to extend the constant-flux solution to lower values of $\beta _{\rm e}$, the helicity barrier explicitly violates the constant-flux solution, placing it into an entirely different class of turbulent dynamics to that which can be obtained via dynamic phase alignment. Importantly, the effects of the helicity barrier are thus manifest even on the largest perpendicular scales accessible to the system, having dramatic implications for the turbulent heating properties thereof.

Figure 10. The phase angle (3.19) between fluctuations of the electrostatic potential $\phi$ and parallel magnetic vector potential ${A_{\parallel }}$, as a function of perpendicular wavenumber, for the three simulations labelled ‘comparison’ in table 1 (solid lines). Simulations CF and HB both have $\rho _{\rm i} = 0.1 L$ and finite $d_{\rm e}$, whereas simulation ULB has $d_{\rm e} =0.1 L$ but $\rho _{\rm i} \rightarrow 0$. The dashed lines show the theoretical scaling (3.20), whereas the dotted black line shows the expected scaling in the ultra-low-beta regime. We do not expect exact agreement between (3.20) and (3.19) because the former is a result derived using a ratio of fluxes (see § 3.2) whereas for the latter we are plotting a ratio of amplitudes (or, equivalently, energies).

4 Summary and discussion

The findings presented in this paper serve as a detailed illustration of the sensitivity of imbalanced turbulence to small changes in its characteristic physical parameters. Starting from equations derived in the low-beta asymptotic limit of gyrokinetics (isothermal KREHM; see § 2), we showed that the requirement for the simultaneous conservation of both free energy and helicity implies the presence of a critical electron beta (see § 3.2) given by

(4.1)\begin{equation} \beta_{\rm e}^{\text{crit}} = \frac{2Z}{1+\tau/Z} \frac{m_{\rm e}}{m_{\rm i}} \frac{1}{\sigma_{\varepsilon}^2}, \end{equation}

where $Z$ is the ion charge, $\tau$ is the equilibrium-temperature-ratio between ions and electrons, with $m_{\rm e}/m_{\rm i}$ their mass ratio, and $\sigma _{\varepsilon }$ is the ratio of the injection rates of cross-helicity and energy at large scales (injection imbalance). This theoretical prediction is well-supported by simulations of imbalanced turbulence in isothermal KREHM across a wide range of $\beta _{\rm e}$ and $\sigma _{\varepsilon }$ (see figure 8). Systems situated close to either side of (4.1) in parameter space exhibit dramatically different turbulent dynamics, as evident from even the most cursory glance at the real-space turbulence snapshots shown in figure 6 (showing two simulations with values of $\beta _{\rm e}$ that differ by a factor of four, but are otherwise identical). Specifically, in systems with $\beta _{\rm e}$ below (4.1), the free energy injected at the largest perpendicular scales is able to undergo a constant-flux, Alfvénic cascade to smaller scales $k_{\perp } \rho _{\rm e} \lesssim 1$ where it dissipates, and in so doing depositing the majority of the turbulent free energy into electron heating. Systems with $\beta _{\rm e}$ exceeding (4.1), on the other hand, are unable to support a constant-flux solution, inevitably forming a helicity barrier. This prevents all but the balanced portion of the injected free energy ($\approx 2 {\varepsilon ^{-}}$) from cascading past the ion Larmor scale $k_{\perp } \rho _{\rm i} \sim 1$, resulting in a build-up of turbulent free energy at larger scales $k_{\perp } \rho _{\rm i} \lesssim 1$. Fluctuations on these scales eventually form fine parallel structures which, in our system, dissipate via parallel hyperviscosity: in a more comprehensive system, they would excite high-frequency ion-cyclotron waves (ICWs), leading to perpendicular ion heating (Squire et al. 2022, 2023). Thus, the constant-flux and helicity-barrier states, demarcated by the critical beta (4.1), offer entirely different propositions for turbulent heating: the majority of the injected turbulent free energy is converted into electron heating, in the former case, or ion heating, in the latter.

Assuming that the existence of the helicity barrier, or otherwise, plays a central role in determining turbulent heating, these results have clear implications for observations of imbalanced Alfvénic turbulence. For highly imbalanced plasmas with a modest plasma beta, we would expect to observe dominant ion heating, and for the spectral slopes of the electromagnetic fields to exhibit a steep ‘transition range’ scaling $\sim k_{\perp }^{-4}$ bracketing $k_{\perp } \rho _{\rm i} \sim 1$, a distinctive feature of helicity-barrier-mediated turbulence (Meyrand et al. 2021; Squire et al. 2022, 2023). Conversely, for plasmas with either a small imbalance (at large solar radii) or a very low plasma beta ($\beta _{\rm e}$ approaching $m_{\rm e}/m_{\rm i}$), we would expect to observe more electron heating, and for the steep transition range scaling to be replaced by the much shallower $\sim k_{\perp }^{-7/3}$ scaling predicted from KAW turbulence (see § 3.1). Given that much of the solar wind typically has $m_{\rm e}/m_{\rm i} \ll \beta _{\rm e} \lesssim 1$ (Bruno & Carbone 2005), (4.1) would suggest that it should usually display signatures of helicity-barrier-mediated turbulence. This is consistent with observations: ions are typically hotter than electrons (Cranmer et al. 2009), with significant power in ICWs around $k_{\parallel } \rho _{\rm i} \sim 1$ (Huang et al. 2020; Bowen et al. 2020b), whereas spectra observed by PSP usually show the aforementioned steep transition range scalings (Bowen et al. 2020a; Duan et al. 2021; Bowen et al. 2024).

Nevertheless, there may be regions of the solar wind in which the helicity barrier does not operate. Equation (4.1) predicts that a constant-flux cascade is always possible if $\beta _{\rm e} \leqslant m_{\rm e}/m_{\rm i}$, irrespective of imbalance, as well as the fact that the helicity barrier will not form at sufficiently low imbalance. This means that the helicity barrier is unlikely to operate in regions of the solar wind that are strongly magnetised or have low imbalance, or indeed some combination of the two. That being said, the values of $\sigma _{\varepsilon }$ below which the helicity barrier will not operate are quite low. Indeed, our simulations showed helicity-barrier formation at values of the injection imbalance as low as $\sigma _{\varepsilon } = 0.2$ (see figure 8), with lower values limited by numerical resolution, rather than by some fundamental constraint on the dynamics (see § 3.2 or the last paragraph in § 3.4). It is worth noting in this context, however, that the injection imbalance $\sigma _{\varepsilon }$ is not the same as the often-measured normalised cross-helicity $\sigma _c$, so careful work is required in order to extract an exact correspondence between values of $\sigma _{\varepsilon }$ used in simulations and the dynamics observed in astrophysical systems.

4.1 Limitations of the isothermal approximation

An important limitation of this current study concerns the possible role of electron kinetics, which we have here neglected in order to isolate the effects of including finite electron inertia. The inclusion thereof has two important physical consequences that could potentially alter the results above: (i) it introduces electron Landau damping as another channel for turbulent electron heating, allowing energy to be transferred to small scales in velocity space (in addition to in wavenumber space); (ii) it modifies the conservation of helicity viz., instead of (2.16) being everywhere conserved, it is able to be injected and/or removed by higher-order velocity moments of the kinetic distribution function. Both of these effects are most significant around $k_{\perp } d_{\rm e} \sim 1$, and so we expect little change for values of $\beta _{\rm e}$ significantly above (4.1). For those close to (4.1), however, the dynamics could be significantly modified, e.g. the helicity barrier may not form because the helicity ceases to be globally conserved or, should the helicity barrier persist, electron Landau damping could provide another source of dissipation for fluctuations at $k_{\perp } \rho _{\rm i} \lesssim 1$. This latter effect would be significant, as the amplitude reached by fluctuations on these scales plays a central role in determining the amount of perpendicular ion heating by ICWs (Squire et al. 2022, 2023), and thus the effective fraction of the injected energy earmarked for electron heating. In either case, the reintroduction of electron kinetics should have the effect of shifting the critical beta (4.1) towards larger values (i.e. moving the solid black curve in figure 8 to the right), giving rise to more electron, and less ion, heating. Confirming these predictions is the subject of ongoing work. In perhaps a preview of what is to come, recent simulations of balanced turbulence by Zhou et al. (2023a) suggest that the advection of energy in velocity space is the dominant mechanism of nonlinear energy transfer at $k_{\perp } d_{\rm e} \sim 1$, which becomes the primary route to dissipation.

There are, of course, other mechanisms for turbulent heating that may be playing a role. A compressive energy cascade, although unable to exchange energy with Alfvénic motions (Schekochihin et al. 2019), is able to cause parallel ion heating (Schekochihin et al. 2009), and will break helicity conservation around $k_{\perp } \rho _{\rm i} \sim 1$, potentially arresting the breakdown of the constant flux solution associated with helicity-barrier formation. In order for either of these effects to be significant, however, the powers in the compressive and Alfvénic cascades would likely need to be comparable, which is generally not observed in the solar wind (Bruno & Carbone 2013; Chen 2016; Chen et al. 2020). Other heating mechanisms that rely on the presence of large-amplitude fluctuations, such as stochastic heating (Chandran et al. 2010) or sub-ion-Larmor-radius KAW turbulence (Arzamasskiy et al. 2019; Isenberg & Vasquez 2019), may also play a role.

Even with these considerations taken into account, the results of this paper represent a robust prediction of two fundamentally different types of turbulence: one of a constant-flux cascade of energy to small scales, the other involving the build-up and dissipation of energy at the largest scales due to the helicity-barrier mechanism. The critical beta (4.1) thus marks a boundary between two dramatically different regimes of turbulent heating, each with clear observational signatures that can further constrain the possible physical processes at play in imbalanced Alfvénic turbulence. This physics could play an important role in many magnetised astrophysical environments where the same symmetry exists (i.e. having imbalance), such as accretion discs, the intracluster medium, and the solar wind context discussed in detail in this work. Therefore, advancing our understanding of the dynamics related to the helicity barrier and the resultant mechanisms of turbulent heating will carry important implications for plasma turbulence across a diverse array of astrophysical contexts.