Linking the two-field dynamics of plasma edge turbulence with the one-field balanced model through systematic unstable forcing at low resistivity

The self-generation and amplification of zonal flows from the interplay with turbulent drift waves are key constituents of particular interest in the investigation of magnetically confined plasmas.^1–5 This is because the zonal flow–drift wave interaction mechanism is thought to have a critical role in the observed level of heat and particle transport perpendicular to the projected magnetic surfaces.^6–10 Numerical simulations are the most direct way for the study of the crucial roles in the complex flow structures and the development of theories.^11,12 However, direct simulations of the entire nonlinear plasma equations remain very expensive computationally. Thus, the selection of important elementary physical processes in the zonal flow and drift wave dynamics using reduced models is necessary.

The use of simplified models based on assumptions for plasma regimes has advantages in improving our understanding of the key features in drift wave–zonal flow interaction, where the most relevant physical mechanism can be identified. Among them, the Hasegawa-Mima (HM) models [also known as the Charney-Hasegawa-Mima (CHM) model in geophysics]^13,24 should be one of the simplest formulations including the most essential physics modeling the adiabatic electrons at zero resistivity.^14,15 It offers a qualitative characterization for the typical observations in realistic plasma flows with desirable analytical tractability. A modified HM model is proposed later in Refs. 2 and 15 with the attractive features of restoring the zonal flows and Galilean invariance under boosts in the poloidal direction and improves the CHM model by introducing a balanced particle response on magnetic surfaces. The major drawback of such models is the lack of instability to maintain the turbulent solutions, and so usually, external forcing is required in order to excite the turbulent drift wave dynamics.

The Hasegawa-Wakatani (HW) models, on the other hand, by introducing the coupled two-field evolution of the electrostatic potential and the density perturbation, contain more complete physics by inherently including the drift wave instability through finite plasma resistivity.^16,17 A modified HW (MHW) model is further proposed in Ref. 18 to reinforce the strong zonal flows by properly removing the zonal contributions to the model for the parallel gradient of the parallel current density. It has the desirable features of Galilean invariance and zonal jets. Formal analysis suggests that the MHW model solution approaches similar zonal structures as observed in the HM model in the corresponding adiabatic limit.² However, with a more precise characterization of the full turbulent field, the exact convergence is not guaranteed in the MHW model with persistent zonal transport. See Fig. 4 in Ref. 17 and Fig. 3 in Ref. 12 for moderate values of the resistivity, as well as Fig. 14 in the present paper for a low resistivity case. A balanced flux treatment for the electron parallel responses is introduced in Refs. 12 and 17 to further improve the model performance. The balanced flux constraint comes from the zero net density fluctuation response averaged on the magnetic surfaces for adiabatic electrons. The modified HM and HW models with the balanced flux have been shown to be more physically satisfying in creating many desirable features, such as the Galilean invariance in the poloidal direction and the stronger persistent zonal jets^5,12,17 with simultaneously reduced particle flux.

Here, we propose a simple HM formulation which systematically incorporates the desirable features observed from the two-field flux-balanced HW model. The strategy is to introduce a forcing effect to the modified HM model through a precise analytical derivation of the drift wave instability in the low resistivity regime. The analysis is carried out through the expansion of the analytic instability solution of the HW model as the adiabaticity approaches infinity (that is, the zero resistivity limit). The two-field model states are shown to converge to the unstable branch solution in the leading-order approximation equivalent to the one-field HM model potential vorticity. A “systematic unstable forced modified Hasegawa-Mima” (“SUF”–MHM) model is then proposed. The leading order instability correction to the HM equation is treated as an explicit external forcing derived from the limit drift instability in the HW model. The unstable forcing is added to the modified HM model to excite fluctuating drift waves in the same fashion as in the HW model. Zonal jets and similar turbulent dynamics can be generated in the next stage from this one-field SUF-MHM model by selective decay effects and secondary instability mechanisms.^4,5 The forcing added to the model is based on the observation that the regime of strongest linear instability is directly linked with the regime with secondary instability of drift waves for nonlinear energy transfer.

The SUF-MHM model is used to mimic the exact drift wave–zonal flow energy feedback loop.^12,19,20 Galilean invariance is automatically satisfied through the model construction. The modeling procedure generates zonal jets through the following flow developments: (i) excitation of drift waves from the consistent primary instability due to resistive particle motion in the HW model; (ii) generation of zonal flows through the nonlinear interaction with drift waves; (iii) effective quench of strong radial particle transport and fluctuations with the formation of zonal jet barrier; and (iv) recovery of the saturated steady state statistics in model variables. Numerical tests are carried out in the low resistivity regime. Similar dynamical evolution is observed from simulations of the one-field SUF-MHM model and the balanced Hasegawa-Wakatani (BHW) model (see Fig. 1). Further statistical agreement is achieved in a quantitative way for flows in regimes with a moderate density gradient. By comparing the similarity and difference in the solutions between the two model simulation results, a better understanding about the energy mechanism can be gained in the roles of nonlinearity and instability generated in the plasma turbulence.

In the structure of this paper, background and basic ideas in the flux-balanced models are introduced first in Sec. II. Then, the growth rate and corresponding eigenstates are derived in a systematic fashion for the linearized two-field HW model in Sec. III. The SUF-MHM model is constructed based on the limit form of growth rate. Detailed numerical simulations of the models are provided in Sec. IV. The importance of the balanced flux correction is further emphasized in Sec. V. A summary discussion is given in Sec. VI.

