A flux-balanced fluid model for collisional plasma edge turbulence: Model derivation and basic physical features

In toroidal magnetically confined plasmas, the level of cross field heat and particle transport in the edge region is in large part set by drift-wave turbulence driven by temperature and density gradients.¹⁷ This turbulence can itself generate zonal flows,⁸ which are known to mediate turbulent transport by shearing the turbulent eddies and absorbing some of the drift wave turbulence energy.^{3,6,8,11,18,28,32} Understanding this process is a critical step towards the goals of minimizing heat transport in magnetically confined plasmas and designing compact, economically viable fusion devices.

A wide variety of models have been used to study the drift wave turbulence–zonal flow dynamics in toroidal plasmas, with a varying degree of physics fidelity;^{1,11,14,15,21,28,29,32,33} interested readers may find a longer list of relevant references in Ref. 8. In the present work, we propose a new fluid model which continues the rich tradition of studies based on reduced fluid models^{1,7,13,16,21–24,30,31} initiated by the pioneering works of Hasegawa and Mima¹⁴ and Hasegawa and Wakatani.¹⁵ The approximations made to derive many of the reduced models, including ours, are often only satisfied for parameter regimes which do not correspond to the values measured in the edge of magnetic confinement fusion plasma experiments. We therefore do not expect numerical results obtained with reduced fluid models to be in quantitative agreement with measurements. However, because of their relative simplicity, reduced models have played an important role in characterizing the mechanisms responsible for the formation and dynamics of zonal flows and continue to be important tools to provide more insights into these important phenomena. In order to put our new model in context, we propose a very brief review of the main reduced models considered in the plasma physics literature.

The Hasegawa-Mima (HM)¹⁴ model is a one-field partial differential equation (PDE) which is the simplest known model containing the drift wave turbulence–zonal flow feedback loop mechanism. It was recognized in the 1990s^10,11 that in the context of magnetic confinement fusion, it is more physically relevant to consider a corrected form of the HM equations, now referred to as the modified Hasegawa-Mima (mHM),⁷ in which the magnetic surface averaged electron density response is subtracted from the original HM (oHM) equation.¹⁰ The mHM model has the desirable feature of being Galilean invariant under boosts in the poloidal direction⁷ and is known to lead to a stronger generation of zonal flows. This modification is also at the heart of the new model we propose.

A limitation of the HM models is that they do not contain a natural instability, so that turbulent forcing must be added externally in order to study the turbulence–zonal flow dynamics. In contrast, Hasegawa and Wakatani¹⁵ have shown that a generalized version of the HM models which includes electron ion friction in the parallel direction naturally contains a drift instability due to the finite plasma resistivity and thus drift wave turbulence induced transport. This is the original Hasegawa-Wakatani (oHW) model, a system of coupled PDEs for the ion vorticity (which is the Laplacian of the electrostatic potential) and the ion density. However, the oHW model has shortcomings too. Just like the oHM, it is not Galilean invariant. Furthermore, zonal flows are not generated in the oHW model,²¹ unless one makes the adiabaticity parameter wave number dependent.²⁴ Numata et al. have proposed a modified Hasegawa-Wakatani (mHW) model obtained by subtracting the zonal components from the resistive coupling term. The mHW model addresses the issues of the oHW mentioned above: its drift wave turbulence can lead to the strong generation of zonal flows, and it can be shown that the model has the desired Galilean invariance. Even if so, the mHW has the weakness that it does not converge to either the oHM model or the mHM model in the limit of zero resistivity and zero dissipation, as we will demonstrate numerically in this article.

We propose a generalized Hasegawa-Wakatani model with an improved treatment of the electron response along the magnetic field lines. Our model addresses the aforementioned limitation of the mHW model by solving for the mHM potential vorticity instead of the vorticity and reduces to the mHM model in the adiabatic, nondissipative limit, as desired, while maintaining the Galilean invariance of the mHW model. We call our model the flux-balanced Hasegawa-Wakatani model (bHW) to emphasize the fact that in this model, the surface-averaged electron density is constant in time when the electrons are adiabatic, as one would physically expect.^7,10,11 Our simple modification to the mHW model leads to major differences in the observed dynamics. Most significantly, the generation of zonal flows is enhanced and the turbulent fluctuations about the zonal mean state are increased. In our new Hasegawa-Wakatani model, we have also chosen to generalize the form of dissipative effects introduced in the oHW model to represent a wider variety of physical processes. As we acknowledged previously, we do not expect quantitative agreement between our fluid model and measurements from the edge of magnetic confinement fusion plasma experiments, but we hope to improve our understanding of the interplay between competing physical effects and of the implications of making simplifying fluid assumptions to model subtle kinetic effects. As an illustration, in this article we make the observation that the simplified linear Landau damping term of the oHW model¹⁵ which we also included in our model, acting mostly on the largest scale modes, may not always act as a pure dissipative effect and can effectively increase the variability in the flow fluctuations. To help the reader throughout the article, we provide a summary of the similarities and differences between the flux-balanced Hasegawa-Wakatani model and the modified Hasegawa-Wakatani model in Table I.

The structure of this paper is as follows. We first briefly review the classical Hasegawa-Mima and Hasegawa-Wakatani models in Sec. II, in order to provide a historical context for our modifications to the Hasegawa-Wakatani model. In Sec. III, we present our new Hasegawa-Wakatani model and derive its conserved quantities. In Secs. IV and V, we study the main differences between the mHW model and our model using numerical simulations; Sec. IV focuses on statistical considerations while Sec. V illustrates our conclusions and observations from Sec. IV with snapshots of results of direct numerical simulations and explores the role of the different dissipation terms. We summarize our work in Sec. VI and present the main linear stability properties of our new balanced model in  Appendix A.

In this section, we review the central features of the classical Hasegawa-Mima (HM) and Hasegawa-Wakatani (HW) models, as well as their modified versions,^10,21 known to excite more realistic and stronger drift-wave–zonal flow dynamics. This short review is meant to provide the motivation for our new model. The presentation is purposefully brief as clear, longer, and more detailed presentations are readily available in the literature.^7,8,17,21

