5.1 Reynolds stress as a proxy for the drive of zonal flows

Zonal flows are linearly stable and are driven by turbulence through the quadratic $E \times B$ non-linearity appearing in the gyrokinetic equation. To better understand their evolution, we analyse the properties of Reynolds stress (Diamond et al. Reference Diamond, Itoh, Itoh and Hahm2005), more exactly the off-diagonal component $\langle \tilde {V}_x\tilde {V}_\chi \rangle$ of the Reynolds stress tensor resulting from the combination of fluctuating $E\times B$ flow components in the radial and poloidal directions. We justify in the following that this Reynolds stress can be considered as a valid proxy of the zonal flow drive (Weikl et al. Reference Weikl, Peeters, Rath, Seiferling, Buchholz, Grosshauser and Strintzi2018). By analysing Reynolds stress we will be able to identify the different possible mechanisms driving zonal flows, their statistical properties and relative importance.

We start by considering an approximate evolution equation for the shearing rate $\omega _{E\times B}$ associated with the $E\times B$ zonal flows. As shown by Parra & Catto (Reference Parra and Catto2009) and Abiteboul et al. (Reference Abiteboul, Garbet, Grandgirard, Allfrey, Ghendrih, Latu, Sarazin and Strugarek2011); Abiteboul (Reference Abiteboul2012), such an equation can be obtained from the radial conservation equation for the total gyrocentre charge density, which, in turn, is derived by taking the appropriate velocity moment and flux-surface average of the gyrokinetic equation. The approach in Abiteboul (Reference Abiteboul2012) has been considered here, but starting from the gyrokinetic equation in the limit of the local (flux-tube) delta-f model rather than the global full-f model considered by Abiteboul (Reference Abiteboul2012). Furthermore, assuming the electrostatic limit, invoking $m_e/m_i\ll 1$, and making use of the quasi-neutrality equation in the long-wavelength approximation (final result thus valid only to second order in $k_\perp \rho _i$), leads to a relation that can be interpreted as a generalised vorticity equation:

(5.1)\begin{equation} \frac{\partial }{\partial t} (\varOmega+\varPi) = \frac{\partial^{2}}{\partial x^{2}} (\mathcal{R} + \mathcal{P}) + \frac{\partial}{\partial x}\mathcal{N}. \end{equation}

One identifies on the left-hand side of (5.1) the generalised vorticity term $\varOmega +\varPi$, composed of the actual vorticity associated with zonal flows and closely related to the zonal flow shearing rate $\omega _{E\times B}$ (by neglecting the $z$ dependence of $B_0$ and $g^{xx}\simeq 1$),

(5.2)\begin{equation} \varOmega = n_{0,i}m_i \frac{\partial^{2}}{\partial x^{2}} \left\langle \frac{g^{xx}}{B_0^{2}}\,\varPhi \right\rangle_{yz} \simeq \frac{n_{0, i}m_i}{B_0}\,\omega_{E\times B}, \end{equation}

as well as a perpendicular pressure term, related to lowest-order finite ion Larmor radius effects,

(5.3)\begin{equation} \varPi =\frac{m_i}{2q_i} \frac{\partial^{2}}{\partial x^{2}} \left \langle \frac{g^{xx}}{B_0^{2}}\, P_{\perp,i} \right \rangle_{yz}, \end{equation}

with the fluctuating part of the perpendicular ion gyrocentre pressure $P_{\perp ,i}$ expressed in terms of the corresponding gyrocentre distribution fluctuation $\delta f_i$,

(5.4)\begin{equation} P_{\perp,i} = \int \mu B_0\delta f_i\,\textrm{d}^{3}v. \end{equation}

On the right-hand side of (5.1) one identifies the second radial derivative of a term $\mathcal {R}=(n_{0,i}m_i/\mathcal {C})\textrm {RS}$, proportional to the Reynolds stress component $\textrm {RS}$ driving zonal flows:

(5.5)\begin{equation} \textrm{RS} = \left\langle \frac{1}{B_0^{2}}\frac{\partial \varPhi}{\partial y} \left( g^{xx}\frac{\partial \varPhi}{\partial x} + g^{xy}\frac{\partial \varPhi}{\partial y} \right) \right\rangle_{yz}. \end{equation}

This Reynolds stress term also has a finite Larmor radius correction term $\mathcal {P}$, again expressed in terms of $P_{\perp ,i}$:

(5.6)\begin{equation} \mathcal{P} = \frac{m_i}{2q_i\mathcal{C}} \left\langle \frac{1}{B_0^{2}} \left( g^{xx} \frac{\partial\varPhi}{\partial x} \frac{\partial P_{\perp,i}}{\partial y} + 2g^{xy}\, \frac{\partial \varPhi}{\partial y} \frac{\partial P_{\perp,i}}{\partial y} + g^{xx} \frac{\partial \varPhi}{\partial y} \frac{\partial P_{\perp,i}}{\partial x} \right) \right\rangle_{yz}. \end{equation}

The last contribution on the right-hand side of (5.1) is the radial derivative of the so-called neoclassical term $\mathcal {N}$, related to both curvature and $\boldsymbol {\nabla } B$ drifts:

(5.7)\begin{equation} \mathcal{N} = -\sum_\textrm{species} \frac{2{\rm \pi}\mathcal{C}}{m} \left\langle \gamma_2 \frac{\partial B_0}{\partial z} \int \textrm{d}v_\parallel \,\textrm{d}\mu\,\frac{mv_\parallel^{2} + \mu B_0}{B_0^{2}} \left( \delta f + \frac{q\,\bar{\varPhi}}{T_0} f_0 \right) \right\rangle_{yz}, \end{equation}

with $f_0$ the Maxwellian background distribution, $\gamma _2{=}g^{xx}g^{yz}-g^{xy}g^{xz}$ and the constant $\mathcal {C}{=}B_0/|\boldsymbol {\nabla } x\times \boldsymbol {\nabla } y|$.

The different terms appearing in (5.1) can be monitored as a diagnostic along a gyrokinetic simulation. For the purpose of verification, this diagnostic was first applied to simple test cases, including the Rosenbluth–Hinton problem (Rosenbluth & Hinton Reference Rosenbluth and Hinton1998) addressing the linear dynamics of zonal fluctuations (in this case, only the linear neoclassical term contributes on the right-hand side of (5.1)), as well as the non-linear decay of an initially single unstable eigenmode (as discussed in § 5.4). For these simple tests, only longer wavelength modes were considered, so as to stay well within the limit $k_\perp \rho _i\ll 1$ assumed for the derivation of relation (5.1).

After this successful initial verification phase, the terms appearing in (5.1) were monitored and compared for the fully developed turbulence simulations studied in this paper, considering the reference gradient case $R/L_{T,i}=6$. To validate Reynolds stress as a good proxy for the drive of zonal flows, the correlation between $\varOmega \sim \omega _{E\times B}$ and the two non-linear driving terms $\partial ^{2}\mathcal {R}/\partial x^{2}$ and $\partial ^{2}\mathcal {P}/\partial x^{2}$ in (5.1) was estimated over these simulations. To this end, the following correlation estimator between two observables $a$ and $b$, functions of the radial variable $x$ and time $t$, was applied:

(5.8)\begin{equation} \textrm{Corr}(a,b) = \frac{\sigma_{x,t}(a, b)}{\sigma_{x,t}(a)\sigma_{x,t}(b)} = \frac{ \langle (a-\langle a\rangle_{x,t}) (b-\langle b\rangle_{x,t}) \rangle_{x,t}} {\sqrt{\langle(a-\langle a\rangle_{x,t})^{2}\rangle_{x,t}} \sqrt{\langle(b-\langle b\rangle_{x,t})^{2}\rangle_{x,t}}}. \end{equation}

Significant correlation estimates were obtained in this way between $\varOmega$ and the Reynolds stress drive term $\partial ^{2}\mathcal {R}/\partial x^{2}$. Correlation values were found to increase further by considering only the longer wavelength contributions, achieved by filtering out $|k_x|\rho _i > 0.5$ and $|k_y|\rho _i>0.5$ Fourier modes from the signals. This is in agreement with the long-wavelength approximation assumed for deriving (5.1). The positive correlation values $\textrm {Corr} = 0.77$ and $0.37$ were obtained in this way between $\varOmega$ and $\partial ^{2}\mathcal {R}/\partial x^{2}$ for the adiabatic and kinetic electron turbulence simulations, respectively, whereas the correlation between $\varOmega$ and $\partial ^{2}\mathcal {P}/\partial x^{2}$ provided the values $\textrm {Corr} = 0.75$ and $0.40$, respectively. These results validate considering the Reynolds stress term $\partial ^{2}\textrm {RS}/\partial x^{2}$ as a proxy for the drive of zonal flow shear $\omega _{E\times B}$.

