First principles of modelling the stabilization of microturbulence by fast ions

The observation that fast ions stabilize ion-temperature-gradient-driven microturbulence has profound implications for future fusion reactors. It is also important in optimizing the performance of present-day devices. In this work, we examine in detail the phenomenology of fast ion stabilization and present a reduced model which describes this effect. This model is derived from the high-energy limit of the gyrokinetic equation and extends the existing"dilution"model to account for nontrivial fast ion kinetics. Our model provides a physically-transparent explanation for the observed stabilization and makes several key qualitative predictions. Firstly, that different classes of fast ions, depending on their radial density or temperature variation, have different stabilizing properties. Secondly, that zonal flows are an important ingredient in this effect precisely because the fast ion zonal response is negligible. Finally, that in the limit of highly-energetic fast ions, their response approaches that of the"dilution"model; in particular, alpha particles are expected to have little, if any, stabilizing effect on plasma turbulence. We support these conclusions through detailed linear and nonlinear gyrokinetic simulations.


Fast ion stabilization of ITG microturbulence
Microturbulence fundamentally limits the confinement time in current and future tokamak experiments [2]. The very gradients that are required to achieve high central densities and temperatures also provide a source of free energy. This free energy drives transport of particles, momentum, and energy, usually far in excess of collisional transport. The ion temper ature gradient (ITG) mode, for instance, limits the core temper ature of tokamaks [3][4][5]. This turbulence occurs on the scale of the thermal ion Larmor radius ρ i .
Gyrokinetics is the reduction of the Fokker-Planck kinetic equation that rigorously handles electromagnetic fields whose fluctuations vary on spatial scales similar to ρ i , but on timescales much slower than the gyro-frequency Ω i [6][7][8]. A multitude of computational tools have been developed to solve the nonlinear gyrokinetic equation [9][10][11][12]. These tools take as inputs the equilibrium magnetic field, density, and temperature profiles, and predict the ensuing turbulent fluctuations and the associated transport fluxes. By using experimentally-determined profiles and comparing the computed fluxes to experimentally-inferred fluxes (usually from powerbalance calculations), and through more detailed comparisons of turbulence characteristics, the validity of the gyrokinetic approach can be verified [13]. This has been widely demonstrated [14][15][16][17].
However, due to the scarcity of computational resources, these matching exercises of necessity entail the use of simulations that neglect various parts of the complex physics of the experimental setup. It was discovered during one of these exercises, in which experiments on the Joint European Torus [18] were analysed, that effects from non-thermal minority ions had to be accounted for. Indeed, the predictions of the gyrokinetic codes in the absence of such ions suggested far greater transport than was observed [19]. This effect-the fast ion stabilization of ITG turbulence-is the subject of this paper.
This effect is important because the presence of energetic ions are essential to sustain fusion-relevant bulk temperatures. Two external heating methods are used to supply the plasma with high-energy ions: neutral beam injection (NBI) [20] and ion cyclotron resonance heating (ICRH) [21]. As well as the fast ions resulting from external heating, another class of fast ions, energetic alpha particles, are generated from the fusion reaction itself. In the future, it is anticipated that alpha particles will provide the majority of the heating, and thus, that the plasma will be energetically self-sustaining (i.e. a burning plasma). In the local gyrokinetic simulations used herein, the practical differences between these several classes of fast ions (as well as, of course, their mass and charge) are in their distribution in phase space.
As the main medium of heat injection, fast ions are a critical component of the fusion plasma and thus worthy of study in their own right. However, above and beyond this, it has been observed that a large fast ion population can improve the plasma confinement [22], leading to elevated core temperatures and/or densities. It has been shown in gyrokinetic simulations [19] that this effect is non-trivial. It is important to understand how and why this stabilization occurs, par ticularly when extrapolating to future net-power-producing fusion reactors, where alpha particles provide most of the heating. Although, like externally driven fast ions, alpha particles will contribute significantly to the plasma pressure in reactors, their density is very small (n α 0.01n e ) compared to NBI and ICRH ions. It is therefore not obvious that they will affect turbulence in a similar way.
Owing to the potentially profound implications of this and other effects of fast ions on microturbulence, there has been widespread study of the topic. Follow-up works [23,24] further characterize the stabilization and examine cases both where it is weak and where it is strong. Other research has sought to determine which properties of the fast ions contribute most strongly to their effects; these properties are typically categorized as dilution (of the main ions), non-Maxwellian velocity space distributions, and electromagnetic effects. Tardini et al and Holland et al [25,26] studied the effect that dilution has on suppressing ITG turbulence, with [27] presenting a reduced model for this effect (expanded further in this work). The active response of a hot Maxwellian impurity was investigated analytically [28] and in simulations [29]. The effects of non-Maxwellian fast ions on turbulence were initially studied using isotropic fast ion velocity-space distributions in [30] and was further generalized for anisotropic fast ions in [31]. The electromagnetic effects of fast ions consist a rich topic, primarily focused on the destabilization of Alfvén eigenmodes (AEs) and energetic particle modes (EPMs) [32,33]. The theory of how electromagnetic fluctuations interact with ITG turbulence is continuing to be developed [34][35][36][37][38][39]. In the present work, an analytic model for the effect of fast ions on ITG turbulence (in contrast to modes which are driven by fast ions themselves) is presented.
We seek to explain the dominant mechanism for fast ioninduced stabilization of ITG turbulence from first principles, focusing on a case where the effect was very strong: discharge 73224 of the Joint European Torus (JET) [40]. The results of a comprehensive linear study are presented in section 2, in which it is examined to what extent the stabilization can be characterized as a change in the linear growth rate. Much about the phenomenon will be learned by doing this because it allows a wider variety of high-resolution simulations to be performed than would be unfeasible with nonlinear turbulence simulations. Then, in section 3, the basic problem is simplified further by approximating the fast ion distribution function in the energetic limit, which yields an approximate analytic solution to the gyrokinetic equation. With the most important elements distilled, this simplified model is inserted into Maxwell's equations, in which an effective parameter model, applicable to linear and nonlinear simulations alike, makes itself evident. The electrostatic and electromagnetic effects are modelled respectively as effective modifications to the temper ature ratio τ = T i /T e and β e = 8πn e T e /B 2 , parameters to which microturbulence is known to be sensitive. This model is then benchmarked against nonlinear gyrokinetic simulations and elaborated upon further.

