Impact of fast ions on density peaking in JET: fluid and gyrokinetic modeling

The effect of fast ions on turbulent particle transport, driven by ion temperature gradient ( ITG ) / trapped electron mode turbulence, is studied. Two neutral beam injection ( NBI ) heated JET discharges in different regimes are analyzed at the radial position ρ t  =  0.6, one of them an L-mode and the other one an H-mode discharge. Results obtained from the computationally ef ﬁ cient ﬂ uid model EDWM and the gyro-ﬂ uid model TGLF are compared to linear and nonlinear gyrokinetic GENE simulations as well as the experimentally obtained density peaking. In these models, the fast ions are treated as a dynamic species with a Maxwellian background distribution. The dependence of the zero particle ﬂ ux density gradient ( peaking factor ) on fast ion density, temperature and corresponding gradients, is investigated. The simulations show that the inclusion of a fast ion species has a stabilizing in ﬂ uence on the ITG mode and reduces the peaking of the main ion and electron density pro ﬁ les in the absence of sources. The models mostly reproduce the experimentally obtained density peaking for the L-mode discharge whereas the H-mode density peaking is signi ﬁ cantly underpredicted, indicating the importance of the NBI particle source for the H-mode density pro ﬁ le.


Introduction
The main ion density peaking is crucial for the performance of a fusion reactor.Obtaining high values of the central density is of particular importance since fusion power scales with the square of the density and the density at the plasma periphery needs to be below the Greenwald limit in order to avoid disruptions.The density profile is determined by a balance between outward diffusion and convection known as a 'pinch', which can be inward or outward, together with possible contributions from the particle sources.
A number of experimental studies have investigated the scaling of the observed density peaking with various key plasma parameters like collisionality and magnetic shear, see [1] and references therein.On the theory side, though much less effort has been devoted to turbulent particle transport compared to heat transport, it is known that the particle pinch driven by ion temperature gradient/trapped electron (ITG/ TE) mode turbulence can support a peaked density profile in the absence of particle sources [2,3].It is also known that increasing collisions reduces the ITG mode driven inward pinch, in qualitative agreement with the experimentally observed collisionality dependence of density peaking [1,3].
Recently, the contribution of the neutral beam injection (NBI) particle source on the density peaking was investigated [4] by performing several dimensionless collisionality scans in various plasma scenarios at JET.To determine the particle transport coefficients a gas puff modulation technique with high quality time-dependent density profile measurements was used.The experimental results indicate that for the JET Original content from this work may be used under the terms of the Creative Commons Attribution 3.0 licence.Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.
H-mode scans the NBI particle fueling contributes significantly to the peaking while for the L-mode scan the turbulent inward pinch is the main contributor.The results were supported by predictive transport simulations using reduced models as well as gyrokinetic simulations [4,5].
Fast ions, produced by the NBI can modify the turbulence and therefore potentially influence the turbulent transport.Previous modeling efforts show that the inclusion of a fast ion species has a stabilizing influence on ITG driven turbulence through a number of effects [6,7] including dilution of the main ions [8], Shafranov shift stabilization [9] and electromagnetic stabilization due to the suprathermal pressure gradients [10,11].How the fast ions affect the electron and main ion particle flux due to the turbulent diffusion and convection and the subsequent density peaking remains an open question.
In the present paper, we analyze two NBI heated JET discharges presented in [4] at the radial position ρ t =0.6 (square root of normalized toroidal flux), one of them an L-mode discharge and the other one an H-mode discharge.The density peaking in the absence of sources is determined by finding the density gradient corresponding to zero particle flux, known as the peaking factor.The analysis is performed using an extended version of the Weiland drift wave model (EDWM) [12], the gyro-fluid model TGLF [13] and linear and nonlinear gyrokinetic simulations using the GENE code [14].The focus is on interpretative comparison between transport models and the role of the NBI particle source by including a high energy (fast) ion species in the modeling.In particular, the dependence of the peaking factor on the fast ion density, temperature and their gradients is investigated and the predicted peaking is compared with experimental results.
The rest of the paper is organized as follows.In section 2 the models used to describe the ITG/TE mode turbulence are presented.Section 3 discusses the interpretative analysis of the discharges as well as comparison with experiments.Finally, conclusions are given in section 4.