The Hasegawa-Wakatani models describe the coupled drift wave–zonal flow interactions with a system of two fields.^16,18 The system is defined on a shearless two-dimensional slab geometry, where the magnetic field is embedded. In convention, the x-coordinate corresponds to the radial direction, and the y-coordinate represents the poloidal direction. The “flux-balanced Hasegawa-Wakatani” (BHW) model improves the original HW models by using the potential vorticity $q=∇2φ−ñ$ with balanced electron response on the magnetic surfaces^12,17 and the density fluctuation n in the following coupled partial differential equations:

where $φ$ is the electrostatic potential, n is the density fluctuation from background density n₀(x), and $u≡∇⊥φ=(−∂yφ,∂xφ)$ is the velocity field. The constant background density gradient $κ=−∇ ln n0$ is defined by the exponential decay profile near the boundary n₀(x). D acts on the two states with the Laplace operator as a homogeneous damping effect.

Drift wave instability is generated from the resistive electron parallel motion through the finite adiabaticity, $α∝1η$⁠, treated as a constant and reciprocal to the resistivity η. The physical quantities $φ$ and n are decomposed into zonal mean states $φ¯,n¯$ and their fluctuations about the mean $φ̃,ñ$ so that

A modified Hasegawa-Wakatani (MHW) model was proposed in Ref. 18 by removing the zonal components in the resistive coupling α(φ – n) to become $α(φ̃−ñ)$ on the right hand side of Eq. (1b). It is shown that this modification is essential for the generation of zonal jets.

As a further correction in the BHW model, the poloidally averaged density $n¯$ along the y-direction is removed from the potential vorticity $q=∇2φ−ñ$⁠. The BHW model offers a more physically relevant formulation with several desirable properties.^12,17 Most importantly, it is shown from rigorous proof and numerical confirmation¹⁷ that in the “low resistivity” limit, η → 0 so that α → ∞, the BHW model converges to the one-field equation as desired

which is called the “modified Hasegawa-Mima” (MHM) model. The modification as compared to the standard HM model is by removing the zonal state $φ¯$ in the definition of potential vorticity q above.

By comparing the one-field MHM formulation (2) with the more complicated two-field BHW model (1), most of the desirable physical features are already modeled in the MHM framework with less computational requirement. However, there exists no internal instability in the MHM model; thus, the solution will simply decay from its initial state if no other external forcing is introduced. It is worthwhile to introduce an equivalent forcing operator simulating the drift wave instability in the two-field system. One direct idea is to introduce the equivalent forcing on each spectral potential vorticity mode with an unstable growth characterizing the drift wave instability in a precise way.

The systematic unstable forced modified Hasegawa-Mima (SUF-MHM) model is then introduced incorporating the basic idea illustrated above

with $F$ being a specific spatial nonlocal operator acting on the nonzonal potential vorticity mode $q̂k$ to model the drift instability [see Eq. (13) for the explicit formulation]. For a precise modeling of the instability effect, it is important to make sure that the added forcing $F$ is free of any adjustable parameters. One of the main tasks of this paper is to offer a systematic derivation of the suitable unstable forcing form, consistent with the two-field BHW model as it approaches the adiabatic limit α → ∞. This is achieved by considering the linearized dynamics of the two-field model and expanding the leading order contribution in the low resistivity limit from the exact linear analysis of the growth rate. It is found that the higher order contributions decay at a much faster rate as α → ∞. We will carry out the detailed derivation in Sec. III and numerical tests of the SUF-MHM model in Sec. IV. Statistically consistent results with similar transient behavior are generated using this simplified model compared with the BHW model for low resistivity regimes.

One important property to point out first for the forcing operator is the maintenance of “Galilean invariance” in the balanced formulation (3) under velocity boost V in the poloidal direction with the transformation

In fact, with the balanced particle response by removing the zonal mean state, the potential vorticity $q′=∇2φ′−φ̃′=q$ is unaltered under the above change in variables. By requiring that the forcing operator $F$ is applied only on the fluctuation modes with k_y ≠ 0, Galilean invariance is guaranteed in this HM model framework with balanced flux (3). This is also automatically satisfied for all the balanced equations.

As a further comment, the same framework with an unstable forcing can also be applied to the Charney-Hasegawa-Mima (CHM) model using the original potential vorticity q = ∇²φ – φ without the balanced response correction in the zonal mean. Galilean invariance is then not valid due to the zonal mean state contribution $φ¯$ in the unbalanced potential vorticity q. We will show in Sec. V A that the unstable forcing leads to only homogeneous turbulence without zonal jets in the unstably forced CHM model (as depicted in the selective decay in Refs. 4 and 21).

First, we display the typical features generated from the one-field SUF-MHM model (3) in comparison with the two-field BHW dynamics (1) by running direct numerical simulations in the same parameter regime. Here, we point out the most representative observations in the typical test case using parameters κ = 0.5 and α = 5 (which as shown below still contains considerable amount of turbulence in the flow field away from the adiabatic α = ∞ limit). Both the simulations start from random initial data with small amplitudes. The detailed numerical setup and complete numerical results with more turbulent features will be discussed thoroughly in Sec. IV.

We display the self-organization of the turbulent states in plasma flow evolution from direct simulation results of both the SUF-MHM and BHW models. The first two rows of Fig. 1 show several snapshots of the ion vorticity ζ = ∇²φ at several typical time instants before the steady state is reached. The one-field SUF-MHM model successfully captures the key physical features at every stage during the evolution of the model. In the starting time with a small amplitude random initial state, nonzonal drift waves are first excited from the drift instability (fluctuating drift wave state, t = 1000); then, the energy in fluctuations begins to transfer to the zonal modes through nonlinear interactions where a competition between the zonal modes and nonzonal drift waves can be observed (coexistence of zonal jets and strong fluctuations, t = 1500); finally, the dominant zonal jets get formed with the nonzonal fluctuations mostly dissipated (zonal jets dominant regime, t > 2000). The SUF-MHM model generates the same representative structures at every dynamical stage in the evolution compared to the two-field BHW model.