In the following, it will be insightful to consider the contributions from different $k_y$ Fourier modes components of the fluctuating fields to the Reynolds stress term $\textrm {RS}$. Relation (5.5) for $\textrm {RS}$ can indeed be written as a sum over $k_y$:

(5.9)\begin{equation} \textrm{RS}(x, t) = \sum_{k_y > 0} \widehat{\textrm{RS}}_{k_y}(x, t), \end{equation}

with the contribution from the $k_y$ Fourier mode $\hat \varPhi _{k_y}(x, z, t) = ({1}/{L_y}) \int _0^{L_y} \varPhi \exp (-{\mathrm i}k_yy)\,\textrm {d}y$ given by

(5.10)\begin{equation} \widehat{\textrm{RS}}_{k_y}(x, t) = 2\,{Re} \left[ \left\langle \frac{1}{B_0^{2}}\,k_y\hat{\varPhi}_{k_y} \left( g^{xx}\textrm{i}\,\frac{\partial \hat{\varPhi}_{k_y}^{\star}}{\partial x} + g^{xy} k_y \hat{\varPhi}_{k_y}^{\star} \right) \right\rangle_z \right], \end{equation}

having invoked the reality condition $\hat \varPhi _{-k_y}=\hat \varPhi _{k_y}^{\star }$. Considering as well the $k_x$ Fourier mode decomposition of $\varPhi$, each of these $k_y$ contributions can also be written as follows:

(5.11)\begin{equation} \widehat{\textrm{RS}}_{k_y}(x, t) = 2\,{Re} \left\{ \sum_{k_x, k_x''} \left\langle \frac{1}{B_0^{2}}\,k_y \left( g^{xx} k_x'' + g^{xy} k_y \right)\hat{\varPhi}_{k_x,k_y} \hat{\varPhi}_{k_x'',k_y}^{\star} \right\rangle_z \exp[\textrm{i}(k_x-k_x'')x] \right\}, \end{equation}

illustrating the drive of zonal modes $(k_x'=k_x-k_x'', 0)$ through non-linear interaction between Fourier modes $(k_x, k_y)$ and $(k_x'', k_y)$, extensively discussed in §§ 5.2 and 5.3.

For the study carried out in this paper, it is essential to understand how different drift modes may non-linearly interact to drive zonal flows via Reynolds stress. To start, we recall in the next subsection the basic mechanism underlying the drive of zonal flows in a simple shearless slab system before considering the more complex case of direct interest to us, i.e. of a sheared toroidal system.

5.2 Modulational instability in shearless slab system

The drive of zonal flows by microturbulence has been studied extensively in the literature considering a simple slab-like plasma confined by a uniform, shearless magnetic field. Such a system was in particular addressed in the original work by Hasegawa & Mima (Reference Hasegawa and Mima1978), where a cold fluid model was assumed for ions and an adiabatic response for electrons. In this model, choosing an orthogonal Cartesian coordinate system $(x, y, z)$, the magnetic field is aligned along $z$, $\boldsymbol {B} = B\boldsymbol {e}_z$ and the background density inhomogeneities along $x$, $\boldsymbol {\nabla } n_0 = (\textrm {d} n_0/\textrm {d}x)\boldsymbol {e}_x$. The corresponding well-known model equation for the non-linear evolution of the electrostatic potential $\varPhi$ associated to the fluctuating fields describes the essentially two-dimensional drift wave turbulence in the $(x, y)$ plane perpendicular to the magnetic field. This model equation led to a first understanding of the generation of zonal flows along $y$, i.e. in the direction both perpendicular to $\boldsymbol {B}$ and the direction of inhomogeneity. The emergence of such large-scale flows can in particular be explained as the result of an anisotropic inverse cascade of energy related to the conservation of energy and enstrophy in the two-dimensional turbulence. The emergence of zonal flows can also be understood at the level of elementary non-linear processes as recalled in the following.

In a shearless slab system, the spatial dependence of linear eigenmodes are given by a single Fourier mode $\varPhi (x, y)\sim \varPhi _{{{\boldsymbol {k}}}}\exp (\textrm {i}{{\boldsymbol {k}}}\boldsymbol {\cdot }{{\boldsymbol {x}}})$, with ${{\boldsymbol {x}}} = x\,{{\boldsymbol {e}}}_x+y\,{{\boldsymbol {e}}}_y$ and ${{\boldsymbol {k}}} = k_x{{\boldsymbol {e}}}_x+k_y{{\boldsymbol {e}}}_y$. The corresponding time dependence is of the form $\sim \exp (-\textrm {i}\omega _{{\boldsymbol {k}}}t)$, with $\omega _{{\boldsymbol {k}}}$ the eigenfrequency of the mode. In the simple Hasegawa–Mima model, one has $\omega _{{\boldsymbol {k}}} = k_y\,v_d/(1+k^{2})$, with ${{\boldsymbol {v}}}_d = -(T_e/eB)(\textrm {d}\log n_0/\textrm {d}x)\, {{\boldsymbol {e}}}_y$ these eigenmodes result from the quadratic non-linearity in the Hasegawa–Mima equation, related to the ${{\boldsymbol {v}}}_{E\times B} = (-\boldsymbol {\nabla }\varPhi \times {{\boldsymbol {B}}})/B^{2}$ drift. The elementary non-linear interaction thus involves a triplet of Fourier modes ${{\boldsymbol {k}}}$, ${{\boldsymbol {k}}}'$, ${{\boldsymbol {k}}}''$ verifying the wave vector matching condition ${{\boldsymbol {k}}} = {{\boldsymbol {k}}}'+{{\boldsymbol {k}}}''$, where each of the modes, e.g. ${{\boldsymbol {k}}}$, is coupled to the two others, ${{\boldsymbol {k}}}'$ and ${{\boldsymbol {k}}}''$ in this case. In case of frequency matching $\omega _{{\boldsymbol {k}}}\simeq \omega _{{{\boldsymbol {k}}}'}+\omega _{{{\boldsymbol {k}}}''}$ and under the condition $k' < k < k''$, one can have a resonant decay of mode ${{\boldsymbol {k}}}$, i.e. a transfer of energy from this mode, into the daughter modes ${{\boldsymbol {k}}}'$ and ${{\boldsymbol {k}}}''$. This basic process is referred to as the resonant three-wave interaction mechanism. One can furthermore show that in the case where mode ${{\boldsymbol {k}}}$ represents a drift wave, i.e. typically with $|k_y| \gg |k_x|$, the decay happens preferentially (meaning with a higher growth rate of decay) if one of the daughter wave vectors, e.g. ${{\boldsymbol {k}}}'$, is (nearly) aligned along the $x$ direction, $|k_x'| \gg |k_y'|$ (Hasegawa et al. Reference Hasegawa, Maclennan and Kodama1979). The $E\times B$ flow associated with the daughter mode ${{\boldsymbol {k}}}'$ is then obviously along $y$, thus explaining the emergence of zonal flows in this direction. Such a mode with vector aligned along the direction of inhomogeneity ${{\boldsymbol {e}}}_x$ is thus referred to as a zonal mode.

Let us still further consider the decay of a pump drift wave ${{\boldsymbol {k}}} = (k_x, k_y)$, $k_y\neq 0$, into a zonal mode ${{\boldsymbol {k}}}'= (k_x', 0)$ and the second daughter wave ${{\boldsymbol {k}}}'' = {{\boldsymbol {k}}}-{{\boldsymbol {k}}}' = (k_x-k_x', k_y)$, itself a drift wave. The non-linear interaction between the two drift waves ${{\boldsymbol {k}}}$ and ${{\boldsymbol {k}}}''$ thus provides the drive to the zonal mode ${{\boldsymbol {k}}}'$ via the Reynolds stress $\textrm {RS}$ discussed in § 5.1, whereas the non-linear coupling between the original drift wave ${{\boldsymbol {k}}}$ and zonal mode ${{\boldsymbol {k}}}'$, leading to the growth of the daughter drift wave ${{\boldsymbol {k}}}''$, actually represents the shearing of the drift mode ${{\boldsymbol {k}}}$ by the zonal flows associated with ${{\boldsymbol {k}}}'$.