Fluid model
The fluid model used to study particle transport is an extended version of the EDWM [12] and can treat an arbitrary number of species in a multi-fluid description.Each species is described by a set of fluid equations, consisting of a continuity, a parallel momentum and an energy equation, for the perturbations in density, parallel velocity and temperature, respectively.The model includes ion finite-Larmor-radius (FLR) effects and ion parallel motion together with electromagnetic and collisional effects.The free electrons are considered adiabatic and electron FLR effects are neglected due to the scaling of the Larmor radius with mass.The density perturbations for ions and trapped electrons are coupled through a quasi-neutrality condition and the resulting eigenvalue equation can be reduced to a set of coupled algebraic equations by using a strong ballooning approximation.Particle and heat fluxes are calculated assuming linear relations between field quantities together with a modified mixing length estimate for the turbulence saturation level [2].The fluxes are computed numerically by adding the separate transport contributions from ITG and TE modes at each wavenumber, assuming isotropic turbulence.The magnetic geometry in the model is circular, shifted ( a s ) including effects of elongation.
Previous versions of the model include FLR effects to the first order in the ion continuity and ion energy equations.When considering fast ions, with substantially larger Larmor radii, the higher order terms in the expansion are no longer small enough to be neglected.Therefore, following [15], we have expanded the model to contain full FLR effects in the continuity and energy equations for the main and fast ions.

Gyro-fluid model
The trapped-gyro-Landau-fluid model (TGLF) [13] is the latest generation gyro-Landau-fluid model.Like its predecessors, it includes kinetic effects such as gyro-averaging and Landau damping, though TGLF unifies trapped and passing particles to a single set of equations.Also, it treats dynamic electrons, main and fast ions and impurities as well as electron-ion collisions and electromagnetic effects, in a shaped Miller geometry [16].It is a quasi-linear model where a system of moments of the gyrokinetic equation are solved for the linear eigenmodes.The quasi-linear fluxes are calculated for the full spectrum of linear modes, 0.1<k θ ρ s <24, where k θ is the poloidal mode number, ρ s =c s /ω ci is the ion gyro radius, = c T m s e is the ion sound speed and ω ci is the ion cyclotron frequency.Here ρ s normalized to deuterium.The saturated potential is then modeled to fit a nonlinear gyrokinetic database.In this paper, we are using the spectral shift saturation model SAT1, presented in [17,18], which includes multiscale effects.We run the model standalone at a single radius with an input parameter file obtained from a JETTO [19] interpretative run.

Gyrokinetic model
The gyrokinetic simulations in this study were performed with GENE [14], which solves the nonlinear gyrokinetic Vlasov equations together with Maxwell's equations in order to find the distribution functions of the species, the electrostatic potential and the parallel components of the magnetic vector potential and magnetic field.It allows for an arbitrary number of kinetic species as well as collisions and electromagnetic effects.The collisions are modeled using a linearized Landau-Boltzmann collision operator.Magnetic fluctuations, A P , were included in all simulations while the low beta justifies to neglect the B P fluctuations (see table 1 for the specific beta values used in the simulation).For the magnetic background field either a numerical equilibrium from the EFIT code [20] or an analytic circular model are used.
Both linear and nonlinear simulations are performed in a flux-tube domain.The computational parameters used in the linear simulations have a resolution of 32×15 in the parallel and radial direction with 32 grid points in the parallel velocity direction and 32 magnetic moments.An initial value solver is typically used, in the cases where sub-dominant modes are presented an eigenvalue solver is used.For the nonlinear simulations a simulation domain in the perpendicular plane of [192,64].In the parallel direction 32 grid points are used and in the parallel velocity direction 32 grid points and 16 magnetic moments.The simulations are typically run up to a simulation time of t=300a/c s , where a is the minor radius of the last closed flux surface.The resolution and simulation domain are checked through convergence tests.