Second, for comparing the statistical consistency at equilibrium, the last row of Fig. 1 compares the equilibrium energy spectra achieved from SUF-MHM and BHW model simulations. The equilibrium statistical spectrum is computed by averaging the energy $k2|φ̂k|2$ in each mode along a long time series after the statistical steady state is reached. For a more detailed calibration of the statistics, both the radially averaged spectrum including the nonzonal fluctuating modes (by taking the summation of all the radial modes with same absolute wavenumber k) and zonal spectrum (with only zonal modes k_y = 0) are compared for the two models. Good agreement in both spectra between the two models is achieved for all the scales. This shows that the quantitative skill of the SUF-MHM model is correctly generating the model statistics in each scale but with the much simpler MHM model structure.

In this section, we construct the Hasegawa-Mima model for α ≫ 1 with an unstable forcing systematically derived from the linear instability analysis of the two-field Hasegawa-Wakatani model. This linear instability generated from the resistive drift waves leads to the excitation of nonzonal fluctuations from the initial state with little energy. Specifically, we are interested in the limit performance of the model as the adiabaticity α approaches the low resistivity regime. In this way, the two-field HW model can be decoupled into a single field HM model with an additional instability forcing representing unstable drift waves in the leading order approximation. In addition, the excitation of the linear unstable nonzonal fluctuating modes is closely related to the nonlinear energy transfer mechanism to zonal modes.

Drift wave instability is due to the nonadiabatic resistive electron motion. We consider purely fluctuating states $(φ̃,ñ)$ with zero background mean flow profile, $φ¯≡0,n¯≡0$⁠, in order to focus on the linear instability from resistive drift waves. The HW models (1) yield the linearized system if we drop the nonlinear terms from the original equations

We formulate the system based on the fluctuating vorticity $ζ̃=∇2φ̃$ and the fluctuating density $ñ$⁠, and the nonlinear coupling terms, $∇⊥φ·∇ζ$ and $∇⊥φ·∇n$⁠, are neglected in the linearized formulation. The above system can be viewed as the dominant dynamics in the starting transient state when the state values are small and nonlinear interactions have not taken over to add a major effect.

We assume that the linear solutions of fluctuating states with nonzonal modes k_y ≠ 0 are taking the following single-mode forms (the subscript k for the single mode variables is neglected for simplicity)

where ω ≡ ω(k) is the wave frequency for the corresponding wavenumber. The dispersion relation can be found by plugging in the above single mode solution. The system decouples into independent subsystems for each single wavenumber since we do not consider the nonlinear terms in the linearized system. The linearized coefficients then form the 2 × 2 system for each wavenumber

Nontrivial solution $(φ̂,n̂)≠0$ of the linearized HW model (4) yields the equation for the wave dispersion relation

with $k2=kx2+ky2$ being the wavenumber square and the dispersion relation $ω*≡ω*(k;κ)$ for the one-field HM model drift waves

The background density gradient κ only contributes to the HM dispersion relation ω_* without drift wave instability. The particle resistivity parameter α adds instability into the system. The homogeneous damping operator with strength D acts as a stabilizing effect of the system, acting strongest on the small-scale fluctuating modes.

In general, the quadratic equation (5) gives two complex roots, $ω±=ωr±+iωi±$⁠, where ω_r and ω_i are the corresponding real and imaginary components of the eigenvalues. At the resulting wave frequency, $exp (−iωt)=exp (−iωrt)eωit$⁠, the real part ω_r represents the wave dispersion and the imaginary part ω_i characterizes the growth rate (for the positive value) or the damping rate (for the negative value) due to the linear instability effect. The two eigenvalues ω^± correspond to the two branches of the eigenmodes representing the characteristic directions for unstable growth or stable damping, that is,

In the above equation, $(φ̂+,n̂+)$ represents the unstable eigen-direction where the energy grows exponentially in time as $exp (ωi+t)$ and $(φ̂−,n̂−)$ represents the stable eigen-direction where the energy gets quickly dissipated in the exponential rate $exp (−|ωi−|t)$⁠. Next, we compute exact solutions for the eigenvalues (5) and eigenmodes (6) based on the values of the adiabaticity parameter α.

In the above analysis, we first provide the general formulas for the linear instability in drift waves with combined effect of dissipations. Specifically, D is fixed at a small value D = 5 × 10⁻⁴ in Sec. IV as used also in Refs. 12 and 17. In the absence of the dissipation effect D = 0, it is more straightforward to compute the dispersion relation from (5) for nondissipative drift waves

Immediately, we can observe that in the high resistivity limit α = 0, the drift waves become nondispersive without any instability as ω = 0; and in the limit with no background density gradient κ = 0 so that ω_* = 0, the wave dispersion frequency is purely imaginary as $ω=−iα(1+k−2)$⁠, that is, always stable with a negative growth rate.

Now, we calculate the explicit solutions for the eigenvalues and eigenvectors for the nondissipative case of (5) and (6). By directly solving the quadratic equation with nonzero parameters α ≠ 0 and κ ≠ 0, the eigenvalues of the system for wavenumber k can be written explicitly as

where for simplicity in the representation, we introduce the parameter $Γ≡Γ(k;κα)$⁠, only dependent on the ratio $κα$ of the two model parameters

