2.1 The two-field gyrofluid model

The model (which retains no dissipation process) is written as equations for the electron gyrocentre number density $N_e$ and the parallel component of the magnetic potential $A_{\|}$. When electron inertia is neglected, it takes the form

(2.1)\begin{gather} \partial_t N_e +[\varphi,N_e]-[B_z,N_e]+\frac{2}{\beta_e}\boldsymbol{\nabla}_{\|} \varDelta_{{\perp}} A_{\|}=0, \end{gather}

(2.2)\begin{gather}\partial_t A_{\|} + \boldsymbol{\nabla}_{\|} (\varphi-N_e-B_z)=0. \end{gather}

Here, $\varDelta _{\perp } = \partial _{xx} + \partial _{yy}$ is the Laplacian in the plane transverse to the ambient field and $[f,g]= \partial _x f \partial _y g-\partial _y f \partial _x g$ the canonical bracket of two scalar functions $f$ and $g$. Furthermore, $\varGamma _n$ denotes the (non-local) operator $\varGamma _n(-\tau \varDelta _{\perp })$ which in Fourier space reduces to the multiplication by the function $\varGamma _n(\tau k_{\perp }^2)$, defined by $\varGamma _n(x) = {\rm I}_n(x) e^{-x}$, where ${\rm I}_n$ is the modified Bessel function of first type of order $n$. For a scalar function $f$, the parallel gradient operator $\boldsymbol {\nabla }_{\parallel }$ is defined by $\boldsymbol {\nabla }_{\parallel } f=-[A_{\parallel } , f]+ \partial _z f$.

The equations are written in a non-dimensional form, using the following quantities as unities: $\varOmega _i^{-1}$ for time, $\rho _s$ for the space coordinates (and thus $\rho _s^{-1}$ for the wavenumber components), $B_0$ for the parallel magnetic fluctuations $B_z$, $n_0$ for the electron gyrocentre density $N_e$, $T_e/e$ for the electric potential $\varphi$ and $B_0\rho _s$ for the parallel magnetic potential $A_{\|}$.

Introducing the operators $L_1 = {2}/{\beta _e} +(1+2\tau )(\varGamma _0-\varGamma _1)$, $L_2 = 1 + (1-\varGamma _0)/\tau - \varGamma _0 +\varGamma _1$, $L_3 = (1-\varGamma _0)/\tau$ and $L_4 = 1-\varGamma _0+\varGamma _1$, the parallel magnetic fluctuations are given by $B_z = M_1 \varphi$, with $M_1 = L_1^{-1} L_2$, and the electrostatic potential by $\varphi =-M_2^{-1} N_e$, where $M_2 = L_3 + L_4L_1^{-1}L_2$ is positive definite. Thus, $B_z$ and $\varphi$ can both be expressed in terms of $N_e$.

At the linear level, the phase velocity $v_{{\rm ph}}= \omega /k_z$ is given by

(2.3)\begin{equation} v_{{\rm ph}}^2\equiv\left(\frac{\omega}{k_z}\right)^2 = \frac{2}{\beta_e} k_{{\perp}}^2 \frac{1 -{\hat M}_1 + {\hat M}_2} {{\hat M}_2}, \end{equation}

where the caret refers to the Fourier symbol of the operator. The associated operator $V_{{\rm ph}}$ is given by

(2.4)\begin{equation} V_{{\rm ph}}= s (-\varDelta_{{\perp}})^{1/2} (1 -M_1 + M_2)^{1/2}M_2^{{-}1/2}, \end{equation}

where $s= (2/\beta _e)^{1/2}$ is the equilibrium Alfvén velocity in sound speed units. In physical space, the eigenmodes, which can be viewed as generalized Elsasser potentials, are given by $\mu ^\pm =\varLambda \varphi \pm s A_{\|}$, where $\varLambda = (-\varDelta _{\perp })^{-1} (1+M_2-M_1)^{1/2} M_2^{1/2}$.

In the absence of dissipation processes (needed in the turbulent regime), the system (2.1)–(2.2) preserves the energy ${\mathcal {E}}$ and the generalized cross-helicity ${\mathcal {C}}$:

(2.5)\begin{gather} {\mathcal{E}} = \frac{1}{2} \int \left( \frac{2}{\beta_e} |\boldsymbol{\nabla}_{{\perp}} A_{\|} |^2 + \frac{4\delta^2}{\beta_e^2}|\varDelta_{{\perp}} A_{\|} |^2 - N_e(\varphi -N_e-B_z) \right) \,{\rm d}^3 {x}, \end{gather}

(2.6)\begin{gather}{\mathcal{C}} ={-}\int N_e A_{\|}\,{\rm d}^3 {x}. \end{gather}

In terms of $\mu ^\pm$, the two invariants are rewritten ($D=(-\varDelta )^{1/2}$)

(2.7)\begin{gather} {\mathcal{E}} = \frac{1}{4} \int \left\{(D \mu^+)^2 + (D \mu^-)^2 \right\} \,{\rm d}^3x, \end{gather}

(2.8)\begin{gather}{\mathcal{C}} = \frac{1}{4} \int V_{{\rm ph}}^{{-}1}\{D \mu^+)^2 - (D \mu^-)^2\} \,{\rm d}^3x. \end{gather}

2.2 Numerical set-up

Equations (2.1)–(2.2) with $\beta _{\rm e}=2$ and $\tau = 1$ were solved in an isotropic periodic box of size $L$, using a spectral method for the space variables and a third-order Runge–Kutta scheme for the time stepping. An eight-order hyper-dissipation operator $\nu _{\perp }(\varDelta _{\perp })^4 + \nu _z \partial ^8_{z}$ acting on $N_e$ and on $A_{\|}$ is supplemented in the equations for the corresponding quantities. For convenience, the parameters $\nu _{\perp }$ and $\nu _{\|}$ are referred to as transverse (or perpendicular) and parallel viscosities. The system is driven by freezing the amplitude of the modes whose transverse wavenumber stands in the first spectral shell and the parallel one corresponds to $k_z (L / 2{\rm \pi} )= \pm 1$, as in Siggia & Patterson (Reference Siggia and Patterson1978). A main interest of this procedure is that it permits simulations with a high imbalance and a moderate nonlinearity parameter. In such a regime, the injection rates of energy and of GCH are not prescribed, in contrast with the situations where the system is driven by an external random force or by a negative damping (Meyrand et al. Reference Meyrand, Squire, Schekochihin and Dorland2021).

3.1 Spectral transition zone and arrest of the GCH cascade

In addition to the imbalance given by the ratio $I={\mathcal {E}}^+ / {\mathcal {E}}^+$ of the energies of the forward- and backward-propagating waves, another quantity governing the dynamics is provided by the nonlinearity parameter $\chi$, defined as the ratio of the inverse nonlinear time to the Alfvén frequency. At early times, when the maintained modes are dominant, a simple phenomenological argument (Miloshevich et al. Reference Miloshevich, Laveder, Passot and Sulem2020) leads to $\chi \approx (k_{\perp 0}/k_{z 0}) |B_{\perp }(k_{\perp 0})| = |B_{\perp }(k_{\perp 0})|$, where $B_{\perp }$ is the transverse magnetic fluctuations (measured in units of the ambient field), and the index zero refers to the maintained modes for which $k_{\perp 0} = k_{z 0}$ in the (anisotropic) units of the gyrofluid model. At later times, in the isotropic simulation box, the root-mean-square value of the magnetic fluctuation is used as an estimation of $\chi$. In this paper, we concentrate on the case where the nonlinearity parameter keeps moderate values (initially $\chi \leq 0.5$), unless otherwise stated. Several simulations, referred to as R ₁ , ... R ₈,  were performed in the above framework. Their parameters are specified in Table 1.

Table 1. Parameters of the runs, together with the values $\chi _0$ and $\chi _{\rm final}$ of the nonlinear parameter at the beginning and at the end of the simulation, estimated by the square root of the magnetic energy at this time. The colour code refers to the colour of the lines in figures 7–11.

