The role of radial electric fields in linear and nonlinear gyrokinetic full radius simulations

The pivotal role played by radial electric fields in the development of turbulence associated with anomalous transport is examined by means of global gyrokinetic simulations. It is shown that the stabilizing effect of E×B flows on ion temperature gradient (ITG) modes is quadratic in the shearing rate amplitude. For a given shearing rate it leads to an increase in the critical gradient. The electric fields (zonal flows) self-generated by ITG modes interact in a nonlinear way and it is shown that a saturated level of both the zonal flow and ITG turbulence is reached in the absence of any collisional mechanism being included in the model. The quality of the global nonlinear simulations is verified by the energy conservation which is allowed by the inclusion of nonlinear parallel dynamics. This demonstrates the absence of spurious damping of numerical origin and thus confirms the nonlinear character of zonal flow saturation mechanism.


Introduction
The development of turbulence from underlying micro-instabilities is one of the mechanisms most widely held responsible for anomalous transport in magnetically confined plasmas. Various types of drift-wave like instabilities have been considered, e.g. ion temperature gradient (ITG) modes, trapped electron modes (TEM), trapped ion modes (TIM), electron temperature gradient (ETG) modes and kinetic ballooning modes (KBM, or Alfvén-ITG). Fluid, gyrofluid, or gyrokinetic theories have been developed, most of them within the electrostatic approximation. Various approaches, differing in their degree of sophistication, have been used to solve the models, ranging from local dispersion relations to ballooning approximations to flux tubes to full radius computations.

29.2
While linear theory is relatively well established for most of the instabilities cited above, the nonlinear evolution of these modes is a much more difficult problem. The challenge is not only of a technical nature but also to understand the physics mechanisms at play. Several theoretical works [1]- [16] have shown that one of the pivotal ingredients of the physics is the existence of radial electric fields that are self-generated by the turbulence. These are the E × B 'zonal flows' (ZF). One mechanism for the ZF generation is by modulational instability of the ITG modes. The importance of the ZF is in that they can have a stabilizing effect on the ITG modes. The mechanism usually invoked is that the ZF, if they are sheared, can tear apart the turbulent eddies of the ITG, hence the alias 'sheared E × B flows'. Therefore a nonlinear feedback loop links the ITG modes amplitude and the ZF amplitude.
Since ZF have k = 0 they are not subject to Landau damping. The question of ZF damping is crucial [5] because it affects the ZF amplitude, which in turn determines the ITG turbulence level and therefore the anomalous heat flux. Various linear or nonlinear mechanisms for the damping of ZF have been considered. Some models have examined Kelvin-Helmolz instabilities of the ZF [10], while some others consider collisionality [7,8]. In nonlinear numerical simulations it is therefore important to verify that the discretization does not introduce fake damping or drive mechanisms. This issue is addressed in this paper.
Radial electric fields can also be imposed by external means, e.g. creating sheared toroidal rotation or with an applied electrode bias [17,18]. These external electric fields could play an important role in the formation of transport barriers. There is a generally observed correlation between the appearance of radial electric fields, reduction of turbulence and improved energy confinement [19]. Experimental evidence of the self-regulating role of ZF has recently been observed in both tokamaks [20,21] and stellarators [22].
In this paper we examine two aspects of the role of E × B flows. First, considering the radial electric field as an equilibrium quantity, we consider ITG, TIM and TEM instabilities. We show that an applied E × B flow is generally stabilizing for ITG modes with a quadratic dependence on the shearing rate, whereas E × B flows can destabilize trapped particle modes. Second, we show in a nonlinear full radius gyrokinetic simulation how the system reaches a quasi-stationary state with finite level of both ITG and ZF amplitudes. The energy conservation property is checked and shows the absence of spurious damping or drive of numerical origin. Since our model does not include any linear damping of ZF (Landau or collisions) this points to the essentially nonlinear mechanism for the saturation of ZF.

