Impact of poloidal convective cells on momentum flux in tokamaks

Calculating particle, momentum and heat fluxes is of utmost importance in magnetized fusion plasmas. Indeed the confinement time is mainly determined by fluxes across the magnetic surfaces. These fluxes can be computed as velocity momenta of the distribution function for each species. At the microscopic level, the motion of a charged particle in a strong guide field can be separated into a cyclotron and gyrocenter motion. The velocity of particle gyrocenters perpendicular to the magnetic field can itself be split into two components. The first one is the E × B drift velocity, while the second one is the magnetic drift, due to field curvature and gradient of the magnetic field modulus. Hence gyrokinetic fluxes, defined here as fluxes of gyrocenters, contain contributions from the $E\times B$ and magnetic drifts. Exact conservation laws can be derived in this framework [1]. Another approach, may be more intuitive, has been used in the past. It consists in separating the fluxes in collisional and turbulent parts. In toroidal fusion devices, collisional processes are amplified by trajectory effects. The resulting fluxes are called 'neoclassical'. Neoclassical and turbulent fluxes are often calculated separately and added up. Without turbulence, neoclassical fluxes match fluxes due to curvature drift, because the contribution of the E × B drift velocity is small [2, 3]. Cases exist however where the electric drift matters, like the neoclassical flux of toroidal momentum [4, 5], or impurity transport when poloidal asymmetries are large (e.g. due to the centrifugal force) [6]. Conversely, turbulent fluxes are usually computed by keeping only the fluctuating E × B velocity [7]. Hence a common methodology consists in using a code to calculate neoclassical fluxes (i.e. without turbulence), and a gyrokinetic code to compute turbulent fluxes, usually without direct contributions of the magnetic drift to fluxes. Both outputs are then added. This procedure can be justified under conditions of scale separation [8].

Recently several gyrokinetic codes have been upgraded to compute both collisional and turbulent transport, thanks to the implementation of accurate collisional operators. Results have been published for heat [9, 10], momentum [11–16] and particle [17] transport channels. These recent developments have raised the question of defining neoclassical and turbulent fluxes in a proper and practical way, and thus the issue of the interplay between collisional and turbulent processes. One way consists in defining the neoclassical flux as due to magnetic drift only, while turbulent fluxes are due to $E\times B$ drift velocity fluctuations. This procedure finds some justification in the fact that it coincides with the conventional definition of a neoclassical flux in absence of turbulence, under the conditions given above, while the turbulent flux agrees with previous results that discarded the magnetic drift contribution to the flux. However this approach leads to inconsistencies, as shown in the references aforementioned. One reason is turbulent scattering in the velocity space, which plays a role similar to Coulomb collisions [15]. Another reason is the development of large scale poloidal asymmetries of the mean $E\times B$ drift velocity due to turbulence. These asymmetries contribute to the neoclassical flux [9], whereas they are small in conventional neoclassical theory. These effects are sometimes called 'synergistic'. To avoid any misunderstanding and lengthy debates, fluxes are split in the present paper into contributions due to $E\times B$ and curvature drifts, without using the words 'turbulent' and 'neoclassical'. We note that another possible definition would consist in separating fluxes due to large and small scales. Though the latter definition seems now to prevail, previous numerical works used the former one. Given these preliminaries, the question addressed in this paper is the following. Is the flux due to magnetic drift significant in a collisionless turbulent state? Since this is still a broad subject, it is limited here to toroidal momentum transport. This topic has triggered much attention due to the observation of intrinsic toroidal rotation in tokamaks, i.e. rotation without external torque (see [18] for an overview).

In several previous works, only the $E\times B$ drift contribution to the parallel (or toroidal) momentum flux was calculated [14, 19–21]. However it was stressed in [22] that the contribution of the magnetic drift to the flux of momentum could be significant in non linear regime. Moreover two research teams recently computed both magnetic and $E\times B$ drift contributions to the momentum flux. Both found that these two contributions are anti-correlated [13, 16], leading in some cases to an almost vanishing total flux. A mechanism is proposed here to explain the correlation between the contributions of magnetic and $E\times B$ drifts to momentum transport, based on turbulent generation of poloidal convective cells. Poloidal convective cells are zonal structures with zero toroidal wavenumber, and finite poloidal wave numbers. We will make the difference with zonal flows by restricting the definition of zonal flows to structures with zero toroidal and poloidal wave numbers. The chain of causes starts with the well established correlation between $E\times B$ radial fluxes of parallel and poloidal momentum due to mean shear flow or any other mechanism responsible for an up–down asymmetry of turbulence intensity. Turbulence is ballooned in a tokamak, i.e. is more intense on the low field side. The resulting asymmetry of the turbulent stress tensor generates poloidal convective cells, i.e. flows that are zonal in the toroidal rotation, but not in the poloidal direction. It appears that these flows are weakly damped at low frequency, and drive poloidal asymmetries of the distribution function, which contribute to a non zero flux of parallel momentum due to magnetic drift. Hence it appears that the magnetic drift component of the flux of parallel momentum is tied to the $E\times B$ radial flux of poloidal momentum. Anti-correlation of the poloidal and parallel components of the turbulent Reynolds stress leads to anti-correlated radial fluxes of parallel momentum due to magnetic and $E\times B$ drifts. This scheme is illustrated in figure 1. A schematic view of the circulation pattern in a tokamak poloidal plane is also shown in figure 2. The interplay between low wavenumber convective cells and turbulent transport has already been mentioned in the context of internal transport barrier formation [23].

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

A simplified calculation is presented here, where the fluid parallel velocity and its gradient are small. In other words only the residual stress tensor is calculated. It appears that the sign of the correlation between the poloidal and parallel components of the $E\times B$ Reynolds stress translates in a correlation of the curvature and $E\times B$ drift contributions to the radial flux of parallel momentum. This general result comes from a simple relationship between the momentum flux due to curvature and the turbulent Reynolds stress (summarized in figure 2). Although poloidal convective cells do not appear explicitly in this relationship, they are instrumental to this mechanism. The overall process is illustrated by a quasi-linear calculation of stress tensors and fluxes based on linear ion temperature gradient (ITG) driven modes in the hydrodynamic limit. This calculation predicts similar amplitudes for the $E\times B$ and magnetic drift components of the parallel momentum flux, but a sign of correlation that depends on plasma parameters. Typically for modes drifting in the ion diamagnetic direction, positive correlation is found for weak drive, and anti-correlation for strong drive. Since the hydrodynamic limit is valid only in the strong drive limit, this result can be considered as encouraging in view of previous numerical calculations, which found anti-correlation [13, 16].

The paper is organized as follows. General expressions of radial fluxes of parallel momentum are derived in section 2. The contribution from the magnetic drift is detailed in section 3. The $E\times B$ and curvature driven fluxes of parallel momentum are compared and discussed in section 4. A conclusion follows. Most technical demonstrations can be found in the appendices, and may be skipped on first reading.

Momentum flux is the same as the Reynolds stress tensor up to a mass density Nm, where N is the unperturbed density and m the ion mass. Both names will be used indistinctly throughout the paper. We calculate separately the contributions from the E × B drift and magnetic drift velocities to the flux of parallel momentum, i.e.

$$\begin{eqnarray}&&{{\boldsymbol{\Pi }}}_{\parallel }^{E}=\displaystyle \frac{1}{N}\displaystyle \int {{\rm{d}}}^{3}{\bf{v}}{{Fv}}_{\parallel }{{\bf{v}}}_{E}\end{eqnarray} \tag{ 1 }$$

and

$$\begin{eqnarray}&&{{\boldsymbol{\Pi }}}_{\parallel }^{D}=\displaystyle \frac{1}{N}\displaystyle \int {{\rm{d}}}^{3}{\bf{v}}{{Fv}}_{\parallel }{{\bf{v}}}_{D},\end{eqnarray} \tag{ 2 }$$

where F is the distribution function, ${v}_{\parallel }$ is velocity that is aligned with the unperturbed magnetic field, ${{\bf{v}}}_{E}$ is the E × B drift velocity and ${{\bf{v}}}_{D}$ is the magnetic drift velocity due to curvature and gradient of the unperturbed magnetic field. The analysis is restricted to electrostatic turbulence, i.e.

$$\begin{eqnarray}&&{{\bf{v}}}_{E}={\bf{b}}\times \displaystyle \frac{{\rm{\nabla }}\phi }{B},\end{eqnarray} \tag{ 3 }$$

where ϕ is the electric potential, ${\bf{b}}=\tfrac{{\bf{B}}}{B}$ and ${\bf{B}}$ is the (unperturbed) magnetic field. Also we use a simplified expression of the magnetic drift velocity, valid at low values of the plasma β and normalized gyroradius ${\rho }_{* }$,

$$\begin{eqnarray}&&{{\bf{v}}}_{D}=\displaystyle \frac{m}{{eB}}\left(\displaystyle \frac{{v}_{\perp }^{2}}{2}+{v}_{\parallel }^{2}\right){\bf{b}}\times \displaystyle \frac{{\rm{\nabla }}B}{B},\end{eqnarray} \tag{ 4 }$$

where m is the mass, e the charge, and ${v}_{\perp }$ is the modulus of the perpendicular velocity.

A simplified geometry of circular concentric magnetic surfaces is considered here. The spatial coordinates are $(r,\theta ,\varphi )$, where r is the minor radius, θ and φ the poloidal and toroidal angles (see figure 2) and the magnetic field is

$$\begin{eqnarray}&&{\bf{B}}={B}_{0}\displaystyle \frac{r}{q(r)}{\rm{\nabla }}(\varphi -q(r)\theta )\times {\rm{\nabla }}r,\end{eqnarray} \tag{ 5 }$$

where ${R}_{0}$ the major radius, $q(r)$ the safety factor, and B₀ the magnetic field on the magnetic axis. The inverse aspect ratio $r/{R}_{0}$ is a small parameter throughout the paper.

2.1. Turbulence parametrization

2.1.1. Electric potential fluctuations

Since an evaluation of the turbulent Reynolds stress tensor is needed, it is necessary to postulate a spatial structure of the fluctuations of the electric potential. We follow a rationale that is close in spirit to the ballooning representation [24–29]. In a magnetized plasma with strong guide field, turbulent structures tend to align with the magnetic field. Given the structure of the magnetic field equation (5), it is convenient to use the variable $\varphi -q(r)\theta$ instead of the toroidal angle. The periodicity in the toroidal direction allows using a a Fourier series

$$\begin{eqnarray}&&\phi (r,\theta ,\varphi ,t)=\displaystyle \sum _{n}{\tilde{\phi }}_{n}(r,\theta ,t)\exp [{\rm{i}}{n}(\varphi -q(r)\theta )],\end{eqnarray} \tag{ 6 }$$

where each Fourier component ${\tilde{\phi }}_{n}(r,\theta ,t)$ can be written as a Fourier integral

$$\begin{eqnarray}&&{\tilde{\phi }}_{n}(r,\theta ,t)={\int }_{-\infty }^{+\infty }\displaystyle \frac{{\rm{d}}{{\rm{\Theta }}}_{k}}{2\pi }{\tilde{\phi }}_{n}({{\rm{\Theta }}}_{k},\theta ,t)\exp [{\rm{i}}{{\rm{\Theta }}}_{k}{nq}(r)]\end{eqnarray} \tag{ 7 }$$

valid locally near a reference magnetic surface. The periodicity condition in the poloidal direction imposes ${\tilde{\phi }}_{n}({{\rm{\Theta }}}_{k}+2\pi ,\theta +2\pi ,t)={\tilde{\phi }}_{n}({{\rm{\Theta }}}_{k},\theta ,t)$. In linear ballooning theory, the function ${\tilde{\phi }}_{n}({{\rm{\Theta }}}_{k},\theta ,t)$ is found to be separable in ${{\rm{\Theta }}}_{k}$ and θ for each n and pulsation ω, i.e. ${\tilde{\phi }}_{n}({{\rm{\Theta }}}_{k},\theta ,t)={{\rm{\Phi }}}_{n\omega }({{\rm{\Theta }}}_{k}){\hat{\phi }}_{n\omega }(\theta ){{\rm{e}}}^{-{\rm{i}}\omega t}$. Separability is not always rigorously demonstrated, but is a good proxy. The amplitude ${{\rm{\Phi }}}_{n}$ is called envelope, usually a localized around an angle ${\theta }_{k}$ called ballooning angle, with a narrow width of the order of ${({L}_{p}{\rho }_{i})}^{-1/2}$, where L_p is a gradient length and ${\rho }_{i}$ the ion thermal gyroradius radius. The ballooning angle ${\theta }_{k}$ measures the departure from up–down symmetry on a given field line. Non zero values of ${\theta }_{k}$ result from symmetry breaking mechanisms, as will be seen later on. This leads to a slow radial dependence of ${\tilde{\phi }}_{n}(r,\theta ,t)$ in r. Hence 'fast' radial variations of ${\tilde{\phi }}_{n}(r,\theta ,t)$ come from the radial dependence of the safety factor $q(r)$ in the exponential argument of equation (6), with a characteristic scale $1/\left(n\tfrac{{\rm{d}}{q}}{{\rm{d}}{r}}\right)$, whereas 'slow' radial variations come from the amplitude ${\tilde{\phi }}_{{\bf{k}}}$, with a mesoscale length ${({L}_{p}{\rho }_{i})}^{1/2}$. We insist here on the fact that a 'slow' radial scale is not a mean gradient length—the only constraint is that it should be larger than the 'fast' scale $1/\left(n\tfrac{{\rm{d}}{q}}{{\rm{d}}{r}}\right)$.

Calculating a mode envelope in turbulent regime is questionable since a global mode has presumably no enough time to form before turbulent decorrelation occurs, except very close to the threshold. Nevertheless simulations show that an instantaneous ballooning angle can still be identified [16]. We therefore adopt the following representation in non linear regime

$$\begin{eqnarray}&&\phi (r,\theta ,\varphi ,t)=\displaystyle \sum _{{\bf{k}}}{\tilde{\phi }}_{{\bf{k}}}(\theta ,t)\exp \{{\rm{i}}{n}[\varphi -q(r)(\theta -{\theta }_{k})]\}.\end{eqnarray} \tag{ 8 }$$

