The effects of anisotropic pressure on plasma displacement and its deviation from flux surfaces

A significant impact of pressure anisotropy on the plasma displacement associated with ideal Magnetohydrodynamic (MHD) stability is found. A heliotron plasma, such as a large helical device plasma, is analyzed. Simulations are performed using the equilibrium solver Anisotropic Neumann Inverse Moments Equilibrium Code and the ideal MHD stability code TERPSICHORE. Both codes provide a treatment of the pressure anisotropy by the bi-Maxwellian model. The ratio of hot particle pressure to total pressure / has been scanned over. Other simulation parameters have been chosen such that the simulations represent an experimentally relevant condition with an external, off-axis heating scheme. The radial location of the peak of the plasma displacement of the n = 1, m = 2 mode number has been compared to the radial location of the ι = 0.5 resonant surface. This comparison shows that this difference in location increases monotonically for increasing / . These results provide insight on the effect of external heating schemes on MHD stability.


Introduction
In both tokamak and stellarator type fusion devices, typically multiple heating schemes are employed to achieve control over the plasma temperature. In the Large Helical Device (LHD), these heating schemes include Neutral Beam Injection (NBI), Ion Cyclotron Resonance Heating (ICRH), and Electron Cyclotron Resonance Heating (ECRH) [1]. The deposited energy generally has a dominant component either * Author to whom any correspondence should be addressed.
Original Content from this work may be used under the terms of the Creative Commons Attribution 4.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. parallel or perpendicular to the magnetic field lines, as a consequence of the physical mechanisms underlying these heating schemes. In a low magnetic field configuration of the LHD, it has been shown that the ratio of parallel to perpendicular kinetic energy can achieve a value of 4 due to heating with NBI [2]. An anisotropy in plasma pressure parallel and perpendicular to the magnetic field lines has shown to have a pronounced effect on the Magnetohydrodynamic (MHD) stability inside the LHD [3][4][5]. Currently, a quantitative analysis of the plasma displacement for experimentally relevant mode numbers does not exist for LHD. The purpose of this work is to provide insight into the effects that pressure anisotropy has on MHD stability and in particular how this affects plasma displacement. Moreover, this work focuses on low-n modes, which comprise the most experimentally relevant mode numbers [1,6].