Gyrokinetic global model
We consider low β magnetic configurations with prescribed equilibrium profiles of safety factor q(ψ), density n 0 (ψ), electron and ion temperatures T e (ψ) and T i (ψ), to which can also be added an electrostatic potential Φ 0 (ψ). The symmetry of the configuration is either axisymmetric or helical. All profiles are taken as function of the poloidal magnetic flux ψ (respectively helical flux). The radial coordinate is defined as s = ψ/ψ s , where ψ s is the edge poloidal flux (respectively helical flux).
The gyrokinetic model is used to describe electrostatic perturbations that satisfy the usual ordering ω/Ω ∼ k /k ⊥ ∼ eδφ/T e ∼ ρ L /L n ∼ ρ L /L T ∼ O( g ), where ρ L is the ion Larmor radius, Ω is the ion cyclotron frequency, L −1 n = |∇ ln n 0 |, L −1 29.3 f = f 0 + δf and the electrostatic potential as Φ 0 + δφ the equations read [23], after linearization: where the brackets indicate Larmor averaging. We consider small Mach numbers of the E × B flow and thus terms of order (v E /v thi ) 2 and B (v E /v thi ) have been neglected. The system of equations is closed by the quasi-neutrality condition Without loss of generality the perturbations are written as where θ is the poloidal angle, ϕ is the toroidal angle and χ is a straight field line poloidal coordinate. Taking advantage of the fact that perturbations have k k ⊥ , the parameters s 0 , m 0 and ω 0 are chosen such that the transformed quantitiesδf andδφ have much smoother poloidal and time dependences. This proves to considerably improve the computational accuracy [24]. For electrons the Larmor radius is neglected. Passing electrons are assumed to respond adiabatically. Non-adiabatic trapped electrons are considered as drift-kinetic.
For nonlinear simulations in this paper we shall neglect magnetic shear and magnetic curvature effects. Electrons are assumed to respond adiabatically along the magnetic field lines. The corresponding equations are where n i is the gyro-averaged ion density andφ is the magnetic surface average of the potential. Note that the parallel nonlinearity is retained in our model. This implies that energy conservation is satisfied: 29.4 The turbulence-driven average radial energy flux is While the nonlinearly driven flux is the most interesting physical quantity in the context of anomalous transport studies, the energy conservation property of the equations can be used as a strong test of the quality of the numerical simulation [25]. If satisfied, it shows that there is no spurious damping or drive mechanism of numerical origin. The equations are discretized using quasi-particles (gyrocentre tracers) and finite elements on a magnetic coordinate system [24]. In [25] several techniques are described to improve the quality of the numerical scheme. In particular it is shown that energy conservation can be satisfied with this PIC-δf scheme.
The stabilization of ITG modes by the shearing effect of E × B flows is often expressed by the criterion |ω E×B | > γ 0 , where ω E×B is the shearing rate and γ 0 is the maximum linear growth rate in the absence of flow. The shearing rate is best described by the expression [3] ω E×B = ∆ψ ∆ϕ in which ∆ψ/|∇ψ| and R∆ϕ are the radial and toroidal correlation lengths of turbulence, respectively. Note that strictly speaking the turbulence suppression criterion [3] is when |ω E×B | is comparable to the inverse turbulence correlation time. For application purposes, however, this is often replaced by the linear growth rate. The last expression in equation (15) is valid in the circular, axisymmetric, large aspect ratio limit and ρ is the minor radius. It shows that the shearing rate is not simply the shear in the poloidal E × B velocity but there is a contribution from the magnetic shear combined with the value of the poloidal E × B velocity.

