2.1. Equations and conservation laws

A description of the Alfvén wave dynamics from the MHD to the electron scales, in an electron–proton plasma, is provided by the two-field gyrofluid model which involves the gyrokinetic scaling corresponding to a strong spectral anisotropy and weak nonlinearity. In the absence of dissipation and driving, it involves two equations for the electron gyrocentre number density $N_e$ and the parallelFootnote ¹ component of the magnetic potential $A_{\|}$, in the formFootnote ²

(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 L_e A_{\|} -\left [\varphi,\frac{2\delta^2}{\beta_e}\varDelta_{\perp} A_{\|} \right ] + \left [B_z,\frac{2\delta^2}{\beta_e}\varDelta_{\perp} A_{\|}\right ] + \boldsymbol{\nabla}_{\|} \left (\varphi-N_e-B_z \right )=0. \end{gather}

In the above equations, the parallel magnetic fluctuations $B_z$ and the electron gyrocentre number density $N_e$ are related to the electric potential $\varphi$ by $B_z= M_1 \varphi$ and $N_e= -M_2 \varphi$, where $M_1$ and $M_2$ are Fourier multiplier operators written as $M_1 = L_1^{-1}L_2$ and $M_2= L_3 + L_4 L_1^{-1}L_2$ (here, $\delta ^2 = m_e/m_i$ denotes the electron to proton mass ratio), the operators $L_i$ being defined as

(2.3)\begin{gather} L_1 = \frac{2}{\beta_e} +(1+2\tau)(\varGamma_0-\varGamma_1), \end{gather}

(2.4)\begin{gather}L_2 = 1 +\frac{1-\varGamma_0}{\tau} - \varGamma_0 +\varGamma_1, \end{gather}

(2.5)\begin{gather}L_3 = \frac{1-\varGamma_0}{\tau} -\delta^2\varDelta_{\perp}, \end{gather}

(2.6)\begin{gather}L_4 = 1-\varGamma_0+\varGamma_1, \end{gather}

(2.7)\begin{gather}L_e = 1-\frac{2 \delta^2}{\beta_e} \varDelta_{\perp}. \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 })$ associated with the Fourier multiplier $\varGamma _n(\tau k_{\perp }^2)$, defined by $\varGamma _n(x) = I_n(x) e^{-x}$, where $I_n$ is the modified Bessel function of first type of order n. Furthermore, $\tau$ denotes the ratio of the proton to the electron temperatures at equilibrium. For a scalar function $f$, the parallel gradient operator $\boldsymbol {\nabla }_{\parallel }$ is defined by

(2.8)\begin{equation} \boldsymbol{\nabla}_{\parallel}\, f=-[A_{\parallel} , f]+\frac{\partial f}{\partial z}. \end{equation}

One recovers the fluctuating magnetic field from the expression of $B_z$ given above and ${\boldsymbol {B}}_{\perp } =\boldsymbol {\nabla }_{\perp }\times (A_{\|}\hat {{\boldsymbol {z}}})$.

Characterizing the plasma equilibrium state by the number density $n_0$, the temperatures $T_{i0}$ and $T_{e0}$ and the ambient field $B_0$, the model is written above in a non-dimensional form, using the following units: time is normalized to the inverse ion gyrofrequency $\varOmega _i^{-1} = (eB_0/(m_ic))^{-1}$ (where $e$ is the proton charge and $c$ the speed of light), and lengths by the sonic Larmor radius $\rho _s= c_s/\varOmega _i$, where $c_s=\sqrt {T_{e0}/m_i}$ is the sound speed. Thus wavenumbers and wavevector components are measured in units of $\rho _s^{-1}$. The other normalization factors include $B_0$ for the parallel magnetic fluctuations $B_z$, $B_0\rho _s$ for the parallel magnetic potential $A_{\|}$, the equilibrium number density $n_0$ for the electron gyrocentre density $N_e$ and $T_e/e$ for the electric potential $\varphi$. Furthermore, the parameter $\beta _e= 8 {\rm \pi} n_0 T_{e0} /B_0^2$ is the ratio of the equilibrium electron pressure to the magnetic pressure due to the ambient field.

At the level of the linear approximation, the KAW parallel phase velocity $v_{\textrm {ph}}= \omega /k_z$ is given by

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

where the hat refers to the Fourier symbol of the operator and $s= (2/\beta _e)^{1/2}$ is the Alfvén velocity $v_A$ measured in sound speed units ($s=v_A/c_s$). In figure 1 we plot $v_{\textrm {ph}}/s$ versus the transverse wavenumber $k_{\perp }$, when electron inertia is neglected. For all values of $\beta _e$, the ratio $v_{\textrm {ph}}/s$, which can also be viewed as the wave frequency normalized by the Alfvén wave frequency $v_A k_{\|}$, is asymptotically equal to unity at large scale (dispersionless limit). On the other hand, while for small ($\beta _e = 0.2$) or moderate ($\beta _e=2$) values of the electron beta parameter, $v_{\textrm {ph}}/s$ involves a monotonic transition to a $k_{\perp }$-scaling at sub-ion scales, for larger $\beta _e$ (e.g. 10 or 50), it displays a local minimum at a wavenumber that (in inverse sonic Larmor radius units) increases with $\beta _e$ at fixed $\tau$ and only depends on $\beta _i$ when measured in units of the ion Larmor radius $\rho _i$. A similar behaviour is observed in figures 3 ($\tau = 100$) and 5 ($\tau =1$) of Howes et al. (Reference Howes, Cowley, Dorland, Hammett, Quataert and Schekochihin2006) where this quantity is computed from the gyrokinetic theory and also in the framework of the full linear kinetic theory taken in the quasi-transverse limit. As discussed below, it turns out that the sensitivity of the dispersion to the transverse wavenumber, as measured by the parallel phase velocity $v_{\textrm {ph}}$, has an important effect on the nonlinear dynamics of the cascade. In particular, it has an impact on the generation of a finite-size condensate and on the formation of large-scale coherent structures.

