Modeling Imbalanced Collisionless Alfvén Wave Turbulence with Nonlinear Diffusion Equations

An important contribution to the solar wind (SW) dynamics originates from nonlinear interactions between counter-propagating Alfvén waves. Outgoing waves are emitted at the Sun's surface as a consequence of mechanisms such as reconnection in the chromospheric magnetic network (see, e.g., the recent review by Marsch 2018), which leads to injection of Alfvénic waves at the base of the fast SW (McIntosh 2012). Ingoing waves are created by reflection on density gradients (Chandran & Hollweg 2009; Chandran & Perez 2019) and by velocity shear (Roberts et al. 1992; Breech et al. 2008), or result from parametric decay instability (Viñas & Goldstein 1991; Zanna et al. 2001). The energies carried by these counter-propagating waves are usually unequal (Tu et al. 1989; Lucek & Balogh 1998; Wicks et al. 2013), with a degree of "imbalance" depending on the type of wind (Tu et al. 1990; Bruno et al. 2014, 2017; D'Amicis et al. 2019) and also on the distance from the Sun (Roberts et al. 1987; Marsch & Tu 1990). Accurate in situ observations are now available for the turbulent energy cascades of Alfvén waves (AWs) or, at the sub-ion scales, kinetic Alfvén waves (KAWs; see, e.g., Goldstein et al. 2015; Kiyani et al. 2015 for reviews). A precise understanding of how imbalance affects these cascades is, however, required in order to predict how much turbulence can heat and accelerate the SW plasma (Chandran et al. 2010; Cranmer et al. 2015; Mallet et al. 2019). Imbalance can also affect cosmic-ray scattering efficiency, as for example discussed in Weidl et al. (2015).

While in the MHD range, imbalanced Alfvénic turbulence has been extensively studied both theoretically and numerically (see Chen 2016 for a review), its dynamics at sub-ion scales remains largely unexplored. Due to the computational cost of 3D imbalanced kinetic simulations, one is led to develop asymptotic models, isolating the AW dynamics within a spectral range extending from the MHD scales (larger than ion and sonic Larmor radii) to the sub-ion scales, assuming a proton–electron homogeneous plasma subject to a strong ambient magnetic field. Such a description is provided by a Hamiltonian two-field gyrofluid model retaining ion finite Larmor radius (FLR) corrections, parallel magnetic fluctuations, and electron inertia (Passot et al. 2018). It takes the form of two dynamical equations for the electron-gyrocenter density and the parallel magnetic potential, from which the electrostatic potential and the parallel magnetic fluctuations are easily derived. Numerical simulations of the two-field gyrofluid are in progress. Here we concentrate on a phenomenological reduction of this model (Passot & Sullem 2019) in the form of nonlinear diffusion equations in Fourier space for the transverse spectra of the two conserved quantities, energy and generalized cross-helicity (GCH). Concentrating on the effect of imbalance at scales larger than the electron skin depth d_e, electron inertia is neglected in most of the simulations, which all address the strong turbulence regime.

The model involves drastic simplifications, leaving for future study the influence of effects such as inhomogeneities and wind expansion. In some of the simulations, we nevertheless retain Landau damping, which was shown to affect the sub-ion inertial range of balanced AW turbulence (Passot & Sulem 2015; Sulem et al. 2016), but neglect ion cyclotron damping and heating of the medium, considered for example by Cranmer & van Ballegooijen (2003), in spite of their potential effect on the dynamics.

In this section we provide a purely phenomonological derivation of the diffusion model discussed in Passot & Sulem (2019) for the time evolution of the transverse spectra $E({k}_{\perp },t)$ and ${E}_{C}({k}_{\perp },t)$ of energy and GCH respectively. The latter quantities are related to the energy spectra ${E}^{\pm }({k}_{\perp },t)$ of the counter-propagating waves by $E({k}_{\perp },t)={E}^{+}({k}_{\perp },t)+{E}^{-}({k}_{\perp },t)$ and ${E}_{C}({k}_{\perp },t)=({E}^{+}({k}_{\perp },t)-{E}^{-}({k}_{\perp },t))/{v}_{\mathrm{ph}}({k}_{\perp })$, where the (parallel) Alfvén phase velocity v_ph is a function of ${k}_{\perp }$.

A nonlinear diffusion model suitable for three-wave interactions (dominant for KAW turbulence) and preserving the existence of absolute-equilibrium solutions should take the form (Thalabard et al. 2015)

$$\begin{eqnarray}&&\displaystyle \frac{\partial }{\partial t}\displaystyle \frac{{E}^{\pm }({k}_{\perp })}{2}=\displaystyle \frac{\partial }{\partial {k}_{\perp }}\left[{D}^{\mp }({E}^{+},{E}^{-},{k}_{\perp })\displaystyle \frac{\partial }{\partial {k}_{\perp }}\displaystyle \frac{{E}^{\pm }({k}_{\perp })}{{k}_{\perp }}\right]+{X}^{\pm },\end{eqnarray} \tag{ 1 }$$

where the corrective terms X^± are necessary because the energies of the forward and backward propagating waves are not conserved independently in the presence of dispersion. Due to energy conservation, it is clear that X⁺ = −X⁻. The equation for the GCH spectrum then reads

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{\partial }{\partial t}\displaystyle \frac{{E}_{C}}{2} & = & \displaystyle \frac{2{X}^{+}}{{v}_{\mathrm{ph}}}+\left(\displaystyle \frac{\partial }{\partial {k}_{\perp }}+\displaystyle \frac{1}{{v}_{\mathrm{ph}}}\displaystyle \frac{\partial {v}_{\mathrm{ph}}}{\partial {k}_{\perp }}\right)\\ & & \times \left\{\sum _{r=\pm 1}{\left(-1\right)}^{\tfrac{r-1}{2}}\displaystyle \frac{{D}^{(-r)}}{{v}_{\mathrm{ph}}}\displaystyle \frac{\partial }{\partial {k}_{\perp }}\left(\displaystyle \frac{{E}^{(r)}}{{k}_{\perp }}\right)\right\}.\end{array}\end{eqnarray} \tag{ 2 }$$

The invariance of the GCH, which implies that the right-hand side of (2) should be in a conservative form, allows us to find X⁺ by the condition that only the first term in the round parentheses survives. Furthermore, dimensional analysis of the terms involved in Equation (1) allows us to estimate ${D}^{\mp }\sim {k}_{\perp }^{3}/{\tau }_{\mathrm{tr}}^{\mp }$ in terms of the transfer time ${\tau }_{\mathrm{tr}}^{\mp }={({\tau }_{\mathrm{NL}}^{\pm })}^{2}/{\tau }_{{\rm{w}}}^{\mp }$. The case with Landau damping is treated in Appendix A. In its absence, ${D}^{\mp }({E}^{\mp },{k}_{\perp })=C^{\prime} {k}_{\perp }^{6}\,{v}_{\mathrm{ph}}{E}^{\mp }/{\widetilde{k}}_{\parallel }^{\mp }$, where ${\widetilde{k}}_{\parallel }^{\mp }$ measures the inverse parallel correlation length of ± eddies or wave packets with perpendicular wavenumber k_⊥. Here, we used arguments from Passot & Sulem (2019) to fix the nonlinear time ${\tau }_{\mathrm{NL}}^{\pm }={({k}_{\perp }^{3}{v}_{\mathrm{ph}}^{2}{E}^{\pm })}^{-1/2}$, consistently with the imbalanced strong MHD turbulence (see Lithwick et al. 2007), and also wrote ${\tau }_{{\rm{w}}}^{\mp }={\widetilde{k}}_{\parallel }^{\mp }{v}_{\mathrm{ph}}$. The parameter C' can be scaled out in the absence of Landau damping but in its presence the former affects the location of the transition to the exponential decay (Passot & Sulem 2015). In the following, it is chosen to equal unity. Retaining only strongly local interactions, the model ignores coupling between co-propagating waves since they have nearly the same velocity and can hardly overtake one another and interact. Interactions between such waves with comparable but not quasi-equal wavenumbers were considered by Voitenko & Keyser (2016), but they were shown not to significantly affect the dynamics in the context of the present model (Passot & Sulem 2019).

Assuming that both waves undergo a strong cascade ${\tau }_{\mathrm{NL}}^{\pm }\sim {\tau }_{{\rm{w}}}^{\mp }$, we get ${\widetilde{k}}_{\parallel }^{\pm }={({k}_{\perp }^{3}{E}^{\pm })}^{1/2}$. However in MHD it is expected that in the imbalanced case, the + wave (which by definition is more energetic) will undergo a weaker cascade as suggested by Chandran (2008). This effect can be modeled by changing the definition of the parallel wavenumber ${\widetilde{k}}_{\parallel }^{-}$ of the "–" wave subject to interactions with the "+" wave. This suggests the parameter χ should be introduced along with the following ansatz:

$$\begin{eqnarray}&&{\widetilde{k}}_{\parallel }^{(r)}={\left({k}_{\perp }^{3}{E}^{(r)}\right)}^{1/2}{\left({E}^{+}/{E}^{-}\right)}^{(1-r)\chi /4}.\end{eqnarray} \tag{ 3 }$$

This expression reproduces the different models found in the literature for imbalanced MHD turbulence, depending on the free exponent χ. The value χ = 0 corresponds to the model of Lithwick et al. (2007),¹ where both waves are in a strong turbulent regime and satisfy the critical balance condition. In contrast, χ = 1 reproduces the model of Chandran (2008), where the transfer time of the stronger wave obeys a weak turbulence phenomenology. The value χ = 1/4 corresponds to the model of Beresnyak & Lazarian (2008). Numerical simulations of imbalanced MHD by these authors seem to favor a value of χ larger than zero but also strictly smaller than unity. Differently, in the presence of dispersion, any value of χ smaller than one leads to an unphysical instability at the dissipation scale (Passot & Sulem 2019). The simulations without Landau damping described below (except the purely MHD ones) are thus performed with χ = 1. The stability range is extended to χ ≥ 0.5 when Landau damping is retained.

Finally, the set of reduced equations are cast as