Applied E × B flows
In this section we focus on the way applied E × B flows can influence the ITG growth rate using a model with adiabatic electrons. The question of the non-adiabatic trapped electron response will be examined in the following section. It will be shown that for toroidal-ITG, slab-ITG and helical-ITG modes, the stabilizing effect of the flows is essentially quadratic in the shearing rate. However, marginal stability is reached when |ω E×B | is comparable to γ 0 within a factor of about 2. In [26] a detailed analysis of this effect is made in both axisymmetric and helically symmetric configurations.
To show the generality of this behaviour let us consider three different magnetic configurations.
(1) A circular cross-section tokamak of aspect ratio 5.5 with a q profile given by q(s) = 1.25 + 3s 2 . (2) A helically symmetric heliac configuration with an elongation of 2 and shearless q profile.
In what follows we shall refer to these configurations as 'tokamak', 'heliac' and 'cylinder' for the sake of simplicity. In all these plasmas we consider T e = T i profiles given by the heliac a/L T = 1.3 and in the cylinder a/L T = 4. The toroidal mode number n (respectively helical mode number and longitudinal wavenumber) is chosen so that it corresponds to the most unstable ITG mode of the corresponding type. The respective average values of k θ ρ Li are k θ ρ Li ≈ 0.5 for the tokamak and heliac cases, and k θ ρ Li ≈ 0.8 for the cylinder case. These values are slightly modified by the presence of applied radial electric fields examined in this paper.
Three different types of profiles of radial electric field are considered.
where B c is the magnetic field on magnetic axis, a is the average minor radius and M , k x and ∆ x are dimensionless parameters. While profiles (a) and (b) are very global, profiles (c) resemble the ZF self-generated by the turbulence, as will be shown in section 5. It should be noted that there is nevertheless a fundamental difference between applied E × B flows, which are treated as equilibrium quantities in this section, and the ZF, which are fluctuating quantities. Figure 1 shows that for all the considered cases there is a quadratic dependence of the resulting ITG growth rate on the E × B shearing rate. Note that throughout this paper the definition used for the shearing rate is that of equation (15), and although it was established 29.6 in [3] in the context of turbulence decorrelation by ZF, it will be shown below that the criterion for ITG stabilization in the presence of applied (equilibrium) E × B flows can also be expressed with the same definition for the shearing rate. The circles correspond to the tokamak case while the diamonds are results for the heliac case. The filled circles have shearless v E profiles (type (a) defined above), while for the three other cases linear profiles (type (b) defined above) have been specified. The curve with filled diamonds corresponds to a helical-ITG mode while the open diamonds are slab-like ITG modes in the heliac. For each of these four cases a parabolic fit of the global gyrokinetic results is shown with dashed curves.
We conclude that the effect of the shearing rate of radial electric fields on ITG modes is essentially quadratic, and that this quadratic nature is a generic feature. As a matter of fact all our simulations of ITG modes in the presence of E × B flows show this behaviour when trapped particle effects do not dominate the instability drive. Trapped particle effects seem to modify this simple behaviour and an example will be shown in the next section.
Looking at the results of figure 1 there appears an asymmetry with the sign of the shearing rate. This asymmetry is due to the fact that the eigenmode in the absence of E × B flow is not up-down symmetric. With one sign of the shearing rate the mode is first straightened so that the radial structure of the mode is aligned in the most unfavourable grad-B drift direction, and therefore the growth rate is maximized; increasing the shearing rate further tilts the radial structure and stabilizes the mode, in a similar mechanism as shown in [27,28]. For all the cases shown in figure 1 the ITG modes are stabilized when the shearing rate is approximately equal to the growth rate without flow: |ω crit E×B | ≈ γ 0 seems therefore to hold, at least within a factor of 2. This should not mislead us to assume that the ITG growth rate in the presence of flow is diminished linearly as γ 0 − |ω E×B |: as we have just shown the ITG growth rate in the presence of flows is quadratically dependent on ω E×B .
Neglecting now the effects of magnetic curvature and magnetic shear we consider a cylindrical configuration with constant axial magnetic field. We set a/L T = 4 and η i = 5 and study the effect of the radial wavelength of E × B flows on ITG mode stability. The results (figure 2) show that considering a very global profile of the shearing rate (ω E×B ∼ ρ/a), the open squares in figure 2, is most effective to stabilize the ITG. For profiles given by equation (16) the effect decreases with the radial wavenumber k x of the E × B flow. The filled triangles in figure 2 are for k x = 5, ∆ x = 0.3 and the filled squares are for k x = 10, ∆ x = 0.15.
The dashed curves in figure 2 are parabolic fits to the data for absolute values of the shearing rate smaller than 0.15: there is an extremely good agreement with the results of the global gyrokinetic simulations, once more showing the generality of the quadratic effect of the shearing rate on ITG modes. For large values of the shearing rate, however, we observe a deviation from this quadratic behaviour. A detailed analysis of the stabilizing effect of E × B shearing rate has been made for the case k x = 10. Figure 3 shows the contribution of the E × B flow to the particle-wave rate of power exchange, q i v E · δE f d 6 z, plotted as a function of the shearing rate. The quadratic stabilizing effect of the shearing rate is well evidenced for small enough values. The stabilizing effect is seen to saturate for large values of the shearing rate (|ω E×B |(a/c s ) >≈ 0.2).
In order to better understand the reason for this behaviour we show in figure 4 the radial profile of the shearing rate for the case k x = 10, ∆ x = 0.15, ω E×B (a/c s ) = 0.1 and the radial position of the maximum ITG mode amplitude as a function of the applied shearing rate: a clear radial shift occurs, pushing the mode away from the maximum shearing rate position. It also appears that this radial shift saturates. It is interesting to plot the ITG growth rate not as a  (16) with k x = 5. Filled squares: for a flow profile given by equation (16) with k x = 10. The dashed curves are parabolic fits to the respective data for small absolute values of the shearing rate.  The result is shown if figure 5. The dashed curve is a parabolic fit to the data. It shows that, for a good part, the deviation from a parabolic dependence of the growth rate on the shearing rate can be attributed to the radial shift of the ITG mode towards regions of smaller shearing rate.
Keeping the shape of the E × B profile with k x = 10 we performed a double parameter scan in the ITG and in the E × B shearing rate amplitude. The growth rates for various values of a/L T are shown in figure 6 as a function of the shearing rate. Except for large values of the shearing rate the dependence of the growth rate is quadratic. In figure 6 negative growth rates are also shown. The question of damped modes is important in relation to the process of nonlinear saturation because they act as an energy sink. In [29] a spectral approach was used to investigate the stable ballooning modes, and it was shown that, in addition to normal modes having an exponential decay there is a continuum mode with a 1/t 2 asymptotic behaviour. For our results with negative γ in figure 6 an exponential decay of the mode amplitude was observed in the simulations, indicating that these are normal modes. The continuum mode is either absent, or may appear only when the amplitude has reached below the noise level inherent to the time-evolution, PIC, δf method used. We note in figure 6 a smooth continuation between the unstable and stable sides; a similar smooth continuation was also obtained in [29], albeit in a different parameter scan, for different types of modes (we do not have ballooning modes in figure 6) and obtained with a different method.
The ITG mode growth rates are plotted in figure 7 versus a/L T for different values of the E × B shearing rate. It shows clearly that there is an upshift in the critical gradient with increasing values of |ω E×B |. The value of the upshift depends on the shearing rate. For all values of the shearing rate the ITG growth rates increase linearily with a/L T above the marginal value, with some levelling off at very high growth rates. We note that the growth rate depends almost exclusively on (a/L T ) − (a/L T ) crit . This is compatible with the results of nonlinear simulations [9] which have shown that self-generated E × B flows can virtually suppress the turbulence when the temperature gradient is slightly above the critical value in the absence of flow. These nonlinear simulations also show that the turbulence suppression mechanism works only up to an upshifted critical gradient value. The critical gradient for marginal ITG stability is shown in figure 8 as a function of the E × B shearing rate at s = s 0 = 0.5 (top). It is found that the critical gradient upshift is proportional to the square of the shearing rate for small values of the shearing rate. Then, due to the saturation of the E × B stabilization (see figure 3), the upshift deviates noticeably from this dependence. As noted earlier, see figures 4-5, this is mainly due to the fact that the E × B flow is pushing the radial position of maximum ITG amplitude, s max , towards lower shearing rate regions. When plotted against ω E×B at s = s max (figure 8, bottom), the quadratic dependence of the critical gradient is more evident.

