Gyrokinetic investigation of Alfv\'en instabilities in the presence of turbulence

The nonlinear dynamics of beta-induced Alfv\'en Eigenmodes (BAE) driven by energetic particles (EP) in the presence of ion-temperature-gradient (ITG) turbulence is investigated, by means of selfconsistent global gyrokinetic simulations and analytical theory. A tokamak magnetic equilibrium with large aspect ratio and reversed shear is considered. A previous study of this configuration has shown that the electron species plays an important role in determining the nonlinear saturation level of a BAE in the absence of turbulence [A. Biancalani, et al., J. Plasma Phys. (2020)]. Here, we extend the study to a turbulent plasma. The EPs are found modify the heat fluxes by introducing energy at the large spatial scales, mainly at the toroidal mode number of the dominant BAE and its harmonics. In this regime, BAEs are found to carry a strong electron heat flux. The feed-back of the global relaxation of the temperature profiles induced by the BAE, and on the turbulence dynamics, is also discussed.


Introduction
Magnetically confined plasmas are complex systems in which waves and instabilities at multiple spatial scales coexist and influence each other. Important examples are microinstabilities, mesoscale zonal flows and macroscopic MHD instabilities like Alfvén modes. Microinstabilities, like ion-temperature-gradient (ITG) modes [1], are linearly unstable due to the gradients of plasma temperature, and nonlinearly interact and saturate forming turbulence states. ITGs carry particle and heat fluxes in the direction of the nonuniformity, and therefore they are particularly deleterious to the heat and particle confinement. One of the products of the nonlinear interaction of microinstabilities is the formation of zonal flows (ZF). ZFs are ExB flows (primarily in the poloidal direction) associated with purely radial variations of the electrostatic potential. They can play the role of the dominant turbulence saturation mechanism [2], by breaking the turbulence vortices, and consequently pushing the energy towards higher radial wavenumbers, where the plasma absorption occurs at kinetic scales. Alfvén modes (AM) are eigenmodes of the shear Alfvén waves, i.e. electromagnetic plasma waves propagating as transverse waves along the magnetic field lines. Various types of AMs exist, such as Global Alfvén Eigenmodes (GAE) [3], Toroidicity-induced Alfvén Eigenmodes (TAE) [4] or Beta-induced Alfvén Eigenmodes (BAE) [5,6,7]. AMs can be driven unstable by suprathermal ions, named here energetic particles (EP), present in tokamak plasmas due to external heating mechanisms and to nuclear fusion reactions. AMs can then lead to a redistribution of the EP population [8], which can have consequences, inter alia, on plasma heating.
The difference in temporal and spatial scales has been invoked in the past to justify a separate treatment of AMs and microturbulence. Nevertheless, they can mutually interact either due to ITG modes [43,44]. In this paper, we discuss the dynamics of AMs in the presence of ITG turbulence, as shown by ORB5 global self-consistent simulations [45], and compare with the estimation of analytical theory.
The structure of the paper is the following. The description of the numerical model of ORB5, of the numerical experiment, and of the considered tokamak case, are given in Sec. 2 and Sec. 3. The heat flux of ITGs and BAEs separately, is described respectively in Sec. 4 and Sec. 5. The results of the numerical simulations of AMs and turbulence are shown in Sec. 6, with dedicated subsections on the evolution of the zonal and nonzonal electric fields, on the evolution of the heat fluxes, on the evolution of the temperature profiles, and on the consequent modification of the spectra at low-n, dominated by the BAEs and at high-n, dominated by the ITGs. The difference between the heat fluxes carried by ITGs and BAEs is also evaluated analytically, to support the numerical findings, and is described in Sec. 7. Finally, a summary of the results and a discussion are given in Sec. 8.