The two branches of the eigenvalues are determined by the parameter θ^± = Arg(–1 + 4γi). In the form of the solution (7), for fixed wavenumber k, it first depends linearly on the adiabaticity parameter α in the outside coefficient. While inside the square bracket, the instability feature is determined by the operator defined by Γ, only dependent on the ratio $κα$⁠.

Still, it is useful to get the explicit expressions for all the components of the solutions. A simple calculation gives that

The two branches of the solutions can be discovered by the signs of the sine and cosine functions depending on the signs of k_y. Putting all the expressions together, we derive the entirely explicit formulas for the two branches of the eigenvalues as

with $γ=καkyk2(1+k2)2$⁠. We introduce the additional parameters for dispersion and growth/damping, $ϖ,ς±$⁠, which are only related to the ratio $κα$⁠. For fixed κ, as α → ∞, we have γ → 0, and then, the corresponding parameters approach the limit ϖ → 0 and $ς+→0,ς−→2$⁠.

Correspondingly using the formulas for the two eigenvectors (6), the two branches of unstable and stable eigenmodes can be written explicitly as

Note that the eigen-directions only depend on the ratio $κα$ from the parameter γ. They retain the self-similarity in the convergence process based on only the value of $κα$⁠. We observe that as $κα→0$⁠, the unstable branch approaches the limit $n̂+=φ̂+$ which is exactly the HM model state for the potential vorticity $q̂+=−k2φ̂+−φ̂+$⁠; while the stable branch reaches $n̂−=−k2φ̂−$ which has no contribution to the potential vorticity $q̂−=−k2φ̂−−n̂−≡0$⁠.

Above, we derived the explicit formulas for the eigenvalues and eigenvectors for any parameter values κ and α in linear instability analysis. Here, we are interested in the two-field HW model performance as it approaches the adiabatic limit with low resistivity η → 0 or α → ∞. As the system approaches the low resistivity regime α ≫ 1, the above formula for the dispersion relation ω can be approximated in the leading order [by using the expansion $(1+x)1/2=1+12x−18x3+O(x4)$ twice] near the value γ ∼ 0 using the following expansions:

Putting the above leading order expansions back into the expression for the eigenvalues (8), we find approximation of the eigenvalues up to the order O(γ³) as

Notice that in the parameter $γ=καkyk2(1+k2)2$⁠, the wavenumber dependence part is always bounded, $kyk2(1+k2)2<1$⁠. The parameter ratio $κα$ determines the bound of γ.

By substituting the explicit form of the parameter γ into the above expansions, we find the explicit forms for the wave frequency ω_r together with the growth rate σ⁺ and the damping rate σ⁻ in the leading orders

From the above expression (10) for the wave frequency ω_r, the leading order term just gives the dispersion relation for the HM drift wave ω_*. Then, the second order offers further correction for this dispersion relation from the two-field model. The next order term decays fast according to the parameter ratio $κ5α4$⁠. The unstable growth σ⁺ shows the leading order growth rate along the unstable direction. The next order depends on the parameter ratio $κ4α3$ decaying also in a much faster rate. As α → ∞, the instability in the system vanishes as σ⁺ → 0 converging to the HM model limit. On the other branch with the damping rate σ⁻, the leading order gives an isotropic damping only dependent on the absolute wavenumber value k. As α → ∞, this term becomes especially strong and dominant driving the energy along this direction to zero rapidly.

Then, we consider the corresponding eigenvectors in this leading order expansions. The direct calculation from the previous formulas (9) gives the unstable and stable eigenmodes in the leading order expansion for the corresponding basis with growth and decay

Consistent with our previous intuitive approximation, the leading order expansions of the eigenstates give the exact HM potential vorticity $q̂+=−k2φ̂+−φ̂+$ along the unstable direction, and the stable branch makes no contribution to the potential vorticity $q̂−=0$ in the leading order as α → ∞. The same HM drift wave frequency ω_* is recovered in the leading order term of ω_r. With the more detailed next order expansions in (11), we can also observe convergence of two branches of the two-field BHW model to the one-field HM limit as α grows large. In fact, it has a self-similarity in the leading order $O(κα)$⁠, which predicts the same leading order statistics for models with the same parameter ratio $κα$⁠. The convergence in the eigen-directions shows invariant performance with constant parameter ratio $κα$⁠. We summarize the major results of the above analysis in Table I.

We check the convergence of the expansion formulas with numerical computations. In Fig. 2, we compute the normalized growth rate αω_i directly from the linearized two-field model (4) in comparison with the leading order approximation formula (10). The background density gradient is fixed at κ = 0.5 or κ = 1 and the values of α varies for different resistivity. No damping effect D = 0 is added in this test. The maximum growth rate always takes place along the k_y axis with k_x = 0. The growth decays fast to small values in the smaller scale modes after the peak. As the system goes toward the adiabatic limit α → ∞, the leading-order approximation in Eq. (10) offers an accurate approximation for the linear growth rates. In comparing the two different parameter cases κ = 1 and κ = 0.5, the same structure is produced for the linear growth consistent with the theoretical formula. The larger density gradient case κ = 1 generates stronger instability growth, and it is four times larger than the other case κ = 0.5 proportional to the coefficient κ² according to the first-order expansion.

Now, we consider the inclusion of homogeneous damping effect DΔ on both the vorticity q and density n equations as in (1). By introducing the single-mode states in the previous form with time dependent coefficients

the original linearized equation (4) gives the following form with the dissipation effect added

