The effect of background flow shear on gyrokinetic turbulence in the cold ion limit

TL;DR

This paper derives a fluid model from gyrokinetics in the cold ion limit, analyzing how background flow shear influences turbulence stability and nonlinear mode interactions, with results validated against gyrokinetic simulations.

Contribution

It introduces a fluid system including flow shear derived from gyrokinetics without closure, providing analytic solutions and benchmarks for turbulence studies.

Findings

Flow shear affects PVG mode stability.

Differences in nonlinear coupling between zonal and non-zonal modes.

Analytic results match gyrokinetic simulations.

Abstract

The cold ion limit of the local gyrokinetic model is rigorously taken to produce a nonlinear system of fluid equations that includes background flow shear. No fluid closure is required. By considering a simple slab geometry with magnetic drifts, but no magnetic shear, these fluid equations reduce to the Charney-Hasegawa-Mima model in the presence of flow shear. Analytic solutions to this model are found to study the impact of ExB flow shear on the stability of a single Parallel Velocity Gradient (PVG) driven mode. Additionally, the model is used to investigate the effect of background ExB flow shear on the basic three-mode nonlinear coupling, which reveals differences between zonal and non-zonal modes. These analytic results agree with gyrokinetic simulations and can serve to benchmark the numerical implementation of flow shear and nonlinear coupling.

Figures

Equations

$$\left(\frac{\partial }{\partial t}-{\omega }_{V\perp }x\frac{\partial }{\partial y}\right)$$

$$+{\vec{v}}_{Mi}\cdot \vec{\nabla }{\overline{h}}_i=\frac{1}{B}\frac{\partial ⟨\overline{\varphi }⟩}{\partial y}{F}_{Mi}\left[\frac{1}{l_n}+\left(\frac{m_i{v}^2}{2T_i}-\frac{3}{2}\right)\frac{1}{l_T}+\frac{m_i{v}_{||}}{T_i}{\omega }_{V||}\right],$$

$${\vec{v}}_{Mi}\equiv \frac{v_{||}^2}{{\Omega }_i}\hat{b}\times \vec{\kappa }+\frac{\mu }{{\Omega }_i}\hat{b}\times \vec{\nabla }B$$

$$a_{||i}\equiv -\mu \hat{b}\cdot \vec{\nabla }B,$$

$${\omega }_{V\perp }\equiv -\frac{d}{dx}\left(\frac{{\vec{E}}_0\times \hat{b}}{B}\cdot \vec{\nabla }y\right)=\frac{\partial }{\partial x}\left(\frac{1}{JB}\frac{d{\Phi }_0}{dx}\right),$$

$$Z_{i}e{\int |}_x{d}^{3}v{\overline{h}}_i=\left(\frac{Z_i^2{e}^2{n}_i}{T_i}+\frac{e^2{n}_e}{T_e}\right)\overline{\varphi },$$

