Direct kinetic Alfvén wave energy cascade in the presence of imbalance

Abstract

A two-field Hamiltonian gyrofluid model for kinetic Alfvén waves retaining ion finite Larmor radius corrections and parallel magnetic field fluctuations is used to study direct turbulent cascades from the magnetohydrodynamic to the sub-ion scales. For moderate energy imbalance and weak enough magnetic fluctuations, the spectrum of the transverse magnetic field and that of the most energetic wave display a steep transition zone near the ion scale, while the parallel transfer (and thus the parallel dissipation) remains weak. In this regime, the perpendicular flux of generalized cross-helicity displays a significant decay past the ion scale, while the perpendicular energy flux remains almost constant. A phenomenological model suggests that the interactions between co-propagative waves present at the sub-ion scales can play a central role in the development of a transition zone in the presence of a helicity barrier.

1 Introduction

Alfvén wave (AW) turbulence, which plays an important role in magnetized collisionless plasmas such as the solar wind (Bruno & Carbone 2013, 2016), is often imbalanced in the sense that the energy carried by the waves propagating in the forward and backward directions relative to the ambient field have significantly different magnitudes (Tu, Marsch & Thieme 1989; Lucek & Balogh 1998; Wicks et al. 2013). The degree of imbalance is dependent on the type of wind (Tu, Marsch & Rausenbauer 1990) and also on the distance from the Sun (Roberts et al. 1987; Marsch & Tu 1990; Chen et al. 2020). For example, in a classical model of solar wind, outgoing AWs are emitted at the Sun's surface, as a consequence of reconnection processes in the chromospheric magnetic network, while ingoing waves result from reflection of these waves on density gradients (Perez & Chandran 2013; Chandran & Perez 2019), or from parametric decay instability (Malara, Primavera & Veltri 2001). In the case of low-frequency highly transverse perturbations, corresponding to the assumptions underlying the gyrofluid model derived in Passot, Sulem & Tassi (2018), imbalance is associated with a non-zero value of an invariant referred to as generalized cross-helicity (GCH), which identifies with the cross-helicity at the magnetohydrodynamic (MHD) scales and with the magnetic helicity at sub-ion scales. As cross-helicity cascades to small scales and magnetic helicity to large scales (Miloshevich et al. 2021), the question arises as to the behaviour of the GCH flux arriving from large scales to the transition zone between MHD and sub-ion scales. Two publications recently addressed this issue. One uses a gyrofluid model (in the limit of small values of the electron beta parameter) which isolates the AW dynamics (Meyrand et al. 2021). The other is based on a hybrid kinetic particle-in-cell simulation which enables the coupling of low-frequency quasi-transverse AWs with quasi-parallel high-frequency ion cyclotron waves (ICWs) (Squire et al. 2022). The role of ICWs and of the dissipation associated with the proton cyclotron resonance at the spectral break separating the MHD and the sub-ion ranges was indeed stressed, using the data from the Messenger and Wind spacecraft missions (Telloni, Bruno & Trenchi 2015; Woodham et al. 2018).

In the regimes studied by Meyrand et al. (2021), the nonlinear parameter, defined as the ratio of the inverse nonlinear time to the AW frequency, is typically above unity. The dynamics then rapidly develops a significant transfer of energy towards large parallel wavenumbers, resulting in the breaking of the asymptotics underlying gyrofluid models. Despite this shortcoming, the gyrofluid simulation reveals the existence of a helicity barrier which depletes the energy flux towards the small perpendicular scales, leading to a transfer in the parallel direction at the ion scale, together with a steepening of the energy (and magnetic) spectra. Energy accumulates at large scales and imbalance increases until saturation occurs under the action of a strong parallel dissipation. In contrast to usual fluid turbulence, the energy level that is reached depends on the magnitude of the parallel viscosity coefficient. In this framework, dissipation, which plays a central role in the dynamics, is not associated with any realistic physical phenomenon. The above scenario was then validated by the kinetic simulation of Squire et al. (2022), where the dissipation originates from ion cyclotron resonance. In this simulation, performed at $\beta _e =0.3$, the level of fluctuations remains moderate because cyclotron resonance is rapidly triggered. In fact, an argument based on critical balance, together with the requirement that $k_{\|} d_i=1$ (where $d_i$ is the ion inertial length) at perpendicular scales such that $k_{\perp } \rho _i=1$ ($\rho _i$ denoting the ion Larmor radius), leads to a saturation amplitude scaling like $\beta _i^{1/2}$. As a consequence, when $\beta _e$ is taken of order unity (for comparable ion and electron temperatures), one may expect that high-amplitude fluctuations should develop before a stationary regime can establish. Recent observations using Parker Solar Probe support the results of the above simulation, with evidence of the coincidence between the presence of ICWs and of magnetic spectra displaying an ion transition zone associated with a polarity change, as expected when AWs trigger the formation of ICWs (Bowen et al. 2021).

Transition zones are commonly observed in the solar wind at 1 AU, where the plasma beta parameter is of order unity or larger (see e.g. Huang et al. (2021) and references therein). When the amplitude of the fluctuations is moderate (such as in the fast wind), the above argument suggests that ICWs can hardly be excited at the perpendicular ion scale. In this context, an alternative mechanism should be at work, like for example the effect of co-propagating AW interactions that become relevant at ion scales (Voitenko & De Keyser 2016; Gogoberidze & Voitenko 2020). In order to address this question, it is appropriate to use a two-field gyrofluid model (Passot et al. 2018), where driving and dissipation have been supplemented, which only retains the AWs, but requires that the fluctuation amplitudes remain small and the parallel dynamics sub-dominant (conditions that can be checked a posteriori). We chose to prescribe the ‘energy imbalance’ (referred to as imbalance in the following) by maintaining constant the amplitudes of the largest modes retained in the simulations, a driving process classically used in hydrodynamics turbulence (Siggia & Patterson 1978). It turns out that, in the quasi-stationary phases of the simulation, the GCH and energy injection rates are roughly constant (equal to the total dissipation rates), the ‘injection imbalance’, measured by their ratio, resulting from the dynamics.