Applied E × B flows: trapped particle effects
In this section we address the question of the role of trapped ions and trapped electrons in the presence of E × B flows. In the previous section the modes under study were toroidal-ITG, helical-ITG or slab-like ITG. We now consider a case for which the most unstable mode, in the absence of radial electric field, is a TIM [30]. We include in our model the non-adiabatic response of trapped electrons which are modelled as drift kinetic. The configuration is a tokamak with the following parameters: R 0 = 1.5 m, a = 0.5 m, B 0 = 1 T, hydrogen ions, T i = T e , T 0 = 5 keV, a/L T = 3.33, a/L n = 0.33, s 0 = 0.7, with a q profile such that q(s 0 ) = 1.5 and the magnetic shear at s = s 0 is unity.
We consider an applied radial electric field with a shearless v E profile, v E = M v thi ρ/ρ 0 , where M is the Mach number of the E × B flow at s = s 0 , and study the instabilities as a function of M . (Note that although this is a shearless v E profile the shearing rate is not zero, see equation (15), because of the finite value of v E and the magnetic shear). The results of our 29.11 TIM is first stabilized but then another mode is destabilized to even higher growth rates than the case without flow. The frequency of the mode, of opposite sign, shows that this mode is rotating in the electron diamagnetic direction. We have checked that this mode is absent if the trapped electron response is neglected. We conclude that it is a TEM. The TEM and TIM have very similar eigenmode structures, see figure 10. Both modes have a maximum amplitude on the unfavourable magnetic curvature region and a ballooning-like appearance, which is, by the way, also the case for toroidal-ITG modes. Comparing the TIM and the TEM, which are the most unstable modes at different Mach numbers of the E × B flow, we notice an effect of shearing the radially extended structures ('fingers'). There is also a small outward shift in the radial position same behaviour as the growth rate. Note that the critical temperature gradient of the TEM with flow is lower as compared to the case without flow. We conclude that the inclusion of E × B flows in this case brings an overall destabilization and that this effect is essentially due to trapped electron dynamics.

