Nonlinear growth of zonal flows by secondary instability in general magnetic geometry

TL;DR

This paper develops a theory for the nonlinear growth of zonal flows in magnetized plasma turbulence via secondary instability, applicable to tokamaks and stellarators, and validates it against gyrokinetic simulations.

Contribution

It introduces a general magnetic geometry theory for zonal flow growth, extending previous models and confirming predictions with simulations.

Findings

Predicted growth rates match gyrokinetic simulations.

Theory applicable to both tokamaks and stellarators.

Error scales with the small parameter of the theory.

Abstract

We present a theory of the nonlinear growth of zonal flows in magnetized plasma turbulence, by the mechanism of secondary instability. The theory is derived for general magnetic geometry, and is thus applicable to both tokamaks and stellarators. The predicted growth rate is shown to compare favorably with nonlinear gyrokinetic simulations, with the error scaling as expected with the small parameter of the theory.

Figures

Equations

$$\frac{\partial h}{\partial t}+v_{\parallel }\frac{\partial }{\partial z}h+({\mathbf{v}}_E+{\mathbf{v}}_d)\cdot {\nabla }_{\perp }h+{\mathbf{v}}_E\cdot \nabla F_0=\frac{q}{T_i}\frac{\partial {⟨\varphi ⟩}_{\mathbf{R}}}{\partial t}{F}_0,$$

$$\left\langle A({\bf r})\right\rangle_{\bf R}=\frac{1}{2\pi}\int_{0}^{2\pi}d\vartheta A({\bf R}+\mbox{\boldmath$\rho$}(\vartheta)).$$

$$\frac{1}{n_0}\int d^3\mathbf{v}{⟨h⟩}_{\mathbf{r}}=\frac{q\varphi }{T_i}+\frac{e\tilde{\varphi }}{T_e}.$$

$$\left\langle A({\bf R})\right\rangle_{\bf r}=\frac{1}{2\pi}\int_{0}^{2\pi}d\vartheta A({\bf r}-\mbox{\boldmath$\rho$}(\vartheta)).$$

$$\varphi =\sum _S{\hat{\varphi }}_S\exp (iS),$$

$$\varphi =\sum _{k_{\psi },k_{\alpha }}{\hat{\varphi }}_{\mathbf{k}}\exp (i{k}_{\psi }\psi +i{k}_{\alpha }\alpha ).$$

$$\overline{\varphi }$$

$$=\frac{1}{V^{\prime }}\sum _{k_{\psi },k_{\alpha }}{\int }_0^{2\pi }d\alpha \int \frac{dz}{B}\exp (i{k}_{\psi }\psi +i{k}_{\alpha }\alpha ){\hat{\varphi }}_{\mathbf{k}}$$

$$=\frac{1}{V^{\prime }}\sum _{k_{\psi },k_{\alpha }=0}{\int }_0^{2\pi }d\alpha \int \frac{dz}{B}\exp (i{k}_{\psi }\psi ){\hat{\varphi }}_{\mathbf{k}},$$

$$\frac{\partial \tilde{\varphi }}{\partial t}+\{\overline{\varphi },\tilde{\varphi }\}=0,$$

$$\frac{\partial }{\partial t}(-\overline{{\rho }^2{\nabla }_{\perp }^2\varphi })+\overline{\{\varphi ,-{\rho }^2{\nabla }_{\perp }^2\varphi -{\rho }^2{\nabla }_{\perp }^2\chi \}}$$

$$+\overline{\{{\rho }^2{\nabla }_{\perp }^2\varphi ,\chi \}}+\overline{{\rho }^2{\nabla }_{\perp }^2\{\varphi ,\chi \}}=0,$$

$$\frac{\partial \chi }{\partial t}+\{\varphi ,\chi \}=0,$$

$$\rho =\frac{m{v}_{\mathrm{T}}}{\sqrt{2}qB(\psi ,\alpha ,z)},$$

$$\chi (\mathbf{r})=\frac{1}{q{n}_0}\int d^3\mathbf{v}\left(\frac{m{v}_{\perp }^2}{4}\right)\left(h(\mathbf{r})-\frac{q\varphi (\mathbf{r})}{T_i}{f}_M\right).$$