Figure 1.

Normalized parallel phase velocity $v_{\textrm {ph}} \sqrt {{\beta _e}/{2}}$ (or, equivalently, ratio of the KAW frequency $v_{\textrm {ph}} k_{\|}$ to the Alfvén frequency $v_A k_{\|}$) versus the transverse wavenumber $k_{\perp }$ for various values of $\beta _e$.

In physical space, the eigenmodes, which can be referred to as Elsasser potentials, are given by

(2.10)\begin{equation} \mu^{\pm} = \varLambda \varphi \pm s A_{\|}, \end{equation}

where

(2.11)\begin{equation} \varLambda = D_e^{-1} (1+M_2-M_1)^{1/2} M_2^{1/2} , \end{equation}

with $D_e^2 = (-\varDelta _{\perp }) L_e$. They obey

(2.12)\begin{align} \partial_t \mu^{\pm} &\pm V_{\textrm{ph}} \partial_z \mu^{\pm} +\frac{1}{4} \varLambda^{-1} D_e^{-2} M_3\{[\varLambda^{-1} M_3 (\mu^++\mu^-) , \varLambda^{-1}M_2 (\mu^+ + \mu^-)]\nonumber\\ & +[(\mu^+-\mu^-), \varDelta_{\perp} (\mu^+-\mu^-)]\} \pm \frac{1}{4} D_e^{-2}\varDelta_{\perp} \{ [(\mu^+ -\mu^-), \varLambda^{-1} M_2 (\mu^++\mu^-)]\nonumber\\ &+ [L_e (\mu^+ -\mu^-), \varLambda^{-1} (1-M_1)(\mu^+ +\mu^-)]\}=0. \end{align}

The systems (2.1) and (2.2) or (2.12) preserve the energy ${\mathcal {E}}$ and the GCH ${\mathcal {C}}$, defined as

(2.13)\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) \,\textrm{d}^3 {x} \nonumber\\ \hspace{-9pc}=\frac{1}{4} \int \{ (D_e \mu^+)^2 + (D_e \mu^-)^2\} \,\textrm{d}^3x, \end{gather}

(2.14)\begin{gather}{\mathcal{C}} =-\int N_e \left( 1 - \frac{2\delta^2}{\beta_e}\varDelta_{\perp} \right) A_{\|}\,\textrm{d}^3 {x} = \frac{1}{4} \int \left \{\left ( V_{\textrm{ph}}^{-1/2}D_e \mu^+\right)^2 - \left ( V_{\textrm{ph}}^{-1/2}D_e \mu^-\right)^2 \right\} \,\textrm{d}^3x, \end{gather}

where $V_{\textrm {ph}}$ is the operator in physical space, which in Fourier space corresponds to multiplication by $v_{\textrm {ph}}$. Note that ${\mathcal {C}}$ was defined with an opposite sign in Passot et al. (Reference Passot, Sulem and Tassi2018).

2.2. Normal-field formulation

As is common with non-canonical Hamiltonian systems, the associated Poisson bracket possesses Casimir invariants, corresponding to $C_{\pm } = \int G^{\pm } \,\textrm {d}^3x$, where $G^{\pm } = L_e A_{\|} \pm \delta N_e$ are referred to as normal fields. In terms of these fields, the two-field gyrofluid model can be rewritten in the form

(2.15)\begin{equation} \partial_t G^{\pm} + [\varphi^{\pm}, G^{\pm}] + \partial_z \left ( \varphi^{\pm} \mp \frac{1}{\delta} G^{\pm} \right) =0, \end{equation}

where $\varphi ^{\pm } =\varphi -B_z \pm ({1}/{\delta }) A_{\|}$. In terms of these variables, GCH simply reads

(2.16)\begin{equation} {\mathcal{C}} = -\frac{1}{4\delta}\int \left \{(G^+)^2 - (G^-)^2 \right \} \, \textrm{d}^3x. \end{equation}

Similarly,

(2.17)\begin{equation} {\mathcal{E}} = \frac{1}{2} \int \left \{\frac{1}{2\beta_e} L_e^{-1} |\boldsymbol{\nabla}_{\perp} (G^+ + G^-)|^2 + \frac{1}{4\delta^2}M_3 M_2^{-1} (G^+ -G^-)^2 \right\}\, \textrm{d}^3x. \end{equation}

Note that the system of (2.15) for $G^{\pm }$ becomes degenerate in the limit $\delta \to 0$, as the two equations become identical and reproduce the equation for $A_{\|}$. The equation for $N_e$ corresponds to the first order in a development in $\delta ^2$. As seen below, this formulation nevertheless is of major interest for estimating the energy and GCH spectral transfers estimated in appendix B.

2.3. Fjørtoft argument for an inverse cascade