The vector ${\bf{k}}$ is a label for the couple $(n,{\theta }_{k})$ so that the sum over ${\bf{k}}$ designates a summation over the toroidal wave number n and also an integral over a distribution of angles ${\theta }_{k}$. The value of the angle ${\theta }_{k}$ depends on time, as shown by simulations [16]. To some extent, this time dependence can be represented by such a distribution. Quantities labeled with a index ${\bf{k}}$ will be called 'Fourier' harmonics, since ${\theta }_{k}$ can be seen as the Fourier counterpart of the radial variable ${nq}(r)$. Each harmonic ${\tilde{\phi }}_{{\bf{k}}}(\theta ,t)$ depends also slowly on the radius r. We will use in the following a 'slow' radial variable ${r}_{0}=\varepsilon r$, where ε is a small parameter. All quantities depends on ${r}_{0}$ with a typical scale that is larger than a typical turbulent vortex size, but is smaller than a gradient length. The label ${r}_{0}$ is omitted to simplify the notations, unless specified otherwise. Equation (8) is in principle not acceptable since it is not periodic in θ. Nevertheless it is a reasonable proxy for a turbulence localized on the low field side of a tokamak (i.e. maximum near $\theta =0$) since ${\tilde{\phi }}_{{\bf{k}}}(\theta ,t)$ is then small for $\theta =\pm \pi$. Hence each Fourier component ${\tilde{\phi }}_{{\bf{k}}}(\theta ,t)$ contains the information on the poloidal localization of fluctuations, and therefore poloidal asymmetries. We postulate a Gaussian form for i.e. ${\tilde{\phi }}_{{\bf{k}}}(\theta ,t)$, exact in the linear framework

$$\begin{eqnarray}&&{\tilde{\phi }}_{{\bf{k}}}(\theta ,t)={\tilde{\phi }}_{{\bf{k}}}(t){a}_{{\bf{k}}}\exp \left\{-\displaystyle \frac{1}{2}{\alpha }_{{\bf{k}}}{(\theta -{\lambda }_{{\bf{k}}}{\theta }_{k})}^{2}\right\},\end{eqnarray} \tag{ 9 }$$

where ${\alpha }_{{\bf{k}}}$ and ${\lambda }_{{\bf{k}}}$ are complex numbers, and ${a}_{{\bf{k}}}$ is a normalizing factor (${\alpha }_{{\bf{k}}}$, ${\lambda }_{{\bf{k}}}$ and ${a}_{{\bf{k}}}$ depend on ${r}_{0}$).

2.1.2. Ballooning angle

The value of the ballooning angle ${\theta }_{k}$ is not easy to determine. In the ballooning formalism, it comes from a calculation of the mode envelope. However, as already mentioned, this calculation is questionable for a turbulent state. An estimate can be found by using a rapid distortion theory [30]. Indeed a structure with initially a zero ballooning angle ${\theta }_{k}=0$ will acquire an effective radial wave number due to a shear flow that is equal to ${\theta }_{k}=-\tfrac{1}{{s}_{0}}{V}_{E}^{\prime }t$ after a time t. Here ${V}_{E}^{\prime }$ is the shear flow rate, i.e. the radial derivative of the mean $E\times B$ velocity. Defining a correlation time ${\tau }_{c}$, one gets the following estimate

$$\begin{eqnarray}&&{\theta }_{k}=-\displaystyle \frac{1}{{s}_{0}}{V}_{E}^{\prime }{\tau }_{c}.\end{eqnarray} \tag{ 10 }$$

Other expressions of ${\theta }_{k}$ have been given in the past [26–29, 31], with similar structures as equation (10) when envelopes are spatially localized. The exact expression of ${\theta }_{k}$ has no real importance here, since we are interested in the relative signs and amplitudes of the momentum fluxes. Sources of symmetry breaking different from shear flow will lead to different expressions of ${\theta }_{k}$, but will cause correlated radial and parallel wavenumbers.

2.2. E × B drift stress tensor

We estimate now the radial fluxes of poloidal and parallel momentum due to $E\times B$ drift, which are proportional to the corresponding components of the Reynolds stress

$$\begin{eqnarray}&&{{\rm{\Pi }}}_{r\theta }^{E}(\theta ,t)=\displaystyle \frac{1}{N}{\displaystyle \int }_{0}^{2\pi }\displaystyle \frac{{\rm{d}}\varphi }{2\pi }\displaystyle \int {{\rm{d}}}^{3}{\bf{v}}{{Fv}}_{{Er}}{v}_{E\theta },\end{eqnarray} \tag{ 11 }$$

$$\begin{eqnarray}&&{{\rm{\Pi }}}_{r\parallel }^{E}(\theta ,t)=\displaystyle \frac{1}{N}{\displaystyle \int }_{0}^{2\pi }\displaystyle \frac{{\rm{d}}\varphi }{2\pi }\displaystyle \int {{\rm{d}}}^{3}{\bf{v}}{{Fv}}_{{Er}}{v}_{\parallel }.\end{eqnarray} \tag{ 12 }$$

The distribution function is decomposed in the same way as the potential

$$\begin{eqnarray}&&F(r,\theta ,\varphi ,{v}_{\parallel },\mu ,t)=\displaystyle \sum _{{\bf{k}}}{F}_{{\bf{k}}}(\theta ,{v}_{\parallel },\mu ,t)\exp \{{\rm{i}}{n}[\varphi -q(r)(\theta -{\theta }_{k})]\},\end{eqnarray} \tag{ 13 }$$

where $\mu =\tfrac{{{mv}}_{\perp }^{2}}{2B}$ is the adiabatic invariant. The equations equations (11), (12) can then be recast as

$$\begin{eqnarray}&&{{\rm{\Pi }}}_{r\theta }^{E}(\theta ,t)=-\displaystyle \sum _{{\bf{k}}}{{ \mathcal A }}_{{\bf{k}}}(\theta ){s}_{0}(\theta -{\theta }_{k}){| {v}_{E{\bf{k}}}(t)| }^{2},\end{eqnarray} \tag{ 14 }$$

$$\begin{eqnarray}&&{{\rm{\Pi }}}_{r\parallel }^{E}(\theta ,t)=\displaystyle \sum _{{\bf{k}}}{{ \mathcal A }}_{{\bf{k}}}(\theta ){v}_{E{\bf{k}}}^{* }(t){V}_{\parallel {\bf{k}}}(t),\end{eqnarray} \tag{ 15 }$$

where

$$\begin{eqnarray}&&{V}_{\parallel {\bf{k}}}=\displaystyle \frac{1}{N}\displaystyle \int {{\rm{d}}}^{3}{\bf{v}}{F}_{{\bf{k}}}{v}_{\parallel }\end{eqnarray} \tag{ 16 }$$

is the Fourier component of the parallel fluid velocity and

$$\begin{eqnarray}&&{v}_{E{\bf{k}}}=-{{\rm{i}}{K}}_{\theta }\displaystyle \frac{{\tilde{\phi }}_{{\bf{k}}}}{{B}_{0}}\end{eqnarray} \tag{ 17 }$$

is the Fourier component of the radial component of the $E\times B$ drift velocity. The wavenumber ${K}_{\theta }$ is a reference poloidal wave number defined as ${K}_{\theta }=-\tfrac{{nq}({r}_{0})}{{r}_{0}}$. The poloidal structure function ${{ \mathcal A }}_{{\bf{k}}}(\theta )$ is defined as (see details in appendix A)

$$\begin{eqnarray}&&{{ \mathcal A }}_{{\bf{k}}}(\theta )=\mathop{\mathrm{lim}}\limits_{{\theta }_{k}\to 0}\displaystyle \frac{{| {\tilde{\phi }}_{{\bf{k}}}(\theta ,t)| }^{2}}{{| {\phi }_{{\bf{k}}}(t)| }^{2}}\end{eqnarray} \tag{ 18 }$$

with the normalization

$$\begin{eqnarray*}&&{\int }_{-\infty }^{\infty }{\rm{d}}\theta {{ \mathcal A }}_{{\bf{k}}}^{2}(\theta )=1.\end{eqnarray*}$$

The amplitude ${{ \mathcal A }}_{{\bf{k}}}$ does not depend on time for the structure equation (8). Note that at this stage, the fluxes depend on the poloidal angle and time (plus a slow radial variation in ${r}_{0}$). A rough estimate of ${V}_{\parallel {\bf{k}}}$ is obtained by a rapid distortion argument [30], i.e. ${V}_{\parallel {\bf{k}}}\simeq -{{\rm{i}}{K}}_{\parallel }\tfrac{e}{m}{\tilde{\phi }}_{{\bf{k}}}{\tau }_{{\bf{k}}}$, where ${\tau }_{{\bf{k}}}$ is a correlation time. The components of the Reynolds stress can be expanded as Fourier series in the poloidal direction

$$\begin{eqnarray}&&{{\rm{\Pi }}}_{r\theta }^{E}(\theta ,t)=\displaystyle \sum _{{\ell }=-\infty }^{+\infty }{{\rm{\Pi }}}_{r\theta ,{\ell }}^{E}(t){{\rm{e}}}^{{\rm{i}}{\ell }\theta },\end{eqnarray} \tag{ 19 }$$

$$\begin{eqnarray}&&{{\rm{\Pi }}}_{r\parallel }^{E}(\theta ,t)=\displaystyle \sum _{{\ell }=-\infty }^{+\infty }{{\rm{\Pi }}}_{r\parallel ,{\ell }}^{E}(t){{\rm{e}}}^{{\rm{i}}{\ell }\theta }.\end{eqnarray} \tag{ 20 }$$

Using the expressions equations (23), (24), one finds

$$\begin{eqnarray}&&{{\rm{\Pi }}}_{r\theta ,{\ell }}^{E}(t)=-\displaystyle \sum _{{\bf{k}}}{C}_{r,{\bf{k}}{\ell }}{| {v}_{E{\bf{k}}}(t)| }^{2}{\theta }_{k},\end{eqnarray} \tag{ 21 }$$

$$\begin{eqnarray}&&{{\rm{\Pi }}}_{r\parallel ,{\ell }}^{E}(t)=\displaystyle \sum _{{\bf{k}}}{C}_{\parallel ,{\bf{k}}{\ell }}{\omega }_{c}{\tau }_{{\bf{k}}}{| {v}_{E{\bf{k}}}(t)| }^{2}{\theta }_{k},\end{eqnarray} \tag{ 22 }$$

where ${\omega }_{c}=\tfrac{{{eB}}_{0}}{m}\gt 0$ is the cyclotron pulsation. The coefficients ${C}_{r,{\bf{k}}{\ell }}$ and ${C}_{\parallel ,{\bf{k}}{\ell }}$ measure the distortion of turbulent structures in the radial and parallel directions. They can be formally written as (see appendix A)

$$\begin{eqnarray}&&{C}_{r,{\bf{k}}{\ell }}=\displaystyle \frac{{K}_{r,{\bf{k}}{\ell }}}{{K}_{\theta }{\theta }_{k}}\end{eqnarray} \tag{ 23 }$$

$$\begin{eqnarray}&&{C}_{\parallel ,{\bf{k}}{\ell }}=\displaystyle \frac{{K}_{\parallel ,{\bf{k}}{\ell }}}{{K}_{\theta }{\theta }_{k}}\end{eqnarray} \tag{ 24 }$$

and

$$\begin{eqnarray}{K}_{r,{\bf{k}}{\ell }} & = & -{\rm{i}}\displaystyle \frac{{\int }_{-\infty }^{\infty }{\rm{d}}\theta {{\rm{e}}}^{-{\rm{i}}{\ell }\theta }\tfrac{\partial {\tilde{\phi }}_{n}}{\partial r}{\tilde{\phi }}_{n}^{* }}{{\int }_{-\infty }^{\infty }{\rm{d}}\theta {| {\tilde{\phi }}_{n}| }^{2}}\\ {K}_{\parallel ,{\bf{k}}{\ell }} & = & -{\rm{i}}\displaystyle \frac{{\int }_{-\infty }^{\infty }{\rm{d}}\theta {{\rm{e}}}^{-{\rm{i}}{\ell }\theta }\tfrac{1}{{{qR}}_{0}}\tfrac{\partial {\tilde{\phi }}_{n}}{\partial \theta }{\tilde{\phi }}_{n}^{* }}{{\int }_{-\infty }^{\infty }{\rm{d}}\theta {| {\tilde{\phi }}_{n}| }^{2}}\end{eqnarray} \tag{ 25 }$$

The wave number ${\ell }$ characterizes fluctuation poloidal symmetries. The useful values of ${\ell }$ are low, typically ${\ell }=0,{\ell }=\pm 1,...$, and therefore different from the much larger turbulent poloidal wave numbers. The coefficients ${K}_{r,{\bf{k}}{\ell }}$ and ${K}_{\parallel ,{\bf{k}}{\ell }}$ can be seen as effective radial and parallel wavenumbers for a given set $({\bf{k}},{\ell })$.

2.2.1. Wave numbers

Calculations in appendix A show that the coefficients ${C}_{r,{\bf{k}}{\ell }}$ and ${C}_{\parallel ,{\bf{k}}{\ell }}$ are close to their ${\ell }=0$ components for small ballooning angles ${\theta }_{k}\ll 1$, i.e.

$$\begin{eqnarray}&&{C}_{r,{\bf{k}}{\ell }}={C}_{r{\bf{k}}}+o({{\ell }}^{2}{\theta }_{k}^{2}),\end{eqnarray} \tag{ 26 }$$

$$\begin{eqnarray}&&{C}_{\parallel ,{\bf{k}}{\ell }}={C}_{\parallel {\bf{k}}}+o({{\ell }}^{2}{\theta }_{k}^{2}),\end{eqnarray} \tag{ 27 }$$

where ${C}_{r{\bf{k}}}$ and ${C}_{\parallel {\bf{k}}}$ are short notations for ${C}_{r,{\bf{k}}0}$ and ${C}_{\parallel ,{\bf{k}}0}$, i.e. ${\ell }=0$ components of the coefficients ${C}_{r,{\bf{k}}{\ell }}$ and ${C}_{\parallel ,{\bf{k}}{\ell }}$. These coefficients are related to the parameters used for the turbulent field representation

$$\begin{eqnarray}&&{C}_{r{\bf{k}}}={s}_{0}\displaystyle \frac{{\mathfrak{R}}({\alpha }_{{\bf{k}}}{\delta }_{{\bf{k}}})}{{\mathfrak{R}}({\alpha }_{{\bf{k}}})},\end{eqnarray} \tag{ 28 }$$

$$\begin{eqnarray}&&{C}_{\parallel {\bf{k}}}=\displaystyle \frac{1}{{q}_{0}{K}_{\theta }{R}_{0}}{| {\alpha }_{{\bf{k}}}| }^{2}\displaystyle \frac{{\mathfrak{I}}({\delta }_{{\bf{k}}})}{{\mathfrak{R}}({\alpha }_{{\bf{k}}})},\end{eqnarray} \tag{ 29 }$$

where