$$\begin{eqnarray}&&\displaystyle \frac{\partial }{\partial t}\displaystyle \frac{E}{2}=\displaystyle \frac{\partial }{\partial {k}_{\perp }}\left\{{k}_{\perp }^{6}V\,\sum _{r=\pm 1}\displaystyle \frac{{E}^{(-r)}}{{\widetilde{k}}_{\parallel }^{(-r)}}\displaystyle \frac{\partial }{\partial {k}_{\perp }}\left(\displaystyle \frac{{E}^{(r)}}{{k}_{\perp }}\right)\right\},\end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial }{\partial t}\displaystyle \frac{{E}_{C}}{2}=\displaystyle \frac{\partial }{\partial {k}_{\perp }}\left\{{k}_{\perp }^{6}\sum _{r=\pm 1}\displaystyle \frac{{\left(-1\right)}^{\tfrac{r-1}{2}}{{VE}}^{(-r)}}{{v}_{\mathrm{ph}}{\widetilde{k}}_{\parallel }^{(-r)}}\displaystyle \frac{\partial }{\partial {k}_{\perp }}\left(\displaystyle \frac{{E}^{(r)}}{{k}_{\perp }}\right)\right\}.\end{eqnarray} \tag{ 5 }$$

In the absence of Landau damping V = v_ph (which is the only quantity in the model that incorporates the kinetic effects). It is constant at wavenumbers that are small compared with the smallest of the inverse ion Larmor radius ${\rho }_{i}^{-1}={(\sqrt{2\tau }{\rho }_{s})}^{-1}$ (where $\tau ={T}_{0i}/{T}_{0e}$ is the ion to electron temperature ratio at equilibrium) and the inverse sonic Larmor radius ${\rho }_{s}^{-1}={(\sqrt{{m}_{i}/{T}_{{oe}}}{{\rm{\Omega }}}_{i})}^{-1}$ (where m_i is the ion mass and Ω_i the ion gyrofrequency), used as the wavenumber unit (with ${{\rm{\Omega }}}_{i}^{-1}$ taken as the time unit). It grows linearly (v_ph ≈ α k_⊥) at smaller scales down to ${d}_{e}={(2/{\beta }_{e})}^{1/2}\delta {\rho }_{s}$ (with δ² denoting the electron to proton mass ratio), where saturation occurs due to electron inertia.

We treat Landau damping (see Appendix A) by adding dissipative terms −γE and −γE_C in Equations (4)–(5), respectively, where $\gamma ({k}_{\perp },{\widetilde{k}}_{\parallel })$ is the AW Landau dissipation rate (dominated by electrons) of transverse and parallel wavenumbers k_⊥ and ${\widetilde{k}}_{\parallel }=\max ({\widetilde{k}}_{\parallel }^{+},{\widetilde{k}}_{\parallel }^{-})$. It also affects (mostly through the ions) the transfer time (Passot & Sulem 2015) and thus modifies the function V entering Equations (4)–(5), which becomes $V={v}_{\mathrm{ph}}^{2}/({v}_{\mathrm{ph}}+\mu \sqrt{2\tau })$ where μ is a numerical constant (in the following, μ = 3).

Equations (4)–(5) were solved numerically using a finite difference scheme both in k_⊥ (with a logarithmic discretization) and t (using a forward Euler scheme with adaptive time step), modifying a code developed for weak turbulence of gravitational waves (Galtier et al. 2019). We also introduced different grids for the fields and the fluxes, with linear interpolation on the fields. Their coupling ensures better stability.

Simulations were performed for β_e = 2, τ = 1 (case I) typical of the SW at 1 au, and β_e = 0.04, τ = 10 (case II), more suitable for regions closer to the Sun (Roytershteyn et al. 2019). For small β_e, d_e significantly exceeds the electron Larmor radius ${\rho }_{e}=\sqrt{2}\delta {\rho }_{s}$, thus permitting electron inertia to be retained, while electron FLR corrections are neglected. Simulations with and without Landau damping were performed, the latter regime being of interest for comparison with analytical predictions. In this case, hyperdiffusive terms of the form ν k_⊥⁸E(k_⊥, t) and $\nu {k}_{\perp }^{8}{E}_{C}({k}_{\perp },t)$, are supplemented in Equations (4)–(5), with a coefficient ν depending on the resolution and the parameters. The system is driven at large scales (injection wavenumber k₀), dissipation taking place at small scales (dissipation or pinning wavenumber k_d between 10¹ and 10⁴). Driving is performed either through boundary conditions by prescribing the spectra E^± or the fluxes η and of energy and GCH at the smallest wavenumber k_min = k₀ (within the range 10⁻³ to 10⁻²), or differently through injection terms of the form $\varepsilon f({k}_{\perp })$ and $\eta f({k}_{\perp })$ supplemented in the right-hand side of Equations (4)–(5), f denoting a function with a compact support located near the injection wavenumber such that ${\int }_{0}^{\infty }f({k}_{\perp }){{dk}}_{\perp }=1$. Initial conditions consist of a spectral bump in the MHD range near the smallest retained wavenumber.

