Impact of poloidal convective cells on momentum flux in tokamaks

Radial fluxes of parallel momentum due to E × B and magnetic drifts are shown to be correlated in tokamak plasmas. This correlation comes from the onset of poloidal convective cells generated by turbulence. The entire process requires a symmetry breaking mechanism, e.g. a mean shear flow. An analytical calculation shows that anti-correlation between the poloidal and parallel components of the turbulent Reynolds stress results in anti-correlation of the fluxes of parallel momentum generated by E × B and curvature drifts.


Introduction
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É 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][12][13][14][15][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É 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É 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É 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É B drift contribution to the parallel (or toroidal) momentum flux was calculated [14,[19][20][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É 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É 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É 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É 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É 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].
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É B Reynolds stress translates in a correlation of the curvature andÉ 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 Figure 1. Schematic chart that illustrates the link between correlated radial fluxes of parallel momentum due toÉ B and magnetic drifts, and the correlation between the poloidal and parallel components of the turbulence stress tensor. amplitudes for theÉ 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É 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.

Fluxes of parallel momentum
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.
where F is the distribution function, v  is velocity that is aligned with the unperturbed magnetic field, v E is the E×B drift velocity and 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.
where f is the electric potential, = b B B and 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 r * , where m is the mass, e the charge, andv is the modulus of the perpendicular velocity. A simplified geometry of circular concentric magnetic surfaces is considered here. The spatial coordinates are q j ( ) r, , , where r is the minor radius, θ and j the poloidal and toroidal angles (see figure 2) and the magnetic field is where R 0 the major radius, ( ) q r the safety factor, and B 0 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][25][26][27][28][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 j q -( ) q r instead of the toroidal angle. The periodicity in the toroidal direction allows using a a Fourier series is found to be separable in Q k and θ for each n and pulsation ω, Separability is not always rigorously demonstrated, but is a good proxy. The amplitude F n is called envelope, usually a localized around an angle q k called ballooning angle, with a narrow width of the order of r -( ) L p i 1 2 , where L p is a gradient length and r i the ion thermal gyroradius radius. The ballooning angle q k measures the departure from up-down symmetry on a given field line. Non zero values of q k result from symmetry breaking mechanisms, as will be seen later on. This leads to a slow radial dependence of f q n come from the radial dependence of the safety factor ( ) q r in the exponential argument of equation (6), with a characteristic scale ( ) 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 ( ) 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 The vector k is a label for the couple q ( ) n, k so that the sum over k designates a summation over the toroidal wave number n and also an integral over a distribution of angles q k . The value of the angle q 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 k will be called 'Fourier' harmonics, since q k can be seen as the Fourier counterpart of the radial variable ( ) nq r . Each harmonic f q ( ) t , k depends also slowly on the radius r. We will use in the following a 'slow' radial variable e = r r 0 , 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 contains the information on the poloidal localization of fluctuations, and therefore poloidal asymmetries. We postulate a Gaussian form for i.e. f q where a k and l k are complex numbers, and a k is a normalizing factor (a k , l k and a k depend on r 0 ).

Ballooning angle
The value of the ballooning angle q 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 q = 0 k will acquire an effective radial wave number due to a shear flow that is equal to after a time t. Here ¢ V E is the shear flow rate, i.e. the radial derivative of the meanÉ B velocity. Defining a correlation time t c , one gets the following estimate Other expressions of q k have been given in the past [26][27][28][29]31], with similar structures as equation (10) when envelopes are spatially localized. The exact expression of q 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 q k , but will cause correlated radial and parallel wavenumbers.

E×B drift stress tensor
We estimate now the radial fluxes of poloidal and parallel momentum due toÉ B drift, which are proportional to the corresponding components of the Reynolds stress ò ò The distribution function is decomposed in the same way as the potential is the adiabatic invariant. The equations equations (11), (12) can then be recast as is the Fourier component of the parallel fluid velocity and is the Fourier component of the radial component of theÉ B drift velocity. The wavenumber q K is a reference poloidal wave number defined as = - The amplitude  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 k  is obtained by a rapid distortion argument [30] where t k is a correlation time. The components of the Reynolds stress can be expanded as Fourier series in the poloidal direction Using the expressions equations (23), (24), one finds measure the distortion of turbulent structures in the radial and parallel directions. They can be formally written as (see appendix A) can be seen as effective radial and parallel wavenumbers for a given set ℓ ( ) k, .

Wave numbers Calculations in appendix A show that the coefficients
are close to their = ℓ 0 components for small ballooning angles q  . These coefficients are related to the parameters used for the turbulent field representation R R . The symbol R (resp. I) designates the real (resp. imaginary) part of a complex number.
In view of the field structure equation (9), the real part of a k must be positive, i.e. R a ( )  0 k , since modes are spatially localized in the poloidal direction. As expected, the average radial and parallel wave numbers are correlated and proportional to q 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 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 are even functions of θ. Here calculations are restricted to order one in q 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 P r 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][34][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 P q ℓ

Flux of parallel momentum due to magnetic drift
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 thé E B Reynolds stress, and therefore with theÉ B radial flux of parallel momentum.

Generation ofÉ 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) where N is the unperturbed density,  the gyroaverage operator and - . It is reminded that the inverse aspect ratio r R 0 0 is a small parameter in the present work. To simplify the calculation, we use the long wavelength limit of the gyroaverage operator   1, i.e. the polarization term is kept, but not the finite Larmor radius (FLR) effects. For non zonal modes f á ñ = y 0, this gives an equation for the potential vorticity f r f -î 2 2 , i.e.
where G  is the ion parallel flux, and * 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 where ¹ ℓ 0 since zonal flows are excluded. One can show the following variant of the Taylor identity [39] (see appendix C) 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 ( where the Reynolds stress tensor has been expanded in a Fourier series and^ℓ K , 2 has to be understood as an operator that acts on P q w However it appears that their damping rate is small at low frequencies w w t  , hence favoring their onset and sustainment. It appears that the w 1 scaling of the frequency spectrum plays an important role in the generation of parallel flow via the magnetic drift. Obviously the limit w = 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 w 1 infrared singularity via resonance broadening.

The value of
 plays also an important role. First, it is assumed that r ℓ K i , 2 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 r p 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 0 correspond to scales larger than the distance between resonant surfaces, itself of the order of the ion gyroradius. A more accurate calculation for r p 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]). as discussed in [12]. Numerical simulations suggest that r p is rather in the range of 0.1 [13].

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 f w ℓ 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 For a ballooned turbulence, the Reynolds stress is ballooned too. If fluctuations are strongly localized on the low field side, the prefactor q ( ) cos 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 k , with ℓ (see equation (26)). Given the simplicity of equation (44), one may actually wonder whether a more direct calculation of P r 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 An integration by part allows a reformulation in terms of the parallel gradient of the parallel flux 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].
This expression can be compared with theÉ B flux of momentum equation (77) It appears that these two contributions to the total flux of parallel momentum are anti-correlated when C r k and C k  are of the same sign (for positive safety factor q 0 ). This is also the condition for anti-correlation of thé E 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 k and C k  , and the corresponding signs of fluxes.