$$\begin{eqnarray}&&{s}_{0}={\displaystyle \frac{r}{q}\displaystyle \frac{{\rm{d}}{q}}{{\rm{d}}{r}}\right|}_{r={r}_{0}}\end{eqnarray} \tag{ 30 }$$

is the magnetic shear calculated at $r={r}_{0}$, and ${\delta }_{{\bf{k}}}={\lambda }_{{\bf{k}}}-1$. The symbol ${\mathfrak{R}}$ (resp. ${\mathfrak{I}}$) designates the real (resp. imaginary) part of a complex number.

In view of the field structure equation (9), the real part of ${\alpha }_{{\bf{k}}}$ must be positive, i.e. ${\mathfrak{R}}({\alpha }_{{\bf{k}}})\geqslant 0$, since modes are spatially localized in the poloidal direction. As expected, the average radial and parallel wave numbers are correlated and proportional to ${\theta }_{k}$, which measures the strength of the symmetry breaking mechanism. However a close inspection of equations (28), (29) indicates that the proportionality coefficients are not the same. For instance, in the limit ${\mathfrak{R}}({\alpha }_{{\bf{k}}})\gg {\mathfrak{I}}({\alpha }_{{\bf{k}}})$ (non propagative localized mode), it appears that ${C}_{r{\bf{k}}}$ is of the sign of ${\mathfrak{R}}({\delta }_{{\bf{k}}}){s}_{0}$, while ${C}_{\parallel {\bf{k}}}$ is of the sign of ${\mathfrak{I}}({\delta }_{{\bf{k}}}){K}_{\theta }{q}_{0}$. Hence the sign and amplitude of the correlation between the radial and parallel wave numbers is by no way trivial. This relationship has been discussed in [32], in the context of turbulent momentum transport in slab geometry. The components ${C}_{r,{\bf{k}}{\ell }}$ and ${C}_{\parallel ,{\bf{k}}{\ell }}$ are even functions of ${\ell }$, so that ${{\rm{\Pi }}}_{r\theta }^{E}(\theta ,t)$ and ${{\rm{\Pi }}}_{r\parallel }^{E}(\theta ,t)$ are even functions of θ. Here calculations are restricted to order one in ${\theta }_{k}$.

Since we are interested in signs, a rapid distortion theory may not be accurate enough. A quasi-linear theory provides a more precise value of the time and poloidal average of ${{\rm{\Pi }}}_{r\parallel }^{E}$. The calculation is done in appendix B, and confirms the estimate equation (22). The structure of the later is in line with previous calculations of the residual stress (see [33–35], and overviews [19, 21]). The structure of the flux of poloidal momentum equation (21) is well-known [36] and was used abundantly in the context of transport barrier formation. We note that ${{\rm{\Pi }}}_{r\theta ,{\ell }}^{E}$ and ${{\rm{\Pi }}}_{r\parallel ,{\ell }}^{E}$ are anti-correlated when ${C}_{r,{\bf{k}}{\ell }}$ and ${C}_{\parallel ,{\bf{k}}{\ell }}$ are of the same sign.

The expression equation (2) of the radial flux of parallel momentum due to magnetic drift shows that a finite flux can only be due to up–down symmetries of the distribution function. The rationale here is as follows: because turbulence is ballooned, the stress tensor is ballooned too, thus leading to poloidal asymmetries of the poloidal flow, and therefore of the electric potential. The resulting distorted distribution function is correlated with the $E\times B$ Reynolds stress, and therefore with the $E\times B$ radial flux of parallel momentum.

3.1. Generation of $E\times B$ poloidal convective cells

Damping of time dependent and poloidally asymmetric flows is not primarily due to collisions, but rather wave particle resonant effects. An estimate can be obtained by solving the linear gyrokinetic equation (see [37, 38] for overviews on gyrokinetic theory)

$$\begin{eqnarray}&&\displaystyle \frac{\partial G}{\partial t}+{{\bf{v}}}_{E}\cdot {\rm{\nabla }}G+{{\bf{v}}}_{D}\cdot {\rm{\nabla }}G+{v}_{\parallel }{{\rm{\nabla }}}_{\parallel }G={F}_{M}\left(\displaystyle \frac{\partial }{\partial t}+{{\bf{v}}}_{* }\cdot {\rm{\nabla }}\right)\displaystyle \frac{e}{T}({ \mathcal J }\cdot \phi )\end{eqnarray} \tag{ 31 }$$

coupled to the Poisson equation

$$\begin{eqnarray}&&-{\rm{\nabla }}\cdot \left(\displaystyle \frac{{Nm}}{{B}^{2}}{{\rm{\nabla }}}_{\perp }\right)\phi +\displaystyle \frac{{{Ne}}^{2}}{T}(\phi -{\langle \phi \rangle }_{\psi })=e\displaystyle \int {{\rm{d}}}^{3}{\bf{v}}{{ \mathcal J }}^{-1}\cdot (F-{F}_{M}),\end{eqnarray} \tag{ 32 }$$

where N is the unperturbed density, ${ \mathcal J }$ the gyroaverage operator and ${{ \mathcal J }}^{-1}$ its inverse, ${{\bf{v}}}_{* }=\tfrac{T}{{{eB}}_{0}}\tfrac{\partial \mathrm{ln}{F}_{M}}{\partial r}$ the kinetic diamagnetic velocity, ${F}_{M}$ is the unshifted Maxwellian distribution function built with the density N and temperature T, and $G=F-{F}_{M}+{F}_{M}\tfrac{e\phi }{T}$ is the non adiabatic part of the distribution function. The distribution functions F and G are functions of the gyrocenter-center position, the adiabatic invariant μ and parallel velocity ${v}_{\parallel }$. The parallel non linear term has been neglected. The mirror force is also ignored since the quantity of interest is the ion parallel flux, to which passing particles mostly contribute. This nevertheless sets a lower bound in frequency, namely the trapped ion bounce frequency of the order of $\tfrac{{v}_{{Ti}}}{{q}_{0}{R}_{0}}\sqrt{{r}_{0}/{R}_{0}}$ . It is reminded that the inverse aspect ratio ${r}_{0}/{R}_{0}$ is a small parameter in the present work. To simplify the calculation, we use the long wavelength limit of the gyroaverage operator ${ \mathcal J }\to 1$, i.e. the polarization term is kept, but not the finite Larmor radius (FLR) effects. For non zonal modes ${\langle \phi \rangle }_{\psi }=0$, this gives an equation for the potential vorticity $\phi -{\rho }_{i}^{2}{{\rm{\nabla }}}_{\perp }^{2}\phi$, i.e.

$$\begin{eqnarray}&&\quad N\left(\displaystyle \frac{\partial }{\partial t}+{{\bf{v}}}_{E}\cdot {\rm{\nabla }}\right)(1-{\rho }_{i}^{2}{{\rm{\nabla }}}_{\perp }^{2})\displaystyle \frac{e\phi }{T}+{{\rm{\nabla }}}_{\parallel }{{\rm{\Gamma }}}_{\parallel }+{\rm{\nabla }}\cdot (N{{\bf{v}}}_{E}+N{{\bf{v}}}_{* p})=0,\end{eqnarray} \tag{ 33 }$$

where ${{\rm{\Gamma }}}_{\parallel }$ is the ion parallel flux, and ${{\bf{v}}}_{* p}$ the diamagnetic velocity. The vorticity equation (33) has been derived by assuming that the electric potential wavelength is smaller than the gradient lengths of density and temperature. The axisymmetric electric potential n = 0 is Fourier expanded in poloidal angle and time

$$\begin{eqnarray}&&{\phi }_{n=0}(\theta ,t)=\displaystyle \sum _{{\ell }\omega }{\phi }_{{\ell }\omega }\exp \{{\rm{i}}({\ell }\theta -\omega t)\},\end{eqnarray} \tag{ 34 }$$

where ${\ell }\ne 0$ since zonal flows are excluded. One can show the following variant of the Taylor identity [39] (see appendix C)

$$\begin{eqnarray}&&\displaystyle \frac{1}{{B}_{0}}{\int }_{0}^{2\pi }\displaystyle \frac{{\rm{d}}\varphi }{2\pi }({{\bf{v}}}_{E}\cdot {\rm{\nabla }}){{\rm{\nabla }}}_{\perp }^{2}\phi ={{\rm{\nabla }}}_{\perp }^{2}{{\rm{\Pi }}}_{r\theta }^{E},\end{eqnarray} \tag{ 35 }$$

where

$$\begin{eqnarray}&&{{\rm{\nabla }}}_{\perp }^{2}{{\rm{\Pi }}}_{r\theta }^{E}=\displaystyle \frac{1}{{r}_{0}}\displaystyle \frac{\partial }{\partial {r}_{0}}\left({r}_{0}\displaystyle \frac{\partial }{\partial {r}_{0}}{{\rm{\Pi }}}_{r\theta }^{E}\right).\end{eqnarray} \tag{ 36 }$$

It is reminded that ${r}_{0}$ is the reference radius in the neighborhood of which the calculations are performed. Keeping FLR effects would replace one component of the drift velocity ${{\bf{v}}}_{E}$ by ${{\bf{v}}}_{E}+{{\bf{v}}}_{* p}$ in the stress tensor ${{\rm{\Pi }}}_{r\theta }^{E}$, i.e. ${v}_{{Er}}{v}_{E\theta }$ becomes $({v}_{{Er}}+{v}_{* {pr}}){v}_{E\theta }$. This does not change the structure of the calculation as long as pressure fluctuations are proportional to potential fluctuations. Implications are discussed in section 4.2. A subsidiary small parameter ${\rho }_{\pi }=\tfrac{{\rho }_{i}}{{\lambda }_{\pi }}\lt 1$ is introduced at this point, which is the ratio of the thermal Larmor radius to the typical radial wavelength of the turbulent Reynolds stress tensor. Only terms of order ${\rho }_{\pi }^{2}$ are retained in the following.

The divergence of the parallel flux in the vorticity equation (33) can be calculated by computing the axisymmetric components $n=0,{\ell }\ne 0$ of the perturbed distribution function, which is a linear solution of the gyrokinetic equation (31), i.e.

$$\begin{eqnarray}&&{G}_{{\ell }\omega }={F}_{M}\displaystyle \frac{\omega }{\omega -{\ell }\tfrac{{v}_{\parallel }}{{q}_{0}{R}_{0}}+{\rm{i}}{0}^{+}}\displaystyle \frac{e{\phi }_{{\ell }\omega }}{T}.\end{eqnarray} \tag{ 37 }$$

It appears that the ${\ell }$ component of the parallel flux divergence is

$$\begin{eqnarray}&&{\left\{{{\rm{\nabla }}}_{\parallel }\left[\displaystyle \int {{\rm{d}}}^{3}{\bf{v}}{v}_{\parallel }F\right]\right\}}_{{\ell }}=-{\rm{i}}\omega ({\sigma }_{{\ell }\omega }+{\rm{i}}{\nu }_{{\ell }\omega })\displaystyle \frac{e{\phi }_{{\ell }\omega }}{T},\end{eqnarray} \tag{ 38 }$$

where ${\sigma }_{{\ell }\omega }\omega$ is a frequency shift and ${\nu }_{{\ell }\omega }\omega$ is the flow damping rate (${\nu }_{{\ell }\omega }\omega \gt 0$). The explicit expressions of ${\sigma }_{{\ell }\omega }$ and ${\nu }_{{\ell }\omega }$ are

$$\begin{eqnarray}&&{\sigma }_{{\ell }\omega }=P.P.\,\displaystyle \frac{1}{{(2\pi )}^{1/2}}{\int }_{-\infty }^{+\infty }{\rm{d}}\zeta {e}^{-\tfrac{{\zeta }^{2}}{2}}\displaystyle \frac{-{\ell }\tfrac{{v}_{{Ti}}}{{q}_{0}{R}_{0}}\zeta }{\omega -{\ell }\tfrac{{v}_{{Ti}}}{{q}_{0}{R}_{0}}\zeta +{\rm{i}}{0}^{+}}\end{eqnarray} \tag{ 39 }$$

and

$$\begin{eqnarray}&&{\nu }_{{\ell }\omega }={\left(\displaystyle \frac{\pi }{2}\right)}^{1/2}\displaystyle \frac{\omega }{| {\ell }| {\omega }_{t}}{e}^{-\tfrac{{\omega }^{2}}{2{{\ell }}^{2}{\omega }_{t}^{2}}},\end{eqnarray} \tag{ 40 }$$

where ${\omega }_{t}=\left|\tfrac{{v}_{{Ti}}}{{q}_{0}{R}_{0}}\right|$ is a thermal transit frequency. The last term in the vorticity equation (33) cancels out since the divergence of the diamagnetic current term is linear so that its toroidal average is zero. The other term is a particle flux which is equal to zero for an ITG turbulence. The vorticity equation (33) can then be solved by combining the Taylor identity equation (35) with the parallel flow divergence equation (38), thus providing the Fourier components of the electric potential ${\phi }_{{\ell }\omega }$

$$\begin{eqnarray}&&\displaystyle \frac{{\phi }_{{\ell }\omega }}{B}=-\displaystyle \frac{{K}_{\perp ,{\ell }}^{2}{\rho }_{i}^{2}}{1+{\sigma }_{{\ell }\omega }+{\rm{i}}{\nu }_{{\ell }\omega }}\displaystyle \frac{{{\rm{\Pi }}}_{r\theta ,{\ell }\omega }^{E}}{-{\rm{i}}\omega },\end{eqnarray} \tag{ 41 }$$

where the Reynolds stress tensor has been expanded in a Fourier series

$$\begin{eqnarray}&&{{\rm{\Pi }}}_{r\theta }^{E}(\theta ,t)=\displaystyle \sum _{{\ell }\omega }{{\rm{\Pi }}}_{r\theta ,{\ell }\omega }^{E}\exp \{{\rm{i}}({\ell }\theta -\omega t)\}\end{eqnarray} \tag{ 42 }$$

and ${K}_{\perp ,{\ell }}^{2}$ has to be understood as an operator that acts on ${{\rm{\Pi }}}_{r\theta ,{\ell }\omega }^{E}$, ${K}_{\perp ,{\ell }}^{2}{{\rm{\Pi }}}_{r\theta ,{\ell }\omega }^{E}=-{{\rm{\nabla }}}_{\perp }^{2}{{\rm{\Pi }}}_{r\theta ,{\ell }\omega }^{E}$. It is reminded that only the leading term in ${\rho }_{\pi }^{2}\simeq {K}_{\perp ,{\ell }}^{2}{\rho }_{i}^{2}$ is kept.

The expression of the electric potential Fourier components equation (41) is of central importance since it provides the structure of time-dependent poloidal convective cells which are driven by turbulent eddies. This result calls for several comments:

• 1.  
  It appears that ${\phi }_{{\ell }\omega }^{* }={\phi }_{-{\ell }-\omega }$, as expected. However, even if the stress tensor is up–down symmetrical at lowest order in ballooning angle, i.e. ${{\rm{\Pi }}}_{r\theta ,-{\ell }\omega }^{E}={{\rm{\Pi }}}_{r\theta ,{\ell }\omega }^{E}$ (a consequence of the parity of ${C}_{r,{\bf{k}}{\ell }}$ and ${C}_{\parallel ,{\bf{k}}{\ell }}$ with ${\ell }$), this is not the case for the potential, which develops up–down asymmetries ${\phi }_{-{\ell }\omega }\ne {\phi }_{{\ell }\omega }$, due to damping.
• 2.  
  Poloidal convective cells are Landau damped due to their finite poloidal wave number, as expected. However it appears that their damping rate is small at low frequencies $\omega \ll {\omega }_{t}$, hence favoring their onset and sustainment. It appears that the $1/\omega$ scaling of the frequency spectrum plays an important role in the generation of parallel flow via the magnetic drift. Obviously the limit $\omega =0$ is an artifact. In fact two cut-off frequencies appear, that correspond to the approximations made in this calculation: the trapped ion bounce pulsation, and the ion collision frequency. The highest is the ion bounce pulsation, which remains nevertheless lower than the ion transit frequency in the limit of low aspect ratio. Collisions regularize the $1/\omega$ infrared singularity via resonance broadening.
• 3.  
  The value of ${K}_{\perp ,{\ell }}{\rho }_{i}\simeq {\rho }_{\pi }$ plays also an important role. First, it is assumed that ${K}_{\perp ,{\ell }}^{2}{\rho }_{i}^{2}$ is positive, i.e. we consider here the case of radial oscillations of the Reynolds stress. Secondly, the mechanism presented here requires values of ${\rho }_{\pi }$ that are not too small, i.e. corrugations of the Reynolds stress with a spatial scale of a few gyroradii [11, 13, 16]. This may push the present model to its validity limit, since it was assumed that the variation with r₀ correspond to scales larger than the distance between resonant surfaces, itself of the order of the ion gyroradius. A more accurate calculation for ${\rho }_{\pi }\simeq 1$ is intricate as it requires a more precise procedure to solve the vorticity equation, as done for zonal flows ([40] and references therein, and also for convective cells driven by turbulence at low magnetic shear [23]).
• 4.  
  The ordering of gyrokinetic theory is not broken since the amplitude of the normalized potential $\tfrac{e{\phi }_{{\ell }\omega }}{{T}_{i}}$ remains small, of the order of ${\rho }_{* }$. This can be verified by using equation (41) , and noting that the turbulence stress tensor scales as ${\rho }_{* }^{2}{v}_{T}^{2}$, frequencies as $\tfrac{{v}_{T}}{R}$ so that $\tfrac{e{\phi }_{{\ell }\omega }}{{T}_{i}}\simeq {\rho }_{\pi }^{2}{\rho }_{* }$. Macroscopic flows ${\rho }_{\pi }\simeq {\rho }_{* }$ would lead to very small amplitudes of the potential turbulence on the electric potential $\tfrac{e{\phi }_{{\ell }\omega }}{{T}_{i}}\simeq {\rho }_{* }^{3}$ with a negligible effect, as discussed in [12]. Numerical simulations suggest that ${\rho }_{\pi }$ is rather in the range of 0.1 [13].

3.2. Magnetic drift contribution to the flux of parallel momentum

The banana-plateau component of the neoclassical flux of parallel momentum is known to exhibit a Pfirsch-Schlüter scaling and thus vanishes in the collisionless regime [4, 5, 41] (see [42] for an overview). Therefore one is left with the cross-correlation between the magnetic drift and the perturbed distribution function due to the electric potential fluctuations ${\phi }_{{\ell }\omega }$ driven by turbulent eddies. This cross-correlation is responsible for a finite flux of parallel momentum. The time dependent radial flux of parallel momentum due to magnetic drift reads

$$\begin{eqnarray}{\langle {{\rm{\Pi }}}_{r\parallel }^{D}\rangle }_{\theta } & = & \displaystyle \frac{1}{N}{\displaystyle \int }_{0}^{2\pi }\displaystyle \frac{{\rm{d}}\varphi }{2\pi }{\displaystyle \int }_{0}^{2\pi }\displaystyle \frac{{\rm{d}}\theta }{2\pi }\displaystyle \int {{\rm{d}}}^{3}{\bf{v}}{{Fv}}_{{Dr}}{v}_{\parallel }\\ & = & \displaystyle \sum _{{\ell }=-\infty }^{+\infty }\displaystyle \int {{\rm{d}}}^{3}{\bf{v}}{v}_{\parallel }{F}_{{\ell }}{v}_{{Dr}{\ell }}^{* },\end{eqnarray} \tag{ 43 }$$

where ${v}_{{Dr}}$ is the radial component of the magnetic drift velocity. This flux can be calculated using the expression equation (41) of the amplitudes of poloidal convective cells (see appendix D). A compact expression is found for its time average

$$\begin{eqnarray}&&{\langle {{\rm{\Pi }}}_{r\parallel }^{D}\rangle }_{\theta ,t}=-{q}_{0}{\langle \cos (\theta ){\rho }_{i}^{2}{{\rm{\nabla }}}_{\perp }^{2}{{\rm{\Pi }}}_{r\theta }^{E}\rangle }_{\theta ,t}.\end{eqnarray} \tag{ 44 }$$

For a ballooned turbulence, the Reynolds stress is ballooned too. If fluctuations are strongly localized on the low field side, the prefactor $\cos (\theta )$ can be replaced by 1 in equation (44). This can also be demonstrated by using the structure of the stress tensor equation (21) and the expression of ${C}_{r,{\bf{k}}{\ell }}$ with ${\ell }$ (see equation (26)).

Given the simplicity of equation (44), one may actually wonder whether a more direct calculation of ${{\rm{\Pi }}}_{r\parallel }^{D}$ is possible, and indeed it is. Since the parallel velocity plays a subdominant role in the magnetic drift, essentially because the resonant velocity goes like the pulsation, the curvature driven flux of momentum equation (43) can be reformulated as

$$\begin{eqnarray}&&{\langle {{\rm{\Pi }}}_{r\parallel }^{D}\rangle }_{\theta ,t}=-\displaystyle \frac{T}{{{NeB}}_{0}{R}_{0}}{\displaystyle \int }_{0}^{2\pi }\displaystyle \frac{{\rm{d}}\varphi }{2\pi }{\displaystyle \int }_{0}^{2\pi }\displaystyle \frac{{\rm{d}}\theta }{2\pi }\sin (\theta ){{\rm{\Gamma }}}_{\parallel }.\end{eqnarray} \tag{ 45 }$$

An integration by part allows a reformulation in terms of the parallel gradient of the parallel flux

$$\begin{eqnarray}&&{\langle {{\rm{\Pi }}}_{r\parallel }^{D}\rangle }_{\theta ,t}=-{q}_{0}\displaystyle \frac{T}{{{NeB}}_{0}}{\displaystyle \int }_{0}^{2\pi }\displaystyle \frac{{\rm{d}}\varphi }{2\pi }{\displaystyle \int }_{0}^{2\pi }\displaystyle \frac{{\rm{d}}\theta }{2\pi }\cos (\theta ){{\rm{\nabla }}}_{\parallel }{{\rm{\Gamma }}}_{\parallel }.\end{eqnarray} \tag{ 46 }$$

Using the vorticity equation (33), and the Taylor identity equation (35), one finds equation (44). It is quite remarkable that this result does not depend on the details of the poloidal convective cells, which act as mediators. An example of turbulence self-organization via the generation of poloidal convective cells is shown in figure 3. This figure comes from simulations run with the GYSELA gyrokinetic code [43]. It turns out that poloidal convective cells play also a role in the interplay between turbulent and neoclassical impurity transport [17].

Zoom In Zoom Out Reset image size

4.1. General expressions of momentum fluxes fluxes

The expression equation (44) of the radial flux of parallel momentum due to the magnetic drift can be expressed in the Fourier space using equation (21)

$$\begin{eqnarray}&&{\langle {{\rm{\Pi }}}_{r\parallel }^{D}\rangle }_{\theta ,t}=-{q}_{0}\displaystyle \sum _{{\bf{k}}}{C}_{r{\bf{k}}}{K}_{\perp ,{\ell }}^{2}{\rho }_{i}^{2}{\langle {| {v}_{E{\bf{k}}}| }^{2}\rangle }_{t}{\theta }_{k}.\end{eqnarray} \tag{ 47 }$$

This expression can be compared with the $E\times B$ flux of momentum equation (77)

$$\begin{eqnarray}&&{\langle {{\rm{\Pi }}}_{r\parallel }^{E}\rangle }_{\theta ,t}=\displaystyle \sum _{{\bf{k}}}{\omega }_{c}{\tau }_{{\bf{k}}}{C}_{\parallel {\bf{k}}}{\langle {| {v}_{E{\bf{k}}}| }^{2}\rangle }_{t}{\theta }_{k}.\end{eqnarray} \tag{ 48 }$$

It appears that these two contributions to the total flux of parallel momentum are anti-correlated when ${C}_{r{\bf{k}}}$ and ${C}_{\parallel {\bf{k}}}$ are of the same sign (for positive safety factor ${q}_{0}$). This is also the condition for anti-correlation of the $E\times B$ poloidal and parallel components of the Reynolds stress. The final answer is clearly sensitive to processes that determine the radial and parallel wavenumbers, which are presumably determined by non linear effects. Nevertheless it is interesting to calculate the linear values of ${C}_{r{\bf{k}}}$ and ${C}_{\parallel {\bf{k}}}$, and the corresponding signs of fluxes.

4.2. Application to ITG driven turbulence

An instructive example is the case of an unstable linear toroidal ITG driven mode. General expressions of ${C}_{r{\bf{k}}}$ and ${C}_{\parallel {\bf{k}}}$ are derived in the hydrodynamic limit in the appendix E. In the long wavelength limit ${K}_{\theta }^{2}{\rho }_{i}^{2}\ll 1/{s}_{0}^{2}$ and strong magnetic shear $| {s}_{0}| \gg 1$, explicit expressions of the coefficients ${C}_{r{\bf{k}}}$ and ${C}_{\parallel {\bf{k}}}$ can be found

$$\begin{eqnarray}&&{C}_{r{\bf{k}}}={c}_{r{\bf{k}}}{s}_{0},\end{eqnarray} \tag{ 49 }$$

$$\begin{eqnarray}&&{C}_{\parallel {\bf{k}}}={c}_{\parallel {\bf{k}}}\displaystyle \frac{{\rho }_{i}}{{R}_{0}}{({K}_{\theta }{\rho }_{i}{s}_{0})}^{2}{q}_{0}{s}_{0}.\end{eqnarray} \tag{ 50 }$$

The numbers ${c}_{r{\bf{k}}}$ and ${c}_{\parallel {\bf{k}}}$ depend on the normalized frequencies $\tfrac{{\omega }_{{\bf{k}}}}{{\omega }_{d}}$ and $\tfrac{{\gamma }_{{\bf{k}}}}{{\omega }_{d}}$ only, where ${\omega }_{{\bf{k}}}$ and ${\gamma }_{{\bf{k}}}$ are the real and imaginary parts of the linear frequency in the plasma frame (frame where the local radial electric field is zero), and

$$\begin{eqnarray}&&{\omega }_{d}=-2\displaystyle \frac{{K}_{\theta }T}{{{eB}}_{0}{R}_{0}}\end{eqnarray} \tag{ 51 }$$

is a thermal magnetic field curvature frequency. More precisely

$$\begin{eqnarray}&&{c}_{r{\bf{k}}}=\displaystyle \frac{1}{2}\displaystyle \frac{1}{{\left(1+\tfrac{{\gamma }_{{\bf{k}}}^{2}}{{\omega }_{{\bf{k}}}^{2}}\right)}^{1/2}}({\lambda }_{-}-\mathrm{sgn}({\omega }_{{\bf{k}}}{\omega }_{d}{s}_{0}))\end{eqnarray} \tag{ 52 }$$

$$\begin{eqnarray}&&{c}_{\parallel {\bf{k}}}={2}^{1/2}\displaystyle \frac{1}{| {q}_{0}| }{\left|\displaystyle \frac{{\omega }_{{\bf{k}}}}{{\omega }_{d}{s}_{0}}\right|}^{3/2}{\left(1+\displaystyle \frac{{\gamma }_{{\bf{k}}}^{2}}{{\omega }_{{\bf{k}}}^{2}}\right)}^{1/2}{\lambda }_{-}^{1/2}\end{eqnarray} \tag{ 53 }$$

and

$$\begin{eqnarray}&&{\lambda }_{-}={\left(1+\displaystyle \frac{{\gamma }_{{\bf{k}}}^{2}}{{\omega }_{{\bf{k}}}^{2}}\right)}^{1/2}-\mathrm{sgn}({\omega }_{{\bf{k}}}{\omega }_{d}{s}_{0}).\end{eqnarray} \tag{ 54 }$$

The coefficient ${\lambda }_{-}$ is always positive. These expressions are valid for any sign of ${q}_{0}$ and ${s}_{0}$. We now restrict the discussion to the generic case ${q}_{0}\gt 0$ and ${s}_{0}\gt 0$. These considerations lead to the following estimates of the momentum fluxes

$$\begin{eqnarray}&&{\langle {{\rm{\Pi }}}_{r\parallel }^{D}\rangle }_{\theta ,t}\simeq {s}_{0}{q}_{0}\displaystyle \sum _{{\bf{k}}}{c}_{r{\bf{k}}}{({K}_{\perp ,1}{\rho }_{i})}^{2}{\langle {| {v}_{E{\bf{k}}}| }^{2}\rangle }_{t}{\theta }_{k}\end{eqnarray} \tag{ 55 }$$

and

$$\begin{eqnarray}&&{\langle {{\rm{\Pi }}}_{r\parallel }^{E}\rangle }_{\theta ,t}\simeq -{s}_{0}{q}_{0}\displaystyle \sum _{{\bf{k}}}{c}_{\parallel {\bf{k}}}{\omega }_{t}{\tau }_{{\bf{k}}}{({K}_{\theta }{\rho }_{i}{s}_{0})}^{2}{\langle {| {v}_{E{\bf{k}}}| }^{2}\rangle }_{t}{\theta }_{k}.\end{eqnarray} \tag{ 56 }$$