If we choose the dispersion relation ω exactly as the solution of the nondissipative system (7), the coefficients are only subject to the damping effect with all the other terms canceled. The damping effect can be easily eliminated by introducing the damping contribution on the original expansion formula as

Clearly in the unstable branch with growth rate σ⁺ in (10), the damping operator acts as the balancing effect of the linear instability induced in the first order of σ⁺. For each wavenumber, the instability is withheld by the damping when

The marginal stability boundary is determined by the contour when equality in the above relation is reached

From the instability (12) with the damping effect, the linearly unstable regime with positive growth rates is constrained on the “localized” nonzonal fluctuating modes in the largest scales. This corresponds to the first excited drift wave base state. Then, the secondary instability^5,21 from the nonlinear interaction transfers energy from the transient fluctuating modes to the zonal directions, where dominant zonal jets are created.

For illustrating the combined effect with dissipation and instability, Fig. 3 shows the linear growth rate including a weak damping D = 5 × 10⁻⁴ as well as the two typical values of the density gradient κ = 0.5 and κ = 1. The leading-order growth σ⁺–Dk² is compared with the two-field model exact solution. With the dissipation effect, the unstable growth vanishes at large wavenumbers and instability with positive growth is localized inside the metastable boundary defined in (12).

The linear growth rates with damping in the spectral space can be viewed equivalently as the injection of energy in the starting transient state through drift wave instability. In the next stage, the nonlinear effect will take over and balance the energy growth especially among nonzonal fluctuating modes. The mostly unstable nonzonal modes (with k_x = 0) get an equivalent damping effect from the third-order interactions, while energy is transferred to the zonal modes through the nonlinear interactions. The entire energy mechanism will be described through the statistical characterization in Sec. IV C.

The previous discussion offered an explicit form for the growth rate calculated from the expansion of the two-field HW model in the low resistivity regime. The leading order eigenstate gives the Hasegawa-Mima model at the zero resistivity limit, and the leading order growth rate introduces drift wave instability to the state variables. Now, the idea is to start with the leading order MHM model without instability and find the proper equivalent external forcing form to add to the MHM model that gives the same energy mechanism as the drift wave instability introduced from the HW system.

From the idea in Eq. (3), the internal drift wave instability can be directly modeled as an unstable forcing exerted on each nonzonal potential vorticity mode in the form $γ̂kq̂k$⁠. At the low resistivity limit, α → ∞ (or more generally κ/α → 0), linear stability analysis in Eq. (11) states that the unstable branch solution converges to the MHM model state $n̂+=φ̂+$⁠; while the potential vorticity contribution in the stable branch just vanishes $q̂−≡0$ in the leading order term. The higher order corrections are shown to decay in a much faster rate compared to the leading term. This observation implies that as the system goes near the low resistivity regime, the potential vorticity q in the two-field BHW model gradually becomes aligned with the unstable direction with the assigned unstable growth rate σ⁺. The unstable forcing strength then can be mimicked by the unstable branch of the explicit solution derived from the leading order instability (10), while it is reasonable to neglect the contribution from the stable branch at this limit due to its fast time decay rate σ⁻.

Therefore, the systematic unstable forced modified HM (SUF-MHM) equation proposed in Eq. (3) can be rewritten explicitly in the following form:

The same flux-balanced potential vorticity $q=∇2φ−φ̃$ is introduced with special treatment for the zonal modes as in the MHM model. The explicit antidamping effect from the leading order expansion of the growth rate injects energy to each nonzonal spectral mode to simulate the contribution from the unstable resistive drift waves

The growth $γ̂k$ from the unstable forcing is used to mimic the excitation of drift waves in the linearized two-field HW dynamics in the low resistivity regime, where the two-field HW model gets gradually aligned with the HM potential vorticity state. Importantly, we make it independent of any adjustable parameters so that no further tuning of the model is necessary. Zonal modes with k_y = 0 are not forced, and the nonlinear interaction terms act to transfer the increased fluctuating energy to zonal modes and downscale to get dissipated. The energy instability mechanism to excite drift wave turbulence in the starting transient state is modeled in a quantitative way. The systematic unstable forcing for the HM model achieves the same statistical behavior in key model quantities at the statistical steady state as shown in Fig. 1 and Sec. IV.

An important quantity we would like to model from the one-field SUF-MHM model is the “zonal particle flux” $Γ¯=ñũ¯=−ñ∂yφ̃¯$ that quantifies the total zonal transport of particle density. In the HM electrostatic potential function φ₀, the zonal particle flux vanishes at the adiabatic limit φ₀ = n₀ from the unstable branch of the model

The dominant particle flux should come from the next order expansion term in Eq. (11). Using the expansion for the unstable eigenmode, we have the approximation for the density fluctuation in each single wavenumber mode

The above relation in fact represents the balance between the dispersive drift waves and the resistive particle feedback in a short time scale. Taking the contribution from the second component into account, the particle flux can be approximated by the first order correction for each spectral mode as

with $n2=nx2+ny2$ and the summation taken over all permitted indexes m + n = k. Specifically, we can compute the total particle flux by taking the summation about the zero mode m + n = 0

Consistent with the two-field model case, the total particle flux approximation (15) shows always a positive particle transport toward the boundary direction from the unstable branch solution. As the system approaches the adiabatic limit α → ∞, the particle flux becomes weaker and finally will vanish at the zero resistivity.

Direct numerical simulations are carried out for the one-field SUF-MHM model to check the performance of the explicit leading order model (13) in comparison with the two-field BHW model (1). The equations are solved on a doubly periodic square domain of size L along each side. The lowest wavenumber becomes Δk = 2π/L. The variables of interest (φ, n, ζ) get the following spectral representations under Galerkin projection on Fourier modes