$$\{f,g\}=\frac{1}{B}(\hat{\mathbf{b}}\times \nabla f)\cdot \nabla g.$$

$${\varphi }_p={\varphi }_{p0}\exp (i{S}_p)+c.c.,$$

$${\chi }_p={\chi }_{p0}\exp (i{S}_p)+c.c.,$$

$$\frac{\partial {\tilde{\varphi }}_s}{\partial t}+\{{\overline{\varphi }}_s,{\tilde{\varphi }}_p\}=0,$$

$$\frac{\partial }{\partial t}(-{\overline{{\rho }^2{\nabla }_{\perp }^2\varphi }}_s)+\overline{\{{\varphi }_p,-{\rho }^2{\nabla }_{\perp }^2{\varphi }_s-{\rho }^2{\nabla }_{\perp }^2{\chi }_s\}}+\overline{\{{\rho }^2{\nabla }_{\perp }^2{\varphi }_p,{\chi }_s\}}$$

$$+\overline{\{{\varphi }_s,-{\rho }^2{\nabla }_{\perp }^2{\varphi }_p-{\rho }^2{\nabla }_{\perp }^2{\chi }_p\}}+\overline{{\rho }^2{\nabla }_{\perp }^2\{{\varphi }_p,{\chi }_s\}}$$

$$+\overline{\{{\rho }^2{\nabla }_{\perp }^2{\varphi }_s,{\chi }_p\}}+\overline{{\rho }^2{\nabla }_{\perp }^2\{{\varphi }_s,{\chi }_p\}}=0,$$

$$\frac{\partial {\chi }_s}{\partial t}+\{{\varphi }_p,{\chi }_s\}+\{{\varphi }_s,{\chi }_p\}=0.$$

$${\varphi }_s=\exp (i\bar{S}-i\omega t)\sum _n{c}_n\exp (in{S}_p),$$

$${\chi }_s=\exp (i\bar{S}-i\omega t)\sum _n{d}_n\exp (in{S}_p),$$

$${\overline{\varphi }}_s={\overline{c}}_0\exp (i\bar{S}-i\omega t),$$

$${\tilde{\varphi }}_s=\exp (i\bar{S}-i\omega t)\left[{\tilde{c}}_0+\sum _{n\ne 0}{c}_n\exp (in{S}_p)\right],$$

$${\tilde{c}}_0=0,$$

$$c_1=i\frac{{\Omega }_{\varphi }}{\omega }j{\overline{c}}_0,$$

$$c_{-1}=-i\frac{{\Omega }_{\varphi }}{\omega }{j}^{\ast }{\overline{c}}_0,$$

$$c_n=0,|n|\ge 2,$$

$${\Omega }_{\varphi }\equiv \frac{(\hat{\mathbf{b}}\times \nabla \bar{S})\cdot \nabla S_p}{B}|{\varphi }_{p0}|=\frac{\partial S_p}{\partial \alpha }\frac{\partial \bar{S}}{\partial \psi }|{\varphi }_{p0}|,$$

$$\overline{-{\rho }^2{\nabla }_{\perp }^2{\varphi }_s}$$

$$b_{\psi }={\rho }^2{\left(\frac{\partial \bar{S}}{\partial \psi }\right)}^2|\nabla \psi {|}^2.$$

$$\omega {\bar{c}}_0\overline{b_{\psi }}=$$

$$b={\rho }^{2}B\frac{\partial S_p}{\partial \alpha }\frac{\partial \bar{S}}{\partial \psi }.$$

$${\omega }^2\overline{b_{\psi }}=-2\overline{{\Omega }_{\varphi }^2{b}_{\psi }\left(1+\frac{{\chi }_{p0}}{{\varphi }_{p0}}+\frac{{\chi }_{p0}^{\ast }}{{\varphi }_{p0}^{\ast }}\right)},$$

$$-i\zeta d_0=\left(c_{-1}{r}_{+}-c_1{r}_{-}+d_1{j}_{-}-d_{-1}{j}_{+}\right),$$

Full text

1 Introduction