Nonlinear interaction of ITG and zonal flows
Turning now to the global nonlinear gyrokinetic model, equations (7)-(10), we focus our attention to the role of ZF on the saturation mechanism of turbulence. As mentioned in the introduction the ITG modes and the ZF modes (m = 0, n = 0) form a coupled nonlinear system that offers  Figure 12 shows the time evolution of the perturbed gyrokinetic ion density for some selected Fourier mode components. In particular, the Fourier component corresponding to the linearly most unstable ITG mode (m = 24, n = 1) is plotted in red, whereas the ZF contribution (m = 0, n = 0) is plotted in green. The perturbation contains a broad Fourier spectrum of modes, which are not shown here for the sake of clarity of the figure. As expected, the early stage of the simulation is characterized by the exponential growth of the linearly most unstable ITG mode. Then the ZF amplitude starts to pick up, with an instantaneous growth rate that exceeds the ITG linear growth rate by a factor up to about 2. This phase thus corresponds to the nonlinear generation of ZF by ITG modes. At the time of maximal ITG amplitude, t = 1.2 × 10 −4 s, the ZF amplitude is high enough to have a stabilizing effect on the ITG modes, in agreement with the studies presented in section 3. This reduces the ITG amplitude. The nonlinear evolution is also characterized by a significant broadening of the spectrum of ITG modes: the linearly most unstable ITG mode is no longer the dominant mode. On a slow timescale, t = 1.2-4×10 −4 s, the ITG amplitude decreases. For t > 4 × 10 −4 s both ITG and ZF amplitudes are nearly constant. Figure 13 (top) shows that the average turbulent-driven energy flux, equation (14), also reaches a quasi-steady-state, after a phase of rapid growth followed by saturation. In figure 13 (bottom) we show the energy balance relation, ∆(E f + E k )/E f , see equations (11)-(13), as a function of time. The energy conservation is satisfied to within 20% of the field energy for the whole duration of the simulation. There is no indication of divergence and no numerical noise accumulation. In particular, in the nonlinear, quasi-steady-state phase of the time evolution, energy is conserved to better than 10%. Satisfying this property is an important verification of the quality of the numerical simulation. However, the implication is also crucial for the nature of the physics processes: in our model the ZF are linearly completely undamped, neither Landau nor collisionally. The energy conservation indicates that there is no spurious damping of numerical origin. But it is well known [31] that in the absence of any ZF damping the dynamic evolution should exhibit one pulse of turbulence followed by ZF saturation at finite amplitude while the turbulence level is completely quenched. This is clearly not what we observe in our gyrokinetic simulations. Therefore, one is left to explain the existence of a quasi-steady-state with finite amplitude of both ITG modes and ZF through an essentially nonlinear ZF saturation mechanism. The evolution of the ZF radial profile during the fast growth phase (t < 1.2 × 10 −4 s) corresponds to a shearing rate profile that is maximum at s = 0.3, which is also the position of the ITG mode maximum amplitude. The region of finite shearing rate covers the radial extent of the ITG, i.e. s ≈ 0.2-0.4. During the nonlinear phase of the simulation (t > 1.2 × 10 −4 s) there is a definite broadening of the region with finite shearing rate, s ≈ 0.1-0.55. This broadening corresponds to that of the ITG fluctuations which are seen to expand in the radial direction as compared to the linear eigenmodes. As an illustration of a typical shape of ZF profile we show in figure 14 the instantaneous profile of ZF shearing rate ω E×B (s) at time t = 3.5 × 10 −4 s. Note that the position of maximum shearing rate is fluctuating back and forth around s = 0.3. The value of the maximum shearing rate is also fluctuating around a value close to the linear growth rate of the most unstable ITG mode (the dashed line in figure 14).
The quasi-stationary phase of the simulation is therefore characterized by nearly constant volume-average amplitudes, but locally fluctuating, of both ITG turbulence and ZF. The ITG and ZF thus form a self-regulated dynamical system, and while the ZF does not contribute directly to the anomalous heat flux it is instrumental in bringing the ITG amplitude to a finite level [1]- [10].