We first illustrate our main results on a large-resolution simulation (run $R_1$) that displays extended inertial ranges. At the end of this simulation, the system did not yet reach a stationary state. Nevertheless, as mentioned in § 3.2, the observed dynamics is similar to that found in a lower-resolution simulation (run $R_2$), carried out until saturation is reached. For the $R_1$ simulation characterized by an initial imbalance $I = 100$ and an initial value $\chi _0 = 0.47$ of the nonlinearity parameter, performed in a computational box of size $L=11 \times 2 {\rm \pi}$, figure 1 displays the transverse spectra (integrated over $k_z$) $E^\pm (k_{\perp })$ of the forward- and backward-propagating waves (figure 1a), and $E_{B_{\perp }}(k_{\perp })$, $E_{B_z}(k_{\perp })$ of the transverse and parallel magnetic fluctuations (figure 1b), averaged over the time interval [14 656, 14 760] corresponding to $1.5$ Alfvén crossing times (see figure 3 for a characterization of this interval relatively to the full simulation). The spectra $E^+(k_{\perp })$ and $E_{B_{\perp }}(k_{\perp })$ display a transition zone at scales slightly larger than the sonic Larmor radius (corresponding to $k_{\perp } = 1$), with a spectral index close to $-3$, located between a sub-ion range with a spectral index $-2.7$ and a MHD range whose extension is strongly limited by the numerical constraints which make difficult the simulation in a box that would ideally cover a more extended range at large scale, while also retaining the small-scale dynamics. The transition zone is seen to extend up to the wavenumber where the spectrum $E_{B_z}(k_{\perp })$ changes from a quasi-flat range to a decay range displaying the same spectral exponent as $E_{B_{\perp }}(k_{\perp })$. The large-scale $B_z$ fluctuations are much smaller than those of $B_{\perp }$ as expected when only AWs are excited at large scale (see e.g. Kobayashi et al. Reference Kobayashi, Sahraoui, Passot, Laveder, Sulem, Huang, Henri and Smets2017).

Figure 1. Transverse energy spectrum $E^\pm$ of the forward- and backward-propagating waves (a) and energy spectrum of the transverse and parallel magnetic fluctuations $E_{B_{\perp }}$ and $E_{B_z}$ (b), averaged on three Alfvén crossing times $L/v_A$, for run $R_1$.

Figure 2 shows the perpendicular fluxes of energy and GCH, averaged over the same time interval as the spectra displayed in figure 1. One observes that the energy flux is quasi-constant, even if the temporal fluctuations are comparable to the mean value (not shown) at scales larger than the transition one. Differently, the GCH flux undergoes a conspicuous decay at the transition range, a consequence of the helicity barrier whose existence was pointed out by Meyrand et al. (Reference Meyrand, Squire, Schekochihin and Dorland2021) in similar simulations (although in the asymptotic regime $\beta _e \to 0$), performed in the large-$\chi$ regime where the energy flux to small transverse scales is also significantly inhibited by the barrier.

Figure 2. Energy and GCH fluxes in the transverse direction for run $R_1$ averaged in the same time interval as in figure 1.

Because of numerical constraints, run $R_1$ was performed in three steps, associated with increased resolutions ($N=420^3$, $560^3$ and $784^3$ mesh points) and decreased (equal) parallel and perpendicular viscosities. This procedure permitted us to analyse the sensitivity of the parallel and perpendicular dissipations of energy (figure 3a) and GCH (figure 3b). We observe that the perpendicular dissipations are essentially not affected by the transition from one step to the following one, as expected in a usual turbulence cascade: reduction of the viscosity does not affect the dissipation and only shifts the dissipation range to smaller scales. In contrast, we observe that this is not the case for the parallel dissipation of energy and even more for GCH, which both decrease with the viscosity, confirming the absence of a standard turbulent cascade in the direction of the ambient field, at least in the regime of small or moderate nonlinear parameters considered in the present study.

Figure 3. Parallel energy (a) and GCH (b) dissipations ${\mathcal {D}}_z({\mathcal {E}})$ and ${\mathcal {D}}_z({\mathcal {C}})$, respectively, for run $R_1$. The run was started using equal perpendicular and parallel viscosity coefficients $\nu =5 \times 10^{-8}$ using a resolution of $420^3$ collocation points (red line); a first restart was then performed with $\nu = 5 \times 10^{-9}$ and $560^3$ collocation points (green line), and another one with $\nu = 3 \times 10^{-10}$ and $784^3$ collocation points (blue line). The insets show the corresponding perpendicular energy and GCH dissipations ${\mathcal {D}}_{\perp }({\mathcal {E}})$ and ${\mathcal {D}}_{\perp }({\mathcal {C}})$. Due to lack of space, the time labels are not specified for the insets that cover the same time interval as the main graphs.

