Gyrokinetic particle-in-cell simulations of electromagnetic turbulence in the presence of fast particles and global modes

Global simulations of electromagnetic turbulence, collisionless tearing modes, and Alfven Eigenmodes in the presence of fast particles are carried out using the gyrokinetic particle-in-cell codes ORB5 (E. Lanti et al, Comp. Phys. Comm., ${\bf 251}$, 107072 (2020)) and EUTERPE (V. Kornilov et al, Phys. Plasmas, ${\bf 11}$, 3196 (2004)) in tokamak and stellarator geometries. Computational feasibility of simulating such complex coupled systems is demonstrated.


Introduction
Creation and control of burning plasmas is an ultimate goal of the magnetic fusion worldwide effort. Such plasmas will become experimentally accessible in the foreseeable future when the ITER and SPARC facilities start their operation. One of the characteristic features of the burning plasmas is the intrinsic richness of their physics featuring complex couplings and interactions of the microscopic processes (turbulence) with the macroscopic ones (MHD and Alfvénic modes). These interactions can affect properties of the turbulence as well as the evolution of macroscopic instabilities. Such couplings may become especially strong in burning plasmas where fast particles are abundant since they can drive the macroscopic modes unstable. A global approach is needed to assess the physics combining macroscopic modes and turbulence.
Global simulations of the electromagnetic turbulence can be particularly difficult due to the variety of the physics involved, but also as a consequence of the numerical problems, such as the cancellation problem [1,2]. Involving all the relevant time scales ranging from the electron motion to slow MHD dynamics (e.g. growth and nonlinear evolution of the tearing instability) and electromagnetic turbulence adds to the computational difficulty of the global gyrokinetic simulations in the electromagnetic regime. Addressing this problem in its full complexity for the reactor-scale plasma will require the emerging exascale computing. However, first steps in this direction can already be made now using the existing high-performing computing systems and employing the available codes. This has been the main goal of the PRACE computing project EMGKPIC. In this paper we report on the basic approach and the main results of the project. Because of the limited space, the results will be described in a brief manner leaving details to follow-up publications.
In the project, the global gyrokinetic particle-in-cell codes ORB5 [3] and EUTERPE [4] have been used to simulate electromagnetic turbulence in realistic tokamak and stellarator geometries. Demonstrating feasibility of global turbulence simulations using the electromagnetic gyrokinetic PIC codes was the first goal of the project. The global setup permits combining electromagnetic turbulence with global modes, such as tearing instabilities and Alfvénic Eigenmodes in the presence of the energetic particles. Assessing such combinations and affordability of their simulations was another goal of the project. A particularly strong emphasis has been put on the simulations of electromagnetic turbulence which is known to be notoriously challenging in a global setup. Combining it with the macroscopic physics (Alfvénic and MHD) is less difficult if the turbulence problem is solved. The following "case studies" have been identified: (i) global electromagnetic turbulence in tokamak geometry (e. g. ASDEX-Upgrade) (ii) tearing instabilities and their combination with the electromagnetic turbulence (iii) nonlinear Alfvénic modes in the presence of fast particles and their combination with the electromagnetic turbulence (iv) global electromagnetic turbulence in stellarator geometry (Wendelstein 7-X) The ORB5 and EUTERPE codes share the equations solved, the basic discretization principles, and many aspects of technical implementation (see Ref. [5] for details). Despite these similarities, they remain separate projects with a different set of capabilities. Thus, it is only EUTERPE which can address stellarator geometries. On the other hand, tokamak simulations using ORB5 are normally more efficient due to the tokamak axisymmetry explicitly employed in the code. Therefore in the project, all tokamak simulations have been performed using ORB5 and all stellarator simulations using EUTERPE.
In this paper, we report global electromagnetic simulations solving gyrokinetic equations and employing the numerical schemes described in detail in Ref. [5]. According to the "case studies" selected and described above, the first basic component considered in these simulations is electromagnetic turbulence nonlinearly evolving and saturating through zonal flow excitation or relaxation of the plasma profiles. In the finite-beta regime, the Kinetic Ballooning Mode (KBM) turbulence is believed to play an important role [6,7,8,9,10]. We identify the KBM regime performing parameter scans with respect to the plasma beta. We compare the low-beta Ion-Temperature-Gradient driven (ITG) turbulence with the KBM turbulence at a higher beta in the same magnetic configuration and for the same plasma profile shape. The second basic component of the global physics considered here is an MHD perturbation, namely the collisionless tearing mode, coupled to electromagnetic turbulence. We provide an example of this instability evolving separately and in the presence of the background turbulence. A basic feasibility of such simulations is demonstrated. The third component of the global physics is associated with macroscopic nonlinear Alfvénic Eigenmodes destabilized by the fast particles. Frequency evolution of a Toroidal Alfvén Eigenmode in the presence of electromagnetic turbulence is considered. Finally, the electromagnetic turbulence in non-axisymmetric stellarator geometry is addressed and feasibility of stellarator simulations is demonstrated. These examples encompass the physics content computationally available for future in-depth studies using ORB5 and EUTERPE.
The structure of this paper is as follows. In Sec. 2, we consider electromagnetic turbulence in tokamak plasmas. In Sec. 3, the evolution of the tearing mode in the presence of electromagnetic turbulence is addressed. In Sec. 4, the coupled nonlinear system including Alfvén Eigenmodes, fast particles, and turbulence is considered. In Sec. 5, the stellarator simulations are described. Finally, conclusions are made in Sec. 6.