For ${\omega }_{{\bf{k}}}\simeq {\gamma }_{{\bf{k}}}\simeq {\omega }_{d}$, the two numbers ${c}_{r{\bf{k}}}$ and ${c}_{\parallel {\bf{k}}}$ are of the same order of magnitude. The two fluxes are therefore comparable when ${K}_{\perp ,1}{\rho }_{i}\simeq {\rho }_{\pi }\simeq {K}_{\theta }{\rho }_{i}{({\omega }_{t}{\tau }_{{\bf{k}}})}^{1/2}$. This condition appears as reasonable since the intensity wave number spectrum typically peaks at ${K}_{\theta }{\rho }_{i}\simeq 0.1$, provided ${\omega }_{t}{\tau }_{{\bf{k}}}\simeq o(1)$ [44, 45].

Regarding the signs, it appears that for modes that drift in the ion diamagnetic direction $\mathrm{sgn}({\omega }_{{\bf{k}}}{\omega }_{d})=1$, ${c}_{r{\bf{k}}}$ is negative for a weak drive ${\gamma }_{{\bf{k}}}\ll {\omega }_{{\bf{k}}}$, and positive for strong drive ${\gamma }_{{\bf{k}}}\geqslant {\omega }_{{\bf{k}}}$. Since ${c}_{\parallel {\bf{k}}}$ is always positive, this means that ${C}_{\parallel {\bf{k}}}$ and ${C}_{r{\bf{k}}}$ are positively correlated near threshold, and anti-correlated far from threshold. Consequently the two fluxes or parallel momentum are positively correlated for low drive ${\gamma }_{{\bf{k}}}\ll {\omega }_{{\bf{k}}}$ and anti-correlated for strong drive ${\gamma }_{{\bf{k}}}\geqslant {\omega }_{{\bf{k}}}$. For modes drifting in the electron diamagnetic direction, $\mathrm{sgn}({\omega }_{{\bf{k}}}{\omega }_{d})=-1$, anti-correlation always occurs. For ITG modes, drift in the ion diamagnetic direction is expected, so that anti-correlation is expected only far enough from the instability threshold. It is not clear whether this finding agrees or not with previous numerical findings for momentum transport [13, 16]. Nevertheless it is stressed that the hydrodynamic limit that is being used here, is a rather poor approximation that becomes correct only well above the instability threshold, i.e. for ${\gamma }_{{\bf{k}}}\gg {\omega }_{{\bf{k}}}$. This is precisely the regime where an agreement with simulations is found, i.e. anticorrelation. In that regard some further analysis of the numerical simulations would be helpful. An encouraging observation though [16] is that the $E\times B$ flux of parallel momentum is anti-correlated with the ballooning angle ${\theta }_{k}$, which suggests positive ${c}_{\parallel {\bf{k}}}$ in view of equation (56). Let us note that adding the FLR effects mentioned in section 3.1 may change the correlation sign. Indeed changing the stress tensor ${v}_{{Er}}{v}_{E\theta }$ into $({v}_{{Er}}+{v}_{* {pr}}){v}_{E\theta }$ is roughly equivalent to multiplying ${c}_{r{\bf{k}}}$ in equation (52) by a factor $\left(1-\tfrac{{\omega }_{* p}}{{\omega }_{{\bf{k}}}}\right)$, where ${\omega }_{* p}$ is the pressure diamagnetic frequency. Since the later tends to be larger than the mode pulsation ${\omega }_{{\bf{k}}}$ for realistic parameters, this may change the sign of the coefficient ${c}_{r{\bf{k}}}$, and therefore the relative signs of the E × B and magnetic drift components of the momentum flux. As a final note, it is possible (if not likely) that the two terms ${c}_{r{\bf{k}}}$ and ${c}_{\parallel {\bf{k}}}$ are determined by non linear processes and not well captured by a linear analysis.

It is shown here that the $E\times B$ Reynolds stress tensor generates poloidal asymmetries of the plasma flow due to turbulence ballooning. These poloidal convective cells are weakly damped at low frequency. Their radial scale is dictated by the turbulent Reynolds stress, and their poloidal wave numbers are small. These cells drive up–down asymmetries of the distribution function, which are responsible for a non-zero radial flux of parallel momentum due to the geodesic component of the particle magnetic drift. The entire process requires a symmetry breaking mechanism, for instance a mean shear flow. Since the turbulent Reynolds can be seen as a flux of momentum, it appears that the two components of the radial flux of parallel momentum due to curvature and $E\times B$ drift are correlated. This general result comes from a simple relationship between the momentum flux due to magnetic drift and the turbulent Reynolds stress. Although poloidal convective cells do not appear explicitly in this relationship, they play an essential role in this mechanism. An analytic calculation shows that anti-correlation between the components of the turbulent Reynolds stress results in anti-correlation of the two contributions to the flux of parallel momentum that come from $E\times B$ and magnetic drifts. A quasi-linear calculation of all quantities, based on ITG linear stability indicates that positive correlation is expected near threshold, and anti-correlation for strong drive. Hence no firm conclusion can be drawn as to the relevance of this mechanism to explain the numerical results. Nevertheless the hydrodynamic limit that has been used is a rather poor representation of the turbulence near threshold. In fact it is valid only well above the stability threshold, i.e. strong drive. This is encouraging since it is the case where anti-correlation is found, in agreement with numerical findings. It is likely that poloidal convective cells generated by turbulence have other consequences on turbulent transport and turbulence. Indeed they may participate in turbulence self-regulation via vortex shearing processes similar to the well documented effect of zonal flows on turbulence.

The mode structure equation (9) lead to a turbulence intensity that reads

$$\begin{eqnarray}&&{| {\tilde{\phi }}_{{\bf{k}}}(\theta ,t)| }^{2}={| {\phi }_{{\bf{k}}}(t)| }^{2}{{ \mathcal A }}_{{\bf{k}}}(\theta ),\end{eqnarray} \tag{ 57 }$$

where ${{ \mathcal A }}_{{\bf{k}}}(\theta )$ is a form factor

$$\begin{eqnarray}&&{{ \mathcal A }}_{{\bf{k}}}(\theta )={\left(\displaystyle \frac{{\bar{\alpha }}_{{\bf{k}}}}{\pi }\right)}^{1/2}\exp \{-{\bar{\alpha }}_{{\bf{k}}}{(\theta -{\bar{\lambda }}_{{\bf{k}}}{\theta }_{k})}^{2}\}\end{eqnarray} \tag{ 58 }$$

and

$$\begin{eqnarray}{\bar{\alpha }}_{{\bf{k}}} & = & \displaystyle \frac{1}{2}({\alpha }_{{\bf{k}}}+{\alpha }_{{\bf{k}}}^{* })={\mathfrak{R}}({\alpha }_{{\bf{k}}})\\ {\bar{\lambda }}_{{\bf{k}}} & = & \displaystyle \frac{{\alpha }_{{\bf{k}}}{\lambda }_{{\bf{k}}}+{\alpha }_{{\bf{k}}}^{* }{\lambda }_{{\bf{k}}}^{* }}{{\alpha }_{{\bf{k}}}+{\alpha }_{{\bf{k}}}^{* }}=\displaystyle \frac{{\mathfrak{R}}({\alpha }_{{\bf{k}}}{\lambda }_{{\bf{k}}})}{{\mathfrak{R}}({\alpha }_{{\bf{k}}})}\end{eqnarray} \tag{ 59 }$$

The amplitude ${a}_{{\bf{k}}}$ in equation (9) has been chosen such that

$$\begin{eqnarray}&&{\int }_{-\infty }^{\infty }{\rm{d}}\theta {| {\tilde{\phi }}_{{\bf{k}}}(\theta ,t)| }^{2}={| {\tilde{\phi }}_{{\bf{k}}}(t)| }^{2}\end{eqnarray} \tag{ 60 }$$

namely

$$\begin{eqnarray}&&{a}_{{\bf{k}}}={\left(\displaystyle \frac{{\bar{\alpha }}_{{\bf{k}}}}{\pi }\right)}^{1/4}\exp \left\{\displaystyle \frac{{| {\alpha }_{{\bf{k}}}| }^{2}}{2({\alpha }_{{\bf{k}}}+{\alpha }_{{\bf{k}}}^{* })}{({\lambda }_{{\bf{k}}}-{\lambda }_{{\bf{k}}}^{* })}^{2}\right\}.\end{eqnarray} \tag{ 61 }$$

The poloidal Fourier components of the radial and parallel wave numbers equation (25) and equation (26) can then be recast as

$$\begin{eqnarray}{K}_{r,{\bf{k}}{\ell }} & = & {K}_{\theta }{s}_{0}{\int }_{-\infty }^{\infty }{\rm{d}}\theta {{\rm{e}}}^{-{\rm{i}}{\ell }\theta }(\theta -{\theta }_{k}){{ \mathcal A }}_{{\bf{k}}}(\theta ),\\ {K}_{\parallel ,{\bf{k}}{\ell }} & = & \displaystyle \frac{1}{{q}_{0}{R}_{0}}{\int }_{-\infty }^{\infty }{\rm{d}}\theta {{\rm{e}}}^{-{\rm{i}}{\ell }\theta }(\theta -{\lambda }_{{\bf{k}}}{\theta }_{k}){{ \mathcal A }}_{{\bf{k}}}(\theta ).\end{eqnarray} \tag{ 62 }$$

For a strongly ballooned turbulence ${\alpha }_{{\bf{k}}}\gg 1$, and small values of the ballooning angle ${\theta }_{k}\ll 1$ , one gets at first order in ${\theta }_{k}$

$$\begin{eqnarray}&&{K}_{r,{\bf{k}}{\ell }}={K}_{r,0}={K}_{\theta }{s}_{0}\displaystyle \frac{{\mathfrak{R}}({\alpha }_{{\bf{k}}}{\delta }_{{\bf{k}}})}{{\mathfrak{R}}({\alpha }_{{\bf{k}}})}{\theta }_{k},\end{eqnarray} \tag{ 63 }$$

$$\begin{eqnarray}&&{K}_{\parallel ,{\bf{k}}{\ell }}={K}_{\parallel ,0}=\displaystyle \frac{{| {\alpha }_{{\bf{k}}}| }^{2}}{{q}_{0}{R}_{0}}\displaystyle \frac{{\mathfrak{I}}({\delta }_{{\bf{k}}})}{{\mathfrak{R}}({\alpha }_{{\bf{k}}})}{\theta }_{k},\end{eqnarray} \tag{ 64 }$$

where ${\delta }_{{\bf{k}}}={\lambda }_{{\bf{k}}}-1$. Hee ${\mathfrak{R}}(z)$ and ${\mathfrak{I}}(z)$ indicate the real and imaginary parts of a complex number z.

We start from the expression of the E × B drift contribution to the momentum flux, and readily get its time average equation (12)

$$\begin{eqnarray}&&{\langle {{\rm{\Pi }}}_{r\parallel }^{E}(\theta ,t)\rangle }_{t}=\displaystyle \frac{1}{N}\displaystyle \sum _{{\bf{k}}\omega }\displaystyle \int {{\rm{d}}}^{3}{\bf{v}}{\tilde{F}}_{{\bf{k}}\omega }{v}_{\parallel }{\tilde{v}}_{E{\bf{k}}\omega }^{* }.\end{eqnarray} \tag{ 65 }$$

To calculate the distribution function versus the potential, we use a ballooning representation. The electric potential is written in the form

$$\begin{eqnarray}&&\phi (r,\theta ,\varphi ,t)=\displaystyle \sum _{p=-\infty }^{+\infty }{\tilde{\phi }}_{{\bf{k}}\omega }(\theta +2p\pi ,t)\exp \{{\rm{i}}{n}[\varphi -q(r)(\theta +2p\pi -{\theta }_{k})]-{\rm{i}}\omega t\},\end{eqnarray} \tag{ 66 }$$

where ${\theta }_{k}$ is the ballooning mode. The single term p = 0 is kept for strongly ballooned fluctuations. A similar expansion is used for the non adiabatic part of the distribution function, i.e. ${\tilde{G}}_{{\bf{k}}\omega }={\tilde{F}}_{{\bf{k}}\omega }+\tfrac{e}{{T}_{i}}{\tilde{\phi }}_{{\bf{k}}\omega }$. The gyrokinetic equation (31) reads

$$\begin{eqnarray}&&(\omega -{K}_{\parallel }{v}_{\parallel }-{\omega }_{D}){\tilde{G}}_{{\bf{k}}\omega }={F}_{M}(\omega -{\omega }_{* })\displaystyle \frac{e}{{T}_{i}}({ \mathcal J }\cdot {\tilde{\phi }}_{{\bf{k}}\omega }).\end{eqnarray} \tag{ 67 }$$

Here ${K}_{\parallel }$ is an operator

$$\begin{eqnarray}&&{K}_{\parallel }=-{\rm{i}}\displaystyle \frac{1}{{q}_{0}{R}_{0}}\displaystyle \frac{\partial }{\partial \eta },\end{eqnarray} \tag{ 68 }$$

where $\eta =\theta -{\theta }_{k}$ is a shifted poloidal angle. The kinetic magnetic drift frequency ${\omega }_{D}$ is defined as

$$\begin{eqnarray}&&{\omega }_{D}(\eta )=\left(\displaystyle \frac{{v}_{\perp }^{2}}{2{v}_{{Ti}}^{2}}+\displaystyle \frac{{v}_{\parallel }^{2}}{{v}_{{Ti}}^{2}}\right){\omega }_{d}(\eta ),\end{eqnarray} \tag{ 69 }$$

where

$$\begin{eqnarray}&&{\omega }_{d}(\eta )=-\displaystyle \frac{2{K}_{\theta }{T}_{i}}{{{eB}}_{0}{R}_{0}}[\cos ({\theta }_{k}+\eta )+{s}_{0}\eta \sin ({\theta }_{k}+\eta )].\end{eqnarray} \tag{ 70 }$$

${v}_{{Ti}}=\sqrt{{T}_{i}/{m}_{i}}$ is the thermal ion velocity, and ${K}_{\theta }=\tfrac{-{{nq}}_{0}}{{r}_{0}}$ is the poloidal wavenumber—${q}_{0}$ is the safety factor at the reference radius ${r}_{0}$. The definition of the kinetic diamagnetic frequency ${\omega }_{* }$ is the usual one

$$\begin{eqnarray}&&{\omega }_{* }={\omega }_{* n}+{\omega }_{* T}\left(\displaystyle \frac{{v}_{\perp }^{2}}{{v}_{{Ti}}^{2}}+\displaystyle \frac{{v}_{\parallel }^{2}}{{v}_{{Ti}}^{2}}-\displaystyle \frac{3}{2}\right),\end{eqnarray} \tag{ 71 }$$