Variations to the simple Hasegawa–Mima model have been considered in the literature. In particular, the enhanced Hasegawa–Mima model (Krommes & Kim Reference Krommes and Kim2000; Gallagher et al. Reference Gallagher, Hnat, Connaughton, Nazarenko and Rowlands2012) accounts for the fact that the adiabatic electron response is inhibited for magnetic surface-averaged fluctuations, i.e. modes ${{\boldsymbol {k}}}$ with $k_y=0$, in other words zonal modes, which leads them to having a reduced effective inertia compared with standard drift waves with $k_y\neq 0$. This effect results in an amplification of the decay rates of drift waves into zonal modes and, thus, to an enhancement of corresponding energy transfer.

A further refinement to the basic driving mechanism of zonal flows is obtained by accounting for the fact that given an initial large-amplitude drift mode $\boldsymbol {k}$, decaying into a zonal mode $\boldsymbol {k}'$, both triplet interactions $[\boldsymbol {k}, \boldsymbol {k}', \boldsymbol {k}'' = \boldsymbol {k}-\boldsymbol {k}']$ and $[\boldsymbol {k}, -\boldsymbol {k}', \boldsymbol {k}''' = \boldsymbol {k}+\boldsymbol {k}']$ may be simultaneously resonant, i.e. $\omega _{\boldsymbol {k}} \simeq \omega _{\boldsymbol {k}'}+\omega _{\boldsymbol {k}''}$ and $\omega _{\boldsymbol {k}} \simeq -\omega _{\boldsymbol {k}'}+\omega _{\boldsymbol {k}'''}$. This coupled pair of three-wave interactions leads to an effective four-wave interaction involving modes $\boldsymbol {k}$, $\boldsymbol {k}'$, $\boldsymbol {k}''$ and $\boldsymbol {k}'''$, referred to as the modulational instability mechanism (Gallagher et al. Reference Gallagher, Hnat, Connaughton, Nazarenko and Rowlands2012).

An important point to emphasise is that, in both the case of the simple resonant three-wave interaction or the more general modulational instability, resonant coupling between the initial large-amplitude (pump) drift mode ${{\boldsymbol {k}}}=(k_x, k_y\neq 0)$ and a zonal mode ${{\boldsymbol {k}}}'=(k_x', k_y'=0)$ is established via either one sideband ${{\boldsymbol {k}}}''={{\boldsymbol {k}}}-{{\boldsymbol {k}}}'$ or, respectively, two sidebands ${{\boldsymbol {k}}}''$ and ${{\boldsymbol {k}}}'''={{\boldsymbol {k}}}+{{\boldsymbol {k}}}'$, where all of these Fourier modes are linearly decoupled from each other. As a result, in the situation where in addition to the pump drift mode ${{\boldsymbol {k}}}$, the zonal mode ${{\boldsymbol {k}}}'$ itself already has a finite initial amplitude and a well-defined phase, the other daughter waves ${{\boldsymbol {k}}}''$ (and ${{\boldsymbol {k}}}'''$) are free to adapt their phases to ensure a resonant interaction and, thus, an efficient energy transfer from the pump to the zonal mode, resulting in the further amplification of this zonal mode. In fact, a finite-amplitude zonal mode can stimulate the decay of multiple non-zonal modes, thus leading to coherent (as a result of the frequency matching involved) and therefore correlated contributions from these non-zonal modes to the Reynolds stress drive of the zonal mode.

Under realistic conditions, multiple elementary mechanisms such as those discussed previously (resonant three-wave interaction/modulational instability) happen both simultaneously and successively in time, leading to an expanding set of Fourier modes, and ultimately to a fully developed turbulent spectrum.

5.3 Self-interaction mechanism in sheared toroidal system

As in the shearless slab geometry, in the case of a tokamak system, i.e. based on a sheared axisymmetric toroidal magnetic geometry, which is of main interest to our study, the modulational instability involving resonant three-wave interactions remains an essential driving mechanism of zonal flows (Chen et al. Reference Chen, Lin and White2000). In a tokamak, however, one must distinguish another form of the non-linear interaction leading to the drive of zonal modes. This mechanism, referred to as self-interaction, is specific to systems presenting magnetic shear and is explained in the following.

As already discussed in § 2, parallel boundary conditions lead to the linear coupling of $k_x$ Fourier modes. In particular, according to (2.16), the electrostatic potential field $\varPhi _{L}$ of a linear microinstability eigenmode with fixed $k_{x0}$ and $k_y\neq 0$ is composed of Fourier modes $\hat {\varPhi }_{k_{x0} + p2{\rm \pi} k_y\hat {s},k_y}$, $p\in \mathbb {Z}$. The spatial dependence of the corresponding eigenmode structure is thus given by

(5.12)\begin{equation} \varPhi_{L}(x, y, z) = \tilde{\varPhi}_{k_{x0},k_y}(x,z)\exp(\textrm{i}k_yy) + \textrm{c.c.}, \end{equation}

where $\textrm {c.c.}$ denotes the complex conjugate (considered here to ensure that $\varPhi _L$ is real-valued, essential for computing the quasi-linear estimate in (5.14)), and with the complex-valued $(x,z)$-dependent envelope given by

(5.13)\begin{equation} \tilde{\varPhi}_{k_{x0},k_y}(x,z) = \sum_{p=-\infty}^{+\infty} \hat{\varPhi}_{k_{x0} + p2{\rm \pi} k_y\hat{s},k_y}(z) \exp[\textrm{i}(k_{x0} + p2{\rm \pi} k_y\hat{s})x]. \end{equation}

