The generations of zonal flow (ZF) and density (ZD) and their feedback on the resistive drift wave turbulent transport are investigated within the modified Hasegawa-Wakatani model. With proper normalization, the system is only controlled by an effective adiabatic parameter, *α*, where the ZF dominates the collisional drift wave (DW) turbulence in the adiabatic limit $\alpha >1$. By conducting direct numerical simulations, we found that the ZF can significantly reduce the transport by trapping the DWs in the vicinities of its extrema for $\alpha >1$, whereas the ZD itself has little impact on the turbulence but can only assist ZF to further decrease the transport by flattening the local plasma density gradient.

## I. INTRODUCTION

In the edge region of magnetically confined plasmas, the drift wave (DW) turbulence driven by the density and/or temperature gradients plays an important role in determining the level of cross field particle and heat transport. Suppressing the edge DW turbulence is crucial to reduce the losses of particles and energy from the bulk of plasma and thus improve conditions for nuclear fusion. It has been shown that the zonal flow (ZF) can be a great candidate in regulating the anomalous transport^{1,2} since it drives no transport (from $E\xd7B$ drift) and thus can be a safe repository for the energy released by gradient driven turbulence. The ZF here is referred to as the poloidally and toroidally symmetric ($k\theta =kz=0$) and zero/low-frequency vortex modes with finite radial scale (*k _{r}* finite). Moreover, we focus on the two-dimensional slab model for the edge plasmas and thus ignore the toroidal effect, which can lead to the geodesic acoustic mode (GAM or finite frequency ZFs) in the high-q tokamak edge region.

^{3,4}(Note that the suppression of turbulence by GAM can be different.

^{5,6})

The generation of ZF can be seen as a “secondary” instability in DW turbulence via Reynolds' stresses, where the “primary” instability occurs at a linear stability threshold with a zero background flow and causes growth of DWs. The secondary instability can be the modulational instability,^{7,8} an inverse energy cascade,^{9} coherent hydrodynamic instability,^{10} or a resonant-type interaction between the ZF and modulations of the small scale DW turbulence.^{11} Once the ZFs are fully developed, they can largely suppress the DW turbulence through a widely assumed effect of shearing the eddies.^{12–14}

Taking into account the fact that the realistic models for the DW/ZF system are complicated, it is useful to consider a simple non-trivial model, namely, modified Hasegawa-Wakatani^{15} (MHW) model to study the feedback of ZF on the DW turbulence. Compared with the original Hasegawa-Wakatani^{16} (oHW) model, the MHW removes the zonal-averaged components from the Boltzmann response, as the zonal mode has no dependence on the direction parallel to the magnetic field. As a result, the ZFs are not dissipated by the adiabatic term and can be dominant. It is worth noting that even in the oHW model, the generation of large ZF is also possible and can be important for an extremely large adiabatic parameter,^{17} which, however, makes the simulations time-consuming because of the small linear instability.

In the MHW model, a zonal density (ZD), in addition to the ZF, will also be generated due to the non-zero non-linear advection of density resulting from the non-adiabatic response of electrons. Although the generation of ZF and its feedback on the DW turbulence have been well studied, the effect of ZD on turbulent transport has been overlooked, which can be important by modulating the equilibrium plasma density.^{18} Moreover, there is no detailed study of the role of ZF in the transport of radial particle flux within the MHW model yet. Therefore, in this paper, we will examine the effects of both the ZF and ZD on the radial particle flux. To identify their impacts separately, we will not only consider the MHW model but also employ two of its modifications by artificially removing the ZF and ZD, respectively.

By introducing proper normalization, the MHW model is only controlled by the effective adiabatic parameter, $\alpha \u223ckz2vth2/\nu ei|\omega |$, where *v _{th}* is the electron thermal velocity,

*ν*is the electron-ion collisional rate, and $|\omega |$ is the characteristic DW frequency. It was shown that the DW/ZF system has a bifurcation behavior

_{ei}^{15}depending on

*α*: in the adiabatic limit $\alpha >1$, the DW/ZF system is dominated by ZF, whereas, in the hydrodynamic limit $\alpha <1$, turbulence is dominant. This reduction of the ZF effect is related to the enhancement of turbulent viscosity and the reduction of residual flux when the plasma response passes from the adiabatic to hydrodynamic limits. The underlying physics of the shear flow collapse and DW turbulence enhancement in the hydrodynamic limit was investigated in detail in Ref. 19 from the point of view of Reynolds work, wave energy and momentum fluxes, potential vorticity mixing, and predator-prey model. It is worth pointing out that in Ref. 19 the transition of the adiabatic parameter from $\alpha >1$ to $\alpha <1$ is linked to the approaching to the density limit, which is accompanied by a weakening of the edge shear lays and degradation of particle confinement. From direct numerical simulations (DNS), we find that ZD itself has little impact on the turbulence but can only assist ZF to reduce the transport. Therefore, in this paper, we mainly focus on the adiabatic limit $\alpha >1$, where the ZF dominates the turbulence. Keep in mind that $\alpha >1$ means that the standard adiabatic parameter $\alpha \u0302\u2261kz2vth2/\nu ei$ is greater than the characteristic frequency of DW.

## II. MODELS

The models we will consider are based on the MHW equations, which capture the essential characteristics of turbulent transport in the DW/ZF system. If we consider a slab geometry with a constant equilibrium magnetic field $B=B0ez$ (where $ez$ is a unit vector in the *z*-direction), an exponentially decaying background plasma density $n0\u221d\u2009exp(\u2212x/Ln)$ in the *x* (radial)-direction, parallel electron resistivity, and cold ions, the MHW equations are in the quasi-two-dimensional (2D) form

Here, all physical quantities have been normalized as $x/\rho s\u2192x,\kappa \omega cit\u2192t,\u2009e\varphi /Te\kappa \u2192\varphi ,\u2009n1/n0\kappa \u2192n$, where $\varphi $ and *n*_{1} are, respectively, the electrostatic potential and density fluctuations, $\rho s=Cs/\omega ci$ is the ion sound gyroradius, $Cs=Te/M$ is the ion sound speed, *T _{e}* is electron temperature,

*M*is ion mass, $\omega ci=eB0/Mc$ is the ion cyclotron frequency,

*e*is the elementary charge,

*c*is the speed of light, and $\kappa =\u2212\rho s\u2202\u2009ln(n0)/\u2202x\u2261\rho s/Ln$ is a constant as the energy source for DW instabilities. The effective adiabatic parameter is redefined as $\alpha \u2261\alpha \u0302/\omega ci\kappa $ ($\omega ci\kappa $ is the characteristic frequency of DW), which controls the essential physical behavior. The only difference between the MHW and oHW models is in the resistive coupling term: in the MHW model, it is determined only by the non-zonal components $f\u0303=f\u2212\u27e8f\u27e9$, where $\u27e8f\u27e9\u2261\u222b0Lyfdy/Ly$ denotes the integration along the poloidal line at a given radial location; whereas, in the oHW model, all the components contribute to it. The non-linear advection terms are expressed in the Poisson brackets ${a,b}=(\u2202a/\u2202x)(\u2202b/\u2202y)\u2212(\u2202a/\u2202y)(\u2202b/\u2202x)$. The dissipation terms of the form $\u22072A$ with coefficients

*D*and

*ν*are added to the equations for numerical stability.

In the oHW model, there are two different limits depending on *α*: the adiabatic limit as $\alpha \u2192\u221e$ for collisionless plasma, where the density fluctuations become enslaved to the electrostatic potential fluctuations through the Boltzmann relation and oHW equations are reduced to Hasegawa-Mima equation;^{20} and the hydrodynamic limit for $\alpha \u21920$, in which the equations are decoupled similar to the 2D Navier-Stokes equations. In the MHW model, these two limits correspond to, respectively, the ZF-dominated regime and the turbulence-dominated regime,^{15} where the ZF absorbs almost all the kinetic energy from the DWs in the saturated state in the former case, while the turbulence is more isotropic in the latter. In the following, we mainly focus on the adiabatic limit, i.e., $\alpha >1$, where the ZF dominates in the system.^{15}

From Eqs. (1) and (2), we can see that, in the MHW model, the ZF (ZD) can be generated via the non-linear advection of vorticity (density)

where $V(x)\u2261\u2202\u27e8\varphi \u27e9/\u2202x$ is the ZF velocity, $N(x)\u2261\u27e8n\u27e9$ is the ZD, and $v\u0303=ez\xd7\u2207\varphi \u0303$. Equation (3) shows that the non-adiabatic response of electrons is important in the generation of ZD. We are particularly interested in the impacts of the ZF and ZD on the particle flux in the radial direction

where $\u222bfdx\u2261\u222b0Lx\u222b0Lyfdxdy/LxLy$. We will also monitor the total energy and the energy stored in the zonal-components

Moreover, in order to distinguish the effects of the ZF and ZD on Γ_{n}, we will not only use the oHW and MHW models but also consider the modifications of MHW equations by artificially removing the ZF (we will call it MHW w/o ZF) and ZD (referring as MHW w/o ZD) components, respectively.

## III. SIMULATION RESULTS

Equations (1) and (2) will be numerically solved by using a pseudo-spectral Fourier code in a square box with doubly periodic boundary condition in Dedalus.^{21} The size of the box is $Lx=Ly=10\pi $ so that the lowest wavenumber is $\Delta k=0.2$, while the number of the modes is chosen as 256 × 256. The time integration algorithm is the fourth-order Runge-Kutta method with a time step $\Delta t=10\u22123$. We mainly focus on the non-linearly saturated particle flux with fixed dissipation coefficients of $\nu =D=10\u22124$. For the oHW and MHW models, we start the simulations from one small perturbation, while in the cases of the MHW model w/o ZF and ZD, we employ a different initial condition with zero zonal-averaged components.

In all the cases, we observe that the fluctuations first grow linearly due to the resistive DW instability, where the modifications will not affect the growth of the DW in the linear regime (e.g., see Fig. 1 for *α* = 2). Note that, in the late linear stage, the simulations of oHW and MHW are dominated by the mode $(kx,ky)=(0.2,1)$, while those of MHW w/o ZF/ZD are governed by $(kx,ky)=(0.2,0.6)$. The growth rates from the simulations are, respectively, 0.025 and 0.057, which agree with the analytic result of resistive DW in the adiabatic limit, $\gamma =ky2k2/(1+k2)3\alpha $ in our normalized units with $k2=kx2+ky2$.

However, when the amplitudes of our normalized DWs become the order of unity, they will undergo secondary instabilities and thus generate the ZF and/or ZD (as an example, the growth rate of ZF can be $\u223c0.1$ in our normalized unit^{22,23} for *α* = 2, i.e., $\gamma zf\u223c0.1\kappa \omega ci$). The ZF and/or ZD will, in turn, affect the DW, and the system will be saturated eventually when the particle flux Γ_{n} approximately balances the dissipation, $D\alpha =\alpha \u22121\u222b[\alpha \u0302(n\u0303\u2212\varphi \u0303)]2dx$, due to the parallel resistivity. From Fig. 1 we see that, the saturated particles fluxes are quite different, where the time-averaged Γ_{n} at the saturated state from *t* = 400 to *t* = 1000 are 0.36 (oHW), $7.4\xd710\u22123$ (MHW), $2.2\xd710\u22122$ (MHW w/o ZD), and 0.4 (MHW w/o ZF), respectively. This is consistent with the standard recognition that the ZF can largely suppress the anomalous transport, whereas it illustrates that the ZD plays an important role in the MHW model (as seen from the plots of MHW and MHW w/o ZD) but not in the oHW model (e.g., see the comparison between the oHW and MHW w/o ZF). Therefore, we conclude that the ZD can reduce the transport only via its synergy with the ZF. This conclusion has also been confirmed by the simulations for other adiabatic parameters $\alpha >1$. For example, for *α* = 4, the time-averaged Γ_{n} in the non-linear saturated stage are, respectively, 0.21 (oHW), $2.6\xd710\u22123$ (MHW), $3.8\xd710\u22123$ (MHW w/o ZD), and 0.24 (MHW w/o ZF). Note that the larger adiabatic parameter means that the electron response to the electrostatic potential is more adiabatic and, thus, the overall particle flux is smaller.^{19}

However, on one hand, in the hydrodynamic limit $\alpha \u226a1$, even the MHW model is strongly turbulent and thus both the ZF and ZD are small. On the other hand, the generation of ZD largely depends on the non-adiabatic response of electrons, $h\u0303k,\omega \u2261(n\u0303\u2212\varphi \u0303)k,\omega =i(\omega \u2212\kappa ky)\varphi \u0303k,\omega /(\alpha \u2212i\omega )$, which is small for larger $\alpha \u226b1$. Therefore, it is expected that the ZD effect on the transport begins to drop when *α* is above some threshold values (e.g., the ZD reduces the time-averaged particle flux by a factor of 3 for *α* = 2 but only by a factor of 1.5 for *α* = 4). Therefore, the role of ZD in the suppression of transport can be important only for a large but finite range of *α* in the MHW model. However, searching for such a range of adiabatic parameter needs time-consuming simulations to scan the parameter space and is beyond the scope of this paper.

To see the different spatial behavior of the ZD in the oHW and MHW models, we plot the saturated density in Fig. 2. (The different spatial structures of the electrostatic potential are similar and can be found, e.g., in Ref. 15.) It shows that the ZD is more prominent in the MHW, where the potential energy stored in the ZD, $Enz$, is approximately half of the total potential energy, *E _{n}*, whereas in the oHW model, $Enz/En\u22480.1$ and so is $Ekz/Ek$. This demonstrates that although the ZF/ZD can still be generated in the oHW model, their effects are negligibly small. We notice that in the MHW w/o ZF, the saturated density is similar to Fig. 2(a), but the ZD stores less energy $Enz/En\u22480.03$. Moreover, Fig. 2 illustrates that the saturated level of $n\u0303$ (and thus $\varphi \u0303\u2248n\u0303$ in the adiabatic limit) is largely reduced by the ZF. However, this suppression of turbulence is not simply due to the concentration of energy to the ZF since it was also observed in the resistive DW system with a constant externally applied ZF.

^{24}In Sec. IV, we will show that the reduction of transport is because of the trapping of turbulence near the extrema of ZF.

## IV. ZONAL FLOW AND DENSITY IN THE MHW MODEL

In the last section, we showed that the ZD can assist the ZF in suppressing the saturated Γ_{n} in the MHW model, and we will examine the underlying physics in detail in this section. We will focus on a linearized problem of DW with quasi-static ZF velocity and ZD in simple sinusoidal profiles

where *q _{v}* and

*q*are the products of integer numbers and $\Delta k$ so that

_{n}*V*(

*x*) and

*N*(

*x*) are periodic, and $\Delta \phi $ is a constant phase shift needed to be determined from the simulation results. Then the governing equation for $\varphi \u0303=\varphi \u0302(x)\u2009exp(\u2212i\omega t+ikyy)$ is

where $\kappa eff=1\u2212N\u2032$, and $\omega \u0303=\omega \u2212kyV$. Here, $\omega \u2261\omega r+i\gamma $ is complex with *γ* being the growth rate.

We first consider a long-lasting problem of linear resistive DW with ZF only, i.e., $N0=0$. For the inhomogeneous profile of ZF, the fluctuations can be trapped near the extrema^{24–26} of *V*(*x*), e.g., see Fig. 3(a). These trapped structures move along the positive (negative) y-direction near the maximum (minimum) of *V*(*x*) and are reflected back and forth from the “slopes” of ZF in the x-direction. (It is worth noting that the propagation direction of trapped structures depends on the sign of equilibrium magnetic field so that if we flip the sign of magnetic field, the propagation direction of trapped structures will also change.) In fact, the trapped structures near the extrema of *V*(*x*) can be interpreted as eigenmode solutions of Eq. (7) with zero boundary conditions, which is a product of traveling wave in y-direction and standing wave in x-direction. Near the maximum and minimum of *V*(*x*), the eigenmode solutions are growing and decaying, respectively. Therefore, if the ZF is externally applied, DWs will be trapped only near the maximum of *V*(*x*), while those near the minimum will disappear. However, if the ZF is self-generated as in our simulations, trapped DWs also appear near the negative peak of *V*(*x*) as shown in Fig. 3, because fast self-generated ZF needs strong turbulence as an energy supplier.^{27}

The eigenmode solution can be interpreted within the eikonal approximation, which assumes that the DW wavelength is much smaller than that of ZF. Although this assumption limits the range of both DW and ZF parameters, this approach can provide some insights of the characteristics of non-linear interaction of DW and ZF. In this approach, the wave packet near the maximum of *V*(*x*) could be considered as an effective “particle,” dynamic of which is described by the “Hamiltonian”

and the canonical variables *x* and *k _{x}* (e.g., see Refs. 1, 11, 28, and 29). Here, we ignored the growth rate

*γ*in

*ω*and the last term in the bracket $R\u2261i(\omega \u0303\u2212\kappa eff)/(\alpha \u2212i\omega \u0303)$ in Eq. (7) for the adiabatic limit. To have an eigenmode (trapped) solution of the wave packet while keeping $\omega (kx,x)$ as constant, the motion of the particle should be bounded by two turning points corresponding to

*k*= 0. This is only possible in the vicinity of the extrema of

_{x}*V*(

*x*) in the form of Eq. (6).

Equation (7) is solved numerically both near the maximum and minimum of *V* for *α* = 2, where the results of growth rates and contours of the equipotential are plotted in Fig. 4. The size of the computation box $L\u0303x$ is chosen that there is only one extrema of *V*(*x*) for given *q _{v}*, while we vary

*k*in the simulations. In fact, the trapped fluctuations in Fig. 3(a) both near the maximum and minimum of

_{y}*V*are superpositions of different modes (

*k*), where the mode amplitude linearly decreases to zero when $ky\u21921$ in the former case, but it follows a slower decreasing curve till $ky\u223c3$ in the latter. The averaged wavenumber, $k\xafy\u2261\u222bkyEkdk/\u222bEkdk$ where $Ek=|\varphi \u0303k|2$, of the fluctuations near the maximum and minimum of

_{y}*V*(

*x*) are, respectively, $k\xafy\u22480.4$ and $k\xafy\u22481$, which match the visualization in Fig. 3(a). These averaged wavenumbers are used in the numerical simulation of Eq. (7).

From Fig. 4(a), we see that the growth rates of the modes near the maximum of *V*(*x*) in the top panel only slightly change (either increase or decrease) concerning *V*_{0}. However, increasing *V*_{0} can cause the trapping of the DWs. Particularly, for $V0=0$, the plasma is homogeneous and the DW is not localized at all, whereas, for large $V0\u2a9e1$, the DWs are localized as seen in Fig. 4(b). This is consistent with the simulation observations that the trapping of fluctuations begins once the large ZFs, $V0\u2a9e1$, are generated. Moreover, the linear solutions illustrate that the real frequency of these eigenmodes is positive approaching $kyV0$ for $V0\u2a9e1$ as indicated in Eq. (8). For *V*(*x*) having a negative peak, the least stable eigenmode is found and plotted in the bottom panel. We note that the real frequency for these modes also approaches $kyV0$ but with negative sign, agreeing with the numerical observations for Fig. 3(a). It is also worth noting that the ZF does not tilt the eddies of the contour of equipotential as shown in Fig. 4(b), which is consistent with the case, where the DWs are trapped by the inhomogeneity of the equilibrium plasma density.^{30} Therefore, we conjecture that the suppression of the particle flux by the ZF is mainly due to the trapping of the DWs rather than by shearing the eddies.

We then consider the impact of ZD on the linear DW and thus the non-linear turbulent transport by taking $N0\u22600$. As a result, the structure of the Hamiltonian $\omega (kx,x)$ in Eq. (8) and thus the trapped structures will be changed by $N\u2032$ through *κ _{eff}*, where the fluctuations accumulate near the extrema of $\omega (kx=0,x)$. These trapped structures have been observed in the simulations, e.g., see Fig. 3(b). From Fig. 3(b), we see that the maximum (minimum) of

*V*coexists with the minimum (maximum) of

*κ*, which results in a relatively smaller growth rate of DW compared with the case of $N0=0$ and thus indicates a reduced transport. From Eq. (3), we see that ZD evolves faster than ZF due to the viscous damping since ZD has larger effective wavenumber

_{eff}^{26}as seen from Fig. 3(b). This has been observed in the simulations, e.g., see Fig. 5(a). In such a process, ZD will exchange energy with DWs, which alters the structure of $\omega 0(x)$ in Eq. (8) and thus the trapping of DWs. As a result, trapping DWs can penetrate radially and accumulate at the locations, where the maxima of ZF and ZD coexist (and thus smallest

*κ*), e.g., see Fig. 3(c) compared with Fig. 3(b). To see how the ZD affects the energy division in the turbulent and zonal components, we plot $Ekdr=Ek\u2212Ekz,\u2009Endr=En\u2212Enz,\u2009Ekz/Ek$, and $Enz/En$ in Fig. 5. It shows that the presence of ZD can slightly increase $Ekz/Ek$ that the ZF is more dominant and thus reduces the fluctuations of DWs ($Ek,ndr$). The shifts of $Ekz/Ek$ can be also understood by considering a decreased

_{eff}^{15}

*κ*near the maximum of

_{eff}*V*.

## V. DISCUSSIONS AND CONCLUSIONS

In this paper, the generations of ZF and ZD and their roles in suppression of the radial particle flux in the resistive DWs turbulence are examined based on the MHW model and its modifications. (Here, we focus on the adiabatic limit for the resistive DWs, where the ZF can largely regulate the turbulence.) Numerical simulations of these models are conducted, from which we observe that when the amplitudes of DWs become the order of unity as a result of growth from the linear state, the secondary instabilities will lead to the generation of the ZF and/or ZD, which will, in turn, suppress the DW turbulence and thus reduce the radial particle flux. We find that the ZF itself can significantly reduce the particle flux as shown in Fig. 1, whereas the ZD can reduce the transport only through its synergy with ZF.

The underlying physics has been examined, and it is shown that the generation of large ZF can trap the DWs in the vicinities of its extrema despite the sheared flow. We conjecture that such trapping of fluctuations is the reason for the reduction of the radial particle flux rather than the conventional picture of shearing of eddies by ZF. First of all, Fig. 4 (as well as Ref. 30) shows that the velocity shear itself cannot largely change the growth rate of the eigenmode solution of resistive DW and does not tilt the eddies. Therefore, the conventional picture of shearing eddies by ZF does not apply here. Second, the trapping of DWs in this paper is due to the inhomogeneity of ZF, which grows from DWs in the non-linear stage. When ZF is strong enough, the DW fluctuations will be transited into the trapping regime. During such a transition, the conventional picture of shearing eddies can work, where the strong ZFs can break big DW eddies into small ones. However, this is only because our ZF is self-generated when DWs are large. If the ZF is externally applied, the fluctuations will also be trapped as shown in Fig. 4 (see also, e.g., Ref. 25). In such a case, eigenmode solutions will grow from small fluctuations without breaking of large eddies, where the eigenmode with the largest growth rate as shown in Fig. 4 will finally dominate. Therefore, the suppression of the turbulence is not due to the shearing of DW eddies but through the trapping of fluctuations. Last but not least, we notice that the turbulent level is largely suppressed by the ZF, either externally applied or self-generated.^{24} Therefore, the suppression of turbulent transport is also not simply due to the concentration of fluctuation energy to ZFs.

On the other hand, the ZD can flatten the local effective plasma density gradient (instability driver) in the regions, where the DWs are trapped [notice that the flattening of ZD on the background density is small $N/n0\u223c2\kappa \u226a1$ as shown in Fig. 2(b)], which results in a smaller linear growth rate and thus reduction of turbulence. However, in the absence of ZF, the system is strongly turbulent that the impact of ZD on the DW instability and eventual non-linear transport is negligible. Our conclusion of the ZD impact on the transport holds for large but finite *α* in the MHW model. When $\alpha \u2192\u221e$, the ZD is negligible due to the adiabatic response of electrons, while for small $\alpha \u2009\u2a9e\u20090.1$ in the hydrodynamic limit, both the ZF and ZD become unimportant since the turbulence is dominant. Therefore, there should exist a range of *α* that the impact of the ZD on the particle flux is important in the MHW model. Based on this characteristic, it is interesting to discuss the ZD effect in the oHW model. Recalling that ZD can affect the particle flux only via its synergy with ZF, it can be important only when ZF is dominant for an extremely large adiabatic parameter, $\alpha \u2009\u2a9e\u20091000$, in the oHW model.^{17} However, for such large *α*, the ZD is negligibly small and thus hardly makes any impact. Therefore, it is conceivable that ZD has little impact on the turbulent transport for the whole *α* space in the oHW model.

Note that both the oHW and MHW models employ the Boussinesq approximation and assume a constant logarithm equilibrium plasma density gradient. However, considering that we are studying the ZD impact on the particle transport downhill to the plasma density gradient, a relaxation of these assumptions may be important. As an example, the resistive DW instability with non-Boussinesq approximation was considered in Ref. 31 when dealing with the edge plasma with a very steep background density gradient. A more refined model and a systematic study are needed and will lead to interesting directions for future works.

## ACKNOWLEDGMENTS

The authors thank H. Zhu for bringing to the authors' attention the Dedalus and an input script. This work was supported by the U.S. Department of Energy, Office of Science, Office of Fusion Energy Sciences under Award No. DE-FG02-04ER54739 at UCSD.

## DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding author upon reasonable request.