$$(\frac{\partial }{\partial t}$$

$$+\frac{1}{B}\sum _{{\vec{k}}^{\prime }}\left({\vec{k}}^{\prime }\times {\vec{k}}^{\prime \prime }\right)\cdot \hat{b}\left(g_i^{\prime }+\frac{Z_{i}e{F}_{Mi}}{T_i}{\varphi }^{\prime }{J}_0\left(k_{\perp }^{\prime }{\rho }_i\right)\right){\varphi }^{\prime \prime }{J}_0\left(k_{\perp }^{\prime \prime }{\rho }_i\right)$$

$$=-\frac{Z_{i}e{F}_{Mi}}{T_i}\left[v_{||}\hat{b}\cdot \vec{\nabla }\left(J_0\left(k_{\perp }{\rho }_i\right)\varphi \right)+i{\vec{k}}_{\perp }\cdot {\vec{v}}_{Mi}{J}_0\left(k_{\perp }{\rho }_i\right)\varphi -a_{||i}\frac{m_i{v}_{||}}{T_i}{J}_0\left(k_{\perp }{\rho }_i\right)\varphi \right]$$

$$+i\frac{k_y}{B}{J}_0\left(k_{\perp }{\rho }_i\right)\varphi F_{Mi}\left[\frac{1}{l_n}+\left(\frac{m_i{v}^2}{2T_i}-\frac{3}{2}\right)\frac{1}{l_T}+\frac{m_i{v}_{||}}{T_i}{\omega }_{V||}\right],$$

$$\int d^{3}v{g}_i{J}_0\left(k_{\perp }{\rho }_i\right)=\frac{e{n}_i}{T_e}\left\{Z_i\frac{T_e}{T_i}\left[1-I_0\left(k_{\perp }^2{\rho }_{thi}^2\right)\exp \left(-k_{\perp }^2{\rho }_{thi}^2\right)\right]+1\right\}\varphi ,$$

$$K_x$$

$$Z$$

$${\frac{\partial g_i}{\partial t}|}_{K_x}$$

$$+$$

$$=$$

$$+$$

$${\frac{\partial g_i}{\partial t}|}_{K_x}$$

$$+$$

$$=$$

$$+$$

$$\delta n\equiv \int d^{3}v{g}_i=\frac{e{n}_i}{T_e}\left(1+k_{\perp }^2{\rho }_S^2\right)\varphi .$$

$$\left({\frac{\partial }{\partial t}|}_{K_x}+u_f\frac{\partial }{\partial Z}\right)\delta n$$

$$+i\frac{e{n}_i}{T_e}\left(k_x{\rho }_S{\omega }_{Mx}+k_y{\rho }_S{\omega }_{My}\right)\varphi =0,$$

$${\omega }_{Mx}$$

$${\omega }_{My}$$

$$({\frac{\partial }{\partial t}|}_{K_x}$$

$$+\frac{1}{2B}\sum _{{\vec{K}}^{\prime }}\left({\vec{K}}^{\prime }\times {\vec{K}}^{\prime \prime }\right)\cdot \hat{b}\left(k_{\perp }^{\prime 2}-k_{\perp }^{\prime \prime 2}\right){\rho }_S^2{\varphi }^{\prime }{\varphi }^{\prime \prime }+i\left(k_x{\rho }_S{\omega }_{Mx}+k_y{\rho }_S{\omega }_{My}\right)\varphi =0.$$

$$\left({\frac{\partial }{\partial t}|}_{K_x}+u_f\frac{\partial }{\partial Z}\right)\delta u_{||}$$

$$+\frac{Z_{i}e}{m_i}\hat{b}\cdot \vec{\nabla }Z\left(\frac{\partial \varphi }{\partial Z}+\frac{\partial \ln B}{\partial Z}\varphi \right)-i{k}_y\frac{{\omega }_{V||}}{B}\varphi =0.$$

$$({\frac{\partial }{\partial t}|}_Y$$

$$+\frac{1}{B}\left({\rho }_S^2\vec{\nabla }\left({\nabla }_{\perp }^2\varphi \right)\times \vec{\nabla }\varphi \right)\cdot \hat{b}+\left({\omega }_{Mx}{\rho }_S\frac{\partial \varphi }{\partial x}+{\omega }_{My}{\rho }_S\frac{\partial \varphi }{\partial Y}\right)=0$$

$$\left({\frac{\partial }{\partial t}|}_Y+u_f\frac{\partial }{\partial Z}\right)\delta u_{||}$$

$$+\frac{Z_{i}e}{m_i}\hat{b}\cdot \vec{\nabla }Z\left(\frac{\partial \varphi }{\partial Z}+\frac{\partial \ln B}{\partial Z}\varphi \right)-\frac{{\omega }_{V||}}{B}\frac{\partial \varphi }{\partial Y}=0,$$

$${\frac{\partial }{\partial t}|}_{K_x}$$

$$+i\left(k_x{\rho }_S{\omega }_{Mx}+k_y{\rho }_S{\omega }_{My}\right)\varphi =0$$

$${\frac{\partial }{\partial t}|}_{K_x}\delta u_{||}$$

$$\varphi =\hat{\varphi }{\left(1+k_{\perp }^2{\rho }_S^2\right)}^{-1}\exp \left(-i{\int |}_{K_x}dt\frac{k_x{\rho }_S{\omega }_{Mx}+k_y{\rho }_S{\omega }_{My}}{1+k_{\perp }^2{\rho }_S^2}\right)$$

$$\left(1+k_{\perp }^2{\rho }_S^2\right){\frac{{\partial }^2\hat{\varphi }}{\partial t^2}|}_{K_x}$$

Full text

1 Introduction

The presence of sheared flow can significantly alter the turbulence in a magnetized plasma. Sheared flow is thought to be important for tokamak experiments [1, 2, 3] as well as many astrophysical phenomena [4, 5, 6, 7, 8]. However, its effect is not straightforward to understand. A gradient in the flow parallel to the magnetic field is able to drive turbulence through what is called the Parallel Velocity Gradient (PVG) instability [9, 10]. On the other hand, shear in the flow perpendicular to the magnetic field line has been shown to stabilize turbulence and reduce its ability to transport heat and particles [11, 12].

Perpendicular flow shear in tokamaks is particularly important as it may enable the improved confinement regime H-mode [3]. H-mode is a robust plasma phenomenon that improves the energy confinement time of the device by roughly a factor of two and is currently viewed as essential for the success of the ITER experiment [13] as well as a future power plant [14, 15]. Thus, a full understanding of perpendicular flow shear is desirable to ensure that a prospective design can properly exploit H-mode.

To this end, nonlinear gyrokinetic computer codes have recently been employed to study this problem with unprecedented realism and accuracy [16, 17, 18, 19]. However, such codes, which can be hundreds of thousands of lines of code [20, 21, 22], are challenging to write, run, and interpret. Thus, analytic solutions that hold in clearly defined parameter regimes can be very valuable. They can illuminate aspects of the underlying physics and serve as benchmarks to verify the code and its execution is correct.

Such analytic efforts, of course, predate gyrokinetic codes and the supercomputers they require. The PVG instability was originally discovered using simple fluid equations in a slab geometry without magnetic shear [9]. Subsequent work rigorously derived PVG from kinetic theory [23] and included the effect of magnetic shear [10]. A more recent analytic calculation [24] took the cold ion limit and calculated how PVG can be stabilized by perpendicular shear flow. References [25, 26, 27] have used fluid equations, rigorously derived from electromagnetic gyrokinetic theory, to understand the interplay between the density gradient, temperature gradient, parallel flow shear, perpendicular flow shear, and magnetic shear. Lastly, several works [28, 29] study how perpendicular flow shear enables gyrokinetic simulations to exhibit “subcritical” turbulence, which is a nonlinear instability that is sustained despite the fact that the modes are linearly stable.

As evidenced by all of these separate works, there are several different physical effects at play that are challenging to include simultaneously and rigorously. First, the effect of perpendicular flow shear is to take the coherent structures in the plasma turbulence and, with time, gradually shear them to higher wavenumbers. This introduces an explicit time dependence that makes the evolution of any turbulence mode more complex. Previous efforts detailed above (with the exception of reference [28]) have dealt with this time dependence by assuming that the solution is an exponential with a slowly-varying growth rate. This is useful for developing an intuition for the dynamics of the system, but is fundamentally an approximation that is only valid when the mode growth rate is much larger than the perpendicular flow shearing rate. This approximation breaks down when the perpendicular flow is comparable to the growth rate, which is necessarily the case to stabilize turbulence [11]. Second, as the modes are sheared to high wavenumbers, finite gyroradius effects become increasingly important. These are generally ignored (excepting references [24, 28]), but are important as they can determine how the mode is ultimately stabilized.

Third, for fusion devices we are fundamentally interested in properties of the nonlinearly saturated state of turbulence, rather than characteristics about its linear growth. Studying how modes nonlinearly interact can illuminate the physics of this. Without perpendicular flow shear, the usual three-wave nonlinear coupling and three-wave instability analysis only involves time-independent equations, which lead to the usual instabilities that grow exponentially in time. But, as for the linear problem described in the previous paragraph, introducing perpendicular flow shear causes the equations to gain an explicit time-dependence. Thus, the nonlinear problem gains complicated time dependence that can exhibit transient growth [25, 27]. Interestingly, we will see that there is a surprising correspondence between the linear and nonlinear analysis. Only recently have nonlinear dynamics been considered analytically in the presence of sheared flows. Reference [30] uses a fluid model (without finite gyroradius effects) to investigate the stability of zonal flows. Additionally, reference [31] uses the cold ion limit to study the nonlinear dynamics streamers driven by the PVG instability.

In this work, we take the cold ion limit to explore the impact of parallel and perpendicular flow shear. This limit, though not physically motivated, is attractive because it enables simple and exact results. It was used most prominently in deriving the Charney-Hasegawa-Mima equation from fluid models [32] and was quickly employed to find nonlinear analytic results [33]. Since then, the cold ion limit of gyrokinetics has been taken [34, 35] in slab geometry and has potential applications ranging from the scrape-off layer of tokamak plasmas [36] to the solar wind [37].

In section 2, we start from the full nonlinear local gyrokinetic equation in general geometry and use the cold ion limit to rigorously derive a simple fluid model that exactly governs its behavior. No fluid closure is needed because the cold ion limit naturally provides it. From this point onwards, we consider a simple slab-like geometry without magnetic shear and show that our fluid equations reduce to the Charney-Hasegawa-Mima model [32] in the presence of background flow. Ignoring magnetic shear reduces the realism of the model in space, but enables investigations of the full time evolution of a single physical Fourier mode driven by PVG and stabilized by perpendicular flow. This is done in section 3 using a Fourier representation that is very similar to that used in gyrokinetic codes, making comparison relatively simple. Next, in section 4, we focus on the nonlinear coupling term and study how three-wave coupling is affected by perpendicular flow shear. The results of both sections 3 and 4 are compared against gyrokinetic simulations as they are presented. Section 5 provides some concluding remarks.

2 Derivation of the model

We begin with the electrostatic, collisionless, local (flux-tube) $\delta f$ gyrokinetic model in a general geometry [38, 39, 40]. The electrons are assumed to respond adiabatically to the motion of a single ion species. We will set the background flow in the center of our domain to be zero, but allow for the flow to vary in the $\vec{\nabla}x$ direction, which is perpendicular to the direction of the magnetic field $\hat{b}$. We refer to the gradient of the component of the flow parallel to the magnetic field line as “parallel flow shear,” while the gradient in the component perpendicular to the magnetic field is called “perpendicular flow shear.” Including both of these effects enables study of the PVG instability as well as stabilization by $\vec{E}\times\vec{B}$ flow shear. All other gradients of the plasma quantities are also assumed to vary in only in the $\vec{\nabla}x$ direction. The binormal spatial coordinate $y$ is defined such that $\vec{\nabla}y\times\vec{\nabla}x$ is in the $\hat{b}$ direction. We will denote the real-space coordinate along the field line as $z$. Note that this coordinate system is not orthogonal in general and permits the Jacobian to be different from one. However, we will only encounter $J\equiv|(\vec{\nabla}y\times\vec{\nabla}x)\cdot\hat{b}|^{-1}$, which is closely related to the coordinate system Jacobian. Thus, it is appropriate for both a slab representation [25] (i.e. a Cartesian $y$, $x$, and $z$ such that $J=1$) or a toroidal ballooning representation [39] (i.e. $y$ as the binormal angle $\alpha$, $x$ as the poloidal flux $\psi$, and $z$ as the poloidal angle $\theta$ such that $J=B^{-1}$). Note that we have defined the right-handed coordinate system as $(y,x,z)$ in order to be consistent with the ballooning representation.

In this context, the real-space ion gyrokinetic equation is given by

[TABLE]

where $t$ is the time coordinate and the velocity-space coordinates are the velocity parallel to the magnetic field $v_{||}$ and the magnetic moment $\mu\equiv v_{\perp}^{2}/(2B)$ (which is defined with the perpendicular velocity $v_{\perp}$). The magnetic drifts are given by

[TABLE]

and the linear parallel acceleration (i.e. the mirror effect) is

[TABLE]

where $\vec{\kappa}\equiv\hat{b}\cdot\vec{\nabla}\hat{b}$ is the magnetic field curvature vector. The background gradients are the shear in the perpendicular flow

[TABLE]

the shear in the parallel flow $\omega_{V||}$, the ion density gradient scale length $l_{n}^{-1}\equiv-J^{-1}d\ln\left(n_{i}\right)/dx$, and the ion temperature gradient scale length $l_{T}^{-1}\equiv-J^{-1}d\ln\left(T_{i}\right)/dx$. The shear in the parallel flow is $\omega_{V||}\equiv-J^{-1}(RB_{\zeta}/B)d(V/R)/dx$ in toroidal geometry and $\omega_{V||}\equiv-J^{-1}(\hat{b}\cdot\hat{\zeta})dV/dx$ in slab geometry, where $V$ is the flow velocity, $R$ is the major radius, and $\hat{\zeta}$ is the direction of the flow. Additionally, $B$ is the magnitude of the magnetic field, $\vec{E}_{0}=-\vec{\nabla}\Phi_{0}$ is the background electric field, $\Phi_{0}=\Phi_{0}(x)$ is the background electrostatic potential, $Z_{i}$ is the ion charge number, $e$ is the elementary charge, $F_{Mi}\equiv n_{i}\left(m_{i}/2\pi T_{i}\right)^{3/2}\exp\left(-m_{i}v^{2}/2T_{i}\right)$ is the background ion Maxwellian distribution function, $n_{i}$ is the background ion number density, $m_{i}$ is the ion mass, $\Omega_{i}\equiv Z_{i}eB/m_{i}$ is the ion gyrofrequency, $T_{i}$ is the background ion temperature, and $v$ is the particle speed coordinate. The unknowns in this equation are $\overline{h}_{i}\equiv\overline{\delta f}_{i}+\left(Z_{i}eF_{Mi}/T_{i}\right)\overline{\phi}$, the nonadiabatic portion of the fluctuating ion distribution function in real-space, and $\overline{\phi}$, the fluctuating electrostatic potential in real-space. Note that $\overline{\delta f}_{i}$ is the fluctuating ion distribution function in real-space and $\left\langle\ldots\right\rangle\equiv(2\pi)^{-1}\oint_{0}^{2\pi}d\varphi\left(\ldots\right)$ is the particle gyroaverage over gyro-angle (taken at constant particle guiding center). Since equation (1) has two unknowns, $\overline{h}_{i}$ and $\overline{\phi}$, we also require the real-space quasineutrality equation to close the system:

[TABLE]

where the species index indicates either $i$ for ions or $e$ for electrons and $\left.\int\right|_{x}d^{3}v$ indicates that the integral must be taken at constant particle position (not the guiding center position).

Since we are working in the local flux-tube limit of $\delta f$ gyrokinetics, the background gradients are fixed constants and it is appropriate to use periodic boundary conditions in the directions perpendicular to the magnetic field lines. Thus, it is convenient to perform a Fourier analysis of $\overline{h}_{i}$ and $\overline{\phi}$ in the radial and binormal directions (e.g. $\overline{\phi}=\sum_{k_{x}}\sum_{k_{y}}\phi\exp\left(ik_{x}x+ik_{y}y\right)$). This has the advantage of converting the averages over gyro-angle into Bessel functions [39, 41]. The only caveat is that, because the nonlinear term includes a product of $\overline{h}_{i}$ and $\overline{\phi}$, the Fourier transform includes a convolution (which involves three modes). In Fourier-space, equation (1) becomes [39, 41]

[TABLE]

where the summation is performed over the full $k^{\prime}_{x}\in\left(-\infty,\infty\right)$, $k^{\prime}_{y}\in\left(-\infty,\infty\right)$ plane, the nonlinear coupling condition of the convolution is $\vec{k}^{\prime\prime}=\vec{k}-\vec{k}^{\prime}$, and the prime and double prime symbols indicate which wavenumber is used in evaluating the quantity. Here $k_{x}$ and $k_{y}$ are the $x$ and $y$ Fourier wavenumbers, $g_{i}(\vec{k})\equiv h_{i}(\vec{k})-\left(Z_{i}eF_{Mi}/T_{i}\right)J_{0}\left(k\rho_{i}\right)\phi(\vec{k})$ is the complementary ion distribution function in Fourier-space, $J_{0}\left(\ldots\right)$ is the $0^{\text{th}}$ order Bessel function of the first kind, and $\rho_{i}\equiv\sqrt{2\mu B}/\Omega_{i}$ is the ion gyroradius.

In Fourier-space, the quasineutrality equation becomes

[TABLE]

where $\rho_{thi}\equiv\sqrt{T_{i}/m_{i}}/\Omega_{i}$ is the ion thermal gyroradius and the integral is now taken at constant guiding center position. Such a Fourier representation facilitates comparison to (and understanding of) gyrokinetic codes, most of which use such a representation.

Next, we rename the dummy variables $\vec{k}^{\prime}\rightarrow\vec{k}^{\prime\prime}$ and $\vec{k}^{\prime\prime}\rightarrow\vec{k}^{\prime}$ in equation (6), then sum the result with the original equation (6) to make the $(\vec{k}^{\prime}$, $\vec{k}^{\prime\prime})$ symmetry of the nonlinear term explicit. By changing to the double-shearing coordinate system [42, 43], we can elegantly treat both the effects of flow shear and magnetic shear. This coordinate system accounts for the fact that magnetic shear causes an eddy to be radially sheared as you move along a field line, while perpendicular flow shear causes an eddy to be radially sheared as time progresses. A more detailed explanation is given in reference [25]. In Fourier-space, this coordinate transform is given by

[TABLE]

where $u_{f}\equiv\omega_{V\perp}l_{s}$ is the velocity at which the unsheared eddy appears to move along the field line. The parameter $l_{s}$ is the global magnetic shear scale length, which is $l_{s}^{-1}\equiv d(\hat{b}\cdot\vec{\nabla}y)/dx$ in slab coordinates or $l_{s}^{-1}\equiv dq/dx$ in ballooning coordinates (where $q$ is the safety factor). In these new coordinates, the gyrokinetic equation becomes

[TABLE]

where the summation is still over the full plane and $\vec{K}^{\prime\prime}=\vec{K}-\vec{K}^{\prime}$. Note that, since the time derivative is now taken at constant $K_{x}$, the quantities of $k_{x}$ appearing in $\vec{k}_{\perp}$ and the Bessel functions have gained an explicit time dependence.

At this point, as in reference [33], we will take the cold ion limit $T_{i}\ll Z_{i}T_{e}$. Note that we maintain $k_{\perp}\rho_{S}\sim 1$, where $\rho_{S}\equiv\sqrt{Z_{i}T_{e}/m_{i}}/\Omega_{i}$ is the sound gyroradius. In this limit, $J_{0}\left(k_{\perp}\rho_{i}\right)\rightarrow 1$, while the quasineutrality condition retains the lowest order polarization drift. Thus, the gyrokinetic equation becomes

[TABLE]

and quasineutrality becomes

[TABLE]

It is important to note the finite sound gyroradius effect that survives in the quasineutrality equation. This is due to the ion polarization drift and will have important consequences later in this work. Taking velocity-space moments of equation (11), gives a fluid model. We find that no closure is needed because, in the cold ion limit, the right-hand side is zero for all except the density and parallel velocity moments. Thus, $\int d^{3}v~{}v_{||}^{a}\mu^{b}g_{i}=0$ for integers $a\geq 2$ and $b\geq 1$. The density moment is

[TABLE]

where

[TABLE]

contain the effect of the magnetic drifts and density gradient, $\kappa_{x}\equiv\partial\vec{r}/\partial x\cdot\vec{\kappa}$ and $\kappa_{y}\equiv\partial\vec{r}/\partial y\cdot\vec{\kappa}$ are the components of the magnetic field line curvature, and $c_{S}\equiv\sqrt{Z_{i}T_{e}/m_{i}}$ is the sound speed. Substituting quasineutrality (i.e. equation (12)) gives an equation for the evolution of $\phi$:

[TABLE]

Lastly, the parallel velocity moment, $\delta u_{||}\equiv n_{i}^{-1}\int d^{3}v~{}v_{||}g_{i}$, of equation (11) is

[TABLE]

Note that, in the cold ion limit, the effect of the ion temperature gradient vanishes entirely from the model. In real-space, equations (16) and (17) are

[TABLE]

and

[TABLE]

where $\vec{\nabla}_{\perp}\equiv\vec{\nabla}-\hat{b}(\partial/\partial z)$. Note that in real-space, the double shearing coordinate system [25] is given by equation (9) and $Y\equiv y-x(z+u_{f}t)/l_{s}$, which means that $\partial/\partial x=\left.\partial/\partial x\right|_{Y}-(Z/l_{s})\partial/\partial Y$.

Equations (16) and (17) form a closed system that governs the electrostatic, collisionless local gyrokinetic model in the cold ion limit. This fluid model is tremendously simpler than the six-dimensional, integro-differential system of gyrokinetics. Yet, as long as the limit of $T_{i}\ll Z_{i}T_{e}$ is satisfied, the two are equivalent.

This cold ion model is simple enough to enable analytic results for the linear dynamics in a slab with magnetic shear, as was done by Waelbroeck, et al. [24]. However, in this work we will sacrifice realism in space in order to enable a more realistic treatment in time. Instead of a slab geometry with magnetic shear, we will use a slab without magnetic shear (but maintain simple magnetic drifts). Accordingly, instead of having to solve differential equations in space, as was done by Waelbroeck, we will solve them in time. This will permit an investigation of the full time evolution of modes under the effect of finite $\omega_{V\perp}$ as well as their nonlinear interaction.

Without magnetic shear, we can let $l_{s}\rightarrow\infty$. Moreover, to prevent the coordinate system from diverging, we must also set $u_{f}=0$. This is physically motivated because, without magnetic shear, the unsheared eddy no longer has any apparent motion along the field line. Thus, we can ignore the last term in equations (8) and (9) as well as the terms containing $u_{f}$ in equations (16) and (17). Importantly, ignoring magnetic shear allows us to apply a standard periodic boundary condition in the parallel direction, rather than the more complex twist-and-shift boundary condition [44]. However, this simplification is also a limitation as it eliminates mode coupling through the parallel boundary condition. Such coupling occurs in toroidal devices, but will not be included in our model. Lastly, we will assume that $dB/dZ=0$, so that $\hat{b}$ becomes a direction of symmetry. We will maintain magnetic drifts, but also assume that they are constant in $Z$. This is applicable to geometries like purely toroidal field lines or straight field lines with a perpendicular gradient in the field strength. With these simplifications, it becomes useful for notational simplicity to adopt a Cartesian orthonormal coordinate system such that $|\vec{\nabla}x|=|\vec{\nabla}y|=|\vec{\nabla}Z|=1$, $\hat{b}\cdot\vec{\nabla}Z=1$, and all the metric coefficients are $1$. Additionally, we can Fourier analyze in the parallel direction and replace $d\phi/dZ$ with a parallel wavenumber $ik_{||}\phi$. Doing so is only possible because all equilibrium quantities have become independent of $Z$. In this geometry, our model becomes

[TABLE]

and

[TABLE]

which is equivalent to the Charney-Hasegawa-Mima model [32] with the addition of background flow. We note that in the slab limit of most toroidal gyrokinetic codes $\omega_{Mx}=0$ and $\omega_{My}=-c_{S}/l_{n}$.

3 Single mode slab results

By restricting our analysis to the evolution of just a single mode, the nonlinear terms vanish from the $\phi$ evolution equation and the parallel velocity moment (i.e. equations (20) and (21)). Qualitatively similar systems have been analyzed in the past. Some used computational approaches to investigate realistic geometries [11, 42], while others obtained analytic results for a simplified fluid model [29]. The physics discussed in these studies are useful in interpreting the results of this section, which will use the cold ion limit and above simplifications to enable rigorous analytic results.

Taking the time derivative of the equation (20) and substituting equation (21) as well as the ansatz

[TABLE]

produces

[TABLE]

This ordinary differential equation determines the full time evolution of a single mode in the presence of parallel and perpendicular velocity shear. While it appears simple, due to the perpendicular flow shear dependence of $k_{x}=K_{x}+k_{y}\omega_{V\perp}t$ (see equation (8)), $k_{\perp}^{2}$ is quadratic in time and $k_{x}$ is linear. Additionally, we note that the ansatz of equation (22) is important to keep in mind because, even with $\omega_{V||}=k_{||}=0$ the mode will still evolve due to the variation of the $1+k_{\perp}^{2}\rho_{S}^{2}$ factor as it is advected by perpendicular flow shear.

While equation (23) is complicated, it still has an analytic solution given by

[TABLE]

where $\prescript{}{2}{F}_{1}$ is the Gaussian hypergeometric function [45],

[TABLE]

and

[TABLE]

are the coefficients of the polynomials appearing in equation (23). Here $C_{j}$ for any Greek letter $j$ is an integration constant that can be calculated from the initial conditions. This analytic solution has been verified against the numerical solution of equation (23) as shown in figure 1. Additionally, it has been verified by comparison with the local version of the gyrokinetic code GENE [46, 47] as shown in figure 2. GENE is one of the most commonly used gyrokinetic codes with the capacity to solve the full nonlinear gyrokinetic system of equations in multiple different geometries. Among many capabilities, it can model collisions, electromagnetic fluctuations, and global effects, although these features are not used in this work.

All GENE simulations in this work use $Z_{i}=1$ and an ion temperature of $T_{i}=10^{-4}T_{e}$ to ensure that the cold ion limit is well satisfied. To be consistent with the analytic model, all simulations also set the magnetic shear equal to zero and use the standard slab geometry in GENE. Because this slab model does not include magnetic drifts, $\omega_{Mx}$ is forced to be zero in all simulations, but $\omega_{My}$ is varied using the density gradient. Additionally, for all comparisons it was necessary to set the flux surface averaged value of the electrostatic potential to zero when GENE calculates the adiabatic electron response. This is because, without magnetic shear, the standard parallel boundary condition used by most local gyrokinetic codes causes every field line to close on itself, so flux surfaces are not formed [44]. All simulations use the same resolution of 32 grid points in the parallel direction, 32 parallel velocity grid points, and 24 magnetic moment grid points. To ensure that the “wavevector-remap” scheme converges to a fairly smooth and continuous mode evolution, 128 radial wavenumbers were used [48]. Since we will only ever initialize a limited number of individual modes to be finite, the number of binormal wavenumbers does not have an effect, so very low values between 2 and 10 were used.

Unfortunately, hypergeometric functions do not provide much physical insight. Thus, as was done in reference [28], we will find further analytic results by investigating equation (23) in certain limits. We will focus on $|\omega_{V\perp}t|\ll 1$ and $|\omega_{V\perp}t|\gg 1$. The first limit can be used to understand the effect of weak perpendicular flow shear or the initial effect of arbitrary flow shear. The second limit can be used to understand the behavior of a mode after a long time as it is advected to large $k_{x}$ or the behavior of a mode with large $k_{x}$ as it is advected towards $k_{x}=0$. One could also find solutions in the weak perpendicular flow shear limit of $\omega_{V\perp}\ll\gamma$, where $\gamma$ is the mode growth rate. We will see that this is closely related to the $|\omega_{V\perp}t|\ll 1$ limit.

3.1 The $|\omega_{V\perp}t|\ll 1$ limit

By looking in the limit of short time, we can obtain analytic results for situations without perpendicular flow shear as well as the early time behavior of situations with flow shear. To lowest order (i.e. $\omega_{V\perp}=0$), we see that $k_{x}=K_{x}=\text{const}$ and $k_{\perp}=K_{x}^{2}+k_{y}^{2}=\text{const}$, so the time dependence of the wavenumber coefficients vanish. Thus, equation (23) becomes a simple 2nd order ordinary differential equation with constant coefficients. This is solved by an exponential, which can be substituted into equation (22) to find the full solution (without the hat on $\phi$) to be

[TABLE]

where the numerical subscript represents the quantity’s order in the $|\omega_{V\perp}t|$ expansion and

[TABLE]

are the mode frequency and growth rate respectively. Here the superscript $\varnothing$ indicates the initial condition of the quantity (e.g. $k_{x}^{\varnothing}\equiv k_{x}(t=0)=K_{x}$, $k_{\perp}^{\varnothing 2}=K_{x}^{2}+k_{y}^{2}$). When the mode growth rate is much greater than the perpendicular flow shear (i.e. $\omega_{V\perp}\ll\gamma$), this solution actually holds for all time. We can simply replace $K_{x}\rightarrow k_{x}$ and $k_{\perp}^{\varnothing}\rightarrow k_{\perp}$ to capture their time dependence because they vary slowly in the $\omega_{V\perp}\ll\gamma$ expansion. Note that, when the magnetic drifts and finite sound gyroradius effects are ignored, equation (37) reduces to the typical PVG growth rate [10, 25].

For instability, the growth rate $\gamma$ must be real, giving the following condition on the parallel velocity gradient:

[TABLE]

This shows that the only possible instability in this limit is PVG. As is intuitive, in the cold ion limit neither ion temperature gradient (ITG) nor drift wave instabilities can exist. Without parallel flow shear, equation (35) simply governs stable drift-sound waves as they oscillate. On the other hand, we see that a sufficiently large parallel velocity gradient will always overcome the damping and lead to instability, so long as $k_{||}$ and $k_{y}$ are finite. To maximize the growth rate, a mode should have a parallel wavenumber of $k_{||}c_{S}=k_{y}\rho_{S}\omega_{V||}/2$, and as large of a $k_{y}$ as possible. Given this value of $k_{||}$, the stability criterion becomes

[TABLE]

Note that, due to finite sound gyroradius effects in equation (37), the growth rate has a maximum possible value of $\gamma=|\omega_{V||}|/2$, rather than increasing without bound when $k_{y}\rightarrow\infty$ as is the usual result [10, 25].

The dependence of the instability growth rate on $K_{x}$ is complex as shown in figure 3. For example, when the magnetic drifts and density gradient are zero such that $\omega_{My}=0$, it is best for $K_{x}\rho_{S}=0$ to minimize the finite gyroradius damping. However, when $\omega_{My}\neq\omega_{Mx}=0$ the maximum growth rate occurs at $K_{x}\rho_{S}=\pm\sqrt{2\omega_{My}^{2}-(1+k_{y}^{2}\rho_{S}^{2})\omega_{V||}^{2}}/\omega_{V||}$. As can be seen from equation (37), when both $\omega_{Mx}$ and $\omega_{My}$ are non-zero the growth rate spectrum even becomes asymmetric about $K_{x}\rho_{S}=0$.

If the magnetic drifts, density gradient, and finite gyroradius effects are neglected, these results match equation (30) of reference [25] and equation (22) of reference [10] when $T_{i}=0$. This is true even though magnetic shear is included in the analysis of reference [25]. Additionally, the effect of finite gyroradius effects is in agreement with references [24, 31].

To next order in $|\omega_{V\perp}t|\ll 1$, the differential equation for the mode becomes

[TABLE]

Perpendicular flow shear appears only in the inhomogeneous terms through $\left(k_{\perp}^{2}\rho_{S}^{2}\right)_{1}=2K_{x}k_{y}\rho_{S}^{2}\omega_{V\perp}t$ and $\left(k_{x}\rho_{S}\right)_{1}=k_{y}\rho_{S}\omega_{V\perp}t$. This equation can be solved analytically. The solution to the homogeneous equation has the same form as the lowest order solution of equation (35), while a particular solution to the inhomogeneous equation is

[TABLE]

where we note the hat on $\hat{\phi}_{1}$. Thus, we see that the dominant effect of weak flow shear (or alternatively the first effect of flow shear to appear) is a quadratic correction to the lowest order exponential behavior. The lowest and next order analytic solutions are compared to the exact numerical solution of equation (23) in figure 4. Note that care must be taken in converting from $\hat{\phi}_{1}$ to $\phi_{1}$ because equation (22) must be expanded to next order in $|\omega_{V\perp}t|\ll 1$.

While equation (41) is fairly complicated, it can be further simplified to give an expected result about the impact of perpendicular flow shear. Specifically, we are primarily interested in the behavior after a few e-folding times, rather than the details of the transients at the very beginning of the evolution. This limit is reached by assuming that $\gamma t\gg 1$ while retaining $|\omega_{V\perp}t|\ll 1$, which implies that $\omega_{V\perp}\ll\gamma$. Thus, the dominant terms are those that are quadratic in time, so the factor of $\left(1\pm\gamma t\right)$ becomes $\pm\gamma t$. Moreover, we are interested in unstable modes (i.e. $\gamma$ is real), for which the $\exp\left(\left(-i\omega+\gamma\right)t\right)$ term will dominate. In these limits, we can write $\phi_{0}+\phi_{1}=C_{\gamma}(1+\epsilon_{\alpha}t^{2})\exp\left(\left(i\omega+\gamma\right)t\right)$, where

[TABLE]

is a complex constant that can be calculated from equation (41) and the expansion of equation (22). To next order in $|\omega_{V\perp}t|\ll 1$, $|\phi_{0}+\phi_{1}|^{2}=|C_{\gamma}|^{2}(1+2~{}\text{Real}(\epsilon_{\alpha})t^{2})\exp\left(2\gamma t\right)$. Thus, the sign of the real part of $\epsilon_{\alpha}$ indicates if perpendicular flow shear will enhance or stabilize the growth of the instability. By manipulating the coefficient, we find that the instability will be enhanced by flow shear if and only if

[TABLE]

We write the condition in this particular form to facilitate comparison with equation (37), the instability growth rate. We see that this condition implies that the instability will be enhanced if the finite flow shear correction to the growth rate is positive. This is an intuitive result and is consistent with the discussion of the $\omega_{V\perp}\ll\gamma$ limit that followed equation (37). If flow shear is moving the mode to a wavenumber that is more strongly driven (relative to damping), then it should begin to grow faster. Note that converting from $\hat{\phi}_{1}$ to $\phi_{1}$ using an expanded version of equation (22) does not end up having an effect on equation (43). This is because it only introduces terms that are either linear in $t$ or imaginary.

The solution we have derived in this section will only be valid for a short time. Our intuition is that a mode will only be unstable for a finite time before flow shear advects them to high radial wavenumbers, where damping is more effective. However, this hypothesis will be rigorously investigated in the next section.

3.2 The $|\omega_{V\perp}t|\gg 1$ limit

We can also obtain analytic solutions to equation (23) in the long time limit. This reveals the ultimate fate of the mode as it is advected to high $k_{x}$ by flow shear, where it can be stabilized by finite sound gyroradius effects. Alternatively, it can be used to understand the behavior of a high $k_{x}$ mode as it is advected towards $k_{x}=0$. In this limit, equation (23) becomes

[TABLE]

This equation is solved by a polynomial that, when substituted into equation (22), gives

[TABLE]

where

[TABLE]

It is important to note that $\tilde{\gamma}$ is a complex number. Additionally, since the time dependence is polynomial rather than exponential as before, the transformation from $\hat{\phi}$ to $\phi$ has been included according to equation (22). The finite gyroradius factor has introduced a factor of $t^{-2}$, while the exponential phase factor has changed the sign of $\tilde{\omega}$. From this solution, we can calculate the condition for an unbounded solution as $t\rightarrow\infty$. This condition is simply $\text{Real}(\tilde{\gamma})>3/2$, which is equivalent to

[TABLE]

or

[TABLE]

at the most unstable parallel wavenumber of $k_{||}c_{S}=k_{y}\rho_{S}\omega_{V||}/2$. If $\omega_{V||}$ is greater than this value, then $\phi\rightarrow\infty$ as $t\rightarrow\infty$, indicating that the stabilizing mechanisms are not sufficiently strong to restrain the PVG instability. This condition is verified against the numerical solution to equation (23) in figure 5. Such a condition is interesting as it indicates what mechanisms are able to saturate PVG turbulence. For weak $\omega_{V\perp}$, the nonlinear term can be expected to limit the modes. For strong $\omega_{V\perp}$, the flow shear on its own is capable of damping the modes, without the nonlinear term.

Additionally, equation (49) can be cast into the form of the approximate “$\omega_{V\perp}\approx\gamma$” quench rule [11], which states that turbulence is quenched when the perpendicular flow shear becomes approximately equal to the maximum linear growth rate $\gamma$ in the absence of perpendicular flow shear. Ignoring the magnetic drifts for simplicity, we can solve equation (37) for $\omega_{V||}$ using the most unstable parallel wavenumber of $k_{||}c_{S}=k_{y}\rho_{S}\omega_{V||}/2$. Then we can substitute this $\omega_{V||}$ in order to rewrite equation (49) as

[TABLE]

For a typical wavenumber of the most unstable mode, $\vec{K}\rho_{S}=(0,0.3)$, we find the condition $\left|\omega_{V\perp}\right|>2.5\left|\gamma\right|$ for turbulent stabilization. Moreover, the fact that stabilization becomes impossible as $k_{y}\rightarrow 0$ makes sense because perpendicular flow shear only affects a mode in proportion to its $k_{y}$ (as seen in equation (8)).

Notice that equation (48) is fairly similar to equation (38), the condition for instability at $t=0$. The primary difference is that the stabilizing perpendicular flow shear term appears and the stability caused by the $\omega_{Mx}$ magnetic drift term slightly weakens. Interestingly, this indicates that for large $K_{x}$, small (but non-zero) $\omega_{V\perp}$, and/or large $\omega_{My}/\omega_{Mx}$, the mode can be initially stable, but unbounded at long time after being advected by flow shear. Moreover, due to the finite sound gyroradius effect contained in equation (22), the mode can actually decay quadratically, but then have an unbounded limit. In fact, as shown in figure 1, the time dependence can be even more complex. Initially, the mode is stable according to equation (37) and its amplitude decreases due to finite sound gyroradius damping as its radial wavenumber changes (i.e. equation (22)). Shortly thereafter, other effects of flow shear become large enough to drive the mode unstable according to equation (41). Then, at intermediate times, both components of the magnetic drifts are important, which reduces the mode growth sufficiently for the finite sound gyroradius damping to win again. Lastly, at long times, the $\omega_{My}$ term becomes negligible, allowing the parallel flow shear to slightly overpower $\omega_{Mx}$ and the finite sound gyroradius damping according to equation (49). Thus, the mode grows polynomially without limit.

3.3 Zonal modes

Lastly, for zonal modes (i.e. $k_{y}=0$), we can directly solve equation (23) and find a simple solution without needing to take any limits. This is because the effects of both parallel and perpendicular flow shear are proportional to $k_{y}$, so they vanish. We are left with a simple 2nd order ordinary differential equation with constant coefficients, just like in section 3.1. Thus, the solution can be found by setting $k_{y}=0$ in equation (35) to get

[TABLE]

We see that the zonal modes are not driven nor damped linearly [49] and oscillate if they have a finite $k_{||}$ and/or there is a magnetic drift in the $x$ direction. Like other gradient-driven instabilities, the PVG cannot drive the zonal modes directly because its effect is proportional to $k_{y}$. Instead the zonal modes must be driven by the nonlinear interaction of pairs of non-zonal modes (the focus of section 4.1).

4 Nonlinear two-dimensional slab results

To investigate the nonlinear dynamics of turbulence, we will let $k_{||}=0$ in order to study two-dimensional turbulence. Strictly speaking, this is inconsistent with our assumption of adiabatic electrons [50]. However, the dynamics should be similar to a very small value of $k_{||}$ and setting $k_{||}$ equal to exactly zero facilitates benchmarking against gyrokinetic codes. In this section, we will focus on how perpendicular flow shear alters the fundamental nonlinear drive of small amplitude “driven” modes by large “pump” modes. The analysis will be kept general for benchmarking purposes, but the primary physics application is zonal flow dynamics. Zonal flows cannot be directly driven by background gradients, so they only grow due to the nonlinear coupling of non-zonal modes. Nevertheless, in many if not most gyrokinetic simulations, the amplitude of the zonal flows eventually become much larger than the non-zonal modes [51, 52]. Thus, despite the fact that zonal modes do not directly cause transport, the mechanism by which they regulate non-zonal modes is of much interest.

With $k_{||}=0$, the evolution equation for $\phi$ (i.e. equation (20)) can be formulated as

[TABLE]

where the coupling condition is $\vec{K}+\vec{K}^{\prime}+\vec{K}^{\prime\prime}=0$. Note that we have redefined $\vec{K}^{\prime}\rightarrow-\vec{K}^{\prime}$ and $\vec{K}^{\prime\prime}\rightarrow-\vec{K}^{\prime\prime}$ and used the reality condition that $\phi(-\vec{K})=\phi^{\ast}(\vec{K})$ (where the $\ast$ superscript indicates the complex conjugate). We see that the dependence on the parallel velocity moment and the parallel velocity gradient has dropped out. This is not too troubling because the parallel flow shear does not enter into the nonlinear term, which is the primary focus of our analysis. We note that equation (52) is a generalization of the Charney-Hasegawa-Mima equation [32] to include the perpendicular flow shear and simple magnetic drifts.

Substituting the hat notation from equation (22) yields

[TABLE]

where we have defined

[TABLE]

Now, as in a three-wave resonant decay calculation [33], we will consider three Fourier modes with $\vec{K}_{a}+\vec{K}_{b}+\vec{K}_{c}=0$ such that they nonlinearly couple. We will consider two separate cases: two pump modes such that $|\phi_{b}|\sim|\phi_{c}|\gg|\phi_{a}|$ and a single pump mode such that $|\phi_{b}|\gg|\phi_{c}|\sim|\phi_{a}|$. Here and henceforth the Latin letter subscript(s) specifies the Fourier modes of the wavenumber arguments. The first case is straightforward to solve analytically and can model the initial drive of the zonal modes from zero amplitude. The second case is more complex, but is important for understanding how dominant zonal modes couple to non-zonal modes.

4.1 Two pump mode case

Linearizing equation (53) in the ratio of the amplitude of the driven mode “a” to the amplitudes of the two pump modes “b” and “c” produces

[TABLE]

where the mode subscripts have been dropped from $\theta$ because their order has no effect. Immediately, we see that the pump modes have the solutions of $\hat{\phi}_{b}=\hat{\phi}_{b}^{\varnothing}$ and $\hat{\phi}_{c}=\hat{\phi}_{c}^{\varnothing}$. This result is consistent with our previous analysis, as it satisfies equation (23) when $k_{||}=0$. With a constant $\hat{\phi}_{b}$, equation (22) tells us that

[TABLE]

and the same is true for mode “c”. Hence, we see that the phase of each pump mode oscillates with a non-constant frequency due to the gradients in the plasma. The only effect that alters the magnitude of the modes is the perpendicular flow shear, which enters through the factor $(1+k_{\perp}^{\varnothing 2}\rho_{S}^{2})/(1+k_{\perp}^{2}\rho_{S}^{2})$. This advects the mode in $k_{x}$ and modifies the finite gyroradius effects that originated in the quasineutrality equation. However, if a pump mode is zonal the wavenumber is not sheared, so that the mode maintains a constant amplitude. Lastly, we note that, because we have linearized our differential equations, our analytic calculations will not exhibit pump depletion (i.e. losses due to the energy transferred to the driven mode). However, pump depletion is negligible as long as the amplitude of the driven mode remains much smaller than the pump modes.

Since $\hat{\phi}_{b}$ and $\hat{\phi}_{c}$ are known, equation (57) can be integrated directly to find the evolution of the driven mode. Substituting the form of $\Lambda_{abc}$ gives

[TABLE]

Thus, we see that for pump modes with fixed amplitudes $\hat{\phi}_{b}$ and $\hat{\phi}_{c}$, the driven mode “a” does not grow exponentially with time, it can only grow polynomially. This is consistent with the numerical results of reference [52], which used a realistic toroidal geometry. Equation (61) also shows that, if the three modes are oscillating out of phase (i.e. $\theta\neq 0$), then the solution will be oscillatory. However, the presence of perpendicular flow shear (i.e. $\omega_{V\perp}\neq 0$) causes $\theta$ to vary with time. Thus, $\theta(t)$ can pass through $\theta=0$, allowing the mode to grow temporarily. When $\omega_{V\perp}\neq 0$ and $\theta\approx 0$, $k_{x}$ varies linearly in time for all non-zonal modes, which means that temporarily $\left|\phi_{a}\right|$ can grow quickly (e.g. if $k^{2}\rho_{S}^{2}\ll 1$, $\left|\phi_{a}\right|$ can grow as $t^{3}$). The condition of $\theta\approx 0$ is equivalent to the frequency matching condition that is required for mode growth in three-mode coupling without background flow shear [33]. For example, when $\omega_{V\perp}=0$, the wavenumbers in equation (61) become constant and $|\phi_{a}|$ can grow linearly with $t$, but only if $\theta=0$.

In the long time limit, the nonlinear drive (i.e. the right side of equation (61)) for a non-zonal mode “a” is proportional to $t^{-1}$. This means that $\left|\phi_{a}\right|$ decays away to zero as $t^{-3}$. However, the situation of most interest is when mode “a” is zonal. When $k_{ay}=0$, the coupling condition requires the other two modes to have the same value of $k_{y}$, so they get sheared at the same rate. Thus, the numerator of equation (61) varies linearly in time instead of quadratically and the drive of a zonal mode decays more quickly than a non-zonal one (specifically as $t^{-2}$ instead of $t^{-1}$). Lastly, when one of the pump modes is zonal, the drive actually gets stronger with time (specifically it scales linearly with $t$).

These analytic solutions are again verified by comparison with the gyrokinetic code GENE as shown in figures 6 and 7. The simulation is initialized with $\hat{\phi}_{b}^{\varnothing}/(B\rho_{S}^{2})=\hat{\phi}_{c}^{\varnothing}/(B\rho_{S}^{2})=c_{S}/L_{r}$, while all other modes are zero. The wavenumbers are chosen to be $\vec{K}_{b}\rho_{S}=(-0.0086,0.015)$ and $\vec{K}_{c}\rho_{S}=(0.01,0.025)$, while $\omega_{V\perp}=-0.002c_{S}/L_{r}$, $\omega_{Mx}=0$, and $\omega_{My}=c_{S}/L_{r}$. Figure 6 shows that the pump mode “b” decreases because flow shear advects it to large $k_{\perp}^{2}$, where the finite sound gyroradius effects more effectively average over it. The pump mode “c” increases briefly because flow shear initially advects it to a lower value of $k_{\perp}^{2}$, but then it starts to decay after passing through $k_{x}=0$. The only other mode that grows is the one with $\vec{K}_{a}=-(\vec{K}_{b}+\vec{K}_{c})$, which is driven by the nonlinear coupling and is shown in figure 7.

Finally, note that the calculation in this section for the initial drive of zonal modes can be extended to more than just three modes. This is because, during the initial linear phase of the drive, the evolution of all the non-zonal modes is determined by the linear physics (i.e. the drive from the background gradients). Moreover, the zonal modes do not directly couple with each other nonlinearly because their wavevectors are parallel to each other such that $\vec{K}^{\prime}\times\vec{K}^{\prime\prime}=0$ in equation (54). This means that the non-zonal modes can be solved without any nonlinearity and the evolution of each zonal mode can be calculated from equation (61) with a summation over all non-zonal modes out in front. Because the drive of zonal modes decays more quickly with time, their drive will tend to be more dominated by the modes near $k_{x}=0$ (compared to the nonlinear drive of non-zonal modes).

4.2 One pump mode case

In this section, instead of having two large amplitude pump modes, we will include just one: $|\phi_{b}|\gg|\phi_{a}|\sim|\phi_{c}|$. This is appropriate for modeling parametric instability decay. While the previous section studied the initial growth of a zonal mode from small amplitude, this section will focus on the opposite situation. In the quasi-steady-state of nonlinearly saturated turbulence, the amplitudes of the zonal modes are typically much larger than that of the non-zonal modes [51, 52]. Thus, we would like to study how a single large-amplitude zonal mode regulates small-amplitude non-zonal modes. However, we will keep the analysis general to allow for a non-zonal pump mode.

Linearizing equation (53) in the ratio of the driven mode amplitudes (here modes “a” and “c”) to the pump mode amplitude (here mode “b”) gives

[TABLE]

The pump mode solution is again given by equation (22) with $\hat{\phi}_{b}=\hat{\phi}_{b}^{\varnothing}$, which is compared to GENE simulations in figure 8. We see good agreement, except for the subtle step-like behavior in the GENE evolution (which can also be seen in figure 6). This is an artifact of the “wavevector-remap” scheme used to implement flow shear and will converge to a continuous evolution in the limit of high radial wavenumber resolution [48].

Note that strictly speaking equations (62) and (64) should have a second term. This is because, due to the reality condition of $\phi(-\vec{K})=\phi^{\ast}(\vec{K})$, the pump mode is actually two modes — one at $\vec{K}_{b}$ and one at $-\vec{K}_{b}$. Thus, there is a second set of three modes involving a new mode $\vec{K}_{d}$ that satisfies the coupling condition $\vec{K}_{a}-\vec{K}_{b}+\vec{K}_{d}=0$. However, this second term can be ignored if it does not lead to instability or if mode “d” is not included in the simulation domain.

Calculating the evolution of the driven modes “a” and “c” requires considerable mathematics. Solving equation (62) for $\hat{\phi}^{\ast}_{c}$ and substituting it into the complex conjugate of equation (64) yields

[TABLE]

where

[TABLE]

and we note that $\left|\phi_{b}\right|^{2}$ can be time dependent. The analogous equation for $\hat{\phi}_{c}$ can be found by swapping the subscripts “a” and “c”. These equations are identical to those in reference [33], except for the presence of the flow shear $\omega_{V\perp}$ introduces $G_{abc}$ and makes $\theta$, $\Lambda_{abc}$, and $\left|\phi_{b}\right|^{2}$ depend on time.

Equation (65) is straightforward to solve numerically as it can be written as a 2nd order, ordinary differential equation with polynomial coefficients. This is done for many different “a” modes in figure 9, which are compared against GENE simulations. We see that some modes are unstable for the entire simulation (e.g. $\vec{K}_{a}\rho_{S}=(-0.15,0.1)$), while others only ever oscillate (e.g. $\vec{K}_{a}\rho_{S}=(-0.35,0.3)$). Certain modes start off unstable, but then the perpendicular flow shear advects them to wavenumbers that are stable (e.g. $\vec{K}_{a}\rho_{S}=(-0.25,0.2)$). However, GENE and our numerical solution do not agree well for several modes. This is because the evolution of these modes is not governed by simple three-wave resonant decay. For these cases, there is a fourth mode in the simulation domain that is unstable, which complicates our analysis as discussed at the beginning of this subsection. These cases have been indicated in figure 9 with gray and should not necessarily agree well. All the other modes match the numerical prediction. The most significant source of error was the “wavevector-remap” scheme used to model perpendicular flow shear, which is only exact in the limit of infinite $k_{x}$ resolution [48]. This error was minimized by using a $k_{x}$ resolution four times denser than the points shown in figure 9. Note the $\vec{K}_{a}\rho_{S}=(-0.1,0.2)$ mode shows good agreement, but the numerical solution is only plotted for a short time because the numerical solution encountered convergence issues.

Equation (65) does not have an analytic solution in general because the polynomials are of very high degree (i.e. the effective coefficient for the second derivative term is degree 10, the first derivative term is degree 9, and the $\hat{\phi}_{a}$ term is degree 6). Nevertheless, we see that if the nonlinear coupling is weak, then $\hat{\phi}_{a}=\text{const}$ is a solution. This is an important and nontrivial solution because the mode will still evolve through the dependences in equation (22). As in the linear analysis of section 3, further analytic results can be found by investigating equation (65) in the $|\omega_{V\perp}t|\ll 1$ and $|\omega_{V\perp}t|\gg 1$ limits.

4.2.1 The $|\omega_{V\perp}t|\ll 1$ limit

To lowest order in the short time limit, equation (65) becomes

[TABLE]

where the numerical subscript indicates the quantity’s order in the $|\omega_{V\perp}t|\ll 1$ expansion. To lowest order, all the coefficients are constant in time and simply equal to their values without perpendicular flow shear. Surprisingly, this differential equation has the same form as the one we solved in section 3.1. Using that solution, we find

[TABLE]

where

[TABLE]

are the mode frequency and growth rate respectively. Inspecting this solution, we see that the mode can be nonlinearly unstable (i.e. have a real growth rate) only if $(k_{b}^{2}-k_{c}^{2})(k_{a}^{2}-k_{b}^{2})>0$, which is equivalent to $k_{a}\leq k_{b}\leq k_{c}$ or $k_{c}\leq k_{b}\leq k_{a}$. Thus, to lowest order we find the traditional nonlinear coupling present in a resonant three wave decay process [33]. This solution is verified against GENE simulations in figure 10. Figure 10 shows that, after a short transient in the GENE simulations, there is good agreement for the mode growth rate. To see the impact of flow shear we must go to next order in the expansion.

To next order, the differential equation for the driven mode becomes

[TABLE]

where

[TABLE]

Flow shear enters through $\theta_{1}$ and $\left(\Lambda_{abc}\Lambda_{cab}\left|\phi_{b}\right|^{2}\right)_{1}$, which depend linearly on time, and $G_{abc1}$, which is constant. This equation can be solved analytically. The homogeneous solution has the same form as the lowest order solution $\hat{\phi}_{a0}$, while the particular solution is

[TABLE]

Thus, we see that the dominant effect of weak flow shear (or alternatively the first effect of flow shear to appear) is a quadratic correction to the lowest order exponential behavior. This is similar to the effect of flow shear on the PVG-driven instability that we found in section 3.1.

However, we notice that this solution, which is fairly complex, can be simplified. Specifically, we are primarily interested in the behavior after a few e-folding times (i.e. $\gamma t\gg 1$ while retaining $|\omega_{V\perp}t|\ll 1$), rather than the details of the transients at the very beginning of the evolution. This limit is equivalent to $\omega_{V\perp}\ll\gamma$ and the dominant terms are those that are quadratic in time, giving

[TABLE]

As in section 3.1, for unstable modes (i.e. $\gamma$ is real), the $\exp\left(\left(i\omega+\gamma\right)t\right)$ term will dominate and the real part of its coefficient indicates if flow shear will enhance or stabilize its growth. Therefore, the instability will be enhanced by perpendicular flow shear if and only if $\left(\Lambda_{abc}\Lambda_{cab}\left|\phi_{b}\right|^{2}\right)_{1}>\theta_{0}\theta_{1}/2$. In other words, the instability is enhanced if the finite flow shear correction to the growth rate (i.e. equation (70)) is positive. This result should be expected according to our previous discussions of the $\omega_{V\perp}\ll\gamma$ limit. If flow shear is moving the mode to a wavenumber that is more strongly driven (relative to damping), then it should begin to grow faster.

4.2.2 The $|\omega_{V\perp}t|\gg 1$ limit

To lowest order in the long time limit, the coefficients in equation (65) reduce to

[TABLE]

We see that the nonlinear coupling term is small, which means that we can return to equations (62) and (64) and solve them on their own. As should be expected the nonlinear coupling term

[TABLE]

is small there too, so that the lowest order solution is simply

[TABLE]

This limit is verified against the full numerical solution in figure 11. This means that the driven modes are not driven at all and simply evolve in isolation according to the $t^{-2}$ dependence contained in equation (22). Thus, they either decay as $t^{-2}$ as they move towards $k_{x}\rightarrow\pm\infty$ or grow as $t^{2}$ as flow shear brings them towards $k_{x}=0$ from large $|k_{x}|$. To next order, the coupling term enters equations (62) and (64), which can be solved to find

[TABLE]

where $t^{\varnothing}$ is the time at which $C_{\lambda}$ and $C_{\mu}$ are determined. The expression for $\hat{\phi}_{c1}$ is the same except the “a” and “c” indexes must be swapped.

From equations (81) and (82), we see that the dominant effect of flow shear at long times is the $t^{-2}$ dependence from the ion polarization drift factor in equation (22). The dominant effect of the nonlinear coupling on the amplitude of a decaying mode is given by equations (22) and (83), which scales as $t^{-3}$.

4.2.3 Zonal pump modes

When the pump mode is zonal, the solution in the $|\omega_{V\perp}t|\ll 1$ limit remains largely the same. However, the behavior in the $|\omega_{V\perp}t|\gg 1$ limit does not. In contrast to a non-zonal pump mode, the nonlinear coupling is not small for a zonal pump mode. Even at long times the nonlinear interaction is still important. To lowest and next order, the coefficients in equation (65) are

[TABLE]

and for the equation to balance we require that $\partial/\partial t\sim\omega_{V\perp}$. This means that once again we must solve equation (65) with constant coefficients (and $G_{abc}=0$). Thus, like in section 4.2.1, the lowest order solution is given by

[TABLE]

where

[TABLE]

The important distinction, which we have made explicit in equation (87), is that the lowest order nonlinear coupling term is necessarily negative. Thus, the nonlinear coupling makes the mode oscillate and modifies the phase of the mode, which can have important consequences on which other modes it can resonate with. However, it does not affect the scaling of the mode amplitude. To lowest order, the dominant effect of flow shear is the $t^{-2}$ factor from the ion polarization drift in equation (22). This scaling is verified in figure 12 for a mode that is decaying in amplitude as flow shear advects it to $k_{x}\rightarrow\infty$. The oscillation frequency predicted by equation (87) was also found to agree.

To the next nontrivial order, $O\left(\omega_{V\perp}^{2}t^{2}\right)$, the differential equation gains time dependent coefficients. It becomes identical to equation (71), except the first order coefficients are proportional to $t^{-2}$ instead of $t$ (and $G_{abc1}=0$). The inhomogeneous solution to this is

[TABLE]

where $E_{I}\left(z\right)\equiv-\int_{-z}^{\infty}dt\exp\left(-t\right)/t$ is the exponential integral function. As before, we can expand in $\gamma t\gg 1$ to get the behavior after a few e-folding times. We find

[TABLE]

which has a purely imaginary exponent and is proportional to $t^{-3}$. Thus, once again, our solution does not affect the amplitude of the decaying mode. This means that, factoring in the $t^{-2}$ dependence from equation (22), the effect of the nonlinear coupling term on the amplitude of a driven mode decays away more rapidly than $t^{-5}$.

Comparing to the previous subsection, we find a surprising result. The long time effect of the nonlinear coupling on the mode amplitude is much weaker when the pump mode is zonal. Specifically, the coupling from a non-zonal pump scales as $t^{-3}$, while the coupling from a zonal pump scales more weakly than $t^{-5}$. This suggests that coupling between non-zonal modes has the dominant nonlinear effect on the amplitude of a mode when it is at large $k_{x}$. Intuitively, it seems like the opposite should be true because the zonal pump does not vary much with time, while a non-zonal pump does. Nevertheless, the coupling with zonal modes still has the dominant effect on the phase and oscillation frequency of a mode and this effect does not get weaker when the mode is at high $k_{x}$.

5 Conclusions

In this work we have taken the cold ion limit of local $\delta f$ gyrokinetics to derive an exact yet simple fluid model (i.e. equations (16) and (17)). This model was applied to a slab geometry that included magnetic drifts, but no magnetic shear. The resulting equations were found to reduce to the Charney-Hasegawa-Mima model in the presence of background flow. Then, an analytic solution (i.e. equation (24)) was found for the full time evolution of a single Fourier mode driven unstable by a parallel velocity gradient, but stabilized by perpendicular flow shear. Studying this solution revealed that quite complicated behavior is possible (e.g. figure 1). Additionally, we calculated simple criteria governing the initial stability (i.e. equation (39)) and long time stability (i.e. equation (49)) of a mode. These two criteria were found to be somewhat different.

The fluid model was also used to study the basic three-mode nonlinear coupling in two-dimensional turbulence. The initial nonlinear drive of zonal flows was investigated using two large pump modes (i.e. equation (61)) and the analytic results were verified against nonlinear gyrokinetic simulations (i.e. figures 6 and 7). Next, the physics occurring in the quasi-steady-state of nonlinearly saturated turbulence was studied using one large pump mode. Specifically, we focused on how a single large zonal mode drives and regulates non-zonal modes. We found that coupling with a zonal mode has an important effect on the phase and oscillation frequency of a mode that does not get weaker as flow shear advects the mode from high $k_{x}$ or to high $k_{x}$. However, coupling with non-zonal modes actually has the dominant effect on the amplitude of a mode. Again, the results were verified against gyrokinetic simulations (i.e. figures 8 through 12).

Acknowledgments

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. Simulations were performed with the support of EUROfusion and MARCONI-Fusion. This work was supported by a grant from the Swiss National Supercomputing Centre (CSCS) under project ID s863. This work was supported in part by the Swiss National Science Foundation.