The present work concerns the dynamics of zonal flows in magnetic confinement fusion devices, and in particular the three dimensional magnetic geometries of stellarators. The magnetic field lines in these devices trace out nested topologically toroidal surfaces, and zonal flows are special linearly stable modes, extending across these surfaces, that suppress turbulence and thereby improve confinement [1]. In the absence of a nonlinear turbulent drive, they oscillate and damp, tending asymptotically to a non-zero “residual” level. This linear behavior has already been studied in the stellarator context [2, 3, 4, 5]. The nonlinear mechanism of growth, however, is less understood. This growth can be estimated by “secondary” instability, and an analogous nonlinear decay mechanism is given by the “tertiary” instability. The goal of the present paper is to clarify this growth mechanism, and a future paper is planned to investigate nonlinear decay.

In the theory of plasma turbulence, the notion of secondary instability, generally speaking, provides a useful tool for understanding the nonlinear processes that determine fluctuation amplitudes. It is essentially an instability of an instability. That is, when a linearly unstable plasma wave, called a primary instability, grows to sufficiently large amplitude, nonlinear processes arise that can drive a new mode. The new mode, called a secondary instability, must eventually grow faster than the primary mode, since it has a growth rate proportional to the amplitude of the primary mode, and thus could ultimately cause saturation of the primary mode, by drawing energy away from it. On the other hand, the secondary mode can also function more indirectly, by driving otherwise stable modes called zonal flows, that are involved in mode saturation by “back-reaction” or shearing. This second scenario is of great importance in ion temperature gradient (ITG) turbulence in magnetic fusion devices, which is the motivation for the present work.

The secondary instability that drives zonal flows has been derived using fluid theory in two dimensions [6, 7], and related instabilities have also been calculated that assume scale separation between the zonal flow and primary mode (e.g. [8]). Some works on zonal flow growth have modeled the effect of magnetic geometry [9, 10], by deriving local fluid approximations of linear theory that, however, ultimately do not account for the variation of the linear modes within a flux surface. In the present work, we treat the ITG mode as a given parameter of the theory, so that its non-local coupling to the zonal mode may be fully calculated. It is demonstrated here that this non-local coupling strongly affects zonal flow growth. The derived theory is applicable to axisymmetric geometries (tokamaks) and also the fully three-dimensional geometry of stellarators. Although the theoretical zonal flow growth rate could be computed fully analytically by using an analytical model of the ITG mode, an exact result is obtained by instead numerically simulating the ITG mode and using the result as an input of the theory; we will take this second approach here.

The outline of the paper is as follows. In Sec. 2 we introduce the equations and notation, and the general representation needed for three dimensional magnetic geometry. In Sec. 3 we present the complete details of the derivation of the secondary mode growth rate. By directly calculating the effect of the infinite hierarchy of coupled Fourier modes, we demonstrate its equivalence to the truncated “four-wave” approximation, agreeing with expectations from previous work [7]. Because the derivation is asymptotically exact, it should agree with gyrokinetic simulations in the limit that the scale of the modes exceeds the ion Larmor radius. Therefore, in Sec. 3.1 we test this, and find that the difference between the theoretical and numerical growth rates indeed scales with the small parameter of the theory. Although the theory is valid for stellarator geometries, we choose for the sake of simplicity to make this first numerical demonstration in circular tokamak geometry. In Sec. 4 we summarize and make some final remarks.

2 Equations

Let us consider the nonlinear electrostatic ion gyrokinetic equation,

[TABLE]

where $h({\bf R}_{\perp},z,\varepsilon,\mu)$ is the non-adiabatic part of the ion distribution function, which depends on the gyro-center coordinate ${\bf R}_{\perp}$ and the independent velocity coordinates are $\mu=m_{i}v_{\perp}^{2}/(2B)$ and $\varepsilon=m_{i}v^{2}/2=m_{i}(v_{\parallel}^{2}+v_{\perp}^{2})/2$ where $m_{i}$ is the ion mass. We employ the following mostly standard conventions: ${\bf v}_{E}=\hat{\bf b}\times\boldsymbol{\nabla}\left\langle\phi\right\rangle_{\bf R}/B$ is the gyro-averaged $E\times B$ velocity; ${\bf v}_{d}=\Omega_{c}^{-1}\hat{\bf b}\times(v_{\parallel}^{2}\hat{\bf b}\cdot\boldsymbol{\nabla}\hat{\bf b}+v_{\perp}^{2}/2\boldsymbol{\nabla}\ln B)$ is the magnetic drift velocity, ${\bf B}=\hat{\bf b}B$ is the equilibrium magnetic field, $\Omega_{c}=qB/m_{i}$ is the ion cyclotron frequency; $q=Ze$ is the charge of an ion; $T_{i}$ and $T_{e}$ are the background ion and electron temperatures; $F_{0}=n_{0}(\pi v_{\scriptsize{\mathrm{T}}}^{2})^{-3/2}\exp(-v^{2}/v_{\scriptsize{\mathrm{T}}}^{2})$ is the background Maxwellian distribution function, where the ion thermal velocity is $v_{\scriptsize{\mathrm{T}}}=\sqrt{2T_{i}/m_{i}}$. The gyro-average is defined