with the spatial variables $x=(x,y)∈[−L/2,L/2]×[−L/2,L/2]$ and the corresponding spectral wavenumbers $k=2πLn,n=(nx,ny)∈ℤ2$⁠. The standard pseudospectral code^12,17 to discretize the model with a 3/2-rule for de-aliasing the nonlinear terms is applied on the square domain with the size of L = 40 and resolution of N = 256. A fourth-order explicit-implicit Runge-Kutta scheme is used to integrate the time steps. Only a weak dissipation effect D = 5 × 10⁻⁴ is introduced in both electrostatic potential and the density (the same values used in Refs. 12 and 17). In addition, to stabilize the truncated numerical system, a hyperviscosity νΔ^2sq is added with the strength of ν = 7 × 10⁻²¹ and order of s = 4. The stiff hyperviscosity operator is integrated in an implicit scheme, while the explicit scheme is applied for all the other terms. The time integration step is kept small Δt = 5 × 10⁻³ throughout all tests. Two values for the background density gradient κ = 0.5 and κ = 1 are considered. The value for adiabaticity is fixed α = 5 which can still generate quite turbulent features in the plasma flows for both tested density gradients κ.

In all the numerical tests, the initial state is set from a random Gaussian field with a small amplitude. In this way, we are able to investigate the roles of linear instability and nonlinear interaction for the self-generation of dominant zonal structures starting from the small homogeneous state with little energy. We focus on the model skill in creating zonal jets from the SUF-MHM model with the balanced flux correction for q. Notice that without instability, the solution of the HM models will simply decay in time purely due to the damping effect without any external excitation. The results will be compared with the two-field BHW model simulations starting from the same initial configuration. Also, we are interested in the regeneration of the statistical energy transport as observed in the two-field BHW model for plasma edge turbulence using this simplified one-field system.

First, we illustrate the general flow evolution in time starting from the initial state with little energy. In Fig. 4, the time-series of the total energy and energy among zonal modes are compared in the same plot. Starting from the nearly zero initial energy, instability accumulates for a while before the rapid growth. This first stage of energy growth is among the nonzonal fluctuating drift waves due to the forcing effect simulating the internal instability. The generation of zonal jets happens in the next stage from the nonlinear transfer of energy between scales. It is confirmed from the observation that the rise of zonal energy always lags behind the total energy in time, implying the secondary instability taking place later after the energy in nonzonal fluctuations from linear instability gets accumulated. Besides, more turbulent features with the higher level of total energy are observed with the larger value of the density gradient κ.

For a detailed illustration of the energy exchange mechanism between different scales, we compute the energy spectra $Ek=k2|φ̂k|2$ for each wavenumber measured at several different time instants. First, the radial averaged spectra are plotted by taking the summation of all the radial modes with the same absolute value $k=kx2+ky2$⁠; second, the zonal spectra plot the energy in each zonal mode (k_x, 0). The radially averaged spectra have an emphasis for the energy in the nonzonal fluctuating modes, while the zonal spectra track the zonal jet structure. Figure 5 compares the energy spectra from the SUF-MHM and BHW models with κ = 0.5. The energy spectra can be compared with the corresponding flow snapshots shown in Fig. 1. The system starts with a flat energy spectrum. Then, the linear instability (modeled as unstable forcing in the SUF-MHM model) first takes place as the most dominant effect. The energy in the most unstable fluctuating modes (near k ∼ 1) begins to rise, while the energy in zonal modes (along the k_y = 0 axis) stays small. In the second stage, the nonlinear interaction effect takes over. The excited energy in the nonzonal fluctuating modes is transferred to the zonal direction through the nonlinear coupling mechanism. In the final stage, the energy cascades downscales and gets dissipated through the strong damping effect on small scales, while the small-scale energy in zonal modes is almost unchanged. Both the SUF-MHM and BHW models evolve in a similar fashion at every stage of the process, confirming the same energy mechanism generated from the explicit forcing added to the SUF-MHM equation.

For completeness, we also show the dynamical evolution of flow solutions from the two models in the more turbulent regime using κ = 1 in Fig. 6. Similar model performances with more turbulent dynamics are observed as in the previous case with κ = 0.5. The energy starts to grow first at the most linearly unstable modes and then transports up and down scales through the nonlinear effects. Energy in zonal modes always grows after the fluctuating energy gets accumulated large enough to trigger nonlinear energy exchanges. The snapshots of the relative vorticity ζ = ∇²φ are also plotted for the two models accordingly. Similar structures in the evolution are again observed. However, in this more turbulent regime, the one-field model maintains relatively stronger fluctuating vortices in the final state due to the strong downscale energy cascade and the lack of feedbacks from the stable branch.

Furthermore, we compare the equilibrium energy spectra achieved from both the two-field BHW model and the one-field SUF-MHM model in this more turbulent regime κ = 1. Figure 7 plots the equilibrium spectra in company with the results for the κ = 0.5 case shown in Fig. 1. In the stronger turbulent regime, there is larger difference in the smaller scale modes in the spectra between the two models. The better agreement in the zonal mode spectra implies that the additional energy in the small-scales of the SUF-MHM model comes from the induced fluctuating vortices in the flow field. The one-field model generates stronger energy in the small-scale modes (representing the smaller scale vortices observed in the snapshots in Fig. 6). Still, the essential energy structure is captured in the simpler SUF-MHM model. The large scale overall structures are captured by the forced one-field model with accuracy, and the decaying slopes in the inertial regime agree with each other through the two sets of models.