Simulation results
In this paper we analyze two NBI heated JET discharges, 79811 (L-mode with Carbon wall) and 87424 (H-mode with ILW) with experimental density and temperature profiles shown in figure 1.The profiles are time averaged over 1 second and include error bars evaluated at fixed diagnostic coordinates.The main plasma parameters for discharge 87424 are B 0 =3.4T (at magnetic axis), I p =2.5 MA, P NBI =22.4MW, q 95 =4.0,Z eff ≈1.2 and T e /T i ∼1 and for discharge 79811 B 0 =3.3T, I p =2.0 MA, P NBI =5.8MW, q 95 =5.7,Z eff ≈2.2 and T e /T i ∼1.The electron temperature and density were measured with Thomson scattering and the ion temperature and rotation for discharge 79811 were measured with charge exchange recombination spectroscopy (CXRS) using the Carbon 6+ line.The ion temperature is not available for the ILW discharge 87424 since the impurity level is to low for CXRS to measure.Based on similar JET discharges it is assumed that T i ∼T e .
The analysis is performed at the radial position r = 0.6 t , which is within the region ρ t =0.5 to ρ t =0.8 analyzed by the experimental technique in [4].First, in section 3.1 we consider kinetic deuterium ions and electrons with the discharge parameters given in table 1.We include collisions and electromagnetic effects in the simulations and, for simplicity, neglect the effect of impurities.In TGLF the effect of rotation is included using the spectral shift model while in the simulations with GENE and EDWM, rotation is not included.For these two discharges, the effect of rotation is small (see figure 3).Unless otherwise stated, a circular geometry is used in all models.Then, in section 3.2 we include fast deuterium ions in the simulations.The fast ions are included as a kinetic   s a q r q d d is the magnetic shear, q is the safety factor and B 0 is the magnetic field at the magnetic axis.79811 87424 species with a Maxwellian background distribution function.This is justified by recent GENE results [21] which showed a lack of sensitivity to the NBI fast ion distribution.