Electromagnetic turbulence in tokamak plasmas
First, we consider tokamak geometry with the aspect ratio A = 10, cirular cross-sections, the safety factors q(ρ) = 0.8 + 0.8ρ 2 and q(ρ) = 1.1 + 0.8ρ 2 (two cases compared), where ρ is the radius of the circular flux surface, and the temperature and density profiles: Here, s = ψ/ψ a , ψ is the poloidal magnetic flux, ψ a is the poloidal magnetic flux at the plasma edge, s 0 = 0.5, κ n = 0.3, and ∆ T = ∆ n = 0.208. Two different temperature gradients will be considered corresponding to κ T = 1.0 and κ T = 2.0. The ubiquitous nonlinear generation of the small scales in the phase space (filamentation) is controlled with the Krook operator, see Ref. [11] for further details. The machine size is determined by L x = 2r a /ρ s = 360 with r a the minor radius and ρ s the characteristic bulk-ion gyroradius, the ion-to-electron mass ratio is m i /m e = 200. The character of the turbulence is defined by the plasma β: the electromagnetic-ITG regime for small β changing into the KBM regime when β increases. This dependence is shown in Fig. 1 where the linear growth rate is plotted as a function of β. The growth rate is computed from the linear evolution of the perturbed heat flux. One sees a characteristic ITG-KBM transition for different plasma temperature and safety factor profiles. The nonlinear evolution of the radial heat flux is shown in Fig. 2 for the both regimes. In Fig. 2, one sees that the simulation enters the nonlinear phase after the linear growth in both (a) the ITG regime and (b) the KBM regime. The latter case shows an oscillatory relaxation dynamics resulting from the combination of the profile flattening by the turbulence and an energy source, provided here by the Krook operator. This profile relaxation process is also indicated by the evolution of the electrostatic potential shown in Fig. 4. One sees here how the finger-like structures develop and propagate outwards resulting in a flattening of the temperature profile shown in Fig. 5(b). At later times, the temperature gradient restores to some extent leading to the next relaxation cycle. In contrast, the ITG regime shows a much weaker relaxation of the temperature profile, see Fig. 5(a). The electrostatic potential, shown in Fig. 3, develops the characteristic pattern of the turbulent eddies decorrelated by the zonal flow. The evolution of the density profile is shown in Fig. 6.
As a next step, consider a realistic ASDEX-Upgrade equilibrium, similar to the "NLED-AUG" case [12], down-scaled to L x = 300, with the safety factor profile shown in Fig. 7(a); temperature and density profiles shown in Figs. 11 and 12. Similar to the previous case, the ITG-KBM transition can be seen in Fig. 7(b) where the growth rate is shown as a function of plasma β computed using the linear phase of the perturbed energy flux evolution. The perturbed energy flux is shown in Fig. 8(a) for the ITG regime and in Fig. 8(b) for the KBM regime. One can see that the energy flux has entered the nonlinear phase of its evolution. The values of the flux measured in the gyro-Bohm units are comparable in both regimes. The evolution of the electrostatic The profile relaxation is much weaker for the ITG turbulence than in the KBM case. potential is shown in Fig. 9 for the ITG regime, and in Fig. 10 for the KBM regime. The temperature and density nonlinear profile evolution is plotted in Figs. 11 and 12. One can see that the temperature profile relaxation caused by the KBM turbulence is not very strong for the shaped ASDEX-Upgrade geometry, in contrast to the circular crosssection tokamak case shown in Fig. 5(b). In future, the effect of the plasma shaping on the electromagnetic turbulence in tokamaks will be addressed in more detail. In Fig. 12(a), one can see a density peaking corresponding to a turbulent particle pinch in the ITG regime. In contrast for the KBM turbulence, the particle flux is outward, see Fig. 12(b). Note that there is a weak particle pinch also for the large-aspect ratio tokamak case in the ITG regime. It can be barely seen in Fig. 6(a). The particle flux reverses its sign in the KBM regime in the both tokamak geometries considered here. In future work, the question of the particle fluxes driven by the electromagnetic turbulence will be addressed in more detail.