Next, we compute the zonal particle transport $Γ=∫ũñ$ using the approximation formulas (14) and (15) for the SUF-MHM model. Figure 8 gives the time-series and snapshots of the approximated particle flux in the leading order expansion at several measured time instants. Starting from the near-zero value in the initial state, the particle flux jumps to large values as the nonzonal drift wave states are excited. The zonal particle flux always reaches its strongest value before the zonal jets are completely formed. Finally, the dominant zonal structure blocks the strong zonal transport of the particle density. Again, the weaker turbulent case with κ = 0.5 shows a quite similar particle flux structure as the two-field BHW model (see Fig. 3 of Ref. 12) The stronger turbulence case κ = 1 with larger linear instability gets the small-scale structures maintained in time and more zonal particle transport in the one-field model case.

In the SUF-MHM model approximation, still we observe that the total particle flux is not entirely quenched at the final steady state. This may be a result of the insufficient modeling of the nonlinear coupling with the stable branch solution. As a result, stronger small vortices are maintained in the final flow field (specifically, for the more turbulent case κ = 1, as shown in the larger energy among small-scale modes shown in Fig. 7). The difference becomes smaller as the system gradually approaches the adiabatic limit.

The statistical higher-order moment feedback generated from the SUF-MHM model is computed here, which can offer an illustration for the drift wave–zonal flow interaction mechanism observed from the numerical results in Secs. II C and IV A. After the growth of energy among the unstable modes in the starting transient state, the growth in energy will finally get saturated when the nonlinear interactions take over as a dominant effect. It is found that the third-order moment feedbacks to the statistical energy equation usually have a central role in the balanced energy mechanism.^12,20 Specifically, we show that the one-field SUF-MHM model with the balanced flux correction can effectively enhance the zonal feedbacks from the higher-order moments.

To characterize the statistical energy exchange between different scales, we look at the statistics $EEk$ in the first two moments defined in each spectral mode $φ̂k$ of the electrostatic potential

where we use the pointed bracket, $⟨·⟩$⁠, to denote the statistical ensemble-averaged solutions. The statistical energy equation for the SUF-MHM model (13) can be derived by multiplying $φ̂k*$ on both the sides of the equation for each projected mode

where $QF,k=−⟨φ̂k*(∇⊥φ·∇q)k⟩+c.c.$ contains the important third-order statistical feedbacks to the system. Importantly, the higher-order terms due to the nonlinear interactions play the central role in mediating the growing unstable modes and driving the system to the final equilibrium. In fact, if the third moments on the left hand side of the above equation are purely ignored, the internal instability will lead to fast growth in energy among the unstable modes and fast decay in the other over-damped modes. The higher-order moments then transfer the growing energy in the unstable subspace to the stable one to get dissipated. Unfortunately, it is usually expensive to compute the third moments directly since it requires the inclusion of all the triad modes across the entire spectrum.²⁰ On the other hand though, in statistical equilibrium, the time derivative on the left hand side vanishes. Thus, the higher-order feedback can be calculated easily based on the equilibrium statistics in the lower moments on the right hand side.

In Fig. 9, negative values in the third moments show the effective damping to stabilize the linear instability, while the positive values represent the in-flow of energy in the linearly over-damped modes. The positive third moment feedback takes place along the zonal modes near k_y = 0, orthogonal to the most linearly unstable modes. For the unstable subspace, the third moments quench the unstable forcing effect on the fluctuating drift waves k_y ≠ 0 through the negative third moment feedback stabilizing the flow field. This characterizes the nonlinear energy mechanism in transferring the energy from the linearly unstable drift wave modes to the linearly stable zonal subspace. The large positive third moments at the zonal modes near k_y = 0 implies the transfer of energy to the zonal modes in generating strong zonal jet structures.

In the final test case, we consider the situation with a larger aspect ratio L_x: L_y = 5:1 for an extended x domain size in the SUF-MHM model simulations. A larger number of zonal jets are generated, and they interact with each other. Thus, this computational geometry will introduce richer phenomena for the time evolution of the jets. Figures 10 and 11 (as in Ref. 12 for the BHW model results with a different value of α = 0.5) give the corresponding plots for the simulation results. In the process to the creation of zonal flows, similar dynamical structures as the square domain case are observed. The most linearly unstable fluctuating modes first get excited. The energy transports to the zonal modes in the final state to form zonal jets. In the larger aspect ratio test, a larger number of zonal jets are generated. Compared with the zonal spectra in the square domain case, a larger number of zonal modes with higher level of energy appear here. At the same time, we can observe the reorganization of the fluctuating vortices gradually to turbulent zonal jets. The snapshots of the ion vorticity at several sampled time instants in the extended domain also show the similar self-reorganization of zonal jets from the homogenous initial state. In both the test cases, we can observe first the generation of drift waves from the explicitly added forcing. Then, quickly energy among these nonzonal drift waves transfers to the zonal states. Final zonal jet structures are developed.

Furthermore, we show the time evolution of the zonal mean flow $v¯=∂xφ¯$ in Fig. 12 for the two tested cases. Here, we can observe directly the formation and interaction between multiple jets. With the elongated x domain size, a larger number of zonal jets are generated. The zonal jets frequently merge and regenerate in time displaying much more complicated dynamics, and they do not have a characteristic spacing. This implies the multiscale structures formed by the groups of jets with a larger α = 5 in SUF-MHM consistent with the direct simulation results for the two-field BHW model with a large aspect ratio.¹² 