Magnetic equilibria and MHD stability
In MHD theory the velocity space distribution is often approximated by a Maxwellian distribution. This distribution, however, assumes an isotropic pressure, which is a sufficiently accurate model for the thermal particles in the plasma. The hot particles, however, can exhibit a strong pressure anisotropy as a consequence of external heating schemes. A modified version of the bi-Maxwellian distribution function [7] can be used to model this plasma anisotropy and serves as an extension to the Maxwellian distribution. The bi-Maxwellian distribution function which is used to model the hot particles, is given by, where N (s) is a density-like factor, m h is the mass of the hot ions, µ is the magnetic moment given by where v ⊥ is the velocity component perpendicular to the field lines, B is the magnetic field strength, B C is the critical magnetic field, E is the total particle energy and T ∥ and T ⊥ are the parallel and perpendicular components of the temperature, respectively. The parameter B C in this equation controls the location at which particles are deposited. It has been shown that B · ∇F h = 0, meaning that the bi-Maxwellian forms a zeroth order solution to the Fokker-Planck equation [8].
In order to model MHD stability in fusion devices, an equilibrium solver called Anisotropic Neumann Inverse Moments Equilibrium Code (ANIMEC) is employed [9]. This code forms an extension to the well-known Variational Moments Equilibrium Code [10] and models thermal and hot particle components separately, making use of equation (1) for the latter. This code calculates a magnetic equilibrium by minimizing the following energy where B is the local magnetic field strength, Γ is the adiabatic index and p ∥ (s, B) is the pressure parallel to the field lines. In order to simplify the mathematics pertaining to the equilibrium calculations, the adiabatic index has been set to zero. The implications of setting the adiabatic index to zero have been detailed in [4]. Note that s is the normalized flux surface label which relates to the normalized radial coordinate ρ as s = ρ 2 . In addition to ANIMEC, an ideal MHD stability code called TERPSICHORE [11], has been employed. This code has been extended to include energy calculations that take into account plasma anisotropy. The fully-interacting Kruskal-Oberman (KO) energy principle [12] is used in this research, which is a specific formulation of the MHD energy balance. The MHD energy balance evaluated in TERPSICHORE, can be summarized by where W P denotes the potential energy in the plasma, W V is the vacuum energy, ω is a complex frequency and W K is the plasma kinetic energy. The operator < · > indicates a volume average over the entire simulated plasma volume. The potential energy can be further divided into the following terms where W C 2 is a stabilizing term associated with field line bending and plasma compression, W BI describes pressure-driven instabilities and W J describes current-driven kink modes. These energy terms are described in more detail in [3]. The KO energy principle in combination with the modified bi-Maxwellian distribution function (1) forms an MHD stability model that considers the effects of hot particles. The ANIMEC and TERPSICHORE codes complement each other, allowing for MHD stability calculations, including mode structure and plasma displacement calculations, for different levels of pressure anisotropy. The total normalized pressure < β > inside a fusion device can be separated into a thermal and hot part, such that, < β > = < β th > + < β h >. The pressure anisotropy in a plasma can be characterized by separating < β > into a parallel and perpendicular component. The following definitions of the normalized pressure as well as its parallel and perpendicular components are employed In these equations p ∥ and p ⊥ are the parallel and perpendicular components of the total pressure p. In ANIMEC, the pressure anisotropy can be influenced by specifying the ratio of perpendicular to parallel hot particle temperature as a function of s, (T ⊥ /T ∥ )(s), as well as the critical magnetic field B C . The ANIMEC code has been used to calculate magnetic equilibria for a low magnetic field configuration of the LHD. Typical results of these calculations are shown in figure 1. This figure shows the parallel and perpendicular components of the pressure, respectively, corresponding to the hot particle component, in a poloidal cross section. These simulations have been performed at a normalized total pressure of < β > ≈ 3%. These figures clearly show that the parallel and perpendicular hot particle pressures are generally not flux surface quantities. This is in contrast to the thermal pressure which, following from isotropic MHD theory, is a flux surface quantity.

Simulation set-up
In order to realistically model the effect of heating schemes on MHD stability, the ANIMEC code takes as input thermal and hot pressure profiles. For this research, a thermal pressure profile of p th ∝ 1 − s and a hot pressure profile of p h ∝ s(1 − s) are used. This hot particle pressure profile is an off-axis profile which therefore serves as a model for off-axis heating schemes such as NBI, ICRH and ECRH.
A scan has been performed over the ratio < β h > / < β >. This ratio will be scanned over the following set of values, {0, 1/10, 1/3, 1/2, 2/3, 9/10}. For these simulations we fix < β > ≈ 3%. A critical magnetic field of B C = 2.3T has been used which corresponds to particle deposition at the low field side. In addition, we have set (T ⊥ /T ∥ )(s) = 4, so as to model perpendicular pressure dominant plasmas. We have modeled perpendicular pressure dominant plasmas given that ANIMEC does not produce converged magnetic equilibria for parallel pressure dominant plasmas at < β > ≈ 3% and for < β h > / < β >> 1/2. When modeling parallel dominant plasmas, corresponding to a hot particle temperature ratio of (T ⊥ /T ∥ )(s) = 1/4, the parallel pressure and its gradient increase in magnitude for increasing < β h > / < β >. Large parallel pressures result in a large plasma energy as shown in (2), which makes it increasingly difficult for ANIMEC to produce a self-consistent equilibrium. For < β h > / < β >> 1/2 a self-consistent equilibrium cannot be produced, irrespective of the size of the time step. For lower ratios of < β h > / < β > the parallel pressure dominant case shows similar behavior as the perpendicular pressure dominant case. For this parameter scan, TERPSICHORE is used to analyze the n = 1 mode family [13][14][15]. This mode family has been chosen as subject for this analysis, given that the most dominant MHD instabilities in LHD experiments belong to the n = 1 mode family [6].