where

$$\begin{eqnarray}{\omega }_{* n} & = & -\displaystyle \frac{{K}_{\theta }{T}_{i}}{{{eB}}_{0}{L}_{{n}_{i}}}\\ {\omega }_{* T} & = & -\displaystyle \frac{{K}_{\theta }{T}_{i}}{{{eB}}_{0}{L}_{{T}_{i}}}\end{eqnarray} \tag{ 72 }$$

are the density and temperature diamagnetic frequencies and ${L}_{{n}_{i}}$ and ${L}_{{p}_{i}}$ are the density and pressure gradient lengths calculated at the reference radius ${r}_{0}$. The gyroaverage ${ \mathcal J }$ is fairly well represented by a Bessel function with an argument ${K}_{\perp }{\rho }_{c}$, i.e. ${J}_{0}({K}_{\perp }{\rho }_{c})$, where ${\rho }_{c}=\tfrac{{m}_{i}{v}_{\perp }}{{{eB}}_{0}}$ is the kinetic ion gyroradius and

$$\begin{eqnarray}&&{K}_{\perp }^{2}={K}_{\theta }^{2}[1+{s}_{0}^{2}{\eta }^{2}].\end{eqnarray} \tag{ 73 }$$

The solution of the gyrokinetic equation (67) involves an integro-differential operator that relates ${\tilde{G}}_{{\bf{k}}\omega }$ to ${\tilde{\phi }}_{{\bf{k}}\omega }$ [46, 47]. A formal solution can be written in a Wentzel-Kramers-Brilloin (WKB) sense by dividing the rhs of equation (67) by the resonant term $\omega -{K}_{\parallel }{v}_{\parallel }-{\omega }_{D}$

$$\begin{eqnarray}&&{\tilde{F}}_{{\bf{k}}\omega }=-{F}_{M}\left\{1-\displaystyle \frac{\omega -{\omega }_{* }}{\omega -{\omega }_{D}-{K}_{\parallel }{v}_{\parallel }+{\rm{i}}{0}^{+}}\right\}\displaystyle \frac{e{\tilde{\phi }}_{{\bf{k}}\omega }}{T}.\end{eqnarray} \tag{ 74 }$$

We assume that the spectral turbulence intensity is of the form

$$\begin{eqnarray}&&{| {\tilde{\phi }}_{{\bf{k}}\omega }| }^{2}={{ \mathcal A }}_{{\bf{k}}}(\theta ){| {\phi }_{{\bf{k}}}| }^{2}\displaystyle \frac{1}{\pi }\displaystyle \frac{{\rm{\Delta }}{\omega }_{{\bf{k}}}}{{(\omega -{\omega }_{{\bf{k}}})}^{2}+{\rm{\Delta }}{\omega }_{{\bf{k}}}^{2}}\end{eqnarray} \tag{ 75 }$$

which gives

$$\begin{eqnarray}&&{\langle {{\rm{\Pi }}}_{r\parallel }^{E}(\theta ,t)\rangle }_{t}=\displaystyle \frac{1}{N}\displaystyle \sum _{{\bf{k}}\omega }\displaystyle \int {{\rm{d}}}^{3}{\bf{v}}{F}_{M}{v}_{\parallel }{| {v}_{E{\bf{k}}}| }^{2}\displaystyle \frac{{{eB}}_{0}}{{K}_{\theta }{T}_{i}}(\omega -{\omega }_{* })\displaystyle \frac{{\rm{\Delta }}{\omega }_{{\bf{k}}}}{{(\omega -{\omega }_{{\bf{k}}})}^{2}+{\rm{\Delta }}{\omega }_{{\bf{k}}}^{2}}\delta (\omega -{K}_{\parallel }{v}_{\parallel }-{\omega }_{D}).\end{eqnarray} \tag{ 76 }$$

In the hydrodynamic limit $\tfrac{{\omega }_{D}}{\omega }\simeq \tfrac{{K}_{\parallel }^{2}{v}_{\parallel }^{2}}{{\omega }^{2}}\ll 1$, one gets the result

$$\begin{eqnarray}&&{\langle {{\rm{\Pi }}}_{r\parallel }^{E}\rangle }_{\theta ,t}=\displaystyle \sum _{{\bf{k}}}{C}_{\parallel {\bf{k}}}\displaystyle \frac{{\omega }_{c}{\rm{\Delta }}{\omega }_{{\bf{k}}}}{{\omega }_{{\bf{k}}}^{2}+{\rm{\Delta }}{\omega }_{{\bf{k}}}^{2}}{\langle {| {v}_{E{\bf{k}}}| }^{2}\rangle }_{t}{\theta }_{k}.\end{eqnarray} \tag{ 77 }$$

This result is identical to equation (22) with ${\tau }_{{\bf{k}}}=\tfrac{{\rm{\Delta }}{\omega }_{{\bf{k}}}}{{\omega }_{{\bf{k}}}^{2}+{\rm{\Delta }}{\omega }_{{\bf{k}}}^{2}}$. Strictly speaking there is also a contribution to the $E\times B$ flux of parallel momentum that comes from the axisymmetric perturbations of the potential ${\phi }_{{\ell }\omega }$ equation (41) and the distribution function response equation (37) . However this contribution is of second order in ballooning angle ${\theta }_{k}$ and will not be retained here. Hence the $E\times B$ flux of momentum is in this peculiar case 'turbulent', i.e. produced by small scale fluctuations.

The purpose of this appendix is to demonstrate the identity

$$\begin{eqnarray}&&\displaystyle \frac{1}{{B}_{0}}{\int }_{0}^{2\pi }\displaystyle \frac{{\rm{d}}\varphi }{2\pi }({{\bf{v}}}_{E}\cdot {\rm{\nabla }}){\rm{\Omega }}={{\rm{\nabla }}}_{\perp }^{2}{{\rm{\Pi }}}_{r\theta }^{E},\end{eqnarray} \tag{ 78 }$$

where ${\rm{\Omega }}={{\rm{\nabla }}}_{\perp }^{2}\phi$ is the vorticity, and ${{\bf{v}}}_{E}={\bf{b}}\times \tfrac{{\rm{\nabla }}\phi }{B}$ the $E\times B$ drift velocity. The demonstration is restricted to a geometry of concentric circular surfaces with large aspect ratio . We use the mode structure equation (8)

$$\begin{eqnarray}&&\phi (r,\theta ,\varphi ,t)=\displaystyle \sum _{{\bf{k}}}{\tilde{\phi }}_{{\bf{k}}}(r,\theta ,t)\exp \{{\rm{i}}{n}{\chi }_{{\bf{k}}}\},\end{eqnarray} \tag{ 79 }$$

where ${\chi }_{{\bf{k}}}=\varphi -q(r)(\theta -{\theta }_{k})$, and the dependence on r of ${\tilde{\phi }}_{{\bf{k}}}(r,\theta ,t)$ is made explicit for clarity throughout this section. The average over ϕ implies that the lhs of equation (78) appears as a sum over the index ${\bf{k}}$ of operators acting on ${\tilde{\phi }}_{{\bf{k}}}$, namely

$$\begin{eqnarray}&&{\int }_{0}^{2\pi }\displaystyle \frac{{\rm{d}}\varphi }{2\pi }\displaystyle \frac{{\bf{b}}}{{B}_{0}^{2}}\cdot ({\rm{\nabla }}\phi \times {\rm{\nabla }}{\rm{\Omega }})=\displaystyle \sum _{{\bf{k}}}\displaystyle \frac{{\bf{b}}}{{B}_{0}^{2}}\cdot [{({\rm{\nabla }}\phi )}_{{\bf{k}}}\times {({\rm{\nabla }}{\rm{\Omega }})}_{{\bf{k}}}^{* }],\end{eqnarray} \tag{ 80 }$$

where

$$\begin{eqnarray}&&{({\rm{\nabla }}\phi )}_{{\bf{k}}}=\displaystyle \frac{\partial {\tilde{\phi }}_{{\bf{k}}}}{\partial r}{\rm{\nabla }}r+\displaystyle \frac{\partial {\tilde{\phi }}_{{\bf{k}}}}{\partial \theta }{\rm{\nabla }}\theta +{\rm{i}}{n}{\rm{\nabla }}{\chi }_{{\bf{k}}}{\tilde{\phi }}_{{\bf{k}}}\end{eqnarray} \tag{ 81 }$$

and ${({\rm{\nabla }}{\rm{\Omega }})}_{{\bf{k}}}$ is given by a similar expression. This suggests a change of variables for each ${\bf{k}}$, namely

$$\begin{eqnarray}&&{r}_{0}=r,\end{eqnarray} \tag{ 82 }$$

$$\begin{eqnarray}&&{\chi }_{{\bf{k}}}=\varphi -q(r)(\theta -{\theta }_{k}),\end{eqnarray} \tag{ 83 }$$

$$\begin{eqnarray}&&\eta =\theta .\end{eqnarray} \tag{ 84 }$$

The variable η plays the role of a coordinate along the field line, r₀ represents the slow variation of the field in the radial electric field, and ${\chi }_{{\bf{k}}}$ is a coordinate transverse to the field. We use the following ordering

$$\begin{eqnarray}&&\displaystyle \frac{1}{{R}_{0}}\displaystyle \frac{\partial {\tilde{\phi }}_{{\bf{k}}}}{\partial \eta }\ll \displaystyle \frac{\partial {\tilde{\phi }}_{{\bf{k}}}}{\partial {r}_{0}}\ll \displaystyle \frac{{{nq}}_{0}}{{r}_{0}}{\tilde{\phi }}_{{\bf{k}}}\end{eqnarray} \tag{ 85 }$$

which allows neglecting the slow variation of the field along the field line. We note that the unit vector along the magnetic field can be written,

$$\begin{eqnarray}&&{\bf{b}}=\displaystyle \frac{{\rm{\nabla }}{\chi }_{{\bf{k}}}}{| {\rm{\nabla }}{\chi }_{{\bf{k}}}| }\times {\rm{\nabla }}{r}_{0},\end{eqnarray} \tag{ 86 }$$

where

$$\begin{eqnarray}&&| {\rm{\nabla }}{\chi }_{{\bf{k}}}| =\displaystyle \frac{{q}_{0}}{{r}_{0}}{[1+{s}_{0}^{2}{(\eta -{\theta }_{k})}^{2}]}^{1/2}.\end{eqnarray} \tag{ 87 }$$

For a strongly ballooned turbulence $| {\rm{\nabla }}{\chi }_{{\bf{k}}}| \simeq \tfrac{{q}_{0}}{{r}_{0}}$ . Ignoring the derivatives with respect to η, one finds the following identity

$$\begin{eqnarray}&&{\int }_{0}^{2\pi }\displaystyle \frac{{\rm{d}}\varphi }{2\pi }\displaystyle \frac{{\bf{b}}}{{B}_{0}^{2}}\cdot ({\rm{\nabla }}\phi \times {\rm{\nabla }}{\rm{\Omega }})=\displaystyle \frac{1}{{B}_{0}^{2}}\displaystyle \frac{{q}_{0}}{{r}_{0}}\displaystyle \frac{\partial }{\partial {r}_{0}}\left\{\displaystyle \sum _{{\bf{k}}}({\rm{i}}{n}{\tilde{\phi }}_{{\bf{k}}}{{\rm{\Omega }}}_{{\bf{k}}}^{* })\right\}.\end{eqnarray} \tag{ 88 }$$

Using ${{\rm{\nabla }}}_{\perp }^{2}={\rm{\nabla }}\cdot {\rm{\nabla }}-{({\bf{b}}\cdot {\rm{\nabla }})}^{2}$, and neglecting again the derivatives along the field lines. Each Fourier harmonics of the vorticity can be written as

$$\begin{eqnarray}&&{{\rm{\Omega }}}_{{\bf{k}}}=\displaystyle \frac{1}{{r}_{0}}\displaystyle \frac{\partial }{\partial {r}_{0}}{r}_{0}\displaystyle \frac{\partial }{\partial {r}_{0}}({r}_{0}{\tilde{\phi }}_{{\bf{k}}})+2{\rm{i}}{n}({\rm{\nabla }}{\chi }_{{\bf{k}}}\cdot {\rm{\nabla }}{r}_{0})\displaystyle \frac{\partial }{\partial {r}_{0}}{\tilde{\phi }}_{{\bf{k}}}-{n}^{2}{| {\rm{\nabla }}{\chi }_{{\bf{k}}}| }^{2}{\tilde{\phi }}_{{\bf{k}}},\end{eqnarray} \tag{ 89 }$$

where it has been used that the variation of ${\tilde{\phi }}_{{\bf{k}}}$ with ${r}_{0}$ is faster than the radial variation of $({\rm{\nabla }}{\chi }_{{\bf{k}}}\cdot {\rm{\nabla }}{r}_{0})$. The largest term is the third one in the rhs of equation (89) , but it does not contribute to equation (88) for parity reasons. The largest contribution therefore comes from the second term of the rhs of equation (89) , i.e.

$$\begin{eqnarray}&&{\int }_{0}^{2\pi }\displaystyle \frac{{\rm{d}}\varphi }{2\pi }\displaystyle \frac{{\bf{b}}}{{B}_{0}^{2}}\cdot ({\rm{\nabla }}\phi \times {\rm{\nabla }}{\rm{\Omega }})\simeq {{\rm{\nabla }}}_{\perp }^{2}\left[-\displaystyle \sum _{{\bf{k}}}{s}_{0}(\eta -{\theta }_{k}){\left|\displaystyle \frac{{{nq}}_{0}}{{r}_{0}}\displaystyle \frac{{\tilde{\phi }}_{{\bf{k}}}}{{B}_{0}}\right|}^{2}\right],\end{eqnarray} \tag{ 90 }$$

where

$$\begin{eqnarray}&&{{\rm{\nabla }}}_{\perp }^{2}=\displaystyle \frac{1}{{r}_{0}}\displaystyle \frac{\partial }{\partial {r}_{0}}{r}_{0}\displaystyle \frac{\partial }{\partial {r}_{0}}.\end{eqnarray} \tag{ 91 }$$

The expression within the brackets of the rhs of equation (90) is the $r\theta$ component of the $E\times B$ stress tensor, as can be verified from equation (14). Equation (90) demonstrates the Taylor identity equation (78).

The purpose of this section is to compute the radial flux of parallel momentum due to magnetic drift equation (43). The geodesic component ${v}_{{Dr}}$ of the magnetic drift velocity reads

$$\begin{eqnarray}&&{v}_{{Dr}}={v}_{D}\displaystyle \frac{1}{2{\rm{i}}}({\delta }_{{\ell },1}-{\delta }_{{\ell },-1}),\end{eqnarray} \tag{ 92 }$$