[TABLE]

where \mbox{\boldmath\rho}(\vartheta)=\hat{\bf b}\times{\bf v}/\Omega_{c} and ${\bf v}=v_{\parallel}\hat{\bf b}+v_{\perp}[\hat{\bf e}_{1}\cos(\vartheta)+\hat{\bf e}_{2}\sin(\vartheta)]$, and $\hat{\bf e}_{1}\times\hat{\bf e}_{2}=\hat{\bf b}$. The electrostatic potential $\phi$ is determined by the quasi-neutrality constraint, assuming adiabatic electrons:

[TABLE]

Here, the non-zonal electrostatic potential is denoted $\tilde{\phi}=\phi-\overline{\phi}$ where the zonal average is defined in Eqn. 7 below, and the angle average at particle position ${\bf r}$ is given by

[TABLE]

Note that for simplicity trapped electrons have been neglected; the modified Boltzmann electron model is formally only valid in this limit.

2.1 Representing turbulence in general magnetic geometry

In order to keep the theory as general as possible, we would like to allow a representation that is slightly more general than the usual flux Fourier series applied in flux tube geometry. The eikonal representation is used [11], i.e. the fluctuations are decomposed as a sum of modes, as follows:

[TABLE]

The coefficients $\hat{\phi}_{S}=\hat{\phi}_{S}(\psi,\alpha,z)$, vary slowly in space, where $z$ denotes the coordinate of distance along the field line, $2\pi\psi$ is the toroidal magnetic flux and $\alpha$ is the final flux coordinate such that the equilibrium magnetic field is written ${\bf B}=\boldsymbol{\nabla}\psi\times\boldsymbol{\nabla}\alpha$. The eikonal functions $S=S(\alpha,\psi)$ also vary slowly in space, but the local perpendicular wavenumber ${\bf k}_{\perp}=\boldsymbol{\nabla}S=k_{\psi}\boldsymbol{\nabla}\psi+k_{\alpha}\boldsymbol{\nabla}\alpha$ is large, i.e. ${\bf k}_{\perp}\sim{\cal O}(1/\rho)$, so that the factor $\exp(iS)$ is a rapidly varying phase factor. The sum in Eqn. 5 can equivalently be written as a sum over the functions $k_{\psi}=\partial S/\partial\psi$ and $k_{\alpha}=\partial S/\partial\alpha$, i.e..

[TABLE]

Note that scale separation, i.e. $k_{\perp}\sim 1/\rho_{i}\gg L$, where $L$ is the variation scale of the equilibrium, ensures these forms are equivalent since $\boldsymbol{\nabla}_{\perp}\phi\approx\sum_{\bf k}i{\bf k}_{\perp}\hat{\phi}_{\bf k}$. It also implies that the wavenumbers $k_{\alpha}$ and $k_{\psi}$ are real, because otherwise $\phi$ would grow exponentially over a macroscopic distance and quickly violate the assumption of small fluctuations that underlies gyrokinetics, $q\phi/T_{i}\ll 1$. Thus, Eqn. 6 constitutes a local Fourier series, but the point we emphasize here is that the waveneumbers $k_{\alpha}$ and $k_{\psi}$ are free to also vary slowly in the $\alpha$ and $\psi$ coordinates. This fact can matter when considering turbulence in larger domains such as a flux surface or a full toroidal volume.

The flux surface average is (see e.g. [12])

[TABLE]