Figure 4 displays two-dimensional colour plots of the time-averaged parallel (figure 4a) and perpendicular (figure 4b) GCH dissipations, defined by $D_z^C(k_{\perp }, k_z) = \nu _z k_z^8 E_C(k_{\perp }, k_z)$ and $D_{\perp }^C(k_{\perp }, k_z) = \nu _{\perp } k_{\perp }^8 E_C(k_{\perp }, k_z)$, respectively, where both the viscosities $\nu _{\perp }$ and $\nu _z$ are here equal to $3 \times 10^{-10}$ (and where $E_C(k_{\perp }, k_z)$ stands for the GCH power spectrum, after angle-averaging in the transverse plane). In addition to the usual small-scale dissipation in both parallel and perpendicular directions, one clearly observes the existence of a significant parallel dissipation of GCH, which is larger near a transverse wavenumber close to $k_{\perp }=1$ than at the smallest scales.

Figure 4. Parallel (a) and perpendicular (b) dissipation of GCH in the $(k_{\perp },k_z)$ plane (logarithmic coordinates) for run $R_1$.

These observations can be understood by considering the contributions of the $+$ and $-$ waves, separately. Figure 5 shows the two-dimensional colour plots of the time-averaged parallel dissipations $D_z^+(k_{\perp }, k_z) = \nu _z k_z^8 E^+(k_{\perp }, k_z)$ and $D_z^-(k_{\perp }, k_z) = \nu _z k_z^8 E^-(k_{\perp }, k_z)$ of the energies associated with the forward- and backward-propagating waves. A dissipation spot is visible near $k_{\perp } = 1$ in the case of the $+$ waves only. In contrast, comparable dissipations of the forward- and backward-propagating waves take place at small scales, in both the parallel and perpendicular directions. Indeed, the (weak) perpendicular transfer $\epsilon ^{-}$ of $E^{-}$ is almost not affected by the GCH barrier, while $\epsilon ^{+}$ is reduced past the barrier to a value close to $\epsilon ^{-}$ since at these scales imbalance has almost disappeared (see also Meyrand et al. (Reference Meyrand, Squire, Schekochihin and Dorland2021), in a case where the helicity barrier affects the energy transfer more significantly than in the present simulation). In order to reach a stationary state, the transfer in the $z$ direction and the resulting parallel dissipation thus have to become sizable near the ion scale for the $+$ waves. At small scale, the dissipations of the two types of waves are comparable and thus almost cancel out in the GCH dissipation, making the parallel dissipation of GCH to be dominant at transverse wavenumbers close to $k_{\perp } =1$, as observed in figure 4.

Figure 5. Parallel dissipation of the energies of the $+$ (a) and $-$ (b) waves in the $(k_{\perp },k_z)$ plane for run $R_1$.

A synthetic view of the dissipations in terms of the transverse wavenumber is obtained when integrating on $k_z$ the dissipations $D_z^+$, $D_z^-$, $D_{\perp }^+$, $D_{\perp }^-$ associated with the energy of the backward- and forward-propagating waves, and the analogous dissipations $D_z^C$, $D_{\perp }^C$ for the GCH. The resulting functions $d_z^+(k_{\perp })$, $d_z^-(k_{\perp })$, $d_{\perp }^+(k_{\perp })$, $d_{\perp }^-(k_{\perp })$, $d_z^C(k_{\perp })$, $d_{\perp }^C(k_{\perp })$ are displayed in figures 6(a) (dissipation of energies) and 6(b) (dissipation of GCH), under the same time averaging as in figures 4 and 5. The perpendicular energy dissipation is concentrated at small scales, with comparable magnitudes for both directions of propagation. Differently, for the parallel dissipation, while the energy of the backward-propagating waves is dissipated at small scales only, the dissipation of the forward-propagating waves, and thus the GCH, also displays a maximum near the transition zone.

Figure 6. (a) Parallel (red) and perpendicular (blue) dissipations of the energy of the $+$ wave, and parallel (green) and perpendicular (magenta) dissipations of the energy of the $-$ wave, integrated on $k_z$, versus $k_{\perp }$ for run $R_1$. (b) Parallel (red) and perpendicular (blue) dissipations of GCH, integrated on $k_z$, versus $k_{\perp }$ for the same run.