where

$$\begin{eqnarray}&&{v}_{D}=-\displaystyle \frac{\mu {B}_{0}+{{mv}}_{\parallel }^{2}}{{{eB}}_{0}{R}_{0}}.\end{eqnarray} \tag{ 93 }$$

We compute the momentum flux by using again the linear solution equation (37) of the gyrokinetic equation (31). The frequency Fourier components of the parallel momentum flux equation (43) read

$$\begin{eqnarray}{\langle {{\rm{\Pi }}}_{r\parallel \omega }^{D}\rangle }_{\theta } & = & -{q}_{0}\displaystyle \sum _{{\ell }=\pm 1}{\displaystyle \int }_{-\infty }^{+\infty }\displaystyle \frac{{\rm{d}}\zeta }{{(2\pi )}^{1/2}}{e}^{-\tfrac{{\zeta }^{2}}{2}}\displaystyle \frac{-{\ell }\tfrac{{v}_{{Ti}}}{{q}_{0}{R}_{0}}\zeta }{\omega -{\ell }\tfrac{{v}_{{Ti}}}{{q}_{0}{R}_{0}}\zeta +{\rm{i}}{0}^{+}}\\ & & \times (1+{\zeta }^{2})\displaystyle \frac{1}{2}({\delta }_{{\ell },1}+{\delta }_{{\ell },-1})\displaystyle \frac{-{\rm{i}}\omega {\phi }_{{\ell }\omega }}{B}.\end{eqnarray} \tag{ 94 }$$

Using equation (38) (or equivalently equations (40), (39)), the integral over the velocity can be expressed as a function of ${\sigma }_{{\ell }\omega }$ and ${\nu }_{{\ell }\omega }$, namely

$$\begin{eqnarray}&&{\langle {{\rm{\Pi }}}_{r\parallel \omega }^{D}\rangle }_{\theta }=-{q}_{0}\displaystyle \sum _{{\ell }=\pm 1}\left[1+\left(1+\displaystyle \frac{{\omega }^{2}}{{\omega }_{t}^{2}}\right)({\sigma }_{{\ell }\omega }+{\rm{i}}{\nu }_{{\ell }\omega })\right]\displaystyle \frac{1}{2}({\delta }_{{\ell },1}+{\delta }_{{\ell },-1})\displaystyle \frac{-{\rm{i}}\omega {\phi }_{{\ell }\omega }}{B}.\end{eqnarray} \tag{ 95 }$$

We note that it is essential to deal with time dependent perturbations to get a finite flux. Static perturbations resonate at zero parallel velocity ${v}_{\parallel }=0$ and therefore do not contribute to the flux of parallel momentum. This is also the basic reason why the banana-plateau neoclassical flux of parallel momentum is zero. In the present case, this cancellation is prevented by the increase of the potential amplitude at low frequency, as shown by the expression of ${\phi }_{{\ell }\omega }$ equation (41). In other words, it is the time derivative of the electric potential that matters. Equation (95) combined with equation (41) leads to the following expression of the time Fourier transform of the momentum flux

$$\begin{eqnarray}&&{\langle {{\rm{\Pi }}}_{r\parallel \omega }^{D}\rangle }_{\theta }={q}_{0}\displaystyle \frac{1+\left(1+\tfrac{{\omega }^{2}}{{\omega }_{t}^{2}}\right)({\sigma }_{{\ell }\omega }+{\rm{i}}{\nu }_{{\ell }\omega })}{1+{\sigma }_{{\ell }\omega }+{\rm{i}}{\nu }_{{\ell }\omega }}{K}_{\perp ,1}^{2}{\rho }_{{i}}^{2}{{\rm{\Pi }}}_{r\theta ,1\omega }^{E},\end{eqnarray} \tag{ 96 }$$

where the up–down symmetry of the turbulent Reynolds stress ${{\rm{\Pi }}}_{r\theta ,-1\omega }^{E}={{\rm{\Pi }}}_{r\theta ,1\omega }^{E}$ has been used (consequence of the parity of ${C}_{r,{\bf{k}}{\ell }}$ and ${C}_{\parallel ,{\bf{k}}{\ell }}$ with ${\ell }$). The next step consists in taking the zero frequency limit of equation (96), which is equivalent to a time average. For frequencies lower than the transit frequency, $\omega \ll {\omega }_{t}$, the frequency shift is close to unity, ${\sigma }_{{\ell }\omega }=1+o({\omega }^{2}/{\omega }_{t}^{2})$, while the damping rate ${\nu }_{{\ell }\omega }\simeq \omega /{\omega }_{t}$ vanishes. The expression equation (96) becomes equation (44).

Assuming an electron adiabatic response, and using the solution of the gyro-kinetic equation (74), the following electroneutrality equation is obtained for a non axisymmetric mode $n\ne 0$

$$\begin{eqnarray}&&\left\{\tau +1-\displaystyle \frac{1}{N}\displaystyle \int {{\rm{d}}}^{3}{\bf{v}}{F}_{M}{{ \mathcal J }}^{-1}\cdot \displaystyle \frac{\omega -{\omega }_{* }}{\omega -{K}_{\parallel }{v}_{\parallel }-{\omega }_{D}}{ \mathcal J }\cdot \right\}{\tilde{\phi }}_{{\bf{k}}\omega }=0,\end{eqnarray} \tag{ 97 }$$

where $\tau =\tfrac{{T}_{i}}{{T}_{e}}$ is the ratio of the ion to electron temperature at $r={r}_{0}$. This equation, which describes reasonably well toroidal ITG modes, can be easily extended to non adiabatic electrons [46] (see also [7] and references therein). In the hydrodynamic limit

$$\begin{eqnarray}&&\displaystyle \frac{{\omega }_{D}}{\omega }\sim \displaystyle \frac{{K}_{\parallel }^{2}{v}_{{Ti}}^{2}}{{\omega }^{2}}\sim {K}_{\perp }^{2}{\rho }_{i}^{2}\ll 1\end{eqnarray} \tag{ 98 }$$

the electro-neutrality condition equation (97) can be expanded at first order in $\tfrac{{\omega }_{D}}{\omega }$, $\tfrac{{K}_{\parallel }^{2}{v}_{{Ti}}^{2}}{{\omega }^{2}}$ and ${K}_{\perp }^{2}{\rho }_{i}^{2}$, where ${\rho }_{i}=\tfrac{{m}_{i}{v}_{{Ti}}}{{{eB}}_{0}}$ is the thermal ion gyroradius. The equation that rules ${\tilde{\phi }}_{{\bf{k}}\omega }$ now reads [47]

$$\begin{eqnarray}&&\left[-\displaystyle \frac{{\omega }_{t}^{2}}{{\omega }^{2}}\displaystyle \frac{{\partial }^{2}}{\partial {\eta }^{2}}-{K}_{\theta }^{2}{\rho }_{i}^{2}(1+{s}_{0}^{2}{\eta }^{2})+{\rm{\Lambda }}(\omega ,\eta )\right]{\tilde{\phi }}_{{\bf{k}}\omega }=0,\end{eqnarray} \tag{ 99 }$$

where

$$\begin{eqnarray}&&{\rm{\Lambda }}(\omega ,\eta )=\displaystyle \frac{\tau \omega +{\omega }_{* n}}{{\omega }_{* p}-\omega }+\displaystyle \frac{{\omega }_{d}(\eta )}{\omega }\end{eqnarray} \tag{ 100 }$$

is the η dependent local dispersion relation. The transit frequency is defined ${\omega }_{t}=\left|\tfrac{{v}_{{Ti}}}{{q}_{0}{R}_{0}}\right|$ and ${\omega }_{* p}$ is the pressure diamagnetic frequency

$$\begin{eqnarray}&&{\omega }_{* p}=-\displaystyle \frac{{K}_{\theta }{T}_{i}}{{{eB}}_{0}{L}_{{p}_{i}}},\end{eqnarray} \tag{ 101 }$$

where ${L}_{{p}_{i}}$ is the pressure gradient length. The next step consists in expanding the magnetic drift frequency near $\eta =0$ up to $o({\eta }^{2})$ and $o({\theta }_{k}^{2})$. After regrouping the various terms, one finds

$$\begin{eqnarray}&&\left[-\displaystyle \frac{{\omega }_{t}^{2}}{{\omega }^{2}}\displaystyle \frac{{\partial }^{2}}{\partial {\eta }^{2}}+\left(\displaystyle \frac{{\omega }_{d}}{\omega }\left({s}_{0}-\displaystyle \frac{1}{2}\right)-{K}_{\theta }^{2}{\rho }_{i}^{2}{s}_{0}^{2}\right){(\eta -{\eta }_{{\bf{k}}})}^{2}+{\rm{\Lambda }}\right]{\tilde{\phi }}_{{\bf{k}}\omega }=0,\end{eqnarray} \tag{ 102 }$$

where Λ is now a local dispersion relation independent of η and

$$\begin{eqnarray}&&{\rm{\Lambda }}(\omega )=\displaystyle \frac{\tau \omega +{\omega }_{* n}}{{\omega }_{* p}-\omega }+\displaystyle \frac{{\omega }_{d}}{\omega }-{K}_{\theta }^{2}{\rho }_{i}^{2}-\displaystyle \frac{1}{4}\displaystyle \frac{{\omega }_{d}^{2}}{{\omega }^{2}}\displaystyle \frac{{({s}_{0}-1)}^{2}}{\left[\tfrac{{\omega }_{d}}{\omega }\left({s}_{0}-\tfrac{1}{2}\right)-{K}_{\theta }^{2}{\rho }_{i}^{2}{s}_{0}^{2}\right]}{\theta }_{k}^{2},\end{eqnarray} \tag{ 103 }$$

where ${\omega }_{d}$ is the local magnetic drift frequency. Strictly speaking ${\omega }_{d}=-2\tfrac{{K}_{\theta }T}{{{eB}}_{0}{R}_{0}}(1-{\theta }_{k}^{2})$, but the analysis is restricted here to low values of ${\theta }_{k}$. Calculations are run at first order in ${\theta }_{k}$, so that the magnetic drift is the one defined in equation (51). The complex number ${\eta }_{{\bf{k}}}$ is proportional to the ballooning angle ${\theta }_{k}$

$$\begin{eqnarray}&&{\eta }_{{\bf{k}}}=-\displaystyle \frac{1}{2}\displaystyle \frac{{\omega }_{d}}{\omega }\displaystyle \frac{({s}_{0}-1)}{\left[\tfrac{{\omega }_{d}}{\omega }\left({s}_{0}-\tfrac{1}{2}\right)-{K}_{\theta }^{2}{\rho }_{i}^{2}{s}_{0}^{2}\right]}{\theta }_{k}.\end{eqnarray} \tag{ 104 }$$

The smoothest solution of this differential equation is

$$\begin{eqnarray}&&{\tilde{\phi }}_{{\bf{k}}\omega }={\phi }_{0}\exp \left\{-\displaystyle \frac{1}{2}{\alpha }_{{\bf{k}}}{(\eta -{\eta }_{{\bf{k}}})}^{2}\right\}\end{eqnarray} \tag{ 105 }$$

with the conditions

$$\begin{eqnarray}&&\displaystyle \frac{{\omega }_{t}^{2}}{{\varpi }_{{\bf{k}}}^{2}}{\alpha }_{{\bf{k}}}^{2}=\displaystyle \frac{{\omega }_{d}}{{\varpi }_{{\bf{k}}}}\left({s}_{0}-\displaystyle \frac{1}{2}\right)-{K}_{\theta }^{2}{\rho }_{i}^{2}{s}_{0}^{2},\end{eqnarray} \tag{ 106 }$$

$$\begin{eqnarray}&&\displaystyle \frac{{\omega }_{t}^{2}}{{\varpi }_{{\bf{k}}}^{2}}{\alpha }_{{\bf{k}}}=-{\rm{\Lambda }}({\varpi }_{{\bf{k}}}).\end{eqnarray} \tag{ 107 }$$

The two equations (106), (107) provide the values of ${\alpha }_{{\bf{k}}}$ and of the complex mode frequency ${\varpi }_{{\bf{k}}}$. Moreover, the solution equation (105) is acceptable only if ${\mathfrak{R}}({\alpha }_{{\bf{k}}})\gt 0$, to guarantee mode spatial localization. The mode frequency ${\varpi }_{{\bf{k}}}$ is written as ${\varpi }_{{\bf{k}}}={\omega }_{{\bf{k}}}+{\rm{i}}{\gamma }_{{\bf{k}}}$ , where ${\omega }_{{\bf{k}}}$ is the real part, and ${\gamma }_{{\bf{k}}}$ the growth rate. The later is positive above the instability threshold. We will focus on that case. A limit of interest is ${K}_{\theta }{\rho }_{i}\ll 1$, and $\tfrac{{\omega }_{t}^{2}}{{\omega }^{2}}{\alpha }_{{\bf{k}}}\ll {\rm{\Lambda }}$. The dispersion relation then becomes ${\rm{\Lambda }}({\varpi }_{{\bf{k}}})=0$, i.e.

$$\begin{eqnarray}&&\tau {\varpi }_{{\bf{k}}}^{2}-({\omega }_{d}-{\omega }_{* n}){\varpi }_{{\bf{k}}}+{\omega }_{* p}{\omega }_{d}=0\end{eqnarray} \tag{ 108 }$$

which yields

$$\begin{eqnarray}&&{\omega }_{{\bf{k}}}=\displaystyle \frac{1}{2\tau }({\omega }_{d}-{\omega }_{* n}),\end{eqnarray} \tag{ 109 }$$

$$\begin{eqnarray}&&{\gamma }_{{\bf{k}}}=\displaystyle \frac{1}{2\tau }{[4\tau {\omega }_{* p}{\omega }_{d}-{({\omega }_{d}-{\omega }_{* n})}^{2}]}^{1/2}.\end{eqnarray} \tag{ 110 }$$

It is recovered that modes drift linearly in the ion diamagnetic direction for flat density profiles ${\omega }_{* n}\simeq 0$, while they rotate in the electron diamagnetic direction for strong density gradients [47]. We now turn to the quantities of interest, i.e. the radial and parallel wavenumber given by equations (28), (29). Since equation (107) depends on the detail of the local dispersion relation, we will use exclusively equation (106) and leave ${\varpi }_{{\bf{k}}}$ arbitrary, except at the very end, when making estimates for ITG modes. It appears that ${\delta }_{{\bf{k}}}={\lambda }_{{\bf{k}}}-1$ is given by the following relation