In the construction of both the BHW model (1) and the MHM model (2), the crucial role in using the balanced potential vorticity q is emphasized. The profound change introduced due to this simple correction in the potential vorticity function has been analyzed in detail from the selective decay principle⁴ and the secondary instability analysis.⁵ Here, we further characterize the important effect in introducing the balanced flux correction by applying this unstable forcing to the CHM model. We also show from the two-field HW model simulations that small-scale fluctuations are difficult to be quenched completely at low resistivity without adopting the balanced flux form.

To show the important role in the balanced flux correction, we consider the classical CHM model with the systematic unstable forcing which arises from the HW model with α ≫ 1

The potential vorticity q defined above is without the balanced correction to remove the zonal mean state in the second component. In this model, as we have shown in the secondary instability results,⁵ the transfer of energy to form the zonal structure is effectively weakened. It is found that only homogeneous turbulence can be generated from forcing the CHM model without the balanced flux correction.

In the left panel of Fig. 13, the time-series of the total energy and energy in zonal modes are compared using both the SUF-MHM model with the balanced flux $q=∇2φ−φ̃$ and the SUF-CHM model using q = ∇²φ – φ. For both the models, total energy grows in the first stage from the unstable forcing on the fluctuating modes. However, zonal modes are excited in the SUF-MHM model with a large amount of energy induced in the zonal modes. In contrast, the zonal state energy excited in the CHM model remains at a negligible level during the entire evolution time. This illustrates the lack of skill in the CHM model in properly generating the zonal jet structures from the nonlinear energy exchange mechanism.^4,19,21 These results with unstable forcing are in agreement with those for long time selective decay of the CHM model which is necessarily homogeneous,^4,21 while for the SUF-MHM model, the selective decay states are necessarily anisotropic with zonal jets of different wavelengths.⁴ 

For further comparison, Fig. 13 also shows snapshots of the flow vorticity field using the original CHM model simulations in the same tested parameter regime as done previously for the SUF-MHM model. This can be compared with the SUF-MHM model results shown in Fig. 1. Fluctuating drift waves are induced in the starting state in a similar fashion from the unstable forcing. Next, the flow breaks into complete homogeneous turbulence showing no clear zonal structures. Again, it confirms that the CHM model lacks the skill in properly transporting the excited nonzonal drift waves to the zonal direction to form zonal jets.

In the last part, we comment about the important role of the balanced flux model in guaranteeing the exact convergence to the MHM model at the low resistivity limit using a large value of α = 5 (see Fig. 4 of Ref. 17 for the performance with a slower converging rate using a moderate value of α). Figure 14 compares the long time performance of the BHW and MHW models in the numerical simulations. The MHW model gets much stronger energy in the fluctuations and stronger particle flux at the same time. By checking the snapshots at a long simulation time, the BHW model finally reduces to an almost zonal state with tiny zonal particle flux, while fluctuations are maintained in the MHW model with large persistent particle flux. The above results provide some explanation for the more turbulent structures that we observed from the previous SUF-MHM model results. There, the instability is added for all the time as an equivalent forcing. Thus, the mechanism for the zonal mean state to balance the drift wave growth is not fully modeled.

We proposed a systematic Galilean invariant unstable forcing form for the modified one-field Hasegawa-Mima model as a link to the drift wave instability naturally induced from the two-field Hasegawa-Wakatani system with finite resistivity. Instead of introducing the unstable forcing just empirically, the explicit forcing effect is derived in a precise way based on the leading-order growth rate from the HW model instability analysis at low resistivity. The direct unstable forcing on the potential vorticity in the HM model is validated based on the observation that the unstable branch solution of the HW model will converge to the exact HM model state at the low resistivity limit in the leading order, while the higher order terms decay at a much faster rate. We use the new SUF-MHM model to study key mechanisms in the drift wave–zonal flow interacting dynamics through direct numerical simulations. It is shown first in a qualitative way that persistent zonal jets are automatically generated from the excited drift waves due to the systematic unstable forcing, just following the same nonlinear interaction mechanism as observed in the two-field BHW model. Then, with a quantitative comparison for the statistics in equilibrium state variables, the SUF-MHM model displays significant skill in producing the same energy spectrum from the two-field BHW model. The self-generation of zonal jets is also illustrated through the important roles of the third-order moment feedback in the statistical dynamics (see Figs. 1 and 7). Furthermore, much richer dynamics at low resistivity with interacting jets and multiscale structures are discovered in simulations on an extended radial computational domain geometry (see Fig. 12 and Figs. 10 and 11 in Ref. 12).

In addition, the role of the balanced flux correction on the magnetic surfaces introduced in Refs. 12 and 17 is further confirmed in the low resistivity regime. First, we show the decay to homogeneous turbulence using the unstably forced CHM model without the balanced flux in the potential vorticity. Second, the direct simulation of the two-field MHW model shows its lack of skill in completely quenching the nonzonal fluctuation with large particle flux at the low resistivity limit. These numerical results confirm the theoretical analysis using the selective decay principle and secondary instability on a drift wave base state discussed in Refs. 4 and 5. In the next stage, we plan to consider the development of low-order statistical models for uncertainty quantification^20,22 based on the energy transfer mechanism analyzed here as a further direction. Using the statistical model reduction strategies developed in Ref. 20 with successful applications in the geophysical model,²³ it has the potential to show that the crucial statistical responses in the plasma turbulence can be captured in a more realistic setting with significantly fewer degrees of freedom.²⁴ 

The research of A.J.M. was partially supported by the Office of Naval Research under Grant No. N00014-19-S-B001. D.Q. was supported as a postdoctoral fellow on the grant.