Global linear gyrokinetic simulation of energetic particle-driven instabilities in the LHD stellarator

Energetic particles are inherent to toroidal fusion systems and can drive instabilities in the Alfvén frequency range, leading to decreased heating efficiency, high heat fluxes on plasma-facing components, and decreased ignition margin. The applicability of global gyrokinetic simulation methods to macroscopic instabilities has now been demonstrated and it is natural to extend these methods to 3D configurations such as stellarators, tokamaks with 3D coils and reversed field pinch helical states. This has been achieved by coupling the GTC global gyrokinetic PIC model to the VMEC equilibrium model, including 3D effects in the field solvers and particle push. This paper demonstrates the application of this new capability to the linearized analysis of Alfvénic instabilities in the LHD stellarator. For normal shear iota profiles, toroidal Alfvén instabilities in the n  =  1 and 2 toroidal mode families are unstable with frequencies in the 75 to 110 kHz range. Also, an LHD case with non-monotonic shear is considered, indicating reductions in growth rate for the same energetic particle drive. Since 3D magnetic fields will be present to some extent in all fusion devices, the extension of gyrokinetic models to 3D configurations is an important step for the simulation of future fusion systems.


Introduction
Instabilities driven by energetic particle (EP) components are of interest for magnetic fusion concepts since they can lead to decreased heating efficiency, high heat fluxes on plasmafacing components, and decreased ignition margins for reactor systems. Since 3D magnetic field perturbations will be present to some extent in all toroidal configurations, the analysis of EP instabilities in 3D systems is an important goal for fusion simulations. To address this, the GTC global gyrokinetic particle-in-cell (PIC) model [1] has been adapted to the VMEC 3D equilibrium model [2], and 3D effects included in the field solvers and particle push. Initial applications of this model have been made [3] to EP instabilities in several stellarators (LHD, W7-X) and pedestal turbulence in tokamaks with 3D fields [4]. Other gyrokinetic models that have been developed for both stellarators and tokamaks include the EUTERPE [5] and GENE [6] models. EUTERPE is a particle-based approach, while GENE is a continuum model that solves the 5D kinetic equations of all species. Additionally, the MEGA model [7] is a hybrid MHD-kinetic approach that couples a particle description for the fast ions with a full MHD model for the thermal plasma component. MEGA is applicable to EP instabilities in both tokamaks and stellarators. Another hybrid model is FAR3D [8], which couples reduced MHD equations for the thermal plasma with a Landau closure model for the fast ions and is designed for 3D systems. These models are all based to varying degrees on the gyrokinetic approach [9][10][11][12], which incorporates both the guiding center dynamics of particle trajectories and the effects arising from the finite helical Larmor orbits that center upon the guiding center trajectory. GTC solves the gyrokinetic equation using particle-based methods; the feasibility of the PIC method for gyrokinetics was initially demonstrated by W. E. Lee [13]. The specific gyrokinetic approach with adiabatic electrons, as described below, was formulated [14] to avoid high frequency modes and the time step limitation related to the electron Courant condition. The gyrokinetic orderings ) are applicable to most plasma components and regimes of interest for magnetic confinement systems. Gyrokinetics constitutes the most advanced first principles model that is also feasible to apply to global energetic particle instabilities in magnetically confined plasmas. The gyrokinetic PIC method used by GTC couples particle stepping in fluctuating fields with self-consistent electromagnetic field solutions based on Poisson's and Ampere's laws (based on retaining the potential, φ, and parallel vector potential, A ∥ ). For particle based gyrokinetic models, the small electron mass presents a numerical difficulty for simultaneously treating the dynamics of ions and electrons in simulations. A fluid-kinetic electron model currently implemented in GTC overcomes this difficulty by expanding the electron drift kinetic equation using the electron-ion mass ratio as a small parameter. The model accurately recovers low frequency plasma dielectric responses and faithfully preserves linear and nonlinear wave-particle resonances. Maximum numerical efficiency is achieved by overcoming the electron Courant condition and suppressing tearing modes and high frequency modes, thus effectively suppressing electron noise. In GTC the parallel vector potential is separated into adiabatic and nonadiabatic components, similar to the mixed variables (symplectic/Hamiltonian) pullback transformation [15] used in EUTERPE to avoid the so-called 'cancellation' problem. GTC can address kinetic issues specific to 3D configurations, such as multiple trapping regions, particles that transition back and forth between trapped and passing, and orbit trajectories that are more non-local than in similar axisymmetric tokamak systems.
In this paper the application of GTC for the linear analysis of energetic particle instabilities that have been observed in the LHD stellarator is demonstrated. A parameter survey indicates that Alfvén modes similar to those observed in LHD resonate with injected beam ions and are predicted to be unstable. Predicting the onset and effects of these instabilities is of significant importance due to their impact on heating efficiency, plasma-facing component heat loads, and possible diagnostic use. The importance of non-axisymmetric effects in all toroidal devices motivates the development of comprehensive new gyrokinetic simulation methods that can apply across the full range of symmetry-breaking perturbation levels. The paper is organized as follows. In section 2, we describe the gyrokinetic model, the LHD profiles and equilibria that are used, and discuss the resonance conditions that allow the neutral beam ions to destabilize the shear Alfvén eigenmodes. Next in section 3, the gyrokinetic results are presented for both normal and non-monotonic shear discharges in LHD; these include variations of growth rates and real frequencies with beam parameters, the relation of the frequencies obtained with the shear Alfvén continuum gap structure, and the Alfvén mode structures. Finally, in section 4 the conclusions are presented.