Density peaking without fast ions
In general, the particle flux of species s is given by is the normalized density gradient of species s and a is the minor radius.The first term inside the bracket is the diffusive contribution and the second is the convection.In steady state the density peaking in circular geometry can be expressed as where the first term is the peaking factor considered in this work and the second term is the contribution from the particle source.
Figure 2 shows the eigenvalue spectrum and particle flux spectrum for both discharges obtained with TGLF.The linear growth rate (γ), frequency (ω) and electron particle flux (Γ e ) are given as a function of wavenumber (k θ ρ s ) for the most dominant mode.Here, ρ s is normalized to hydrogen.The growth rate has a maximum around k θ ρ s =0.3 and the horizontal dashed lines, corresponding to the E × B shear rate, illustrate that the importance of rotation is small.The negative frequency indicates that the mode is drifting in the ion diamagnetic direction corresponding to an ITG mode.The magnitude of Γ e is normalized by the number of k θ ρ s -modes in the simulation, and is thus expressed in arbitrary units.We observe the prevalent outward flux from the smallest wavenumbers and an inward flux at higher wavenumbers (for 79811).For these discharges there is negligible contributions of Γ e at electron scales (ES) and it is therefore sufficient to only include ion scale (IS) turbulence where k θ ρ s <0.7.Note that the corresponding condition is k θ ρ s <1 should one normalize to deuterium.
To determine the electron peaking factor (PF e ) we perform a scan in the normalized electron density gradient (a L n e ), henceforth referred to only as the density gradient, to find the value corresponding to zero electron particle flux.In the scan the ion density gradient is adapted to preserve quasineutrality.The electron particle flux given in gyro-Bohm units is illustrated in figure 3(a) as a function of electron density gradient and calculated with TGLF using the parameters in table 1.Note that Γ e is a sum of the particle flux from multiple wavenumbers between 0 and 0.7, henceforth referred to as the ion scale spectrum.The resulting peaking factor, obtained in a pure ITG regime, is illustrated with vertical dashed lines in the figure.The particle flux spectrum, when the density gradient is equal to PF e , is depicted in figure 3(b).Here, one can see the balance of outward and inward transport from the different wavenumbers.
We start our model comparison by comparing the peaking factors from TGLF and linear GENE with nonlinear GENE, including effects of non-circular geometry, for both discharges (simulation setups described in sections 2.2 and 2.3).TGLF uses Miller geometry while in GENE the numerical equilibrium from EFIT is used.As shown in figure 4, the peaking factors from TGLF reproduce those of nonlinear GENE quite well for both the L-mode (79811) and the H-mode (87424).Also, when comparing linear GENE for a single wavenumber with nonlinear GENE simulations, the best match is for k θ ρ s =0.3 and 0.4 for 79811 and 87424, respectively.Therefore, for the remainder of this section we restrict ourselves to linear GENE and EDWM for single wavenumbers in addition to TGLF, which includes a wavenumber spectrum at ion scales.Furthermore, a circular cross section is used, since EDWM is more limited in terms of geometry.
In figure 5 the peaking factor from all three models, EDWM, TGLF and GENE, is compared with the experimentally obtained density peaking.All models are found to reproduce the density peaking in the L-mode (79811) discharge, which demonstrates that the main contribution to the peaking is from the turbulent inward pinch, in concurrence with the experimental results of [4].For the H-mode (87424) discharge, on the other hand, all of the models underpredict the density peaking, with PF e negative or close to zero, while the experimental peaking is even larger than the experimental peaking for the L-mode.This indicates that the NBI fueling is important for the steady state density profile in the H-mode discharge, in agreement with [4].
To rule out that this discrepancy is not merely due to uncertainties in the temperature gradient (a/L T ), we perform a sensitivity analysis of PF e with respect to a/L T .The result is presented in figure 6.For EDWM and GENE the result is shown for k θ ρ s =0.4,though the trend is the same at other wavenumbers.Here the electron and ion temperature gradients are set equal and varied around their experimental value.All three models predict an increase in PF e with a/L T .With EDWM and TGLF the increasing trend is weak while with GENE it is more pronounced.This increase is due to a change in electron temperature gradient, changing only the ion temperature gradient tends to reduce PF e .We note that even an increase in temperature gradient by 20% is in this case insufficient to reproduce the experimental density peaking of 87424.Therefore, any uncertainties in a/L T are unlikely to alter the conclusions in [4] regarding the importance of the particle source for the H-mode discharge.
In simulations, we observe (see for instance figure 5) that the L-mode (79811) is peaked while the H-mode (87424) has a rather flat density profile.The key local parameters that differ between the discharges are magnetic shear (ŝ), safety   factor (q), plasma beta (β) and collisionality (ν) (see table 1).To quantify the effect of each parameter we replace step by step the values of ŝ, q, β and ν in 87424 to those of 79811 using TGLF.This corresponds to an increase in magnetic shear and safety factor and a decrease in beta and collisionality.The result is displayed in figure 7 where there is a negligible change due to ŝ and q while with β and ν replaced the peaked profile of 79811 is essentially reproduced in the H-mode (87424).This is consistent with previously known results regarding the effect of these parameters on density peaking (see for instance [1,22] and references therein).