Conclusions
The full radius gyrokinetic simulations presented in this paper prove to be a useful tool in the understanding of the effects of applied E×B flows on ITG modes and of the complex mechanisms involved in the interactions of ITG modes with ZF. To sum up the main results obtained: • The effect of applied E × B on the growth rate of ITG modes is quadratic in the shearing rate (see figures [1][2][3][4][5][6]. Note that this is in line with fluctuation-flow models [32,33]. • For large and localized shearing rates the stabilizing effect seems to saturate (see figure 3), largely because the E × B flow appears to push the ITG mode to radial positions where the shearing rate is lower (see figures 4 and 5). • The behaviour of the ITG growth rate as a function of the ion temperature gradient (see figure 7) for different values of the E × B shearing rate yields an upshift of the critical temperature gradient that is quadratic in the shearing rate (see figure 8). • The inclusion of trapped particle effects can lead to a more complex behaviour: sometimes the effect of sheared E × B flows is destabilizing (see figures 9-11). • Energy conserving nonlinear full radius gyrokinetic simulations, figures 12-14, in which a quasi-stationary state is reached with finite values of both ITG turbulence and ZF in the absence of any linear damping of ZF, confirm the existence of a nonlinear saturation mechanism of ZF. Note that this is also in line with fluctuation-flow models [31]- [33].
In this paper we have considered separately the equilibrium E × B flows and self-generated E × B flows (ZF). The next step will be to consider the presence of both. In other words, it will be interesting to study the dynamical evolution of ZF and ITG modes in the presence of an externally applied E × B flow. This could be of interest for the study of transport barrier formation.