Description of the GTC gyrokinetic model
GTC is a global gyrokinetic, full torus, electromagnetic particle-in-cell model [16] based on Boozer magnetic coordinates [17]. Computational efficiency is gained by modifying these coordinates to approximately follow field lines. The applicability of GTC to fast ion destabilized global Alfvénic instabilities in tokamaks has been extensively demonstrated [18][19][20]. The implementation of GTC primarily used in this paper can be classified as a gyrokinetic model with adiabatic electrons. As will be described below, the energetic particle and thermal ion species are treated using gyrokinetics, while an adiabatic fluid description is used for the electrons. This includes most of the physics expected to be of importance for Alfvén instabilities. Several tests have been made including the non-adiabatic gyrokinetic electron terms for cases given in this paper and this leads to about a 10% reduction of growth rate due to electron Landau damping effects, but no significant change in real frequency or mode structure; however, including such effects increases the computational requirements and will be left for future research. The gyrokinetic equation for the thermal and energetic ions (σ is the species index) is given below: This is supplemented by the equations of motion for the position X of the gyrocenters, and the parallel velocity v ∥ : where m σ , Z σ , Ω σ are the mass, charge number, and cyclotron frequency of species σ, μ is the magnetic moment, b 0 is the unit vector b 0 = B 0 /|B 0 |, and the E × B and magnetic drift velocities are given as: With the curvature and grad-B drifts given by In our model we exclude the compressional component of the magnetic field perturbation by assuming δB ∥ = 0 and representing the perturbed magnetic field as: In equations (1)-(7) the perturbed electrostatic φ potential, magnetic fields (δB), and vector potentials A ∥ are gyroaveraged quantities, evaluated at the gyrocenter's position.
Equation (1) is solved by evolving δ = w f f 0 / weights synchronized with the particle trajectories using the following weight evolution equation.
In our electromagnetic simulations we use the fluid-kinetic adiabatic electron model [16] based on the separation between adiabatic and non-adiabatic electron response. In this model the non-zonal component of parallel electric field in written as The effective potential δφ eff in the lowest adiabatic order Here we have used the Clebsch representation of the magnetic field δ ψ δψ α δα where ψ is the poloidal flux function, α = q(ψ)θ − ζ is a field line label, and θ, and ζ are the poloidal and toroidal angles of Boozer magnetic coordinates [17]. In equation (10) δn e is obtained by solving the electron continuity equation The electron current is calculated from the Amperé's law The perturbed electrostatic potential is calculated from the gyrokinetic Poisson's equation [21]  ( ) Where δn is the perturbed gyrocenter's density and δφ ∼ is the second-gyroaveraged potential defined as The perturbed magnetic field potential functions are obtained from Faraday's law Where δφ δφ δφ = − ind e ff . These equations (1)-(19) form a closed nonlinear system to lowest order in the electron-ion mass expansion. The electron non-adiabatic terms and kinetic equation have been presented in [16] and will not be given here since, as mentioned above, the calculations of this paper are based on the adiabatic fluid model for the electrons.
In order to test the use of this model for stellarators, parameters are chosen to approximately match an LHD regime where Alfvénic activity was observed [22], although some simplifications have been made in the profiles. Specifically, the thermal ion density, ion temperature, and electron temper ature profiles are taken as constant in order to null out the drives for other instabilities caused by core density and temper ature gradients. The electron density profile is determined from the quasi-neutrality condition n e = n ion + n fast-ion for the three species (electrons, ions, fast ions) included in the calculation. The fast ions and thermal ions are treated using gyrokinetics,  while the electrons are incorporated using an adiabatic fluid model [16]. The LHD major radius is 3.7 m; the magnetic field on axis is 0.62 T; ion and electron temperatures are 1 keV; the central electron density is 0.884 × 10 13 cm −3 ; the plasma and beam species were hydrogen. The fast ion component is model led as a Maxwellian distribution with a constant temper ature versus flux surface. GTC also includes options for slowing-down models of beam distributions; these will be considered in future research on EP instabilities in stellarators. The computational parameters for these calculations were: 60 radial grid points, 128 toroidal grid points, 200 poloidal grid points, 40 particles per cell for ions and fast ions, 20 particles per cell for electrons, uniform marker temperatures, and 4-point gyro-orbit stencils for ions and fast ions. The time step for LHD is limited to about 1/10 of that for a similar axisymmetric system. These resolution and time step limits for 3D systems lead to significant computational requirements and currently limit the extent of parameter/profile surveys. Another modelling issue is that stellarators generally tend to have higher fast ion losses through the last closed flux surface   than tokamaks; several methods have been tested in the simulations for taking these escaping fast ions into account. For the results given in this paper, as ions escape through the last surface, their δf weights are set to zero. For the LHD cases given in this paper, about 40% of the initial fast ion markers and 7% of the thermal ion markers are lost through either the outer or inner radial boundaries. These leave the simulation domain at early times and do not present an obvious limitation to the simulation time. They do, however, reduce the marker resolution near the boundary regions, and reduce numerical efficiency by the retention of markers that do not contribute. Resolution has been tested in a few cases by increasing the initial particles per cell up to 100 without significant changes in the results, indicating particle counts are adequate for the linear analysis presented in this paper. Techniques that reinsert escaping ions back at another location at the same magnetic field value and such that they drift back into the simulation domain are under development. The calculations reported here are based on version 0706 of GTC. The primary changes from versions used in the earlier applications of this model to stellarators are the use of a Gaussian drop-off in the fields for the edge and magnetic axis boundary conditions instead of a linear extrapolation, and the zero weight/no-reinsertion method described above for treating escaping fast ions.
In order to reduce noise levels and target specific instabilities, a Fourier mode filter is used. The filter takes effect between the field solve and particle steps and involves a fast Fourier transform of the field data, followed by a nulling out of components not included in the filter, and then an inverse fast Fourier transform of the fields before they are passed to the particle trajectory step. For simplicity, the calculations given in this paper are based on one toroidal mode with 8 poloidal modes for the filter. Specifically, for n = 1, m = 1-8 are used; for n = 2, m = 1-8 are used, etc. Previous calculations [3] have also included the toroidal field period coupled modes, e.g. n = 1, n = −9, n = 11, but have not indicated for LHD that significant changes in stability properties result from including the higher order modes, due to its relatively high aspect ratio ( R 0 / a ~ 6) and number of field periods (N fp = 10).