Fast ion effects
Next, fast deuterium ions are included in the simulations with the parameters in table 2. Since the fast ions are already present in the plasma and thus have affected the plasma parameters in table 1, we keep the electron density and its gradient fixed.In order to preserve quasi-neutrality we reduce the main ion density according to the fast ion concentration.We investigate the effect of the fast ions on the growth rate of the most dominant unstable mode and on the electron and main ion peaking factor.Fast particle driven modes were not present in these discharges and are not considered in this work.Also, the equilibrium plasma parameters such as the safety factor, magnetic shear and Shafranov shift as well as the magnetic geometry, will remain constant.
The effect on the mode growth rate by the inclusion of fast ions is shown in figure 8 for both discharges.These growth rate spectra are calculated with TGLF, though the trends are the same for the other models.The inclusion of fast ions reduces the mode growth rate slightly, which reduces the turbulence fluctuation level in concurrence with previous results [10].For these discharges, however, the reduction in growth rate is quite small.
Without fast ions the electron and ion peaking factors are equal since we do not include any other species in the simulations.With the inclusion of fast ions they differ and one can consider the effect on the ion and electron peaking factors separately.Figure 9 illustrates the effect on the electron particle flux, (a) and (b), and the main ion particle flux, (c) and (d), due to their inclusion.The peaking factors are obtained with TGLF without and including fast ions.In the scan in electron (a) or ion (c) density gradient, the fast ion density gradient is kept fixed.First, we note the similarities between the electrons and ions in terms of trends.Both display a reduction of the peaking factor by the inclusion of fast ions, as can be seen in figures 9(a) and (c).Also, there is a net positive ion and electron particle flux for the density gradient corresponding to the peaking factor without fast ions, as can be seen in figures 9(b) and (d).In order to reestablish the balance of convective and diffusive flux, the zero flux density gradient is reduced.This tendency is the same for both ions and electrons, although it is more pronounced for the ions since their density is reduced by the inclusion of fast ions and are therefore more affected by effects due to dilution.
This reduction of the peaking factor is supported by the ion particle flux spectra in figure 10 from nonlinear GENE simulations including effects of a non-circular geometry.These spectra are taken at = a L 0.014 without and including 10% fast ions, respectively, for 87424 and both correspond to approximately zero flux.In TGLF, the reduction is similar when a Miller geometry is used.PF i is reduced from 0.045 (in figure 4) to −0.13 with 10% fast ions.Now we investigate further which fast ion parameters most affect the density peaking in the three models.First, the scalings of PF i with fast ion density is presented in  Effect on PF e by replacing the values of ŝ, q, β and ν in discharge 87424 to those in discharge 79811 (see table 1), calculated with TGLF.wavenumbers, though.Next, the sensitivity of the peaking factor to the fast ion density and temperature gradients is investigated.For the NBI heated discharges studied here,  L L T n f f (see table 2).For ICRH heated discharges, however, the relation between the fast ion gradients tends to be the opposite,  L L n T f f .Figure 12 shows the dependence of PF i on the fast ion density gradient, obtained with TGLF, for different fast ion temperature gradients for discharge 87424 (H-mode) covering both the NBI-like and ICRH-like regions.In this analysis other effects from ICRH are neglected.As shown in the figure, the decrease in PF i is enhanced by larger fast ion density or temperature gradients.We also note the tendency of lower PF i when the fast ions are ICRH-like.A similar trend on the effect on the mode growth rate by the inclusion of NBI and ICRH heated fast ions was found in [7].
We end this section by highlighting the effect on the PF e when the fast ions are included with the experimental density   of 4.9% and 3.0% in discharge 79811 and 87424, respectively, in comparison to the experimental density peaking (see figure 13).We note the greater reduction in PF e for experimental parameters similar to those of discharge 79811 than those of 87424, where the difference is negligible.

Conclusions
In the present paper we studied the effect of fast ions on particle transport due to ITG/TE mode turbulence using the fluid model EDWM, the gyro-fluid model TGLF as well as linear and nonlinear simulations using the gyrokinetic code GENE.The analysis was performed with experimental data from two NBI heated JET discharges, one L-mode and one H-mode at the radial position ρ t =0.6.The zero flux density gradient (peaking factor) was determined by performing scans in the electron (or ion) density gradient.In the scan the ion (or electron) density gradient was adapted in order to preserve quasi-neutrality while keeping the fast ion density gradient fixed.In EDWM and linear GENE we considered single wavenumbers while in TGLF the ion scale wavenumber spectrum was used in the analysis.In general the results of the reduced models reproduced the gyrokinetic results qualitatively.All three models predicted a peaked density profile for the L-mode while the H-mode was predicted to have a rather flat density profile, in the absence of sources.This difference in zero flux peaking between the two discharges was found to be mainly due to lower values of beta and collisionality in the L-mode discharge.Nonlinear GENE and TGLF simulations predicted a peaking factor of 62% and 88% of the experimental density peaking for the L-mode and 4.3% and 4.9% for the H-mode when effects of a non-circular geometry were taken into account.
The experimental fast ion density at ρ t =0.6 is 3.0% and 4.9% in the L-mode and H-mode discharge, respectively.All three models predicted a slight reduction of the linear growth rate of the mode and a reduction of the peaking factor by the inclusion of a fast ion species.The reduction increased with increasing fast ion density, temperature or corresponding gradients.When the experimental fast ion parameters were used the reduction in peaking was noticiable in the L-mode and negligible in the H-mode due to the lower values of fast ion density, temperature and corresponding gradients in the H-mode discharge.
For the L-mode discharge, the simulations reproduced most of the experimental density peaking.Hence, the main contribution to the peaking comes from the turbulent inward pinch as opposed to the NBI particle source, in agreement with [4].For the H-mode discharge, the simulated peaking factor was much lower than the experimentally obtained density peaking, even more so by the inclusion of fast ions.This shows that the NBI particle fueling plays an important role for the steady state density profile in the H-mode discharge, also in agreement with [4].
Note that, even though inclusion of fast ions reduces the zero flux peaking factor, the density peaking is also determined by the source, see equation (2).In order to determine their effect on the experimental density peaking one also needs to evaluate their effect on the diffusion and convection coefficient separately and not only on their ratio (PF).A stabilization of the turbulence by the inclusion of fast ions can thus in principle act to increase the density peaking rather than decrese it.Future work will aim to investigate this as well as the influence of impurities and rotation on density peaking.Grant Agreement No 633053.The views and opinions expressed herein do not necessarily reflect those of the European Commission.
The simulations were performed on resources provided by the Swedish National Infrastructure for Computing (SNIC) at PDC-HPC and the MARCONI (CINECA, Bologna, Italy) supercomputer systems.2).