For simplicity, all the models we will discuss in this article are considered for a shearless slab geometry,^2,7 in which the toroidal magnetic surfaces are imagined to be flattened into planes parallel to the y and z axes, as shown in Fig. 1, where (x, y, z) is the Cartesian coordinate system used to describe the geometry, x representing the radial distance, which can be viewed as a flux surface label, and y and z playing the roles of the poloidal and toroidal angles, respectively. The magnetic field is assumed to be solely in the z direction, $B=B0∇z$⁠, and B₀ is constant and uniform. The equilibrium density depends on the radial variable, $n0(x)$⁠, and we will treat the density profile within the framework of the standard “local approximation,”²⁷ in which $n0′(x)/n0(x)=constant$⁠.^21,31 The electron temperature is uniform throughout the plasma, and the ion temperature is assumed to be small compared to the electron temperature: $Ti/Te≪1$⁠. In all the models discussed in this article, the physical quantities will be assumed to be uniform in the z-direction, corresponding to the fact that in strongly magnetized plasmas, the dynamics is highly anisotropic, with much slower variations of the physical quantities along the magnetic field than across the magnetic field. The problem is therefore reduced to a two-dimensional problem in which all quantities only depend on x, y, and the time t, and the computational domain Ω is a rectangle whose sides have lengths L_x and L_y. Finally, the boundary conditions on the perturbed quantities (n, φ), which are the quantities the models solve for, are periodic in both x and y.³¹ 

The Hasegawa-Mima models, both original and modified, are equations for the ion vorticity obtained by combining the ion continuity equation and the ion momentum equation, together with the assumption of an adiabatic electron response to close the system of equations.^7,14 The only difference between the oHM and mHM models lies in the treatment of that adiabatic electron response. More than a decade after the introduction of the oHM model, it was indeed recognized that one should subtract the zonal mean of the electrostatic potential in the equation relating the adiabatic electron density and the electrostatic potential in order to prevent unphysical net radial transport of electrons.^7,9,10,24

The oHM and mHM equations can be formulated under the same framework by defining a switch parameter j with j = 0 for the oHM model and j = 1 for the mHM model, which accounts for the different treatment of the adiabatic electrons. The unified equation is^13,17

In Eq. (1), $J(φ,q)=∂xφ∂yq−∂yφ∂xq$ is the Jacobian associated with the advection term $vE·∇q$⁠, where $vE$ is the $E×B$ velocity, t is the time normalized to the ion cyclotron frequency $ωci=eB0/m$⁠, and x and y are normalized in terms of the hybrid ion thermal Larmor radius $ρs=ωci−1(Te/mi)1/2=miTe/eB0, δj0$ is the Kronecker delta, which is equal to 1 if j = 0 and 0 otherwise, $q=∇2φ−(φ̃+δj0φ¯)$ is the potential vorticity, with $φ=eϕ/Te$ the normalized electrostatic potential, where e is the charge of the electron and ϕ the electrostatic potential, and $κ=−d ln n0/dx$⁠. A bar over a quantity f represents the zonally averaged mean of that quantity, which only depends on x, and a tilde represents the fluctuation component of f, obtained by removing the mean from f

where L_y is the length of the domain in the y direction.

The mHM modification to the oHM model is simple, only appearing in the definition of the potential vorticity q. Yet it has important physical implications. First, unlike the oHM model, the mHM model is Galilean invariant under boosts in the poloidal direction,⁷ which is a desired property for our shearless slab geometry. Second, it can be shown, using linear theory, that in the absence of mean flow, the drift wave dispersion relation is identical in both models and given by $ω=kyκ/(1+k2)$⁠, where k_y is the y component of the wavevector k, and k the magnitude of k. However, in the presence of a constant and uniform background mean flow in the y direction, $vE=v¯ŷ$⁠, the dispersion relation is modified in different ways in the two models⁷ 

We see that at small scales, $k≫1$⁠, the dispersion relations of the oHM and mHM models agree and are given by $ω=ω*+kyv¯$⁠, where $ω*=κky/(1+k2)$ is the drift wave frequency. However, at larger scales, i.e., scales comparable to ρ_s corresponding to $k∼1$⁠, the dispersion relations differ. In the mHM model, the mean flow leads to a simple Doppler shift in the dispersion relation, as expected, but in the oHM model, the Doppler shift is reduced by the factor $k2/(1+k2)$⁠.

It is clear from the dispersion relations above that in the absence of mean flow or in the presence of a steady, uniform mean flow in the poloidal direction, the Hasegawa-Mima models do not have instabilities. This does not remain true if we assume the presence of a radially varying zonal mean flow. As we will show in an article currently in progress and more focused on the HM models, a drift instability can grow in these conditions and break up the zonal jet structure. Here too, the oHM and mHM models differ, and we find that the zonal jets are more likely to be broken up by this instability in the oHM model than in the mHM model. This provides a hint for the most important physical difference between the oHM and mHM models from the point of view of the present article as well as applications to magnetic confinement fusion, namely, the well-known result that the proper treatment of the electron adiabatic response in the mHM model leads to much stronger zonal jet structures when a random forcing term is added to the equation in order to mimic an instability leading to turbulent behavior.^5,7,11 This crucial observation provides the motivation for the mHW model and our new bHW model, which will be the focus of the remainder of this article.

The Hasegawa-Mima models capture essential features of the drift wave–zonal flow dynamics when a forcing term is included in the equation but does not include any internal drift instability to drive the drift wave–zonal flow feedback loop in the absence of forcing. The Hasegawa-Wakatani models^21,31 address this limitation by including electron-ion friction, which relaxes the slaving relation between the electron density and the electrostatic potential, leading to a drift wave instability.³⁰ Because the one-to-one correspondence between density and potential is lost, the HW models are two-field models for the ion vorticity and the ion density fluctuation. As for the HM models, one can make the distinction between the original HW model (oHW)³¹ and the modified HW model (mHW),²¹ which differ in their treatment of the parallel current. The mHW model accounts for the fact that zonal modes, with wavenumber k_y = 0, do not contribute to the parallel current, while the oHW does not. As for the HM models, we can write the two HW models in a unified form, with the switch parameter j such that j = 0 for the oHW model and j = 1 for the mHW model

where $ζ=Δφ$ is the ion vorticity and $n=n1/n0$ is the relative density fluctuation, with $N=n0+n1$ the total ion density, and where all the quantities have been normalized in the same way as for the HM models. The term $μΔζ$ is an approximate model for collisional ion viscosity perpendicular to the magnetic field, where $μ=μ̃/(ρs2ωci)$ with $μ̃=3Tiνii/(10miωci2)$ the kinematic ion viscosity coefficient, T_i the ion temperature, ν_ii the ion-ion collision frequency, and m_i the ion mass. We note that strictly speaking, Eqs. (2a) and (2b) are not the mHW model, in the sense that the mHW model has different dissipation terms than the oHW model. Specifically, in the mHW model, dissipation for the ion vorticity is written as $DΔ2ζ$ and dissipation for the ion density is written as $DΔ2n$⁠, with D the same unspecified constant in both equations. We do not dwell further on this distinction, since these terms are not given any physical justification in Ref. 21, and seem to only be included in order to guarantee numerical stability.