Tearing instability and electromagnetic turbulence
We consider a tokamak with circular cross-sections, the aspect ratio A = 10, the machine size L x = 200, and the safety factor [13,14]: with ρ the radial coordinate, r a the minor radius, q a = 3.5, and ν = 1. This safety factor profile has a q = 2 resonance at s = 0.58 and is unstable with respect to the tearing mode. In gyrokinetic simulations, the tearing instability drive is included via the shifted Maxwellian distribution function for the electrons with the parallel-velocity shift determined by the ambient parallel current, similar to Ref. [13]. In Fig. 13(d), the evolution of the tearing instability is shown for the parallel magnetic potential at the s = 0.55 flux surface. One sees that the mode enters the nonlinear phase. The simulation is quite expensive computationally because of a relatively low growth rate of the collisionless tearing instability and the physical importance of the collisionless electron skin depth which is resolved in the simulations. The reduced mass ratio m i /m e = 200 and β = 0.4% are used. The mode considered in Fig. 13 is driven only by the profile of the ambient parallel current. The plasma Figure 14. Evolution of the electrostatic potential in tokamak plasmas during the collisionless tearing instability. Flat plasma temperature and density are considered. One sees how an island develops in time with the X-points clearly visible. Figure 15. Evolution of the electrostatic potential in tearing-unstable tokamak plasma with finite temperature and density gradients. One sees that the microturbulence is excited in addition to the tearing component. Also, a global mode, possibly of an Alfvénic nature, is excited at later times. temperature and density profiles are flat. The snapshots of the n = 1 harmonic of the perturbation are shown in Fig. 13 for the linear stage at ω ci t = 0.25 × 10 6 : (a) the parallel magnetic potential and (c) the electrostatic potential. One sees that a typical tearing mode structures has developed. In Fig. 13(b), the radial derivative of the parallel magnetic potential is plotted. This quantity is closely related to the tearing instability parameter ∆ . As expected, it has a sharp maximum at the resonant flux surface.
In Fig. 14, the evolution of the electrostatic potential is shown in the poloidal crosssection which includes all Fourier harmonics. One can see the tearing mode linearly excited and growing into the nonlinear island structure. The X-points are clearly visible in the nonlinear stage. Toroidal mode numbers −30 ≤ n ≤ 30 are kept in this simulation and the diagonal poloidal filter is used with the half-width ∆m = 5. For comparison, evolution of the electrostatic potential in the case of finite temperature and density gradients is shown in Fig. 15. Here, the plasma profiles are given by Eqs. (1) and (2) with s 0 = 0.5, κ n = 0.1, κ T = 0.5, and ∆ T = ∆ n = 0.2. One sees that the microturbulence develops in addition to the tearing mode excited in the early phase of the simulations via the choice of the initial weights of the markers. At the later stage, a global mode, possibly of an Alfvénic nature, appears. We will study this case in detail in our future work. Here, we just note that it could give an example of a self-consistent combination Figure 16. Evolution of the TAE n = 2 frequency for the case of flat bulk plasma temperature. Here, the windowed Fourier transform of the electrostatic potential is shown for the toroildal mode number n = 2 and the poloidal mode numbers −5 ≤ m ≤ −2. One sees that the frequency of the mode remains at the toroidicityinduced gap. The continuum is computed using the slow-sound approximation. Figure 17. Simulations using the same parameters as in Fig. 16 except the bulkplasma temperature which has now a weak gradient κ T = 0.4, see Eq. (4). One sees that the frequency evolves along the continuum branch. The same quantity is shown as in Fig. 16 (the Fourier transform of the electrostatic potential).
of the MHD tearing instability, microturbulence, and, eventually, an Alfvénic mode and a zonal flow. Such interlinked physical systems [15] can be addressed only within a global nonlinear framework and represent the main target of our work.

