Impact of Energetic-Particle-Driven Geodesic Acoustic Modes on Turbulence

The impact on turbulent transport of geodesic acoustic modes excited by energetic particles is evidenced for the ﬁrst time in ﬂux-driven 5D gyrokinetic simulations using the GYSELA code. Energetic geodesic acoustic modes (EGAMs) are excited in a regime with a transport barrier in the outer radial region. The interaction between EGAMs and turbulence is such that turbulent transport can be enhanced in the presence of EGAMs, with the subsequent destruction of the transport barrier. This scenario could be particularly critical in those plasmas, such as burning plasmas, exhibiting a rich population of suprathermal particles capable of exciting energetic modes.

Understanding turbulent transport is crucial in numerous plasma physics frameworks, ranging from plasma laboratories such as nuclear fusion devices [1] to astrophysical systems such as the solar tachocline [2] or the atmospheres [3].In this Letter, we focus on the turbulent transport in toroidal nuclear fusion devices (tokamaks), where accurate predictions are essential on the route towards the steadystate production of energy.Together with turbulence, energetic particles (EPs) constitute a ubiquitous component of current and future tokamaks, due to both nuclear reactions and heating systems.EPs are characterized by energies larger than the thermal energy.Whereas the impact of turbulence on EP transport has been analyzed and found to be weak [4], the effect of EPs on turbulence has not been much studied so far (see, e.g., Ref. [5]) and represents the aim of this Letter.This analysis is done via the excitation by EPs of a class of modes naturally existing in tokamaks: the geodesic acoustic modes (GAMs) [6], which are the oscillatory component of large scale E Â B zonal flows.The EP-driven GAMs are called energetic geodesic acoustic modes (EGAMs).These modes have been predicted theoretically [7,8], observed experimentally [9,10], and reported very recently numerically in the absence of turbulence [11] in gyrokinetic simulations with the 5D GYSELA code [12].The motivation of the present Letter relies upon fluid simulations where the turbulence level was controlled by GAMs in the transitional regime between core and edge [13].In addition, experimental evidence of the role of GAMs in the edge-turbulence suppression has been reported for the first time during the analysis of the L-H transition in the ASDEX Upgrade tokamak [14].However, in the context of coreturbulence suppression, the role of GAMs is less evident for several reasons.First, these modes are Landau damped in the core plasma.Second, since they are nonlinearly generated by turbulence, their external control has proven difficult.Last, their frequency !GAM is close to the characteristic turbulence frequency !turb , which means that the shearing rate provided by GAMs might be large compared to the autocorrelation time.In that respect, theoretical predictions of the shearing effect provided by GAMs are not straightforward.The first two difficulties are overcome in this Letter by the presence of EPs providing a continuous excitation of the mode.The possibility of efficiently exciting EGAMs in (global) flux-driven gyrokinetic simulations opens the way to the analysis of turbulence regulation by external ways as well as the impact of EPs on turbulence, and sheds light on the third above-mentioned difficulty.In this Letter, we provide first numerical evidence of the excitation of EGAMs in a fully developed ion temperature gradient (ITG) turbulence as well as the impact of these EGAMs on the turbulent transport.This is done by a source specifically designed to increase the population of EPs.Three main results are presented: (1) EGAMs can be excited in the presence of ITG turbulence in a steady-state regime, (2) the resulting oscillating radial electric shear is not able to suppress the core turbulence, on the contrary (3) a complex interaction between EGAMs and turbulence is evidenced, leading to a radial propagation of turbulence and the destruction of an existing transport barrier.These results are illustrated in Fig. 1, where we display the dependence on time and minor radius of the turbulent diffusivity governed by fluctuations of the radial E Â B drift velocity v Er .Three phases are clearly visible, and will be detailed in the remainder of this Letter: (A) the energetic particle source is applied to a statistical steady-state turbulence regime, (B) a transport barrier is triggered, and (C) EGAMs and turbulence coexist and interact with each other.
GYSELA solves the standard gyrokinetic equation in the electrostatic limit for the ion guiding-center distribution function in a simplified magnetic topology consisting of toroidal flux surfaces with circular poloidal cross sections.In this model, the electrons are considered adiabatic.The self-consistent system is composed by the gyrokinetic equation coupled to the quasineutrality condition, providing the ion guiding-center distribution function F and the electrostatic potential : À hi À 1 n eq r ?: n eq B 2 r ?¼ n i À n eq n eq : (2) In Eqs. ( 1) and (2) and hereafter, all quantities are normalized.Three normalizing (constant) quantities only are needed, namely density n 0 , temperature T 0 , and the magnetic field on the magnetic axis B 0 .In this framework, temperatures are normalized to T 0 , velocities are normalized to the thermal velocity v T0 ¼ ðT Here, m i is the ion mass and e the elementary charge.In Eq. ( 2), B is the normalized modulus of the dimensionless magnetic field B ¼ ðR 0 =RÞb, where b is the unit vector parallel to the magnetic field and R 0 the major radius on the magnetic axis.R is the major radius: R ¼ R 0 þ r cos, r being the minor radius and the poloidal angle.The electrostatic potential is normalized to the thermal energy ( standing for e=T 0 ).The subscript eq stands for equilibrium quantities, the brackets hÁi represent a flux-surface average, and n i is the normalized ion guiding-center density.
The right-hand side of Eq. ( 1) accounts for collisions CðFÞ [15] and heat sources S. The source term is essential for the study of neoclassical and turbulent transports on energy confinement times, away from the initial state.The EGAM instability belongs to the more general class of bump-on-tail instabilities, where the drive comes from a positive slope of the distribution function in velocity space.Therefore, it is required to invert the negative slope of the initial Maxwellian distribution function so as to obtain an out-of-equilibrium state.The system will then tend to relax towards the thermodynamical equilibrium unless it is forced by additional external sources aiming at inverting continuously the population.To this end, the source S is made of two components: S ¼ S th þ S EP , with S th a source of thermal energy [16] and S EP an EP heating, hence, leading to a positive slope in energy.S th is localized in the inner region while S EP is spread across the mid radial position r mid ¼ ðr min þ r max Þ=2.The source S EP is built by using projections onto the Laguerre and Hermite polynomial basis, which enables a separation between injection of energy, parallel momentum, and vorticity with a convenient choice of real coefficients [17].The coefficients are chosen so as to inject only parallel energy.For symmetry reasons, the source is written as where S EP;0 is the amplitude of the source and S r its radial envelope ( R rdrS r ¼ 1).stands for =T s? , with the adiabatic invariant and T s? the normalized transverse temperature of the source, set to 1 in the present simulations.S AE includes the decomposition onto the specified polynomial basis where H h and L l are the Hermite and Laguerre polynomials of degree h and l, respectively.Here, is an arbitrary normalized velocity v 0 divided by the normalized parallel temperature of the source T sk , and . Both v 0 and T sk are critical parameters in view of exciting EGAMs (see below).This source mimics the effects of two tangential neutral beam injectors, oriented in the co-and countercurrent directions so that no external momentum injection exists.
Two simulations are presented, characterized by the following dimensionless parameters: the collisionality is ?¼ 0:02 (banana regime), the profile of the safety factor q is parabolic with relatively low magnetic shear (0 < rd r logq < 0:4) and such that qðr mid Þ % 2:7.The normalized Larmor radius is ?¼ 0 =a ¼ 1=150 (with a the minor radius of the poloidal cross section).The initial density and temperature profiles are characterized by R 0 =L n ¼ 2:2 and R 0 =L T ¼ 6:5, respectively, with L n;T ¼ Àðd r logfn; TgÞ À1 the density and temperature gradient lengths.L n remains constant during the simulations.Both simulations are identical in the time window 0 < t < t init ¼ 2250, where only the source S th is used, with an input power of 4 MW.At t init , turbulence has reached a statistically steady state.The source S EP is then switched on, leading to an additional input power of 2 MW (beginning of phase A).In the following, the two simulations are referred to as simulation with EPs (v 0 ¼ 2, T sk ¼ 0:5) and simulation without EPs (v 0 ¼ 0, T sk ¼ 1).In both simulations, the parallel and transverse injected energies are the same.
In the simulation with EPs, the source S EP modifies the distribution function by depleting the population of particles around v k % AEv 0 and creating two bumps on the tail of the distribution function.This modification leads to a reduction of ITG activity, which is explained as follows.Considering only resonant modes k k ¼ 0, the resonance for ITG modes is !À n dT E ¼ 0, where n is the toroidal wave number, E the normalized kinetic energy, and dT ¼ q rR 0 d the normalized precession frequency, with d a dimensionless frequency [18].Following this simplified framework, the resonance occurs at E ¼ E ! !=ðn dT Þ, which is estimated by considering that turbulence develops around the diamagnetic frequency !? %k T=BL n [16], such that E ! % ðT= d ÞR 0 =L n for resonant modes.Since the chosen value of R 0 =L n is close to v 2 0 =2, the source S EP leads to a reduction of the number of resonant particles for the ITG instability in this simulation.In the simulation with EPs, the depletion around E ! occurs mainly in the outer region ( > 0:5).This is illustrated in Fig. 2(a), where we plot the time evolution of the distribution function at v k ¼ v 0 normalized to its value at t ¼ t init .The observed reduction of turbulence for > 0:5 (phase B of Fig. 1) then appears consistent with the linear stabilization of ITG modes due to this depletion mechanism.
In Fig. 2(b) we plot the time evolution of the slope of the distribution function at the resonant velocity v res ¼ q! EGAM R, with !EGAM the frequency of the EP mode [7,8,11].After the introduction of S EP , the slope is clearly inverted, mainly in the region > 0:5.It is the positiveness of @ v k F eq at the resonant velocity v res in the simulation with EPs which drives EGAMs.This is evidenced by analyzing the contribution of the up-down poloidally asymmetric components, i.e., sin, to the axisymmetric modes, i.e., n ¼ 0, of the electrostatic potential.The time Fourier transform of the imaginary part of the mode 1;0 for t > t init is shown in Fig. 3 for both simulations.We observe a peak at ! EGAM ¼ 0:4! GAM (and also the second harmonic at 2! EGAM ) which clearly dominates the spectrum in the simulation with EPs (solid red line).In the simulation without EPs (dotted blue line), the spectrum does not exhibit this peak.The imaginary part normalized to the modulus j 1;0 j is also plotted (shadowed region) for the simulation with EPs.This curve clearly shows that the imaginary part dominates over the real part.Finally, analysis of the non axisymmetric modes proves that the EGAM frequency is embedded in the turbulence spectrum (not shown here).Therefore, analytical predictions of the effect of EGAMs on turbulence is not straightforward and the gyrokinetic simulations presented in this Letter prove especially useful.
An estimate of the E Â B diffusivity is given by  where Q EÂB ¼ hv Er pi is the radial heat flux associated to v Er and p is the pressure.The time evolution of EÂB is plotted in Fig. 4 for both simulations with (solid and dotted black lines) and without (dashed blue line) EPs in the region 0:5 < < 0:8.In the simulation with EPs, the expression ( 5) is decomposed into the contributions of axisymmetric and nonaxisymmetric modes.Three observations are made from this figure: (1) in the absence of EPs, the modification of turbulent transport remains within the temporal fluctuations and is therefore not noticeable, (2) the contribution of axisymmetric modes in the simulation with EPs exhibits large amplitude oscillations around zero at the EGAM frequency, leading to a vanishing contribution when averaged over some EGAM cycles, and (3) the E Â B diffusivity, mainly due to nonaxisymmetric modes, increases during the EGAM activity in phase C and oscillates also at the EGAM frequency, indicating a strong interaction between turbulence and EPs via EGAMs.As shown in Fig. 4, the increase of the E Â B diffusivity is due to the nonaxisymmetric modes and represents, therefore, an increase of the turbulent transport triggered by the excitation of EGAMs.
The complex interaction between turbulence and EGAMs is captured in Fig. 5, where the evolution of the oscillating part of R 0 =L T is plotted at two different time windows.The oscillating part of R 0 =L T corresponds to R 0 =L T À hR 0 =L T i t , where hÁi t is a time average.This figure corresponds to the first instants of the EGAM excitation at the end of the phase B (lower panel) and to the phase where turbulence and EGAMs coexist, during phase C (top panel).During the transport barrier, the inner region exhibits avalanchelike behavior, as observed in Fig. 5 (lower panel), with fronts propagating outward and vanishing at % 0:5.The propagation velocity is of the order of v aval % 0:8v ?, with v ? the ion diamagnetic velocity.In the same figure, we observe static oscillations in the outer region, characterized by horizontal traces.By static we mean that there is no front propagation; i.e., the beginning of the EGAM oscillations does not exhibit avalanchelike behavior for > 0:5.However, in phase C (top panel of Fig. 5), we see that there is not a single propagation velocity; i.e., both outward propagating fronts and static oscillations coexist and the outer radial region is also characterized by an avalanchelike behavior.This means that the inner avalanches propagate also in the outer region in the presence of EGAMs and the energy flows along the whole radial dimension by means of radially elongated structures.This result indicates, in particular, that a radial electric shear oscillating at a frequency !$ !turb is likely not an efficient way to suppress turbulence.
In conclusion, we have presented, for the first time, evidence of the modification of turbulence by energetic particles.This has been done in flux-driven gyrokinetic simulations, which represent a powerful tool in view of elucidating the complex interaction between turbulence and energetic particles via the excitation of EGAMs.We have shown that the presence of energetic particles in a transport barrier efficiently excites EGAMs, leading to an oscillating radial electric shear that does not suppress the turbulence.On the contrary, the effect on turbulent transport leads to the destruction of the preexisting transport barrier.The complex interaction has been evidenced by the coupling between turbulence-driven avalanches and EGAM-driven static oscillations, resulting in a regime where turbulent transport (modulated at the EGAM frequency) and EGAMs coexist in steady state.
This work was granted access to the HPC resources of CINES under the allocation 2012052224 made by GENCI (Grand Equipement National de Calcul Intensif).This work was supported by EURATOM and carried out within the framework of the European Fusion Development Agreement.The views and opinions expressed herein do not necessarily reflect those of the European Commission.

χ
FIG. 1 (color online).Color map of the E Â B diffusivity ( EÂB ) with the three phases analyzed in this Letter.

FIG. 2 (
FIG. 2 (color online).Time evolution of the distribution function at v k ¼ v 0 (top panel) and the derivative in velocity space at the EGAM resonant velocity (lower panel).

FIG. 3 (
FIG. 3 (color online).Time Fourier transform of the sin component of 1;0 for both simulations with (solid red line) and without (dotted blue line) EPs.For the simulation with EPs, the imaginary part of 1;0 =j 1;0 j is plotted (shadowed region).

FIG. 4 (
FIG. 4 (color online).Time evolution of the E Â B diffusivity for both simulations.For the simulation with EPs, the respective contributions to EÂB of axisymmetric modes (dotted black line) and of nonaxisymmetric modes (solid black line) are distinguished.For the simulation without EPs, no distinction is made (dashed blue line).

FIG. 5 (
FIG. 5 (color online).Oscillating part of R 0 =L T in two time intervals of the simulation with EPs: end of phase B (bottom) and during phase C (top).