Simulation results
As an example, the flux surface-averaged pressure profile of the thermal particles as well as the profile of the parallel and perpendicular components of the hot particle pressure, as calculated by ANIMEC, are shown in figure 2 for The flux surface-averaged thermal pressure profile p th and the parallel and perpendicular parts of the hot particle pressure profile p h∥ and p h⊥ , respectively. These profiles have been calculated by ANIMEC for a pressure ratio of < β h > / < β > ≈ 1/2. Table 1. The radial location of the ι = 0.5 surface, the n = 1, m = 2 mode peak, their difference in location ∆s, the pressure gradient at the iota surface and the eigenvalue for each simulation. < β h > / < β > ≈ 1/2. In figure 3 the mode structure of the plasma displacement ξ n,m (s) corresponding to the five most unstable mode numbers is shown for the pressure ratios, < β h > / < β > = 0, 1/10, 1/3, 1/2, 2/3, 9/10, respectively. The rotational transform ι has been plotted to identify the surfaces at which mode resonances occur. The ι = 0.5 surface is shown to indicate the radial location of the n = 1, m = 2 mode number combination. At this flux surface, a resonant mode is expected and has been observed for all pressure ratios. Figure 3 shows, however, that the peak of the resonant mode does not match exactly the location of the ι = 0.5 surface. We denote the location of the ι = 0.5 surface by s ι , the location of the n = 1, m = 2 mode peak by s p and we define the difference between these values ∆s ≡ s p − s ι . These quantities are presented, for each simulation, in table 1. Additionally, the pressure gradient at the location of the resonant ι = 0.5 surface and the growth rates pertaining to each simulation are given. Observing figure 3, it is apparent that for < β h > / < β > ⩽ 1/3 the n = 1, m = 2 mode numbers represent the dominant instability, while the n = 1, m = 1 mode number becomes dominant for < β h > / < β > ⩾ 1/2. This change in behavior is most likely due to the fact that the pressure gradient at the ι = 0.5 surface changes sign when the hot pressure ratio in increased from < β h > / < β > = 1/3 to < β h > / < β > = 1/2, which has a stabilizing effect on the n = 1, m = 2 mode number. It should also be noted that for < β h > / < β > ⩾ 2/3, corresponding to figures 3(e) and (f ), the plasma displacement pertaining to the n = 1, m = 2 mode number is small in amplitude compared to the dominant mode number. This implies that the physical relevance of the last two rows in table 1 is limited and will therefore be ignored in the following discussion.
Both figure 3 and table 1 show that the introduction of hot particles to the plasma has a profound effect on the rotational transform. For the top four rows of table 1, we observe that the difference ∆s in location between the mode peak and the ι = 0.5 surface, increases monotonically for increasing < β h > / < β >. It is worth pointing out that this is true, despite the fact that s ι and s p have a more complicated dependence on < β h > / < β >. These simulations show that ∆s increases significantly for relatively low hot particle pressures of < β h > / < β >⩽ 1/3, which have been shown to be experimentally relevant in the case of heating through NBI [6]. Experimental results for the LHD obtained by Motojima et al, have shown that for an NBI heated plasma with high β, the n = 1, m = 2 resonant mode disappeared when increasing β [1]. Similar behavior is found in figure 3 when increasing the hot pressure ratio from < β h > / < β > = 1/3 to < β h > / < β > = 1/2.  Figures (a) through (f ) correspond to the following ratios of plasma pressure < β h > / < β > = 0, 1/10, 1/3, 1/2, 2/3, 9/10, respectively. The rotational transform profile is shown in red. The red dashed lines indicate the flux surface coordinate s where the ι = 0.5 surface is located.

Conclusion
In conclusion, in this work we have performed an analysis on the plasma displacement in an LHD plasma with an anisotropy in pressure. We have used the ANIMEC code to calculate magnetic equilibria for plasmas with a pressure anisotropy. Furthermore, the ideal MHD stability code TERPSICHORE was used to calculate the plasma displacement. Using a pressure distribution that models off-axis heating through, for example, NBI, we have scanned the ratio < β h > / < β > in order to analyze the mode structure of unstable modes. The results show that the introduction of hot particles to the plasma has a profound effect on the mode structure as well as on the rotational transform profile. Most notably, we have found that the difference in location between the peak of the n = 1, m = 2 mode and the ι = 0.5 surface increases monotonically for increasing < β h > / < β >. These results are applicable in experimentally relevant regimes and can therefore have a significant impact on MHD mode suppression techniques.