The parameter $α≡Tekz2/n0e2ηωci$⁠, where η is the resistivity parallel to the field lines, is often referred to as the adiabaticity parameter.^4,21,24 In the collisionless limit, $α→∞$⁠, the electrons have an adiabatic response along the field lines, whereas for small α electron-ion friction decouples the density and electrostatic potential. In order to make the mathematical connection with the HM models in the asymptotic collisionless limit $α→∞$ and $μ→0$⁠, we can derive an equation for the potential vorticity $q=∇2φ−n$ valid for both models, by assuming that the oHW and mHW models both only include the viscous dissipation $μΔζ$ as the dissipation mechanism (these dissipative terms will not matter in the end since we will take the limit $μ→0$⁠). Subtracting Eq. (2b) from Eq. (2a), we find

We see that in the collisionless limit $μ→0$⁠, Eq. (3) has the same form as the equivalent equation (1) for the HM models, with the HM potential vorticities $∇2φ−(φ̃+δj0φ¯)$ replaced with the HW potential vorticity $q=∇2φ−n$⁠. In the collisionless limit, $α→∞$⁠, Eqs. (2a) and (2b) give n = φ in the oHW model, so the oHW potential vorticity becomes $q=∇2φ−φ$⁠, and the oHW model coincides with the oHM model as desired.

However, in the mHW model, $α→∞$ leads to $ñ=φ̃$⁠, leaving $n¯$ unspecified. The potential vorticity converges to $q→∇2φ−φ̃−n¯$ and the mHW equations (2a) and (2b) do in general not coincide with the mHM model in the zero resistivity limit. Physically, the mHW model does not guarantee the absence of a net radial transport of electrons when the electrons are adiabatic. Throughout this article, we will return to this important conclusion, which is confirmed by our numerical simulations, as it is a key difference between the mHW model and the new bHW model we introduce in Sec. III. Even if the mHW model does not converge to the mHM model in the collisionless limit, it has been shown that the improved treatment of the parallel current as compared to the oHW model leads to stronger zonal mean structures than in the oHW model,²¹ just as the mHM model does compared to the oHM model.

 Appendix A contains a detailed discussion of the linear stability of resistive drift waves in the HW models. At this point, it suffices to say that for $α≠0, κ≠0$⁠, and μ = 0, all the modes are unstable and the growth rates are smaller for small-scale modes, $k≫1$⁠.

We are now ready to introduce the new bHW model, which addresses a major shortcoming of the mHW model, namely the fact that it does not converge to the mHM model in the collisionless limit. In this new model, we also consider generalized dissipation effects, which allow us to study separately the roles of the different dissipation terms, which are often included for the mere sake of numerical stability, and whose simple forms are rarely justified physically.

Let us start with the dissipation terms. Starting from the mHW model in Eq. (2), we relax the ad hoc assumption that the viscosity coefficients are identical for the ion vorticity and the density fluctuation. Furthermore, we include the model Landau damping term considered in the original Hasegawa-Wakatani article when the drift wave frequency $ω*$ is comparable to $kzvTi$⁠, where $vTi$ is the ion thermal velocity.³¹ The equations become

Here, $μΔζ$ is the same ion collisional viscosity term as in (2a), and $DΔn$ can be viewed as a dissipation term due to the collisional diffusion of the electrons perpendicular to the magnetic field.⁴ Finally, $Cω*φ$⁠, with $C=Ti/(ωciTe)$⁠, is the simple model for ion Landau damping introduced in Ref. 31. Observe that Hasegawa and Wakatani suggest a cutoff for this term, which they set to zero when $ω*$ is large. For simplicity, we do not impose such a cutoff, noticing that the absence of cutoff will only be noticeable when k_x is small and k_y is moderate, since $ω*$ is small for large k, as well as for small k_y. Furthermore, even if $ω*$ can be larger than $kzvTi$ when k_x is small and k_y moderate, the ion viscosity term, $μΔζ$⁠, is still likely to dominate in these circumstances, since it contains a term proportional to $ky4$⁠.

We would like to stress the fact that to the best of our knowledge, the forms of our dissipation terms $μΔζ, DΔn$⁠, and $Cφ$ cannot be derived rigorously from kinetic theory for magnetic confinement plasmas. However, our choice to consider distinct coefficients μ and D and to introduce different forms of dissipation allows us to analyze the relative importance of these effects and determine which effects would need to be modeled more accurately in a higher fidelity model. For example, as we will discuss in more detail in Sec. V in the light of our numerical results, and in our linear stability analysis in  Appendix A, we find that Landau damping as modeled in Eq. (4) has a stabilizing effect on the largest scales and at the same time increases the variability of the small-scale fluctuations.

We now turn to the key modification to the mHW model we make in the bHW model. As we did in Sec. II, if we subtract Eq. (4b) from Eq. (4a), the terms with the adiabaticity parameter cancel, and we obtain an equation for the potential vorticity in the mHW model $qm≡ζ−n=Δφ−n$ which reads

As discussed in Sec. II, the issue with the equation above is that in the limit $μ,D,C→0$ and $α→∞$ it violates the balanced electron response on the magnetic surfaces in the mHM sense, due to the inclusion of a non-zero zonal mean state $n¯$ in the potential vorticity. We propose a simple modification to the mHW model to address this limitation, namely, to replace the original potential vorticity $qm$ with the corrected form $qb=ζ−ñ$⁠. The resulting flux-balanced equation for the new potential vorticity $qb$ does not depend on the zonal mean density $n¯$ explicitly. The immediate and desired implication of this modification is that in the collisionless limit, $α→∞$ gives the slaving relation $ñ→φ̃$ from Eq. (4b), and the potential vorticity $qb$ converges to the mHM potential vorticity $q→∇2φ−φ̃$ unlike the mHW potential vorticity. As a result, in the adiabatic limit and in the absence of the dissipative terms, Eq. (5) when expressed for $qb$ is identical to Eq. (1) of the mHM model.

We now have all the elements to introduce the flux-balanced Hasegawa-Wakatani (bHW) equations for the balanced potential vorticity $qb=∇2φ−ñ$ and the relative particle density $n=n¯+ñ$ as