The Fjørtoft (Reference Fjørtoft1953) argument originally referred to a simple method to predict the existence of an inverse cascade by noting the impossibility of simultaneous direct cascades of two ideal invariants whose spectra differ by a power of the wavenumber. Other versions of the argument can be found in Nazarenko (Reference Nazarenko2011) (and references therein). We first revisit the use of the Fjørtoft (Reference Fjørtoft1953) argument to support the existence of an inverse GCH cascade at sub-ion scales (Schekochihin et al. Reference Schekochihin, Cowley, Dorland, Hammett, Howes, Quataert and Tatsuno2009). In this range, the two-field model simplifies and reduces, when neglecting electron inertia, to electron reduced MHD (ERMHD) (Schekochihin et al. Reference Schekochihin, Cowley, Dorland, Hammett, Howes, Quataert and Tatsuno2009; Boldyrev et al. Reference Boldyrev, Horaites, Xia and Perez2013), under the conditions $\tau k_{\perp }^2 \gg 1$, $\tau \sim 1$:

(2.18)\begin{gather} \partial_t A_{\|} + \left(1+\frac{1}{\tau}\right) \boldsymbol{\nabla}_{\|} \varphi=0, \end{gather}

(2.19)\begin{gather}\partial_t \varphi - \frac{\dfrac{2\tau}{\beta_e}}{1+\dfrac{\beta_e}{2}(1+\tau)} \boldsymbol{\nabla}_{\|}\varDelta_{\perp} A_{\|} =0. \end{gather}

In this regime, the invariants read

(2.20)\begin{gather} {\mathcal{E}} = \frac{1}{\beta_e}\int \left (|{\boldsymbol{B}}|^2 + \frac{2}{\beta_e} \frac{1}{1+\tau} B_z^2 \right )\,\textrm{d}^3x, \end{gather}

(2.21)\begin{gather}{\mathcal{C}} = \left ( 1 + \frac{2}{\beta_e}\frac{1}{1+\tau} \right ) \int A_{\|}B_z\,\textrm{d}^3x. \end{gather}

In the incompressible limit where the beta parameter tends to infinity, the GCH invariant reduces, in the quasi-transverse limit, to the generalized magnetic helicities of EMHD (Biskamp et al. Reference Biskamp, Schwarz, Zeiler, Celani and Drake1999), or of extended MHD when the ion velocity and electron inertia are taken to zero ((35) of Abdelhamid, Lingam & Mahajan (Reference Abdelhamid, Lingam and Mahajan2016) and (29) of Miloshevich, Lingam & Morrison (Reference Miloshevich, Lingam and Morrison2017)).

The magnetic field ${\boldsymbol {B}}$ can be written ${\boldsymbol {B}} = \boldsymbol {\nabla }_{\perp } \times (A_{\|} {\hat {\boldsymbol {z}}}) + B_z {\hat {\boldsymbol {z}}}$. Introducing ${\boldsymbol {A}}_{\perp }$ such that $B_z = {\hat {\boldsymbol {z}}}\boldsymbol {\cdot } (\boldsymbol {\nabla }_{\perp } \times {\boldsymbol {A}}_{\perp })$, one has $\int A_{\|} B_z\,\textrm {d}^3x = (1/2)\int {\boldsymbol {A}}\boldsymbol {\cdot } {\boldsymbol {B}}\,\textrm {d}^3x$ with ${\boldsymbol {A}} = {\boldsymbol {A}}_{\perp } + A_{\|} {\hat {\boldsymbol {z}}}$ (see e.g. (F6) of Schekochihin et al. (Reference Schekochihin, Cowley, Dorland, Hammett, Howes, Quataert and Tatsuno2009)), together with $\int |{\boldsymbol {B}}|^2\,\textrm {d}^3x = \int |\boldsymbol {\nabla }_{\perp } \times {\boldsymbol {A}}|^2\,\textrm {d}^3x$.

In the large $\beta _e$ limit, ${\mathcal {E}}= ({1}/{\beta _e})\int |{\boldsymbol {B}}|^2\,{\textrm {d}x}$ and, like in incompressible MHD without ambient field (see e.g. Pouquet et al. (Reference Pouquet, Rosenberg, Stawarz and Marino2019) for a recent review), assuming that the helicity is of a given sign and that it is maximal, one gets $E_c(k_{\perp })= (1/k_{\perp }) E(k_{\perp })$, which leads to conjecture the existence of an inverse cascade of magnetic helicity, by generalizing the argument developed by Fjørtoft (Reference Fjørtoft1953) for two-dimensional incompressible turbulence.

In the case of finite $\beta _e$, the energy and GCH spectra cannot be directly related to each other. Nevertheless, when turbulence is not too strong, the first and third terms in the energy given by the first line of (2.13) are approximatively in equipartition (this property was checked to be accurately satisfied in the present simulations). In the ERMHD regime, where the coefficients $M_1$ and $M_2$ are scale-independent, it is easily seen that ${\mathcal {E}} \sim \int |{\boldsymbol {B}}|^2\,\textrm {d}^3x$, which then permits use of the above argument. In summary, the phenomenological argument for the existence of an inverse GCH cascade at finite $\beta _e$ requires both equipartition between total magnetic and internal energies and single-sign maximal helicity, while in the large $\beta _e$ limit the first condition is not necessary.

