Multi-scale analysis of global electromagnetic instabilities in ITER Pre-Fusion-Power Operation plasmas

Global electromagnetic gyrokinetic simulations are performed with the Particle-in-Cell code ORB5 for an ITER Pre-Fusion-Power-Operation (PFPO) plasma scenario, with half-field (2.65 T) and half-current (7.5 MA). We report on a 'multi-scale' analysis of the discharge, considering eigenmodes and instabilies across three scale-lengths. Although the scenario will nominally have neutral beam heating with particles injected with 1 MeV, Alfv\'en eigenmodes are investigated in the absence of such source, and Reversed Shear (RSAE), Toroidal (TAE) and Elliptical (EAE) Alfv\'en eigenmodes are found with weak damping for moderately low toroidal mode numbers ($10 \leq n \leq 35$). At higher toroidal mode numbers ($40 \leq n \leq 70$), unstable Alfv\'enic modes have been observed close to rational surfaces and are labelled as Beta-induced Alfv\'en eigenmodes (BAE)/Alfv\'enic Ion Temperature Gradient (AITG) modes, since their frequency is associated with the BAE gap and they are driven by the bulk plasma on the Alfv\'enic continuum. These modes are unstable in the absence of energetic particles, and adding a species of energetic particles (with an isotropic 1 MeV slowing down distribution) has negligible impact on their growth rate. At higher toroidal mode numbers ($150 \lessapprox n<220$), low frequency microscale instabilities are observed.


Introduction
Energetic particles (EPs), such as Neutral Beam Injection (NBI) or alpha particles generated by fusion reactions, can have a significant impact on the stability of electromagnetic instabilities in tokamak plasmas. In turn, these instabilities can lead to energetic particle transport, or losses, affecting the efficiency of the plasma heating, or potentially increasing wall heat loads. In addition, the bulk plasma can also interact with these electromagnetic instabilities, not only contributing to the damping, but in some cases also driving these instabilities.
In previous work, we studied the effect of alpha particle drive on Alfvén eigenmodes in a Q = 10 scenario for ITER. In this paper, we consider a broad range of electromagnetic instabilities potentially present in the pre-fusion-power-operation (PFPO) phase of ITER. Specifically, we consider a scenario with a half-field, half-current, hydrogen plasma, and 1 MeV neutral beam injection (NBI). We consider the Alfvén eigenmodes which may potentially be destabilized by the NBI ions, and also smaller scale instabilities, namely meso-scale Alfvénic modes such as Beta-induced Alfvén eigenmodes/Alfvénic Ion Temperature Gradient modes (BAE/AITG), as well as microscale instabilies such as those associated with (electromagnetic) turbulence.
In this work, we use the global electromagnetic gyrokinetic particle-in-cell code ORB5 [1,2] for the simulations, which are also compared to Alfvén continuum calculations using the linear global gyrokinetic eigenvalue solver LIGKA [3].
This paper is organized as follows: §2 describes the theoretical model used for the calculations, and the numerical implementation of the numerical tool ORB5; §3 outlines the scenario being studied in this work; §4 presents the results, split into three sections: linear Alfvén eigenmodes which might be driven unstable by EPs, unstable AEs driven by the bulk plasma, and instabilities associated with microturbulence. Finally, §5 presents a summary and outlook of the paper.

Physical model
In this work, the numerical results presented are obtained using the ORB5 code [1]. This is a global electromagnetic gyrokinetic particle-in-cell (PIC) code which uses markers to sample the 5D phase space (R, v , µ), with the equations independent of the gyrophase, and the magnetic moment (µ) of a marker constant in the absence of collisions.
The perturbed part of the distribution function is then evolved according to the gyrokinetic Vlasov equation (where subscript s denotes a plasma species s) according to the equations for the particle equations of motioṅ and with the gyroaveraged potential φ = φ(R + ρ)dα/(2π) with ρ the gyroradius of the particle and α the gyro-phase. These equations are coupled to the field equations, the linearized gyrokinetic quasineutrality equation For the parallel Ampère's law, the perturbed magnetic potential A is split in to the symplectic part A (s) , which is found from and the Hamiltonian part A (h) , solved using the mixedvariable parallel Ampère's law, with n 1s = d 6 Zδf s δ(R + ρ − x) the perturbed gyrocenter density, ρ s = √ m s T s /(q s B) the thermal gyroradius, q s the particle charge, d 6 Z = B * dRdv dµdα the phase space volume, j 1s = d 6 Zv δf s δ(R + ρ − x) the perturbed parallel gyrocenter current. Note that due to the scale separation between ion and electron Larmor radii, electrons are treated differently from ion species, and terms with ρ e are neglected.
The particular method of using the mixedvariables formulation, and then applying a pullback procedure to the particles is described in reference [2], where its advantages for solving electromagnetic equations are discussed.
A set of straight field line coordinates are used, with radial coordinate s = ψ/ψ edge (where ψ is the poloidal flux), toroidal angle ϕ, and poloidal angle (where θ is the geometric poloidal angle, and q(s) the safety factor) ‡. The field equations are solved on a basis of cubic finite elements and using Fourier methods in the two angular directions. This allows the use of Fourier filtering, which helps to reduce noise by excluding non-physical modes which are far from being field aligned [4]. The local poloidal Fourier filter will retain modes with poloidal mode numbers m n (s) = nq(s) ± ∆m at each radial point for each toroidal mode n, where a typical values of ∆m is 5.
The distribution function f (R, v , µ, t) is separated into a time-independent background part (F 0 (Ψ = s 2 , E = 1 2 v 2 + µB, v )), and a time dependent part δf (R, v , µ, t). The time dependent part of the distribution function is represented by markers, using the so-called particle-in-cell (PIC) Lagrangian method, with each "marker" representing a small 5D volume of phase-space, and carrying a particle weight, w(t). These markers are evolved using the equations of motion for the gyrocenters, and with a weight equation to evolve the weights. The fields are solved from the charge and current densities, calculated by depositing the charge/current from the markers using the gyrocenters for electron markers, and gyrorings (using in this work a 4-point average) for the ion markers.
The system of particles and fields is evolved using a 4th-order Runge-Kutta timestepping scheme. ‡ θ * in ORB5 is often labelled elsewhere as χ.
Linear simulations can be performed by linearizing equations for the particle trajectories and weights, equations 2, 3, and 4.

Scenario description
We consider a half-field, half-current ITER PFPO [5,6] scenario, modelled using ASTRA [7] and made available at ITER through the Integrated Modelling & Analysis Suite (IMAS) [8]. The major radius, R = 6.2 m; the minor radius, a = 2.0 m; the magnetic field on axis B 0 = 2.65 T; the bulk ions are H, with some additional Beryllium and Neon impurities (< 1%). The electron temperature profile is peaked, with the electron temperature on axis T e,s=0 ≈ 8 keV, and the electron temperature at the top of the pedestal (T e,ped. top ≈ 3 keV). The main ion temperature profile has very similar core and pedestal top values, but the peaking profile is slightly different, with the steeper gradient region closer to the core than in the case of the electron temperature. Impurity ion species are taken to have the same temperature profiles as the main ion species. The density profiles of the ions and electrons (linked via quasineutrality and the impurity density profiles) are almost flat, but marginally peaked. The on-axis electron density n e,s=0 ≈ 4 × 10 19 m −3 . The temperature and density profiles for the background species are shown in figure 2(a).
The energetic particles in the scenario are from a negative-ion neutral beam injection (NBI, or N-NBI) source with a single birth energy of 1 MeV. In this work, when considered, the energetic particles, are modelled using an isotropic slowing down distribution [9] with birth energy 1 MeV. If we calculate the equivalent temperature by integrating moments of this distribution function, we find that this corresponds to 242 keV or T EP /T e = 30.9 on axis, with the ratio decreasing to about 24 at the top of the pedestal. The profile of the equivalent temperature for the NBI is shown in figure 3. While this neglects the anisotropy present in a realistic distribution function (also considered in the ASTRA model), we leave an anisotropic analytical distribution function [10], or a numerical treatment of the NBI from ASTRA's Fokker-Plank solver to future work. In fact, we leave most of the modelling of drive from NBI ions to future work.
The profiles and equilibrium considered were taken from IMAS at ITER, the scenario in question being shot 101006, run 50. The equilibrium used was generated by running the CHEASE MHD equilibrium code [11] reading from IMAS, and outputting an file using the CHEASE-ORB5 interface. § Due to the steep gradients at the edge of the plasma, we consider at most the region of the plasma 0 ≤ s ≤ 0.9, which corresponds to the outermost interior ring plotted in figure 1(a).
Finally, while all simulations presented here are electromagnetic and performed with kinetic electrons, we increase the electron mass in order to reduce the numerical cost. In previous work [12], it was found that the mode damping was only weakly affected by increasing the electron mass, but that the numerical savings allowed more confidence in the numerical timestep convergence. Therefore, unless otherwise but converges to its own solution. As a result, the q-profile obtained in the output of CHEASE is not identical to the qprofile from the IDS, as written by ASTRA.  specified, the ratio of m H /m e shall be taken as 400.

Results
By changing the toroidal mode number(s) in the Fourier filter, it is possible to perform simulations of reduced systems, and thereby separate the physics scales of interest: of particular interest for linear studies. Unless stated otherwise, the simulations presented retain only a single toroidal mode number in the filter, and are initialized with a density perturbation (a perturbation on the particle weights depending only on the spatial coordinates) corresponding to the same toroidal mode number. The radial and poloidal structure of the perturbation may vary, but is typically either radially broad and approximately field aligned (as is used in this work), or radially narrower, with e.g. a pair of poloidal harmonics designed to initialize the TAE, EAE, etc, in which case it will be radially located at the corresponding rational/half-rational surface of the TAE/EAE.

Alfvén eigenmodes
We start by considering the range of Alfvén eigenmodes which might reasonably be expected to be driven unstable by the energetic particles, focussing on 10 < n < 35.
First, we show n = 12, for which we perform a simulation with reduced mass ratio m H /m e = 400.
In figure 4, we show the poloidal cross section of the electrostatic potential at t =40 000 ω −1 ci . For the frequency evolution, we show in figure 6 how we fit the frequency, taking the example of s = 0.4685 (the location of the peak in figure 5).
In the upper panel, we plot the evolution of the envelope of the toroidal mode |Φ n,s (t)|. Fitting this from t >20 × 10 3 ω −1 ci , we obtain a growth rate of γ =−1.95 × 10 −3 ω A . We can now normalize the signal by multiplying (Φ n,s (t)) by e −γt . We then perform a Discrete Fourier transform of this normalized signal (again, for t >20 × 10 3 ω −1 ci ), the results of which are in the lower panel. In this case, we find (ω, γ) = (0.222, −1.95 × 10 −3 ) ω A , or f = 51.4 kHz (γ/ω = −0.88 %). Such weak damping is sufficiently small that we might possibly expect these modes to be driven unstable when energetic particles are added [12].
If we perform the procedure described above on all radial points (independently), then we can stack the lower panels of figure 6 over radius, giving us a radial spectrogram, shown in figure 7. In this figure, we overlay the Alfvén continuum from LIGKA ¶ [13], both with and without the pressure upshift (BAE gap). This technique is designed to best identify frequencies, and we note that, since every radial position is normalized in a different way (to account for the growth/decay Taking the real part here is optional when dealing with a complex signal such as the toroidal/poloidal Fourier coefficients of the field. In such a case, we can refrain from taking the real part, and then the resulting analysis will allow us to differentiate between clockwise and anti-clockwise rotating modes. In cases where we analyse a real signal, however, the positive and negative frequency components in the output will necessarily be complex conjugates of each other, and the sign therefore has no meaning. ¶ See discussion in Appendix A regarding differences in the safety factor profile used for ORB5 and LIGKA. An alternative method to what we have described above, which can also isolate growth rates and frequencies from multiple components, is the DMusic algorithm [14], in which multiple complex frequencies are fit to the signal. However, where we have a single dominant component at a given radius, we appreciate the fast numerical speed and high intuitivity of the method that we have described. Considering the analysis described here, we can say that the least damped modes for the n = 12 case are a BAE/RSAE at s ≈ 0.45; and a TAE at s ≈ 0.65. In addition, we also observe a weak EAE in the core, and another high frequency mode (high frequency RSAE or NAE (Triangularity-induced Alfvén eigenmode)), located around s = 0.15 and s = 0.24 respectively.
Next, we consider a simulation with n = 16. Similar to the analysis for n = 12, we show in figure 8 the radial profiles of the poloidal harmonics of the electrostatic potential Φ. Here, we also mark with vertical lines the rational surfaces (solid) and halfrational surfaces (dashed). In figure 9, we again show the radial spectrogram of the analysis. Here we see some TAEs (s = (0.37, 0.52)), with an RSAE at the minimum of q around s = 0.46. We also observe a very weak high frequency mode in the EAE gap at s ≈ 0.25. Again, the observed damping is weak, with the measured fit damping rates −0.005 < γ/ω A ≤ 0 across a broad range of the plasma, in particular the region 0.4 ≤ s ≤ 0.6.
Next we consider n = 20, for which we increase  8  9  10  11  12  13  14  15  16  17  18  19  20   21  22  23  24  25  26  27  28  29  30  Log. fit scaling Lorentzian fit (N s , N θ * , N ϕ ) = (1024, 512, 128). The radial profile of the electrostatic potential is show in figure 10, showing a peak of two poloidal harmonics (m = 20 and 21), peaked at approximately s = 0.45. We note that, even though this modes structure strongly reminds us of a TAE, we do not have a half-rational surface at s = 0.45, located between two rational surfaces (s = 0.350, s = 0.537), but not on a half-rational surface. We show in figure 12 that the minimum of the safety factor in this region is very close to the half-rational value of a TAE. From the spectrogram in figure 11, we note that the frequency of this mode is close to the bottom of the TAE gap, similar to behaviour also observed in the ramp-up phase of tokamaks [15,16] We therefore label this mode an TAE-like RSAE. We observe the

Unstable mesoscale modes
Further increasing the toroidal mode number, we observe a transition to a very localized instability close to the q max at s = 0.27. We therefore present analysis for simulations performed on a reduced radial domain, 0.1 ≤ s ≤ 0.4, which retain the physical effects, but allows us to reduce the numerical cost for such simulations (both in terms of reducing the required radial resolution as well as the q max , and therefore the poloidal resolution required to resolve m ≈ nq max . We therefore run a simulation for n = 50 with, (N s , N θ * , N ϕ ) = (384, 480, 240). This simulation has been performed with the number of markers, N p = (64, 64, 16, 16) × 10 6 , for H, e − , Be, Ne.
In figure 17, we can see that the dominant mode is single poloidal harmonic (m = 56) peaked at s = 0.26.  We also see a secondary peak of m = 55 peaked at the rational surface at s = 0.187. In figure 18, we see that the mode is growing with a clearly identified growth rate and frequency, (ω, γ) = (0.1658, 6.7 × 10 −3 ) ω A , or f =38.4 kHz (γ/ω = 4.0 %).
If we extend this case to include other toroidal mode numbers, we find that after a fixed amount of time, the electrostatic potential on outboard midplane can be visualized in figure 19. Here we see a pattern with a near periodicity of ∆n of approximately 8. We can understand this pattern by noting that the local q max = 1.1191, and therefore the proximity of the rational surfaces to the q max is periodic with 1/(1 − q max ) = 8.4. This also explains why a peak in intensity is observed for n = 59, since the rational fraction 66/59 is the closest approximation of q max in  this range of n in the denominator.
In figure 20, we show the frequency and growth rate for these cases, both located at the q max (a), and at the location of the mode peak (b). Here we can further see the almost periodic behaviour in the growth and frequency.
Regarding the mode drive, we show in figure 21 the spectrogram, but we draw particular attention to the kinetic spectrum overplotted from LIGKA. This is obtained from the local kinetic model of LIGKA, and we have plotted in yellow the points where the imaginary part of the continuum frequency is positive. Since the mode that we observe corresponds with tip of this unstable continuum, and since the mode is close to the position where the q max is close to the value of a rational fraction, we label this mode   For the explanation of the overplotted Alfvén continuum, see figure 7. The Fourier analysis has been performed on the complex spatial Fourier coefficients, and therefore the sign of the frequency can be determined. In this case, the LIGKA data has been mirrored in the (ω) = 0 plane.
a BAE [17], or Alfvénic Ion Temperature Gradient (AITG) mode [18,19], although we note that the mode also has similarities with the RSAE. We find (ω, γ) = (−0.152, 6.15 × 10 −3 ) ω A . Wanting to be sure that the artificial electron mass is not responsible for the mode growth in the simulations, we perform also a simulation with . We therefore have confidence that this mode is driven unstable by the bulk plasma in the model considered, and not that this is a numerical instability. For a detailed discussion of the effect of the electron mass and the time step on the numerical convergence, we refer to Reference [12].

Effect of energetic particles on BAE/AITG
modes An interesting question is whether energetic particles might have any effect on such a mode such as this bulk-plasma-driven mode.
Still with the m H /m e and ∆t =0.5 ω −1 ci , we therefore add a simplified population of NBI particles, using an isotropic slowing down distribution with an injection energy of 1 MeV. The effective temperature of this distribution depends weakly on radius and is shown in figure 3. We consider a simplified density profile for the NBI, shown in figure 22, constructed with a tanh function.
We therefore state that energetic particles, at least those considered here, have no effect on the linear behaviour of these BAE/AITG modes.

Microscale instabilities
Considering now the microscales, we perform a study of large toroidal mode numbers, looking for ion gyroradius-scale micro-instabilities (such as those associated with turbulence).
We consider the radial domain 0.4 ≤ n ≤ 0.7, the region with the steepest electron temperature gradient and consider toroidal mode numbers in the range n > 120, increasing the toroidal resolution to keep N ϕ > 3n, ideally ≥ 4n, but also choosing convenient values for parallelization.
We begin below by showing the case of n = 180, with (N s , N θ * , N ϕ ) = (384, 1152, 576). In the late  linear phase, we find a peak at s = 0.467, from where we measure low or zero frequency, and a growth rate of γ = 4.15 × 10 5 s −1 .
In the case of a larger mode number, n = 216, we find that the mode is linearly stable. For a smaller mode number, n = 160, we find a much weaker growth rate, which allows us to measure more accurately the linear (non-zero) frequency, finding (5000 ≤ tω −1 ci ≤ 54300) ω = −4.7 kHz, γ = 5.9432 × 10 4 s −1 .
We have therefore identified this range of toroidal mode numbers, n < 216, with a peak of γ for n ≈ 180, which driven unstable by bulk-plasma gradients.

Summary and outlook
In conclusion, we have shown a multi-scale nature of instabilities in this ITER PFPO scenario. By means of global electromagnetic gyrokinetic simulations with ORB5, we have identified three classes of eigenmode/instability, namely the weakly damped 'macroscale' Alfvén eigenmodes in the range n = [10 → 35], a range of 'mesoscale' unstable bulk-plasma driven Alfvénic instabilities in the range n = [45 → 60], and a range of unstable 'microscale' instabilities in the range n = [150 → 220].
We have shown that, for an isotropic slowing down distribution function of NBI energetic particles with E birth = 1 MeV, the growth rates of these mesoscale modes appear insensitive to the addition of energetic particles, although a more realistic treatment of the EP distribution function is left for future work.
We have identified a number of effects (that we refer to as 'multi-scale') which would be of strong interest to consider in a unified setting. Since there are partially different radial locations key to the different single-scale phenomena, the nonlinear interplay of these different scales would interesting. Considering the numerical demands of dealing with large mode numbers, and time scale separation and radial domains, this would be very demanding, however we think that we have laid out the scales and phenomena which such an effort should retain. Finally, we note that the label 'macroscale' is often used to refer to MHD instabilities, for example fishbone and sawtooth physics, which are not included in this study, but are also relevant for this scenario. We note that the safety factor profile, q(s), is of major importance for the Alfvénic modes, and any sawtooth activity will lead to a change in this profile. It is therefore of particular interest to consider time-dependent scenarios, which should especially be considered with reduced models [20]. The study of these MHD phenomena (in isolation) is of interest from (gyro-)kinetic codes or (kinetic-) MHD models, especially since these might lead to a redistibution of the energetic particle profile. profile As noted previously, when we overlay the Alfvén continuum from LIGKA, we must note that the safety factor profile is not consistent between LIGKA and ORB5, since these two codes interface to different  Figure A1: Comparison of the safety factor profiles (ORB5, from CHEASE; LIGKA, from HELENA) in the core of the plasma. The "ORB5" line corresponds to the inset plot in figure 1(b). equilibrium solvers, HELENA [21] and CHEASE respectively. Since these equilibrium solvers have different boundary conditions at the magnetic axis, there is not an exact match of the q-profile. We have taken care to match the q-profile as closely as possible in the core region, especially to match the region relevant for §4.2. We therefore plot a comparison of the safety factor profiles in figure A1.