Alfvén Eigenmodes in the presence of fast ions and electromagnetic turbulence
In this Section, we address the issue of coupling between the turbulence and Alfvén Eigenmodes destabilized by the fast ions. We consider a circular cross-section tokamak with the aspect ratio A = 10 and the machine size L x = 350. The magnetic geometry is that of the ITPA benchmark [16,17]. The simulation includes the toroidal Fourier modes −40 ≤ n ≤ 40 and all the poloidal modes inside a diagonal filter with the halfwidth ∆m = 5. All particle species (bulk ions, electrons, and fast ions) are kinetic and nonlinear. The reduced mass ratio m i /m e = 200 and β = 1.72% are used. For the fast particles, Deuterium ions are selected with the temperature T f /T i = 40 and the density n f /n i = 0.005. Here, T i is the bulk-ion temperature and n i is the bulk-ion density. The Maxwell distribution function is chosen for the fast ions. The fast-ion temperature Figure 18. Electrostatic potential plotted as a function of the flux surface label s and the toroidal mode number n. Only the relevant part of the spectrum is shown. One sees that there is a mode activity in the case with the finite temperature gradient (b). In contrast, the case of the flat temperature gradient (a) does not have the mode activity. A consequence of the mode activity in (b) is a higher amplitude of the perturbation in the nonlinear phase which can be related to the characteristic frequency evolution observed in Fig. 17.
In Fig. 16, the frequency evolution is shown for the case of a flat bulk-ion profile for the toroidal Fourier harmonic n = 2. It corresponds to the dominant TAE instability for the parameters considered. There are other TAEs, for example at the toroidal mode number n = 6, which are however less unstable. One can see that the frequency remains at the toroidicity-induced gap for all times of the nonlinear evolution. The shear-Alfvén continuum is computed using the CONTI code [18,19] within the slowsound approximation.
In contrast, the frequency of the same TAE evolves along the continuum branch when the bulk-plasma species has a temperature gradient, see Fig. 17. Here, we apply a rather shallow bulk-plasma temperature profile given by the expression: with ρ the minor radius, κ T = 0.4, ∆ T = 0.04, ρ 0 = 0.5, and ∆ρ = 0.4. All other physical and computational parameters remain unchanged. A possible explanation for this frequency evolution is a higher amplitude of the Fourier harmonics corresponding to the drift instabilities. An indication for such instabilities excited by the temperature profile considered can be seen in Fig. 18. Here, the elecrostatic potential is plotted as a function of the flux surface label and the toroidal mode number. One can see a coherent instability developing at n ≈ 20 when the bulk-plasma temperature has a gradient. Clearly, more work needs to be done in future in order to fully understand the physical mechanism behind these observations.