In contrast, when considering the Hall reduced MHD equations for dispersive Alfvén waves ((E19) and (E20) of Schekochihin et al. (Reference Schekochihin, Cowley, Dorland, Hammett, Howes, Quataert and Tatsuno2009) or (2.24) and (2.25) of Passot & Sulem (Reference Passot and Sulem2019)), which correspond to the regime $\tau \ll 1$ and $\tau k_{\perp }^2 \ll 1$, $N_e ({\boldsymbol {k}})\approx - k_{\perp }^2 \varphi ({\boldsymbol {k}})$ and $B_z({\boldsymbol {k}}) \approx (1+ 2/\beta _e)^{-1} k_{\perp }^2\varphi {(\boldsymbol {k}})$, one has

(2.22)\begin{gather} {\mathcal{E}} = \frac{1}{2} \int \left [ \frac{2}{\beta_e}|\boldsymbol{\nabla}_{\perp} A_{\|} |^2 + \left (1 + \frac{k_{\perp}^2}{1+ \dfrac{\beta_e}{2}}\right) \left(1+ \frac{2}{\beta_e}\right)^2 \frac{1}{k_{\perp}^2} B_z^2 \right ]\,\textrm{d}^3x, \end{gather}

(2.23)\begin{gather}{\mathcal{C}} = \left( 1 + \frac{2}{\beta_e}\right) \int A_{\|} B_z\,\textrm{d}^3 x =\frac{1}{2}\left( 1 + \frac{2}{\beta_e}\right) \int {\boldsymbol{A}}\boldsymbol{\cdot} {\boldsymbol{B}}\,\textrm{d}^3x. \end{gather}

Two regimes are then to be distinguished. When formally taking the limit $k_{\perp } \to \infty$ (which in the Hall reduced MHD model requires en extremely small value of $\tau$) and assuming equipartition of the magnetic and internal energies, one writes the total energy in the form

(2.24)\begin{equation} {\mathcal{E}}= \frac{1}{\beta_e}\int|{\boldsymbol{B}}_{\perp}|^2\,\textrm{d}^3x = \left ( \frac{1 + \dfrac{2}{\beta_e}}{2+ \dfrac{2}{\beta_e}}\right)\frac{2}{\beta_e} \int |{\boldsymbol{B}}|^2\,\textrm{d}^3x. \end{equation}

The same argument as in the case of ERMHD can then be used to conjecture the existence of an inverse GCH cascade. As expected, this argument does not apply in the reduced MHD regime corresponding to the limit $k_{\perp } \to 0$. Using in this case the strongest assumption of local-in-scale energy balance, one gets

(2.25)\begin{equation} B_z({\boldsymbol{k}}) = \sqrt{\frac{2}{\beta_e}} \frac{1}{1+ \dfrac{2}{\beta_e}}k_{\perp} |{\boldsymbol{B}}_{\perp}({\boldsymbol{k}})|, \end{equation}

and thus

(2.26)\begin{equation} {\mathcal{C}} = \sqrt{\frac{2}{\beta_e}} \int |{\boldsymbol{B}}_{\perp}|^2\,\textrm{d}^3 x. \end{equation}

In this case, ${\mathcal {E}}$ and ${\mathcal {C}}$ are proportional and no energetic condition prevents the existence of a simultaneous direct cascade of the two invariants. From the above argument, we can expect that an inverse cascade of GCH develops when turbulence is driven at sub-ion scales and that this cascade progressively slows down when approaching the non-dispersive scales.

4.1. Spectra and fluxes

When choosing $\beta _e=2$ and $k_f = 13$, an inverse cascade of GCH quickly develops, giving rise to a spectral bump propagating towards large transverse scales. This is illustrated in figure 2(a). In this case, the inverse cascade extends over more than one decade before reaching the weakly dispersive scales. The observed dynamics is qualitatively similar to that obtained in numerical simulations of EMHD (Kim & Cho Reference Kim and Cho2015) and also of incompressible MHD in the absence of an ambient field (Müller et al. Reference Müller, Malapaka and Busse2012; Linkmann & Dallas Reference Linkmann and Dallas2017). At the final time of the simulation, the energy has almost reached the lowest modes and a separate simulation is needed to address the situation where the cascade crosses $k_{\perp }=1$ and penetrates into the weakly dispersive range. It is discussed in § 6 with run $R_{1.3}^2$.

Figure 2.

Perpendicular (a,c) and parallel (b,d) energy and GCH spectra (a,b) and GCH fluxes (c,d) for run $R_{13}^2$. The integration time that extends up to $t=300$ (in $\varOmega _i^{-1}$ units) is indicated in the the colour bars, where $+$GCH refers to positive values and $-$GCH to the modulus of negative values (dashed lines visible at the smallest scales) of the GCH spectra. For clarity, the axes of the energy and GCH spectra have been shifted. The same convention is adopted for subsequent graphs.

The perpendicular GCH flux displayed in figure 2(c) shows at late times a remarkably constant range, associated with an inverse cascade, in the spectral region between the peak and the forcing. Interestingly, a self-similar inverse cascade develops in the parallel direction as well, as shown in figure 2(b), also associated with a negative and constant GCH flux (figure 2d). The existence of this parallel inverse cascade has important consequences on the types of structures that form in physical space, as discussed in § 8. When considering energy fluxes (not shown), we note that the perpendicular one is significant at wavenumbers larger than $k_f$ but very small for $k_{\perp } < k_f$. Energy is transferred to large scales together with the cross-helicity, but only the latter undergoes a true inverse cascade. In the parallel direction, the energy flux, although not displaying a constant range, is decaying more slowly as $k_{\perp }$ decreases away from $k_f$.