Application to ITG driven turbulence
An instructive example is the case of an unstable linear toroidal ITG driven mode. General expressions of C r k and C k  are derived in the hydrodynamic limit in the appendix E. In the long wavelength limit r q K s 1 i 2 2 0 2  and strong magnetic shear | | s 1 0  , explicit expressions of the coefficients C r k and C k  can be found The numbers c r k and c k  depend on the normalized frequencies  , the two numbers c r k and c k  are of the same order of magnitude. The two fluxes are therefore comparable when r r r w t p q . This condition appears as reasonable since the intensity wave number spectrum typically peaks at r [44,45]. Regarding the signs, it appears that for modes that drift in the ion diamagnetic direction w w = ( ) sgn 1 d k , c r k is negative for a weak drive g w k k  , and positive for strong drive g w  k k . Since c k  is always positive, this means that C k  and C r 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 g w k k  and anticorrelated for strong drive g w  k k . For modes drifting in the electron diamagnetic direction, , 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 g w k 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É B flux of parallel momentum is anti-correlated with the ballooning angle q k , which suggests positive c 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 Er pr E is roughly equivalent to multiplying c r k in equation (52) by a factor -w w ( ) * 1 p k , where w * p is the pressure diamagnetic frequency. Since the later tends to be larger than the mode pulsation w k for realistic parameters, this may change the sign of the coefficient c r 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 k and c k  are determined by non linear processes and not well captured by a linear analysis.

Conclusion
It is shown here that theÉ 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É 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É 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 anticorrelation 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.

Appendix A. Mode structure
The mode structure equation (9) lead to a turbulence intensity that reads The amplitude a k in equation (9) has been chosen such that The poloidal Fourier components of the radial and parallel wave numbers equation (25) and equation (26) can then be recast as ò ò q q q q q q l q q   The solution of the gyrokinetic equation (67) involves an integro-differential operator that relates w G k to f w k [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 w w We assume that the spectral turbulence intensity is of the form 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 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 f w ℓ equation (41). In other words, it is the time derivative of the electric potential that matters. Equation with ℓ). 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, w w t  , the frequency shift is close to unity, s w w = +