Stellarator turbulence
In this Section, we apply the EUTERPE code to study the electromagnetic turbulence in a stellarator plasma. We consider a Wendelstein 7-X [20] configuration falling into the class of variants which exploit the rotational transform ι = 5/6 and the corresponding edge islands for divertor operation. Furthermore, the specific case chosen here is characterized by a very high magnetic mirror field, on the magnetic axis (B max − B min )/(2 < B >) ≈ 0.25, which has to be compared to the magnetic mirror of 0.1 in the Wendelstein 7-X standard case. On the one hand, this magnetic mirror leads to a small bootstrap current, on the other hand, it prevents the formation of a vacuum-field magnetic well. In fact, this configuration even possesses a vacuum-field magnetic hill, which means that, when going from the magnetic axis to the plasma boundary, the enclosed volume grows more strongly than the enclosed magnetic flux. Ideal MHD theory shows that such magnetic configurations are plagued by the lack of ideal MHD stability. In Ref. [21], a similar W7-X variant was found to be unstable against low-mode-number ideal MHD modes at even low plasma pressure. For the configuration studied here a free-boundary plasma equilibrium was computed using the ideal MHD equilibrium code VMEC [22]. For a plasma volume of V plasma = 25 m 3 , a uniform number density of n 0 = 1.4 × 10 19 m −3 and the temperature profile shown in Fig. 20(a) correspond to β = 1.52% at the center of the simulation volume, i. e. at the normalized toroidal flux s = 0.5. The characteristic crosssections illustrating the non-axisymmetric geometry of the stellarator and the rotational transform are shown in Fig. 19. The reduced mass ratio m i /m e = 200 is used. The resulting turbulent evolution of the heat flux is shown in Fig. 20(b). One can see that the simulation clearly enters the nonlinear phase.
The electrostatic potential at different times is shown in Fig. 21 as a a function of the flux surface label and the toroidal mode number. Here, one can see that a linear instability develops at early times. The low-mode-number part of the spectrum (which includes the zonal flow) is excited when the turbulence evolves. At later times, the temperature profile is modified via the turbulent relaxation. In Fig. 22, the electrostatic potential is shown in the poloidal cross-section corresponding to the toroidal angle ϕ = 0. One sees that the perturbation covers only a small fraction of the poloidal angles in the linear regime. It spreads over the whole poloidal domain in the nonlinear phase. In future, we will study global electromagnetic turbulence in stellarator plasma in detail.

Conclusions
In this paper, we have considered the electromagnetic turbulence in tokamak plasmas. The ITG-KBM transition has been identified and the relaxation of the profiles in the case of the KBM turbulence has been observed. For the large-aspect-ratio tokamak, the mode saturation in the KBM regime appears to be due to the flattening of the temperature profile caused by an outward transport of the turbulent finger-like structures. We have considered the electromagnetic turbulence in a down-scaled ASDEX-Upgrade plasma where a similar ITG-KBM transition has been observed as well. The profile relaxation is weaker in the ASDEX-Upgrade for the parameters considered. The multiscale physics has been addressed coupling the electromagnetic turbulence to the collisionless tearing instability and Alfvénic modes destabilized by the fast ions. Electromagnetic turbulence has also been addressed in the non-axisymmetric stellarator geometry.
A typical computational effort for the cases considered here requires large jobs on Marconi (128 nodes for many days of the simulation duration) or Marconi100 (64 nodes). In future, exascale computing systems will be needed for the reactor-size turbulence simulations and for massive parameter scans using the global codes. The spatial resolution can also become an issue when the collisionless reconnection is of importance requiring the collisioneless electron skin depth to be resolved. Extending the simulations to the realistic mass ratio does not normally pose a problem. Of course, such simulations are computationally more expensive because of a smaller time step needed. This is the main reason to use a reduced mass ratio throughout this paper. A number of simulations at the realistic mass ratio has also been performed for various scenarios and will be reported elsewhere. Since the collisionless electron skin depth is affected when the electron mass is decreased, the radial resolution requirements can become more challenging for the realistic mass ratio when the tearing mode is present.
The main goal of this paper has been to present a set of the characteristic physics applications which can be used as a starting point for future more detailed numerical experiments with ORB5 and EUTERPE. Here, we demonstrate that the simulations of such kind are possible using existing global gyrokinetic particle-in-cell codes on the available HPC systems.
awarding us access to Marconi100 at CINECA, Italy. We acknowledge PRACE for awarding us access to Joliot-Curie at GENCI@CEA, France.