The model
The investigation of the dynamics of AMs, EPs, and turbulence, requires to treat all 3 species (thermal ions, thermal electrons, and energetic ions) with a kinetic model. Among the main reasons, we mention here the importance of the wave-particle resonances which are essential for determining the ion and electron Landau damping, and for the AM drive due to the EPs. Due to the low frequencies of the modes of interest, with respect to the ion cyclotron frequency, we are allowed to reduce the complexity of the kinetic model from 6D in phase space to 5D, by averaging out the fast cyclotron motion of the particles around the magnetic field lines: this reduced model is called gyrokinetics (GK). Global simulations allow to investigate the intrinsically multiscale dynamics of a system made of microturbulence, meso-scale zonal flows, and macro-scale AMs. By means of global simulations, we can investigate the interplay of the several instabilities at the different positions where they develop.
In this paper, we use the GK framework for solving numerical simulations, and for making comparisons with analytical theory, which helps interpreting the results of the simulations. The main numerical tool used for the investigations described here is the global GK particle-in-cell code ORB5. The model equations of ORB5 are the gyrocenter trajectories, and the two equations for the fields.
The gyrocenter trajectories are: The phase-space coordinates are Z = (R, p , µ), i.e. respectively the gyrocenter position, canonical parallel moment p = mU + (e/c) A G and magnetic moment µ = mv 2 ⊥ /(2B). The timedependent fields are named φ and A , and they are respectively the perturbed scalar potential and the parallel component of the perturbed vector potential. In our notation, on the other hand, A is the equilibrium vector potential. The Jacobian is given by the parallel component of B * = B + (c/e)p ∇ × b, where B and b are the equilibrium magnetic field and magnetic unit vector. The summation is over all species present in the plasma. The gyroaverage operator is defined by: where α here is the gyroangle and ρ L = ρ L (α, µ) is the Larmor radius. The gyroaverage operator reduces to the zeroth Bessel function J 0 (k ⊥ ρ Li ) if we transform into Fourier space. The gyroaverage is calculated for all ion species. For electrons, ρ L → 0, therefore φ G = φ(R) (see [37] for more detail). In other words, we take into account finite Larmor radius of ions, and we neglect them for electrons. The quasineutrality equation is: Here f and f M are the total and equilibrium (i.e. independent of time) distribution functions, the integrals are over the phase space volume, with dV being the real-space infinitesimal and dW = (2π/m 2 )B * dp dµ the velocity-space infinitesimal. Finally, the Ampère equation is: Electromagnetic particle-in-cell gyrokinetic codes are affected by a numerical problem called the "cancellation problem", which arises from the fact that some terms in the Ampère equation are solved analytically, and some terms are evaluated by means of a marker discretisation, which introduces a statistical error. A detailed description of the cancellation problem and possible mitigation techniques can be found in [46]. In particular, the "pull-back" mitigation technique used for the simulations discussed in this paper, is described in [47].
This study of AMs in the presence of turbulence is performed by running a first simulation with turbulence only, and switching on the EP effects in a second simulation, namely a restart. This is done in the following way: -firstly, a simulation with 3 species is initialized, where the EP density is n EP /n e = 0.01, and the EP profiles are identical to those of the thermal species; -secondly, a "restart" simulation is performed, i.e. starting from the last time step of the previous one, and with modified EP profiles taken from Ref. [34]. In particular, we consider a flat EP temperature profile, with the EP temperature 10 times higher than the thermal species at the reference position at mid-radius. The EP density gradient is 10 times higher than the thermal gradients at the reference position. Moreover, EPs are allowed to redistribute in phase space, i.e. they are allowed to follow perturbed orbits, like the thermal species. This defines fully nonlinear simulations.
The quasineutraility is imposed by ORB5 at every restart, by modifying the electron profile in order to have n e = n i + n EP at each radial position. As a consequence, due to the larger density gradient of the EPs species after the "switch", the electron density is slightly steeper. The effect of this "switch" on the linear dynamics of the ITGs has been found to be negligible.
Perturbed modes with toroidal mode number in the range 20 ≤ n ≤ 30 are initialized, while all toroidal modes in the range 0 ≤ n ≤ 40 are simulated. For each toroidal mode n, a radial dependent filter is applied, allowing for poloidal mode numbers in the range m = n q(s) ± δm. A spatial grid of (ns, nchi, nphi) = (256, 384, 192) points, a time step of dt=5 Ω −1 i (Ω i being the ion cyclotron frequency), and a number of markers of (n i , n e , n EP ) = (2, 10, 2)e8 respectively for the thermal ions, electrons and EP are used for the turbulence simulations. Unicity boundary conditions are applied at potentials at the axis and Dirichlet at the edge. No collisions are considered here. We consider a Krook operator which acts like a source for the thermal species, by restoring the initial plasma profiles, and at the same imposes an artificial damping to all modes, with value γ K = 1.0 · 10 −4 Ω i .