where the integral in $z$ is taken over one toroidal transit, and $V^{\prime}(\psi)=\int d\alpha\int dz/B(\psi,\alpha,z)$. Note that we have used the fact that, due to scale separation, the non-zero $k_{\alpha}$ terms will average to zero upon integration over $\alpha$. The integration over $\alpha$ still appears due to the fact that $\hat{\phi}$, $B$, etc. vary (smoothly) in $\alpha$.

3 Derivation of the global secondary instability in general geometry

The secondary instability of Rogers and Dorland [6] is a nonlinear mechanism for the generation of zonal flows by ion-scale turbulence. Two key assumptions make it possible to obtain an analytic result for this mode: (1) the Larmor scale is small $k_{\perp}^{2}\rho_{i}^{2}\ll 1$, and (2) the ion motion parallel to the magnetic field is slow, \mbox{k_{\parallel}}v_{\scriptsize{\mathrm{T}}}\ll\omega. The second assumption is called called “quasi-two-dimensional” (or “pseudo-three-dimensional” [13]) because parallel electron motion is retained, but the ion dynamics is only two-dimensional [14, 15]. That is, the two-dimensional motion of the ions is calculated in a three-dimensional domain, i.e. a magnetic geometry with shear, etc, with coupling along the field line being provided by the electrons, zeroing out the flux-surface-average density fluctuations. For the secondary instability that drives zonal flows, it is thus important to include information about the 3D mode structure of the primary mode. This is justified in the strongly ballooning limit, in which case parallel ion motion is slow, but nevertheless determines the the structure of the mode along the field line, i.e. its “ballooning structure” (see e.g. [16]). In this limit, the ballooning structure does not strongly affect the primary mode growth, but we will see that it does strongly affect the secondary instability. This is essentially because the global secondary instability has a strong zonal component, and it is therefore a non-local instability that depends on the properties of the turbulence across the entire flux surface. Intuitively, a zonal mode is driven more effectively by fluctuations that occupy a large fraction of a flux surface. This idea has been mentioned by previous authors [17, 18], but has not been investigated in detail.

We will now derive an expression for the global secondary instability. We first neglect the linear terms from the gyrokinetic equation. Then, assuming $k_{\perp}^{2}\rho_{i}^{2}\ll 1$, we may expand the Bessel functions and derive a fluid system (this is a 3D version of the zeroth-order system of [19]). Note that the fluid system is derived as an exact asymptotic limit of the gyrokinetic system, i.e. there is no truncation performed in the fluid moment hierarchy. The system is

[TABLE]

where we recall that $\tilde{\phi}=\phi-\overline{\phi}$, and note that the thermal Larmor radius varies in space, i.e.

[TABLE]

and $\chi$ (proportional to the ion gyrocenter perpendicular temperature) is defined

[TABLE]

For arbitrary functions $f$ and $g$, the Poisson bracket is defined

[TABLE]

We consider the primary mode as a given solution of the linearized gyrokinetic equation, and denote this solution by $\phi_{p}$, $\chi_{p}$. We emphasize that the $z$-dependence of the primary mode is determined by parallel ion motion, which therefore must be retained in the linearized gyrokinetic equation, but can be neglected in the system above, which is derived to describe the secondary instability. The primary mode has the form

[TABLE]

were c.c. means “complex conjugate,” $\chi_{p0}(\psi,\alpha,z)$ and $\phi_{p0}(\psi,\alpha,z)$ are complex functions. The zonal average of the primary is assumed to be zero $\partial S_{p}/\partial\alpha\neq 0$, i.e. it has no component with $k_{\alpha}=0$. In axisymmetric geometry, this assumption is trivial because $k_{\alpha}$ is everywhere the same, i.e. all field lines are the same. In non-axisymmetric geometry, $k_{\alpha}$ of a global mode is free to vary slowly over a surface, but we will assume it is nowhere zero. This is justified because a mode is locally stable wherever $k_{\alpha}=0$, contradicting the assumption that it is globally unstable.