LHD test case
Two LHD cases are considered here. The first case is motivated by an LHD discharge [22] where Alfvén activity was observed with toroidal modes numbers n = 1 and 2. This case had β~3% , magnetic axis at R 0 = 3.7 m, and an iota profile with normal shear (increasing with radius). The shape of the LHD outer flux surface is shown in figure 1(a) with the colors indicating the magnetic field strength (red = higher, blue = lower). The second case has a non-monotonic (reversed) shear region in the iota profile near the center; such profiles have been produced in LHD [23] by using neutral beam current drive and appropriate plasma start-up programming. The iota profiles for these cases are given in figure 1(b), with the reversed shear region indicated for the second profile.
In figure 1(a) and in subsequent figures, the variable denoted as r a / is the flux surface label equal to (ψ/ψ edge ) 1/2 , where ψ is the toroidal magnetic flux. The variation in flux surface shape as the toroidal angle is incremented within a field period is shown in figure 2. Here the direction of increasing toroidal angle is counter-clockwise (VMEC convention).
As there are no direct measurements of the fast ion density profile, a set of model profiles (normalized to the central electron density) as given in figure 3 are used. For the normal shear case, the fast ion profile model consists of a centrally flattened region with an exponentially decaying region on the outside (solid lines). This profile shape has been chosen specifically to select out the toroidal Alfvén eigenmodes that reside in the outer gaps, which are expected to be the ones that were observed [22]. For the non-monotonic shear case, a centrally peaked sequence of profiles have been used (dotted lines) that match onto the profiles used in the normal shear case on the outside. These are used to provide instability drive in the central reversed shear region of the plasma.