Magnetic equilibrium and plasma profiles
The magnetic equilibrium and plasma profiles used here are the same as in Ref. [34] for the study of AMs, and in Ref. [45] for the study of AMs and turbulence. The major and minor radii are R 0 = 1.0 m, and a = 0.1 m, and the toroidal magnetic field at the axis is B 0 = 3.0 T. An equilibrium with circular concentric flux surfaces is considered. The safety factor has a value of 1.85 at the axis, it decreases from ρ=0 to ρ=0.5, where the minimum value is located  (q(ρ = 0.5)=1.78), and then it raises to the edge, where it reaches the maximum value (q(ρ = 1)=2.6) (see Fig. 1). Here ρ is a normalized radial coordinate defined as ρ = r/a. A reference radial position is chosen at ρ = ρ r = 0.5, corresponding to s = 0.525, where the flux radial coordinate s is defined as s = ψ pol /ψ pol (edge), ψ pol being the poloidal magnetic flux. The ion and electron temperature profiles are the same: T e (ρ) = T i (ρ). The temperature at the reference radius is chosen to have a value of T e (ρ = ρ r ) corresponding ρ * = ρ s /a = 0.00571 (with ρ s = T e /m i /Ω i being the sound Larmor radius). The electron thermal to magnetic pressure ratio of β e = 8π n e T e (ρ r )/B 2 0 = 5 · 10 −4 . The ion species (thermal and energetic) have the mass of the deuterium. The value of the electron mass is chosen as m i /m e = 200. This value is chosen in order to have quicker numerical simulations. This value is found to be at convergence for the saturation levels of the BAE, and for the linear dynamics of the ITG. The dynamics of the BAE at very low EP concentration depends on the value of the electron mass due to the electron Landau damping, and in particular the BAE in the absence of EPs is found to be slightly above marginal stability, with a growth rate which is small but finite. In this paper, we show the results of simulation where a krook operator is used, which adds an artificial damping to all modes, with a value of γ K = 1.0 · 10 −4 Ω i . This value is sufficiently high to make BAEs linearly stable in the absence of EPs (independently on the electron mass).
The equilibrium considered has some simplifications with respect to present-day tokamaks: namely a higher aspect ratio, and circular concentric flux surfaces. This allows to capture the physics of a simple case, before going to more accurate experimental cases. At the same time, an effort of comparing the results of the code ORB5 with experimental measurements of EP-driven modes is being done in separate papers ( [48,49]).

ITG linear dynamics
The radial location is near the s = 0.525 reference surface, where the peak of the temperature gradient is. The EPs are found not to be sensibly changing the ITG dynamics for concentration lower than n EP /n e = 0.05, which is the regime of interest in this paper.
Here we show the spectrum of linear electromagnetic simulations of ITG, without sources/sinks, i.e. without krook operator. The most unstable mode is around n = 30, with growth rate: In order to characterise the ITG heat flux, a linear electromagnetic simulation of a single ITG, with n = 26, is considered. The EP concentration is n EP /n e = 0.01 The heat flux is measured, and it is found to grow exponentially. The growth rate of the heat flux is measured to be γ h = 0.83e − 3.
We write the volume averaged heat flux as: where s is the species index, and n is the toroidal mode number. We measure for this ITG with n = 26 the following values respectively for the thermal ions, electrons, and EPs: α IT G,e = 0.014 (9) α IT G,EP = 0.002 (10)  where α is in units of n s T s c s ρ * T 2 e /e 2 , with n s being the density, and T s being the temperature measured in eV.
The ratios of the heat fluxes and the electron heat flux, which can be calculated as the coefficients of the thermal and energetic ion heat transport normalized to the coefficient of the electrons, are: Although the values of these ratios for the thermal ions and EPs slightly depend on the electron mass, nevertheless they are found to always be respectively higher and lower than unity. This characterises the ITG in the regime of interest.