In figure 3, we show the Elsasser spectra $E^{\pm }(k_{\perp })$ (figure 3a) and the time evolution of the total Elsasser energies ${\mathcal {E}^{\pm } =\sum _{k_{\perp }} E^{\pm }}(k_{\perp })$ and of the positive/negative GCHs ${\mathcal {E}}_C^{\pm } =\sum _{k_{\perp }} E_C^{\pm }(k_{\perp })$ (figure 3b). The inverse GCH cascade is only associated with an inverse transfer of $E^+$. This behaviour is specific to cases where the ratio $\epsilon ^+/\epsilon ^-$ is small or moderate, as also seen in EMHD simulations by Kim & Cho (Reference Kim and Cho2015). In this run, a non-negligible transfer of $E^-$ is nevertheless visible, which might develop into a different type of inverse cascade involving longer spatial and temporal scales. Simulations with a spectral range extending to much larger scales and spanning a longer time interval would be necessary to characterize this effect more precisely. Both $\mathcal {E}^+$ and $\mathcal {E}_C^+$ grow at long time as power laws with exponents similar but not identical to those reported in Kim & Cho (Reference Kim and Cho2015) in the case of EMHD. Furthermore, we verified that, while keeping the parameters of run $R_{13}^{2}$, but injecting only one type of wave (infinitely imbalanced driving), the other type of wave is generated by interactions of driven modes. While the energy of the dominant mode is transferred to the large scales through the propagation of a spectral bump whose envelope obeys a $k_{\perp }^{-3/2}$ scaling law, the subdominant wave undergoes a self-similar cascade with a $k_{\perp }^{-1}$ spectrum (not shown), indicating a behaviour of both spectra similar to that observed in EMHD simulations with maximal helicity injection (Kim & Cho Reference Kim and Cho2015).

Figure 3.

Perpendicular Elsasser spectra (a) and time evolution of the total Elsasser energies ${\mathcal {E}}^{\pm }$ and GCHs ${\mathcal {E}}_C^{\pm }$ (b) for run $R_{13}^2$.

4.2. Shell-to-shell transfers

In the following, we analyse the shell-to-shell transfers, i.e. the amount of energy and GCH transferred per unit time from one spectral shell in the transverse plane to another, after integration on the parallel wavenumbers. The transfer between parallel wavenumbers after integration on the transverse wavenumbers is also considered. These quantities provide useful information concerning the degree of (spectral) locality of the nonlinear interactions. Their expressions are derived in appendix B and explicitly given in (B 21) and (B 22). It is to be noted that in this context $R_{13}^2$ is generic and the statements about the non-local nature of transfers apply to the other simulations. However, as shown in the following subsection, other features may differ.

Figure 4 shows the GCH shell-to-shell transfers both in the transverse plane (figure 4a,c) and in the longitudinal direction (figure 4b,d), at an early time ($t=5$; figure 4a,b) and at a time close to the end of the simulation ($t=100$; figure 4c,d). Transfers in the perpendicular directions are computed from (B 21) and (B 22) of appendix B. Similar expressions are used for the longitudinal transfers, where cylindrical shells are replaced by slabs.

Figure 4.

Perpendicular (a,c) and parallel (b,d) transfers of GCH for run $R_{13}^2$ at (a,b) an early ($t=5$) or (c,d) a late ($t=100$) time of the simulation.

At early time, inverse transfers of GCH appear to be strongly non-local, i.e. off-diagonal, which indicates transfer between non-neighbouring shells. We observe the characteristic cross-signature which demonstrates the interactions of the forcing with other shells. The modes located on the upper vertical and left horizontal axes that intersect at the forcing wavenumber are associated with an inverse transfer close to the forcing. On the contrary, the non-local transfers between the forcing and the wavenumbers that are much larger are consistently associated with direct transfers.

The inverse GCH cascade is mainly driven by the non-local interactions between the forcing and the propagating bump, which can be seen as the orange (respectively cyan) areas located in the regions where $k_{\perp }$ (respectively $p_{\perp }$) is slightly smaller and larger than $k_f$ with $p_{\perp }$ (respectively $k_{\perp }$) between $5$ and $10$.

In the parallel direction, there is no significant inverse transfer but interacting triads involving the forcing and its first harmonic are clearly visible (see also the parallel spectrum at the first output in figure 2b). The oblique lines ($k_{\|} + p_{\|} = k_f$ and circular permutations) and those parallel to the axes ($k_{\|} = k_f$ or $p_{\|}=k_f$) correspond to the interactions between triads including the driving mode. No other interaction is significant at this early time.

The transfers at late times display more local features, in both the parallel and perpendicular directions. The perpendicular inverse transfer of GCH still involves non-local interactions between the spectral bump and the forcing, while energy (not shown) has predominantly local features. It is noteworthy to mention that early bump propagation is always dominated by the non-local cross-type interaction between the forcing and the bump. The parallel GCH transfer associated with the direct cascade towards scales smaller than the forcing is local, while the inverse parallel cascade involves both strongly local and non-local features.

5.1. Varying the amplitude of the fluctuations