These linearly coupled $k_x$ Fourier modes, all having same $k_y$, drive the zonal modes $(k_x', 0)$, with $k_x'=p'2{\rm \pi} k_y\hat {s}$, $p'\in \mathbb {Z}$, forming a set of harmonics. Indeed, any two Fourier modes $\hat {\varPhi }_{k_{x0} + p2{\rm \pi} k_y\hat {s}, k_y}$ and $\hat {\varPhi }_{k_{x0} + p''2{\rm \pi} k_y\hat {s},k_y}$ composing the physical eigenmode will drive, via three Fourier mode coupling, the zonal mode $\hat {\varPhi }_{p'2{\rm \pi} k_y\hat {s},0}$ with $p'=p-p''$ (see (5.11)). Note that this drive of zonal modes is via the same quadratic non-linearity in the gyrokinetic equation related to $E\times B$ drifts as that driving zonal modes through the modulational instability, but in this case involving Fourier modes that are already linearly coupled to each other.

Assuming that the relative phases between the Fourier modes $\hat {\varPhi }_{k_{x, 0} + p2{\rm \pi} k_y\hat {s},k_y}$ remains set by the linear coupling, even during the non-linear turbulent evolution (to what extent this assumption holds is validated in § 5.5), the phases of the associated contributions to Reynolds stress driving the zonal modes is fixed. This translates in direct space to an essentially fixed periodic radial dependence (with period ${\rm \Delta} x_\textrm {MRS} = 1/k_y\hat {s}$ corresponding to the distance between MRSs) of the contribution to Reynolds stress through this self-interaction mechanism from a given $k_y$ eigenmode, the overall magnitude of this contribution obviously varying in time as the amplitude of the eigenmode evolves. This is a critical difference compared with the drive of zonal modes ${{\boldsymbol {k}}}'$ through the modulational instability discussed in § 5.2, where the relative phases between the pump mode ${{\boldsymbol {k}}}$ and sidebands ${{\boldsymbol {k}}}''$, ${{\boldsymbol {k}}}'''$ are free to adapt to enable a coherent (i.e. in phase) drive of a given zonal mode.

Figure 9 plots the quasi-linear estimates of this self-interaction contribution to the Reynolds stress drive term $\partial ^{2}\textrm {RS}/\partial x^{2}$ from the linear eigenmode with $k_{x0}=0$, $k_{y}\rho _i = 0.28$ and $R/L_{T_i} =6$, whose ballooning structures are given in figure 2, considering both the case of kinetic and adiabatic electron response. These results are obtained by inserting the corresponding eigenmode structure (5.12) into relations (5.10) and (5.11), leading to the Reynolds stress contribution denoted by $\widetilde{\textrm {RS}}_{k_{x0},k_y}(x)$:

(5.14)\begin{align} \widetilde{\textrm{RS}}_{k_{x0},k_y}(x) & = \textrm{RS}[\varPhi_L] = \widehat{\textrm{RS}}_{k_y}[\tilde{\varPhi}_{k_{x0}, k_y}] \nonumber\\ & = 2\,{Re} \left[ \left\langle \frac{1}{B_0^{2}}\,k_y\tilde{\varPhi}_{k_{x0}, k_y} \left( g^{xx}\textrm{i}\,\frac{\partial \tilde{\varPhi}_{k_{x0}, k_y}^{\star}}{\partial x} + g^{xy} k_y \tilde{\varPhi}_{k_{x0}, k_y}^{\star} \right) \right\rangle_z\right]\nonumber\\ & = 2\,{Re} \left[ \sum_{p'=-\infty}^{+\infty} \exp({\textrm{i}p'2{\rm \pi} k_y\hat{s}x}) \sum_{p=-\infty}^{+\infty} \left\langle \frac{1}{B_0^{2}}\,k_y \left( g^{xx} k_x'' + g^{xy} k_y \right)\hat{\varPhi}_{k_x,k_y} \hat{\varPhi}_{k_x'',k_y}^{\star} \right\rangle_z\right]. \end{align}

In the last equality of relation (5.14), $k_x = k_{x0} + p2{\rm \pi} k_y\hat {s}$, $k_x'' = k_{x0} + p''2{\rm \pi} k_y\hat {s}$ and $p'' = p - p'$. This relation also clearly points out how this contribution to Reynolds stress from a given $k_y$ eigenmode through self-interaction is periodic with period ${\rm \Delta} x_\textrm {MRS} = 1/k_y\hat {s}$.

Figure 9.

Solid lines indicate the quasi-linear estimate of $\partial ^{2}\widetilde{\textrm {RS}}_{k_{x0},k_y}/\partial x^{2}$ normalised to its maximum value as a function of $x$ for the linear eigenmode with $k_{x0}=0$, $k_y\rho _i=0.28$ and either kinetic (blue) or adiabatic (red) electrons. Dashed lines denote the time-averaged $\langle \partial ^{2}\widehat{\textrm {RS}}^{\textrm {si}}_{k_y}/\partial x^{2}\rangle _t$ normalised with respect to the RMS in time of the total contribution $\partial ^{2}\widehat{\textrm{RS}}_{k_y}/\partial x^{2}$, for $k_y\rho _i=0.28$, in turbulence simulations considering either kinetic (blue) or adiabatic (red) electrons. Here $k_{y,\min }\rho _i=0.035$, $R/L_{T_i}=6$ and other parameters are given in table 1.

The overall amplitude of the quasi-linear estimates shown in figure 9 are naturally irrelevant. Furthermore, as these contributions to Reynolds stress have period ${\rm \Delta} x_\textrm {MRS} = 1/k_y\hat {s}$, only one such period is shown. Note how the radial profile of the quasi-linear estimate of $\partial ^{2}\widehat{\textrm{RS}}_{k_y}/\partial x^{2}$ is narrow in the case of kinetic electrons and localised around the MRS at $x=0$, clearly related to the fine structures of the eigenmode at MRSs and the associated broad tail in ballooning representation. As expected, the corresponding radial profile is much broader in the case of adiabatic electrons, as fine structures at MRSs are essentially absent in this case.

5.4 Evidence of zonal flow drive by modulational instability and self-interaction in reduced simulations

In this subsection, we consider reduced non-linear simulation setups in tokamak geometry to clearly illustrate the two basic mechanisms driving zonal flows discussed in §§ 5.2 and 5.3. For these reduced simulations, the same physical parameters as summarised in table 1 are considered, however with particular initial conditions defined as follows. An unstable linear eigenmode (hereby called the pump mode with $k_y=k_{y,\textrm {pump}}\neq 0$) is initialised with an amplitude a few times (${\sim }5$) less than the corresponding value in the fully saturated turbulence simulation discussed in § 5.5. The purpose of the simulation is to study how this single linear eigenmode drives zonal modes via either the modulational instability or self-interaction. The eigenmode with $k_{x0}=0$ and $k_{y}\rho _i=0.28$ was chosen for this pump mode, as it is among the most linearly unstable and also contributes significantly to the non-linear fluctuation spectra in the fully developed turbulence simulations (see figures 3 and 7b). In addition, zonal Fourier modes (with $k_y=0$) are initialised to amplitudes 10–12 orders of magnitude less than the pump mode to provide a necessary seed for possible modulational instabilities. All Fourier modes not part of the pump and zonal modes are initialised to zero. With this initial setup, only $k_y$ modes that are harmonics of $k_{y,\textrm {pump}}$ can possibly develop non-linearly ($k_y = pk_{y,\textrm {pump}}$, $p\in \mathbb {N}$). For these reduced simulations, we therefore set $k_{y, \min }=k_{y,\textrm {pump}}$ and actually only considered $N_{k_y}=8$ Fourier modes. The flux-tube width in the $x$ direction was set to $L_x=M/\hat {s}k_{y, \min }$ with $M=32$, so that it remains the same as in full turbulence simulations and the resulting fine $k_x$ spectrum allows for a detailed analysis of $k_x$-dependence of zonal mode growth. The remaining numerical resolutions are kept the same as in table 1. This system is then allowed to evolve until higher $k_y$ Fourier harmonics start to develop amplitudes similar to the fundamental $k_{y,\textrm {pump}}$. These steps ensure that the dominant non-linear interactions mainly involve only $k_y = 0$ and $k_{y,\textrm {pump}}$. This reduced non-linear setup thus clearly isolates the contribution to the drive of zonal modes from a single $k_y$ mode, while in the much more complex case of a standard fully developed turbulence simulation, multiple $k_y$ contributions may act simultaneously.

Figure 10 plots the evolution of the pump mode and the zonal modes driven by it. The solid black line represents the time trace of the $z$-averaged amplitude of the most dominant Fourier mode $\langle |\hat {\varPhi }_{k_{x0},k_{y,\textrm {pump}}}|\rangle _z(t)$ composing the pump mode. Note that the other linearly coupled Fourier modes composing the pump mode are not shown. Coloured lines represent the zonal mode amplitude $\langle |\hat {\varPhi }_{k_x,k_y=0}(t)|\rangle _z$ for each $k_x$. These results are shown for reduced non-linear simulations considering either adiabatic electrons (figure 10a) or kinetic electrons (figure 10b).

Figure 10.

Lin–log plots. Coloured lines represent the time evolution of the $z$-averaged amplitudes $\langle |\hat {\varPhi }_{k_x,k_y=0}|\rangle _z$ of zonal modes in reduced non-linear simulations initialised with a single large-amplitude $k_{y, \textrm {pump}}\neq 0$ eigenmode, considering either (a) adiabatic or (b) kinetic electrons. The colourbar maps the value of $k_x$ in units of $2{\rm \pi} k_{y, \textrm {pump}}\hat {s}$. Modes with $k_x=p'2{\rm \pi} k_{y, \textrm {pump}}\hat {s}$, $p'\in \mathbb {N}$, are plotted with dashed lines. The solid black line represents the time trace of the pump Fourier mode amplitude $\langle |\hat {\varPhi }_{k_x=0,k_{y,{\textrm {pump}}}}|\rangle _z$. The dotted black line represents an evolution proportional to $(\textrm {e}^{2\gamma t}-1)$, where $\gamma$ is the linear growth rate of the pump mode.

In the initial stage, the pump eigenmode grows exponentially with corresponding linear growth rates, i.e. $\gamma R/v_{\textrm {th},i}=0.25$ in the case of adiabatic electrons and $\gamma R/v_{\textrm {th},i} = 0.47$ in the case of kinetic electrons.

Zonal Fourier modes $\hat {\varPhi }_{k_x',k_y=0}$ with $k_x'=p'2{\rm \pi} k_{y,\textrm {pump}}\hat {s}$, $p'\in \mathbb {Z}$, (see dashed coloured time traces in figure 10) are driven by the large amplitude pump eigenmode through the self-interaction mechanism, i.e. are driven by multiple quadratic non-linearities, each involving two exponentially growing $k_x$ Fourier components of the pump, $\hat {\varPhi }_{k_x,k_{y,\textrm {pump}}}(t)$, $\hat {\varPhi }_{k_x'',k_{y,\textrm {pump}}}(t) \sim \textrm {e}^{\gamma t}$, where $k_x^{(\prime \prime )}=k_{x0}+p^{(\prime \prime )}2{\rm \pi} k_{y,\textrm {pump}}\hat {s}$, $p^{(\prime \prime )}\in \mathbb {Z}$, and $p-p''=p'$. One thus has for $k_x'=p'2{\rm \pi} k_{y,\textrm {pump}}\hat {s}$:

(5.15)\begin{align} \hat{\varPhi}_{k_x',k_y=0}(t) - \hat{\varPhi}_{k_x',k_y=0}(0) &\sim \int_0^{t} \hat{\varPhi}_{k_x,k_{y,\textrm{pump}}}(t')\, \hat{\varPhi}^{*}_{k_x'',k_{y,\textrm{pump}}}(t')\, \textrm{d}t' \sim \int_0^{t} \textrm{e}^{2\gamma t'}\,\textrm{d}t' \nonumber\\ &\sim (\textrm{e}\rm{}^{2\gamma t}-1). \end{align}

Note that for $t \ll 1/\gamma$, these modes thus start by growing linearly in time. In the initial stage of the simulations shown in figures 10(a) and 10(b), i.e. for $t \lesssim 28 R/v_{\textrm {th},i}$ and $t \lesssim 15 R/v_{\textrm {th},i}$ in the case of adiabatic and kinetic electrons, respectively, the Fourier modes driven by self-interaction, $\hat {\varPhi }_{p'2{\rm \pi} k_{y,\textrm {pump}}\hat {s},k_y=0}$, therefore dominate the zonal $k_x$ spectrum. Black dotted lines in these figures indicate the time evolution ${\sim }(\textrm {e}^{2\gamma t}-1)$, providing a good fit to the initial evolution of these particular zonal modes. As discussed in § 5.3, these modes are driven by the Reynolds stress resulting from the self-interaction of the linear eigenmode with $k_y = k_{y,\textrm {pump}}$, which in direct space is periodic in $x$ with period ${\rm \Delta} x_\textrm {MRS} = 1/k_{y,\textrm {pump}} \hat {s}$ and aligned with corresponding MRSs. This is clearly reflected in figures 11(a) and 11(b), plotting the effective zonal flow shearing rate $\omega _\textrm {eff}$ as a function of $x$ and time $t$, again for simulations with either adiabatic or kinetic electrons. At least in the initial stage of the simulations, these plots indeed present a periodic radial variation of $\omega _\textrm {eff}$ with period $1/k_{y,\textrm {pump}} \hat {s}$ and perfectly aligned with the MRSs of the pump mode.

Figure 11.

Effective zonal flow shearing rate $\omega _\textrm {eff}$ as a function of radial position $x$ and time $t$ for the same reduced non-linear simulations as in figure 10. Cases with either (a) adiabatic or (b) kinetic electrons are shown.

Zonal Fourier modes $\hat {\varPhi }_{k_x',k_y=0}$ with $k_x'\neq p'2{\rm \pi} k_{y,\textrm {pump}}\hat {s}$, $p'\in \mathbb {Z}$ (see solid coloured time traces in figure 10) are however driven by the pump eigenmode via the modulational instability mechanism, i.e. mainly through the quadratic non-linearities involving the large-amplitude Fourier component $\hat {\varPhi }_{k_{x0},k_{y,\textrm {pump}}}$ of the pump and either one of the initially low amplitude sideband modes $\hat {\varPhi }_{k_{x0}\pm k_x',k_{y,\textrm {pump}}}$. As a result, these particular zonal modes grow as ${\sim } \textrm {e}^{\gamma _\textrm {mod} t}$, where $\gamma _\textrm {mod}$ stands for the growth rate of the modulational instability. Given that $\gamma _\textrm {mod} \sim |\varPhi _{\textrm {pump}}|$ (Hasegawa et al. Reference Hasegawa, Maclennan and Kodama1979) and that the pump amplitude itself grows exponentially, $|\varPhi _{\textrm {pump}}| \sim \textrm {e}^{\gamma t}$, these zonal modes effectively end up growing super-exponentially.

One can also see in figure 10(b) that all plotted zonal modes driven by self-interaction, $\hat {\varPhi }_{p'2{\rm \pi} k_{y,\textrm {pump}}\hat {s},k_y=0}$, $p'=1, \ldots , 7$, are from the start driven up to similarly large amplitudes in the case of kinetic electrons, as opposed to the adiabatic case in figure 10(a), where these same modes are driven more weakly. This is explained by the different linear eigenmode structures of the pump in the adiabatic and kinetic electron cases (see the corresponding ballooning representations in figure 2). With kinetic electrons, the Fourier mode components $\hat {\varPhi }_{k_x,k_{y,\textrm {pump}}}(z)$, where $k_x=k_{x0}+ p2{\rm \pi} k_{y,\textrm {pump}}\hat {s}$ and $p\in \mathbb {Z}^{*}$ (set of integers excluding zero), which compose the tail of the ballooning representation, have much higher relative amplitudes compared with the main Fourier component $\hat {\varPhi }_{k_{x0},k_{y,\textrm {pump}}}(z)$ than in the case with adiabatic electrons. As a result, the self-interaction drive of a zonal mode $\hat {\varPhi }_{p'2{\rm \pi} k_{y,\textrm {pump}}\hat {s},k_y=0}$, $p' \in \mathbb {Z}^{*}$, dominated by the non-linear interaction between $\hat {\varPhi }_{k_{x0},k_{y,\textrm {pump}}}(z)$ and $\hat {\varPhi }_{k_{x0} + p'2{\rm \pi} k_{y,\textrm {pump}}\hat {s},k_{y,\textrm {pump}}}(z)$, is stronger with kinetic than with adiabatic electrons.

The weaker drive of the zonal modes by the self-interaction mechanism in the case with adiabatic electrons also explains why these particular zonal modes end up being overwhelmed by those driven by the modulational instability mechanism (reflected by the absence of fine structures at MRSs on $\omega _\textrm {eff}$ in figure 11(a) for $t \gtrsim 28R/v_{\textrm {th},i}$), whereas they remain significant in the case with kinetic electrons even after saturation of the modes driven by the modulational instability (reflected by the persistent fine structures on $\omega _\textrm {eff}$ in figure 11(b) even for $t \gtrsim 15R/v_{\textrm {th},i}$).

In the following two subsections, evidence of the self-interaction and modulational instability mechanisms are provided in the case of fully developed turbulence simulations, and their relative importance are compared between the cases of adiabatic and kinetic electron responses. Such simulations involve a fully saturated spectra of multiple $k_y$ modes than the case in the just considered reduced non-linear simulations. We find that the results obtained in the reduced non-linear simulations follow, in general, to the full turbulence scenarios, but with certain differences.

5.5 Evidence of self-interaction in turbulence simulations

In this subsection, the self-interaction mechanism is analysed in fully developed turbulence simulations with physical and numerical parameters given in table 1. We first test whether, in this non-linear system, the relative phase differences between the linearly coupled $k_x$ Fourier modes remain the same as in corresponding linear eigenmode, as was assumed earlier in § 5.3. It is found that the relative phase difference is only partially preserved, and the consequent effect on the Reynolds stress driving zonal flows via self-interaction is illustrated.

To find the non-linear modification of an eigenmode, first the absolute value of the time-averaged ballooning structure of the electrostatic potential $|\langle \hat {\varPhi }_{b,\textrm {nl}}(\chi ,t)/\hat {\varPhi }_{b,\textrm {nl}}(\chi _0=0,t)\rangle _{t}|$ for an eigenmode with $k_{x0}=0$ and $k_y\rho _i=0.28$ in the full turbulence simulation is plotted along with the corresponding linear ballooning structure $|\hat {\varPhi }_{b,\textrm {lin}}(\chi )/\hat {\varPhi }_{b,\textrm {lin}}(\chi _0=0)|$, in figure 12(a). See (2.18) for the definition of the ballooning structure $\hat {\varPhi }_{b}(\chi )$. The non-linear result with kinetic electrons starts significantly deviating from the linear result for $|\chi |\gtrsim {\rm \pi}$. That is, the ‘giant tails’ in the ballooning structure are reduced compared with the linear case. Nonetheless, the non-linear ballooning structure with kinetic electrons retain more significant tails compared with the case with adiabatic electrons. One also finds that, with kinetic electrons, the relative phase along the ballooning structure $\delta \phi _\textrm {nl}(\chi ,t)=\phi [\hat {\varPhi }_{b,\textrm {nl}}(\chi ,t)/\hat {\varPhi }_{b,\textrm {nl}}(\chi _0=0,t)]$ remains approximately constant in time (in contrast to adiabatic electron results) and equal to its linear value $\delta \phi _\textrm {lin}(\chi )=\phi [\hat {\varPhi }_{b,\textrm {lin}}(\chi )/\hat {\varPhi }_{b,\textrm {lin}}(\chi _0=0)]$; here, $\phi [A]=\textrm {arg}(A)$ denotes the phase or argument of the complex number $A$. This is true, in particular, for $|\chi |\lesssim 4{\rm \pi}$, as demonstrated by the solid blue line in figure 12(b) representing ${\rm \Delta} \phi (\chi ,t)=\delta \phi _\textrm {nl}(\chi ,t)-\delta \phi _\textrm {lin}(\chi )$ for $\chi =2{\rm \pi}$ in the kinetic electron case. Although there are jumps of $2{\rm \pi}$, the line closely adheres to the linear phase difference. To be quantitative, one considers the estimate $\textrm {MOD}({\rm \Delta} \phi (\chi ,t))$, where $\textrm {MOD}(A)=(\langle |\textrm {mod}_{2{\rm \pi} }(A)|^{2}\rangle _t)^{1/2}$, $\textrm {mod}_{2{\rm \pi} }(A) \equiv A-2{\rm \pi} \times \textrm {round}(A/2{\rm \pi} )$, $A\in \mathbb {R}$, and the function round provides the nearest integer. One notes that, lower the value of $\textrm {MOD}({\rm \Delta} \phi )$, more strongly is the relative phase fixed by the linear coupling. Further, for uniform random values of ${\rm \Delta} \phi$ between $-{\rm \pi}$ and ${\rm \pi}$, one obtains $\textrm {MOD}({\rm \Delta} \phi )=0.58{\rm \pi}$. It is found that, for the case with kinetic electrons, $\textrm {MOD}({\rm \Delta} \phi (\chi =2{\rm \pi} ,t))=0.34{\rm \pi}$. At $\chi =4{\rm \pi}$, as illustrated by the dotted blue line in figure 12(b), the phase imposed by linear coupling appears to be weakened but still present, giving $\textrm {MOD}({\rm \Delta} \phi (\chi ,t))=0.50{\rm \pi}$. The corresponding results with adiabatic electrons, denoted by the solid and dashed red lines for $\chi =2{\rm \pi}$ and $4{\rm \pi}$, respectively, clearly do not retain constant phase differences. The respective values of $\textrm {MOD}({\rm \Delta} \phi (\chi ,t))$ are $0.57{\rm \pi}$ and $0.61{\rm \pi}$.

Figure 12.

(a) Solid lines denote the absolute value of the time-averaged ballooning structure of the electrostatic potential $|\langle \hat {\varPhi }_{b,\textrm {nl}}(\chi ,t)/\hat {\varPhi }_{b,\textrm {nl}}(\chi _0=0,t)\rangle _{t}|$ for $k_{x0}=0$ and $k_y\rho _i=0.28$ in the non-linear simulation described in table 1, with $k_{y,\min }\rho _i=0.035$ and $R/L_{T,i}=6$. Dotted lines denote the corresponding ballooning representation in linear simulations (same as in figure 2). Inset: enlarged view near $\chi =0$. Red and blue colours represent adiabatic and kinetic electron simulations respectively. (b) Phase difference ${\rm \Delta} \phi (\chi ,t)=(\phi [(\hat {\varPhi }_{b,\textrm {nl}}(\chi ,t)/\hat {\varPhi }_{b,\textrm {nl}}(\chi _0=0,t))(\hat {\varPhi }_{b,\textrm {lin}}(\chi _0=0)/\hat {\varPhi }_{b,\textrm {lin}}(\chi ))])$ plotted as a function of time. Here, $\hat {\varPhi }_{b,\textrm {nl}}$ and $\hat {\varPhi }_{b,\textrm {lin}}$ denote the ballooning representation of the electrostatic potential for the same eigenmode considered in (a), in non-linear and linear simulations, respectively. Solid and dotted lines represent $\chi =2{\rm \pi}$ and $4{\rm \pi}$, respectively.

In other words, in a full turbulence simulation with adiabatic electrons, for a linearly unstable eigenmode with given $k_{x0}$ and $k_y$, linear coupling between the Fourier modes $(k_x=k_{x0}+2{\rm \pi} p k_y\hat {s},k_y)$ where $p\in \mathbb {Z}$ have become subdominant compared with non-linear coupling effects. Thus, each $(k_x,k_y)$ Fourier mode can be considered linearly decoupled as in a shearless slab system and, hence, modulational instability is the only driving mechanism of zonal flows in turbulent simulations with adiabatic electrons. On the other hand, with kinetic electrons, for a linearly unstable physical mode with given $k_{x0}$ and $k_y$, linear coupling of Fourier modes is preserved to some extent but reduced to only the three Fourier modes with $k_x=k_{x0}+2{\rm \pi} pk_y\hat {s}$ and $|p|=0,\pm 1$. Only these Fourier modes have their relative phases largely fixed by their linear dynamics (and not for all $p\in \mathbb {Z}$, as assumed in § 5.3). Hence, the self-interaction mechanism is effectively reduced to driving the zonal Fourier modes with $k'_x=\pm 2{\rm \pi} k_y\hat {s}$. Furthermore, the self-interacting contribution to Reynolds stress $\widetilde{\textrm {RS}}_{k_{x0},k_y}(x)$ (defined in (5.14)) from the particular physical eigenmode, becomes essentially sinusoidal in $x$, spanning the full distance between corresponding MRSs with a period of ${\rm \Delta} x_\textrm {MRS}=1/k_y\hat {s}$, and with fixed phase/sign in time. Note that the radial width of the fine structures do not have a dependence on electron–ion mass ratio in turbulence simulations as they are already broadened to span the distance between MRSs as a result of non-linear coupling effects. This is unlike the situation in linear simulations where such a dependence exists (Dominski et al. Reference Dominski, Brunner, Görler, Jenko, Told and Villard2015).

The net contribution $\textrm {RS}^\textrm {si}$ from the self-interaction mechanism to Reynolds stress can be obtained by summing $\widetilde{\textrm {RS}}_{k_{x0},k_y}(x)$ over all $k_y$'s and all ballooning angles (measured by $k_{x0}$; see (2.19)):

(5.16)\begin{equation} \textrm{RS}^\textrm{si}(x)=\sum_{k_y>0}\,\widehat{\textrm{RS}}^\textrm{si}_{k_y}(x), \end{equation}

where

(5.17)\begin{equation} \widehat{\textrm{RS}}^\textrm{si}_{k_y}=\sum_{k_{x0}=-{\rm \pi} k_y\hat{s}}^{{\rm \pi} k_y\hat{s}}\widetilde{\textrm{RS}}_{k_{x0},k_y}(x). \end{equation}

In figure 9, the self-interaction contribution $\partial ^{2}\textrm {RS}^\textrm {si}/\partial x^{2}$ from a particular $k_y$ in full turbulence simulation is plotted. The dashed lines denote the time average of $\partial ^{2}\widehat{\textrm{RS}}^\textrm {si}_{k_y}/\partial x^{2}$ for $k_y\rho _i=0.28$, normalised by the RMS in time of total $\partial ^{2}\widehat{\textrm{RS}}_{k_y}/\partial x^{2}$, in simulations with kinetic (blue) and adiabatic (red) electrons; the total Reynolds stress contribution $\widehat {\textrm {RS}}_{k_y}$ from a given $k_y$ being defined in (5.11). This diagnostic simultaneously measures the relative importance of self-interaction drive of zonal flows with respect to the total contribution to Reynolds stress drive from the considered $k_y$, as well as how good its phase/sign is fixed in time. One notes that, because for a given $k_y$, the eigenmode with $k_{x0}=0$ is the most unstable/dominant, similar results are obtained for $\widetilde{\textrm {RS}}_{k_{x0}=0,k_y}(x)$ as well, instead of $\widehat{\textrm{RS}}^\textrm {si}_{k_y}$ in this diagnostic. From figure 9, one can conclude that, in full turbulence simulations with kinetic electrons, each $k_y$ contribution $\partial ^{2}\widehat{\textrm{RS}}^\textrm {si}_{k_y}/\partial x^{2}$ from self-interaction is significant. Furthermore, it is a sinusoidal with a period of $1/k_y\hat {s}$ and maintains the same sign at each radial position (although with randomly varying amplitude in time, as will be shown at the end of § 5.6). In particular, $\partial ^{2}\widehat{\textrm{RS}}^\textrm {si}_{k_y}/\partial x^{2}$ is zero at MRSs, negative to the left and positive to the right. Whereas with adiabatic electrons, the self-interaction contribution is weak, does not have a fixed sign, and essentially averages out to zero over time. Further, with kinetic electrons, these significant time-averaged self-interaction contributions, localised at MRSs of each $k_y$, align at LMRSs, driving time-averaged zonal flows at these radial positions. However, away from LMRSs, the time-averaged self-interaction contributions $\partial ^{2}\widehat{\textrm{RS}}^\textrm {si}_{k_y}/\partial x^{2}$ from each $k_y$ are spatially misaligned and cancel each other out, resulting in relatively negligible time-averaged zonal flow levels.

In figure 13, the time average of the total $\partial ^{2}\widehat{\textrm{RS}}_{k_y}/\partial x^{2}$ is plotted as a function of both position $x$ and wave number $k_y$. Although between LMRSs, $\partial ^{2}\widehat{\textrm{RS}}_{k_y}/\partial x^{2}$ follows essentially the same spatial orientation (phase) as $\partial ^{2}\widehat{\textrm{RS}}^\textrm {si}_{k_y}/\partial x^{2}$ in figure 9, near LMRSs, there is a reversal of spatial phase. We suspect this to be a result of the back reaction of zonal modes driven by self-interaction on turbulence. At the moment, we postpone further analysis on this to a later time.

Figure 13.

(a) Each $k_y$ contribution of $\langle \partial ^{2}\textrm {RS}/\partial x^{2}\rangle _{t}$ plotted as a function of the radial coordinate $x$. (b) Total $\langle \partial ^{2}\textrm {RS}/\partial x^{2}\rangle _{t}$ plotted versus $x$. Dashed line denote the standard deviation. Results correspond to the non-linear simulation described in table 1, with $k_{y,\min }\rho _i=0.035$, $R/L_{T,i}=6$ and kinetic electrons.

Recall from earlier in this subsection that in kinetic electron turbulence simulations, the self-interaction contribution from a given $k_y\neq 0$ is effectively reduced to driving the zonal Fourier modes with $k'_x=\pm 2{\rm \pi} k_y\hat {s}$. The combined self-interacting contributions from the various $k_y$ therefore drive zonal Fourier modes with $k_x=2{\rm \pi} p\hat {s}k_{y, \min }$, where $p\in \mathbb {Z}$. This can be seen as peaks in the $k_x$ spectra of the time-averaged shearing rate $\langle \omega _\textrm {eff}\rangle _{t}$, shown with a blue dashed line in figure 14 and corresponding to the time-constant structures at LMRSs separated by ${\rm \Delta} x_{\textrm {LMRS}}=1/\hat {s}k_{y, \min }$ in figure 5. Given that, in the case of the reference turbulence simulation with kinetic electrons, the highest contribution to heat flux (see figure 7b) and the $|\varPhi |^{2}$ amplitude spectra is at $k_y\rho _i\sim 0.21$ and that, for a given $k_y$, eigenmodes with $k_{x0}=0$ are the most dominant, the zonal mode with $k_x\rho _i\sim 2{\rm \pi} \hat {s}\times 0.21\sim 1.06$ is driven strongly via the self-interaction mechanism and, hence, has the highest value in the $k_x$ spectra of shearing rate. Similarly, those $k_y$ with lower heat flux and amplitude contributions contribute less towards the self-interaction mechanism and this explains the number and position of peaks in the $k_x$ spectra of the shearing rate with kinetic electrons in figure 14.

Figure 14.

The $k_x$ spectra of the effective shearing rate $\omega _\textrm {eff}$ in turbulence simulations with adiabatic (red) and kinetic (blue) electrons. Corresponding parameters are given in table 1, with $k_{y, \min }\rho _i=0.035$ and $R/L_{T,i}=6$. Solid lines and dashed lines represent time average of the amplitude and amplitude of the time average of $\omega _\textrm {eff}$, respectively. All traces are normalised by corresponding maximum linear growth rate $\gamma _{\max }$. Vertical dotted lines indicate positions where $k_x=2{\rm \pi} p\hat {s}k_{y, \min }$; $p\in \mathbb {Z}$.

Solid red line in figure 14 denotes the $k_x$ spectra of $\omega _\textrm {eff}$ in the simulation with adiabatic electron response. These zonal Fourier modes are driven predominantly by modulational instability. In the next subsection, we test whether, consequently, as a result of the frequency matching involved in the modulational instability mechanism, those zonal modes having high shearing rate amplitudes end up playing a major role in pacing the dynamics of much of the $(k_x,k_y\neq 0)$ Fourier modes in adiabatic electron simulations.

5.6 Evidence of modulational instability in turbulence simulations: bicoherence and correlation analysis

In this subsection, we verify, in fully developed turbulence simulations, that indeed with adiabatic electron response zonal flows are driven predominantly via modulational instability, as found in reduced simulations in § 5.4 and as discussed in § 5.5. We test this by doing a bicoherence-like analysis. We simultaneously show that, as a consequence, this leads to correlated contributions from the various $k_y$'s to the Reynolds stress drive of zonal modes. We furthermore show that in kinetic electron simulations, the self-interaction contributions to Reynolds stress from the various $k_y$ are random in time and decorrelated with each other.

As discussed earlier, modulational instability mechanism involves resonant three-wave interactions that require frequency matching. The strength of a particular resonant interaction between the three Fourier modes, ${\boldsymbol {k}}=(k_x,k_y)$, the zonal mode ${\boldsymbol {k}'}=(k'_{x},0)$ and daughter mode ${\boldsymbol {k}''}={\boldsymbol {k}}-{\boldsymbol {k}'}=(k_x-k'_x,k_y)$, can thus be measured via a bicoherence-type analysis, essentially involving the time average of the triplet product

(5.18)\begin{equation} T({\boldsymbol{k}};+{\boldsymbol{k}'})=\hat{\varPhi}_{\boldsymbol{k}}(t) \hat{\varPhi}^{*}_{\boldsymbol{k}'}(t)\hat{\varPhi}^{*}_{\boldsymbol{k}''}(t), \end{equation}

where $\hat {\varPhi }_{\boldsymbol {q}} (t) \sim \exp [-i(\omega _{\boldsymbol {q}} t + \phi _{\boldsymbol {q}})]$ is the complex time-dependent amplitude of the electrostatic field with any Fourier mode ${\boldsymbol {q}}$, having a frequency $\omega _{\boldsymbol {q}}$, initial phase $\phi _{\boldsymbol {q}}$ and evaluated at $z=0$. If the Fourier modes $[\boldsymbol {k}, \boldsymbol {k}', \boldsymbol {k}'' = \boldsymbol {k}-\boldsymbol {k}']$ are frequency matched in the simulation time, i.e. $\omega _{\boldsymbol {k}} = \omega _{\boldsymbol {k}'} + \omega _{\boldsymbol {k}''}$, then $\langle T({\boldsymbol {k}};+{\boldsymbol {k}'})\rangle _t\neq 0$, where $\langle \cdot \rangle _t$ denotes the time average over the simulation time.

Here, we remark that, for a given set of Fourier modes $[\boldsymbol {k}, \boldsymbol {k}', \boldsymbol {k}'' = \boldsymbol {k}-\boldsymbol {k}']$ under frequency matching condition, different resonant decay mechanisms could be possible, each with a specific phase difference ${\rm \Delta} \phi =\phi _{\boldsymbol {k}}-\phi _{\boldsymbol {k}'}-\phi _{\boldsymbol {k}''}$ associated with it. Hence, a non-zero time-averaged $\langle T({\boldsymbol {k}} ;+{\boldsymbol {k}'})\rangle _t$ over the simulation time is indicative of a particular resonant decay mechanism being persistent throughout the simulation.

Now, a normalised measure of the strength of the resonant three-wave interaction can be calculated by the following estimate:

(5.19)\begin{equation} b_{N}({\boldsymbol{k}} ;+{\boldsymbol{k}'}) = \frac{|\langle T({\boldsymbol{k}} ;+{\boldsymbol{k}'})\rangle_t|}{\langle |T({\boldsymbol{k}} ;+{\boldsymbol{k}'})|\rangle_t}, \end{equation}

defined as the bicoherence between the Fourier triplet $[\boldsymbol {k}, \boldsymbol {k}', \boldsymbol {k}'' = \boldsymbol {k}-\boldsymbol {k}']$. Note $0\leq b_{N}\leq 1$. Further, $b_{N}({\boldsymbol {k}} ;+{\boldsymbol {k}'})\simeq 1$ indicates a fully resonant three-wave interaction between ${\boldsymbol {k}}$, ${\boldsymbol {k}''}$ and the zonal mode ${\boldsymbol {k}'}$, whereas $b_{N}\simeq 0$ indicates a non-resonant process. As modulational instability is a simultaneous resonant interaction between both triplets $[\boldsymbol {k}, \boldsymbol {k}', \boldsymbol {k}'' = \boldsymbol {k}-\boldsymbol {k}']$ and $[\boldsymbol {k}, -\boldsymbol {k}', \boldsymbol {k}''' = \boldsymbol {k}+\boldsymbol {k}']$, we define total bicoherence

(5.20)\begin{equation} B_{N}({\boldsymbol{k}} ;{\boldsymbol{k}'}) = (b_N({\boldsymbol{k}} ;+{\boldsymbol{k}'})+ b_N({\boldsymbol{k}} ;-{\boldsymbol{k}'}))/2, \end{equation}

such that a value of $B_{N}\simeq 1$ indicates a fully resonant modulational instability mechanism. In general, values of $B_N$ closer to zero are indicative of non-resonant interactions whereas larger values of $B_{N}$ are characteristic of modulational instability mechanisms.

Figure 15(a) shows $B_N({\boldsymbol {k}} ;{\boldsymbol {k}'})$ as a function of ${\boldsymbol {k}}=(k_x,k_y)$, in the reference simulation with adiabatic electrons, for the fixed zonal mode with ${\boldsymbol {k}'}=(k'_{x}=0.13\rho _i^{-1},0)$ which has the highest contribution to the zonal shearing rate $k_x$ spectra (see the solid red line in figure 14). One can see that most ${\boldsymbol {k}}=(k_x,k_y)$ have high values of $B_N\gtrsim 0.3$, reflecting the dominance of resonant three-wave interaction processes. However, in figure 15(b), the corresponding kinetic electron simulation result for the zonal mode with ${\boldsymbol {k}'}=(k'_{x}=0.26\rho _i^{-1},0)$ having highest contribution to the shearing rate $k_x$ spectra shows much smaller values of $B_N$, indicating much weaker resonant three-wave interactions. Similar differences are seen between adiabatic and kinetic electron results for other zonal modes ${\boldsymbol {k}'}=(k'_{x},0)$ as well.

Figure 15.

The bicoherence level $B_N$, as defined by (5.20), in (a) adiabatic and (b) kinetic electron turbulence simulations, for the zonal modes with $k'_{x}\rho _i=0.13$ and $0.26$, respectively. Results correspond to simulations with parameters as given in table 1, with $k_{y, \min }\rho _i=0.035$ and $R/LT_i=6$.

Note that, for a given zonal mode ${\boldsymbol {k}'}=(k_x',0)$, under wave vector matching condition, one has the following approximate estimate between associated triple products and the Reynolds stress contributions driving the particular zonal mode:

(5.21)\begin{equation} \sum_{k_x}T({\boldsymbol{k}} ;{\boldsymbol{k}'}) \sim \hat{\varPhi}_{\boldsymbol{k}'}^{*}\widehat{\textrm{RS}}_{k_y}, \end{equation}

where $\widehat{\textrm{RS}}_{k_y}$ has been defined in (5.11). Now, if $\langle T({\boldsymbol {k}} ;{\boldsymbol {k}'})\rangle _t$ (and the corresponding total bicoherence $B_N({\boldsymbol {k}} ;{\boldsymbol {k}'})$) for multiple Fourier modes ${\boldsymbol {k}}=(k_x,k_y)$ are simultaneously significant for the same dominant zonal mode ${\boldsymbol {k}'}$ (as is the case in figure 15(a) for adiabatic electron simulation), then based on (5.21) one can conclude that the Reynolds stress contributions $\widehat{\textrm{RS}}_{k_y}$ from the various $k_y$ tend to be in phase with the particular zonal mode. This implies a coherent and, thus, correlated contributions (drives) from these various $\widehat{\textrm{RS}}_{k_y}$ ($\partial ^{2}\widehat{\textrm{RS}}_{k_{y}}/\partial x^{2}$). To verify this explicitly, we define an effective correlation function $C_{\textrm {RS}}$ measuring the average correlation between all pairs of $[\partial ^{2}\widehat{\textrm{RS}}_{k_{y,1}}/\partial x^{2}$, $\partial ^{2}\widehat{\textrm{RS}}_{k_{y,2}}/\partial x^{2}]$ for $k_{y,1}\neq k_{y,2}$:

(5.22)\begin{equation} {C_{\textrm{RS}}}[\,f]= \left.\sum_{\substack{k_{y,i}, k_{y,j} \\ k_{y,j}>k_{y,i}} } \,\frac{\textrm{Cov}[\,\hat{f}_{k_{y,i}},\hat{f}_{k_{y,j}}]}{\sigma[\, \hat{f}_{k_{y,i}}]\sigma[\,\hat{f}_{k_{y,j}}]} \right/ \sum_{\substack{k_{y,i}, k_{y,j} \\ k_{y,j}>k_{y,i}}} 1 . \end{equation}

Here $f=\partial ^{2}\textrm {RS}/\partial x^{2}$, $\hat {f}_{k_{y}}$= $\partial ^{2}\widehat{\textrm{RS}}_{k_y}(x)/\partial x^{2}$, covariance $\textrm {Cov}[a,b]=(\sigma ^{2}[a+b]-\sigma ^{2}[a]-\sigma ^{2}[b])/2$ and variance $\sigma ^{2}[a]=\langle |a-\langle a\rangle _t|^{2}\rangle _t$, with $\langle \cdot \rangle _t$ representing average over simulation time. Note that ${C_{\textrm {RS}}}\in [0,1]$, with 1 corresponding to perfect correlation between Reynolds stress drive from all $k_y$ and 0 corresponding to total decorrelation between them. In figure 16, $C_\textrm {RS}$ is plotted as a function of the radial position $x$, for adiabatic and kinetic electron turbulence simulations. We can find that ${C_{\textrm {RS}}}[\partial ^{2} \textrm {RS}/\partial x^{2}]$ for adiabatic electron case (solid red line) only shows ${\sim }7\,\%$ correlation. However,the important point to note is that this is roughly an order of magnitude higher than that for kinetic electron simulation (solid blue line). Thus, this diagnostic verifies that in adiabatic electron turbulence simulations where modulational instability mechanism is stronger, Reynolds stress drive $\partial ^{2}\widehat{\textrm{RS}}_{k_{y}}/\partial x^{2}$ from the various $k_y$ are indeed more correlated with each other, when compared with kinetic electron simulations where modulational instability mechanism is weaker.

Figure 16.

Correlation between the $k_y$ modes of $\partial ^{2}\textrm {RS}/\partial x^{2}$, as defined by (5.22), as a function of $x$ in turbulence simulations with kinetic (solid blue line) and adiabatic (solid red line) electrons. Corresponding parameters are given in table 1, with $k_{y, \min }\rho _i=0.035$ and $R/L_{T,i}=6$. Correlation between the $k_y$ modes of self-interacting contribution $\partial ^{2}\textrm {RS}^\textrm {si}/\partial x^{2}$ in the simulation with kinetic electrons is shown with a dashed blue line.

The same correlation diagnostic is now performed on the self-interacting part of Reynolds stress in kinetic electron simulations. In figure 16, the correlation between the $k_y$ contributions to $\partial ^{2}\textrm {RS}^\textrm {si}/\partial x^{2}$, measured by ${C_{\textrm {RS}}}[\partial ^{2} \textrm {RS}^{\textrm {si}}/\partial x^{2}]$, is found to be nearly zero. This implies that each $k_y$ contribution to $\partial ^{2}\textrm {RS}^\textrm {si}/\partial x^{2}$ at a given radial position (while, in general, having a fixed sign) has an amplitude that is completely independent and decorrelated with each other. The self-interaction drive of zonal modes from each $k_y$ thus acts like random kicks. Consequently, if the decorrelated (incoherent) contributions from the self-interaction mechanism becomes significant, it may disrupt the (coherent) modulational instability interactions. We interpret this as being the cause for lower bicoherence values of $B_N$ in figure 15(b), as well as the lower level of $\rm C_{RS}[\partial ^{2} \textrm {RS}/\partial x^{2}]$ in figure 16 with kinetic electrons as compared with adiabatic electrons.

To summarise, in this section, the two different mechanisms driving zonal flows, namely modulational instability and self-interaction, have been discussed in detail. We have shown that, as a result of the stronger linear $k_x$ couplings within an eigenmode in kinetic electron simulations, the associated drive of zonal flows via self-interaction is much more significant than in the case of adiabatic electron simulations. This has been illustrated in both reduced and full turbulence simulations. Furthermore, in turbulence simulations, it was shown that the contributions to zonal flow drive from each $k_y$ mode, via the self-interaction mechanism, are random and decorrelated with each other in time. In the next section, we explain how such a decorrelated drive can lead to a decrease in the fluctuating zonal flow levels with decreasing $k_{y,\min }\rho _i$.