We will now assume that the linear mode is sufficiently localized (in $\alpha$ and $z$) such that it does not significantly interact with itself nonlinearly (as modes that are very extended in ballooning space can do). For a global eigenmode (a mode that is everywhere on the surface evolving with the same complex frequency), this means that the mode structure in physical space is close to a single eikonal, i.e. the structure is the same in ballooning (the “covering space” [20]) and physical space. Alternatively it may be convenient to not require the primary mode to be a global mode of the system, but rather a linearly growing structure that is locally wave-like (i.e. a wave packet). In either case, this localization will enable us to focus on the secondary mode as a saturation mechanism, and rule out nonlinear self-interaction.

The secondary mode is denoted by $\phi_{s}$ and $\chi_{s}$, where $\phi_{s}\ll\phi_{p}$ and $\chi_{s}\ll\chi_{p}$. Now taking $\phi=\phi_{p}+\phi_{s}$ and $\chi=\chi_{p}+\chi_{s}$ and expanding the nonlinearity, Eqns. 10-12 yield

[TABLE]

The above differential system is linear in $\phi_{s}$ and $\chi_{s}$ with coefficients that are periodic (in their fast spatial variation) due to the dependence $\exp(iS_{p})$, and so Floquet theory applies. Despite the smooth three-dimensional variation that is present, the fast variation of the coefficients is only in one dimension, namely in the direction of $\boldsymbol{\nabla}S_{p}$. The direction that is binormal to $\boldsymbol{\nabla}S_{p}$ and $\hat{\bf b}$ can be considered ignorable, so we may assume a simple wave-like dependence, i.e. $\exp(i\lambda)$, such that $\hat{\bf b}\cdot\boldsymbol{\nabla}\lambda=\boldsymbol{\nabla}S_{p}\cdot\boldsymbol{\nabla}\lambda=0$. We have another free eikonal function corresponding to the Floquet exponent, which we will denote $\nu(\psi,\alpha)$; its gradient is parallel to that of $S_{p}$. In order to ensure that the secondary mode has a non-zero zonal component, we would like the overall factor $\exp(i\lambda+i\nu)$ to have a non-zero zonal average. This is ensured by requiring $\partial\nu/\partial\alpha=-\partial\lambda/\partial\alpha$. Thus we combine these eikonal functions into a single function $\bar{S}(\psi)=\lambda+\nu$.

We may now formally write the secondary mode in terms of a series in harmonics of $\exp(iS_{p})$:

[TABLE]

where we note that $\bar{S}=\bar{S}(\psi)$, and the $c_{n}$ and $d_{n}$ are functions of $z$, $\psi$ and $\alpha$. The secondary mode $\phi_{s}$ of Eqn. 21 can be split into zonal and non-zonal components

[TABLE]

whereas the primary mode has no zonal component, $\phi_{p}=\tilde{\phi}_{p}$. Now, substituting Eqns. 16, 23, and 24 into Eqn. 18, we obtain an equation for each component of $\tilde{\phi}_{s}$

[TABLE]

where $j=\phi_{p0}/|\phi_{p0}|$, and

[TABLE]

Thus, only three components of $\phi_{s}$, namely the principal component $\overline{c}_{0}$ and the two sidebands $c_{\pm 1}$, are needed for the calculation. For the temperature component of the secondary, Eqn. 22, the series does not truncate by itself, and so we must retain the infinite hierarchy of coefficients $d_{n}$. Fortunately, these may be computed using continued fractions; for brevity, we put this part of the calculation in A.

Finally we turn to Eqn. 19, which will yield the dispersion relation. Note that the first term may be simplified because $\tilde{c}_{0}=0$:

[TABLE]

where we define

[TABLE]

The remaining terms of Eqn. 19 consist of coupling terms between the sidebands, $c_{\pm 1}$ and $d_{\pm 1}$, and the primary mode. Manipulating these terms, Eqn. 19 yields the following

[TABLE]

where $j=\phi_{p0}/|\phi_{p0}|$, $r=\chi_{p0}/|\phi_{p0}|$, and

[TABLE]

We now substitute the expressions for $d_{\pm 1}$ (Eqns. 47-48) and $c_{\pm 1}$ (Eqns. 26-27) into Eqn. 32. All of the resulting terms are then proportional to $\overline{c}_{0}$, and this constant can be divided out, leaving the dispersion relation for the mode. After some algebra, we find the simple result

[TABLE]