Dynamics of AMs without turbulence: heat flux
As an example, we consider here a fully nonlinear simulation with only n = 5. No other mode is allowed to develop (therefore there is no turbulence). The EP concentration is n EP /n e = 0.01. The mode is found to grow linearly, saturate at around t = 1.3e4 Ω −1 i , and damp due to the EP radial redistribution (see Ref. [34] for a detailed analysis of the BAE dynamics in this equilibrium). The evolution in time of the scalar potential and of the heat flux are shown respectively in Fig We write the volume averaged heat flux as in Eq. 7. In the linear phase, we measure for this  This is a simulation with only n = 5, fully nonlinear, and with n EP /n e = 0.01, representing a BAE.
BAE with n = 5 the following values respectively for the thermal ions, electrons, and EPs: α BAE,i = 0.023 (13) α BAE,e = 0.035 (14) α BAE,EP = 0.083 (15) Note that, with respect to the ITG (see Sec. 4), a BAE drives a lower ion heat flux, a similar electron heat flux, and a much higher EP heat flux (1 order of magnitude higher). The ratios of the heat fluxes and the electron heat flux, which can be calculated as the coefficients of the thermal and energetic ion heat transport normalized to the coefficient of the electrons, are: These characteristic values are found to be stable during the linear phase, and they change after the nonlinear saturation (see Fig. 5-bottom). They are found to be at convergence with the electron mass for m i /m e = 200. Note that the EP species carries a heat flux of the same order of magnitude of the thermal species, in these simulations. In order to understand this, we can calculate the value of the ion diamagnetic frequency, for a mode with n=5 sitting at s=0. 4: The diamagnetic frequency is about two orders of magnitude lower than the mode frequency [34]. So, for our case, we have that an energetic species with the same concentration of the thermal species, is expected to carry a radial flux of about 2 orders of magnitude higher than the thermal species, due to the higher temperature (1 order of magnitude) and higher density gradient (1 order of magnitude). The EP species considered here has a concentration of 2 orders of magnitude lower than the thermal species, which explains the fact that the heat flux is of the same order of magnitude of the thermal heat flux. The reason why a BAE can carry a significant amount of heat flux of the thermal species, with respect to other AMs, is due to its relatively low frequency, which increases the importance of the wave-particle resonances [7]. Regarding the thermal ions, the dominant resonance is expected to be with the transit frequency. Regarding the electrons, which typically have a much higher thermal velocity, the resonances are expected to occur mainly with barely trapped electrons, similarly to what happens for EGAMs [50]. The strong electron transport due to AMs found here with ORB5 simulations, is consistent with earlier theoretical predictions [51] (see also Ref. [52] for a treatment in non-toroidal magnetic equilibria) and experimental observations [53].
6 Nonlinear dynamics of AMs with turbulence 6

.1 Evolution of the fields
Here, the results of the nonlinear simulations with turbulence and AMs are shown. The restart with EPs switched on is performed at t ≃ 4.9e4 Ω −1 i . Before discussing the results of the simulations where the EPs have their nominal temperature, i.e. T EP (s r )/T e (s r ) = 10, we have run a simulation where the EP switch is performed only by increasing the density profile to κ n = 10, but keeping the EP temperature equal to the thermal species. This test is done to study the effect of the modification in the density profiles alone. No sensible change is observed. Therefore, we can state that the effect of density gradient alone of the minority with thermal temperature, is negligible. We can now study the effect of the EPs with nominal temperature i.e. T EP (s r )/T e (s r ) = 10. In Fig. 6-a, the evolution in time of the maximum of the radial electric field is shown, for a simulation where the EPs are switched on at the restart, and compared with a simulation where the EP are not switched on at the restart. One can see that the effect of the EPs is to excite both the nonzonal and zonal components of the electric field, which start growing, then saturate around t ≃ 6e4 Ω −1 i , and then damp in time. The radial structure shows that the BAE mode with n=5, m=9 is excited by the EPs in the second part of the simulation, when the EPs are switched on (see Fig. 6-b). This BAE is not clearly visible before the EPs are switched on. This mode is identified as a BAE due to the polarization, i.e. a clear n=5, m=9 signature, and the frequency, which is slightly higher than the linear BAE frequency, i.e. two orders of magnitude higher than the frequency of the ITG with n=5 (see Sec. 4). This frequency is observed at each radius in the domain where the linear BAE and the linear ITG are observed. Therefore, this BAE is shown to be dominant in amplitude on the ITG modes in this time window (see Fig. 6-a). The saturation level of the BAE for this value of krook is E r ≃ 1.0 · 10 5 V /m (lowering the value of the krook corresponds to decreasing the artificial damping and reaching a higher saturation level). This value is found to be the same as in the absence of turbulence (see Ref. [34]).