Alfvén resonance conditions
The resonance conditions for wave-particle interactions in stellarators depend on the frequency, mode number and rotational transform profiles of the device through the relation [24]: ω µ ω where m, n are the mode numbers of the instability, µ ν , = equilibrium mode numbers, N fp = number of field periods (=10 for LHD), and j = 0, ±1 is a coupling parameter. ω θ and ω φ are the poloidal/ toroidal drift frequencies, which may be calculated by following orbit trajectories in the 3D equilibrium fields.
This has been evaluated for the normal shear LHD equilibrium and parameters of the observed [22] Alfvén instabilities, leading to the results shown in figure 4. The resonance lines cross over most of the regions of phase space encountered by passing particles, confirming that tangentially injected beams should excite such instabilities.

Normal shear discharge
The GTC gyrokinetic model starts with an initial field perturbation and integrates particle trajectories and electromagnetic fields (φ and A ) in time to follow the growth of EP driven Alfvén instabilities. The characteristic behaviour of an unstable Alfvén frequency mode is shown in figure 5, where perturbations oscillating in the Alfvén range frequency grow exponentially with all modes growing at close to the same rate. This example is for the normal shear case with n = 1, T fast = 120 keV, n fast (0)/n e (0) = 0.0185. The growth rate can be inferred from the slope of the curves in figure 5(b) and the frequency from a Fourier transform of the data in figure 5(a).
In some cases, especially for stellarators, several eigenmodes can be present, each with unique frequencies and growth rates. This leads initially to a modulational waveform; however, if followed long enough, one mode dominates. For simplicity, this paper will restrict its analysis to the time intervals where a single mode dominates.
Another useful diagnostic for displaying the instability growth and radial extent is shown in figure 6; here the RMS averaged evolution of the potential is shown as a function of radius and time for the normal shear case with n fast (0)/n e (0) = 0.0185. Figure 6(a) is for n = 1, T fast = 120 keV while figure 6(b) is for n = 2, T fast = 60 keV. As can be seen, the n = 2 has a narrower radial extent than the n = 1 case.
As the fast ion density is increased, the Alfvén instability drive increases, leading to larger growth rates. An example of the variation of the growth rate and frequency with fast ion density for an n = 1 TAE instability with T fast = 100 keV is   given in figure 7. Due to the increased simulation time needed to resolve smaller growth rates, it has not been feasible to determine the marginal stability threshold (growth rate = 0) with this model.
In figure 8 the effects of varying the fast ion temperature are given for n = 1 and n = 2 at n fast (0)/n e (0) = 0.022. In the LHD experiment, beams are injected at 180 keV. The energy moment of a slowing-down distribution with 180 keV birth energy is equal to that of the Maxwellian distribution used here at about 100 keV. While fast ion destabilized Alfvén modes involve wave-particle resonances, the gyrokinetic results show very broad peaks in the variation with beam energy. This is due to the fact that sideband couplings induce secondary resonances at other velocities, which can encompass a wider range of energies. 3D stellarator equilibria offer significantly more mode coupling combinations (i.e. sideband couplings) than tokamaks. Also, since the gyrokinetic model includes all of the various trapped particle populations that are present in 3D systems, there are many other resonant frequencies beyond the usual passing particle transit resonance that can be involved. The n = 2 results show more structure than the n = 1, with several maxima present, possibly due to interactions with different particle resonances. Calculations were also carried out for n = 3 and 4, but these mode families did   The frequency ranges shown in figure 8(b) can be related to shear Alfvén continua obtained from the STELLGAP code [25] with acoustic coupling effects [26] included. Continuum plots for n = 1 and n = 2 are shown in figure 12. Here the slow-sound approximation [27] has been used to simplify the plots. The dashed black lines indicate the frequency ranges of the data in figure 8(b) and indicate that the unstable modes reside in the upper part of the m = 1,2 gap for n = 1 and the upper part of the m = 2,3 gap for n = 2. The frequencies obtained from these gyrokinetic calculations (f ~ 75-82 kHz for n = 1) are somewhat higher than seen experimentally (f ~ 60-70 kHz for n = 1). This is likely due to the use of a flat ion density profile in the gyrokinetic model calculations.
Recently reported [7] reconstructions of the experimental profiles have indicated a hollow ion density profile with higher ion densities near the edge (leading to a lower Alfvén velocity) than assumed here.