$$\begin{eqnarray}&&{\delta }_{{\bf{k}}}=-\displaystyle \frac{\tfrac{{\omega }_{d}}{2{\varpi }_{{\bf{k}}}}({s}_{0}-1)}{\tfrac{{\omega }_{d}}{{\varpi }_{{\bf{k}}}}\left({s}_{0}-\tfrac{1}{2}\right)-{K}_{\theta }^{2}{\rho }_{i}^{2}{s}_{0}^{2}}=-\displaystyle \frac{1}{2}\displaystyle \frac{{\omega }_{d}{\varpi }_{{\bf{k}}}}{{\omega }_{t}^{2}}\displaystyle \frac{1}{{\alpha }_{{\bf{k}}}^{2}}({s}_{0}-1)\end{eqnarray} \tag{ 111 }$$

from which its imaginary part ${\mathfrak{I}}({\delta }_{{\bf{k}}})$ can be deduced

$$\begin{eqnarray}&&{\mathfrak{I}}({\delta }_{{\bf{k}}})=-\displaystyle \frac{1}{2}\displaystyle \frac{{K}_{\theta }^{2}{\rho }_{i}^{2}{s}_{0}^{2}({s}_{0}-1)}{{\left|\tfrac{{\omega }_{d}}{{\varpi }_{{\bf{k}}}}\left({s}_{0}-\tfrac{1}{2}\right)-{K}_{\theta }^{2}{\rho }_{i}^{2}{s}_{0}^{2}\right|}^{2}}\displaystyle \frac{{\gamma }_{{\bf{k}}}{\omega }_{d}}{{| {\varpi }_{{\bf{k}}}| }^{2}}\end{eqnarray} \tag{ 112 }$$

thus leading to

$$\begin{eqnarray}&&{K}_{\parallel ,0}=-\displaystyle \frac{1}{2{q}_{0}{R}_{0}}{K}_{\theta }^{2}{\rho }_{i}^{2}{s}_{0}^{2}({s}_{0}-1)\displaystyle \frac{1}{{\mathfrak{R}}({\alpha }_{{\bf{k}}})}{\left|\displaystyle \frac{{\varpi }_{{\bf{k}}}}{{\alpha }_{{\bf{k}}}{\omega }_{t}}\right|}^{2}\displaystyle \frac{{\gamma }_{{\bf{k}}}{\omega }_{d}}{{\omega }_{t}^{2}}{\theta }_{k}.\end{eqnarray} \tag{ 113 }$$

Using the relation ${\omega }_{d}=-2| {q}_{0}| {K}_{\theta }{\rho }_{i}{\omega }_{t}$, one gets the expression of ${C}_{\parallel {\bf{k}}}$

$$\begin{eqnarray}&&{C}_{\parallel {\bf{k}}}=\mathrm{sgn}({q}_{0})\displaystyle \frac{{\rho }_{i}}{{R}_{0}}{({K}_{\theta }{\rho }_{i}{s}_{0})}^{2}({s}_{0}-1)\displaystyle \frac{1}{{\mathfrak{R}}({\alpha }_{{\bf{k}}})}{\left|\displaystyle \frac{{\varpi }_{{\bf{k}}}}{{\alpha }_{{\bf{k}}}{\omega }_{t}}\right|}^{2}\displaystyle \frac{{\gamma }_{{\bf{k}}}}{{\omega }_{t}}.\end{eqnarray} \tag{ 114 }$$

The calculation of ${K}_{r,0}$ is somewhat more delicate because of compensation effects. Using the rhs of equation (111), one finds

$$\begin{eqnarray}&&{C}_{r{\bf{k}}}=\displaystyle \frac{{K}_{r,0}}{{K}_{\theta }{\theta }_{k}}=-\displaystyle \frac{1}{2}\displaystyle \frac{{\omega }_{d}}{{\omega }_{t}^{2}}{s}_{0}({s}_{0}-1)\displaystyle \frac{{\mathfrak{R}}\left(\tfrac{{\varpi }_{{\bf{k}}}}{{\alpha }_{{\bf{k}}}}\right)}{{\mathfrak{R}}({\alpha }_{{\bf{k}}})}.\end{eqnarray} \tag{ 115 }$$

It appears readily that

$$\begin{eqnarray}&&{C}_{r{\bf{k}}}=-\displaystyle \frac{1}{2}{s}_{0}({s}_{0}-1)\displaystyle \frac{{\omega }_{d}{\omega }_{{\bf{k}}}}{{\omega }_{t}^{2}}\displaystyle \frac{{\sigma }_{{\bf{k}}}}{{| {\alpha }_{{\bf{k}}}| }^{2}},\end{eqnarray} \tag{ 116 }$$

where

$$\begin{eqnarray}&&{\sigma }_{{\bf{k}}}=1-\displaystyle \frac{{\mathfrak{I}}({\alpha }_{{\bf{k}}})}{{\mathfrak{R}}({\alpha }_{{\bf{k}}})}\displaystyle \frac{{\gamma }_{{\bf{k}}}}{{\omega }_{{\bf{k}}}}.\end{eqnarray} \tag{ 117 }$$

The condition equation (106) can be used to provide a relation between ${\mathfrak{R}}({\alpha }_{{\bf{k}}})$ and ${\mathfrak{I}}({\alpha }_{{\bf{k}}})$, i.e.

$$\begin{eqnarray}&&2{\mathfrak{R}}({\alpha }_{{\bf{k}}}){\mathfrak{I}}({\alpha }_{{\bf{k}}})=\displaystyle \frac{{\gamma }_{{\bf{k}}}{\omega }_{{\bf{k}}}}{{\omega }_{t}^{2}}\left[\displaystyle \frac{{\omega }_{d}}{{\omega }_{{\bf{k}}}}\left({s}_{0}-\displaystyle \frac{1}{2}\right)-2{K}_{\theta }^{2}{\rho }_{i}^{2}{s}_{0}^{2}\right]\end{eqnarray} \tag{ 118 }$$

from which one gets the following result

$$\begin{eqnarray}&&{\sigma }_{{\bf{k}}}=1-\left[\displaystyle \frac{{\omega }_{d}}{2{\omega }_{{\bf{k}}}}\left({s}_{0}-\displaystyle \frac{1}{2}\right)-{K}_{\theta }^{2}{\rho }_{i}^{2}{s}_{0}^{2}\right]\displaystyle \frac{{\gamma }_{{\bf{k}}}^{2}}{{[{\mathfrak{R}}({\alpha }_{{\bf{k}}})]}^{2}{\omega }_{t}^{2}}.\end{eqnarray} \tag{ 119 }$$

An explicit expression of ${[{\mathfrak{R}}({\alpha }_{{\bf{k}}})]}^{2}$ can be found by solving equation (106) , i.e.

$$\begin{eqnarray}{[{\mathfrak{R}}({\alpha }_{{\bf{k}}})]}^{2} & = & \displaystyle \frac{1}{2{\omega }_{t}^{2}}\left[{\omega }_{{\bf{k}}}{\omega }_{d}\left({s}_{0}-\displaystyle \frac{1}{2}\right)-{K}_{\theta }^{2}{\rho }_{i}^{2}{s}_{0}^{2}\{{\omega }_{{\bf{k}}}^{2}-{\gamma }_{{\bf{k}}}^{2}\}\right]\\ & & +\displaystyle \frac{1}{2{\omega }_{t}^{2}}\left\{{\left[{\omega }_{{\bf{k}}}{\omega }_{d}\left({s}_{0}-\displaystyle \frac{1}{2}\right)-{K}_{\theta }^{2}{\rho }_{i}^{2}{s}_{0}^{2}\{{\omega }_{{\bf{k}}}^{2}-{\gamma }_{{\bf{k}}}^{2}\}\right]}^{2}\right.\\ & & +{\left.{\gamma }_{{\bf{k}}}^{2}{\left[{\omega }_{d}\left({s}_{0}-\displaystyle \frac{1}{2}\right)-2{K}_{\theta }^{2}{\rho }_{i}^{2}{s}_{0}^{2}{\omega }_{{\bf{k}}}\right]}^{2}\right\}}^{1/2}.\end{eqnarray} \tag{ 120 }$$

It results from equation (119) that ${C}_{r{\bf{k}}}$ changes sign when

$$\begin{eqnarray}&&\displaystyle \frac{{\gamma }_{{\bf{k}}}^{2}}{{[{\mathfrak{R}}({\alpha }_{{\bf{k}}})]}^{2}{\omega }_{t}^{2}}\geqslant \displaystyle \frac{1}{\tfrac{{\omega }_{d}}{2{\omega }_{{\bf{k}}}}\left({s}_{0}-\tfrac{1}{2}\right)-{K}_{\theta }^{2}{\rho }_{i}^{2}{s}_{0}^{2}}\end{eqnarray} \tag{ 121 }$$

provided the rhs is positive. The positivity condition is fulfilled for modes drifting in the ion diamagnetic direction $\mathrm{sgn}({\omega }_{d}{\omega }_{{\bf{k}}})$ for low wave numbers and positive magnetic shear. This change of sign typically occurs when one moves from the situation of weak drive (near threshold) ${\gamma }_{{\bf{k}}}\ll {\omega }_{{\bf{k}}}$ to a situation of strong drive ${\gamma }_{{\bf{k}}}\gg {\omega }_{{\bf{k}}}$.

To make these expressions more explicit, we concentrate on the limit of a strong magnetic shear $| {s}_{0}| \gg 1$ and low wavenumbers ${K}_{\theta }^{2}{\rho }_{i}^{2}\ll 1$. Using the relations equations (118) (120), explicit expressions of ${\mathfrak{R}}({\alpha }_{{\bf{k}}})$, ${\mathfrak{I}}({\alpha }_{{\bf{k}}})$ can be found

$$\begin{eqnarray}&&{[{\mathfrak{R}}({\alpha }_{{\bf{k}}})]}^{2}=\displaystyle \frac{1}{2}\displaystyle \frac{| {\omega }_{{\bf{k}}}{\omega }_{d}{s}_{0}| }{{\omega }_{t}^{2}}{\lambda }_{+}\end{eqnarray} \tag{ 122 }$$

$$\begin{eqnarray}&&{[{\mathfrak{I}}({\alpha }_{{\bf{k}}})]}^{2}=\displaystyle \frac{1}{2}\displaystyle \frac{| {\omega }_{{\bf{k}}}{\omega }_{d}{s}_{0}| }{{\omega }_{t}^{2}}{\lambda }_{-}\end{eqnarray} \tag{ 123 }$$

from which ${\sigma }_{{\bf{k}}}$ and ${| {\alpha }_{{\bf{k}}}| }^{2}$ can be derived

$$\begin{eqnarray}&&{\sigma }_{{\bf{k}}}=1-\mathrm{sgn}({\omega }_{{\bf{k}}}{\omega }_{d}{s}_{0}){\lambda }_{-},\end{eqnarray} \tag{ 124 }$$

$$\begin{eqnarray}&&{| {\alpha }_{{\bf{k}}}| }^{2}=\displaystyle \frac{| {\omega }_{{\bf{k}}}{\omega }_{d}{s}_{0}| }{{\omega }_{t}^{2}}{\left(1+\displaystyle \frac{{\gamma }_{{\bf{k}}}^{2}}{{\omega }_{{\bf{k}}}^{2}}\right)}^{1/2},\end{eqnarray} \tag{ 125 }$$

where

$$\begin{eqnarray}&&{\lambda }_{\pm }={\left(1+\displaystyle \frac{{\gamma }_{{\bf{k}}}^{2}}{{\omega }_{{\bf{k}}}^{2}}\right)}^{1/2}\pm \mathrm{sgn}({\omega }_{{\bf{k}}}{\omega }_{d}{s}_{0}).\end{eqnarray} \tag{ 126 }$$

The numbers ${\lambda }_{+}$ and ${\lambda }_{-}$ are always positive, and satisfy the useful relationships

$$\begin{eqnarray}&&{\lambda }_{+}{\lambda }_{-}=\displaystyle \frac{{\gamma }_{{\bf{k}}}^{2}}{{\omega }_{{\bf{k}}}^{2}},\end{eqnarray} \tag{ 127 }$$

$$\begin{eqnarray}&&{\lambda }_{+}+{\lambda }_{-}=2\left(1+\displaystyle \frac{{\gamma }_{{\bf{k}}}^{2}}{{\omega }_{{\bf{k}}}^{2}}\right).\end{eqnarray} \tag{ 128 }$$

The expression equation (125) of ${| {\alpha }_{{\bf{k}}}| }^{2}$ can be obtained directly from equation (106), thus providing a cross-check . Plugging equations (127), (128) in the expressions of ${C}_{\parallel {\bf{k}}}$ and ${C}_{r{\bf{k}}}$ equations (113), (116), one obtains the following relations

$$\begin{eqnarray}&&{C}_{r{\bf{k}}}={c}_{r{\bf{k}}}{s}_{0},\end{eqnarray} \tag{ 129 }$$

$$\begin{eqnarray}&&{C}_{\parallel {\bf{k}}}={c}_{\parallel {\bf{k}}}\displaystyle \frac{{\rho }_{i}}{{R}_{0}}{({K}_{\theta }{\rho }_{i}{s}_{0})}^{2}{q}_{0}{s}_{0},\end{eqnarray} \tag{ 130 }$$

where ${c}_{r{\bf{k}}}$ and ${c}_{\parallel {\bf{k}}}$ are numbers that depend on the normalized frequencies $\tfrac{{\omega }_{{\bf{k}}}}{{\omega }_{d}}$ and $\tfrac{{\gamma }_{{\bf{k}}}}{{\omega }_{d}}$ only, namely

$$\begin{eqnarray}&&{c}_{r{\bf{k}}}=\displaystyle \frac{1}{2}\displaystyle \frac{1}{{\left(1+\tfrac{{\gamma }_{{\bf{k}}}^{2}}{{\omega }_{{\bf{k}}}^{2}}\right)}^{1/2}}({\lambda }_{-}-\mathrm{sgn}({\omega }_{{\bf{k}}}{\omega }_{d}{s}_{0})),\end{eqnarray} \tag{ 131 }$$

$$\begin{eqnarray}&&{c}_{\parallel {\bf{k}}}={2}^{1/2}\displaystyle \frac{1}{| {q}_{0}| }{\left|\displaystyle \frac{{\omega }_{{\bf{k}}}}{{\omega }_{d}{s}_{0}}\right|}^{3/2}{\left(1+\displaystyle \frac{{\gamma }_{{\bf{k}}}^{2}}{{\omega }_{{\bf{k}}}^{2}}\right)}^{1/2}{\lambda }_{-}^{1/2}.\end{eqnarray} \tag{ 132 }$$