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. (Reference Meyrand, Squire, Schekochihin and Dorland2021), 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 (Reference Loureiro and Boldyrev2018) and Milanese et al. (Reference Milanese, Loureiro, Daschner and Boldyrev2020). 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}$.Footnote ¹ 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. (Reference Meyrand, Squire, Schekochihin and Dorland2021). 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 Reference Elsässer1950), 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 Reference Zocco and Schekochihin2011; Loureiro et al. Reference Loureiro, Dorland, Fazendeiro, Kanekar, Mallet, Vilelas and Zocco2016; Adkins et al. Reference Adkins, Schekochihin, Ivanov and Roach2022):

(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. Reference Schekochihin, Cowley, Dorland, Hammett, Howes, Quataert and Tatsuno2009).

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. Reference Schekochihin, Cowley, Dorland, Hammett, Howes, Plunk, Quataert and Tatsuno2008, Reference Schekochihin, Cowley, Dorland, Hammett, Howes, Quataert and Tatsuno2009), as in hydrodynamic turbulence (see, e.g., Alexakis & Biferale Reference Alexakis and Biferale2018, 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. Reference Schekochihin, Cowley, Dorland, Hammett, Howes, Quataert and Tatsuno2009). 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. Reference Teaca, Verma, Knaepen and Carati2009) 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 (Reference Williamson1980) 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. Reference Meyrand, Squire, Mallet and Chandran2024). 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. Reference Meyrand, Squire, Schekochihin and Dorland2021); 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.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 Reference Velli, Grappin and Mangeney1989). 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 (Reference Kolmogorov1941) 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 Reference Davidson2013). 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 Reference Schekochihin2022). 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 Reference Goldreich and Sridhar1995; Cho & Lazarian Reference Cho and Lazarian2004; Schekochihin et al. Reference Schekochihin, Cowley, Dorland, Hammett, Howes, Quataert and Tatsuno2009), 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 Reference Boldyrev2006; Chandran, Schekochihin & Mallet Reference Chandran, Schekochihin and Mallet2015; Mallet & Schekochihin Reference Mallet and Schekochihin2017 in the former case, or Boldyrev & Perez Reference Boldyrev and Perez2012; Meyrand & Galtier Reference Meyrand and Galtier2013; Zhou, Liu & Loureiro Reference Zhou, Liu and Loureiro2023b, 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. Reference Chen, Bale, Bonnell, Borovikov, Bowen, Burgess, Case, Chandran, Dudok de Wit and Goetz2020).

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. (Reference Lithwick, Goldreich and Sridhar2007) 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 Reference Beresnyak and Lazarian2008, Reference Beresnyak and Lazarian2009a, Reference Beresnyak and Lazarian2010; Mallet & Schekochihin Reference Mallet and Schekochihin2011; Schekochihin Reference Schekochihin2022), evidence for it holding throughout the inertial range is less clear. Previous work (see, e.g., Perez & Boldyrev Reference Perez and Boldyrev2010; Mallet & Schekochihin Reference Mallet and Schekochihin2011) 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. (Reference Lithwick, Goldreich and Sridhar2007), 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 (Reference Schekochihin2022). 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 Reference Cho2011; Kim & Cho Reference Kim and Cho2015; Voitenko & De Keyser Reference Voitenko and De Keyser2016). 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. Reference Schekochihin, Cowley, Dorland, Hammett, Howes, Quataert and Tatsuno2009; Meyrand & Galtier Reference Meyrand and Galtier2010; Loureiro & Boldyrev Reference Loureiro and Boldyrev2018), 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. Reference Meyrand, Squire, Schekochihin and Dorland2021; Squire et al. Reference Squire, Meyrand, Kunz, Arzamasskiy, Schekochihin and Quataert2022), 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 Reference Loureiro and Boldyrev2018; Milanese et al. Reference Milanese, Loureiro, Daschner and Boldyrev2020), 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. (Reference Meyrand, Squire, Schekochihin and Dorland2021) (see also Alexakis & Biferale Reference Alexakis and Biferale2018), 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. Reference Meyrand, Squire, Schekochihin and Dorland2021). 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. Reference Meyrand, Squire, Schekochihin and Dorland2021). 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. Reference Bowen, Mallet, Bale, Bonnell, Case, Chandran, Chasapis, Chen, Duan and Dudok de Wit2020a; Duan et al. Reference Duan, He, Bowen, Woodham, Wang, Chen, Mallet and Bale2021; Bowen et al. Reference Bowen, Bale, Chandran, Chasapis, Chen, Dudok de Wit, Mallet, Meyrand and Squire2024). 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. (Reference Meyrand, Squire, Schekochihin and Dorland2021) and Squire et al. (Reference Squire, Meyrand and Kunz2023) 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. Reference Zhou, Liu and Loureiro2023b), 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.Footnote ² 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.Footnote ³ 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 Reference Loureiro and Boldyrev2018; Milanese et al. Reference Milanese, Loureiro, Daschner and Boldyrev2020). 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. (Reference Milanese, Loureiro, Daschner and Boldyrev2020) 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 (Reference Loureiro and Boldyrev2018) and Milanese et al. (Reference Milanese, Loureiro, Daschner and Boldyrev2020). It is clear that in the latter case, $\cos \alpha _{k_{\perp }}$ decreases with perpendicular wavenumber, as demonstrated by Milanese et al. (Reference Milanese, Loureiro, Daschner and Boldyrev2020), 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.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. Reference Squire, Meyrand, Kunz, Arzamasskiy, Schekochihin and Quataert2022, Reference Squire, Meyrand and Kunz2023), 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. (Reference Zhou, Liu and Loureiro2023a) 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. Reference Schekochihin, Kawazura and Barnes2019), is able to cause parallel ion heating (Schekochihin et al. Reference Schekochihin, Cowley, Dorland, Hammett, Howes, Quataert and Tatsuno2009), 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 Reference Bruno and Carbone2013; Chen Reference Chen2016; Chen et al. Reference Chen, Bale, Bonnell, Borovikov, Bowen, Burgess, Case, Chandran, Dudok de Wit and Goetz2020). Other heating mechanisms that rely on the presence of large-amplitude fluctuations, such as stochastic heating (Chandran et al. Reference Chandran, Li, Rogers, Quataert and Germaschewski2010) or sub-ion-Larmor-radius KAW turbulence (Arzamasskiy et al. Reference Arzamasskiy, Kunz, Chandran and Quataert2019; Isenberg & Vasquez Reference Isenberg and Vasquez2019), 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.