It turns out that the effect of the helicity barrier is enhanced when the large-scale imbalance is increased, resulting in a steeper transition zone, in agreement with solar wind observations reported in Huang et al. (Reference Huang, Sahraoui, Andrés, Hadid, Yuan, He, Zhao, Galtier, Zhang and Deng2021). This is exemplified in figure 7(a) which displays the time-averaged spectrum $E^+(k_{\perp })$ for three runs with different initial imbalance: $I=10$ (run $R_3$), $I=100$ (run $R_2$) and $I=1400$ (run $R_8$). Time averages are performed in the intervals [20 000, 25 000], [30 000, 40 000] and [7064, 7096], respectively. Compared with run $R_1$, they have a similar value of $\chi$ but are integrated in a box of size $L=5.5\times 2{\rm \pi}$, with resolutions of $280^3$ (for $I = 10$, 100) and $720^3$ (for $I = 1400$) mesh points. When considering the ratio $E^+(k_{\perp })/E^-(k_{\perp })$ (figure 7b), we observe that it becomes close to unity for $k_{\perp } \gtrsim 2$, whatever the value of the imbalance at large scales.

Figure 7. Influence of the imbalance on the perpendicular $E^+(k_{\perp })$ spectrum (a) and on the ratio $E^+(k_{\perp })/ E^-(k_{\perp })$ (b): $I = 10$ (run $R_3$), $I=100$ (run $R_2$), $I=1400$ (run $R_8$).

3.2 Temporal evolution of the dissipations and of the ideal invariants

While at zero imbalance the energy quickly reaches its saturation value (with nevertheless significant temporal oscillations), when the imbalance is increased, the time needed to reach saturation can become very large. This is shown in figure 8 which also indicates that, for a small initial value of $\chi$, the typical saturation time increases with $I$ (runs $R_2$, $R_3$, $R_4$).Footnote ¹

Figure 8. Time evolution (averaged over three Alfvén crossing times) of the energy (a) and GCH (b), normalized to their initial values, for runs $R_2$, $R_3$ and $R_4$, with similar values of the nonlinear parameter $\chi$ (whose time evolution is displayed in the inset in b) and different imbalances: $I=100$ (red), $I=10$ (green), $I=1$ (orange), respectively; a fourth run ($R_5$) with $I=100$ and a larger $\chi$ is also displayed in magenta.

Interestingly, for run $R_2$ with an energy imbalance $I=100$, the injection imbalance, given by the GCH to energy injection ratio, is relatively small, with a value close to $0.1$ at saturation, as seen in figure 9(a). Here, the injection rates are given by the slope of the cumulative injections, estimated as the sums of the energy or the GCH and of their respective cumulative dissipations. Figure 9(b) displays the energy and GCH fluxes for run $R_2$ in the steady state. They look qualitatively similar to those of run $R_1$, and in particular confirm the persistence of the GCH barrier in the saturated regime, in spite of the small value of the injection imbalance. In contrast, an energy barrier is hardly visible in this regime.

Figure 9. (a) Cumulative injections of energy (solid line) and GCH (dashed line). (b) Energy (solid line) and GCH (dashed line) fluxes of run $R_2$, averaged on the time interval $[30\,000, 40\,000]$.

Furthermore, as illustrated in figure 10, the ratio of perpendicular to parallel energy and GCH dissipations at saturation decreases as $I$ increases. When saturation is reached with a small value of $\chi$, perpendicular transfer and dissipation still strongly dominate the parallel ones. This is no longer the case when $I$ is large, as illustrated with run $R_5$. Note that in this simulation, the initial value of $\chi$ was increased in order to shorten the time needed to reach saturation. For such a large value of $\chi$, the parallel dissipation rapidly becomes comparable to the perpendicular one. It is to be noted that this situation is not specific to imbalanced turbulence but also occurs for balanced simulations at large $\chi$ (not shown), indicating that the breakdown of the asymptotics is not limited to imbalanced runs. Moreover, it is conspicuous that the ratio of the perpendicular to parallel dissipation is always smaller for the GCH than for the energy, pointing out that the barrier acts predominantly on the GCH.