When, in the absence of Landau damping, Equations (4)–(5) are integrated in the strong turbulence regime with prescribed boundary values E^± (k₀) and ${k}_{\perp }^{8}$ dissipation terms, we observe a direct transfer of energy and GCH, with the development of power-law spectra that progressively extend to the small scales, up to the moment when the dissipation scale is reached. Although the MHD-range spectrum that forms during this transient is very close to the ${k}_{\perp }^{-5/3}$ stationary solution, the sub-ion spectra first develop a balanced regime (as a result of a faster transfer of energy than of GCH) with a spectrum steeper than ${k}_{\perp }^{-8/3}$ before converging, while imbalance develops, toward an approximately −7/3 steady state when imbalance is weak (not shown). Nonlinear diffusion models for finite-capacity systems (no ultraviolet divergence of the energy) often exhibit such transient anomalous power-law spectra (Thalabard et al. 2015). The case of weak KAW turbulence has recently been addressed by David & Galtier (2019).

When the boundary conditions in k₀ consist of prescribing energy and GCH η fluxes, the imbalance between E⁺ and E⁻ at large scales is found to strongly depend on the small-scale extension of the dispersive range, as seen in Figure 1 (top). This is to be compared with Figure 1 (bottom) where the phase velocity is assumed constant (nondispersive MHD). Both χ = 1 (for comparison with the dispersive case) and the more realistic value χ = 0.25 are displayed. In both panels, the spectra are plotted at times at which a stationary regime with constant (positive) energy and GCH fluxes have established. In these simulations (case I with δ = 0), we varied the dissipation wavenumber by changing the hyperdiffusivity coefficient. In the MHD regime, E^± behave as power laws until the pinning scale is approached. Conversely, in the dispersive case, the spectra approach each other exponentially (consistent with the linear variation of $\phi ({k}_{\perp })=(1/2)\mathrm{ln}({E}^{+}({k}_{\perp })/{E}^{-}({k}_{\perp }))$ when $\alpha \eta /\epsilon \ll 1$ (see Appendix B). It follows that, in contrast with usual hydrodynamic turbulence, large-scale quantities depend on the dissipation scale ${k}_{d}^{-1}$. This effect is much stronger in the presence of dispersion where ${E}^{+}({k}_{0})/{E}^{-}({k}_{0})$ varies exponentially with k_d, while it scales like a power law in MHD when keeping χ = 1 (see Equations (17) and (13), respectively). For MHD with χ = 0.25, the spectra differ by less than 1% from the case χ = 0 for which there is no dependence on k_d (see Equation (20)). To understand heuristically the behavior in the dispersive case, we must note that (i) the dissipative scales for E and E_C are the same since the equation for the GCH spectrum is linear in E_C, and (ii) the E_C spectrum steepens at the dispersive scale even more than the energy spectrum since $| {E}_{C}({k}_{\perp })| \leqslant E({k}_{\perp })/{v}_{\mathrm{ph}}({k}_{\perp })$. As a consequence, in order for the GCH dissipation $\nu \int {k}_{\perp }^{8}{E}_{C}{{dk}}_{\perp }$ to match the injected GCH rate η prescribed at k₀, the magnitude of E_C must be larger, and consequently the imbalance at large scales is enhanced compared with the MHD problem in the same setting.

Zoom In Zoom Out Reset image size

Difference between the dispersive (GYRO) and nondispersive (MHD) cases is also seen on Figure 2, which displays the imbalance ${E}^{+}({k}_{0})/{E}^{-}({k}_{0})$ at the outer wavenumber ${k}_{0}={k}_{\min }$ (located in the MHD range) versus the parameter $a\equiv {v}_{\mathrm{ph}}({k}_{0})\eta /\epsilon \approx \sqrt{2/{\beta }_{e}}\eta /\epsilon$ (MHD) or $a\equiv \alpha \eta /\epsilon$ (GYRO), with α defined in Appendix B, in cases I and II. Graphs corresponding to different β_e collapse on the same curve, with an almost perfect agreement in MHD. Note that, in the GYRO simulations, the assumption ${v}_{\mathrm{ph}}=\alpha {k}_{\perp }$ is only approximate, especially close to the ion scale. We furthermore observed that changing ε while keeping η/ε constant has no effect.

Zoom In Zoom Out Reset image size

Influence of the degree of imbalance in the presence of Landau damping is considered in Figure 3 in case I with δ = 0, where, for comparison, a simulation involving hyperdiffusion is also presented. The top panel displays the energy spectra E^±, while the bottom panel shows the transverse magnetic spectrum ${E}_{{B}_{\perp }}$ (see Appendix C). With hyperdiffusion, the ${k}_{\perp }^{-5/3}$ MHD spectrum is continued at sub-ion scales by the classical ${k}_{\perp }^{-7/3}$ range, the degree of imbalance significantly decreasing only near the pinning wavenumber. When Landau damping is retained, the degree of imbalance decreases with the scale, the more so when χ is closer to 1 (visible in the top panel when comparing the runs with ${E}^{+}({k}_{0})/{E}^{-}({k}_{0})=100$ for χ = 0.6 and 0.5). In the sub-ion range, the steepening of ${E}_{{B}_{\perp }}$ displayed in the bottom panel increases with the degree of imbalance (a more pronounced effect when χ is smaller). Such a steepening is often observed as a transition range in the SW at 1 au, depending on the fluctuation power (Bruno et al. 2014), and sometimes associated with proton Landau damping (Sahraoui et al. 2010) and imbalance degree or Alfvénicity (Bruno et al. 2014, 2017; D'Amicis et al. 2019). The present model suggests that both effects are to act simultaneously. An alternative mechanism for steepening of the spectrum related to reconnection is suggested by Vech et al. (2018).

Zoom In Zoom Out Reset image size

Figure 4 displays ${E}_{{B}_{\perp }}$ in case II with δ² = 1/1836. For this relatively small β_e, an even stronger spectral steepening is observed for a large imbalance. Furthermore, the multiplicative factor in the estimate of E(k_⊥) in terms of ${E}_{{B}_{\perp }}$ (k_⊥) is responsible for a steepening of the latter at scales smaller than d_e, an effect visible in particular in the simulation with hyperdiffusivity, χ = 1, and ${E}^{+}({k}_{0})/{E}^{-}({k}_{0})=100$. A ${k}_{\perp }^{-11/3}$ spectrum, classical for balanced IKAW turbulence (Chen & Boldyrev 2017; Passot et al. 2018; Roytershteyn et al. 2019), is still observed for this level of imbalance. Spectra, obtained with Landau damping and electron inertia, are also displayed in this figure: for χ = 0.6 with ${E}^{+}({k}_{0})/{E}^{-}({k}_{0})=1$ or ${E}^{+}({k}_{0})/{E}^{-}({k}_{0})=100$, showing that the spectra get steeper as the imbalance increases, and for χ = 1 with ${E}^{+}({k}_{0})/{E}^{-}({k}_{0})=100$, showing that the steepening is more pronounced as χ is decreased. Such a steepening (although more moderate) is also observed in 3D fully kinetic simulations (Grošelj et al. 2018).

Zoom In Zoom Out Reset image size

We consider a situation where a stationary turbulence is affected by small-scale injection of KAWs resulting, e.g., from magnetic reconnection (Chaston et al. 2005; Liang et al. 2016; Shi et al. 2019 and references therein). We focus on the case where the injection is imbalanced. This setting differs from that used for studying a possible inverse cascade of GCH by driving an initially zero solution. Such a cascade, predicted to exist at sub-ion scales using absolute-equilibrium arguments (Passot et al. 2018), and necessarily involving inverse energy transfer, is not captured by the diffusion model and is currently investigated via direct simulations of the gyrofluid model.

In the following, the initial stationary regime is obtained with a large-scale driving at k₀ = 3.2 (larger than k_min) and a k⁸ hyperdiffusivity ν = 10⁻¹⁶ in GYRO and ν = 10⁻¹⁸ in MHD runs. An injection of E⁺ energy only (maximal imbalance) is then provided at a wavenumber k_i = 51.2 at a rate smaller than or equal to the one imposed at large scales. Stationary solutions are obtained by using a k⁻⁸ hypodiffusivity equal to 10⁻⁷ for MHD and 10⁻⁶ for GYRO. The results presented below pertain to the case where the initial turbulence is balanced but later becomes imbalanced as a result of additional forcing at k_i.

Let us first address the regime where the driving takes place in the MHD range for χ = 0. We see in Figure 5 (top) that some imbalance is generated at wavenumbers larger than k_i but not at smaller wavenumbers, indicating the absence of inverse transfer of GCH. This point can be understood by solving Equation (19) (right) for $\phi =(1/2)\mathrm{log}({E}^{+}({k}_{\perp })/{E}^{-}({k}_{\perp }))$ with the initial condition taken at k₀ ≪ k_⊥ < k_i. For the parameters of the simulation, the initial condition for ϕ is of order unity and the constant C turns out to be very large, leading to a fast decrease of ϕ − ψ (where ψ = 2 tanh⁻¹ a) with ${k}_{\perp }$. We thus find that for wavenumbers approaching k₀, the imbalance remains unchanged, equal to the value given by Equation (20). Numerical simulations show that the absence of inverse transfer of GCH still holds for larger values of χ, including in particular 0.25.

Zoom In Zoom Out Reset image size

When injection takes place in the dispersive range (with χ = 1), assuming that αη/ ≪ 1, Equation (16) indicates that ϕ is almost independent of k_⊥, and thus that the imbalance created at k_i extends to large scales, as exemplified in Figure 5 (bottom). This evolution results from the development of a transient negative GCH flux at ${k}_{\perp }\lt {k}_{i}$ due to the small-scale forcing, leading to a decrease of the total GCH flux and to an enhancement of imbalance that saturates when the flux recovers its original value (consistent with the GCH flux ${\eta }_{0}$ injected at k₀). In the SW, dispersive scales are always present, making the possibility of such an inverse transfer relevant. It would be of interest to investigate whether the process by which reconnection events can generate an inverse flux toward larger scales (Franci et al. 2017) can also generate imbalance in the case where KAWs are generated at the reconnection sites. The present modeling could also be useful to analyze the recently predicted cascade reversal at d_e in a reduced two-fluid model (Miloshevich et al. 2018) and a 3D extended magnetohydrodynamic model (Miloshevich et al. 2017).

This Letter provides an analysis of the influence of imbalance between the energies of counter-propagating AWs, on the dynamics of a collisionless plasma. It aims to contribute to the understanding of regimes encountered in the SW, particularly in regions close to the Sun explored by space missions such as Parker Solar Probe or Solar Orbiter.

The main results can be summarized as follows. (i) The imbalance produced by large-scale injection of GCH at a prescribed rate is enhanced by wave dispersion. (ii) A steep range in the spectrum of the transverse magnetic fluctuations, consistent with the transition region reported in SW observations (Sahraoui et al. 2010), develops at the sub-ion scales, under the combined influence of Landau damping and strong imbalance, an effect enhanced when β_e is decreased. We conjecture that the simulation results can be more specifically related to the observations inside the trailing edge, which is characterized by the highest level of Alfvénicity, i.e., imbalance (Bruno et al. 2014). Existence of a shallower spectrum at smaller scales is then expected to originate from SW regions that are less imbalanced and more energetic at these scales. (iii) Under some conditions the system develops an inverse transfer of imbalance when imbalanced forcing takes place at small scales in an already fully developed turbulence.

Future works include the study of the parent two-field gyrofluid models that is in particular expected to address the question of the characteristic nonlinear timescale in imbalanced turbulence, evaluate the assumption of strongly local interactions, and investigate the role of KAW decay instability. The influence on the global dynamics of the coupling of the AWs with the slow modes, important at small β_e as they can generate large-scale parametric decay instabilities, will be studied using an extension of the present gyrofluid including both kinds of waves.

Computations have been done on the "Mesocentre SIGAMM" machine, hosted by Observatoire de la Côte d'Azur.

The dissipation rate $\gamma ({k}_{\perp },{k}_{\parallel })$ is evaluated from Equation (D.21) of Howes et al. (2006), obtained from the linearized gyrokinetic equations in the limit ${\delta }^{2}\tau \ll {\beta }_{i}=\tau {\beta }_{e}\ll 1$ (see definitions in the text). In a non-dimensional form, one has

$$\begin{eqnarray}&&\gamma =\sqrt{\displaystyle \frac{\pi }{2}}\displaystyle \frac{1}{{\beta }_{e}}\left(\displaystyle \frac{{{\rm{\Gamma }}}_{0}(\tau {k}_{\perp }^{2})}{{\tau }^{3/2}}\exp \left[-\displaystyle \frac{{\overline{\omega }}^{2}}{\tau {\beta }_{e}}\right]+\delta \right){k}_{\parallel }{k}_{\perp }^{2}.\end{eqnarray} \tag{ 6 }$$

Here ${\overline{\omega }}^{2}={k}_{\perp }^{2}(1+\tau -{{\rm{\Gamma }}}_{0}(\tau {k}_{\perp }^{2}))/(1-{{\rm{\Gamma }}}_{0}(\tau {k}_{\perp }^{2}))$, where Γ_n(x) = I_n(x)e^−x and I_n is the first type modified Bessel function of order n. While Equation (6) includes both ion and electron Landau damping, at small β_e, the primary contribution comes from electrons, so that Equation (63) of Howes et al. (2006) can also be used. Landau damping also affects the transfer times ${\tau }_{\mathrm{tr}}^{\pm }$ of both counter-propagating waves, due to the temperature homogenization process along the magnetic field lines on the correlation length scale ${k}_{\parallel }^{\pm -1}$. The associate timescale ${v}_{\mathrm{th}}{\widetilde{k}}_{\parallel }^{\pm }$, which explicitly arises in Landau fluid closures (Snyder et al. 1997; Sulem & Passot 2015), being proportional to the thermal velocity v_th of the particles, is very short for the electrons and cannot affect the dynamics. It is in contrast relevant in the case of the ions for which it is given in the present units by ${({\tau }_{H}^{\pm })}^{-1}=\mu \sqrt{2\tau }\,{\widetilde{k}}_{\parallel }^{\pm }$, where μ denotes a numerical constant of order unity. We are thus led to write

$$\begin{eqnarray}&&{\tau }_{\mathrm{tr}}^{\pm }={\tau }_{{NL}}^{\pm }\left(\displaystyle \frac{{\tau }_{{NL}}^{\pm }}{{\tau }_{w}^{\pm }}+\displaystyle \frac{{\tau }_{{NL}}^{\pm }}{{\tau }_{H}^{\pm }}\right),\end{eqnarray} \tag{ 7 }$$

which leads to

$$\begin{eqnarray}&&V=\displaystyle \frac{{v}_{\mathrm{ph}}^{2}}{{v}_{\mathrm{ph}}+\mu \sqrt{2\tau }}.\end{eqnarray} \tag{ 8 }$$

It is possible to relate the flux ratio η/ to the imbalance ${E}^{+}({k}_{0})/{E}^{-}({k}_{0})$ at the outer scale. This can be done by rewriting the spectra in the form

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{E({k}_{\perp })}{{k}_{\perp }}=\rho ({k}_{\perp })\cosh \phi ({k}_{\perp }),\\ \displaystyle \frac{{v}_{\mathrm{ph}}({k}_{\perp }){E}_{C}({k}_{\perp })}{{k}_{\perp }}=\rho ({k}_{\perp })\sinh \phi ({k}_{\perp }).\end{array}\end{eqnarray} \tag{ 9 }$$

This leads to

$$\begin{eqnarray}&&{E}^{\pm }({k}_{\perp })=\displaystyle \frac{1}{2}{k}_{\perp }\rho ({k}_{\perp }){e}^{\pm \phi ({k}_{\perp })}.\end{eqnarray} \tag{ 10 }$$

Two cases are to be distinguished, depending on the value of χ.

B.1. The Case χ = 1

From Equation (4) we see that, in the case χ = 1, solutions with constant fluxes η and obey

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{d}{{{dk}}_{\perp }}{\rho }^{2}({k}_{\perp })=-\displaystyle \frac{2\varepsilon \,{\widetilde{k}}_{\parallel }^{+}({k}_{\perp })}{{k}_{\perp }^{7}{v}_{\mathrm{ph}}({k}_{\perp })},\\ {\rho }^{2}({k}_{\perp })\displaystyle \frac{d}{{{dk}}_{\perp }}\phi ({k}_{\perp })=-\displaystyle \frac{\eta \,{\widetilde{k}}_{\parallel }^{+}({k}_{\perp })}{{k}_{\perp }^{7}},\end{array}\end{eqnarray} \tag{ 11 }$$

where ${\widetilde{k}}_{\parallel }^{+}={k}_{\perp }^{2}\sqrt{\rho /2}\exp (\phi /2)$.

• 1.  
  In the MHD regime where ${v}_{\mathrm{ph}}=\sqrt{2/{\beta }_{e}}$, it is easily shown that, defining a = ηv_ph/,

  $$\begin{eqnarray}&&\phi ={\phi }_{0}+a\mathrm{ln}\rho /{\rho }_{0}.\end{eqnarray} \tag{ 12 }$$

  When substituted into Equation ((11), left), this leads to $\rho ={\left(3-a){e}^{{\phi }_{0}/2}\epsilon \sqrt{{\beta }_{e}}/16\right)}^{2/3}{k}_{0}^{-8/3}{\left({k}_{\perp }/{k}_{0}\right)}^{-8/(3-a)},$ which, after some algebra, prescribes for the MHD regime

  $$\begin{eqnarray}&&\displaystyle \frac{{E}^{+}({k}_{0})}{{E}^{-}({k}_{0})}={\left(\displaystyle \frac{{k}_{d}}{{k}_{0}}\right)}^{\tfrac{16a}{3-a}},\quad \phi ={\phi }_{0}-\displaystyle \frac{8a}{3-a}\mathrm{ln}\displaystyle \frac{{k}_{\perp }}{{k}_{0}},\end{eqnarray} \tag{ 13 }$$

  where k_d denotes the pinning wavenumber. This prediction excellently matches the numerical results presented in Figure 2.
• 2.  
  In the far sub-ion range, v_ph ≈ αk_⊥, where $\alpha =2\sqrt{\tfrac{1+\tau }{{\beta }_{e}(2+(1+\tau ){\beta }_{e})}}$. The relationship between ρ and ϕ derived from Equation (11) reads in this case

  $$\begin{eqnarray}&&\displaystyle \frac{\alpha \eta }{\epsilon }{k}_{\perp }\displaystyle \frac{d}{{{dk}}_{\perp }}\mathrm{ln}\rho =\displaystyle \frac{d\phi }{{{dk}}_{\perp }}.\end{eqnarray} \tag{ 14 }$$

  If we assume $\eta \alpha /\epsilon \ll 1$, ϕ can be approximated by a constant ϕ₀ in Equation ((11), left) which is then solved as

  $$\begin{eqnarray}&&\rho ={\left(\displaystyle \frac{3\epsilon {e}^{{\phi }_{0}/2}\sqrt{\beta }}{10\alpha \sqrt{2}}\right)}^{2/3}{k}_{\perp }^{-10/3}.\end{eqnarray} \tag{ 15 }$$

  Equation (14) becomes

  $$\begin{eqnarray}&&\phi ={\phi }_{0}^{{\prime} }-\displaystyle \frac{10}{3}\displaystyle \frac{\alpha \eta }{\epsilon }{k}_{\perp },\end{eqnarray} \tag{ 16 }$$

  which determines the pinning wavenumber k_d where ϕ vanishes. From Equation (10), it is clear that for small αη/, ϕ will be nearly constant up to the vicinity of k_d where the spectra E^± approach each other exponentially. Therefore, we can obtain a dispersive imbalance relation, whose behavior is very different from the MHD case, namely,

  $$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{{E}^{+}({k}_{0})}{{E}^{-}({k}_{0})}=\exp \left(\displaystyle \frac{20\alpha \eta ({k}_{d}-{k}_{0})}{3\epsilon }\right),\\ {\phi }_{0}^{{\prime} }=\displaystyle \frac{10\alpha \eta {k}_{d}}{3\epsilon }\quad \mathrm{when}\quad {k}_{d}\gg {k}_{0}.\end{array}\end{eqnarray} \tag{ 17 }$$

  To apply this formula to the case where both MHD and dispersive ranges are present, one has to match ϕ at the transition wavenumber k_⊥ using both (13) and (16). But the contribution due to the second term in (13) is negligible and therefore we can simply extend (17) to the full range.

B.2. The MHD Regime with χ = 0

After some algebra, it is easily shown that the equations for ρ and ϕ read

$$\begin{eqnarray}&&\begin{array}{l}\rho \displaystyle \frac{d}{{{dk}}_{\perp }}\left({\rho }^{1/2}({k}_{\perp })\cosh (\displaystyle \frac{\phi ({k}_{\perp })}{2}\right)=-\displaystyle \frac{\epsilon }{{2}^{3/2}{k}_{\perp }^{5}{v}_{\mathrm{ph}}},\\ \rho \displaystyle \frac{d}{{{dk}}_{\perp }}\left({\rho }^{1/2}({k}_{\perp })\sinh (\displaystyle \frac{\phi ({k}_{\perp })}{2}\right)=-\displaystyle \frac{\eta }{{2}^{3/2}{k}_{\perp }^{5}}.\end{array}\end{eqnarray} \tag{ 18 }$$

From here, defining $\tanh \psi /2=\eta {v}_{\mathrm{ph}}/\epsilon$, we derive (A being a constant)

$$\begin{eqnarray}&&\sqrt{\rho }\sinh \displaystyle \frac{\phi -\psi }{2}=A,\quad \displaystyle \frac{d\phi }{{dk}}=\displaystyle \frac{C}{{k}_{\perp }^{5}}{\sinh }^{4}\displaystyle \frac{\phi -\psi }{2},\end{eqnarray} \tag{ 19 }$$

where $C=\tfrac{\epsilon /{v}_{\mathrm{ph}}}{{2}^{1/2}\,{A}^{3}\cosh \psi /2}$. Let us consider the case where energy and GCH are injected at k_⊥ = k_i. Imposing ρ = 0 as boundary condition at ${k}_{\perp }=+\infty$, we find that A = 0. As a result ϕ = ψ for k_⊥ > k_i. In this case, the imbalance becomes independent of ${k}_{\perp }$ in the form

$$\begin{eqnarray}&&\displaystyle \frac{{E}^{+}({k}_{\perp })}{{E}^{-}({k}_{\perp })}={\left(\displaystyle \frac{\epsilon +\eta {v}_{\mathrm{ph}}}{\epsilon -\eta {v}_{\mathrm{ph}}}\right)}^{2}.\end{eqnarray} \tag{ 20 }$$

We also get ${\rho }^{3/2}=\displaystyle \frac{3\epsilon \sqrt{{\beta }_{e}}}{16\,{k}_{\perp }^{4}}{\cosh }^{-1}\tfrac{\psi }{2}$. For ${k}_{\perp }\lt {k}_{i}$, where we can assume $\psi =0$ (in the absence of GCH injection at k₀), we can deduce from Equation ((19), right) (using that the constant C is large as a result of the small value of ρ at k_i) that ϕ also tends to zero, as it is observed in the simulations of the diffusion model (Figure 5, bottom). This explains the absence of propagation of the imbalance to large scales.

The transverse magnetic energy spectrum ${E}_{{B}_{\perp }}$, commonly measured in the SW, can be related, at least approximately, to the total energy spectrum. Writing ${\boldsymbol{k}}=({k}_{\perp },\theta ,{k}_{\parallel })$ in cylindrical coordinates, the total energy can be expressed as ${ \mathcal E }=\int E({k}_{\perp }){{dk}}_{\perp }$, with an energy spectrum given by $E({k}_{\perp })=(1/2)\int ({s}^{2}| {k}_{\perp }{L}_{e}{\widehat{A}}_{\parallel }{| }^{2}+| {k}_{\perp }{L}_{e}{\rm{\Lambda }}\widehat{\varphi }{| }^{2}){{dk}}_{\parallel }{k}_{\perp }d\theta$, equivalent to Equation (2.36) of Passot & Sulem (2019). Here, $\widehat{\varphi }({\boldsymbol{k}})$ and ${\widehat{A}}_{\parallel }({\boldsymbol{k}})$ refer to the Fourier transforms of the electrostatic and parallel magnetic potentials, respectively, ${L}_{e}={(1+2{\delta }^{2}{k}_{\perp }^{2}/{\beta }_{e})}^{1/2}$, and Λ defined in Equation (2.17) of Passot & Sulem 2019) is a function of k_⊥, which tends to 1 as ${k}_{\perp }\to 0$ and is proportional to $1/{k}_{\perp }$ in the sub-ion range. The first term in the integral rewrites ${L}_{e}^{2}{E}_{{B}_{\perp }}$, where ${E}_{{B}_{\perp }}=({s}^{2}/2)\int | {\widehat{B}}_{\perp }({\boldsymbol{k}}){| }^{2}{{dk}}_{\parallel }{k}_{\perp }d\theta$ corresponds to the magnetic energy spectrum, and the second one reduces in the MHD regime to the kinetic energy spectrum. Their difference, referred to as the residual energy spectrum, is observed to remain small if initially zero in direct numerical simulations of the parent gyrofluid. In spite of the nonlinear interactions, the solution can indeed be viewed as a superposition of eigenmodes of both polarizations (which satisfy ${s}^{2}| {k}_{\perp }{L}_{e}{\widehat{A}}_{\parallel }{| }^{2}=| {k}_{\perp }{L}_{e}{\rm{\Lambda }}\widehat{\varphi }{| }^{2}$), and we are thus led to write ${E}_{{B}_{\perp }}({k}_{\perp })\approx (1/2){L}_{e}^{-2}E({k}_{\perp })$.