Evolution of the total heat fluxes
The volume averaged radial heat flux can be studied in a simulation with EPs (see Fig. 7-left). At the moment of the BAE saturation, very similar levels of the heat fluxes of the turbulence simulations are obtained with the simulations which retains only 0 ≤ n ≤ 10. If we filter out the n = 10 mode, slightly lower values of the thermal heat fluxes are obtained. This means that these heat fluxes are not due to the high-n ITG modes, but to the n = 5 and n = 10 modes. The ratios of the heat fluxes with the electron heat flux, are given in Fig. 7-right. We note that the heat fluxes go from an ITG dominated regime, to a BAE dominated regime, and then back to ITG. In particular, in the first time range when the EPs are not yet switched on, we have Γ i /Γ e bigger than 1 (Γ i /Γ e >≃ 3), which goes down to values lower than unity when the EPs are switched on and the heat flux is dominated by the BAE (Γ i /Γ e >≃ 0.7), and then back to values higher than unity when the BAE is damped. Viceversa, Γ EP /Γ e > starts with values higher than unity when they are switched on, indicating a dominant Alfvénic activity, to decrease to vaues lower than unity after the BAE has faded away. When sending simulations with 0 ≤ n ≤ 10 only, we also try to switch the Krook on, like in the simulation with turbulence, but the values are not sensibly changed. If we consider the simulation where also the mode with n = 10 is filtered out, the values are not sensibly changed.

Evolution of the temperature profiles
In this subsection, we show the temperature profile evolution in time in the simulations with and without EPs. The simulation is the same as in the previous section. Sources are applied to the thermal species at the restart in the form of a Krook operator, like in the previous sections. In Fig. 8, the difference of the temperature at a certain time and the temperature at the beginning of the simulation is shown. The original profiles are normalized to their value at the reference radius s r = 0.525. The time is chosen around the peak of the BAE amplitude. Around this time, the temperature profile is found to have small oscillations and small radial structures are also formed.
Note that the ITG turbulence redistributes the temperature of about a factor 2 higher for the thermal ions than for the electrons, whereas the BAE redistributes both ions and electrons of about the same amount. Note also that the amplitude of the ion temperature perturbation of the BAE is about 3 times higher than the ITG.   In Fig. 9, the evolution of the temperature profiles is shown. Note the flattening of the profiles around the location of the BAE, due to the heat flux carried by the BAE.

Spectra in toroidal mode number
Both ITG turbulence and AMs are known to induce cross-field heat transport. In the simulations presented here, both of these effects are co-existing, raising two general questions: What is the relative importance of each of these two mechanisms, and do they interfere in any way?
The toroidal mode number spectrum of the ion heat flux (for the details of the diagnostic, see Ref. [54]) for a simulation where the EPs are switched on can be seen as a the continuous red line in Fig. 10. For comparison, the same spectrum for a simulation without EPs is shown as a dashed red line. A sensible general trend that we find is that the heat fluxes are higher in the simulation with EPs. This is explained by the fact that the EP density profile consitutes an additional source of free energy to the system.

Analysis of the spectrum at low-n
We observe that the nonzonal spectra given in Fig. 10 are dominated by the modes n = 5, n = 10 and n = 15, namely the first, second and third harmonics of the main BAE mode. The contributions of the different poloidal mode numbers to the heat fluxes for the dominant mode, namely the mode with n = 10, is depicted for a simulation without and with EPs, in Fig. 11. In the case of the simulation without EPs, a broad spectrum of poloidal components shows the characteristic polarization of the ITG mode. On the other hand, for the case of the simulation with EPs, the heat flux is dominated by the poloidal mode m = 18. Note also that the peak of the m = 18 component is at the position of the peak of the BAE field. These two are signatures of the second harmonics of the BAE. Therefore, we can state that the ITG dynamics is subdominant in the simulation with EPs, with respect to the BAE, in carrying the heat transport. It is also important to note that the other poloidal components (i.e., m = 18) are increased in the simulation with EPs, in comparison with the simulation without EPs. This means that the BAE second harmonic is efficient in modifying the dynamics of the ITG of n = 10, due to the nonlinear interaction. This is an example of cross-scale interaction, with the AM being the macro-scale mode, and the ITG turbulence being constituted mainly by micro-instabilities.