where $\Omega_{\phi}$ is defined by Eqn. 29. Note that the radial wavenumber of the primary $\partial S_{p}/\partial\psi$ does not enter explicitly into the final result. It is also interesting to consider that when we can neglect the spatial variation of the magnetic field (i.e. $|\boldsymbol{\nabla}\psi|$, $B$, $\partial S_{p}/\partial\alpha$, etc.) the only affect of the geometry on the secondary mode is that the growth rate is determined by the RMS average of the primary mode on the flux surface. Thus, to induce similar secondary (zonal flow) growth rates, a localized mode must grow to a higher peak amplitude than a mode than is more uniformly distributed on a flux surface.

Remarkably, one can obtain the same dispersion relation as Eqn. 34 by performing the truncation $d_{n}=0$ for $|n|\geq 2$, which is commonly done to obtain an analytic result. This explains why such truncation has been found to be an excellent approximation to the full-series numerical gyrokinetic solutions (see e.g. [21]). Note that the unsheared slab result is obtained by taking the functions $\Omega_{\phi}$, $b_{\psi}$, $\phi_{p0}$, and $\chi_{p0}$ to be independent of $\alpha$ and $z$. The resulting expression can be compared with Eqn. (D43) of Ref. [19] (revealing a sign error in the earlier result).

3.1 Numerical simulation of zonal flow growth by secondary instability

Because Eqn. 34 is an exact asymptotic result, it should agree with fully gyrokinetic simulations, in arbitrary geometry, in the applicable limit, $k_{\perp}^{2}\rho^{2}\ll 1$. We test this using the GENE code [22], using a concentric circular model tokamak geometry (“$\hat{s}$-$\alpha$” geometry, without shift $\alpha=0$). This is a significant simplification, as compared with general stellarator geometry, but is a good initial test because it allows for nontrivial variation in the primary mode eigenfunctions along the field line. The parameter $\hat{s}$, representing the global shear in the magnetic field, is varied, along with the background temperature gradient, and the primary mode wavenumber. We can thereby sample parameter space to confirm that theory is able to uniformly describe different physical regimes.

For the model geometry, $|\boldsymbol{\nabla}\psi|$ is constant (independent of $z$), which somewhat simplifies the evaluation of $b_{\psi}$ in Eqn. 34. The expression can then be straightforwardly translated into simulation units and evaluated directly.

We perform direct numerical simulations of both the primary and secondary instabilities. The primary mode solution is used as input for the theoretical expression of Eqn. 34. The theoretical growth rate obtained can then be compared directly with the numerically observed growth rate, obtained by directly simulating the secondary mode, and fitting an exponential to the time dependence the evolving zonal potential $\overline{\phi}(t)$.

The secondary mode is simulated as follows. The primary mode is artificially frozen ($\partial/\partial t\rightarrow 0$), and the linear terms for the secondary mode are also removed, so its dynamics are fully determined by the nonlinear term. As long as parallel ion motion remains negligible (as assumed here) this is equivalent to a simulation without these modifications in the limit of large primary mode amplitude; note that the finite-$k_{\parallel}$ secondary mode (see [23]) has no zonal component, and is outside the scope of our work. Critically, our procedure avoids the need for tedious tuning of the simulation parameters to obtain a converged result. Another advantage of removing the linear terms, both theoretically and numerically, is that the solution obtained is proportional to the amplitude of the primary mode, so can be normalized to yield an easily reproducible result with a clear interpretation. Because very few wavenumbers are required for the mode, the (nonlinear) simulation of the secondary mode is not much more computationally expensive than the (linear) simulation of the primary mode.

We considered several values of $\hat{s}$ in the range of [math] to $1$, several temperature gradients $R/L_{T}$ in the range of $10$-$20$, and several wavenumbers with $k_{x}\rho$ and $k_{y}\rho$ in the range of $0.05$-$0.3$. The variation of the primary mode along the field line is sensitive to magnetic shear, allowing us to test the new prediction of theory, namely that the growth rate is sensitive to this variation. The set of parameters was also chosen to verify that the error, as measured by the difference between the theoretical and numerical growth rates, depends mainly on the size of the small parameter $k_{\perp}^{2}\rho^{2}$. This comparison is complicated by the fact that the small parameter of the theory in fact varies along the field line in the presence of finite magnetic shear. We thus chose to measure the largest value of $k_{\perp}$ as a representative value. As demonstrated by Fig. 1, we find agreement between the numerical and theoretical results, with corrections that scale linearly with the small parameter, as expected.

