Non-Oberbeck–Boussinesq zonal flow generation

Self-organization from turbulent to coherent states is a ubiquitous process in fluids. In particular, much interest and effort has been drawn to the formation of zonal flows (ZFs) [1–3]. These coherent flows arise in atmospheres, in the form of banded cloud structures on Jupiter [4], Saturn's north-polar hexagon [5] or mid-latitude westerlies on earth and in the ocean as stationary jets [6]. In magnetized fusion plasmas ZFs are key players for the reduction of the radial transport of particles and heat and for the transition to improved confinement regimes in tokamaks [7–12].

Reynolds stress is quintessential for ZF generation in all fluids [1–3, 13–17], but in magnetized plasmas also other stresses like the Maxwell [18, 19] or the diamagnetic stress [20, 21] can become significant. Virtually all of the work on ZF theory so far rely on $\newcommand{\flucRrel}[1]{{\delta{#1}}} \newcommand{\re}{{\rm Re}} \renewcommand{\ne}{{n}} \newcommand{\flucR}[1]{{\widetilde{#1}}} \flucRrel{f}$ models [1, 13, 22], which invoke the so called Oberbeck–Boussinesq (or thin layer) approximation [23, 24]. However, the latter breaks down, if the background density varies over more than one order of magnitude or if the relative density fluctuations exceed roughly 10%. This for example prevails in the edge of tokamak fusion plasmas, where experimental measurements typically feature relative density fluctuation levels around the order 0.1 in the edge and up to unity at the last closed flux surface [25–33]. Moreover, typical edge background density gradient (e-folding) lengths reach from $\newcommand{\re}{{\rm Re}} \renewcommand{\ne}{{n}} \newcommand{\rhoN}{{\rho_{s0}}} 50 \rhoN$ in low-confinement to $\newcommand{\re}{{\rm Re}} \renewcommand{\ne}{{n}} \newcommand{\rhoN}{{\rho_{s0}}} 10 \rhoN$ in high-confinement tokamak plasmas [34, 35]. Here, $\newcommand{\re}{{\rm Re}} \renewcommand{\ne}{{n}} \newcommand{\rhoN}{{\rho_{s0}}} \newcommand{\teN}{{T_{e0}}} \rhoN:=\sqrt{\teN m_i }/(e B_0)$ is the drift scale with reference electron temperature $\newcommand{\re}{{\rm Re}} \renewcommand{\ne}{{n}} \newcommand{\teN}{{T_{e0}}} \teN$, ion particle charge e, ion mass m_i and reference magnetic field B₀.

Non-Oberbeck–Boussinesq (NOB) effects on ZF generation are an unresolved issue. However, theoretical and experimental studies of poloidal ZFs in the edge of fusion plasmas indicate that unknown mechanisms beyond the Reynolds stress exist [36] and that steep background density gradients and large relative density fluctuations affect the poloidal ZF dynamics [15, 37, 38]. Moreover, the importance of large relative density fluctuations for toroidal momentum transport, as suggested by theoretical estimates in the strong and weak turbulence regime [39, 40] and experimental measurements in the TORPEX and PANTA device [41, 42], point towards a similar significance for poloidal momentum transport.

In the following we generalize the theory of ZFs to NOB effects. To this end, we decompose the density and electric potential of a full-F gyro-fluid model of a magnetized plasma [43] with the help of a density weighted Favre average [44]. This well known decomposition strategy in compressible fluid dynamics (see e.g. [45]) is here for the first time introduced to plasma physics and enables us to disentangle the density fluctuations from the ZF dynamics, while retaining the relevant physical effects. As a result, we identify novel agents in the poloidal ZF dynamics, which become significant for high relative density fluctuations or steep background density gradients. We confirm the herein proposed NOB mechanism for radial advection of ZFs with the help of numerical simulations of a fully nonlinear model for drift wave-ZF dynamics. The exploited model is based on the specified extension of the Hasegawa–Wakatani model to the full-F framework. Additionally, we show how the ZF dynamics is distributed among the proposed NOB actors and provide scalings with collisionality, reference background density gradient length and the maximum of the relative density fluctuation amplitude.

2.1. ${ \newcommand{\flucRrel}[1]{{\delta{#1}}} \newcommand{\re}{{\rm Re}} \renewcommand{\ne}{{n}} \newcommand{\flucR}[1]{{\widetilde{#1}}} \flucRrel{f}}$ formalism

We start our discussion with a short re-derivation of the conventional ZF equation and Reynolds stress from a cold ion $\newcommand{\flucRrel}[1]{{\delta{#1}}} \newcommand{\re}{{\rm Re}} \renewcommand{\ne}{{n}} \newcommand{\flucR}[1]{{\widetilde{#1}}} \flucRrel{f}$ gyro-fluid model, which couples small relative density fluctuations to the electric potential via $\renewcommand{\ne}{{n}} \newcommand{\re}{{\rm Re}} \renewcommand{\vec }[1]{{\mathbf{#1}}} \vec {E}\times\vec {B}$ advection and linear polarization [46–48]

$$\begin{align} \newcommand{\e}{{\rm e}} \renewcommand{\vec }[1]{{\mathbf{#1}}} \newcommand{\re}{{\rm Re}} \renewcommand{\ne}{{n}} \newcommand{\flucR}[1]{{\widetilde{#1}}} \newcommand{\flucRrel}[1]{{\delta{#1}}} \displaystyle \label{eq:deltafdtne} \qquad \frac{\partial}{\partial t }\flucRrel{\ne} + \vec {\nabla} \cdot \left(\flucRrel{\ne} \vec {u}_E \right) +\frac{1}{L_n B_0} \frac{\partial }{\partial y} \phi = \Lambda_\delta, \nonumber \end{align} \tag{ 1a }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \renewcommand{\vec }[1]{{\mathbf{#1}}} \newcommand{\re}{{\rm Re}} \newcommand{\Ni}{{N}} \renewcommand{\ne}{{n}} \newcommand{\flucR}[1]{{\widetilde{#1}}} \newcommand{\flucRrel}[1]{{\delta{#1}}} \displaystyle \label{eq:deltafdtNi} \qquad \frac{\partial}{\partial t }\flucRrel{\Ni} + \vec {\nabla} \cdot \left(\flucRrel{\Ni} \vec {u}_E \right) +\frac{1}{L_n B_0} \frac{\partial }{\partial y}\phi = 0, \nonumber \end{align} \tag{ 1b }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \renewcommand{\vec }[1]{{\mathbf{#1}}} \newcommand{\re}{{\rm Re}} \newcommand{\OmegaciN}{{\Omega_0}} \newcommand{\Ni}{{N}} \renewcommand{\ne}{{n}} \newcommand{\flucR}[1]{{\widetilde{#1}}} \newcommand{\flucRrel}[1]{{\delta{#1}}} \displaystyle \label{eq:dfpolcold} \qquad \vec {\nabla} \cdot \left(\frac{1}{\OmegaciN} \frac{\vec {\nabla}_\perp \phi}{B_0}\right)= \flucRrel{\ne} - \flucRrel{\Ni}. \nonumber \end{align} \tag{ 1c }$$

Here, $\newcommand{\flucRrel}[1]{{\delta{#1}}} \newcommand{\flucR}[1]{{\widetilde{#1}}} \newcommand{\re}{{\rm Re}} \renewcommand{\ne}{{n}} \flucRrel{\ne}:= \ne/n_G-1$ is the relative electron density fluctuation, $\newcommand{\flucRrel}[1]{{\delta{#1}}} \newcommand{\flucR}[1]{{\widetilde{#1}}} \newcommand{\re}{{\rm Re}} \renewcommand{\ne}{{n}} \newcommand{\Ni}{{N}} \flucRrel{\Ni}:= \Ni/n_G-1$ is the relative ion gyro-center density fluctuation, ϕ is the electric potential and $\newcommand{\re}{{\rm Re}} \renewcommand{\ne}{{n}} \newcommand{\OmegaciN}{{\Omega_0}} \OmegaciN := e B_0/m_i$ is the ion gyro-frequency. The reference background density n_G(x) refers to a constant reference background gradient length $\newcommand{\re}{{\rm Re}} \renewcommand{\ne}{{n}} \newcommand{\neref}{{n_0}} L_n:= -1/\partial_x \ln{(n_G/\neref)}$ with constant reference density $\newcommand{\re}{{\rm Re}} \renewcommand{\ne}{{n}} \newcommand{\neref}{{n_0}} \neref$. For the sake of simplicity the magnetic field B  =  B₀ is assumed constant and the unit vector in the magnetic field direction is $\renewcommand{\ne}{{n}} \newcommand{\re}{{\rm Re}} \renewcommand{\vec }[1]{{\mathbf{#1}}} \newcommand{\bhat}{\hat{\vec b}} \bhat := \vec {B}/B_0 = {\hat {\vec {e}}}_z$. The perpendicular gradient and the $\renewcommand{\ne}{{n}} \newcommand{\re}{{\rm Re}} \renewcommand{\vec }[1]{{\mathbf{#1}}} \vec {E}\times\vec {B}$ drift velocity are defined by $\renewcommand{\ne}{{n}} \newcommand{\re}{{\rm Re}} \renewcommand{\vec }[1]{{\mathbf{#1}}} \newcommand{\bhat}{\hat{\vec b}} \vec {\nabla}_\perp := - \bhat \times (\bhat \times \vec {\nabla})$ and $\renewcommand{\ne}{{n}} \newcommand{\re}{{\rm Re}} \renewcommand{\vec }[1]{{\mathbf{#1}}} \newcommand{\bhat}{\hat{\vec b}} \vec {u}_E:= \bhat \times \vec {\nabla} \phi/B_0$, respectively. The term $\Lambda_\delta$ denotes a closure for the parallel dynamics, which is discussed later in more detail. Taking the time derivative over the polarization equation (1c) yields the $\newcommand{\flucRrel}[1]{{\delta{#1}}} \newcommand{\re}{{\rm Re}} \renewcommand{\ne}{{n}} \newcommand{\flucR}[1]{{\widetilde{#1}}} \flucRrel{f}$ drift-fluid vorticity density equation

$$\begin{align} \newcommand{\e}{{\rm e}} \renewcommand{\vec }[1]{{\mathbf{#1}}} \newcommand{\re}{{\rm Re}} \newcommand{\vordf}{{\mathcal{W}_\delta}} \newcommand{\OmegaciN}{{\Omega_0}} \newcommand{\neref}{{n_0}} \renewcommand{\ne}{{n}} \displaystyle \label{eq:dtvordf} \qquad \frac{\partial}{\partial t } \vordf + \vec {\nabla} \cdot \left(\vordf\vec {u}_E \right) = \OmegaciN \neref \Lambda_\delta, \nonumber \end{align} \tag{ 2 }$$

with the linear $\renewcommand{\ne}{{n}} \newcommand{\re}{{\rm Re}} \renewcommand{\vec }[1]{{\mathbf{#1}}} \vec {E}\times\vec {B}$ vorticity density $\renewcommand{\ne}{{n}} \newcommand{\neref}{{n_0}} \newcommand{\vordf}{{\mathcal{W}_\delta}} \newcommand{\re}{{\rm Re}} \renewcommand{\vec }[1]{{\mathbf{#1}}} \newcommand{\bhat}{\hat{\vec b}} \vordf :=n_0 \nabla_\perp^2 \phi / B_0=$$\renewcommand{\ne}{{n}} \newcommand{\neref}{{n_0}} \newcommand{\vordf}{{\mathcal{W}_\delta}} \newcommand{\re}{{\rm Re}} \renewcommand{\vec }[1]{{\mathbf{#1}}} \newcommand{\bhat}{\hat{\vec b}} \bhat \cdot \vec {\nabla} \times \left(\neref \vec {u}_E\right)$. Now we apply the average over the 'poloidal' y coordinate $\newcommand{\re}{{\rm Re}} \renewcommand{\ne}{{n}} \newcommand{\avgR}[1]{{\langle #1 \rangle}} \avgR{h} := L_y^{-1}\int_0^{L_y} {\rm d}y \hspace{1mm} h$ to equation (2), which is the 2D equivalent of a flux surface average. Reynolds decomposition $\newcommand{\avgR}[1]{{\langle #1 \rangle}} \newcommand{\re}{{\rm Re}} \renewcommand{\ne}{{n}} \newcommand{\flucR}[1]{{\widetilde{#1}}} h = \avgR{h} + \flucR{h}$ and integration over the 'radial' coordinate x result in the $\newcommand{\flucRrel}[1]{{\delta{#1}}} \newcommand{\re}{{\rm Re}} \renewcommand{\ne}{{n}} \newcommand{\flucR}[1]{{\widetilde{#1}}} \flucRrel{f}$ evolution equation for poloidal ZFs [13]

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\OmegaciN}{{\Omega_0}} \renewcommand{\ne}{{n}} \newcommand{\re}{{\rm Re}} \newcommand{\flucR}[1]{{\widetilde{#1}}} \newcommand{\avgR}[1]{{\langle #1 \rangle}} \displaystyle \label{eq:deltafmeanflow} \qquad \frac{\partial }{\partial t} \avgR{u_{y}} = - \frac{\partial}{\partial x} \avgR{\flucR{u_{x} } \flucR{u_{y}} } + \OmegaciN\int_{x_0}^x {\rm d} x \avgR{\Lambda_\delta } . \nonumber \end{align} \tag{ 3 }$$

Here, we introduced $u_{x} := - \partial_y \phi/B_0$, $u_{y} := \partial_x \phi/B_0$ and the anticipated Reynolds stress $\newcommand{\avgR}[1]{{\langle #1 \rangle}} \newcommand{\re}{{\rm Re}} \renewcommand{\ne}{{n}} \newcommand{\flucR}[1]{{\widetilde{#1}}} \mathcal{R} := \avgR{\flucR{u_{x}} \flucR{u_{y}} }$ [49], where $\newcommand{\avgR}[1]{{\langle #1 \rangle}} \newcommand{\re}{{\rm Re}} \renewcommand{\ne}{{n}} \newcommand{\flucR}[1]{{\widetilde{#1}}} \avgR{u_{x} u_{y} } = \avgR{u_{x} } \avgR{u_{y} } + \avgR{\flucR{u_{x}} \flucR{u_{y}} }$ and $\newcommand{\re}{{\rm Re}} \renewcommand{\ne}{{n}} \newcommand{\avgR}[1]{{\langle #1 \rangle}} \avgR{u_{x} } = 0$ was used. In passing we note that we assume that radial boundary conditions give rise to no additional terms in equation (3) and for the remainder of this letter.

2.2. Full-F formalism

In full-F theory the splitting of the gyro-fluid moment variables into fluctuating and background parts is avoided and the quasi-neutrality constraint for electrons and ions is rendered by the nonlinear polarization equation [43]. The cold ion full-F gyro-fluid model [50, 51]

$$\begin{align} \newcommand{\e}{{\rm e}} \renewcommand{\vec }[1]{{\mathbf{#1}}} \newcommand{\re}{{\rm Re}} \renewcommand{\ne}{{n}} \displaystyle \label{eq:fullFdtne} \qquad \frac{\partial}{\partial t }\ne +\vec {\nabla} \cdot \left(\ne \vec {u}_E\right) =\Lambda, \nonumber \end{align} \tag{ 4a }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \renewcommand{\vec }[1]{{\mathbf{#1}}} \newcommand{\re}{{\rm Re}} \newcommand{\Ni}{{N}} \renewcommand{\ne}{{n}} \displaystyle \label{eq:fullFdtNi} \qquad \frac{\partial}{\partial t } \Ni +\vec {\nabla} \cdot \left(\Ni \vec {U}_E\right) =0 , \nonumber \end{align} \tag{ 4b }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \renewcommand{\vec }[1]{{\mathbf{#1}}} \newcommand{\re}{{\rm Re}} \newcommand{\OmegaciN}{{\Omega_0}} \newcommand{\Ni}{{N}} \renewcommand{\ne}{{n}} \displaystyle \label{eq:ffpolcoldlwl} \qquad \vec {\nabla} \cdot \left(\frac{\Ni}{\OmegaciN} \frac{\vec {\nabla}_\perp \phi}{B_0}\right)= \ne - \Ni, \nonumber \end{align} \tag{ 4c }$$

evolves the full electron density $\newcommand{\re}{{\rm Re}} \renewcommand{\ne}{{n}} {\ne}$ and ion gyro-center density $\newcommand{\re}{{\rm Re}} \renewcommand{\ne}{{n}} \newcommand{\Ni}{{N}} \Ni$. In the gyro-center $\renewcommand{\ne}{{n}} \newcommand{\re}{{\rm Re}} \renewcommand{\vec }[1]{{\mathbf{#1}}} \vec {E}\times\vec {B}$ drift velocity by $\renewcommand{\ne}{{n}} \newcommand{\re}{{\rm Re}} \renewcommand{\vec }[1]{{\mathbf{#1}}} \vec {U}_E:= \vec {u}_E + \vec {U}_p$ the ponderomotive correction $\renewcommand{\ne}{{n}} \newcommand{\OmegaciN}{{\Omega_0}} \newcommand{\re}{{\rm Re}} \renewcommand{\vec }[1]{{\mathbf{#1}}} \newcommand{\bhat}{\hat{\vec b}} \vec {U}_p :=-\bhat \times \vec {\nabla} \vec {u}_E^2/(2 \OmegaciN)$ appears. Both, the latter ponderomotive correction and the polarization charge nonlinearity on the left hand side of equation (4c) are crucial for energetic consistency and an exact momentum conservation law [52]. We refer to the parallel closure term Λ later on. In the long wavelength limit we can again reformulate equations (4b) and (4c) into a drift-fluid vorticity density equation

$$\begin{align} \newcommand{\e}{{\rm e}} \renewcommand{\vec }[1]{{\mathbf{#1}}} \newcommand{\re}{{\rm Re}} \newcommand{\vorfF}{{\mathcal{W}}} \newcommand{\OmegaciN}{{\Omega_0}} \renewcommand{\ne}{{n}} \displaystyle \label{eq:dtvorfullF} \qquad \frac{\partial}{\partial t } \vorfF + \vec {\nabla} \cdot \left(\vorfF\vec {u}_E \right) - \OmegaciN \vec {\nabla} \cdot \left(n \vec {U}_p \right) =\OmegaciN \Lambda, \nonumber \end{align} \tag{ 5 }$$

where the nonlinear $\renewcommand{\ne}{{n}} \newcommand{\re}{{\rm Re}} \renewcommand{\vec }[1]{{\mathbf{#1}}} \vec {E}\times\vec {B}$ vorticity density is given by $\renewcommand{\ne}{{n}} \newcommand{\vorfF}{{\mathcal{W}}} \newcommand{\re}{{\rm Re}} \renewcommand{\vec }[1]{{\mathbf{#1}}} \newcommand{\bhat}{\hat{\vec b}} \vorfF := \vec {\nabla} \cdot \left(\ne \vec {\nabla}_\perp \phi/B_0\right) = \bhat \cdot \vec {\nabla} \times \left(n \vec {u}_E\right)$. As above, we obtain the averaged poloidal momentum equation [38, 53, 54]

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\OmegaciN}{{\Omega_0}} \renewcommand{\ne}{{n}} \newcommand{\re}{{\rm Re}} \newcommand{\flucR}[1]{{\widetilde{#1}}} \newcommand{\avgR}[1]{{\langle #1 \rangle}} \displaystyle \label{eq:fullfmeanflow} \qquad \frac{\partial }{\partial t} \avgR{\ne u_{y} }=& -\frac{\partial}{\partial x} \left(\avgR{\ne } \mathcal{R} +\avgR{\flucR{\ne} \flucR{u_{x}}} \avgR{u_{y} } +\avgR{\flucR{\ne} \flucR{u_{x}}\flucR{u_{y}} } \right) \nonumber \\ &+\OmegaciN\int_{x_0}^x {\rm d} x \avgR{\Lambda } . \nonumber \end{align} \tag{ 6 }$$

The divergence of the full-F stress drive terms of equation (6) is related to the averaged radial flux of vorticity density minus the ponderomotive correction via the full-F Taylor identity

$$\begin{align} \newcommand{\bi}{\boldsymbol} \newcommand{\e}{{\rm e}} \newcommand{\vorfF}{{\mathcal{W}}} \newcommand{\OmegaciN}{{\Omega_0}} \renewcommand{\ne}{{n}} \newcommand{\re}{{\rm Re}} \newcommand{\flucR}[1]{{\widetilde{#1}}} \newcommand{\avgR}[1]{{\langle #1 \rangle}} \displaystyle \label{eq:taylorfullF} \qquad \avgR{\flucR{u_x} \flucR{\vorfF}} =&\, \frac{\partial}{\partial x} \big(\avgR{\ne} {\mathcal R} + \avgR{\flucR{\ne} \flucR{u_x}} \avgR{u_y} + \avgR{\flucR{\ne} \flucR{u_x} \flucR{u_y}}\big) \nonumber \\ &+\OmegaciN \avgR{\flucR{U_{p,x}} \flucR{n}} . \nonumber \end{align} \tag{ 7 }$$

The interpretation of equation (6) is problematic since (i) absolute density fluctuations $\newcommand{\flucR}[1]{{\widetilde{#1}}} \newcommand{\re}{{\rm Re}} \renewcommand{\ne}{{n}} \flucR{\ne}$ arise instead of relative density fluctuations $\newcommand{\avgR}[1]{{\langle #1 \rangle}} \newcommand{\flucR}[1]{{\widetilde{#1}}} \newcommand{\re}{{\rm Re}} \renewcommand{\ne}{{n}} \flucR{\ne}/\avgR{\ne}$, (ii) the time evolution of the averaged poloidal momentum $\newcommand{\avgR}[1]{{\langle #1 \rangle}} \newcommand{\re}{{\rm Re}} \renewcommand{\ne}{{n}} \avgR{\ne u_{y} }$ is given in terms of the averaged poloidal velocity $\newcommand{\re}{{\rm Re}} \renewcommand{\ne}{{n}} \newcommand{\avgR}[1]{{\langle #1 \rangle}} \avgR{u_{y} }$ and (iii) background density gradient $\newcommand{\avgR}[1]{{\langle #1 \rangle}} \newcommand{\re}{{\rm Re}} \renewcommand{\ne}{{n}} \partial_x \ln{\avgR{\ne }}$ effects are not obvious. Despite these obstacles equation (6) has been recently used to show that the second term, occasionally misinterpreted as advective, and the cubic term can be comparable to the Reynolds stress related drive $\newcommand{\avgR}[1]{{\langle #1 \rangle}} \newcommand{\re}{{\rm Re}} \renewcommand{\ne}{{n}} \partial_x (\avgR{\ne} \mathcal{R})$ [15, 37, 38].

Thus, we go a step further and utilize a density weighted Favre decomposition instead of the Reynolds decomposition according to $\newcommand{\avgF}[1]{{\left[ \left[#1 \right]\right]}} \newcommand{\re}{{\rm Re}} \renewcommand{\ne}{{n}} \newcommand{\flucF}[1]{{\widehat{#1}}} h := \avgF{h} + \flucF{h}$ and $\newcommand{\avgF}[1]{{\left[ \left[#1 \right]\right]}} \newcommand{\avgR}[1]{{\langle #1 \rangle}} \newcommand{\re}{{\rm Re}} \renewcommand{\ne}{{n}} \avgF{h} :=\avgR{\ne h}/\avgR{\ne }$ [44]. Note that the Favre decomposition reduces to the Reynolds decomposition if the density $\newcommand{\re}{{\rm Re}} \renewcommand{\ne}{{n}} \ne$ is only a function of x. Now we combine the poloidal average of equation (4a) divided by $\newcommand{\avgR}[1]{{\langle #1 \rangle}} \newcommand{\re}{{\rm Re}} \renewcommand{\ne}{{n}} \avgR{\ne }$

$$\begin{align} \newcommand{\e}{{\rm e}} \renewcommand{\ne}{{n}} \newcommand{\re}{{\rm Re}} \newcommand{\avgR}[1]{{\langle #1 \rangle}} \newcommand{\avgF}[1]{{\left[ \left[#1 \right]\right]}} \displaystyle \label{eq:fullflognemeanflow} \qquad \frac{\partial }{\partial t}\ln{\avgR{\ne }} =& - \frac{\partial}{\partial x} \avgF{u_{x}} - \avgF{u_{x}}\frac{\partial}{\partial x} \ln{\avgR{\ne }} \nonumber \\&+ \frac{\avgR{\Lambda }}{\avgR{\ne }} \nonumber \end{align} \tag{ 8 }$$

with equation (6) divided by $\newcommand{\avgR}[1]{{\langle #1 \rangle}} \newcommand{\re}{{\rm Re}} \renewcommand{\ne}{{n}} \avgR{\ne }$ to obtain a ZF evolution equation for the Favre averaged poloidal velocity

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\OmegaciN}{{\Omega_0}} \renewcommand{\ne}{{n}} \newcommand{\re}{{\rm Re}} \newcommand{\avgR}[1]{{\langle #1 \rangle}} \newcommand{\flucF}[1]{{\widehat{#1}}} \newcommand{\avgF}[1]{{\left[ \left[#1 \right]\right]}} \displaystyle \label{eq:fullfmeanflowcombined} \qquad \frac{\partial }{\partial t} \avgF{u_{y}} =& - \frac{\partial}{\partial x} \avgF{\flucF{u_{x}} \flucF{u_{y}} } - \avgF{u_{x}} \frac{\partial}{\partial x} \avgF{u_{y}} \nonumber \\& - \avgF{\flucF{u_{x}} \flucF{u_{y}} } \frac{\partial}{\partial x} \ln{\avgR{\ne }} \nonumber \\& - \avgF{u_{y}} \frac{\avgR{\Lambda }}{\avgR{\ne }} + \frac{\OmegaciN}{\avgR{\ne }}\int_{x_0}^x {\rm d} x \avgR{\Lambda } , \nonumber \end{align} \tag{ 9 }$$

where we used $\newcommand{\avgF}[1]{{\left[ \left[#1 \right]\right]}} \newcommand{\re}{{\rm Re}} \renewcommand{\ne}{{n}} \newcommand{\flucF}[1]{{\widehat{#1}}} \avgF{u_{x} u_{y} } = \avgF{u_{x}} \avgF{u_{y}}+\avgF{\flucF{u_{x}} \flucF{u_{y}} }$. The Favre stress $\newcommand{\avgF}[1]{{\left[ \left[#1 \right]\right]}} \newcommand{\re}{{\rm Re}} \renewcommand{\ne}{{n}} \newcommand{\flucF}[1]{{\widehat{#1}}} \mathcal{F} := \avgF{\flucF{u_{x}}\flucF{u_{y}}}$ can be rewritten into

$$\begin{align} \newcommand{\e}{{\rm e}} \renewcommand{\ne}{{n}} \newcommand{\re}{{\rm Re}} \newcommand{\flucR}[1]{{\widetilde{#1}}} \newcommand{\avgR}[1]{{\langle #1 \rangle}} \newcommand{\avgF}[1]{{\left[ \left[#1 \right]\right]}} \displaystyle \label{eq:Favrestress} \qquad \mathcal{F} = \mathcal{R} - \avgF{\flucR{u_{x}}} \avgF{\flucR{u_{y}}} +\avgR{\flucR{\ne} \flucR{u_{x}}\flucR{u_{y}} } / \avgR{\ne } . \nonumber \end{align} \tag{ 10 }$$

Consequently, the first term $-\partial_x \mathcal{F}$ on the right hand side of equation (9) is the superposition of the conventional Reynolds stress drive $T_1 := -\partial_x \mathcal{R}$, the quadruple fluctuation term $\newcommand{\avgF}[1]{{\left[ \left[#1 \right]\right]}} \newcommand{\re}{{\rm Re}} \renewcommand{\ne}{{n}} \newcommand{\flucR}[1]{{\widetilde{#1}}} T_2 := \partial_x \left(\avgF{\flucR{u_{x}}} \avgF{\flucR{u_{y}}} \right)$ and the triple fluctuation drive $\newcommand{\avgR}[1]{{\langle #1 \rangle}} \newcommand{\flucR}[1]{{\widetilde{#1}}} \newcommand{\re}{{\rm Re}} \renewcommand{\ne}{{n}} T_3 := -\partial_x \left(\avgR{\flucR{\ne} \flucR{u_{x}}\flucR{u_{y}} } / \avgR{\ne } \right)$. The novel second term on the right hand side of equation (9) represents radial advection of poloidal ZFs $\newcommand{\re}{{\rm Re}} \renewcommand{\ne}{{n}} \newcommand{\avgF}[1]{{\left[ \left[#1 \right]\right]}} \avgF{u_{y}}$. Its direction depends on the sign of the averaged radial particle flux $\newcommand{\avgF}[1]{{\left[ \left[#1 \right]\right]}} \newcommand{\avgR}[1]{{\langle #1 \rangle}} \newcommand{\re}{{\rm Re}} \renewcommand{\ne}{{n}} \avgR{\Gamma_{x}}:= \avgR{\ne }\avgF{u_{x}}$, which is typically positive, so that T₄ describes an outward pinch of ZFs. The novel third term $\newcommand{\re}{{\rm Re}} \renewcommand{\ne}{{n}} \newcommand{\avgR}[1]{{\langle #1 \rangle}} T_5 := \mathcal{F} /L_{\avgR{n }}$ on the right hand side of equation (9) is proportional to the inverse of the background density gradient length $\newcommand{\avgR}[1]{{\langle #1 \rangle}} \newcommand{\re}{{\rm Re}} \renewcommand{\ne}{{n}} 1/L_{\avgR{n }}:= -\partial_x \ln{\avgR{\ne }}$. This term is large for small reference background density gradient lengths L_n, or has large radially localized values if the density profile $\newcommand{\re}{{\rm Re}} \renewcommand{\ne}{{n}} \newcommand{\avgR}[1]{{\langle #1 \rangle}} \avgR{n}$ develops into a staircase like pattern [55]. In contrast to the Favre stress drive $-\partial_x \mathcal{F}$, the background density gradient drive T₅ contributes to the ZF generation even if the Favre stress is radially homogeneous $\partial_x \mathcal{F}=0$. Remarkably, the background density gradient drive remains finite in the small relative density fluctuation limit, where the density n is only a function of x and the Favre stress $\mathcal{F}$ resembles the conventional Reynolds stress $\mathcal{R}$.

In order to interpret the dynamics of the background density gradient drive T₅ let us assume for a moment that the turbulent viscosity hypothesis $\newcommand{\re}{{\rm Re}} \renewcommand{\ne}{{n}} \newcommand{\avgF}[1]{{\left[ \left[#1 \right]\right]}} \mathcal{F} := -\nu_T(x)\partial_x \avgF{u_y}$ holds [13, 56]. In this case equation (9) reduces to a simple advection-diffusion equation for ZFs

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\OmegaciN}{{\Omega_0}} \renewcommand{\ne}{{n}} \newcommand{\re}{{\rm Re}} \newcommand{\avgR}[1]{{\langle #1 \rangle}} \newcommand{\avgF}[1]{{\left[ \left[#1 \right]\right]}} \displaystyle \qquad \frac{\partial }{\partial t} \avgF{u_{y}} =& - \left(\avgF{u_{x}}+ V \right)\frac{\partial}{\partial x} \avgF{u_{y}} \nonumber \\& +\frac{\partial}{\partial x} \left(\nu_T\frac{\partial}{\partial x}\avgF{u_y} \right) - \avgF{u_{y}} \frac{\avgR{\Lambda }}{\avgR{\ne }} \nonumber \\& + \frac{\OmegaciN}{\avgR{\ne }}\int_{x_0}^x {\rm d} x \avgR{\Lambda } , \nonumber \end{align} \tag{ 11 }$$

where the background density gradient pinch velocity $\newcommand{\re}{{\rm Re}} \renewcommand{\ne}{{n}} \newcommand{\avgR}[1]{{\langle #1 \rangle}} V:=\nu_T/L_{\avgR{n}}$ appears now in addition to the radial outward pinch velocity $\newcommand{\re}{{\rm Re}} \renewcommand{\ne}{{n}} \newcommand{\avgF}[1]{{\left[ \left[#1 \right]\right]}} \avgF{u_{x}}$. The direction of the additional pinch depends on the sign of the turbulent viscosity $\nu_T$.

Finally, we extend the theory for energy transfer inside the kinetic $\renewcommand{\ne}{{n}} \newcommand{\re}{{\rm Re}} \renewcommand{\vec }[1]{{\mathbf{#1}}} \vec {E}\times \vec {B}$ energy $\newcommand{\avgR}[1]{{\langle #1 \rangle}} \renewcommand{\ne}{{n}} \newcommand{\re}{{\rm Re}} \renewcommand{\vec }[1]{{\mathbf{#1}}} E(t):= m_i \int {\rm d}A\avgR{\ne\vec {u}_E^2}/2$ to the full-F formalism. Here, the Favre decomposition $E = E_0+ E_1$ is pivotal to derive the conservation laws for the zonal (or mean) $\newcommand{\avgF}[1]{{\left[ \left[#1 \right]\right]}} \newcommand{\avgR}[1]{{\langle #1 \rangle}} \newcommand{\re}{{\rm Re}} \renewcommand{\ne}{{n}} E_0(t) := m_i \int {\rm d}A\avgR{\ne}\avgF{u_y}^2/2$ and turbulent part $\newcommand{\flucF}[1]{{\widehat{#1}}} \newcommand{\avgR}[1]{{\langle #1 \rangle}} \renewcommand{\ne}{{n}} \newcommand{\re}{{\rm Re}} \renewcommand{\vec }[1]{{\mathbf{#1}}} E_1(t) := m_i \int {\rm d}A\avgR{\ne}\flucF{\vec {u}}_E^2/2$ of the kinetic $\renewcommand{\ne}{{n}} \newcommand{\re}{{\rm Re}} \renewcommand{\vec }[1]{{\mathbf{#1}}} \vec {E}\times \vec {B}$ energy and supersedes the Reynolds decomposition in the $\newcommand{\flucRrel}[1]{{\delta{#1}}} \newcommand{\re}{{\rm Re}} \renewcommand{\ne}{{n}} \newcommand{\flucR}[1]{{\widetilde{#1}}} \flucRrel{f}$ formalism [19, 21]. With the help of equations (8) and (9) we obtain the conservation laws for the zonal and turbulent kinetic $\renewcommand{\ne}{{n}} \newcommand{\re}{{\rm Re}} \renewcommand{\vec }[1]{{\mathbf{#1}}} \vec {E}\times \vec {B}$ energy

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\OmegaciN}{{\Omega_0}} \renewcommand{\ne}{{n}} \newcommand{\re}{{\rm Re}} \newcommand{\avgR}[1]{{\langle #1 \rangle}} \newcommand{\avgF}[1]{{\left[ \left[#1 \right]\right]}} \displaystyle \label{eq:dtE0} \qquad \frac{\partial}{\partial t} E_0 =& \int {\rm d}A m_i\Big(\avgR{\ne} \mathcal{F} \frac{\partial}{\partial x} \avgF{u_{y}} - \frac{\avgF{u_{y}}^2}{2} \avgR{\Lambda } \nonumber \\ & + \OmegaciN \avgF{u_{y}}\int_{x_0}^x {\rm d} x \avgR{\Lambda } \Big), \nonumber \end{align} \tag{ 12a }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\OmegaciN}{{\Omega_0}} \renewcommand{\ne}{{n}} \newcommand{\re}{{\rm Re}} \newcommand{\avgR}[1]{{\langle #1 \rangle}} \newcommand{\avgF}[1]{{\left[ \left[#1 \right]\right]}} \displaystyle \label{eq:dtE1} \qquad \frac{\partial}{\partial t} E_1 =& \int {\rm d}A m_i\Big(- \avgR{\ne} \mathcal{F} \frac{\partial}{\partial x} \avgF{u_{y}} + \frac{\avgF{u_{y}}^2}{2} \avgR{\Lambda } \nonumber \\ & - \OmegaciN \avgF{u_{y}}\int_{x_0}^x {\rm d} x \avgR{\Lambda } - \frac{e}{m_i}\avgR{\phi \Lambda}\Big) . \nonumber \end{align} \tag{ 12b }$$

This unveils that the Favre stress term $\newcommand{\avgF}[1]{{\left[ \left[#1 \right]\right]}} \newcommand{\avgR}[1]{{\langle #1 \rangle}} \newcommand{\re}{{\rm Re}} \renewcommand{\ne}{{n}} \avgR{\ne} \mathcal{F} \partial_x \avgF{u_{y}}$ is the central mechanism for energy transfer between the zonal and turbulent kinetic $\renewcommand{\ne}{{n}} \newcommand{\re}{{\rm Re}} \renewcommand{\vec }[1]{{\mathbf{#1}}} \vec {E}\times \vec {B}$ energy. As a consequence, density fluctuations (see equation (10)) manifest as an additional transfer channel in the full-F formalism.

Self-sustained drift wave turbulence is maintained by the non-adiabatic parallel coupling of the relative density fluctuations and the electric potential, which can arise due to various mechanisms. Here, we exemplarily consider resistive drift wave turbulence, which arises due to resistive friction between electrons and ions along the magnetic field line. This mechanism enters the 2D gyro-fluid models via the parallel closure terms ($\Lambda_\delta$ or Λ) of the Hasegawa–Wakatani (HW) type as summarized in table 1. Here, we introduced the full-F adiabaticity parameter $\newcommand{\re}{{\rm Re}} \renewcommand{\ne}{{n}} \newcommand{\neref}{{n_0}} \newcommand{\OmegaciN}{{\Omega_0}} \newcommand{\teN}{{T_{e0}}} \newcommand{\e}{{\rm e}} \alpha:= \teN k_\parallel^2 /(\eta_\parallel e^2 \neref \OmegaciN)$ with parallel wavenumber $k_\parallel$ and parallel Spitzer resistivity $\newcommand{\re}{{\rm Re}} \renewcommand{\ne}{{n}} \newcommand{\e}{{\rm e}} \eta_\parallel := 0.51 m_e \nu_e/(\ne e^2)$ [61, 62]. In the electron collision frequency $\nu_e$ the Coulomb logarithm is treated as a constant so that $\newcommand{\re}{{\rm Re}} \renewcommand{\ne}{{n}} \newcommand{\e}{{\rm e}} \eta_\parallel$ has no explicit dependence on $\newcommand{\re}{{\rm Re}} \renewcommand{\ne}{{n}} \ne$. As opposed to this in $\newcommand{\flucRrel}[1]{{\delta{#1}}} \newcommand{\re}{{\rm Re}} \renewcommand{\ne}{{n}} \newcommand{\flucR}[1]{{\widetilde{#1}}} \flucRrel{f}$ models the density dependence in the collision frequency $\newcommand{\re}{{\rm Re}} \renewcommand{\ne}{{n}} \nu_e (\ne) \approx \nu_{e0}$ is completely neglected so that $\newcommand{\re}{{\rm Re}} \renewcommand{\ne}{{n}} \newcommand{\OmegaciN}{{\Omega_0}} \newcommand{\teN}{{T_{e0}}} \alpha_\delta := \teN k_\parallel^2 /(0.51 m_e \nu_{e0} \OmegaciN)$ reduces to a parameter. Only then, the poloidal variations of the adiabaticity parameters vanish ($\newcommand{\re}{{\rm Re}} \renewcommand{\ne}{{n}} \newcommand{\flucR}[1]{{\widetilde{#1}}} \flucR{\alpha_\delta} = \flucR{\alpha} = 0$), and the full-F and $\newcommand{\flucRrel}[1]{{\delta{#1}}} \newcommand{\re}{{\rm Re}} \renewcommand{\ne}{{n}} \newcommand{\flucR}[1]{{\widetilde{#1}}} \flucRrel{f}$ closures coincide in the limit of $\newcommand{\avgR}[1]{{\langle #1 \rangle}} \newcommand{\re}{{\rm Re}} \renewcommand{\ne}{{n}} \avgR{\ne } \approx n_G$ and $\newcommand{\flucRrel}[1]{{\delta{#1}}} \newcommand{\flucR}[1]{{\widetilde{#1}}} \newcommand{\re}{{\rm Re}} \renewcommand{\ne}{{n}} \flucRrel{\ne} \ll 1$.

Table 1. HW closures for $\newcommand{\flucRrel}[1]{{\delta{#1}}} \newcommand{\re}{{\rm Re}} \renewcommand{\ne}{{n}} \newcommand{\flucR}[1]{{\widetilde{#1}}} \flucRrel{f}$ and full-F models.

	Ordinary HW	Modified HW
$\newcommand{\re}{{\rm Re}} \renewcommand{\ne}{{n}} \newcommand{\OmegaciN}{{\Omega_0}} \Lambda_\delta/(\alpha_\delta \OmegaciN)$	$\newcommand{\flucRrel}[1]{{\delta{#1}}} \newcommand{\flucR}[1]{{\widetilde{#1}}} \newcommand{\re}{{\rm Re}} \renewcommand{\ne}{{n}} \newcommand{\teN}{{T_{e0}}} e\phi/\teN-\flucRrel{\ne}$ [57–59]	$\newcommand{\flucRrel}[1]{{\delta{#1}}} \newcommand{\flucR}[1]{{\widetilde{#1}}} \newcommand{\re}{{\rm Re}} \renewcommand{\ne}{{n}} \newcommand{\teN}{{T_{e0}}} e\flucR{\phi}/\teN- \flucR{\flucRrel{\ne}}$ [60]
$\newcommand{\re}{{\rm Re}} \renewcommand{\ne}{{n}} \newcommand{\neref}{{n_0}} \newcommand{\OmegaciN}{{\Omega_0}} {\Lambda}/(\alpha \neref \OmegaciN)$	$\newcommand{\avgR}[1]{{\langle #1 \rangle}} \newcommand{\re}{{\rm Re}} \renewcommand{\ne}{{n}} \newcommand{\teN}{{T_{e0}}} e\phi/\teN -\ln\left(\ne/\avgR{\ne }\right)$	$\newcommand{\flucR}[1]{{\widetilde{#1}}} \newcommand{\re}{{\rm Re}} \renewcommand{\ne}{{n}} \newcommand{\teN}{{T_{e0}}} e\flucR{\phi}/\teN-\flucR{\ln\left(\ne\right)}$

We use the open source library Feltor [63] to numerically solve the full-F gyro-fluid equations (4a)–(4c) with the modified HW parallel closure of table 1. Numerical stability is ensured by adding hyperdiffusive terms of second order $\renewcommand{\ne}{{n}} \newcommand{\re}{{\rm Re}} \renewcommand{\vec }[1]{{\mathbf{#1}}} - \nu \vec {\nabla}_\perp^4 n$ and $\renewcommand{\ne}{{n}} \newcommand{\re}{{\rm Re}} \renewcommand{\vec }[1]{{\mathbf{#1}}} - \nu \vec {\nabla}_\perp^4 N$ to the right hand side of equations (4a) and (4b). Moreover, we append the right hand side of equations (4a) and (4b) by a density source of the form $\omega_S z \Theta (z)$ with $\newcommand{\avgR}[1]{{\langle #1 \rangle}} \newcommand{\re}{{\rm Re}} \renewcommand{\ne}{{n}} z:=g(x) \left(n_G - \avgR{\ne}\right)$ to maintain the initial profile in a small region $x \in \left[0, \,x_b\right]$. Here, we defined the Heaviside function $\Theta(z)$ and $g(x) := \left[1-\tanh{(x-x_b)/\sigma_b}\right]/2$. The corresponding parameters are fixed to $\newcommand{\re}{{\rm Re}} \renewcommand{\ne}{{n}} \newcommand{\rhoN}{{\rho_{s0}}} \newcommand{\csN}{{c_{s0}}} \nu=5\times 10^{-4} \csN \rhoN^3$, $\newcommand{\re}{{\rm Re}} \renewcommand{\ne}{{n}} \newcommand{\OmegaciN}{{\Omega_0}} \omega_S = 0.1 \OmegaciN$, $x_b= 0.1 L_x$ and $\newcommand{\re}{{\rm Re}} \renewcommand{\ne}{{n}} \newcommand{\rhoN}{{\rho_{s0}}} \sigma_b = 0.5 \rhoN$ with cold ion sound speed $\newcommand{\re}{{\rm Re}} \renewcommand{\ne}{{n}} \newcommand{\OmegaciN}{{\Omega_0}} \newcommand{\rhoN}{{\rho_{s0}}} \newcommand{\csN}{{c_{s0}}} \csN:=\rhoN\OmegaciN$. The box with size $\newcommand{\re}{{\rm Re}} \renewcommand{\ne}{{n}} \newcommand{\rhoN}{{\rho_{s0}}} L_x=L_y=128 \rhoN$ is resolved by a discontinuous Galerkin discretization with P  =  3 polynomial coefficients and at least $N_x=N_y=256$ equidistant grid cells. The initial (gyro-center) density fields $\newcommand{\flucRrel}[1]{{\delta{#1}}} \newcommand{\flucR}[1]{{\widetilde{#1}}} \renewcommand{\ne}{{n}} \newcommand{\re}{{\rm Re}} \renewcommand{\vec }[1]{{\mathbf{#1}}} n(\vec {x}, 0)=N(\vec {x}, 0) = n_G(x)\left(1+\flucRrel{n}_0(\vec {x})\right)$ consist of the reference background density profile n_G, which is perturbed by a turbulent bath $\newcommand{\flucRrel}[1]{{\delta{#1}}} \newcommand{\flucR}[1]{{\widetilde{#1}}} \renewcommand{\ne}{{n}} \newcommand{\re}{{\rm Re}} \renewcommand{\vec }[1]{{\mathbf{#1}}} \flucRrel{n}_0(\vec {x})$.

NOB effects on drift wave-ZF dynamics, as it is described by equations (8) and (9) with $\langle \Lambda \rangle=0$, are in this setup studied by varying the adiabaticity parameter (or inverse collisionality) α and the reference background gradient length L_n. In figure 1 we show that L_n crucially determines the time evolution of ZFs in the high collisionality regime with $\alpha=0.0005$. While stationary ZFs emerge for $\newcommand{\re}{{\rm Re}} \renewcommand{\ne}{{n}} \newcommand{\rhoN}{{\rho_{s0}}} L_n=128 \rhoN$, a radial outward pinch of ZFs occurs for a four times smaller reference background density gradient length $\newcommand{\re}{{\rm Re}} \renewcommand{\ne}{{n}} \newcommand{\rhoN}{{\rho_{s0}}} L_n=32 \rhoN$.

Zoom In Zoom Out Reset image size

In this steep gradient and high collisionality regime the ZF signature is no longer solely determined by the conventional Reynolds stress drive, which is illustrated in figure 2. The Reynolds stress drive T₁ is here comparable to the radial advection term T₄, which explains the observed radial outward propagation of ZFs in figure 1.

Zoom In Zoom Out Reset image size

In the following the parametric dependence of each term T_i on the right hand side of equation (9) is investigated. To this end, the contribution of each term T_i on ZF evolution is measured by taking the L² norm, denoted by $\newcommand{\re}{{\rm Re}} \renewcommand{\ne}{{n}} \newcommand{\bi}{\boldsymbol} \big\| h\big\|_{2}$, of the time integrated contribution. Following this, we propose a measure of the relative ZF contribution

$$\begin{align} \newcommand{\bi}{\boldsymbol} \newcommand{\e}{{\rm e}} \renewcommand{\ne}{{n}} \newcommand{\re}{{\rm Re}} \displaystyle \qquad M_i:=\frac{\big\| \int_0^{t_1} {\rm d}t \hspace{0.5mm} T_i\big\|_{2} }{\sum_{j=1}^{5}\big\| \int_0^{t_1} {\rm d}t T_j\big\|_{2}}. \nonumber \end{align} \tag{ 13 }$$

In figure 3(a) we show that the relative contribution M_i of the NOB ZF terms ($T_2, \dots, T_5$) decreases with the reference background density gradient length L_n in the high collisionality regime ($\alpha=0.0005$). The summed up relative contribution of the NOB ZF terms exceeds the one of the conventional Reynolds stress for $\newcommand{\re}{{\rm Re}} \renewcommand{\ne}{{n}} \newcommand{\rhoN}{{\rho_{s0}}} L_n=32 \rhoN$. For steep reference background density gradients the radial advection term T₄ exhibits the largest relative contribution to the ZF dynamics of all the NOB terms.

Zoom In Zoom Out Reset image size

In figure 3(b) the dependence of the relative importance of each term on the adiabaticity parameter α is depicted for a fixed reference background density gradient length $\newcommand{\re}{{\rm Re}} \renewcommand{\ne}{{n}} \newcommand{\rhoN}{{\rho_{s0}}} L_n=32\rhoN$. While the conventional Reynolds stress term is again the dominating ZF contributor in particular for small collisionalities, all the NOB terms, except the background density gradient drive T₅, gain in importance for higher collisionalities. Interestingly, in the small collisionality regime ($\alpha\geqslant0.01$) the background density gradient drive T₅ exceeds all the remaining NOB actors. The quadruple fluctuation drive T₂ is for all studied parameters the smallest contributor to the ZF dynamics.

The dependence of the ZF terms on the time averaged maximum of the relative density fluctuation level $\newcommand{\avgR}[1]{{\langle #1 \rangle}} \newcommand{\flucR}[1]{{\widetilde{#1}}} \newcommand{\re}{{\rm Re}} \renewcommand{\ne}{{n}} \newcommand{\bi}{\boldsymbol} \langle \big\| \flucR{\ne}/\avgR{\ne} \big\|_\infty \rangle_t$ is shown in figure 4. Here, we denote the time average by $\langle h \rangle_t$ and compute the maximum with the help of the supremum norm $\newcommand{\re}{{\rm Re}} \renewcommand{\ne}{{n}} \newcommand{\bi}{\boldsymbol} \big\| h \big\|_\infty$. In figure 4 the conventional Reynolds stress drive T₁ contribution weakens with increasing relative density fluctuation level. The radial advection term T₄ and the triple fluctuation drive T₃ are the dominating NOB ZF contributors for high relative density fluctuations, while the background density gradient drive T₅ can be relevant likewise for small relative density fluctuations.

Zoom In Zoom Out Reset image size

We have generalized the ZF equation (3) to account for NOB effects in equation (9). Most importantly, the former Reynolds stress $\mathcal{R}$ is replaced by the Favre stress $\mathcal{F}$, which adds to its predecessor in case of high relative density fluctuations. The latter is accompanied by two new agents in the NOB ZF equation (9). The first of these radially advects ZFs by the Favre averaged radial drift velocity, which is proportional to the averaged radial particle flux. The second term scales inversely with the background density gradient length and affects the ZF dynamics even if the relative density fluctuations are small or if the Favre stress is radially homogeneous. Thus, this term may be of significance in or during the formation of radial transport barriers, where steep density profiles form with strongly reduced radial particle transport.

Additionally we extended the ordinary and modified HW model to the full-F theory. We simulated the full-F gyro-fluid model with the modified HW closure to numerically corroborate our theoretical results. The simulations successfully reproduced the predicted radial advection of ZFs, which appeared for small reference background density gradient lengths and large averaged radial particle flux. Moreover, our numerical parameter study showed that the NOB ZF drives can be comparable to the Reynolds stress drive in the herein scanned parameter range. In particular the deviation between the Reynolds and Favre stress drive increases with the relative density fluctuation amplitude, collisionality and inversely with the reference background density gradient length. This deviation is mainly reasoned in the triple fluctuation drive. Its importance in steep background density gradient regimes is in qualitative agreement with the theoretical estimate in the strong turbulence regime [38]. A similar dependence as for the Favre stress drive is found for the radial ZF advection mechanism. For the background density gradient drive only a dependence on the reference background density gradient is observed.

The presented results strongly argue in favor of the development and application of full-F gyro-fluid or gyro-kinetic models for simulation of fusion edge plasma turbulence, and in general demonstrate exemplarily the relevance of NOB effects for ZF formation in fluids and plasmas with large fluctuations and inhomogeneities. The latter conditions prevail e.g. during the low- to high-confinement mode transition. Thus, a consistent full-F simulation approach of this phenomenon is crucial to allow for the herein presented NOB ZF mechanisms.

Finally, we emphasize that the relative error between the Favre and Reynolds average of the poloidal velocity, derived to $\newcommand{\avgF}[1]{{\left[ \left[#1 \right]\right]}} \newcommand{\avgR}[1]{{\langle #1 \rangle}} \newcommand{\re}{{\rm Re}} \renewcommand{\ne}{{n}} \newcommand{\flucR}[1]{{\widetilde{#1}}} |\avgF{\flucR{u_y}}/\avgR{u_y}|$, is typically below a few percent. Thus, our proposed NOB ZF theory is also applicable to experimental measurements of the Reynolds averaged poloidal velocity $\newcommand{\re}{{\rm Re}} \renewcommand{\ne}{{n}} \newcommand{\avgR}[1]{{\langle #1 \rangle}} \avgR{u_y}$.

This work was supported by the Austrian Science Fund (FWF) Y398. R.K. was supported with financial subvention from the Research Council of Norway under grant 240510/F20. The computational results presented have been achieved using the Vienna Scientific Cluster (VSC) and the EUROfusion High Performance Computer (Marconi-Fusion). This work has been carried out within the framework of the EUROfusion Consortium and has received funding from the Euratom research and training programme 2014–2018 under grant agreement No 633053. The views and opinions expressed herein do not necessarily reflect those of the European Commission.