When varying the energy injection rate, the nature of the turbulent cascade that develops at relatively early time changes significantly (see figure 5a). At large amplitude, a self-similar cascade starts developing (as clearly seen in run $R_{13s}^2$) but the spectrum rapidly distorts. When the driving is decreased (run $R_{13}^2$), the dynamics is that of a propagating spectral bump, as described in the previous section. Further decreasing the injection rate, the early dynamics is dominated by the development of a parametric decay instability of the pump waves driven by the random forcing (conspicuous in figure 5(a) where two bumps emerge at time $t\approx 110$). Indeed, in contrast with Alfvén waves in ideal MHD, three KAWs can interact resonantly, producing a decay instability (Voitenko Reference Voitenko1998a) (see appendix A for a brief description of this linear instability in the context of the two-field model, in the absence of injection and dissipation). This instability competes with the turbulence dynamics (see Shi et al. (Reference Shi, Li, Xiao and Wang2017) and Fu et al. (Reference Fu, Li, Guo, Li and Roytershteyn2018) for a study in the MHD context). Even though the growth rate becomes larger when the level of the magnetic fluctuations is increased, the transfer to neighbouring modes due to turbulence is more efficient, thus hiding or preventing the instability. When $\epsilon _E$ is small enough, as is the case for run $R_{13w}^2$, the growth rate is smaller, but the turbulence is also less efficient, making possible the appearance of the decay. A more detailed study of the resonance conditions is discussed in the next subsection.

Figure 5.

Snapshots at short (a) and long (b) times of the GCH spectra for runs $R_{13s}^2$ (1), $R_{13}^2$ (2) and $R_{13w}^2$ (3) corresponding to strong, intermediate and weak injection rates, respectively. Final times displayed in (a) are $t=1.6$ (1), $t=12.5$ (2), $t=110$ (3), and in (b) are $t=15$ (1), $t=300$ (2), $t=2850$ (3).

Decreasing $\epsilon _E$ only affects the early-time dynamics. Indeed, as seen in figure 5(b), the late evolution of the two runs with strong and intermediate driving displays similar features associated with a strongly bent spectrum with a bump at large scale, separated from the forcing by a spectral domain where the dynamics tends to settle to an absolute equilibrium.Footnote ⁵ The case of a weak driving has the same tendency but the evolution is much slower. Similar simulations that display decay at early time but with a larger value of $\chi _f$ show in a more conspicuous way an evolution towards a flat spectrum in the range between the driving and the maximum of the spectrum.

Even if the inverse cascades displayed by the simulations previously discussed may suggest the existence of qualitatively different kinds of cascades, a more likely scenario is that there is a continuous transition between a self-similar cascade and a spectral bump propagating to large scales. Inspection of the time evolution of individual energy and GCH transverse Fourier modes indicates that when the inverse cascade reaches a given wavenumber, the corresponding modes first grow and then decrease with a characteristic decay time that becomes shorter when wave effects are enhanced. The relative effect of the strength of the turbulence compared to that of the waves, measured by the nonlinearity parameter at the driving scale $\chi _f$, is then likely to characterize the cascade properties.

The nonlinear parameter $\chi ({\boldsymbol {k}})$ at the wavenumber ${\boldsymbol {k}}$ can be estimated as follows. From (2.2) with $\delta =0$, we estimate the nonlinear time $\tau _{\textrm {NL}}(k_{\perp })$ (defining ${\hat {M}}_3 = 1 -{\hat {M}}_1 + {\hat {M}}_2$) by

(5.1)\begin{equation} \tau_{\textrm{NL}}({\boldsymbol{k}})^{-1} = k_{\perp}^2 {\hat{M}}_3 |\varphi({\boldsymbol{k}})|. \end{equation}

The nonlinear parameter $\chi ({\boldsymbol {k}}) = \tau _{\textrm {NL}}({\boldsymbol {k}})^{-1} / (v_{\textrm {ph}} k_{\|})$ reads (using (2.9))

(5.2)\begin{equation} \chi({\boldsymbol{k}}) = \frac{1}{s} k_{\perp} {\hat{M}}_2^{1/2} {\hat{M}}_3^{1/2} \frac{1}{k_{\|}}|\varphi({\boldsymbol{k}})|. \end{equation}

At least when the nonlinearities are not too strong, one can expect equipartition of the magnetic and potential energy spectra (see e.g. Slepyan (Reference Slepyan2015) and references therein), a property which is observed to hold with good accuracy in most of the simulations considered in this paper (in particular, the ratio of the two spectra at $k_f$ is typically $1\pm 0.03$, the error reaching $4\,\%$ at large $\beta _e$). This leads to a phenomenological estimate that near the injection scale

(5.3)\begin{equation} |\varphi({\boldsymbol{k}})| \approx s {\hat{M}}_2^{-1/2} {\hat{M}}_3^{-1/2} |{\boldsymbol{B}}_{\perp}({\boldsymbol{k}})|, \end{equation}

which allows one to relate the nonlinear parameter $\chi _f$ at the injection wavevector ${\boldsymbol {k}}_f$ (whose perpendicular and parallel components are denoted $k_{\perp f}$ and $k_{\| \,f}$) to the magnitude of magnetic fluctuations $B_{\perp } (k_{\perp f})= \sqrt {k_{\perp } E_{B_{\perp }}(k_{\perp f})}$ (where $E_{B_{\perp }} (k_{\perp f})$ holds for the one-dimensional transverse magnetic spectrum) by