Characterizing the effect of fast ions on the linear ITG mode
We begin by considering the simplified linear gyrokinetic system and ask whether, and how, the presence of fast ions affects the growth of unstable ITG modes that give rise to microturbulence. We accomplish this by analyzing a case in which their effect is particularly strong: JET discharge 73224 [19]. We will find that the presence of fast ions has a nontrivial effect on the ITG mode beyond their global effects on the plasma and their dilution of turbulence driven by the thermal ions 3 . In this section many of the results and analyses presented are a summary of the more extended linear analysis expounded in [41].
When the fluctuations are small enough to treat linearly, ITG may cause the electromagnetic fields to grow exponentially. Thus we find, for example, that the fluctuating electric potential φ ∝ e −iωt , where ω = ω + iγ, ω is the frequency (by convention, positive for waves propagating in the ion diamagnetic direction), and γ (if positive) is the growth rate. In this section, we repeatedly calculate this growth rate because, although eventually the fields grow until the nonlinear interaction between modes competes with the linear physics [42] and the system reaches a saturated turbulent state, the linear growth rate has a direct impact on the strength of the saturated turbulence [43]. Indeed, it is often used for the estimation of saturated field amplitudes in reduced transport models [44,45]. Now in order to study the 'effect of fast ions on the linear growth rate', it is necessary to somewhat artificially isolate particular effects because, without generating fast ions, the discharge would be fundamentally different. In other words, since fast ions are used to heat the plasma, drive the current, and act as a particle source, it is not immediately clear how to compare cases 'with' and 'without' fast ions. In this section, we define this comparison as being between a baseline case that includes fast ions and a hypothetical equivalent case in which fast ions do not participate in the turbulence, but in which global properties of the main ions and electrons are nevertheless the same, except for constraints imposed locally by the absence of fast ion charge.
To be effective at heating, the fast ion pressure ought to be comparable to the thermal pressure. Therefore an additional effect on the magnetic geometry is expected. Our observations indicate that fast ions have little effect on the flux surface shape, but have a non-trivial impact on the safety factor, magnetic shear, and Shafranov shift, which are all known to play a significant role in microturbulence. For more details on the stabilizing role of the fast ion pressure gradient, the reader is directed to [22,24,41]. Henceforth, in order to to isolate the effect of fast ions on the local ITG mode, the equilibrium plasma parameters, including the safety factor, magnetic shear, and flux surface shape, will remain fixed.

Gyrokinetic framework and baseline case
Apart from these global effects, there remains a nontrivial effect of fast ions on the ITG mode as manifest in local gyrokinetics. The simulations presented in this work were performed with the gs2 code [9,46], which solves the gyrokinetic equation [6]: for the perturbed distribution function for several isotropic species s: δf s = Z s eφ ∂F0s ∂E + h s . While we will focus on simulations of isotropic fast ions in this work, we will occasionally discuss the implications of anisotropy. In this case, the gyrokinetic equation is the same, but there is an additional contribution to δf s from the possible μ-dependence of F 0s . The mass and charge of the species are m s and Z s e respectively. The fluctuating fields φ and A (the vector magnetic potential) are represented by a scalar electromagnetic potential χ ≡ φ − v · A/c. The equilibrium distribution of species s is F 0s E, µ, σ , and h s is the non-adiabatic part of the fluctuating distribution function which does not depend on gyrophase ϑ. The equilibrium magnetic field has magnitude B and points in the direction of the unit vector b ≡ B/B. The parallel velocity is defined as v = σ 2 [E − µB (θ)] /m s , with energy E, the exactly-conserved magnetic moment μ, and sign σ = ±1. For a Maxwellian species with a temperature T s , the thermal speed v ts ≡ 2T s /m s , and for a non-Maxwellian species, this represents a characteristic speed defined using an The characteristic Larmor radius is given by ρ s ≡ v ts /Ω s , where Ω s ≡ Z s eB/m s c. The tokamak minor radius a provides an approximate length scale on which the equilibrium and fluctuations parallel to the magnetic field vary, and ρ * ≡ ρ i /a 1. The quantity v D is the magnetic drift velocity of the guiding center and v χ ≡ (c/B) b × ∇ χ R includes the E × B drift along with the drifts and streaming associated with the fluctuating magnetic field. The gyro-average at fixed guiding center R is denoted by φ R = 2π 0 φ (r) dϑ/2π. For simplicity, gradients of equilibrium plasma flows are ignored and the equation is solved in the frame rotating with the plasma. A conservative linearized Fokker-Planck collision operator C [47,48] is employed to model collisions of ions and electrons, and of each species with themselves (collisions between ions of different species are omitted).
Throughout this work, we will focus on JET discharge 73224 around the flux surface with a half-width of r = 0.375a (where a is the half-width of the last closed flux surface). The nominal parameters for this baseline case are based on those from [49] and are listed in tables A1 and A2.
The resolution for the linear simulations of this section are as follows: there are 58 grid points along the field line per poloidal turn (parametrized by poloidal angle θ), and the parallel domain extends to θ = ±22π. The velocity space resolution is 36 grid points in energy and 46 in pitch angle. The time step is 0.03 a/v ti , run until the complex frequency is conv erged to within a factor of 2 × 10 −4 . There are two species of fast ions present in these simulations, each represented as a hightemperature Maxwellian. The local pressure gradient, used to rescale the local geometrical parameters according to the Miller prescription, is calculated consistently depending on the local fast ion pressure gradient.
Unless otherwise stated, all simulations include perpendicular magnetic fluctuations ( A ), with compressive fluctuations (δB ) artificially disabled. This is justified by the relatively low β e , and by the growth rates shown in figure 1. At a critical β e , the kinetic ballooning mode (KBM) becomes unstable and the growth rate significantly increases. Only then is there a discernible difference when including δB . As long as the mode is ITG-like, increasing β e is stabilizing and only the fluctuations of A need to be considered. For a more detailed treatment on the effect of compressive fluctuations on ITG modes, see [38].