The potential vorticity equation (6a) has a form analogous to the vorticity equation in the mHM model, Eq. (1), but unlike the mHM model, it contains a resistive drift instability through the coupling with the density fluctuation $ñ$⁠. The equation for the relative particle density fluctuation (6b) is the same as in the mHW equation (2b), except for the dissipation term $DΔn$⁠.

In this section and Sec. III C, we present some elementary properties of our model. We begin with a direct comparison with the mHW model. A good starting point is to write the equations for the balanced potential vorticity $qb=Δφ−ñ$ according to the mHW model and to the bHW model side by side

with $ũ=−∂φ̃/∂y$ the velocity fluctuation along the x direction and

with $f̃=ñ$ as in the equation above, or $f̃=q̃$ as we will have below. The two differences in the mHW model are (1) the zonal density gradient feedback term $(∂xn¯)∂φ̃/∂y$⁠, which modifies the original background density gradient profile $−κ=d(ln n0)/dx$⁠, and (2) the additional eddy flux feedback, $∂x(ũñ¯)$⁠, due to the transport of the particle density along the x direction. Both terms are zonal mean density feedbacks which can act on the largest scales (to change the zonal jet structure) and can also modify the structures in the small-scale fluctuations (to change the particle flux). It is important to observe that although the zonal mean density $n¯$ does not appear in Eq. (7b), $qb$ depends on $n¯$ in the bHW model through the dependence of $ñ$ on $n¯$⁠, which we will write explicitly shortly. For the same reason, the vorticity $ζ=∇2φ$ depends on $n¯$ in the mHW model.

Let us now turn to the equations for the zonal mean quantities $q¯b(x)$ and $n¯(x)$⁠. If, in order to highlight similarities between the mHW and bHW models, we modify the ad hoc dissipation terms in the mHW model in such a way that the dissipation in the density equation is $DΔn$ in both models, we obtain the following equations for the mean balanced potential vorticity $q¯b=∂2φ¯/∂x2$⁠:

and an equation for the zonal mean density $n¯$ which is identical in both models

Equations (8) and (9) highlight the feedback of the fluctuations on the zonal mean states, through the nonlinear coupling terms $∂x(ũq̃b¯)$ and $∂x(ũñ¯)$⁠. In the mHW model, the density fluctuation feedback $∂x(ũñ¯)$ is also present in the equation for the mean balanced potential vorticity (8a), just like it is in Eq. (7a). Considering the steady-state version of Eq. (9), we see that a strong particle flux $ũñ¯$ can lead to a large zonal mean density gradient $∂n¯/∂x$ which can then become a significant contribution in the term $(∂xn¯−κ)∂φ̃/∂y$ in Eq. (7a) and effectively alter the background density gradient in the mHW model.

A summary of the similarities and differences between the bHW and mHW models is given in  Appendix B, including the equations for the total energy E and enstrophy W which we will discuss in Sec. III C. Once more, in order to draw the appropriate parallels between the two models, we modified the form of the ad hoc dissipation originally proposed in the mHW model²¹ in such a way as to make the dissipation term agree in the density equation of both models.

For the analysis of the fundamental properties of a new model and the verification of numerical codes, it is always helpful to identify dynamical invariants. We define the total energy E and potential enstrophy W²¹ for the bHW model as

where the integrals are over the computational domain Ω, which is a periodic box with sides L_x and L_y, and where we have separated $E¯,W¯$⁠, the energy and enstrophy of the zonal mean state with $qb¯=∂2φ¯/∂x2$⁠, from $Ẽ,W̃$⁠, the energy and enstrophy of the fluctuations about the zonal mean state, with $qb̃=Δφ̃−ñ$⁠. Note that the enstrophy W in the bHW model is defined in terms of the balanced potential vorticity $qb=∇2φ−ñ$⁠. It is easy to verify that the nonlinear terms $J(φ,qb)$ and $J(φ,n)$ in (6) conserve both the energy and the balanced enstrophy, and we obtain the following dynamical equations for the total energy and potential enstrophy:

where $v¯=∂φ¯∂x$ and D_E and D_W come from the dissipation terms in the bHW equations

Comparing the energy and enstrophy equations (11) with the analogous equations for E and W in the mHW model, as given in  Appendix B, we notice the additional term $∫Ωv¯ũñ¯dV$ in the energy equation for the bHW model. This contribution to the total energy from the fluctuations originates from the eddy diffusivity term $(ũñ¯)x$ in the equation for the vorticity in the bHW model, which is absent in the mHW model. This additional term, in which we recognize the advection by the mean velocity $v¯$⁠, represents the zonal flow transport of the particle flux $ũñ¯$⁠.

As in the mHW model, the resistive coupling of $ñ$ and $φ̃$ through the adiabaticity $α(φ̃−ñ)$ gives a negative-definite term in the energy equation and acts as an energy sink. That term does not enter in the evolution equation for the enstrophy W. The dissipation for the total energy D_E is always positive definite. This is not true for the dissipation D_W of the potential enstrophy, whose terms may not all always be positive. Specifically, the first integral on the right hand side contains interactions between the fluctuating potential and density $Δφ̃Δñ,ñφ̃$ due to ion Landau damping and the fact that the damping coefficients μ and D are not equal; the last term in the second integral, $(D−μ)|∇ñ|2$⁠, is always negative if $D<μ$ and can thus act as a source of potential enstrophy instead of a sink. In Sec. V B, we illustrate these different dissipation effects on the flow field and energy spectra through direct numerical simulations.

To close this section, we propose a more mathematically rigorous proof of the convergence of the bHW model to the mHM model in the adiabatic limit $α→∞$⁠. We have already observed that formally, the limit $α→∞$ implies that $φ̃=ñ$⁠, meaning that the bHW potential vorticity is equal to the mHM potential vorticity and that the models are unified in this limit. Here, we rely on the energy and enstrophy equations (11) to find a more precise description of the relation between the bHW model and the mHM model in the adiabatic limit. For simplicity, we assume that $μ=D$ and C = 0, so that the dissipation terms in the energy/enstrophy equations take the simple form $DW=D∫|∇qb|2$ and $DE=D∫(|Δφ|2+|∇n|2)$⁠. We then consider states in statistical equilibrium, i.e., such that the time derivatives for the statistical expectations $⟨E⟩$ and $⟨W⟩$ vanish. Here, $⟨·⟩$ represents an ensemble average, which could in principle be estimated with a Monte-Carlo approach, using a large number of numerical simulations. When statistical equilibrium has been reached, the enstrophy equation in (11) gives us the following estimate for the total particle flux:

where $⟨·⟩eq$ can be viewed as the time average along the stationary trajectory of the functional f by making the usual ergodic hypothesis. This relation shows that there always exists a statistically positive particle flux towards the boundary of the domain and that the total particle flux $Γeq$ depends on the ratio $D/κ$⁠.

Turning to the equilibrium statistical equation for the total energy, we obtain the following estimate from the balance between the flux and dissipation terms:

To derive the inequality above, we used the identity $κΓeq=⟨DW⟩eq$ and assumed that the zonal mean flow can be bounded as $V=maxΩ|1+κ−1v¯|$⁠. The statistical expectation on the right hand side of the inequality is positive and finite in the equilibrium state. Therefore, the expectation $⟨(ñ−φ̃)2⟩eq$ vanishes as $α→∞$ since the right-hand side of the inequality goes to zero. This demonstrates that in the adiabatic limit, the fluctuations $ñ$ and $φ̃$ approach the Boltzmann relation $ñ=φ̃$ in the mean square sense under expectation. Note however that there may still exist a non-zero zonal structure in the density $n¯(x)$ in this limit, in addition to the fluctuation $ñ=φ̃$⁠.

For the remainder of this article, we turn to direct numerical simulations to highlight characteristic features of the bHW model and compare them with the mHW model. We solve the equations on a doubly periodic square domain of size L along each side, so that the smallest wavenumber is $2π/L$⁠, which is also the spacing $Δk$ between any two wavenumbers. We write the quantities $(φ,n,ζ)$ we solve for as the following Fourier series:

where the spatial variables $x=(x,y)∈[−L/2,L/2]×[−L/2,L/2]$ and the corresponding wavenumbers are

where N is the number of modes we keep in our simulations. For all the results shown in this article, we used L = 40 and N = 256. We know from our linear stability analysis (see  Appendix A) that the larger the domain size L, the more unstable modes we will have inside the circle of strong instability. We solve the equations using a pseudo-spectral approach, in which the nonlinear terms are calculated in real space instead of Fourier space. To stabilize the truncated numerical system, hyperviscosities $νΔsqb$ and $νΔsn$ are added in the vorticity and density equations, respectively, with $ν=7×10−21$ and order s = 8. We rely on a fourth-order explicit-implicit Runge-Kutta scheme for time stepping, where the stiff hyperviscosity operator $νΔ8$ is the only term treated implicitly. In agreement with the difference in the formulations of the mHW and bHW models, we numerically solve for the unknowns $(ζ,n)$ when we consider the mHW model and solve for the balanced potential vorticity $qb=∇2φ−ñ$ and n when we consider the bHW model.

We keep the background density gradient fixed at $κ=0.5$ and will vary the value of α in the range $[0.01,10]$⁠. In this section, the dissipation coefficients μ and D are set to $μ=5×10−4,D=5×10−4$⁠, and the ion Landau damping term is turned off by setting C = 0. In contrast, in Sec. V we will consider different values for the parameters $μ,D$⁠, and C to analyze the role of dissipation effects more closely. The parameters we used for all our numerical simulations in this article are summarized in Table II.

Previous work has shown that the randomness of the turbulent flow depends sensitively on the parameters κ and α.²¹ This can be readily understood from the linear stability analysis we present in  Appendix A, since $κ/α$ determines the size of the region of instability in k-space. Hence, if κ is held fixed, a smaller value of α, corresponding to higher plasma resistivity, leads to a larger region of instability and a more energetic and turbulent vorticity field. When the value of α is increased, the vorticity field is regularized into dominant anisotropic zonal jets. This is the phenomenon we now study, comparing the bHW and mHW models.

Because of the turbulent nature of the flows, it is more appropriate to adopt a statistical viewpoint, and instead of looking at physical quantities at a given instant in time, we take time-averages of the solution once the stationary state is reached. We will focus on the total kinetic energy and the relative enstrophy for both models, given by

The relative enstrophy gives more weight to the small scale modes, while the kinetic energy puts more emphasis on the large scale structures. We decompose $φ$ and ζ into their time averages $⟨φ⟩eq, ⟨ζ⟩eq$ and their time fluctuating parts $φ′, ζ′$ according to $φ=⟨φ⟩eq+φ′, ζ=⟨ζ⟩eq+ζ′$⁠, and compute the energy and enstrophy in the time-averaged mean

and the time-averaged kinetic energy and enstrophy of the variances

The results are shown in Fig. 2 for simulations of the bHW equations and the mHW equations, with α varying in the interval $[10−2,10]$⁠. Since one of the main reasons to study the reduced HM and HW fluid models is to better understand the drift wave–zonal flow interplay, we separately plot the kinetic energy and enstrophy contained in zonal modes, i.e., with k_y = 0, and the total kinetic energy and enstrophy obtained by summing over all the resolved modes. We first focus on the plots for the variances, in the left-hand column, to observe that in both models, the kinetic energy and enstrophy of the variances decrease as α increases. This is a clear signature of the transition from a strongly turbulent regime to a zonal flow dominated regime, as we expected. For α small, the flow is strongly fluctuating with large and small scale modes both containing considerable amount of energy. As α increases, stronger zonal jets develop and the variances in kinetic energy and enstrophy decrease. In the regime $α≫1$⁠, the system tends to the HM equations which do not have an internal instability, and the variances reach minimum values.

Still looking at this first column, we highlight differences between the bHW and mHW models. The bHW model always contains a large variance for the kinetic energy in comparison to the mHW model, which has weaker variability. Furthermore, much of the variance in the bHW model is in zonal modes, unlike the mHW model. Our improved treatment of the electron dynamics parallel to the field lines enhances the zonal flow variability as the system transitions to the turbulence dominated regime when α decreases. The plot for the variance of the relative enstrophy, which is a physical quantity which gives more weight to smaller scale structures, helps to explain the different dynamics in the two models. Indeed, we notice that the mHW model has a much larger variance in enstrophy as $α<0.1$ than the bHW model. This is due to the generation of many small and strong vortices, as Fig. 5 in Sec. V illustrates explicitly with a contour plot of the vorticity field. Even if so, the zonal mode variance in the relative enstrophy remains larger in the bHW model. This is a central point of this article: even for very small values of α, the variability of the turbulent fluctuations is in zonal modes in the bHW model, while the mHW equations describe fully homogeneous turbulence in the small α limit.