Analysis of the spectrum at high-n
It is worthy to note that not only the low-n part of the heat flux spectrum, i.e., where the BAEs are dominant, has higher levels in the presence of EPs, but also the higher-n part, i.e., where the ITGs are dominant. This is of interest for the question of how the presence of EPs modifies the ITG turbulence. In this regime, we can state that an EP population driving a BAE linearly unstable, affects the turbulence dynamics by increasing the heat transport. The fluctuation amplitude of the scalar potential φ and density ρ for the different toroidal mode numbers is also shown in Fig. 10 (respectively with blue lines and blue black lines). Differently from the heat flux, the fluctuation amplitude given by the scalar potential is shown not to be sensibly modified by the presence of EPs in the range of toroidal mode numbers of the ITG. Moreover, the fluctuation amplitude given by the perturbed density is found to be decreased by the presence of EPs.
We can also investigate the relative importance of the ion and electron heat fluxes in the domain of high-n. This is shown in Fig. 12 and 13. We note that the ion heat flux is always bigger than the electron heat flux in the radial region 0.5 < s < 0.7, with or without EPs, because the spectrum is dominated by the ITG turbulence. On the other hand, the electron heat flux is bigger than the ion heat flux in the radial region 0.3 < s < 0.5, when the EPs are switched on,  Heat fluxes, and fit with n_ i with EPs Fit with _i = -2.5+-0.8 _e with EPs Fit with _e = -2.5+-0.7 Figure 13: Same as in Fig. 12, but with a measurement done far from the location of the BAE, namely in the region 0.5 < s < 0.7. Here the ion heat flux is always bigger than the electron heat flux, as the ITG turbulence dominates the spectra.
because in this region the spectrum becomes dominated by the BAE.

Analytical estimation of the heat fluxes of ballooning and nonballooning modes
In the previous sections, we have described the results of global gyrokinetic simulations of AMs and ITG turbulence. We have shown that AMs drive a strong electron heat flux, in comparison with ITGs, which have the ion heat flux dominant. AMs and ITGs differ mainly because of their frequency, spatial structure, and polarization. In this Section, we want to investigate how the spatial structure of the modes influences the electron heat flux, by estimating the heat flux analytically. To this aim, we consider two types of modes. The first is an ITG: a mode with ballooning structure, i.e. with higher amplitude on the low-field side of the tokamak and lower amplitude at the high-field side of the tokamak. The second is a BAE: a mode with non-ballooning structure, i.e. with equivalent intensity on the low-field side and high-field side.
Let us recall that the distribution function of each species s is written as where r is the particle position, R the guiding-centre position, φ the electrostatic potential, e s the charge, v the speed, λ = v 2 ⊥ /(v 2 B) the ratio between the magnetic moment µ = m s v 2 ⊥ /(2B) and the kinetic energy m s v 2 /2 = x 2 T s , and the Maxwellian with particle density n s (ψ) and temperature T s (ψ). In this notation, the equation for g a becomes where Ω s = e s B/m s is the gyrofrequency and the derivatives are taken at fixed energy and magnetic moment. J 0 is the zeroth-order Bessel function of the first kind, and corresponds to a gyroaveraging operator in real space. The magnetic field is taken to be B = ∇ψ × ∇α, k ⊥ = k ψ ∇ψ + k α ∇α, where in this section ψ is the toroidal magnetic flux, and α = qθ − ϕ, with θ and ϕ respectively the poloidal and toroidal angles. The diamagnetic and drift frequencies, respectively, are defined by where η s = d ln T s /d ln n s and v ds = Ω −1 s b × v 2 b · ∇b + 0.5v 2 ⊥ ∇B/B denotes the drift velocity. Within this framework, we introduce the flux-surface averaged quasilinear radial heat flux, for species s: The kinetic function g s will be considered in the frequency regime where ω be , and ω tr,e are the bounce and transit frequency. This allows us to solve analytically Eq. (19) for the electrons.
In the following, we will focus on the electron transport, which is found to be dominant for BAEs in the simulations.

Passing electrons
For passing electrons the streaming term is dominant, hence [55] whereψ is so that A = c∇ ψ /(iω). By replacing this expression in Eq. (20), after writing ω = ω r + iγ, one finds In the specific case E → 0, φ =ψ, and for a periodic potential in θ, e.g. φ = cos θ + i sin θ, we have Then, close to marginality, passing particles do not give a contribution to transport.

Trapped electrons
We consider trapped electrons to be mostly electrostatic. We bounce average the electron kinetic equation to obtain where we took J 0 = 1, since we are interested in transport at the ion scale, and the bounce-average is defined in the appendix. By inserting Eq. (25) in (20), one obtains with and We evaluate term by term, using a simple model for the equilibrium field B = B 0 (1 − ǫ cos θ), in the marginal limit γ → 0 + . Results are in the Appendix, and will be used in the following.