Fast ion-induced stabilization of the ITG mode: general observations
The effect of fast ions on the ITG mode growth rate and frequency in this case is shown in figure 2. The baseline case has both NBI fast deuterium and ICRH fast helium-3, with both types together consisting 20% of the positive charge. Then, the growth rates were recalculated for a case with the different types of fast ions removed. A significant increase in the ITG mode growth rate is observed, with a modest change in the frequency spectrum. Finally, both fast ions were removed to generate a case 'without fast ions'. This shows that the effect of NBI is relatively small compared to ICRH. While the growth rate is sensitive to the presence of ICRH fast ions, the ITG mode frequency and eigenfunction are not significantly altered by the presence of fast ions, as was also observed in [49].
A 'mixing-length' estimate for the bulk ion heat flux based on the linear physics is q i ∼ γ/k 2 ⊥ . Using this estimate, our results suggest about a factor of 2 increase in the turbulence amplitude when fast ions are removed. This alone does not account for the strong nonlinear effect in [19] and later in figure 8. These latter results showed about a factor of 10 increase in the thermal ion heat flux. It is in this sense that we say the effect of fast ions on turbulence cannot be explained in the context of linear gyrokinetics.
The decrease in the ITG mode growth rate when fast ions are included as a kinetic species, and the subsequent (disproportionate) decrease in the turbulence amplitude, is the subject of the remainder of this paper. We begin with the simplest explanation for this phenomenon: that only the indirect effects of fast ions are responsible for the stabilization. Namely, we examine cases where the only effects of fast ions considered are their effects on the bulk ion and electron densities.