We now consider the kinetic energy and enstrophy in the mean, in the second column of Fig. 2. The kinetic energy and enstrophy in zonal modes almost overlap with the total kinetic energy and enstrophy, proving that the mean state is always dominantly in the zonal direction. We can distinguish three different regimes depending on the value of α. Starting at around $α∼0.1$⁠, we notice a significant jump in the value of the kinetic energy and the enstrophy in the mean as compared to situations with $α<0.1$⁠, indicating the presence of stronger zonal jets. A close comparison between the plots for the kinetic energy and enstrophy shows a large kinetic energy and a moderate enstrophy for $α∈(0.1,0.5)$⁠, which is a signature of the presence of one or very few zonal jets. For $α>0.5$⁠, the kinetic energy is a bit lower and fairly constant as a function of α, whereas the enstrophy increases. This indicates the presence of multiple zonal jets. These conclusions are confirmed by the snapshots of the vorticity field for each of the α regimes discussed here, which are shown in Fig. 5 in Sec. V.

Besides the kinetic energy and enstrophy in the statistical equilibrium state discussed in Sec. IV A 1, it is also important to investigate the dynamical evolution of the statistical quantities, in particular to better understand the nonlinear energy exchange processes in the model and the model sensitivity to various perturbations. While in many situations we are mostly interested in the statistics of the two lowest order moments, that is, the statistical mean and variance, the higher-order feedbacks to these low order moments play a key role in the dynamics and should be accurately treated, as has been highlighted and illustrated in Refs. 19, 20, and 25 for general systems with quadratic nonlinearity. In this section, we highlight these mechanisms for the bHW and mHW models.

We consider here the statistical average $⟨k̃2|φ̂k|2⟩$ of the kinetic energy in each spectral mode $φ̂k$ of the electrostatic potential. To derive the dynamical equations for $⟨k̃2|φ̂k|2⟩$⁠, we first project the equations for the vorticity in the mHW and bHW models onto each individual spectral mode, to obtain equations for each $φ̂k$⁠. We then multiply by $φ̂k*$ on both sides of the equations and take the statistical ensemble average (see Ref. 20 for details) to obtain the desired equations for $⟨k̃2|φ̂k|2⟩$⁠. For the bHW model, the equations are

Here, c.c. stands for the complex conjugate completion. Importantly, the nonlinear interaction term in the vorticity equation produces the higher-order feedback $QF,kb=−⟨φ̂k*J(φ,Δφ+n¯)k⟩+c.c.$ for the non-zonal modes $ky≠0$⁠; and $QF,kb=−⟨φ̂k*J(φ,Δφ−ñ)k⟩+c.c.$ for the zonal modes k_y = 0. For the mHW model, the equations are

with the higher-order feedback $QF,km=−⟨φ̂k*J(φ,Δφ)k⟩+c.c.$ The removal of the zonal mean density $n¯$ in the term involving the adiabaticity parameter α in both models, as well as in the definition of the potential vorticity in the bHW model, is responsible for the different form the equations take for the zonal modes. We finally note the presence of the cross-covariance term $⟨φ̂kn̂k*⟩$ in both models, representing the interactions between the potential and density.

As emphasized in Refs. 19 and 20, the third-order moments on the left hand side of Eqs. (12) and (13) play the critical role of mediating the growing unstable modes and driving the system to the final equilibrium. As we know from the linear stability analysis, the linear operators on the right hand sides of the statistical equations contain positive eigenvalues, corresponding to the linearly unstable directions. If the third order moments on the left hand sides are not included, the internal instability will lead to fast growth in kinetic energy among the unstable modes and fast decay in the other over-damped modes. When the third-order moments are kept, they transfer the growing energy from the unstable subspace to the stable subspace, where the energy eventually gets dissipated.

We thus see the importance of a rigorous analysis of the contributions from the third-order moments. Unfortunately, it is usually expensive to compute these moments from direct numerical simulations since it requires the inclusion of all the triad modes across the entire spectrum.²⁰ On the other hand, the third-order moment feedbacks in statistical equilibrium are much easier to compute. Indeed, in statistical equilibrium, the time derivatives on the left hand sides of Eqs. (12) and (13) vanish, so the third-order moments can be expressed in terms of the first and second order moments appearing on the right-hand sides. Specifically, for the bHW model, we can write

The equilibrium moments, $⟨|φ̂k|2⟩eq,⟨φ̂k*n̂k⟩eq$⁠, can be readily calculated from the equilibrium mean and variance by averaging along the stationary time trajectory. The equilibrium third-order moments for the mHW model can be computed in an analogous way.

To numerically investigate the differences and similarities between the mHW and bHW models from the point of view of the statistical moments, we choose a regime with small α, namely, $α=0.01$⁠, which corresponds to a physical situation with strong mean flow–drift wave interactions. Figure 3 presents the equilibrium statistics of the bHW and mHW models for both the second-order variances and the third-order moment feedbacks in the spectral domain. We first look at the first column, which shows contours of the logarithm of the second-order variance of each energy mode $k2|φ̂k|2$⁠. The stronger variability for the zonal modes (k_y = 0) in the bHW model is clearly visible. This is consistent with the observation we made in Sec. IV A 1 that even in the small α regime, the bHW model generates energetic zonal flows which are responsible for most of the flow variability. In contrast, we note that the contour plot for the mHW model is radially symmetric in the Fourier domain, indicating the absence of anisotropy in the flow statistics.

We next turn to the third-order moment feedbacks shown in the second and third columns of Fig. 3. In the middle column, the logarithm of the magnitude of the third-order moments is plotted in order to amplify the small values, and we use dashed lines to represent positive values and solid lines for negative values. In the last column, we plot the contours of the third-order moments themselves (i.e., not their logarithm). The red colors correspond to positive values and the blue colors to negative values. We first observe that both the bHW and the mHW models show strong non-zero contributions from the third-order moments, particularly at the largest scales in the system. Comparing the results of the linear stability analysis we present in Fig. 8 of  Appendix A with these figures, we see that the linearly unstable modes are subject to negative third-order moment feedbacks which act as an effective damping mechanism against the growth in energy. At the same time, the linearly stable modes are subject to positive third-order moment feedbacks. This is characteristic of a nonlinear transfer of energy from the unstable modes, whose energy grows in time, to the stable modes, which then dissipate the energy. Note that the third-order moment feedback conserves energy as a whole, so it cannot act as a source or sink of energy, and instead redistributes energy at different scales. Finally, one observes much stronger third-order moment feedbacks for the zonal modes in the bHW model than in the mHW model, which we would expect given our previous results: once again, even in the regime $α≪1$⁠, zonal modes play a central role in the dynamics in the bHW model, while they are essentially absent in the mHW model. These results also highlight the need for accurate approximations of the third-order moment feedback terms when applying model reduction strategies to investigate the statistical properties of a system and for uncertainty quantification. In the present case, the differences between the bHW and mHW models, localized in the vicinity of k_y = 0, lead to major differences in the observed dynamics.