Influence of the space structure on the electron heat flux
We consider two types of modes, one with a maximum on the outboard-side and a minimum at the inboard-side (ballooning, ITG-like) one a simple trigonometric function (non-ballooning, BAE-like) Figure 14: Ratio of the coefficients of the electron heat flux due to temperature gradients for ITG-like modes and BAE-like modes. Here ω r ω (0) The result is the estimation of the electron heat flux of ballooning and non-ballooning modes as a dependence on the frequency. As an example, we consider the heat flux given by the temperature gradient term, χ e,T (the others behave similarly). The ratio of the coefficients of the contributions of non-ballooning and ballooning modes is shown in Fig. 14. One can see that, below a certain frequency given by ω r /ω de (0) = 1.5, the electron heat flux is dominated by ballooning modes ("b") like ITGs, whereas above this frequency, the electron heat flux is dominated by non-ballooning modes ("nb") like BAEs. This shows how the spatial structure of BAEs is important, together with the frequency, in driving an important electron heat flux, in comparison with ITGs.

Conclusions and discussion
The transport of heat and particles is one of the most important problems for the confinement of tokamak plasmas in present day tokamaks and future devices. Historically, the radial transport of energetic particles (EPs) by Alfvén modes (AMs) and the heat transport of the thermal species by turbulence, have been treated separately. Nowadays, it is becoming clearer that these two problems can actually be connected, because their mutual influence can be strong in some regimes. Moreover, it is becoming more feasible to study their selfconsistent interaction by means of numerical simulations, due to the more powerful supercomputers, and to more efficient numerical schemes. A kinetic model must be used for this study, due to the importance of the wave-particle resonances with all species (thermal ions, thermal electrons, and energetic ions). The multi-scale nature of the problem demands a global description, because different modes have different spatial size and localization.
In this paper, we have adopted a global gyrokinetic model and applied it to a simplified tokamak equilibrium, where numerical simulations are sufficiently fast to allow the study of the nonlinear dynamics of AMs in the presence of turbulence. The heat flux of AMs has been studied and compared with that of ITGs. It has been found that, in the selected regime, AMs drive an important electron heat flux, in contrast to ITGs which drive a dominant ion heat flux. This has been found to depend not only on the difference of the frequency, but also on the different spatial structure, by means of analytical estimations.
For relatively low concentration of EPs (1%), with 10 times higher temperature than the thermal species, the radial electric field of BAEs has been observed to grow and saturate at levels one order of magnitude higher than that of ITG turbulence. As a consequence, zonal flows driven by the BAEs via forced-driven excitation have also been observed, at levels 10 times higher than those driven by turbulence alone. The BAEs dominate the heat flux spectra at the low toroidal mode numbers and the main harmonics, injecting energy at the large spatial scales.
The global character of the problem has been studied, and we have found that the location of the BAE influences also the dynamics of turbulence in the part of the spectrum of high toroidal mode numbers. In particular, at the radial location of the BAEs, the electron heat flux becomes dominant in the presence of EPs, confirming the importance of the Alfvénic activity, whereas at radial locations far from the BAEs, the ITGs is found to dominate the heat flux even in the presence of EPs, and the ion heat flux remain the dominant one.
In summary, we have shown that an AM like a BAE can drive a heat flux of the thermal species as efficiently and even more efficiently than ITG turbulence, especially for the electrons. This modifies the temperature profiles of the thermal species, flattening them at the location of the AMs. The spectra of the heat fluxes are modified even in the part of high toroidal mode numbers, dominated by the ITGs, with the electron heat flux increasing in amplitude.
As next steps, we will first investigate how changing the location of the AMs influences the turbulence dynamics, in monotonic safety factor profiles and reversed-shear safety-factor profiles. In particular, one effort will be to isolate the effect of the zonal flows. These are known to suppress the turbulence by enhancing the cascade of energy from large spatial scale to small spatial scales, and it will be important to assess how this can enter the multi-scale interaction of AMs and turbulence. Having demonstrated the feasibility of global self-consistent gyrokinetic simulations in simplified equilibria, we will then relax some of the limitations, to approach cases which are closer to experimentally relevant scenarios. As an ultimate goal, comprehensive theoretical studies of burning plasmas can be envisioned, to develop a predictive capability for ITER and future fusion power plants.