(5.4)\begin{equation} \chi_f = \frac{k_{\perp f}}{k_{\| \,f}} |B_{\perp}(k_{\perp f})|. \end{equation}

In order to have a qualitative relation between the level of magnetic field fluctuation $|B_{\perp }(k_{\perp f})|$ at the driving wavenumber and $\epsilon _E$, it is nevertheless possible to use a simple strong turbulence phenomenology and write, still assuming energy equipartition,

(5.5)\begin{equation} \epsilon_E\approx\frac{2}{\beta_e}|B_{\perp}({k_{\perp f}})|^2 \tau_{\textrm{NL}}(k_{\perp f})^{-1}\approx\frac{2}{\beta_e}|B_{\perp}({k_{\perp f}})|^3 k_{\perp f} v_{\textrm{ph}}(k_{\perp f}). \end{equation}

This prompts the definition of the nonlinearity parameter based on turbulent phenomenology by, up to a numerical constant,

(5.6)\begin{equation} \chi_f^t= \frac{k_{\perp f}}{k_{\| \,f}}\left ( \frac {\beta_e}{2} \frac{\epsilon_E}{k_{\perp f}v_{\textrm{ph}}(k_{\perp f})} \right )^{1/3}. \end{equation}

Assuming that $\tau$ is of order unity, we find, using table 1 in Passot & Sulem (Reference Passot and Sulem2019), that at large scales $v_{\textrm {ph}}\sim s$, and at small scales $v_{\textrm {ph}}\sim s k_{\perp }$ when $\beta _e\lesssim 1$ and $v_{\textrm {ph}}\sim s^2 k_{\perp }$ when $\beta _e$ is large. For a driving acting at sub-ion scales, $\chi _f^t\sim \epsilon _E^{1/3}k_{\perp f}^{1/3}k_{\| \,f}^{-1}\beta _e^{1/2}$ for small or moderate $\beta _e$ and $\chi _f^t\sim \epsilon _E^{1/3}k_{\perp f}^{1/3}k_{\| \,f}^{-1}\beta _e^{2/3}$ for large $\beta _e$.

The expression of the $\chi _f$ parameter suggests that the dynamics will be sensitive to the angle of the injected KAWs with respect to the ambient magnetic field. Indeed, in simulations for which $k_{\| \,f}\ll k_{\perp f}$, which makes $\chi _f$ very large, the dynamics is quasi-two-dimensional. The waves have very long time scales and the inverse cascade is more self-similar. Moreover, in these conditions, the dynamics is influenced by the quasi-conservation of the squared magnetic potential, which is observed to undergo a self-similar inverse cascade (not shown).

Note that for small $\beta _e$, where $v_{\textrm {ph}} \sqrt {{\beta _e}/{2}}$ has a finite limit, the parameter $\chi _f^s$ becomes independent of $\beta _e$ when velocities are measured in units of Alfvén speed instead of sound speed and the length unit is kept unchanged. Indeed, in this limit, after a proper rescaling of the dependent and independent variables, the equations become independent of $\beta _e$ and reproduce the model of Zocco & Schekochihin (Reference Zocco and Schekochihin2011), taken in the isothermal limit. Interestingly, the small-$\beta _e$ limit is already accurately approached when $\beta _e=0.2$. In contrast, at large $\beta _e$, the parameter $\chi _f^t$ increases with $\beta _e$, even when using the Alfvén velocity unit.

The values of $\chi _f^t$ for the different runs considered in this paper are given in table 1, together with $|B_{\perp }(k_{\perp f})|$, which is measured immediately after the transient phase, at the beginning of the simulation. Note that, for some runs, $|B_{\perp }(k_{\perp f})|$ varies significantly during the time evolution. The ratio $\chi _f/\chi _f^t$ has a mean value of $1.78$ with a standard deviation of $0.43$. Both quantities vary in the same way from one run to another at moderate values of $\beta _e$, but not necessarily at larger values. Although relatively large, such uncertainty nevertheless permits using $\chi _f^t$ as a crude prediction of the nonlinear parameter $\chi _f$ which can be only measured a posteriori. It is natural to attempt to isolate the parameter $\chi _f$ to see if it solely controls the onset of self-similar bump transitions, etc. However, since it is not possible to precisely assign $\chi _f$ when initializing the simulations, we would have to initiate multiple runs to be able to select the appropriate desired values. Such ensemble would drain computational resources. Nevertheless, this conjecture that $\chi _f$ is the governing parameter appears to be consistent with the various simulations presented in this section.

5.2. Analysis of a decay event

Another example of the decay instability is encountered when, taking the same injection rate and $\beta _e$ as for run $R_{13}^2$, the driving wavenumber $k_{\perp f}$ is increased to $36$. Figure 6(a) clearly shows that, after the development of self-similar spectra (with exponents $-2$ for energy and $-3$ for GCH, together with a very weak negative GCH flux), a transient regime develops, corresponding to the onset of unstable modes, an effect that does not persist after turbulence has become sufficiently developed. This suggests that, if the decay instability can be a leading mechanism in the regime of weak turbulence (Voitenko Reference Voitenko1998b), its effect cannot be more than a transient in the strong turbulence regime. When present, the decay instability nevertheless significantly affects the global dynamics. As seen in figure 6(b), a sharp increase of the total energy and GCH takes place when the decay instability is acting. The instability can be identified at a time close to $t=7$. Later on, energy and GCH are still growing, but at lower rates, consistent with the existence of an inverse cascade. As expected, increasing $\beta _e$ for example to $\beta _e=50$ (not shown), the other parameters being fixed, $\chi _f$ increases to 1.1 and the decay instability is suppressed. It is also interesting to note that the cascade that develops at early time in this case is self-similar.