Figure 10. Time evolution (averaged over three Alfvén crossing times) of the ratio between perpendicular and parallel dissipations of energy (a) and GCH (b), for the simulations shown in figure 8 (same colour code).

As the nonlinear parameter $\chi$ is increased, the transition zone is less extended and displays a shallower slope (not shown), except possibly when the imbalance is also increased, as in the simulations presented in Meyrand et al. (Reference Meyrand, Squire, Schekochihin and Dorland2021).

3.3 Transition zone with negligible parallel dissipation

The transition zone not only exists in regimes where $\nu _z=\nu _{\perp }$ (as in all the runs presented in the previous sections), but persists when the parallel viscosity $\nu _z$ is decreased while keeping $\nu _{\perp }$ unchanged. Figure 11(a) shows the transverse energy spectra $E^\pm (k_{\perp })$ when $\nu _z$ is reduced by a factor $40$ relative to its value in run $R_2$. Parallel viscosity can in fact be taken arbitrarily small and even zero without significantly affecting the accuracy of the numerical simulation during a significant integration time. As seen in figure 11(b), when reducing $\nu _z$ the perpendicular dissipation remains unchanged, while the parallel one (which is always sub-dominant) decreases further and the contribution of the spectral range close to $k_{\perp } = 1$ still dominates that of the small transverse scales. Regimes with small parallel dissipation are in fact relevant for the solar wind where the dynamics being quasi-transverse (Sahraoui et al. Reference Sahraoui, Goldstein, Belmont, Canu and Rezeau2010), the parallel transfer and consequently the parallel dissipation at the smallest scales are much weaker than the transverse ones. The above results suggest that the existence of the transition zone is not a consequence of the parallel dissipation, at least for small values of the nonlinear parameter. A possible origin of the transition in this case could be the interactions between co-propagating waves that take place at the sub-ion scales, as suggested by Voitenko & De Keyser (Reference Voitenko and De Keyser2016) and Gogoberidze & Voitenko (Reference Gogoberidze and Voitenko2020). This issue is addressed in the next section, using a semi-phenomenological model.

Figure 11. (a) Energy spectra $E^\pm$, averaged over the time interval [15 000, 20 000], for run $R_6$ with a parallel viscosity smaller by a factor of 40 than in run $R_2$. (b) Perpendicular (solid lines) and parallel (dotted lines) energy dissipations (averaged over three Alfvén crossing times) for simulations performed with different parallel viscosities, all other parameters being kept unchanged (red, run $R_2$; blue, run $R_6$; cyan, run $R_7$).

4.1 A phenomenological model including co-propagating wave interactions

We here concentrate on a diffusion model in spectral space which was first derived in Passot & Sulem (Reference Passot and Sulem2019), and then reproduced in Miloshevich et al. (Reference Miloshevich, Laveder, Passot and Sulem2020) using a heuristic approach based on the conservation of the energy and and of the GCH, together with a phenomenological estimate of the transfer times. This model does not capture the interactions between co-propagating waves, not only in the MHD range, where they do not exist, but also in the sub-ion range where co-propagative kinetic AWs can actually interact. In Passot & Sulem (Reference Passot and Sulem2019), following Voitenko & De Keyser (Reference Voitenko and De Keyser2016), we enriched the model by retaining interactions between triad wavenumbers which are comparable but not necessarily asymptotically close. We introduce the nonlinear characteristic frequencies for the mode ${\boldsymbol k}$ propagating in the $\pm$ direction, $\gamma _{\boldsymbol k}^{\pm (\uparrow \uparrow )}= V_{\boldsymbol k}^{(\uparrow \uparrow )} \sqrt {k_{\perp }^3 {E}^\pm (k_{\perp })}$ and $\gamma _{\boldsymbol k}^{\pm (\uparrow \downarrow )}= V_{\boldsymbol k}^{(\uparrow \downarrow )} \sqrt {k_{\perp }^3 {E}^\mp (k_{\perp })}$ (where $V_{\boldsymbol k}^{(\uparrow \uparrow )}$ and $V_{\boldsymbol k}^{(\uparrow \downarrow )}$ estimate the strength of the corresponding interactions), associated with the interactions between co-propagating $(\uparrow \uparrow )$ or counter-propagating $(\uparrow \downarrow )$ waves, respectively. Restricting the discussion to the case where both waves undergo strong nonlinear interactions, we can define global inverse transfer times as $(\tau _{tr, G}^\pm )^{-1} = \gamma _{\boldsymbol k}^{\pm (\uparrow \uparrow )} +\gamma _{\boldsymbol k}^{\pm (\uparrow \downarrow )}$. The ratio $\alpha (k_{\perp })=V_{\boldsymbol k}^{(\uparrow \uparrow )} / V_{\boldsymbol k}^{(\uparrow \downarrow )}$ should be an increasing function of $k_{\perp }$, which is essentially zero in the MHD range and saturates to a finite value in the dispersive range. Identifying the inverse transfer times $\gamma _{\boldsymbol k}^{\pm (\uparrow \downarrow )}$ associated with the counter-propagating waves with the nonlinear frequencies $\sqrt {k_{\perp }^3 v_{{\rm ph}}^2 {E}^\mp (k_{\perp })}$, the new inverse transfer times are then rewritten as

(4.1)\begin{equation} (\tau_{{\rm tr}, G}^\pm)^{{-}1} \approx k^{3/2}_{{\perp}} v_{{\rm ph}} \left(({E}^\mp)^{1/2}+ \alpha(k_{{\perp}})({E}^\pm)^{1/2} \right). \end{equation}

Defining $u^\pm = {E}^\pm / k_{\perp }$, the equation governing the evolution of ${E}^\pm (k_{\perp })$ is rewritten, in the stationary regime, as

(4.2)\begin{equation} \frac{{\rm d}}{{\rm d} k_{{\perp}}} u^\pm(k_{{\perp}})={-}C\frac{\varepsilon\pm\eta v_{{\rm ph}}}{k^5_{{\perp}} v_{{\rm ph}}}\frac{1} {u^\mp(k_{{\perp}})^{1/2} + \alpha(k_{{\perp}}) u^\pm(k_{{\perp}})^{1/2} }, \end{equation}

where the energy and GCH fluxes $\epsilon$ and $\eta$ are a priori functions of $k_{\perp }$ and $C$ is a numerical constant. This equation is similar to (7.34) of Passot & Sulem (Reference Passot and Sulem2019) when $\nu =0$, except that the weighting factor of the co-propagating wave interactions is here taken to be $\alpha (k)$ instead of its square.

The function $\alpha$, which is taken as $\alpha (k_{\perp })=3((1+ck_{\perp }^2)^{1/2}-1)/((1+ck_{\perp }^2)^{1/2}+1)$ in Voitenko & De Keyser (Reference Voitenko and De Keyser2016), with $c$ denoting a slowly varying function of $k_{\perp }$, can also be written as

(4.3)\begin{equation} \alpha(k_{{\perp}})= d\left( \frac{v_{{\rm ph}}(k_{{\perp}})-s}{v_{{\rm ph}}(k_{{\perp}})+s}\right), \end{equation}

where $s = \sqrt {2/\beta _e}$ is the Alfvén velocity (in sound speed units) and $d$ is a free constant parameter. An important remark, related to the existence of the helicity barrier, is that (4.2) does not have a satisfactory solution when $\eta \ne 0$. The linear growth of the phase velocity $v_{{\rm ph}}$ at small scales implies that for $k_{\perp }$ large enough, depending on the sign of $\eta$, one of the spectra ${E}^+$ or ${E}^-$ increases with $k_{\perp }$.

4.2 Numerical integration of the phenomenological model

As discussed in Passot & Sulem (Reference Passot and Sulem2019), (4.2) can be considered either as an initial-value problem when the energy and GCH fluxes $\varepsilon$ and $\eta$ are given and the spectra ${E}^\pm$ prescribed at a wavenumber $k_{\perp }=k_0$ or, alternatively, as a nonlinear eigenvalue problem for $\varepsilon$ and $\eta$, when the spectra are specified at $k_0$ and prescribed to decrease to zero at infinity. Numerical integration of these equations is best performed as an initial-value problem where $k_0$ is taken close to the largest wavenumber of the spectral domain. Indeed, the eigenvalue problem (where ideally $k_0$ is chosen at the lower end of the spectral domain) requires a shooting method and turns out to be extremely difficult to solve numerically, due to its instability. Indeed, under the effect of any small perturbation, the small-scale spectrum takes the form of an absolute equilibrium or goes to zero at a finite value of $k_{\perp }$, depending on the boundary conditions.

In the integrations presented below, we choose $\beta =2$, $\tau =1$, $k_0=200$ and small enough initial conditions such as to obtain a solution where the spectra tend to zero for a value of $k_{\perp }$ slightly in excess of $k_0$. As suggested by the numerical simulations of the gyrofluid equations, in the conditions where $\chi$ is not too large, $\varepsilon$ can be taken constant. We used $\varepsilon =1$ (close to the almost constant energy flux observed in the gyrofluid simulations) and $C=15$. Differently, the gyrofluid simulations show that the flux of GCH undergoes a significant decrease at the barrier.Footnote ² To comply with this observation, we were led to prescribe that $\eta$ decreases faster than $1/v_{{\rm ph}}(k_{\perp })$. For simplicity, we here chose to take $\eta (k_{\perp })=\eta _0/v_{{\rm ph}}^{\gamma }(k_{\perp })$. The numerical factor $\eta _0=0.82$ ensures that the imbalance at large scales tends to a value close to $100$, as in most of the numerical simulations presented in § 3, while the exponent $\gamma =1.1$ leads to a small-scale imbalance which is also compatible with the simulations. A larger value of $\gamma$ would produce a smaller imbalance at small scales. Note that decaying transfer rates also arise in nonlinear diffusion models of turbulence retaining Landau damping (Howes, Tenbarge & Dorland Reference Howes, Tenbarge and Dorland2011).

In figure 12, we display the ${E}^\pm$ spectra (figure 12a,c) and their local spectral exponents (figure 12b,d) for the cases $d=3$ (figure 12a,b) and $d=0$ (figure 12c,d), the other parameters being taken as mentioned above. In the absence of co-propagating wave interactions ($d=0$), the ${E}^+$ spectrum exponent undergoes a transition, at $k_{\perp }$ close to unity, between the Kolmogorov exponent $-5/3$ and a value slightly steeper than the theoretical value $-7/3$, as expected in the presence of imbalance (Passot & Sulem Reference Passot and Sulem2019). When $d=3$ (a value compatible with kinetic theory, used by Voitenko & De Keyser (Reference Voitenko and De Keyser2016)), a clear overshoot of the spectral index is observed, associated with the existence of a transition zone which starts for $k_{\perp }$ slightly smaller than unity, the steeper slope being equal to -3.28 at $k_{\perp }=1.18$. Near this wavenumber, the ${E}^-$ spectrum displays a weak flattening.

Figure 12. Spectra $E^\pm (k_{\perp })$ (a,c) and local slopes (b,d) predicted by the phenomenological model for $d=3$ (a,b) and $d=0$ (c,d), with $C=15$, $\epsilon =1$, $\eta (k_{\perp })=0.82 /v_{{\rm ph}}^{1.1}(k_{\perp })$.

Several remarks can be made concerning the influence of the various parameters of the model, which, at this stage, remains only qualitative. First, the slope of the transition zone is steeper (shallower) as $d$ in increased (decreased). This further indicates the role of the co-propagating wave interactions in the generation of this transition zone.

At larger imbalance (e.g. when taking $\eta _0=0.92$, which leads to a large-scale imbalance close to $I=1000$), the transition zone is slightly steeper, with a slope reaching $3.52$ at $k_{\perp }=1.11$.

In the absence of co-propagating wave interactions, and when the variation of $\eta$ is such that $\eta v_{{\rm ph}}(k_{\perp })$ is constant, the imbalance remains independent of $k_{\perp }$. This is no longer the case when $\eta$ decreases faster than $1/v_{{\rm ph}}(k_{\perp })$. As $\gamma$ is increased from the value $1.1$ chosen in the figure, the imbalance at small scale is strongly reduced but the transition zone takes a different form. While the slope of the ${E}^+$ spectrum is not as steep as for smaller values of $\gamma$, the ${E}^-$ spectrum becomes shallower. The case where $\eta$ decreases faster can be associated with a stronger transfer and/or dissipation in the $z$ direction and thus a larger value of $\chi$. This effect is seen in some simulations but a systematic study has not been performed.