The paper is organized as follows. Section 2 describes the model and the numerical set-up. Section 3 discusses the gyrofluid simulations, showing that, for moderate imbalance, the spectrum of the magnetic fluctuations (and also that of the most energetic wave) displays a steep transition zone near the ion scale, while the parallel transfer (and thus the parallel dissipation) remains weak. The helicity barrier is mostly conspicuous on the perpendicular GCH flux, which displays a significant decay past the ion scale, while the transverse energy flux remains almost constant. In the situation where parallel energy and cross-helicity dissipations are negligible, the effect of co-propagating wave interactions seems to be the only alternative mechanism at the origin of the transition zone within the context of the present model. Section 4 discusses this mechanism using a Leith-type phenomenological spectral model. Section 5 presents the conclusions.

2 Equations and numerical set-up

Let us consider a situation where the plasma dynamics results from weak perturbations of a homogeneous equilibrium state characterized by a density $n_0$, isotropic ion and electron temperatures $T_{0i}$ and $T_{0e}$ (characterized by their ratio $\tau = T_{0i}/T_{0e}$) and subject to an ambient magnetic field of amplitude $B_0$ along the $z$ direction. The various characteristic scales are conveniently measured in terms of the sonic Larmor radius $\rho _s = c_s/\varOmega _i$, where $c_s = \sqrt {T_{0e}/m_i}$ is the speed of sound and $\varOmega _i= eB_0/(mc)$ the ion gyrofrequency. In particular, the ion Larmor radius $\rho _i =v_{th\,i}/\varOmega _i$ (where the ion thermal velocity is given by $v_{th , i}=(2T_{0i}/m_i)^{1/2}$) and the ion inertial length $d_i=v_A/\varOmega _i$ (where $v_A=B_0/(4{\rm \pi} n_0 m_i)^{1/2} =c_s \sqrt {2/\beta _e}$ is the Alfvén velocity) satisfy $\rho _i = \sqrt {2\tau } \rho _s$ and $d_i = \sqrt {2/{\beta _e}} \rho _s$.

The two-field gyrofluid model is derived from the full gyrofluid model of Brizard (1992) under specific scaling assumptions. Such a simple reduced fluid model does not retain physical phenomena like kinetic dissipation, but also misses nonlinear couplings with other modes. Discarding these effects, one here focuses on the nonlinear dynamics of AWs, or kinetic AWs. Fully kinetic simulations will ultimately be necessary to validate the findings. Like the parent model, the model is not restricted to the small-$\beta _e$ regime, unless physics at the electron inertia scale $d_e$ is retained (in the absence of electron finite Larmor radius terms in the two-field equations, the electron Larmor radius must indeed remain smaller than $d_e$). More precisely, leaving aside issues related to kinetic damping, the model reproduces reduced MHD in the limit of large transverse scales, and electron reduced MHD at sub-ion scales, both systems describing AWs or kinetic AWs, uncoupled to slow waves, with no restriction on $\beta _e$ (especially since the parallel magnetic field component is retained). Only at the ion scales, where the system includes a dispersive term associated with the electron pressure gradient, can a coupling between AWs and slow waves become relevant (see Appendix E of Schekochihin et al. 2009). We here neglect the electron inertia, so the model covers a spectral range extending from the MHD to the sub-ion scales, within an anisotropic scaling where the transverse scales are much smaller than the parallel ones.

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 (1978). 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. 2021).

3 The GCH barrier

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. 2020) 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. 2017).

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. (2021) 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. (2021), 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. (2021). 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$).¹

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. (2021).

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. 2010), 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 (2016) and Gogoberidze & Voitenko (2020). 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 Modelling of the transition zone

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 (2019), and then reproduced in Miloshevich et al. (2020) 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 (2019), following Voitenko & De Keyser (2016), 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 (2019) 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 (2016), 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 (2019), (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.² 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 2011).

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 2019). When $d=3$ (a value compatible with kinetic theory, used by Voitenko & De Keyser (2016)), 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.

5 Conclusion

Using a two-field gyrofluid model retaining AWs only, we have observed that, for $\beta _e$ of order unity or larger and a weak or moderate energy imbalance, a transition zone can exist in the absence of a significant parallel energy dissipation, in contrast with the mechanism based on the dissipation induced by ion cyclotron resonance, recently proposed by Squire et al. (2022). This transition zone is found to be steeper at larger imbalance, in line with Parker Solar Probe observations (Huang et al. 2021). In these simulations, a helicity barrier also forms, associated with a decay of the perpendicular CGH flux at the ion scales, while the corresponding energy flux is almost scale-independent. Such simulations are consistent with the possible effect of the interactions between co-propagating AWs. The efficiency of this mechanism is supported by the analysis of a phenomenological model based on diffusion equations in the spectral space, in the spirit of the Leith model (Leith 1967) for hydrodynamic turbulence. It turns out that, in the presence of co-propagative waves, and only in this case, the energy barrier, modelled by a decay of the GCH flux with the transverse wavenumber, leads to a spectral transition zone. As a consequence of weak or negligible parallel dissipation, the injected energy accumulates and a relatively long time is needed for turbulence to reach a steady state. As a consequence, in the context of a more realistic model, one may suspect that before an ICW can be excited, other dissipative effects can come into play, among them Landau damping or stochastic heating. Kinetic simulations for values of the beta parameters of order unity or larger are necessary to address these issues.