Both the bHW and the mHW models converge to an HM-like model without instability in the limit of large α, as we confirmed numerically in Sec. IV A. We showed in Sec. III C that the bHW model converges to the mHM equation (1) (with an additional dissipation term) in the limit $α→∞$⁠, with the balanced potential vorticity converging to its limiting form $qb=∇2φ−ñ→∇2φ−φ̃$⁠. We also explained in Sec. III A that the mHW equations do not guarantee exact convergence to the mHM model due to the incomplete treatment of the zonal mean density $n¯$⁠. In this section, we verify our reasoning with numerical simulations. Specifically, we consider a case in which α is relatively large, α = 5, which is a regime where the linear growth rates are small, so the energy in the small-scale modes is mostly dissipated, and a dominant zonal mode with multiple jets sets in, as is clearly visible in the first column of Fig. 4, which shows contour plots of the ion vorticity ζ at the final time of the simulation, in the top left corner for the bHW model, and in the bottom left corner for the mHW model. To obtain a more quantitative idea of the selective energy decay in each spectral mode, we plot the kinetic energy $k̃2|φ̂k|2$ of each mode at different simulation times in the second column of Fig. 4, with the bHW results at the top and the mHW results at the bottom as before. The energy decay in the small scale modes can be observed for both models at early times, and for both models, the decay saturates at around t = 2500. At the final time, only a dominant zonal mode with three zonal jets remains in the flow field, as can be seen in the contour plots in the first column of Fig. 4. However, the decay rate in the mHW is much slower, with a considerable amount of kinetic energy maintained in the small scales, while the small scale kinetic energy in the bHW model decays quickly to negligible amounts after a short transient state. This fundamental difference between the mHW and bHW models is clearly apparent in the contour plots of the vorticity field in the first column. In the bHW simulation in the upper left corner, the vorticity field is largely dominated by the large scale zonal jet structure, with very limited smaller scale features. On the other hand, the mHW simulation shows many persistent small scale vortices in the vorticity field.

To further verify the convergence to the mHM model, we numerically solve the mHM equations with dissipation and random noise perturbations added at small scales, starting from an initial state given by the final state of the bHW simulation with parameters $κ=0.5,α=5$⁠. The results are shown in the right column of Fig. 4. Since the mHM model does not contain any instability, the fluctuations get dissipated while the zonal jet structure is maintained in time. The top right figure shows the final and converged zonally averaged velocity $v¯=∂xφ¯$ for both the bHW and mHM model simulations. The zonal mean flow clearly converges to the same saturated limit. We have verified that this is true whatever the initial state is chosen for the mHM simulation. Although not shown in the figure, the contour plot for the mHM ion vorticity has strong similarities with the bHW figure, with the same dominant zonal jets, and very few small scale structures, unlike the mHW results.

This concludes our numerical proof of the convergence of the bHW model to the mHM model in the adiabatic limit and of the incomplete convergence of the mHW model in that same limit.

In this section, we illustrate the key results discussed in Sec. IV with plots of the ion vorticity and mean flow obtained with direct numerical simulations for different values of the adiabaticity parameter α. We also study the separate roles of the dissipation terms $μΔζ, DΔn$⁠, and $Cω*φ$ by considering different values for the parameters μ, D, and C than what we had in Sec. IV.

In Sec. IV, we highlighted the critical role of the adiabaticity parameter α in determining the turbulent flow regime in the HW models, with a transition to strong zonal jets as α increases. This transition is clearly visible in Fig. 5, which shows contours of the ion vorticity ζ for the bHW model in the first row and for the mHW model in the second row. These figures are obtained for the parameters $μ=D=5× 10−4$ and C = 0, and for the final time of the simulation, which is much later than the time when statistical equilibrium has been reached. The three contour plots in each row correspond to three different values of the adiabaticity parameter α. From left to right, the values are $α=0.01,α=0.1$⁠, and $α=0.5$⁠. We see that for $α=0.01$⁠, many small-scale vortices are visible on top of the background mean vorticity. As the parameter value increases to $α=0.1$⁠, a stronger single jet pattern begins to form, while still coexisting with many smaller scale fluctuations. Finally, for $α=0.5$⁠, the vorticity field becomes fully dominated by several strong zonal jets. These figures offer an explicit illustration of the statistical results given in Fig. 2.

Figure 5 also confirms another key observation regarding differences between the bHW and mHW models. The bHW model maintains jet structures for a wide range of values of the adiabaticity parameter α, throughout the transition from a regime with dominant zonal jets (large α) to the strong drift wave turbulence regime (small α). For small values of α, the jets are more turbulent and shift in time, but the anisotropic zonal dynamics persists. This is in stark contrast to the mHW model, which loses the jets as $α→0$⁠, as a regime with fully homogeneous turbulence and strong vortices sets in. Because of the persistence of the zonal jets, the particle flux is always smaller in the bHW model than in the mHW model for small values of α. Physically, the unbalanced density flux in the mHW model is responsible for the highly turbulent vorticity and strong particle transport in the limit of small α.

In the large α regime, we highlight another critical difference between the mHW and the bHW models, namely, the fact that in the bHW model, both the zonal mean flow and the fluctuations have a much larger variability than in the mHW model. To demonstrate this, in Fig. 6 we plot the time series of the zonal mean velocity field $v¯=∂φ¯/∂x$ for $α=0.5$⁠, which corresponds to a strong zonal jet regime. The jets generated in the bHW model have large amplitude variations in time, whereas the zonal velocity in the mHW model is mostly steady in time with an almost constant jet amplitude.

Thus far, we have always turned off ion Landau damping by setting C = 0 and set the dissipation coefficients D and μ to be equal, with $μ=D=5×10−4$⁠. In this last section, we take a closer look at the right-hand side of Eq. (6a) by dissociating μ and D and considering finite values for C. Specifically, we study the consequences of increasing the values of μ to $μ=2×10−3$ while keeping the other parameters to their original values and of increasing D to $D=2×10−3$ while keeping the other parameters to their original values. In addition, we test a small value for C, namely, C = 0.01. Since Landau damping mostly acts on the largest scales, we focus in this section on the regime corresponding to $α=0.5,κ=0.5$⁠, where there exist strong but fluctuating large-scale zonal mean modes.