Dilution
On the temporal and spatial scales of interest, the plasma remains quasineutral: s Z s en s = 0. ( When positively-charged impurities are included, n i = n e , and this means that the electromagnetic fields have proportionally less response to the bulk ions. Since the bulk thermal ions are responsible for the instability of the ITG mode, this can lead to a reduction in the vigor of the turbulence. This effect is known as dilution. For the same reason, dilution of the ions implies an enhancement of electron-driven microinstabilities, such as the trapped-electron mode [50] and the electrontemperature-gradient mode (ETG) [46]. In contrast to thermal high-Z impurities, which are known to stabilize ETG [51], the adiabatic response of singly-charged fast ions actually reduces Z Taking the radial derivative of equation (2) gives: In a local simulation, this provides a second independent constraint that must be satisfied. This represents dilution on neighboring flux surfaces, which influences the relationship between the electron and thermal ion density gradients. This also has a nontrivial effect on ITG turbulence and is included when we speak of 'dilution'. When results are thereby labelled, it means that fast ion fluctuations are not included in the gyrokinetic simulations, but their equilibrium effect on the density of the bulk species and on the magnetic drift (through their contribution to ∂β/∂r) are included. In other words, we keep the effects of F 0f and ∇F 0f , but not δf f . Note that, in the high-energy limit (explored in section 3), this is equivalent to the adiabatic approximation, In constructing artificial cases lacking fast ions (in order to isolate and study their effect on the ITG mode and resulting turbulence), there is a choice to be made regarding whether the electron or thermal ion density (and density gradients) are changed to maintain quasineutrality. Fusion products and energetic tails heated from the bulk ion population deplete the thermal ion population. On the other hand, the physical origin of injected and minority-heated fast ions are such that they are accompanied by excess electrons. The difference between these choices for the k y ρ i = 0.57 ITG mode is shown in figure 3. As can be seen, this choice happens to have a relatively small effect on our results, although this may not be true in general. In this section, we choose the convention that it is the electron properties that are changed.

Differentiating classes of fast ions
The different types of fast ions are characterized by whether their strong radial variation is one of particle density ('NBIlike') or energy density ('ICRH-like'). It has previously been shown that these different classes of fast ions respond differently to turbulence when passive [27,52]. It is worth examining to what extent their effect on turbulence differs compared to their respective dilution effects. This contrast is shown in figures 4(a) and (b). The violet spectra each include one of the respective fast ion species, and the dashed black spectra are the cases without the fast ions, for which the electron density and density gradients are changed to maintain quasineutrality. Then, cases identical to the violet ('with fast ions') were run, except that fast ions were not included as a kinetic species; only their effect on the equilibrium n e and L ne . These 'dilution' cases are shown in green. We see that, relative to dilution, NBI-type fast ions are actually destabilizing, whereas ICRH-type fast ions are more strongly stabilizing than dilution only. Therefore, the classification of fast ions, whether L Tf L nf (NBI) or L nf L Tf (ICRH), is critical to predicting and understanding their stabilizing effect. Note that alpha particles are considered NBI-like, but instead of being artificially injected, they are produced by the fusion reaction. Because the fusion source is a strong function of radius, the alpha particle density gradient is sharp. Now the fact that the two cases with fast ions differ from their respective 'dilution only' cases demonstrates that fast ions play a non-trivial role besides mere dilution. However, when both types are present, the (significant) kinetic effects of fast ions may cancel out and one can be left with the illusion that dilution is the dominant effect.
We saw from figure 4 that, relative to dilution, the NBIlike fast ions (strong density gradient) are destabilizing, while the ICRH-like fast ions (strong temperature gradient) are stabilizing. The parameter η f ≡ T f (r)/n f (r) = L nf /L Tf characterizes the relative strengths of the gradients. Since the stabilization is clearly a function of this parameter, its effect is examined in figure 5. Here, the temperature gradient of fast deuterium is modified from the baseline case and the growth rate is compared to the case of pure dilution (black lines). Note that the dilutive effect on the bulk plasma density is not affected when changing the temperature gradient. For both electrostatic and electromagnetic cases, the threshold in the value of η f for when the fast ions become more stabilizing than simple dilution is around 0.7-1.0. This threshold will be explained with an analytic model in section 3.4 after a model perturbed distribution function is obtained.
The effect of fast ions in local gyrokinetic turbulence simulations is often ascribed to electromagnetic effects. Here, some qualitative differences between different kinds of fast ions at different values of β e are catalogued. This is complicated by the fact that β e itself has its own electromagnetic stabilizing effect, and separating which changes to the growth rate are due to the electromagnetic stabilization and which are due to fast ions is not trivial. Therefore, scans in β e for the growth rate of the k y ρ i = 0.57 mode, using several different fast ion parameters, are shown in figure 6. In all these cases, increasing β e is stabilizing. However, increasing the   fast ion pressure gradient (either from the density gradient or the temperature gradient) does not universally decrease the ITG mode growth rate. In particular, note the case of NBIlike fast ions in figure 6(b) (a/L T,f = 0) at low β e : increasing a/L n,f is actually net-destabilizing electrostatically. Dilution cases are not shown in figure 6, but these still have a lower growth rate than their corresponding high-a/L nf cases, and is approximately equivalent in the a/L nf = 0, a/L Tf = 0 case, as expected. Therefore, the stabilization effect of fast ions and dilution itself is somewhat sensitive to β e [53]. The perturbed fast ion parallel current is small, but the direct effect of fast ions on the electrostatic potential could couple to A , causing the behavior shown in figure 6.
We have explored several potential explanations for the fast ion stabilization including: a modified linear growth rate, dilution of bulk ions, changes to magnetic geometry, and magnetic fluctuations. Although these explanations are physically motivated and in aggregate have a non-trivial effect, they are found lacking in describing the full order-of-magnitude stabilization observed in JET discharge 73224. A more fundamental treatment, valid in fully nonlinear turbulence, is thereby motivated. In the next section, a simplified model for the gyrokinetic equation, valid in the high-energy limit, will be derived and this will later be used to develop a reduced model to explain the qualitative effects presented in this section, in addition to the effect of fast ions in nonlinear simulations.

Reduced model for fast ion effect on microturbulence
In the previous section it was demonstrated that the effect of fast ions in stabilizing the ITG mode can go beyond dilution and depends on the details of the equilibrium fast ion phase space distribution. In this section, we examine the leadingorder behavior of the gyrokinetic equation in the high-energy limit. This subsidiary expansion yields an analytic solution that applies rigorously to fast ions with a strong radial dependence. Then, this model for h f is reduced further with additional assumptions based on the ITG mode structure. When this approximate distribution function is inserted into Maxwell's equations, this leads to a physically-transparent effective parameter model. We then discuss this model and benchmark it against linear and nonlinear gyrokinetic simulations.

Energetic expansion of the gyrokinetic equation
In order to obtain a nonlinear model for the effects of fast ions, we will directly expand the gyrokinetic equation, making use of the high-energy nature of the fast ions. Indeed, it is common to perform such expansions for electrons, taking advantage of their small mass. Commonly-used models include: driftkinetic [54][55][56], fluid [57,58], bounce-averaged [59], or adiabatic electron [60] models. The adiabatic electron model, for example, approximates their contribution to the field equations as proportional to the electrostatic potential. In this section, an analogous model, applicable to energetic ions, will be developed.
In the energetic limit of ≡ v ti /v tf 1, the gyrokinetic equation (1) In arriving at equation (4), ion-scale microturbulence is in mind so that the bulk ions are what set the temporal and spatial scales of the turbulent fluctuations so that ∂h f /∂t ∼ ωh f and |∇h f | ∼ h f /ρ i . The ∂h f /∂t term, the nonlinear term, and the first term on the right hand side of equation (1) are smaller than the magnetic drift term by factors of −2 , −3/2 , and −3 , respectively. The radial gradients of the fast ion equilibrium are ordered to be strong such that |∇F 0f | ∼ O −3/2 F 0f /a . To leading order, only the magnetic drift term survives and we recover the adiabatic approximation: h f ≈ 0, which is equivalent to 'dilution' in this limit. To find nontrivial effects, we therefore wish to find a solution for h f correct up to O ( ), which is why the other terms are retained in equation (4) even though they are formally smaller than the magnetic drift term by one power of ε.
As is customary, an eikonal representation is chosen for the fluctuations [61] so that χ =χ (θ) e iS , where b · ∇S = 0 and the ballooning angle θ ∈ (−∞, ∞) has been chosen as the coordinate along the magnetic field line. Hats will be dropped henceforth on h and χ where there is no ambiguity. The perpend icular spatial dependence of the fluctuations is embedded in S, which depends on the poloidal magnetic flux ψ (which labels the flux surface and is a useful 'radial' coordinate) and a field line label α such that B = ∇ψ × ∇α. We perform a Fourier transform in the plane spanned by ∇ψ and ∇α so that the gyroaverage is represented by the Bessel func- This can be solved directly for h f . For the integrating factor, it will be convenient to define: where ω D (θ) ≡ v D · ∇S. The primed variables denote which θ is the independent variable such that ω D ≡ ω D (θ ) (similarly for v , b · ∇θ, χ, J 0 , and z). In equation (6), the lower limit of integration, θ 0 , is defined to be a point where ω D vanishes such that θ 0 < θ if v > 0, and θ 0 > θ if v < 0. Equation (5) can then be rewritten as whose solution is: This solution is largely inspired by that of [34] and can alternatively be derived directly therefrom with the appropriate approximations. The equilibrium distribution can be factored out because F 0f = F 0f E, µ, σ is not a function of θ in these coordinates. Equation (8) is the full solution of the gyrokinetic equation for fast ions correct to O (v ti /v tf ). It is a complete linear model in the sense that it is the furthest one can take the (v ti /v tf ) expansion linearly; the formally next-largest term that would appear in equation (4) is the electromagnetic nonlinear Without further analysis, we can read off one important implication of this solution for h f : there is no zonal reponse of the fast ions. This follows from the fact that the zonal modes are defined as those without any α variation (in the context of gyrokinetic simulations, this is often written as k y = 0). Hence, for such modes ∂S /∂α = 0, and so h f (k y = 0) = 0. Zonal flows (which arise from zonal fluctuations of φ) are well known to be important in the nonlinear saturation of ITG turbulence. Our approximate solution shows that fast ions have no direct impact on the zonal fields. To demonstrate the robustness of this result, simulations have been performed in which h f (k y = 0) = 0 is artificially enforced. There, no discernible difference in the nonlinear fluxes was found; see figure 8.
This result is not without consequence. The saturated level of turbulence is determined by a balance between 'drive' and zonal flows, which interact nonlinearly. Fast ions only directly affect the former and not the latter. Decreasing the drive will allow stronger zonal flows [53] and an overall damping of the turbulence: moreso than if all modes (including the zonal modes) were directly damped. A similar mechanism applies to an adiabatic electron response, which also vanishes for the zonal mode [60], and this makes turbulence remarkably sensitive to the ion-electron temperature ratio [62]. This is a possible explanation for the nonlinear enhancement of the fast ion stabilization. Zonal flows play an important role in the story, but it is precisely because fast ions do not have a zonal response.

Simplifying the model
The rigorous solution, equation (8), still retains too much physics to be a useful reduced model. We can make further approximations which allow us to write h f ∝ χ. When moments of h f are thereby taken, the proportionality factors become response functions, which we will find very useful in interpreting the contribution of fast ions to the fluctuating electromagnetic fields.
Integrate equation (8) directly by parts to obtain the approximate solution: where the approximation is made because the second term is smaller by one power of ε. A difficulty with this approximation is that ω D vanishes at specific θ points. We thus need to assume that χ vanishes sufficiently rapidly away from θ = 0 that the contributions from such resonances are negligible due to the smallness of χ. This is referred to as the 'strongly ballooning' or 'outboard mid-plane localization' approximation. It is valid to the extent that the ITG mode structure peaks at the outboard midplane, which is typically the case. In practice, this is equivalent to approximating the integral θ −∞ ≈ θ θ0+∆θ , stopping just short of the nearest resonance point θ 0 . Alternatively, the solution (9) could also be found by ignoring the (formally small) parallel streaming term in equation (4), promoting ∇F 0f by an additional order in ε, and solving for h f algebraically.
The approximation in equation (9) is supported by the fact that it exhibits the correct energy dependence that we expect from the usual scalings of energetic particle transport in electrostatic and electromagnetic turbulence [30,[63][64][65]. The leading term in equation (9), which we will use for the fast ion response to fluctuating fields, does not contribute to the turbulent transport because it is in-phase with χ. The correction from the second term results in a nontrivial phase factor, which does contribute to turbulent transport with the correct energy dependence.

Model fast ion response function
In this section, the model distribution function of the previous section is reduced further to remove the θ dependence, taking θ = 0 as the only relevant location for ω D and J 0f . In this case, the v j moments of h f become related to numerical parameters R jf , which can be transparently interpreted in the gyrokinetic field equations.
First, consider the relevant Maxwell's equations in the gyrokinetic limit: s Z s δn s = 0 (Poisson's equation/quasineutrality) and ∇ 2 ⊥ A = − (4π/c) s δj s (the parallel component of Ampere's law with δj s ≡ Z s e δf s r v d 3 v being the contribution of each species to the perturbed parallel cur rent). Assume the plasma consists of a thermal ion species, electrons, and fast ions. Considering the low mass and high thermal speed of electrons, their contribution to the perturbed current dominates over the ions. At low β e , the electron contrib ution to the perturbed charge density is approximately their adiabatic response to the electrostatic potential: δn e = (n e e/T e ) φ − φ ψ , where ψ denotes the flux surface average. The field equations become: In equation (11), it was assumed that the thermal ion contribution to the parallel current is small compared to that of the electrons or the fast ions. Equation (11) is not inconsistent with employing the adiabatic electron approximation in equation (10). This is because the leading adiabatic behavior of electrons, even if dominant at low β e , does not contribute to the parallel current, while non-adiabatic corrections do contribute.
With singly-charged bulk thermal ions and a temperature ratio τ ≡ T i /T e , the non-zonal components of quasineutrality (we have already established that the zonal fast ion contribution is negligible) can be written: (12) To capture the gyrokinetic effect of fast ions in response to the fluctuating fields, let us define the response functions: where ω D0 ≡ ω D (θ = 0) (similarly, k ⊥0 and Ω f 0 take their values at θ = 0). Here, the mid-plane localization assumption of the previous section was taken further to assume that v D takes its values at the outboard midplane. This assumption is also applied to v , b · ∇θ, and Ω f , which is equivalent to taking the large aspect ratio approximation and ignoring trapped particle effects. Up until now, the magnetic geometry and equilibrium fast ion distribution are general. However, for the sake of straightforward parametrization, let us make the further assumptions of circular geometry and Maxwellian fast ions: This is the form of the response function that will be used henceforth. Note that R 1f vanishes by odd symmetry in v when F 0f is isotropic. This symmetry is important for the model that follows because the latter depends on the fast ions not coupling the electrostatic and electromagnetic field equations (this coupling is precisely R 1f ). Even if this symmetry is violated, one could still argue that the fast ion current is small (which removes their relevance in Ampere's law) and, when β e is small, v A /c φ, which means that the prefactor on R 1f is small in quasineutrality, equation (10).
The response functions R 0f and R 2f are shown in figure 7. Note that k ⊥ ρ f and η f consist the only nontrivial parameter dependency of the response functions; all other parameters appear as prefactors in equation (14). When |R 0f | 1, this means that dilution is the dominant electrostatic fast ion effect by definition. This allows us to make an important conclusion for alpha particles: unless accompanied by an unphysically strong radial gradient, the prefactor of n α /T α in equation (14) makes the electrostatic kinetic response of alpha particles very weak, in agreement with the results of [30]. We will find that the electromagnetic response R 2f has a relatively weak effect at thermal ion scales, even though it comes with an additional prefactor of T f /T i .

Effective parameter model
Having derived a simple fast ion response function in equation (14), we proceed to interpret it in the context of the gyrokinetic field equations. This is done by generalizing the model presented in [27] in a way that goes beyond dilution and electrostatic turbulence.
Consider that h i is proportional to n i , but otherwise does not depend on the equilibrium ion density, except through the calculation of φ in equation (12). After multiplying equation (12) by n e /n i , it is seen that the same φ is ensured to be calculated from h i for all cases where the quantity Equate this quantity to a case without fast ions where n i = n e , but with an artificially defined τ eff . This defines an effective temperature ratio which mimics effect of fast ions on microturbulence: This is a generalization of the dilution model presented in [27], which did not include the kinetic effect of fast ions  (14) have been normalized out. Note that R 2f is multiplied by an additional factor of T f /T i . approximated by R 0f . It is also analogous to the τ parametrization used in [51] for impurities in ETG. A similar calculation can be performed with the parallel component of Ampere's law, equation (11), which is rearranged thusly: By similar arguments, the same A will be calculated given an electron h e if and only if the bracketed factor is constant. This suggests an effective beta to mimic the electromagnetic response of fast ions: The dimensionless response functions R 0f and R 2f are both of order unity. However, when β e is small, fast ions have little electromagnetic effect except the part of the spectrum where k ⊥ ρ i 1. Furthermore, note that dilution could also have an electromagnetic effect in equation (17). This happens when the electron density, and thereby β e , changes to maintain quasineutrality in the presence of fast ions (as was done in section 2).
It is stressed that nowhere in this derivation was it assumed that the bulk plasma responds linearly to χ, only that the fast ions do. This is justified from the high-energy expansion of the gyrokinetic equation. Therefore, the τ eff model is expected to be valid in fully nonlinear turbulence, and is not strictly a linear model. In cases where the fast ions are not energetic enough to take such an expansion seriously, the only recourse is to perform the fully nonlinear multi-species electromagnetic gyrokinetic simulations. The 'brisk ions' must then be treated as a low-charge nearly-thermal impurity whose behavior is as difficult to predict as the bulk ions. But here we conclude that fast ions, to the extent that they can be classified as such, obey such an expansion and make the reduced model presented here a useful one for investigating their effect on microturbulence. We now compare this model against gyrokinetic simulations.

Benchmarking the reduced model
With the same basic plasma parameters as in section 2, gs2 simulations were performed to calculate the steady-state bulk ion heat flux. In this section, ions are used to balance quasineutrality. Therefore, in these cases, a/L ne = 0.422, and L ni is determined from equation (3). This is done to be consistent with the existing results in the literature and to avoid direct changes to β e . The various cases presented in this section are tabulated in appendix A. For nonlinear simulations simulations, the spectral range in the perpendicular direction goes up to k y ρ i = 2.1 with k y,min ρ i = 0.1, and similarly for k x . Along the field line, there are N θ = 30 grid points in each poloidal turn (of which there are 7 given the k x resolution) of the irrational flux surface. The velocity space resolution is N v × N λ = 18 × 32. The timestep is conservatively kept below the CFL condition [66], and time averages are consistently performed to be the last 60% of the simulation: typically averaging over a period of about 300-500 a/v ti units. For linear simulations, we examine the k y ρ i = 0.4, k x = 0 mode, and this is the mode for which we calculate τ eff corresponding to the nonlinear cases. Figure 8(a) shows the time-trace of the bulk ion heat flux q i for two different simulations: the baseline case that includes fast ions and has a heat flux that approximately matches the experimental power balance, and a case with the NBI and ICRH fast ions removed. The variation of the steady-state heat flux as the thermal ion temperature gradient changes is shown in figure 8(b). From there, one can see that, along with a change in the critical gradient, there is also a strong reduction in the slope. Also shown is another example with fast ions, but here magnetic fluctuations were removed from the simulation. This shows that electromagnetic fluctuations are clearly stabilizing in their own right. Also shown in figure 8(b) is the steady-state heat flux for a case with fast ions, but with their zonal component (h f (k y = 0)) artificially nullified. The fact that this case is nearly indistinguishable from the standard case implies that the fast ions have negligible direct effect on the zonal flows, as discussed in section 3.1.
We wish to use the τ eff model of equation (15) to estimate the expected strength of the fast ion stabilization. To this end, we consult an empirical scaling for the bulk ion heat diffusivity: χ i ∝ τ −3 [62]. Although it has not been derived from first principles, we can use this scaling as a useful indication for the approximate strength of fast ion stabilization based on τ eff . Furthermore, the presence of a carbon impurity in this case interferes with the direct calculation of τ eff (see section 3.6). Nevertheless, for the parameters of JET discharge 73224, one obtains τ eff = 2.9. If the empirical scaling is to be taken seriously, this implies an even stronger stabilization than observed in figure 8. This indicates that the strong stabilization observed in some cases can be at least qualitatively described by the simplified model presented here. In light of this analysis, the sensitivity of microturbulence to the presence of fast ions is no mystery. In fact, the model is oversensitive compared to simulation. The τ eff model is further benchmarked with a collection of additional linear and nonlinear simulations. These cases have the same baseline parameters as those before (tabulated in table A1), but with no carbon impurity and simplified fast ion parameters (table A3). The red circles in figure 9(a) represent the ITG mode growth rate at k y ρ i = 0.4 with kinetic electrons and an adjusted ion temperature, but no fast ions are present, even via dilution. The simulations with fast ions include either 'ICRH-like' (a/L nf = 0, a/L Tf = 5), or 'NBI-like' (a/L nf = 5, a/L Tf = 0), each with a nominal density of n f = 0.15n e and T f = 10T e . The cases NBI2 and ICRH2 each have n f = 0.2n e instead. Another case is 'alpha-like', which have a/L nf = 4.5, a/L Tf = 0.5, n f = 0.0075n e , Z f = 2, m f = 2m i , and T f = 200T e . The case 'BothFI' has two different fast ions species each at half density (n f = 0.075) with the respective gradient length scales and temperatures listed above. For these cases with fast ions, τ eff was calculated at the k y ρ i = 0.4 according to equation (15) and plotted accordingly on the horizontal axis. In most cases, the bulk plasma gradients and magnetic geometry were held fixed, despite changes in the fast ion parameters. The exception to this are the hollow boxes, in which the ion density gradient was changed consistently with the presence of NBI-like fast ions. Figure 9(b) shows the steady-state ion heat flux from the corresponding nonlinear simulations.
The fidelity of the τ eff model is indicated by the proximity of the blue squares to the red circles in figure 9. The growth rates show excellent agreement, while the nonlinear heat flux captures the general trend. For the nonlinear ICRH-like cases, the τ eff model significantly over-predicts the impact of fast ions. One explanation is that these cases are on the verge of being dominated by a fast ion-driven mode, possibly a fast ion-driven ITG mode. This is made clear when the fast ion temperature is increased to T f = 15T i and 20T i . For these modified cases, agreement with the τ eff model is much better, but they are not shown in figure 9 because it is not longer thermal-ITG and the agreement is likely coincidental. Another possible contribution to the τ eff model's under-prediction of the heat flux is that it makes use of the adiabatic electron model in the electrostatic field equation, while the simulations in figure 9 included kinetic electrons. Furthermore, τ eff depends on k ⊥ ρ i , so in nonlinear simulations, we must choose a single representative mode number if we wish to specify a single τ eff parameter.
A prediction of the model that matches particularly well with nonlinear electromagnetic gyrokinetic simulations is that alpha particle fluctuations play little role in stabilizing ITG turbulence. The Maxwellian alpha-particle-like species shown in figure 9 affect the bulk ion heat flux by only about 6%. This can be explained in the τ eff model by the large temperature and small density of these species. Note that high temperature, like those of alpha particles, is exactly where we expect the model to be most accurate. From equation (14), one can see that R 0f → 0 as T f /T i → ∞, which means that dilution is the dominant effect of such fast ions. Even the electromagnetic response is small for this case: β eff ≈ 1.014β e . As demonstrated for the electrostatic case in [30], even accounting for the non-Maxwellian nature of alpha particles has little effect on ITG-driven turbulence. We have neglected the changes on the magnetic geometry caused by the alpha particle pressure gradient. While this is known to have a significant impact on microturbulence, it is beyond the scope of the phenomenon we attempt to isolate in this work.
Surprisingly, in none of the cases shown in figure 9 does β eff depart by more than about 15% from β e . This is evident from equation (17) since R 2f if of order unity, but k ⊥ ρ i ∼ 1 and β e 1. Since fast ions can play a key role in driving electromagnetic modes unstable, one would expect their influence on the ITG mode to be associated with the β e stabilization. Indeed, low k ⊥ ρ i is relevant for Alfvén eigenmodes and energetic particle modes, and is where one might see a significant electromagnetic effect from fast ions. Therefore, according to the first-principles linear model, at ITG-relevant mode numbers, the effect of fast ions is mostly through their contribution to the electrostatic field equation. Their contribution to the electromagnetic field equation is precisely the perturbed fast ion current, which is small compared to that of electrons for the cases studied here. However, this does not rule out, for example, an indirect effect of τ or τ eff on the magnetic fluctuations. Furthermore, if the electron density were increased to maintain quasineutrality (as opposed to decreasing the thermal ion density as was the convention followed in this section), this has a direct change on β e , which can have its own significant effect on the turbulence.
The threshold for stabilization relative to dilution is determined by the sign of R 0f . One can obtain an analytic estimate by taking the k ⊥ ρ f → 0 limit in approximating the Bessel function in equation (14). In this case, one obtains a threshold η f = 1. As k ⊥ ρ f becomes finite and large (as is appropriate for fast ions in thermal ion-scale turbulence), this threshold can be estimated from numerical calculation of the integral in equation (14) (note that all other parameters are multiplicative factors and the threshold only depends on k ⊥ ρ f ). See figure 10 for these calculations. At high k ⊥ ρ f , the threshold approaches η f ≈ 0.70. This is in excellent agreement with figure 5, in which the observed threshold was close to but less than η f = 1.

Limitations of the reduced model
Many assumptions were made on the way to writing equations (15) and (17) that it is worth considering what has been lost.
One can see from figure 9 that the model does poorly in predicting the results of simulations with a relatively large concentration of NBI-like fast ions (the hollow black boxes). The reason for this is because they have a strong density gradient (such that ∇n f ∼ ∇n i ) and thereby have an effect on the density gradients of the bulk plasma via dilution. Note that the effect of local dilution is included in equation (15), but that on neighboring flux surfaces is not included even in principle. To show that this is responsible for the disagreement, other simulations were run where the ion density gradient was not changed, and these cases are also presented in figure 9 (labelled as 'NBI' and 'More NBI'), and there the model performs much better. Unfortunately for the model, it appears that the fast ion-induced change in L ni is the dominant stabilization effect of NBI-like fast ions. Although progress is being made [67], the general theory of how turbulence scales with the density gradient, especially when that of the ions and electrons differ, is not generally known. Note that this does not jeopardize the applicability of the model to alpha particles. Even though the alpha particle gradient scale length is short, their densities are so small that |∇n α | |∇n i | and the overall effect of dilution is weak.
In order to derive equation (15), we had to ignore other impurities, such as thermal carbon or tungsten. These are known to have an impact on turbulence [51,[68][69][70], but their response to ion-scale turbulent fields is not even approximately linear. Therefore, their contribution to equation (12) is difficult to predict a priori. The simulations shown in figure 9 did not have such an impurity, although the simulations in figure 8 did, consistent with the experiment.
While it is undisputed that electromagnetic fluctuations alone can have a strong effect on ITG turbulence, the reduced model presented here indicates that the fast ions play little role in this phenomenon (due to the small contrib ution of the fast ion current relative to that of electrons). Nevertheless, fast ions are known to destabilize Alfvén eigenmodes [32] and geodesic acoustic modes [71], and it is conceivable that this could play a role in interacting with microturbulence [36]. Garcia et al [24] studied such a case where beta-induced Alfvén eigenmodes co-existed with, and had a nontrivial effect on, ITG turbulence. The τ eff and β eff models are not expected to capture the effects of modes that are driven unstable by fast ions, only their approximate effect on thermal ion-driven modes.

Conclusion
This work presented a reduced model, derived analytically from first principles, for the effect of fast ions on ITG and Figure 10. Threshold η f for which fast ions are stabilizing/ destabilizing relative to dilution, as estimated from numerical calculation of equation (14).
other forms of ion-scale turbulence. This model is based on the high-energy limit of the gyrokinetic equation and on the ballooning structure of the ITG mode, and provides physical insight into the sensitivity of microturbulence to the presence of fast ions. Several important linear and nonlinear results were presented, which highlighted the differences between different classes of fast ions as they affect ITG turbulence.
It was found that fast ions with strong density gradients are less stabilizing than those with relatively large temperature gradients, all else being equal. The chief stabilizing effect of the former is dilution of the thermal ions, an effect which is tempered by the destabilizing kinetic response of these 'NBIlike' ions. Fusion-produced alpha particles in burning plasmas fall into this category, and their effect was found to be small, owing to their low density and high energy. The strong stabilization observed with 'ICRH-like' fast ions can be explained in light of the sensitivity of microturbulence to the thermal ion-electron temperature ratio, which acts as a proxy for the contribution of fast ions to Maxwell's equations even in cases when this contribution is not trivial.
To reliably predict the effect of fast ions when the separation of energy is not as extreme as it is for alpha particles, multi-species nonlinear gyrokinetic simulations are required. Nevertheless, the theory presented here provides a useful estimate for the baseline effect, to which more sophisticated physics can later be added. Ideas for expanding this model include predicting the impact of: thermal impurities, changes to the equilibrium thermal plasma density gradients, and nontrivial interaction with energetic particle-driven modes.
The authors would like to thank T. Fülöp, J. Citrin, M.J. Pueschel, and R. Bravenec for helpful discussions, data, and feedback. Simulations were run on the CINECA Marconi cluster. GW was supported by the Vetenskapsrådet  Table A1 lists the local geometrical parameters used for all the simulations in this work. These consist of the local geometric parameters, the thermal temperature gradients, and the Carbon properties (when present). These values do not change among all the different iterations in this work, except when explicitly scanned upon (such as, for example, the scan in a/L Ti in figure 3). Table A2 shows the fast ion parameters, along with the effects they may or may not have on the bulk plasma (depending on the specific case being studied-some are intentionally left nonquasineutral for demonstration purposes). These are the cases used for the linear simulations in section 2, while the nonlinear results for cases B0 and BXi are shown in figure 8. Case B0 is considered the 'baseline' case of this work and, along with table A1 is based on the parameters reported in [49].

Appendix A. Simulation paramters
Finally, in table A3, we tabulate the parameters used in the simulations shown in figure 9. This case is a simplified version of case B0, with simplified fast ions gradients and no carbon impurity. For most cases, the bulk plasma gradients were held fixed, except for the 'Alphas' case, and cases QN1 and QN2, where a/L ni was adjusted for quasineutrality.  Table A3. The plasma species parameters used for the nonlinear cases in section 3. There are no thermal impurities in these simulations. Electron properties were held fixed for these cases: T i /T e = 1 and a/L ne = 0.422. All fast ions are singly-charged deuterium except for the 'Alphas' case. Labels correspond to those in figure 9.  Table A2. The plasma species parameters used for the various cases in section 2. Carbon has the same temperature and temperature gradient as the bulk ions, n c = 0.039n i , a/L nC = 0.422, and T C = T C = T e . In general, ion density and density gradients were held fixed, except for case BXi.