Figure 1 .
Figure 1.Density and temperature profiles of the H-mode discharge 87424 and the L-mode discharge 79811, time averaged for 1 s.

Figure 2 .
Figure 2. TGLF eigenvalue spectra in s /a and electron particle flux spectra in arbitrary units for the discharges 79811 (L-mode) and 87424 (H-mode).The horizontal dashed lines correspond to the E × B shear rate.

Figure 4 .
Figure 4. PF e from linear and nonlinear GENE and from TGLF(Miller).

Figure 5 .
Figure 5. PF e three different models (circular geometry) compared to the experimental density peaking.

6 .
Scaling of e with respect to temperature gradient,

Figure 7 .
Figure 7. Effect on PF e by replacing the values of ŝ, q, β and ν in discharge 87424 to those in discharge 79811 (see table1), calculated with TGLF.

Figure
Figure growth rate spectra for the most dominant mode without including fast ions.

Figure 9 .
Figure 9. Particle flux, calculated with TGLF, without and including fast ions.

Figure 10 .
Figure 10.Ion particle flux spectrum from nonlinear GENE (EFIT) without and including fast ions when = a L 0.014 ni

Figure 11 .
Figure 11.PF i for k θ ρ s =0.4 in EDWM and GENE and ion scale spectrum in TGLF w/o and including fast ions with varying density and temperature.Vertical dashed line corresponds to experimental fast ion density.

Figure 12 .
Figure 12.PF i dependence on fast particle density-and temperature gradient from TGLF.

Figure 13 .
Figure 13.PF e , from three different models compared to the experimental density peaking.The hatched bars are including fast ions with experimental parameters (see table2).

Table 1 .
Discharge parameters at ρ t =0.6.Here, a/L n =−a∇n/n and a/L T =−a∇T/T are the normalized density and temperature gradients computed with respect to ρ t , b

Table 2 .
[23]meters for the NBI generated fast deuterium ions, obtained from PENCIL[23], at ρ t =0.6.The beam power is P NBI =5.8MW and P NBI =22.4MWfor79811 and 87424, respectively.dischargeTf /T i n f /n e a L nf a L Tf The vertical dashed lines in the figure represent the experimental fast ion density.PF i is calculated with the experimental temperature, T f =20.1TiandT f =19.7Tifor79811 and 87424, respectively, and with the lower temperature of T f =10T i .All three models predict the same trends, which is a decrease of PF i with fast ion density and a more pronounced decrease for a higher fast ion temperature.However, EDWM is mainly sensitive to the fast ion density, while TGLF and GENE are sensitive to both.Also, this decrease in peaking is stronger in GENE than in the other two models.In figure11, the result is shown for k θ ρ s =0.4 for EDWM and GENE.The trends are the same at other