Non-monotonic shear discharge
Non-monotonic (reversed) shear rotational transform profiles have been of significant interest in tokamaks. Such profiles lead to new branches of Alfvén modes (the RSAE or Reversed Shear Alfvén Eigenmode) that are typically dominated by a single poloidal mode number. The frequencies of the RSAE modes are more sensitive to the rotational transform profile and show more dynamic behaviour (upward/downward frequency sweeps) in experiments than the TAE modes [28]. RSAEs have also been associated with higher levels of fast ion transport. Non-monotonic shear profiles were formed in LHD [23] using strong neutral beam current drive at low plasma densities. In this case the non-monotonic region refers to a region with negative shear in rotational transform, since the typical stellarator transform profile increases toward the plasma edge (positive shear). The tokamak non-monotonic shear profile has the opposite direction of shear; i.e. a positive shear region superimposed on a dominantly negative shear iota profile. When LHD was operated in this mode, n = 1 and n = 0 GAM (geodesic acoustic mode) activity was observed. The n = 1 signal was characterized by frequency sweeping both upward and downward in frequency, followed by more steady-state frequency lines covering a range from 50 kHz up to 150 kHz. These modes have been analysed using the STELLGAP and AE3D models [29], resulting in stable eigenfrequencies in the ranges seen experimentally.  The GTC model has been applied to an LHD reversed shear profile case similar to those realized experimentally. In order to compare with the earlier normal shear results, the plasma profiles and parameters have been kept the same, except for the fast ion density profile. Both a peaked profile case (shown in figure 3) and an equivalent flat profile case have been used. The peaked profile was tested to determine if placing instability drive in the reversed shear region would produce an RSAE mode; the flat profile was used to allow direct comparison with the normal shear result. For the equivalent flattened profile, the n = 1 growth rate is reduced from that of the normal shear case by about 28% (from 24.9 to 17.9 × 10 3 s −1 ) and the frequency is increased from 79.7 kHz to 115 kHz. The difference in the mode amplitude growth between the normal and reversed shear cases is shown in figure 13(a). In the case of the peaked fast ion profile, the dominant mode remains radially located outside the reversed shear region, as shown in the rms amplitude growth versus time and radius in figure 13(b).
A clearly defined RSAE localized around the minimum in the iota profile did not emerge in the simulation. However, there is a secondary component present around r a / ~ 0.5, that can be seen in figures 13(b) and 14(a) and (b) and may be related to the reversed shear region. In the reversed shear case the primary mode is dominated by m = 2, 3, 4 components for both the peaked and flattened profile cases.
The Alfvén continuum gap for this case is displayed in figure 15 with the frequency and approximate radial extent of the mode indicated by the dashed black line. The mode is predominantly related to the m = 2, 3 gap near r a / = 0.7. There should also be reversed shear Alfvén modes present near r a / ~ 0.4 above the m = 3 and under the m = 4 continuum lines. Finding appropriate fast ion profiles and conditions to excite these modes will be the topic of future research.

Summary
The gyrokinetic GTC model has been adapted to general 3D configurations that can include stellarators, tokamaks with 3D effects, and reversed field pinch helical states. To demonstrate this capability, it is applied here to the LHD stellarator, looking specifically at a low-field, low density regime where Alfvénic activity was observed [22]. Unstable modes that reside near the upper frequencies of the Alfvén gaps are found for n = 1 and 2, but not for n = 3 or 4. Phase space resonance analysis also indicates that tangentially injected beam ions should readily couple to Alfvén modes in the observed frequency ranges. The n = 1 and 2 mode structures have a global characteristic and may be expected to impact fast ion confinement and heating efficiency. The evolution of these modes in some cases shows modulational effects related to multiple competing Alfvén instabilities at separate frequencies. A second LHD application described here is to regimes with non-monotonic shear rotational transform profiles. When compared for similar parameters and plasma profiles, the non-monotonic profile results in about a 28% reduction in the n = 1 growth rate. The dominant mode remains a TAE, but a subdominant coupling to the reversed shear region is apparent in the eigenfunctions. The GTC model is a comprehensive first principles electromagnetic gyrokinetic PIC method, and can include most of the relevant growth and damping effects expected to be of importance for these instabilities. This model can also provide a calibration reference for reduced models of these instabilities. The calculations presented here demonstrate its increasing usefulness for the analysis of Alfvénic instabilities in 3D systems. Future work with this model will add more realism and investigate the nonlinear consequences of these instabilities. For example, experimentally measured profiles for the thermal ion/electron temperature and density will be included; slowing-down beam ions distributions can be used instead of Maxwellian; fast ion density profiles derived from beam deposition models can be factored in; the next order (non-adiabatic) electron kinetic terms will be included; larger mode filter sets can be used (must be accompanied with a higher poloidal grid resolution); and particle reinsertion methods will be developed for 3D systems.