4 Conclusion

We have presented a theory of nonlinear zonal flow growth in general magnetic geometry, and validated it using nonlinear gyrokinetic simulations. The theory is valid in the electrostatic limit, and assumes adiabatic electrons. It applies to both axisymmetric (tokamak) and non-axisymmetric (stellarator) geometries. We note that even in the tokamak context, the role of magnetic geometry has never previously been accounted for, and so the result may be applicable to well-studied tokamak regimes. However, it will be especially interesting to apply the theory to the complex three dimensional geometries of stellarators, or tokamak-related problems with “extreme” magnetic geometry, e.g. in the edge where local magnetic shear is strong.

In the solution, given by Eqn. 34, magnetic geometry appears in two places. First, it enters through the structure of the primary mode, i.e. the variation of the electrostatic potential and temperature perturbation across the magnetic flux surface. Second, it enters through the average of the radial wavenumber, which depends on the relative distribution of flux surface compression on the surface. The first of these effects could be estimated by analytical models of linear instability, but can also be computed numerically using a gyrokinetic code, fully accounting for magnetic geometry, in local or global settings. We note also, that the choice of what model to use for the primary mode ultimately depends on what linear physics is believed to govern the nonlinear state, and therefore a “more complete” primary mode calculation may not actually be a better choice to try to understand the turbulence.

Finally, we would like to reiterate the point that the nonlinear growth mechanism is one process among several that determines the role of zonal flows in turbulence. That is, one must generally also take into account linear decay mechanisms (both collisionless and collisional) and nonlinear decay mechanisms. However, combining an understanding of those decay mechanisms with an estimate for growth, such as the one presented here, we believe it may be feasible, with the help of only linear gyrokinetic simulations, to model the full nonlinear response of zonal flows in fusion experiments with complex magnetic geometries.

5 Acknowledgements

The authors wish to thank T. Bird for discussions that initiated this work, and T. Görler and M. J. Pueschel, for help with the numerical calculation of secondary modes. This research has received funding from the European Research Council under the European Unions Seventh Framework Programme (FP7/2007V2013)/ERC Grant Agreement No. 277870. The gyrokinetic simulations presented in this work used resources of the National Energy Research Scientific Computing Center, a DOE Office of Science User Facility supported by the Office of Science of the U.S. Department of Energy under Contract No. DEAC02-05CH11231.

Appendix A Calculation of $\chi_{s}$

Eqn. 20 leads to expressions for the $d_{n}$. Using the fact that $c_{n}=0$ for $|n|\geq 2$ we write these as

[TABLE]

where we define $j=\phi_{p0}/|\phi_{p0}|$, $j_{+}=j$, $j_{-}=j^{*}$, $r=\chi_{p0}/|\phi_{p0}|$, $r_{+}=r$, and $r_{-}=r^{*}$, and

[TABLE]

We solve Eqn. 38 for positive and negative $n$ separately. For positive $n\geq 3$ we define $q_{n}=j^{*}d_{n}/d_{n-1}$ and find

[TABLE]

Iterating this we obtain

[TABLE]

This is a special continued fraction $x$ of the form $x=1/(y+x)$, which has the solution $x=(-y\pm\sqrt{y^{2}+4})/2$, i.e.

[TABLE]

Now for $n\leq-3$ we define $s_{n}=jd_{n}/d_{n+1}$ and similarly find

[TABLE]

which implies

[TABLE]

The signs in the solutions 42 and 44 for $q_{n}$ and $s_{n}$ can be determined by the condition $|q_{n}|<1$ and $|s_{n}|<1$, which is required for the convergence of the series in Eqn. 22. This leads to the result

[TABLE]

where $s$ is defined as

[TABLE]

with the sign chosen such that $|s|<1$. We may now combine Eqn. 45 with Eqns. 35-37 to obtain the following expressions for $d_{\pm 1}$:

[TABLE]

It is worth noting that the above expressions for $d_{\pm 1}$ are exactly the same as those obtained by simply truncating $d_{n}=0$ for $|n|\geq 2$.