Figure 6.

Transverse energy and GCH spectra (a) and time variation of the total energy and GCH (b) for run $R_{36}^{2}$.

In order to highlight the presence of resonant interactions in run $R_{36}^{2}$, we display in figure 7(a) the colour-scale plot of the two-dimensional spectrum $|{\hat {\mu }}^+({\boldsymbol {k}})|^2$ in the plane $(k_x, k_z)$, showing decay instability at work (at $t \approx 8$). Since the forcing is statistically isotropic in the $(k_x, k_y)$ plane, this plot is sufficient to demonstrate the features of the full $k_x$, $k_y$ and $k_z$ dependence. We see numerous peaks, some of them associated with resonant triads. In fact, due to axisymmetry, each peak is essentially a circle penetrating the plane (see figure 7b). This is important when evaluating the resonance condition since it provides more freedom in choosing the corresponding $k_x$ and $k_y$ wavevector components. Of particular interest are triads with coordinates in the $(k_x,k_z)$ plane given by $(50,-16)$, $(25,20)$ and $(36,36)$ (which corresponds to the forcing). We indeed see that, for this triad, both $k_{1z} + k_{2z} = k_{fz}$ and $\omega _1 + \omega _2 = \omega _f$. In this context, what appears as $k_x$ should be viewed as $k_{\perp }$. Resonance conditions in the transverse plane then consist of five equations (three norms of transverse wavevectors and two resonance conditions) for six unknowns. Choosing arbitrarily one of them, we can easily construct a resonant triad. There is thus an infinite number of such resonant triads which can be viewed as defining a resonant manifold.

Figure 7.

Colour-scale plot of $\ln |{\hat {\mu }}^+|^2$ for run $R_{36}^2$ showing decay instability. (a) Cross-section in a $(k_x, k_z)$ plane. Because the forcing is statistically isotropic in the $(k_x, k_y)$ plane, this plot is sufficient to characterize the full $(k_x, k_y, k_z)$ space (see text). Among the visible peaks, some of them, such as $(50,-16)$, $(25,20)$ and $(36,36)$ (which corresponds to the forcing), are associated with resonant triads. (b) Three-dimensional view. The dominant peaks visible as dots in the $(k_x, k_y)$ plane appear as circles symmetric around the $k_z$ axis. Some circles are not visible since a transparency mask was applied to clearly visualize the most intense peaks.

5.3. Instability of the balanced state

The case where the $\mu ^+$ and $\mu ^-$ waves are driven in a balanced way, i.e. with a zero GCH injection rate, is illustrated with run $R_{13b}^2$, for which $\chi _f=0.26$, in figure 8. Figure 8(a) shows the energy spectra $E^{\pm }$ of the two counter-propagating waves. We observe that the balanced regime is unstable. The spectra of the counter-propagating waves, which identify at early times, separate and their difference increases with time. While positive GCH is transferred towards small scales and dissipated, negative GCH is transferred towards large scales, creating a finite global imbalance. As seen in figure 8(b), after some time, the global imbalance of run $R_{13b}^2$, as measured by the ratio ${\mathcal {E}}^+/{\mathcal {E}}^-$ of the energy of the counter-propagating waves, saturates to a value close to $0.5$, in contrast with the simulation with imbalanced driving for which this ratio keeps increasing with time. Interestingly, when the energy injection rate $\epsilon _E$ is increased by a factor of 256 (run $R_{13sb}^2$), leading to $\chi _f= 1.68$, balanced turbulence becomes stable. Another example where the balanced regime is unstable is provided by $R_{36b}^2$ for which $\chi _f= 0.16$. In this case, we first observe a decay instability before the amplitudes of the counter-propagating waves separate (not shown).

Figure 8.

(a) The $E^{\pm }(k_{\perp })$ spectra for run $R_{13b}^2$ resulting from a balanced driving. (b) Time evolution of ${\mathcal {E}}^-/{\mathcal {E}}^+$ for run $R_{13b}^2$ (red line) and of ${\mathcal {E}}^+/{\mathcal {E}}^-$ for run $R_{13}^2$ (green line). The balanced state is unstable and evolves to a regime where ${\mathcal {E}}^-/{\mathcal {E}}^+ \simeq 2$.

When for the same $\beta _e$, and the same energy injection rate, the wavenumber $k_f$ is reduced to $1.3$ (run $R_{1.3b}^2$), thus increasing the value of $\chi _f$ to $0.85$, the instability is suppressed. It is recovered when reducing $\beta _e$ to $0.2$ (run $R_{1.3b}^{0.2}$), which has the effect of reducing $\chi _f$ to $0.38$. When, keeping $\beta _e=0.2$ and $k_f=1.3$, the energy injection rate is increased by a factor of $16$ (run $R_{1.3sb}^{0.2}$), $\chi _f = 1.68$ and the instability is again suppressed. A more quantitative study would be required to determine a precise instability threshold and to confirm the prediction that $\chi _f$ is the only parameter that controls the instability of the balanced state.