Figure 7 provides a summary of our main results. The top left panel shows time series of the total energy in the system for the different values of μ, D, and C considered in this section. The blue, red, and yellow curves correspond to cases in which there is no Landau damping, and we see that increasing D by a factor of 4 to $D=2×10−3$ only decreases the total energy by a small amount. As expected from our linear stability analysis in  Appendix A, increasing μ by a factor of 4 to $μ=2×10−3$ has a slightly stronger effect, decreasing the total energy to a noticeably lower value, but it does not significantly change the nature of the dynamics. In contrast, ion Landau damping greatly changes the energy time series, even for small values of C, as shown with the purple curve. In this regime with strong zonal jets, Landau damping increases the total energy and leads to much larger fluctuations. This is because our naive model for Landau damping damps the large scales but increases the energy at the intermediate and small scales, as we now explain with a spectral viewpoint.

To better understand the results found for the time series of the total energy, we show the energy spectra of the variance $⟨Etot′2⟩eq$ and the statistical mean $⟨Etot⟩eq$ in the top center and top right panels of Fig. 7, respectively. These figures confirm that increasing D has a small impact on small scale modes and does not significantly modify the largest zonal mean modes of the flow either. The impact of a larger value of μ is more noticeable but approximately uniform across all scales. On the other hand, setting the Landau damping coefficient C to a finite value has a strong effect in changing the zonal flow profile $v¯$⁠. Landau damping effectively removes the energy of the zonal modes at the largest scales and drives the original zonal flow with 3 jets to a configuration with 5 dominant jets in the flow. At the same time, it increases the variance at all scales and particularly strongly at intermediate scales. We illustrate this effect by showing the time-series of the zonal mean flow $v¯=∂xφ¯$ with Landau damping in the second row of Fig. 7. Comparing this figure with the analogous Fig. 6 for the same regime in the absence of Landau damping, we see that the jets are more turbulent when Landau damping is included, meandering in time with stronger fluctuations. Moreover, instead of a persistent 3-jet structure as in the C = 0 case, the flow shows a dominant 5-jet structure, with sudden merging events leading to a 4-jet flow, before the 5 jets reemerge after a short period of time.

In this article, we presented a new two-field fluid model for the study of the drift wave–zonal flow dynamics in magnetically confined plasmas. Our model is inspired by the rich tradition of Hasegawa-Wakatani (HW) models,^21,24,31 which are the simplest known two-field models for magnetized plasmas which contain a drift-wave instability and can generate zonal flows through nonlinear drift wave interactions. The main novelty of our balanced HW (bHW) model as compared to previous HW models is its improved treatment of the electron dynamics parallel to the magnetic field lines, which guarantees a balanced electron density response^7,10 in the limit in which electrons are adiabatic. To achieve this, the bHW model does not solve for the ion vorticity, but instead for a well constructed potential vorticity defined by $qb=∇2φ−ñ$⁠, which we call the balanced potential vorticity and does crucially not include the zonal mean density $n¯$⁠.

Our new bHW model inherits the desirable features of former HW models as well as those of the Hasegawa-Mima model with a modified adiabatic electron response, which we have called mHM.⁷ Like other HW models, the bHW model contains a drift wave instability, and like the mHW model, the bHW model has the desired Galilean invariance along a magnetic flux surface; on the other hand, like the mHM model, the improved electron response in the bHW model leads to enhanced zonal flows and strong fluctuations. Mathematically, the bHW model has the satisfying property of converging to the mHM model in the proper adiabatic limit, which is not the case of previous HW models.

We relied on direct numerical simulations to investigate the main features of the bHW model, with numerical simulations of the recent modified HW model (mHW)²¹ as a point of comparison. We paid close attention to the transition from a strongly drift wave turbulent regime to a much more organized regime with strong zonal structures, which is associated with a decrease in the plasma resistivity, allowing electrons to become more and more adiabatic. We observed that the bHW model has stronger fluctuations than the mHW model in the turbulent regime and throughout the transition to the regime with strong zonal jets. Remarkably, however, the bHW model maintains mainly zonal dynamics even at large plasma resistivity, while in contrast the drift wave turbulent regime in the mHW model is dominated by small scale dynamics. In other words, in the bHW model, the zonal flows are more robust, in the sense that they are observed for any value of the plasma resistivity, but they have a larger variability than in the mHW model in the regime in which they are seen in both models. We also used direct numerical simulations to confirm the convergence to the mHM model in the zero resistivity limit and to study the role of the dissipation terms introduced in our model. We noticed that the simple model for Landau damping introduced by Wakatani and Hasegawa³¹ leads to more turbulent dynamics, with zonal structures which are fluctuating more, and more energy and variability at intermediate and small scales. The significant effect of the Landau damping term on the bHW dynamics provides yet again a clear justification for the development and implementation of more accurate fluid approximations for this effect, as was first done in the pioneering work of Hammett, Perkins, and Dorland.^10,12

This article is meant to be the first in a series of articles in which we study the interplay between drift waves and zonal flows using the new bHW model presented here. In a second article,²⁶ more focused on the analysis of detailed numerical simulations, we investigate the effect of the shape and size of the computational domain on the properties of the drift wave turbulence and zonal structures. In a subsequent article, we will treat the bHW model and the drift wave–zonal flow dynamics it describes as a test bed to explore the applicability of statistical model reduction strategies^20,25 for the efficient description of turbulent magnetized plasmas. The latter is the topic of ongoing research, with progress to be reported in the near future.

A.J.C. would like to thank Jeffrey B. Parker for the suggestions he made during an insightful conversation which led to this work. A.J.M. was partially supported by the Office of Naval Research through MURI N00014-16-1-2161 and DARPA through W911NF-15-1-0636. D.Q. was supported as a postdoctoral fellow on the second grant. A.J.C. was partially supported by the U.S. Department of Energy, Office of Science, Fusion Energy Sciences under Award Nos. DE-FG02-86ER53223 and DE-SC0012398.

In this appendix, we present a detailed linear stability analysis for the bHW model given by Eq. (6). For this analysis, we assume a zero background mean flow, so the linearized equations for the perturbations $(φ1,q1b,n1)$ to the background state are given by

The dispersion relation for the bHW model is found by considering a single mode solution for the perturbations, i.e.,

where we assume $ky≠0$⁠. Substituting the above single mode solution in the linearized equations, we find the following equations for $φ̂1$ and $n̂1$⁠: