Dimits shift, avalanche-like bursts, and solitary propagating structures in the two-field flux-balanced Hasegawa–Wakatani model for plasma edge turbulence

Qi, Di; Majda, Andrew J.; Cerfon, Antoine J.

doi:10.1063/5.0018943

We show that the recently introduced two-field flux-balanced Hasegawa–Wakatani (BHW) model captures the key features of drift-wave turbulent transport mediated by zonal flows observed in more complete and accurate gyrokinetic simulations, such as the existence of a nonlinear upshift of the threshold for drift wave turbulence driven transport, often called the Dimits shift, as well as non-local transport with avalanche bursts and solitary propagating structures. Because of the approximations made in the BHW model, these observations are made for the particle flux instead of the heat flux more commonly studied in ion temperature gradient (ITG) driven turbulence in fluid or gyrokinetic codes. Many of these features are not seen in other Hasegawa–Wakatani models, which confirm the critical role of the electron dynamics parallel to the magnetic field lines. To address questions regarding the role of boundary conditions on the drift-wave zonal flow dynamics, we apply our model to both a channel domain geometry and the more typical doubly periodic geometry. We only observe strong soliton-like solutions in the particle flux for the channel geometry, in the vicinity of the boundaries, where strong velocity shear and density gradients are generated, which are absent in the doubly periodic simulations. Changing the aspect ratio of the simulation domain also has a significant effect. In domains which are elongated in the radial direction, more complex multi-scale dynamics takes place, with multiple zonal jets interacting with each other, and large scale avalanches.

I. INTRODUCTION

The zonal flow–drift wave interaction is known to play a fundamental role in setting the observed level of radial particle and heat transport in magnetically confined plasmas.^1–7 The self-organization of the plasma from a turbulent state to a more laminar state with zonal flows can lead to a sharp reduction of turbulent transport, sometimes even a full suppression.^8,9 Furthermore, regimes with a high level of turbulent transport are not accurately modeled by theories relying on a local diffusion paradigm¹⁰ and are instead better described by approaches allowing for non-local, nondiffusive transport, characterized by avalanche bursts and solitary traveling structures,^10–16 in which zonal flows also play a key role.^13,16–18 While these observations rely on a large amount of numerical evidence, and the importance of these physical processes in modeling and predicting transport in magnetic confinement fusion devices is not disputed, theoretical understanding of the fundamental mechanisms at play and reduced models with accurate yet inexpensive predictive capabilities remain scarce.

Although reduced fluid models can rarely be derived self-consistently from first principles for the values of the plasma parameters measured in magnetic confinement fusion experiments, they have been useful to identify the elementary building blocks required for the emergence of zonal flows from drift wave turbulence,^19–23 the associated reduction of turbulence driven transport,^16,24,25 and the non-local transport characterized by ballistic bursts called avalanches^17,26 and long-time-coherent, soliton-like traveling structures, for which the term “ferdinons” was recently coined.^16,27 In this article, we demonstrate that our recently proposed flux-balanced Hasegawa–Wakatani (BHW) model^22,23 is one of the simplest models which captures all of these characteristic features of turbulent transport in magnetic confinement fusion devices. It is, therefore, a powerful model not only to study the complex interplay between the fundamental physics mechanisms responsible for the observed level of transport, but also to test and benchmark mathematical methods for model reduction as well as uncertainty quantification, which, if successful, may then be tested on the more accurate gyrokinetic model and the significantly more computationally intensive gyrokinetic simulations.

There have been three different Hasegawa–Wakatani (HW) models introduced in the plasma physics community: the original HW (OHW) model,²⁸ the modified HW (MHW) model,²⁰ and our flux-balanced HW model.²² The BHW model differs from the other two by its treatment of the electron response parallel to the magnetic lines, which is constructed so as to ensure, in contrast to the OHW and MHW models, that there is no net radial flux of electrons in the adiabatic electron limit. In this article, we will show that certain characteristics of magnetized turbulent transport discussed in the previous paragraph are only captured by the BHW model, and not the other HW models. This is, in particular, true for what is known as the Dimits shift, namely, the nonlinear upshift in the critical gradient required for the onset of significant turbulent transport.^9,16,25 Our results, thus, show that the electron dynamics is a critical component to the existence of a true Dimits shift, with a sharp transition in the observed level of transport from a turbulent flow regime to a regularized zonal flow regime, which is absent in the other HW models. This result is in agreement with the conclusion St-Onge reached with the Terry–Horton model.²⁴ Observe that the Dimits shift usually refers to the nonlinear upshift for the onset of strong heat flux and is usually considered in the context of the ion temperature gradient (ITG) instability.⁹ Since our BHW model does neither have a temperature gradient nor a heat flux, we will study the Dimits shift for the onset of strong particle transport, driven by a density gradient drift instability, as was done in other models.^24,25 In Appendix A, we highlight the similarities between the BHW model and a recently proposed ITG fluid model in a Z-pinch geometry,¹⁶ and the comparable properties of the particle flux in the BHW model and of the heat flux in that ITG fluid model.

In addition, this article will investigate, within the framework of the BHW model, the role that boundary conditions play on the dynamics. The dominant point of view in the magnetic confinement fusion community has long been that, since turbulent structures are very much smaller than the system size in the directions perpendicular to the magnetic field, the local turbulence properties can be expected to mostly depend on local gradients and are not affected by boundaries in a meaningful way. As a result, a large number of simulations are done with periodic boundary conditions in the directions perpendicular to the magnetic field, both for reduced fluid models and for detailed kinetic models.^{20,24,29–32} However, the non-local nature of the transport characteristics and the large scale structures which emerge as a result of the plasma self-organization raise questions regarding the validity of this approach. Recent results with non-periodic boundary conditions¹⁵ and with domains with different aspect ratios³³ suggest that boundary conditions indeed have an impact on the resulting turbulence and transport. In this article, we will therefore study the emergence of zonal flows and of large scale structures, and the existence of a Dimits shift as well as of non-local transport for two different boundary conditions: a channel domain geometry, where the radial component of the velocity and the density perturbation are set to zero at the inner edge and outer edge of the computational domain,³⁴ and the periodic boundary conditions more commonly used for such study. Furthermore, we also investigate the dependence of our results on the aspect ratio of the computational domain. A computational domain which is elongated in the radial direction includes many more intermediate scales, and their interaction creates richer dynamics with the generation of multiple jets interaction and stronger intermittent bursts in the transient regime.

In the remainder of this article, the flux-balanced Hasegawa–Wakatani model and its fundamental properties are described in Sec. II. The existence of a true Dimits shift and a sharp transition from a turbulence dominated, high flux regime to a regime with strong zonal structures and very small particle flux are studied for the channel domain geometry in Sec. III. The avalanche-like structures in the turbulent transport of particles analyzed in Sec. IV, with a comparison between different HW models and different aspect ratios of the computational domain. In Sec. V, we revisit the simulations of Secs. III and IV with doubly periodic boundary conditions imposed instead of the channel geometry, and again compare the MHW and BHW models. Finally, our results are summarized in Sec. VI.

II. THE FLUX-BALANCED HASEGAWA–WAKATANI MODEL IN PLASMA EDGE TURBULENCE

Below, and for the remainder of the article, we will identify the distinct contributions from the zonal states and from the non-zonal fluctuations by decomposing the electrostatic potential $φ$ and particle density n into the zonal mean states $\bar{φ}, \bar{n}$ and their fluctuations about the zonal mean $\tilde{φ}, \tilde{n}$ as follows:

φ = \bar{φ} + \tilde{φ}, n = \bar{n} + \tilde{n}, \bar{f} (x) = L_{y}^{- 1} \int_{0}^{L_{y}} f (x, y) d y .

Generically, the Hasegawa–Wakatani models^20,22,28 describe the drift wave–zonal flow interactions in a two-field coupled system in a shearless slab geometry. The system supports unstable drift waves driven by a background ion density gradient and a non-adiabatic electron response parallel to the magnetic field. The flux-balanced Hasegawa–Wakatani (BHW) model was recently introduced in Ref. 22 and is conveniently expressed in terms of the flux-balanced potential vorticity $q = \nabla^{2} φ - \tilde{n}$ and the density fluctuation n, as follows:

\frac{\partial q}{\partial t} + u \cdot \nabla q - κ \frac{\partial φ}{\partial y} = μ Δ q,

(1a)

\frac{\partial n}{\partial t} + u \cdot \nabla n + κ \frac{\partial φ}{\partial y} = α (\tilde{φ} - \tilde{n}) + μ Δ n .

(1b)

Here, $φ$ is the electrostatic potential, n is the density fluctuation from the background density profile $n_{0} (x)$ ⁠, and $u \equiv \nabla^{⊥} φ = (- \partial_{y} φ, \partial_{x} φ)$ is the E × B velocity field, where x can be interpreted as the coordinate for the radial direction, and y the coordinate for the poloidal direction, if our shearless slab geometry is interpreted as an approximation of the edge of a magnetic confinement device. The constant background density gradient $κ = - \nabla ln n_{0}$ is a consequence of the assumption of an exponential background density profile in the plasma region of interest: $n_{0} (x) = exp (- κ x)$ ⁠. α is a constant which is often called the “adiabaticity parameter”^20,21 and which is proportional to the inverse of the parallel resistivity. The terms $μ Δ q$ and $μ Δ n$ can be viewed as dissipation due to the collisional diffusion of the electrons perpendicular to the magnetic field²² and have the same dissipation coefficient μ for the vorticity and for the density.³⁵

We stress that in the BHW model (1), the poloidally averaged density $\bar{n}$ is removed from the potential vorticity q. By doing so, we were able to prove rigorously,²² as well as confirm numerically,^22,23 that in the limit of adiabatic electrons, $α \to + \infty$ ⁠, the BHW model converges to the following one-field model:

\frac{\partial q_{MHM}}{\partial t} + \nabla^{⊥} φ \cdot \nabla q_{MHM} - κ \frac{\partial φ}{\partial y} = μ Δ q_{MHM}, q_{MHM} = \nabla^{2} φ - \tilde{φ},

(2)

which is called the modified Hasegawa–Mima (MHM) model,¹⁹ where we draw the attention of the reader to the fact that the zonally averaged potential $\bar{φ}$ is absent in the definition of the MHM potential vorticity: $q_{MHM} = \nabla^{2} φ - \tilde{φ}$ ⁠. In the BHW model, there is therefore no net radial transport of electrons in the limit of adiabatic electrons, which is what one would expect physically.^19,36–38 In contrast, when the modified Hasegawa–Wakatani (MHW) model is expressed in terms of the potential vorticity and the density, the “unbalanced” potential vorticity $q_{MHW} = \nabla^{2} φ - n$ appears, which contains the mean state contribution $\bar{n}$ ⁠. In the formal limit “ $α = + \infty$ ⁠,” the MHW model also agrees with the MHM model, but the convergence of the MHW model to this limit is problematic, in the sense that for α large but finite, the dynamics of the MHW model can be significantly different from that of the MHM model, as we have shown previously in Refs. 22 and 23, with a considerable amount of kinetic energy maintained in the small scale turbulence in the MHW model, which is not fully quenched by the zonal structures. This subtle difference will play a major role in the remainder of this article as it is the reason why a true Dimits shift will be observed in the BHW model, with a sharp transition between a high transport regime and a low transport regime, and not seen in the same way for the MHW model.

A. Boundary conditions for the computational domain geometry

As discussed in the Introduction, the two-dimensional flow is usually defined on a rectangular domain, i.e., $x = (x, y) \in D = [- L_{x}, L_{x}] \times [- L_{y}, L_{y}]$ ⁠, with periodic boundary condition in both x and y directions, which for our model (1) would imply

\begin{matrix} q (x, - L_{y}) = q (x, L_{y}), q (- L_{x}, y) = q (L_{x}, y), \\ n (x, - L_{y}) = n (x, L_{y}), n (- L_{x}, y) = n (L_{x}, y) . \end{matrix}

(3)

The above doubly periodic boundary conditions offer a convenient setup to apply the standard Fourier decomposition to all quantities of interest, which can be computed efficiently with the Fast Fourier Transform (FFT).

However, the emergence of large-scale structures from the nonlinear drift wave dynamics, sometimes comparable to the size of the computational domain, raises questions regarding the role of boundary conditions and the applicability of periodic boundary conditions to this problem, particularly in the “radial” x direction. In this article, we therefore also consider a channel domain geometry,³⁴ with periodic boundary conditions in the “poloidal” y direction and zero radial flow and density fluctuation imposed along the radial edges $x = 0, L_{x}$ of the domain,

\begin{matrix} q (x, - L_{y}) = q (x, L_{y}), u (0, y) = u (L_{x}, y) = 0, \\ n (x, - L_{y}) = n (x, L_{y}), n (0, y) = n (L_{x}, y) = 0. \end{matrix}

(4)

We do not have a rigorous justification derived from first principles for these boundary conditions, whose primary merit is to be an alternative to the usual periodic boundary conditions. Still, these boundary conditions are reasonable in the sense that they would be physically appropriate for a plasma in contact with a solid wall at $x = L_{x}$ ⁠, and such that a strong internal transport barrier is present at x = 0, preventing transport across this flux surface. We also note that these boundary conditions for the density perturbation have been imposed in other reduced fluid models for the study of bursty turbulent transport with avalanches.²⁶

B. Stability analysis for the growth rate in drift wave turbulence

Linear stability analysis of the BHW model (1) highlights the driving mechanism for the generation of non-zonal fluctuations in the starting transient state. We write the solutions $(\tilde{φ}, \tilde{n})$ to the linearized BHW equations as normal modes,

\tilde{φ} = \hat{φ} exp (i (k \cdot x - ω t)), \tilde{n} = \hat{n} exp (i (k \cdot x - ω t)),

where $ω \equiv ω (k)$ is the wave frequency for the corresponding wave vector, and the subscript k for ${\hat{φ}}_{k}$ and ${\hat{n}}_{k}$ has been omitted for the simplicity of the notation. Non-zonal solutions are such that $k_{y} \neq 0$ ⁠.

We first consider the non-dissipative case, with μ = 0. Plugging the normal mode solutions in the linearized BHW equations, we find the following solutions to the dispersion equation:^22,39

ω^{\pm} = \frac{α}{2} \frac{1 + k^{2}}{k^{2}} [\pm R^{\frac{1}{4}} cos \frac{θ}{2} - i (1 \mp R^{\frac{1}{4}} sin \frac{θ}{2})],

(5)

with

\begin{matrix} R (k; \frac{κ}{α}) = 1 + 16 γ^{2}, γ (k; \frac{κ}{α}) = \frac{κ}{α} \frac{k_{y} k^{2}}{{(1 + k^{2})}^{2}} = \frac{ω_{*} k^{2}}{α (1 + k^{2})}, \\ ω_{*} = \frac{κ k_{y}}{1 + k^{2}}, \end{matrix}

where $ω_{*}$ is called the drift wave frequency. From this expression, we can see that the question of whether a mode is stable or not only depends on the physical parameters κ and α through the ratio $κ / α$ ⁠. For $α \neq 0, κ \neq 0$ ⁠, all non-zonal modes have a solution with a positive imaginary part, so they are all unstable to this drift-wave instability. In the limit $α \to + \infty$ ⁠, the imaginary part of the frequency goes to zero for all modes, which means that adiabatic electrons stabilize the instability. In other words, the instability is driven by the combination of the background density gradient and the electron resistivity parallel to the magnetic field lines. The BHW model is thus a model for the study of electrostatic resistive drift wave turbulence, just like the MHW model.²⁰

Let us now turn to the more interesting and physically relevant dissipative case, with $μ > 0$ ⁠. The homogeneous damping operator with strength μ acts as a stabilizing effect for the system, acting most strongly on the small-scale modes. The damping effect can be easily accounted for by considering solutions to the linearized equations with dissipation of the form,

\begin{matrix} {\tilde{φ}}^{\pm} = \hat{φ} e^{- μ k^{2} t} e^{- i ω^{\pm} t} e^{i k \cdot x}, \\ {\tilde{n}}^{\pm} = \hat{n} e^{- μ k^{2} t} e^{- i ω^{\pm} t} e^{i k \cdot x} . \end{matrix}

The cross field collisional diffusion of the electrons acts as a counterbalance to the drift resistive instability. A mode is stable when this damping effect is larger than the linear growth rate. If we write the linear growth rate as $σ^{+} = Im ω^{+}$ ⁠, the criterion can be written as

μ \geq \frac{σ^{+}}{k^{2}} \Leftrightarrow μ \geq \frac{α}{2} \frac{1 + k^{2}}{k^{4}} [- 1 + \frac{1}{\sqrt{2}} {(1 + \sqrt{1 + 16 \frac{κ^{2}}{α^{2}} \frac{k_{y}^{2} k^{4}}{{(1 + k^{2})}^{4}}})}^{1 / 2}] .

(6)

From the criterion (6), we can conclude that in the presence of dissipation in the BHW model, short wavelength modes will be stable, while large wavelength modes will be unstable to the drift resistive instability. This instability corresponds to the first excited drift wave state. In the subsequent, nonlinear phase, a secondary instability³³ takes place through the nonlinear transfer of energy from the transient fluctuating modes to the zonal directions, and dominant zonal jets are created, as we will see in our numerical results.

Note that in the limit $\frac{κ}{α} ≪ 1$ ⁠, a criterion can be derived for the overall linear stability of the system, by Taylor expanding the right-hand side of (6) and retaining only the term which is first order in $\frac{κ^{2}}{α^{2}}$ ⁠,

μ \geq \frac{κ^{2}}{α} max_{k_{x}, k_{y}} \frac{k_{y}^{2}}{{(k^{2} + 1)}^{3}} = \frac{4}{27} \frac{κ^{2}}{α} .

(7)

C. Conserved quantities and their dynamical equations

The nonlinear interactions conserve two important quantities in the BHW model (1). They are the potential enstrophy W and total energy E defined as

W = \frac{1}{2} \int_{D} q^{2} d xdy = \frac{1}{2} \int_{D} {(\nabla^{2} φ - \tilde{n})}^{2} d xdy, E = \frac{1}{2} \int_{D} ({| \nabla φ |}^{2} + n^{2}) d xdy .

(8)

First, the dynamical equation for the total enstrophy W is determined by the potential vorticity equation

\frac{d W}{d t} = κ \int_{D} \tilde{u} \tilde{n} d xdy - μ \int_{D} {| \nabla q |}^{2} d xdy .

(9)

We thus see that the total particle flux,

Γ = \int_{D} \tilde{u} \tilde{n} d xdy,

(10)

acts as a forcing effect exciting drift waves. Second, the dynamical equation for the total energy E can be derived in a similar way based on the vorticity and density equations,

\frac{d E}{d t} = \int_{D} (κ + \bar{v}) (\tilde{u} \tilde{n}) d xdy - α \int_{D} {(\tilde{n} - \tilde{φ})}^{2} d xdy - μ \int_{D} ({| Δ φ |}^{2} + {| \nabla n |}^{2}) d xdy .

(11)

We emphasize the presence of the term $\int \bar{v} (\tilde{u} \tilde{n}) d xdy$ ⁠, which we refer to as the advected particle flux, in this energy equation for the BHW model. This term is absent in the energy equation for the MHW model and represents the only difference between the energy equations in the BHW model and the MHW model.²² It can be understood as the transport of the radial particle flux $\tilde{u} \tilde{n}$ by the poloidal mean velocity $\bar{v} = \partial_{x} \bar{φ}$ ⁠. Using the selective decay principle, we showed in Refs. 40 and 41 that this term plays a major role in explaining the persistent zonal flow $\bar{v}$ over a wide parameter range in the BHW model, blocking direct energy cascade all the way to the smallest scale. Besides, we observe that the negative-definite term $α \int_{D} {(\tilde{n} - \tilde{φ})}^{2} d xdy$ appears in the equation for the evolution of the energy E as an energy sink due to the electron collisional parallel resistivity, while this term is absent from the equation for the evolution of the enstrophy W.

III. DIMITS SHIFT AND STATISTICAL TRANSITION IN THE CHANNEL GEOMETRY

The Dimits shift usually refers to the nonlinear upshift of the critical temperature gradient for the rapid transition from a low transport regime to a strong turbulent transport regime, which is observed in gyrokinetic simulations of ion-temperature-gradient turbulence.^8,9,42 In our BHW model, which does neither have a temperature gradient nor a heat flux, but does have a density gradient and particle transport, an analogous Dimits shift together with a sharp transport regime transition is observed for the total particle flux. In this section, we characterize this transport transition and the corresponding properties of the magnetized plasma turbulence in our channel geometry. These results will be compared with the Dimits shift we obtain for the doubly periodic boundary conditions in Sec. V.

A. Numerical treatment of the boundary conditions

For our numerical simulations of the BHW model, we implemented an efficient pseudo-spectral scheme.²³ For the doubly periodic computational domain (3) with $x = (x, y) \in [- L_{x}, L_{x}] \times [- L_{y}, L_{y}]$ ⁠, the state variables are expanded in Fourier modes according to

φ = \sum_{k} {\hat{φ}}_{k} e^{i k \cdot x}, n = \sum_{k} {\hat{n}}_{k} e^{i k \cdot x}, ζ = \nabla^{2} φ = - \sum_{k} k^{2} {\hat{φ}}_{k} e^{i k \cdot x},

(12)

with spectral wave vector $k = (k_{x}, k_{y}) = (\frac{π}{L_{x}} n_{x}, \frac{π}{L_{y}} n_{y})$ and $k^{2} = {| k |}^{2}$ ⁠, and integer indices $n_{x} = - \frac{N_{x}}{2} + 1, \dots, \frac{N_{x}}{2}$ ⁠, and $n_{y} = - \frac{N_{y}}{2} + 1, \dots, \frac{N_{y}}{2}$ truncated at a large wavenumber. Zonal modes are represented by wavenumbers k_x and k_y = 0, and non-zonal fluctuation modes require $k_{y} \neq 0$ ⁠. Observe that larger domain sizes L_x, L_y will allow for more intermediate scales in the corresponding directions in the spectral domain.

For the channel domain geometry (4), we define the flow solutions on the computational domain with $x = (x, y) \in [0, L_{x}] \times [- L_{y}, L_{y}]$ ⁠. The boundary conditions are periodic in the y direction, and such that there are no density fluctuations and no x-directed flow at the boundaries $x = 0, L_{x}$ of the domain. The latter is satisfied by imposing the following conditions, for any $y \in [- L_{y}, L_{y}]$ ⁠:

\begin{matrix} φ (0, y, t) = φ (L_{x}, y, t) = 0, \\ n (0, y, t) = n (L_{x}, y, t) = 0. \end{matrix}

(13)

Note that strictly speaking, the left and right edges only need to be lines with constant potential (not necessarily zero) for the channel boundary conditions to be satisfied. However, we can set the potential to zero at both edges without loss of generality because the potential is always defined to within a free constant, which we fix in such a way that $φ = 0$ at the left edge, and because the BHW model is Galilean invariant,²² which means that our results will be insensitive to the addition of a large scale profile of the form Vx to the potential profile.

In order to have an efficient pseudo-spectral scheme for the numerical simulations for the channel domain, we use the numerical trick already implemented in Ref. 34, namely, we extend the original state variables f to an odd periodic extension $f^{odd}$ in the x direction:

f^{odd} (x, y) = {\begin{matrix} f (x, y), & 0 \leq x \leq L_{x}, \\ - f (- x, y), & - L_{x} \leq x < 0, \end{matrix}

with $f^{odd} (- L_{x}, y) = f^{odd} (0, y) = f^{odd} (L_{x}, y) \equiv 0$ as required for an odd periodic function. This extended function $f^{odd}$ is doubly periodic on the extended domain $[- L_{x}, L_{x}] \times [- L_{y}, L_{y}]$ and naturally satisfies the channel boundary condition (13) on the constrained half domain $[0, L_{x}] \times [- L_{y}, L_{y}]$ ⁠. In practice, $f^{odd}$ is constructed as a Fourier expansion in the y direction and sine expansion in the x direction,

f^{odd} (x, y) = \sum_{k_{x}} i {\hat{a}}_{k_{x}} (y) e^{i k_{x} x}, {\hat{a}}_{k_{x}} (y) = \sum_{k_{y}} {\hat{f}}_{k_{x}, k_{y}} e^{i k_{y} y},

with ${\hat{a}}_{k_{x}}$ real coefficients from the sine transform. With this approach, the simulations for the channel geometry can be computed using the same pseudo-spectral code as the code for the doubly periodic simulations on the domain $[- L_{x}, L_{x}] \times [- L_{y}, L_{y}]$ ⁠, with the only modification that one must take the odd decomposition of the functions.

For our simulations, the initial state is taken to be Gaussian random fields for the vorticity and the density, with zero mean and very small variance. The basic domain size is $L_{x} = L_{y} = 20$ along the two directions, and an elongated zonal length case L_x = 100 will be compared, which shows richer dynamics. We change the values of the model parameters κ and α to investigate the transitions in the turbulent flow. The collisional diffusion parameter is kept relatively small, either $μ = 1 \times 10^{- 3}$ or $μ = 1 \times 10^{- 4}$ depending on the cases, and diffusionless cases with μ = 0 will also be considered. In addition, hyperviscosities, $- ν Δ^{s} q$ and $- ν Δ^{s} n$ ⁠, are added to each equation, respectively, in our numerical implementations, with small coefficient $ν = 7 \times 10^{- 21}$ and high order s = 8, in order to dissipate the energy in the under-resolved smallest scales. All numerical simulations are based on a standard pseudo-spectral scheme with N = 256 discretization points along each direction. A fourth-order explicit-implicit Runge–Kutta scheme is used for the time integration with the implicit part only for the stiff hyperviscosity term.

B. Sharp transport transition from zonal jet dominated regime to drift wave turbulence dominated regime

In the framework of the BHW model, a Dimits shift is a nonlinear upshift of the threshold for strong particle flux. In other words, it is the presence of a region in parameter space in which one would expect moderate to large particle flux according to linear drift wave analysis, but near-zero flux is in fact observed, because zonal flows quench the drift instability and prevent the radial transport of particles. We now investigate this phenomenon in detail.

In the BHW model (1), the parameters which control the linear drift wave instability to excite small-scale non-zonal fluctuation modes are the adiabaticity parameter α and the density gradient κ. The collisional diffusion parameter μ, on the other hand, acts on both the vorticity and density as a homogeneous damping effect, which is stabilizing. From the linear stability analysis of the BHW model in Sec. II B as well as in Refs. 23 and 39, we know that for given values of κ and μ, there is a critical α value, which we call $α_{cr} (κ, μ)$ ⁠, above which the system is linearly stable. Likewise, for given values of α and μ, there is a critical κ value, which we call $κ_{cr} (α, μ)$ ⁠, below which the system is linearly stable.

In our numerical simulations, for the linearly stable regime, with all unstable modes damped by the dissipation, the flow is regular with no turbulent structures. When one enters the linearly unstable regime, obtained for $α < α_{cr}$ or $κ > κ_{cr}$ ⁠, we see richer dynamical structures induced by the drift wave instability. Still, in a large region above the threshold for the linear instability, the excited drift waves observed in the transient state are completely eliminated through the nonlinear coupling between zonal and fluctuation modes. This nonlinear interaction mechanism is described in detail in Ref. 33 using the secondary instability analysis. In this intermediate regime, we therefore see the generation of dominant, long-lived zonal flows, and near-zero non-zonal turbulent fluctuations, as we will show explicitly in Sec. III C. As a result, there is a near-zero total particle flux $Γ = \int_{D} \tilde{u} \tilde{n} d xdy$ ⁠, mostly determined by random numerical fluctuations. As one further decreases α or increases κ, to second critical parameter values we call $α_{turb} (κ, μ)$ or $κ_{turb} (α, μ)$ ⁠, the secondary instability saturates, and we see the generation of strong turbulent flows which are not quenched by zonal flows; this is the starting point for strong zonal transport. Before these critical values are reached, in the regime $α_{turb} < α < α_{cr}$ or $κ_{cr} < κ < κ_{turb}$ ⁠, the flow converges to regular zonal structures with little turbulence remaining after the starting transient state. The Dimits shift in the BHW model can be represented by the distance $α_{cr} - α_{turb}$ in the adiabaticity parameter, or $κ_{turb} - κ_{cr}$ in the density gradient parameter. The existence of the Dimits shift and a sharp transition in the level of particle transport for $α = α_{turb} < α_{cr}$ or $κ = κ_{turb} > κ_{cr}$ in the BHW model can be clearly seen in Fig. 1, which shows the change in the total particle flux as the parameter κ changes in the range $[0.1, 5]$ with fixed $α = 0.5$ ⁠, and as the parameter α changes in the range $[0.05, 5]$ with fixed κ = 1, for two different values of the collisional diffusion parameter μ: $μ = 10^{- 4}$ and $μ = 10^{- 3}$ ⁠. We observe that, as expected, a larger value of μ pushes the transition point to turbulence to a larger value of κ or a smaller value of α, requiring a stronger instability to induce strong transport. Still the sudden phase transition is seen for both the smaller and the larger value of μ. Finally note that for α and μ fixed, the particle flux continually increases as one increases κ into the drift wave turbulence dominated regime. On the other hand, for κ and μ fixed, the particle flux seems to asymptote as one decreases α toward zero, and would in fact decrease if we plotted the particle flux for even smaller values of α.²³ This is because the growth rate of the linear instability always increases with κ, whereas that growth rate decreases and converges to zero in the limit $α \to 0$ ⁠. The parameter values for α, κ, and μ are picked according to the reference values in Refs. 19 and 28 and the previous simulation tests in Refs. 23 and 34 where the flows display the most interesting structures and statistics. From our numerical simulations, we computed the equilibrium total particle flux in the final statistical steady state by taking the time average of 5000 time steps, which averages out the random, sometimes sharp fluctuations in the time series (see Fig. 3 for the time evolutions of the total fluxes).

FIG. 1.

View large Download slide

Total particle flux $Γ = \int_{D} \tilde{u} \tilde{n} d xdy$ as a function of the parameters α and κ from direct numerical solutions of the BHW model in the channel domain geometry. The curves corresponding to two different values of the collisional diffusion parameter μ, $μ = 1 \times 10^{- 3}$ ⁠, and $μ = 1 \times 10^{- 4}$ ⁠, are compared. The critical values of $κ_{cr}$ and $α_{cr}$ for the onset of the linear instability are listed near the corresponding lines. Notice the logarithmic scale for the y-axis, covering several orders of magnitude.

Table I gives the values of the total particle flux as the parameters κ or α are varied to go from the zonal flow dominated regime to the drift wave turbulence dominated regime, for $μ = 0, 10^{- 4}, 10^{- 3}$ ⁠. We did not fill all boxes in the table because the simulations are computationally expensive, and we decided that the data corresponding to the boxes left blank would not give any more insights into the dynamics. The information in the table complements Fig. 1, from which actual numbers may be hard to extract. We find the following:

For the relatively large value of the diffusion parameter μ, $μ = 1 \times 10^{- 3}$ ⁠, the critical values for the onset of the linear instability are $κ_{cr} = 0.058$ and $α_{cr} = 1.475 \times 10^{2}$ ⁠, while the strong turbulent transport is triggered at the much later stage at around $κ_{turb} \sim 1$ and $α_{turb} \sim 0.5$ ⁠.
For the smaller value of the collisional diffusion parameter, $μ = 1 \times 10^{- 4}$ ⁠, the critical values for the onset of the linear instability are $κ_{cr} = 0.0185$ and $α_{cr} = 1.47 \times 10^{3}$ ⁠, while strong turbulent transport is found for $κ_{turb} \sim 0.4$ and $α_{turb} \sim 1.5$ ⁠.
For the diffusionless case μ = 0, the system is linearly unstable for all values of $α \neq 0, κ \neq 0$ ⁠, but a transition is still observed for the measured level of particle transport, although the transition can be less sharp than for the $μ \neq 0$ cases depending on the parameter values and on the choice of parameter which is varied.

TABLE I.

Total particle flux $Γ = \int_{D} \tilde{u} \tilde{n} d xdy$ with changing values of α and κ. Cases corresponding to different values of the collisional diffusion parameter μ, $μ = 1 \times 10^{- 3}$ ⁠, $μ = 1 \times 10^{- 4}$ ⁠, as well as the diffusionless case μ = 0 are compared.

$α = 0.5$ ⁠, κ	0.1	0.3	0.5	0.8	0.9	1	2
μ = 0	0.0031	0.0155	1.3688	5.1319	…	6.4705	…
$μ = 1 \times 10^{- 4}$	0.0066	0.0316	1.5238	3.9184	…	6.8305	53.8872
$μ = 1 \times 10^{- 3}$	0.0064	…	0.0448	0.1143	0.2320	3.6692	54.4245
κ = 1, α	0.05	0.1	0.6	1	2	3	5
μ = 0	…	11.2053	5.7420	4.0262	1.5184	0.9160	0.2707
$μ = 1 \times 10^{- 4}$	10.0764	10.7077	…	2.1394	0.0093	0.0081	…
$μ = 1 \times 10^{- 3}$	17.1608	16.4841	0.1002	0.0459	0.0484	…	…

$α = 0.5$ ⁠, κ	0.1	0.3	0.5	0.8	0.9	1	2
μ = 0	0.0031	0.0155	1.3688	5.1319	…	6.4705	…
$μ = 1 \times 10^{- 4}$	0.0066	0.0316	1.5238	3.9184	…	6.8305	53.8872
$μ = 1 \times 10^{- 3}$	0.0064	…	0.0448	0.1143	0.2320	3.6692	54.4245
κ = 1, α	0.05	0.1	0.6	1	2	3	5
μ = 0	…	11.2053	5.7420	4.0262	1.5184	0.9160	0.2707
$μ = 1 \times 10^{- 4}$	10.0764	10.7077	…	2.1394	0.0093	0.0081	…
$μ = 1 \times 10^{- 3}$	17.1608	16.4841	0.1002	0.0459	0.0484	…	…

We conclude in this section that the sudden transition in the level of transport characterizing the Dimits shift is successfully captured in the BHW model for a channel domain geometry. As we will show in Sec. V, this is also the case for the more commonly considered doubly periodic geometry. This is a key observation of this article since it is a unique feature of the BHW model among the HW models. Indeed, we will also show in Sec. V that the MHW model does not produce similarly sharp transitions as one decreases α or increases κ. This demonstrates that in reduced fluid models, the Dimits shift sensitively depends on the treatment of the electron dynamics parallel to the magnetic field lines, an observation which was also recently made with the Terry–Horton model.²⁴

C. Detailed analysis of the transport transition

We now take a more detailed look at the dynamical structures in the flow before and after the transition from a low transport regime to a high transport regime. Specifically, we compare the flow solutions and statistics for the different regimes. Figure 2 first plots, for $μ = 10^{- 4}$ ⁠, the time-series for the energy in the non-zonal fluctuation modes with $k_{y} \neq 0$ and the total zonal energy in the modes with k_y = 0 in the zonal flow dominated regime obtained for $κ = 1, α = 2$ ⁠, which is just before the transition, and in the drift wave turbulence dominated regime, obtained for $κ = 1, α = 1$ ⁠, which is shortly after the transition. In both cases, starting from an initial state with small amplitude fluctuations and with little energy, the energy in the fluctuations first increases dramatically due to the drift instability of non-zonal modes. In the zonal flow regime with α = 2, the energy in non-zonal modes is then gradually transferred to the zonal state in the next stage through the secondary energy transfer. The total energy in the zonal modes reaches approximately the same level as the maximum excited fluctuation energy in the starting transient state before the dissipation effect takes over, implying full conversion of the fluctuation energy to the zonal state. In contrast, in the drift wave turbulence dominated regime, observed for α = 1, the non-zonal fluctuations are never completely quenched by the formation of zonal jets. The non-zonal energy is still transferred to drive the zonal modes through the nonlinear coupling, but the instability intermittently injects energy as bursts into the fluctuation modes. In this way, the turbulent fluctuations are maintained and cannot be completely converted to the zonal state. We note that the zonal energy always increases after the burst in the fluctuation energy excited by the linear drift instability, from which we can infer the direction of the nonlinear energy transfer from non-zonal to zonal modes.

FIG. 2.

View large Download slide

Development of the zonal state from initial fluctuations with small amplitude for the case $μ = 10^{- 4}$ ⁠. Two typical parameter regimes with dominant zonal flow (left, $κ = 1, α = 2$ ⁠) and turbulent flow (right, $κ = 1, α = 1$ ⁠) are compared. The upper row compares the time series for the energy in the non-zonal modes $k_{y} \neq 0$ and in the zonal modes with k_y = 0. The lower row shows the development of the zonal mean flow profiles at several times as the zonal flows are generated from the fluctuations.

The emergence of the zonal state is also illustrated in the second row of Fig. 2. We plot the development of the zonal velocity profiles $\bar{v} = \partial_{x} \bar{φ}$ measured at several different time instants. For the simulations with α = 2, corresponding to zonally dominated dynamics, two jets first appear, before the zonal state gradually evolves to a one jet structure. In contrast, a dominant single jet is directly formed in the simulations with α = 1, which correspond to the turbulent regime, with the development of a stronger positively directed flow in the center of the domain and a stronger negatively directed flow near the two edges of the computational domain. This observation is consistent with the selective decay results in Ref. 41 showing a more stable solution with smaller wavenumber, leading to a single dominant large-scale jet.

Next, we illustrate the representative flow field snapshots throughout the flow transition and the corresponding transition for the time series of the enstrophy and the fluxes by comparing the direct simulation results for several different parameter regimes, corresponding to the full range from pure zonal jets to fully turbulent flows. We choose κ = 1 and $μ = 10^{- 3}$ ⁠, which give typical examples of the flow transition, and vary α in order to obtain the desired transition. Snapshots of the ion vorticity and time series of the energetics are shown in Fig. 3. We first confirm with the snapshots that, as the value of α increases, the turbulent flow gradually transitions to regular zonal jets. We can then look at the corresponding effects on the time series for the total enstrophy W, the total particle flux Γ, and the advected particle flux $\int_{D} \bar{v} (\bar{\tilde{u} \tilde{n}}) d xdy$ ⁠. We recall here that the flux terms act as energy source and sink in the equations for the enstrophy (9) and for the energy (11). In the turbulent regime, the enstrophy is excited by the bursty injection from the particle flux. As the value of α decreases, the frequency of intermittent bursts increases, and the solution is kept in a state with high enstrophy for both zonal modes and non-zonal modes. We will return to these features when we discuss avalanche-like bursts in more detail in Sec. IV. For large values of α, the bursts in the fluxes are much less frequent and have much smaller amplitude. The total enstrophy converges to a constant value with most of the enstrophy in the zonal modes: this is the zonal flow dominated regime.

FIG. 3.

View large Download slide

Snapshots of the ion vorticity $ζ = \nabla^{2} φ$ (upper row) and time series of the total enstrophy $W = \int_{D} q^{2} d xdy$ ⁠, as well as the total fluxes $Γ = \int_{D} \tilde{u} \tilde{n} d xdy$ and $\int_{D} \bar{v} \bar{\tilde{u} \tilde{n}} d xdy$ (lower row) for the channel domain geometry, for different values of α and fixed $κ = 1, μ = 1 \times 10^{- 3}$ ⁠. More frequent turbulent intermittent bursts are observed as the flow enters the more turbulent regimes.

Valuable insights can also be obtained by comparing, for values of α corresponding to the different transport regimes, the spectra for the total energy and for the enstrophy, as well as the zonal mean profiles at the end of the simulations, after statistical steady-state has been reached. This is precisely what we display in Fig. 4, separating non-zonal fluctuation modes $k_{y} \neq 0$ and the zonal modes k_y = 0 for the spectra. κ is set to 1, μ is set to $10^{- 3}$ ⁠, and the values of the adiabaticity parameter α are chosen in the interval $[0.05, 1]$ to capture both the turbulent dominated regime and the zonal flow dominated regime. The two regimes can be clearly distinguished from the spectra, with energetic small scales for $α < 0.6$ ⁠, and a sudden drop of both the energy and enstrophy in the small scale modes for $α \geq 0.6$ ⁠. Thus, in the zonal flow dominated regime with relatively weak lineary instability, the excited energy and enstrophy are kept in the large-scale zonal modes and the downscale cascade is blocked. When the linear instability is stronger, turbulent transport is induced by the stronger downscale cascade to create small-scale fluctuations. This characterizes the statistical transition in the turbulent transport.

FIG. 4.

View large Download slide

Spectra for the non-zonal fluctuation modes (⁠ $k_{y} \neq 0$ ⁠, left) and the zonal modes (k_y = 0, right) for different values of α, and κ = 1, and $μ = 10^{- 3}$ ⁠. The kinetic energy $k^{2} | {\hat{φ}}_{k} |^{2} + | {\hat{n}}_{k} |^{2}$ ⁠, and the enstrophy $| {\hat{q}}_{k} |^{2}$ at each scale are compared. A clear transition can be observed in both the zonal and fluctuation spectra depending on the parameter values of α. The bottom rows of the figure compare the zonal structures in the velocity and the density for different values of α. The curves corresponding to the turbulence dominated regime are plotted in solid lines, and the curves corresponding to the zonal flow dominanted regime are plotted in dashed lines.

The two bottom rows of Fig. 4 compare the zonal velocity and zonal density profiles measured at the end of the simulation, in statistical steady-state, for $α = 0.1, 0.3, 0.5, 0.8, 1$ ⁠. In all cases, a strong flow in the positive y direction and a strong density gradient are seen at the center of the computational domain. The zonal profiles before and after the transport transition, however, display different features in the vicinity of the edges of the computational domain. In the turbulent regime with strong transport, stronger zonal flow $\bar{v}$ ⁠, a larger flow shear and a larger density gradient are observed near the boundaries. These features, which were already found and analyzed in our previous study³⁴ of the BHW model in the channel geometry, are not present for doubly periodic boundary conditions and are critical to explain the differences between the channel geometry and the doubly periodic geometry we will find for the particle transport in Secs. IV and V.

IV. AVALANCHE-LIKE STRUCTURES IN DOMAINS WITH DIFFERENT ASPECT RATIOS

Besides the Dimits shift, other salient aspects of drift wave turbulence driven transport which remain incompletely understood are the importance of non-diffusive, non-local transport mechanisms in the observed level of transport, and whether meso-scale ballistically propagating structures observed for various physical quantities such as the heat flux in ITG simulations^10–13,16 are the main non-local transport mechanism. An emerging hypothesis, mostly based on the results of numerical simulations, is that these ballistically propagating structures may be sorted into two categories:¹⁵ avalanches,^10–13 whose signatures are large amplitude, rather disordered radially propagating bursts, and more coherent solitary structures^{12,14–16,27} for which the term ferdinons was recently coined.¹⁶ It is not yet clear how avalanches and soliton-like solutions are related, or whether all soliton-like solutions are ferdinons.¹⁶ In this article, we will not propose a theoretical explanation for these phenomena. However, we will now show that both disordered avalanches and solitary structures are commonly observed for the particle flux in the BHW model, which is therefore a remarkably simple model for their study. Furthermore, our simulations suggest that avalanches and solitary structures are indeed intimately related, in the sense that the coherent solitary structures frequently emerge from the avalanches in the edge of the computational domain, where the velocity shear and density gradients are large. Finally, by comparing the BHW model and the MHW model, we find that the soliton-like solutions do not emerge from the avalanches in the MHW model, which indicates that electron dynamics parallel to the magnetic field plays a fundamental role in this mechanism.

A. Avalanche-like structures in the radial particle flux

We focus here on the same channel domain geometry and the same parameter regime as in Sec. III. First, we show the time evolution of the radial quantities (that is, the zonally y-averaged quantities, varying only along the radial x-direction) to illustrate the representative avalanche-like bursts emerging in the BHW model in the regime past the Dimits threshold, where significant levels of transport are measured. In the first row of Fig. 5, the time evolution of the zonal flow $\bar{v}$ displays a dominant large-scale single jet structure frequently shifting in time. The second row of Fig. 5 shows the time evolution of the radial particle flux transport $\bar{\tilde{u} \tilde{n}}$ ⁠. We observe two major features: (1) a central region, whose extent depends on the value of α, in which we see more incoherent avalanche bursts, very similar of the avalanches reported for the heat flux in gyrokinetic simulations^11,12 and (2) in the vicinity of the edges of the computational domain, coherent radially propagating structures with constant speed. The soliton-like coherent structures near the edges, which may be ferdinons,¹⁶ are traveling at a near-constant speed. In our normalized units, the solitons travel at an approximate speed of 0.072 for $α = 0.1, κ = 1$ ⁠, and a faster speed of 0.125 for $α = 0.5, κ = 1$ ⁠. In a future study, it will be interesting to conduct a detailed analysis of the traveling speed as a function of the key ratio $α / κ$ ⁠. The shifts of the zonal structure and the associated spikes in the zonal flow are strongly correlated with the emergence of the solitary structures from the avalanches. For the smaller α value, $α = 0.1$ ⁠, which corresponds to the more turbulent regime, the heat flux avalanches are stronger in amplitude and have more interacting structures, and the solitary structures emerge more frequently. For the larger α value, $α = 0.5$ ⁠, the heat flux avalanches are weaker, and solitary structures are less frequent. Finally, we note that the solitary structures are most extended for the larger value of α, which has smaller velocity and density gradients near the boundary of the domain.

FIG. 5.

View large Download slide

Time evolution of the zonal mean velocity $\bar{v}$ and the radial particle flux $\bar{\tilde{u} \tilde{n}}$ from BHW model simulations (top and middle rows). Two representative regimes with relatively sparse bursts of zonal transport (⁠ $α = 0.5$ ⁠, left) and more frequent and stronger zonal bursty transport (⁠ $α = 0.1$ ⁠, right) are shown. The PDFs of the radial particle flux $\bar{\tilde{u} \tilde{n}}$ (left) and power spectra from the time series of the radial particle flux (lower right) for different values of the parameter α are compared in the bottom row. As before, solutions corresponding to the turbulence dominated regime are plotted with solid lines, while solutions corresponding to the zonal flow dominated regime are plotted with dashed lines. For all these results, κ = 1 and $μ = 10^{- 3}$ ⁠.

We next turn to the statistical and spectral signatures of this bursty meso-scale transport. The bottom row of Fig. 5 plots the empirical probability density functions (PDF) and the power spectra sampled from time-series of the zonal particle transport $\bar{\tilde{u} \tilde{n}}$ for different values of α. The empirical PDFs are constructed by computing the normalized histograms of the time series of the radial particle flux after time t = 1000. For most values of α, statistical steady-state is established at that point. For the largest values of α, some transient events still occur infrequently at the beginning of that time window, which however do not significantly alter the statistics, and are not seen later on. All the PDFs show large positive skewness, implying a zonal particle transport toward the outer boundary, as expected. Large value events in the radial particle flux are much more common than they would be for a normal distribution with the same mean and variance. For the large transport regime, for $α = 0.1$ ⁠, the PDF has a larger variance, and the power spectrum has larger values at high frequencies. As α increases, the PDFs have more of their mass in the vicinity of zero, and the power spectra decay at a faster rate as the frequency increases, implying weaker zonal transport and fewer occurrences of extreme events. We note that larger values of α have fatter tails, which are due to the existence of more infrequent but larger amplitude, more localized bursts, as seen in the time evolution of the particle flux. Below the Dimits threshold, in the zonal flow dominated regime with near-zero particle flux, these bursts do, however, not appear any longer, once statistical steady-state is reached.

Insights regarding the role of the parallel electron dynamics in the observed avalanches and solitary structures in the BHW model can be gained by plotting the corresponding simulation results for the MHW model in Fig. 6, for the same set of values α, κ, and μ. Zonal jets are also developed in the MHW model case. However, after the zonal mean flow $\bar{v}$ reaches a statistical steady state, the flow has much weaker time variability than in the BHW model. Correspondingly, the radial particle flux in the MHW model has the relatively incoherent avalanche structures present in the BHW in the center of the computational domain, but does not have the solitary structures emerging from the avalanches in the vicinity of the boundaries.

FIG. 6.

View large Download slide

Time evolution of the zonal mean velocity $\bar{v}$ and the radial particle flux $\bar{\tilde{u} \tilde{n}}$ from MHW model simulations (top and middle rows) for the same values of α considered in Fig. 5: $α = 0.5$ (left) and $α = 0.1$ (right). Time series of the total enstrophy, the enstrophy in the zonal modes (top figures of bottom row) and of the total particle flux Γ and of the total advected particle flux $\int_{D} \bar{v} (\bar{\tilde{u} \tilde{n}}) d xdy$ (bottom figures of bottom row) also for these values of α: $α = 0.5$ (left) and $α = 0.1$ (right). For all these results, κ = 1 and $μ = 10^{- 3}$ ⁠.

In the bottom half of Fig. 6, we also show the time series for the total enstrophy W as well as the fluxes Γ and $\int_{D} \bar{v} \bar{\tilde{u} \tilde{n}} d xdy$ from the MHW model, which can be compared with the corresponding BHW results in Fig. 3. We see that in the MHW model, the zonal contribution dominates largely for the enstrophy. Furthermore, for both values of α, the fluxes lack the intermittent bursts we had observed for the time series in the BHW model. In the turbulent regime, the total particle flux thus has a larger mean value in the MHW model, but has less variability and less intermittency. This result supports the following conjecture: the different treatment of the parallel electron dynamics in the BHW and MHW models is responsible for the presence or not of the radially propagating coherent solitary structures near the boundaries of the computational domain, and these solitary structures are responsible for the largest intermittent bursts measured for the fluxes in the BHW model, while the avalanches are not. In this article, we will not present a theory or additional arguments to demonstrate the validity of this hypothesis, but it is the subject of ongoing work, to be reported in a future article.

In order to be complete on this point, it is important to clarify that in a different setup, coherent radially propagating structures have been seen in the MHW model, as reported in Ref. 25. In that work, the radially propagating structures were obtained with a strong imposed background zonal shear flow. However, in our direct simulations of the MHW model, it appears that the drift instability cannot generate and support such a strong zonal structure by itself, leading to the absence of soliton-like transporting structures, in contrast to the BHW model.

B. The effect of different domain aspect ratios

In Secs. III and IV A, our computational domain was always such that $L_{x} = L_{y} = L$ ⁠. As a first investigation of the role of boundary conditions and of Fourier resolution on our results, we consider the effect of a different size of the computational domain on our observations. Specifically, we extend the channel width in the x direction to $[0, a L]$ ⁠, where we call $a = L_{x} / L_{y} > 1$ the aspect ratio. Our past results suggest that by doing so, we may observe the generation of more zonal jets, which may interact with each other, with multi-scale dynamics.²³

1. Multi-scale jet dynamics and large scale avalanche evolution

We first consider a computational domain with aspect ratio a = 5, for the sets of parameters $κ = 1, α = 0.5$ ⁠, corresponding to the strong turbulence regime, and $κ = 0.5, α = 0.5$ ⁠, corresponding to a weaker turbulence regime. The collisional diffusion coefficient is $μ = 10^{- 3}$ for both test cases. The snapshots of the ion vorticity and the time series of the zonal mean flow are compared for these two regimes in Fig. 7. Differences with our results for the computational domain with aspect ratio a = 1, shown in Fig. 3, appear clearly. For both values of α, larger numbers of zonal jets are created in the starting transient stage instead of a dominant single jet structure as before. The multiple jets then interact and merge into more dominant zonal structures. At the same time, the vorticity plots indicate that the flow is more turbulent, with more small-scale vortical structures. The multi-scale jet dynamics is correlated with large-scale dynamics for the particle flux avalanches, also shown in Fig. 7. The avalanches are particularly strong and widespread in the starting transient phase, before gradually shrinking to a narrower region in the center of the computational domain as the zonal jets merge to form stronger and more localized zonal structures. Furthermore, we observe that as before, simple traveling structures are seen emerging from the avalanches in the radial particle flux although they are not as coherent and localized as for the computational domain with aspect ratio a = 1.

FIG. 7.

View large Download slide

Snapshots of the ion vorticity (upper row), and time-series of the zonal mean velocity $\bar{v}$ and the radial particle flux $\bar{\tilde{u} \tilde{n}}$ (middle row) from the BHW model for a computational domain with aspect ratio a = 5. The figures on the left-hand side were obtained for $κ = 1, α = 0.5$ ⁠, and the figures on the right-hand side for $κ = 0.5, α = 0.5$ ⁠, with $μ = 10^{- 3}$ in both cases. Kinetic energy spectra (left), density energy $\frac{1}{2} \int_{D} n^{2} d xdy$ spectra (middle), and enstrophy spectra (right) for different domain aspect ratios a = 5, 4, 3, $κ = 1, α = 0.5$ ⁠, once statistical steady-state was reached (lower row).

We conclude this section with a comparison with simulations from the MHW model for the same values of the parameters α and κ, and the same computational domain. The results for the zonal mean velocity $\bar{v}$ and the radial particle flux $\bar{\tilde{u} \tilde{n}}$ are shown in Fig. 8. In the MHW simulations, we also see initial mergings of zonal jets, but these happen during a short transition period. After this transient, the zonal jets have much less variability than in the BHW model, they do not merge any more, and we do not observe any evidence of strong large scale dynamics and multi-jet interactions. This difference in behavior between the BHW model and the MHW model was already highlighted for doubly-periodic boundary conditions in Ref. 23. As a result of the weakness of the multi-scale interactions, there is little large scale dynamics in the particle flux avalanches as compared to the BHW model, and the avalanches remain spread across most of the computational domain. Finally, as for the simulations with aspect ratio a = 1, in the MHW model we do not see solitary structures emerging from the particle flux avalanches.

FIG. 8.

View large Download slide

Time series of the zonal mean velocity $\bar{v}$ and the radial particle flux $\bar{\tilde{u} \tilde{n}}$ from the MHW model for a computational domain with aspect ratio a = 5. The figures on the left and middle left were obtained for $κ = 1, α = 0.5$ ⁠, and the figures on the middle right and right were obtained for $κ = 0.5, α = 0.5$ ⁠.

2. Statistical invariance with different aspect ratios

In this final subsection, we show the statistical invariance of the simulations for different aspect ratios a of the computational domain once the transient multi-scale evolution we just discussed has taken place, and statistical steady-state has been reached. The bottom row of Fig. 7 displays the spectra of the kinetic energy $\frac{1}{2} \int_{D} | \nabla φ |^{2} d xdy$ ⁠, the density energy $\frac{1}{2} \int_{D} n^{2} d xdy$ and of the enstrophy once statistical equilibrium was reached, for the following aspect ratios of the computational domain: a = 3, 4, 5. A longer x domain length contains larger energy and also stronger particle flux in the transient state. Thus, it generates stronger turbulent dynamics. However in the final statistical equilibrium state, it is observed that the radially averaged steady-state energy and enstrophy spectra are approximately independent of the aspect ratio. This invariance confirms the existence of a consistent statistical structure in the channel domain flows even though much richer dynamics are created with a longer x domain length in the transient phase.

V. DIMITS SHIFT AND AVALANCHE STRUCTURES IN A DOUBLY PERIODIC GEOMETRY

Many studies of turbulent transport in magnetic confinement devices, either with reduced fluid models or the gyrokinetic equations, impose doubly periodic boundary conditions for their simulations. We now investigate how imposing doubly periodic boundary conditions instead of the channel geometry boundary conditions modifies our results. We have already studied the BHW model with doubly periodic boundary conditions extensively in Refs. 22 and 23. In this section, we therefore focus on the Dimits shift and signatures of non-local, non-diffusive transport, which are the central subjects of this article. For all the simulations, the computational domain is $[- L_{x}, L_{x}] \times [- L_{y}, L_{y}] = [- 20, 20] \times [- 20, 20]$ ⁠.

A. Dimits shift and comparison between the BHW and MHW models

We first compare the transition between the low transport regime and the high transport regime in the BHW model and the MHW model. A comprehensive summary of our results is given in Fig. 9, which shows, for each model, the total particle flux Γ as a function of κ for the fixed values $α = 0.5$ and $μ = 10^{- 3}$ ⁠, and as a function of α for the fixed values κ = 1 and $μ = 10^{- 3}$ ⁠, as well as the times series for the total energy and enstrophy and for the particle flux for $α = 0.5$ and κ = 1, and the spectra for the energy in variance $⟨ | f - ⟨ f ⟩ |^{2} ⟩$ ⁠, where we compute the second-order moment of the velocity and density fields $f = \nabla φ, n$ by taking the long time average $⟨ f ⟩ = \frac{1}{T} \int_{0}^{T} f (t) d t$ of the state variables for each scale with the same wavenumber k.

FIG. 9.

View large Download slide

Comparison of the total particle flux from the BHW and MHW model simulations (left column), and the time-series of the total energy and particle flux and the energy spectra for variance from the BHW and MHW model simulations with $α = 0.5$ and κ = 1 (right column). For all these simulations, $μ = 10^{- 3}$ ⁠.

In both models, a transition between a zonal flow dominated regime with low particle transport and a turbulent regime with strong particle transport is observed as we vary the physical parameters α and κ. However, the transition is much less sharp in the MHW model: for a fixed value of α, the particle flux increases with κ as a power law, and for a fixed value of κ, the particle flux increases with decreasing values of α as a power law as well. In other words, the MHW model does not have what one would characterize as a true Dimits shift, marked by a sudden transition in the particle flux,^8,9,24 as obtained in the BHW model, which has near-zero particle flux until the Dimits threshold is reached. This result can be understood in light of our previously published results,^22,23 and in particular, Fig. 4 in Ref. 22. In the MHW model, purely zonal states are never generated, even when the linear instability is weak. Small scales vortices are always present, which drive some amount of turbulent transport. They are responsible for the absence of a true nonlinear upshift for the threshold for strong particle transport in the MHW model. In contrast, when the linear instability is weak, pure zonal states are generated in the BHW model, which fully impede particle transport.

In addition, we identify further differences between the two models in the turbulence dominated regime, as illustrated by the time-series of the total energy and particle flux in the center column of Fig. 9 as well as the energy spectra are compared in detail in the right column of Fig. 9. The time series for the BHW model are characteristic of strong but intermittent turbulent transport, while the time series for the MHW model rapidly converge to a statistical steady state with much weaker fluctuations. In the turbulence dominated regime, transport in the MHW model is strong in averaged value but less bursty and variable than in the BHW model, as we had already found previously for the channel geometry. This is further confirmed by observing the variance energy spectra of the two models shown in the right column of Fig. 9. The variances in the leading modes (that is, the first two largest scales dominated by the variance in zonal modes $k_{x} = 1, 2$ ⁠) have much larger values in the BHW model than in the MHW model case. This is consistent with the much stronger variability in the zonal modes already reported in previous studies of the BHW model.^22,23,34

B. Transport structures for doubly periodic boundary conditions

We now look at the flow and transport features in the BHW model in more detail to compare them with the equivalent results we obtained in the channel geometry. In Fig. 10, we present the flow snapshots and time series of the total enstrophy and total particle flux for $κ = 1, μ = 10^{- 3}$ and three values of α, $α = 0.1, α = 0.5$ ⁠, and α = 1, which can be directly compared with Fig. 3. Just like in the channel domain case, and in agreement with our observations from Sec. V A, as the value of α increases, the flow goes through a transition from a fully turbulent regime with strong non-zonal vorticity and strong and bursty particle transport to a zonal flow dominated regime with nearly zero particle flux. We nevertheless note that for the same values of α, κ, and μ, the flow computed for doubly periodic boundary conditions is more turbulent, with more small scale vortices. This explains why, for fixed κ and μ, the Dimits threshold is found for a larger value of α for the doubly periodic boundary conditions than for the channel geometry.

FIG. 10.

View large Download slide

Snapshots of the ion vorticity field $ζ = \nabla^{2} φ$ (upper row) and time series of the total enstrophy $W = \int_{D} q^{2} d xdy$ ⁠, as well as the total fluxes $Γ = \int_{D} \tilde{u} \tilde{n} d xdy$ and $\int_{D} \bar{v} \bar{\tilde{u} \tilde{n}} d xdy$ (lower row) from the BHW model with doubly periodic boundary condition, for fixed $κ = 1, μ = 1 \times 10^{- 3}$ and different values of α: $α = 0.1$ (left), $α = 0.5$ (center), and α = 1 (right).

The time evolution of the zonally averaged velocity $\bar{v}$ and of the radial particle flux $\bar{\tilde{u} \tilde{n}}$ for the doubly periodic domain case are shown in Fig. 11, for the same sets of parameters as in the previous figure: κ = 1 and $α = 0.1, α = 0.5$ ⁠, and α = 1, $μ = 10^{- 3}$ ⁠. Figure 11 may be analyzed in light of the equivalent plots for the channel geometry presented in Fig. 5. With periodic boundary conditions in the x-direction as well, the zonal jets and the radial particle flux $\bar{\tilde{u} \tilde{n}}$ can cross the boundary and appear on the other side of the computational domain. The radial particle flux is thus observed to travel from both right to left and left to right across the edges of the domain, leading to avalanches throughout the computational domain. The zonal jets have the same variability as in the channel geometry, but we do not see the emergence of the soliton-like solutions which were localized in the regions of strong velocity shear and large density gradient. The strong turbulence regime with $α = 0.1$ displays strong, coherent, and frequent radially propagating fronts, with a very strong maximum particle flux ${max}_{t, x} \bar{\tilde{u} \tilde{n}} = 1895.2$ ⁠; in the weaker turbulence regime with $α = 0.5$ ⁠, the radially propagating fronts are still relatively strong, but appear less frequent and less coherent, and the maximum particle flux is reduced to 303.63; finally in the zonal flow dominant regime with α = 1, the particle flux is much weaker in amplitude, with maximum value 44.07, and much less coherent.

FIG. 11.

View large Download slide

Time evolution of the zonal velocity $\bar{v} = \partial_{x} \bar{φ}$ (left column) and the radial particle flux $\bar{\tilde{u} \tilde{n}}$ (right column) from the BHW model with doubly periodic boundary conditions, for κ = 1, $μ = 1 \times 10^{- 3}$ ⁠, and different values of the adiabaticity parameter α: $α = 0.1$ (top), $α = 0.5$ (middle), and α = 1 (bottom). The maximum values of radial particle flux in the three cases are 1895.2 in $α = 0.1$ ⁠; 303.63 in $α = 0.5$ ⁠; and 44.07 in α = 1.

VI. SUMMARY

We have conducted direct numerical simulations with the flux-balanced Hasegawa–Wakatani (BHW) model²² to study turbulent transport driven by a resistive drift-wave instability and mediated by the emergence of zonal flows. We have considered two different geometries, a channel domain geometry³⁴ and a doubly periodic geometry, and different aspect ratios of the computational domain. For all geometries, we have found the existence of a Dimits shift for the particle flux, i.e., the existence of a region in parameter space in which linear instability is predicted yet near-zero particle flux is measured, and of a sudden increase in the particle flux as the density gradient or the parallel electron resistivity exceeds the values corresponding to the Dimits threshold for our problem. Below the Dimits threshold, the flow is quasi-laminar, in an almost pure zonal state. The energy in the non-zonal modes, which are driven by the resistive-drift instability, is almost fully transferred to the zonal modes by the secondary instability.

For all geometries, we also observe signatures of non-local transport mechanisms once the Dimits threshold is crossed, characterized by avalanches reminiscent of the heat flux avalanches reported in gyrokinetic simulations of the ITG instability.^10–13 These are large amplitude, large and meso-scale disordered radially propagating bursts, which are responsible for the strong intermittency measured in the time series of the particle flux. For domains which are extended in the radial direction, we emphasized the rich multi-scale dynamics we identified in our simulations, where the large scale structure of the avalanches evolves in time in the transient phase during which zonal jets merge and interact with each other. In the channel domain geometry, where strong velocity shear and large density gradients are naturally created by the turbulent dynamics in the vicinity of the radial edges of the computational domain, coherent solitary structures with constant speed of propagation are seen to frequently emerge from the particle flux avalanches. Our results suggest that avalanches and soliton-like solutions have a common origin. The coherent solitary structures are responsible for larger amplitude bursts measured in the total particle transport.

We have made direct comparisons with the modified Hasegawa–Wakatani model (MHW) for both geometries. Particle flux avalanches are also seen in the MHW model, but without significant large-scale dynamics, which is consistent with the absence of multi-jet interactions and jet mergings in that model, except during a short initial transient phase. In the MHW model, zonal flows also reduce the level of particle transport and can fully impede transport for certain choices of the physical parameters, but the transition from low transport to high transport is gradual as one increases the density gradient or parallel resistivity. In that sense, we can say that the MHW model does not have a true Dimits shift, unlike the BHW model. Another major difference with the BHW model is the absence of soliton-like solutions accompanying the avalanches in the particle flux. These differences highlight the key role of the treatment of the electron parallel dynamics in reduced fluid models to capture the key mechanisms of drift turbulence driven transport in magnetically confined devices. This observation motivates a future study in which one would rigorously derive the MHW and BHW models from kinetic theory, and in that process identify the origin of the difference in the electron parallel dynamics, which could come from a different choice of closure or a different choice of orderings, leading to different terms which are neglected. We note that insightful work in that vein was presented in Ref. 21, which analyzed differences in the generation of zonal jets in the OHW and MHW models starting from a reduced kinetic model for trapped ions, but the BHW model was not considered in that article.

Our results demonstrate that the BHW model is one of the simplest model supporting a drift instability and capturing all at once the existence of a Dimits shift, avalanches, and coherent solitary structures. It is therefore a convenient model to complete our theoretical understanding of the Dimits shift, for which much progress was recently made with other reduced fluid models.^16,24,25 When the BHW model is considered in the channel domain geometry, it is also a useful model for the study of soliton-like solutions, and their link with avalanches, which are topics which are even less well studied and understood than the Dimits shift. Specifically, one may wonder if all soliton-like solutions^{12,14–16,27} observed in models of drift wave turbulence are ferdinons^14,16 and if avalanches are made of multiple ferdinons which interact with each other, or if the radially propagating structures in avalanches are inherently different from the soliton-like solutions. Finally, with the rich multi-scale dynamics and abrupt phase transitions supported by the BHW model, it is an ideal testbed to investigate the capability of novel model reduction techniques to construct inexpensive models able to accurately predict the level of turbulent transport in magnetically confined plasmas. We are currently working on all of these questions, and expect to present results in the near future.

ACKNOWLEDGMENTS

This research of A.J.M. was partially supported by the Office of Naval Research N00014-19-1-2286. D.Q. was supported as a postdoctoral fellow on the 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 and the Simons Foundation/SFARI(560651, AB).

DATA AVAILABILITY

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

APPENDIX A: COMPARISON WITH A REDUCED FLUID MODEL FOR CURVATURE-DRIVEN ITG TURBULENCE

In a recent paper,¹⁶ a two-dimensional fluid model for ion-temperature-gradient (ITG) instability in a Z-pinch geometry has been proposed. In the collisionless case, the Z-pinch ITG-driven turbulence model can be formulated in terms of the potential vorticity $q = \nabla^{2} φ - \tilde{φ}$ and the ion temperature disturbance T according to

\frac{\partial q}{\partial t} + J (φ, q) - κ_{T} \frac{\partial Δ φ}{\partial y} = \nabla \cdot J (\nabla φ, T) - \frac{\partial (φ + T)}{\partial y},

(A1a)

\frac{\partial T}{\partial t} + J (φ, T) + κ_{T} \frac{\partial φ}{\partial y} = 0,

(A1b)

with $J (f, g) = \partial_{x} f \partial_{y} g - \partial_{y} f \partial_{x} g$ the Jacobian operator, and κ_T the background temperature gradient. In contrast, the BHW equations (1) expressed in terms of the flux-balanced potential vorticity $q = \nabla^{2} φ - \tilde{n}$ and the particle density disturbance n follow the dynamics given by

\frac{\partial q}{\partial t} + J (φ, q) - κ \frac{\partial φ}{\partial y} = 0,

(A2a)

\frac{\partial n}{\partial t} + J (φ, n) + κ \frac{\partial φ}{\partial y} = α (\tilde{φ} - \tilde{n}) .

(A2b)

Both models (A1) and (A2) are defined on the 2D embedded magnetic surfaces with x as the radial direction and y as the poloidal direction. The Z-pinch ITG model has a reflection symmetry such that the equations and momentum flux are invariant under the reflection

x \mapsto - x, y \mapsto y, φ \mapsto - φ, T \mapsto - T .

This symmetry is naturally satisfied with the channel geometry boundary conditions in the BHW model, where the odd function extension in the x-direction (13) is applied to the physical quantities. Furthermore, the solid boundary walls in the channel geometry form strong mixing barriers that enhance the shear strength near the boundaries (see zonal profiles of $\partial_{x} \bar{v}$ in Fig. 4).

The Z-pinch ITG model includes the temperature gradient instability which is not present in the BHW model. The induced heat flux from the temperature gradient in the ITG model is more directly relevant for the study of the Dimits transition observed in detailed gyrokinetic simulations.⁸ Still, the analysis and numerical results in Ref. 16 reveal many similarities with the results we reported for the channel domain BHW model simulations, such as the Dimits transition occurring as one varies the temperature gradient or the density gradient, respectively, self-generated sheared zonal structures, and more remarkably the generation of localized coherent structures traveling at constant speed, which the authors in Ref. 16 call ferdinons. Specifically, the channel boundary creates strong zonal flow shear driving the ferdinons drifting toward the domain edge, comparable to the analytic results in the Z-pinch ITG model across the shear zones. Particularly strong zonal shear is observed near the boundaries (thus strong zonal transport) in the turbulent regime, while much weaker shear (and suppressed zonal transport) is observed in the zonal flow dominated regime (see Fig. 4 for the zonal profiles).

These observations stress the fact that the particle flux in the BHW model with the channel geometry plays a role comparable to the heat flux in the Z-pinch ITG fluid model. Below we use a table (Table II) to illustrate the representative features observed in the two models.

TABLE II.

Comparison of the structures and characteristics of the channel domain BHW model and periodic domain Z-pinch ITG model.

	Flux-balanced Hasegawa–Wakatani model	Z-pinch ITG-driven model
Model approximation	Finite electron resistivity α, instability from density gradient κ	Adiabatic electron response, instability from temperature gradient κ_T
Limiting model	Converges to the MHM model in Eq. (2) at adiabatic limit $α \to \infty$	Converges to the MHM model at zero temperature gradient $κ_{T} = 0$ and T = 0
Boundary condition	Solid wall boundary condition as Eq. (13) in x direction through odd extension	Doubly periodic boundary condition with reflection symmetry in Eq. (3.4) of Ref. 16
Linear instability		Both the resistive drift instability in Eq. (5) and the ITG instability in (2.22) of Ref. 16 occur in a bounded region dominated by the non-zonal modes along k_x = 0
Conservation law	$\frac{d}{d t} \frac{1}{2} \int q^{2} = κ \int \tilde{u} \tilde{n}$	$\frac{d}{d t} \frac{1}{2} \int T^{2} = κ_{T} \int \tilde{u} \tilde{T}$
Zonal flow equation	$\partial_{t} (\partial_{x}^{2} \bar{φ}) = \partial_{x} (\bar{\tilde{u} \tilde{n}}) + \partial_{x}^{2} (\bar{\tilde{u} \tilde{v}})$	$\partial_{t} (\partial_{x}^{2} \bar{φ}) = \partial_{x}^{2} \bar{(\partial_{x} \tilde{u}) \tilde{T}} + \partial_{x}^{2} (\bar{\tilde{u} \tilde{v}})$
Particle flux $\bar{\tilde{u} \tilde{n}}$	(i) Changes of total particle flux with κ in Fig. 1 and heat flux with κ_T in Fig. 6 of Ref. 16 both indicate the Dimits transition interacting with zonal flows
v.s. heat flux $\bar{\tilde{u} \tilde{T}}$	(ii) Intermittent bursts in the zonal fluxes are observed in turbulent regime in both the particle flux in Fig. 3 and the heat flux in Fig. 11 of Ref. 16
Zonal flow $\bar{v}$ and	(i) Zonal staircases in zonal flow shears arise in both models with similar shapes (see Fig. 4 in the main text and Fig. 8 of Ref. 16)
Zonal shear $\partial_{x} \bar{v}$	(ii) Bursting structures can be observed in shear zones of the zonal mean flow time-series in Fig. 5 compared with Fig. 11 of Ref. 16
Soliton-like ferdinons	Soliton-like structures (ferdinons) are created across the shear flow zone (see Fig. 5 in the main text and Fig. 10 of Ref. 16).

	Flux-balanced Hasegawa–Wakatani model	Z-pinch ITG-driven model
Model approximation	Finite electron resistivity α, instability from density gradient κ	Adiabatic electron response, instability from temperature gradient κ_T
Limiting model	Converges to the MHM model in Eq. (2) at adiabatic limit $α \to \infty$	Converges to the MHM model at zero temperature gradient $κ_{T} = 0$ and T = 0
Boundary condition	Solid wall boundary condition as Eq. (13) in x direction through odd extension	Doubly periodic boundary condition with reflection symmetry in Eq. (3.4) of Ref. 16
Linear instability		Both the resistive drift instability in Eq. (5) and the ITG instability in (2.22) of Ref. 16 occur in a bounded region dominated by the non-zonal modes along k_x = 0
Conservation law	$\frac{d}{d t} \frac{1}{2} \int q^{2} = κ \int \tilde{u} \tilde{n}$	$\frac{d}{d t} \frac{1}{2} \int T^{2} = κ_{T} \int \tilde{u} \tilde{T}$
Zonal flow equation	$\partial_{t} (\partial_{x}^{2} \bar{φ}) = \partial_{x} (\bar{\tilde{u} \tilde{n}}) + \partial_{x}^{2} (\bar{\tilde{u} \tilde{v}})$	$\partial_{t} (\partial_{x}^{2} \bar{φ}) = \partial_{x}^{2} \bar{(\partial_{x} \tilde{u}) \tilde{T}} + \partial_{x}^{2} (\bar{\tilde{u} \tilde{v}})$
Particle flux $\bar{\tilde{u} \tilde{n}}$	(i) Changes of total particle flux with κ in Fig. 1 and heat flux with κ_T in Fig. 6 of Ref. 16 both indicate the Dimits transition interacting with zonal flows
v.s. heat flux $\bar{\tilde{u} \tilde{T}}$	(ii) Intermittent bursts in the zonal fluxes are observed in turbulent regime in both the particle flux in Fig. 3 and the heat flux in Fig. 11 of Ref. 16
Zonal flow $\bar{v}$ and	(i) Zonal staircases in zonal flow shears arise in both models with similar shapes (see Fig. 4 in the main text and Fig. 8 of Ref. 16)
Zonal shear $\partial_{x} \bar{v}$	(ii) Bursting structures can be observed in shear zones of the zonal mean flow time-series in Fig. 5 compared with Fig. 11 of Ref. 16
Soliton-like ferdinons	Soliton-like structures (ferdinons) are created across the shear flow zone (see Fig. 5 in the main text and Fig. 10 of Ref. 16).

References

1.

G. W.

Hammett

,

M. A.

Beer

,

W.

Dorland

,

S. C.

Cowley

, and

S. A.

Smith

,

Plasma Phys. Controlled Fusion

35

,

973

(

1993

).

https://doi.org/10.1088/0741-3335/35/8/006

Google Scholar

Crossref

2.

R. E.

Waltz

,

G. D.

Kerbel

, and

J.

Milovich

,

Phys. Plasmas

1

,

2229

(

1994

).

https://doi.org/10.1063/1.870934

Google Scholar

Crossref

3.

M.

Beer

, “

Gyrofluid models of turbulent transport in tokamaks

,” Ph.D. thesis (

Princeton University

,

1994

).

Google Scholar

4.

W.

Horton

,

Rev. Mod. Phys.

71

,

735

(

1999

).

https://doi.org/10.1103/RevModPhys.71.735

Google Scholar

Crossref

5.

B. N.

Rogers

,

W.

Dorland

, and

M.

Kotschenreuther

,

Phys. Rev. Lett.

85

,

5336

(

2000

).

https://doi.org/10.1103/PhysRevLett.85.5336

Google Scholar

Crossref

PubMed

6.

P. H.

Diamond

,

S.

Itoh

,

K.

Itoh

, and

T.

Hahm

,

Plasma Phys. Controlled Fusion

47

,

R35

(

2005

).

https://doi.org/10.1088/0741-3335/47/5/R01

Google Scholar

Crossref

7.

K.

Itoh

,

S.-I.

Itoh

,

P. H.

Diamond

,

T. S.

Hahm

,

A.

Fujisawa

,

G. R.

Tynan

,

M.

Yagi

, and

Y.

Nagashima

,

Phys. Plasmas

13

,

055502

(

2006

).

https://doi.org/10.1063/1.2178779

Google Scholar

Crossref

8.

A.

Dimits

,

B.

Cohen

,

N.

Mattor

,

W.

Nevins

,

D.

Shumaker

,

S.

Parker

, and

C.

Kim

,

Nucl. Fusion

40

,

661

(

2000

).

https://doi.org/10.1088/0029-5515/40/3Y/329

Google Scholar

Crossref

9.

A. M.

Dimits

,

G.

Bateman

,

M.

Beer

,

B.

Cohen

,

W.

Dorland

,

G.

Hammett

,

C.

Kim

,

J.

Kinsey

,

M.

Kotschenreuther

,

A.

Kritz

et al.,

Phys. Plasmas

7

,

969

(

2000

).

https://doi.org/10.1063/1.873896

Google Scholar

Crossref

10.

G.

Dif-Pradalier

,

P. H.

Diamond

,

V.

Grandgirard

,

Y.

Sarazin

,

J.

Abiteboul

,

X.

Garbet

,

P.

Ghendrih

,

A.

Strugarek

,

S.

Ku

, and

C. S.

Chang

,

Phys. Rev. E

82

,

025401

(

2010

).

https://doi.org/10.1103/PhysRevE.82.025401

Google Scholar

Crossref

11.

J.

Candy

and

R.

Waltz

,

Phys. Rev. Lett.

91

,

045001

(

2003

).

https://doi.org/10.1103/PhysRevLett.91.045001

Google Scholar

Crossref

PubMed

12.

B. F.

McMillan

,

S.

Jolliet

,

T.

Tran

,

L.

Villard

,

A.

Bottino

, and

P.

Angelino

,

Phys. Plasmas

16

,

022310

(

2009

).

https://doi.org/10.1063/1.3079076

Google Scholar

Crossref

13.

T.

Goerler

,

X.

Lapillonne

,

S.

Brunner

,

T.

Dannert

,

F.

Jenko

,

S. K.

Aghdam

,

P.

Marcus

,

B. F.

McMillan

,

F.

Merz

,

O.

Sauter

,

D.

Told

, and

L.

Villard

,

Phys. Plasmas

18

,

056103

(

2011

).

https://doi.org/10.1063/1.3567484

Google Scholar

Crossref

14.

F.

van Wyk

,

E. G.

Highcock

,

A. A.

Schekochihin

,

C. M.

Roach

,

A. R.

Field

, and

W.

Dorland

,

J. Plasma Phys.

82

,

905820609

(

2016

).

https://doi.org/10.1017/S0022377816001148

Google Scholar

Crossref

15.

B. F.

McMillan

,

C. C. T.

Pringle

, and

B.

Teaca

,

J. Plasma Phys.

84

,

905840611

(

2018

).

https://doi.org/10.1017/S0022377818001216

Google Scholar

Crossref

16.

P. G.

Ivanov

,

A. A.

Schekochihin

,

W.

Dorland

,

A. R.

Field

, and

F. I.

Parra

, “

Zonally dominated dynamics and Dimits threshold in curvature-driven ITG turbulence

,” arXiv:2004.04047 (

2020

).

Google Scholar

17.

P.

Beyer

,

S.

Benkadda

,

X.

Garbet

, and

P.

Diamond

,

Phys. Rev. Lett.

85

,

4892

(

2000

).

https://doi.org/10.1103/PhysRevLett.85.4892

Google Scholar

Crossref

PubMed

18.

S.

Benkadda

,

P.

Beyer

,

N.

Bian

,

C.

Figarella

,

O.

Garcia

,

X.

Garbet

,

P.

Ghendrih

,

Y.

Sarazin

, and

P.

Diamond

,

Nucl. Fusion

41

,

995

(

2001

).

https://doi.org/10.1088/0029-5515/41/8/304

Google Scholar

Crossref

19.

R. L.

Dewar

and

R. F.

Abdullatif

,

in Frontiers in Turbulence and Coherent Structures

(

World Scientific

,

2007

), pp.

415

–

430

.

Google Scholar

Crossref

20.

R.

Numata

,

R.

Ball

, and

R. L.

Dewar

,

Phys. Plasmas

14

,

102312

(

2007

).

https://doi.org/10.1063/1.2796106

Google Scholar

Crossref

21.

D. D.

Sarto

and

A.

Ghizzo

,

Fluids

2

,

65

(

2017

).

https://doi.org/10.3390/fluids2040065

Google Scholar

Crossref

22.

A. J.

Majda

,

D.

Qi

, and

A. J.

Cerfon

,

Phys. Plasmas

25

,

102307

(

2018

).

https://doi.org/10.1063/1.5049389

Google Scholar

Crossref

23.

D.

Qi

,

A. J.

Majda

, and

A. J.

Cerfon

,

Phys. Plasmas

26

,

082303

(

2019

).

https://doi.org/10.1063/1.5083845

Google Scholar

Crossref

24.

D. A.

St-Onge

,

J. Plasma Phys.

83

,

905830504

(

2017

).

https://doi.org/10.1017/S0022377817000708

Google Scholar

Crossref

25.

H.

Zhu

,

Y.

Zhou

, and

I.

Dodin

,

Phys. Rev. Lett.

124

,

055002

(

2020

).

https://doi.org/10.1103/PhysRevLett.124.055002

Google Scholar

Crossref

PubMed

26.

Y.

Sarazin

,

X.

Garbet

,

P.

Ghendrih

, and

S.

Benkadda

,

Phys. Plasmas

7

,

1085

(

2000

).

https://doi.org/10.1063/1.873947

Google Scholar

Crossref

27.

Y.

Zhou

,

H.

Zhu

, and

I.

Dodin

,

Plasma Phys. Controlled Fusion

62

,

045021

(

2020

).

https://doi.org/10.1088/1361-6587/ab78f3

Google Scholar

Crossref

28.

A.

Hasegawa

and

M.

Wakatani

,

Phys. Rev. Lett.

50

,

682

(

1983

).

https://doi.org/10.1103/PhysRevLett.50.682

Google Scholar

Crossref

29.

M. A.

Beer

,

S. C.

Cowley

, and

G. W.

Hammett

,

Phys. Plasmas

2

,

2687

(

1995

).

https://doi.org/10.1063/1.871232

Google Scholar

Crossref

30.

P.

Xanthopoulos

,

A.

Mischchenko

,

P.

Helander

,

H.

Sugama

, and

T.-H.

Watanabe

,

Phys. Rev. Lett.

107

,

245002

(

2011

).

https://doi.org/10.1103/PhysRevLett.107.245002

Google Scholar

Crossref

PubMed

31.

B. J.

Faber

,

M. J.

Pueschel

,

J. H. E.

Proll

,

P.

Xanthopoulos

,

P. W.

Terry

,

C. C.

Hegna

,

G. M.

Weir

,

K. M.

Likin

, and

J. N.

Talmadge

,

Phys. Plasmas

22

,

072305

(

2015

).

https://doi.org/10.1063/1.4926510

Google Scholar

Crossref

32.

M. F.

Martin

,

M.

Landreman

,

P.

Xanthopoulos

,

N. R.

Mandell

, and

W.

Dorland

,

Plasma Phys. Controlled Fusion

60

,

095008

(

2018

).

https://doi.org/10.1088/1361-6587/aad38a

Google Scholar

Crossref

33.

D.

Qi

and

A. J.

Majda

,

Chin. Ann. Math., Ser. B

40

,

869

(

2019

).

https://doi.org/10.1007/s11401-019-0165-1

Google Scholar

Crossref

34.

D.

Qi

and

A. J.

Majda

,

Phys. Plasmas

27

,

032304

(

2020

).

https://doi.org/10.1063/1.5136281

Google Scholar

Crossref

35.

A. J.

Majda

and

X. T.

Tong

,

Nonlinearity

28

,

4171

(

2015

).

https://doi.org/10.1088/0951-7715/28/11/4171

Google Scholar

Crossref

36.

W.

Dorland

, “

Gyrofluid models of plasma turbulence

,” Ph.D. thesis (

Princeton University

,

1993

).

Google Scholar

37.

W.

Dorland

and

G. W.

Hammett

,

Phys. Fluids B

5

,

812

(

1993

).

https://doi.org/10.1063/1.860934

Google Scholar

Crossref

38.

A. V.

Pushkarev

,

W. J. T.

Bos

, and

S. V.

Nazarenko

,

Phys. Plasmas

20

,

042304

(

2013

).

https://doi.org/10.1063/1.4802187

Google Scholar

Crossref

39.

D.

Qi

and

A. J.

Majda

,

Phys. Plasmas

26

,

052108

(

2019

).

https://doi.org/10.1063/1.5095406

Google Scholar

Crossref

40.

D.

Qi

and

A. J.

Majda

,

Res. Math. Sci.

7

,

22

(

2020

).

https://doi.org/10.1007/s40687-020-00219-2

Google Scholar

Crossref

41.

D.

Qi

and

A. J.

Majda

,

J. Nonlinear Sci.

29

,

2297

(

2019

).

https://doi.org/10.1007/s00332-019-09544-5

Google Scholar

Crossref

42.

J.

Weiland

and

A.

Zagorodny

,

Rev. Mod. Plasma Phys.

3

,

8

(

2019

).

https://doi.org/10.1007/s41614-019-0029-x

Google Scholar

Crossref

2020

Author(s)

Dimits shift, avalanche-like bursts, and solitary propagating structures in the two-field flux-balanced Hasegawa–Wakatani model for plasma edge turbulence

I. INTRODUCTION

II. THE FLUX-BALANCED HASEGAWA–WAKATANI MODEL IN PLASMA EDGE TURBULENCE

A. Boundary conditions for the computational domain geometry

B. Stability analysis for the growth rate in drift wave turbulence

C. Conserved quantities and their dynamical equations

III. DIMITS SHIFT AND STATISTICAL TRANSITION IN THE CHANNEL GEOMETRY

A. Numerical treatment of the boundary conditions

B. Sharp transport transition from zonal jet dominated regime to drift wave turbulence dominated regime

C. Detailed analysis of the transport transition

IV. AVALANCHE-LIKE STRUCTURES IN DOMAINS WITH DIFFERENT ASPECT RATIOS

A. Avalanche-like structures in the radial particle flux

B. The effect of different domain aspect ratios

1. Multi-scale jet dynamics and large scale avalanche evolution

2. Statistical invariance with different aspect ratios

V. DIMITS SHIFT AND AVALANCHE STRUCTURES IN A DOUBLY PERIODIC GEOMETRY

A. Dimits shift and comparison between the BHW and MHW models

B. Transport structures for doubly periodic boundary conditions

VI. SUMMARY

ACKNOWLEDGMENTS

DATA AVAILABILITY

APPENDIX A: COMPARISON WITH A REDUCED FLUID MODEL FOR CURVATURE-DRIVEN ITG TURBULENCE

References

Citing articles via

Submit your article

Sign up for alerts

Resources

Explore

pubs.aip.org

Connect with AIP Publishing

Dimits shift, avalanche-like bursts, and solitary propagating structures in the two-field flux-balanced Hasegawa–Wakatani model for plasma edge turbulence Free

I. INTRODUCTION

II. THE FLUX-BALANCED HASEGAWA–WAKATANI MODEL IN PLASMA EDGE TURBULENCE

A. Boundary conditions for the computational domain geometry

B. Stability analysis for the growth rate in drift wave turbulence

C. Conserved quantities and their dynamical equations

III. DIMITS SHIFT AND STATISTICAL TRANSITION IN THE CHANNEL GEOMETRY

A. Numerical treatment of the boundary conditions

B. Sharp transport transition from zonal jet dominated regime to drift wave turbulence dominated regime

C. Detailed analysis of the transport transition

IV. AVALANCHE-LIKE STRUCTURES IN DOMAINS WITH DIFFERENT ASPECT RATIOS

A. Avalanche-like structures in the radial particle flux

B. The effect of different domain aspect ratios

1. Multi-scale jet dynamics and large scale avalanche evolution

2. Statistical invariance with different aspect ratios

V. DIMITS SHIFT AND AVALANCHE STRUCTURES IN A DOUBLY PERIODIC GEOMETRY

A. Dimits shift and comparison between the BHW and MHW models

B. Transport structures for doubly periodic boundary conditions

VI. SUMMARY

ACKNOWLEDGMENTS

DATA AVAILABILITY

APPENDIX A: COMPARISON WITH A REDUCED FLUID MODEL FOR CURVATURE-DRIVEN ITG TURBULENCE

References

Citing articles via

Submit your article

Sign up for alerts

Related Content

Resources

Explore

pubs.aip.org

Connect with AIP Publishing

This Feature Is Available To Subscribers Only

Dimits shift, avalanche-like bursts, and solitary propagating structures in the two-field flux-balanced Hasegawa–Wakatani model for plasma edge turbulence