Drift kinetic effects on plasma response in high beta spherical tokamak experiments

The high β plasma response to rotating n=1 external magnetic perturbations is numerically studied and compared with the National Spherical Torus Experiment (NSTX). The hybrid magnetohydrodynamic(MHD)-kinetic modeling shows that drift kinetic effects are important in resolving the disagreement of plasma response between the ideal MHD prediction and the NSTX experimental observation when plasma pressure reaches and exceeds the no-wall limit (Troyon et al 1984 Plasma Phys. Control. Fusion 26 209). Since the external rotating fields and high plasma rotation are presented in the NSTX experiments, the importance of the resistive wall effect and plasma rotation in determining the plasma response is also identified, where the resistive wall suppresses the plasma response through the wall eddy current. The inertial energy due to plasma rotation destabilizes the plasma. The complexity of the plasma response in this study indicates that MHD modeling, including comprehensive physics, e.g. the drift kinetic effects, resistive wall and plasma rotation, are essential in order to reliably predict the plasma behavior in a high beta spherical tokamak device.


Introduction
Two-dimensional (2D) tokamak plasmas, having toroidal symmetry, can be modified by non-axisymmetric (3D) perturbed magnetic fields which are generated by external coils. 3D fields cause a plasma response, producing perturbed plasma current, plasma displacement and magnetic perturbations [1][2][3]. This phenomenon has been extensively studied in various tokamak devices, including DIII-D [5][6][7][8], JET [9], NSTX [10,11], etc. The plasma response is important in predicting the utility of 3D fields, such as the detection of plasma stability with resonant field amplification and control of the edge localized mode using 3D magnetic coils [13,14]. The plasma response can also be utilized to modulate the plasma rotation through the neoclasical toroidal viscosity (NTV) [15][16][17][18]21]. Therefore, understanding the physics and making a reliable prediction of the plasma response can be essential in order to achieve the high performance operation of present tokamak devices and ITER [12].
Previously, both NSTX and DIII-D experiments [5,10] have shown that the amplitude of n = 1 plasma response linearly increases with the normalized pressure β N across the ideal kink stability limit β no−wall N , where n is the toroidal mode number, β N = β(%)/[I p (MA)/a(m)B 0 (T)]. β no−wall N is predicted by ideal MHD modeling [19]. Here β = 2µ 0 p /B 2 0 , p is the volume-averaged plasma pressure, B 0 is the magnetic strength at the plasma center, I p is the plasma current, a is the plasma minor radius, and µ 0 is the magnetic permeability. The conventional plasma response modeling, based on

Nuclear Fusion
Drift kinetic effects on plasma response in high beta spherical tokamak experiments ideal MHD theory, predicts that the plasma response has a peaked amplification at β no−wall N which is inconsistent with the nearly linear increment of the plasma response observed in the experiments. The ideal MHD also fails to predict the plasma stability since the plasmas remain stable in NSTX and DIII-D experiments when β N > β no−wall N . In the DIII-D plasma response experiments, the simulation of plasma response, based on hybrid MHD-kinetic modeling, resolves this inconsistency and shows excellent agreement with the experimental measurements. The study reveals the importance of the drift kinetic effects, produced by thermal ions, in determining the plasma response [20].
In this work, hybrid MHD-kinetic modeling is further applied to understand the observed n = 1 plasma response in the NSTX experiments. Similar to the DIII-D experiments, it is found that the kinetic effects also play an essential role in determining the plasma response in NSTX. A key difference is that the high plasma rotation makes the underlying kinetic contrib ution different from the DIII-D experiments. We also show that the resistive wall effects and plasma inertia are important due to the rotating external field and high plasma flow. The results highlight that a comprehensive modeling of the plasma response is required to predict the plasma behavior in NSTX. Meanwhile, this work further confirms the validity of hybrid MHD-kinetic modeling, which is essential for making a reliable prediction of plasma response in high β tokamak operations. Section 2 introduces the details of the plasma response model. In section 3, the computed n = 1 plasma response in the NSTX experiments is reported and analyzed. Section 4 summarizes this work.

Fluid plasma response formulation
Following previous studies [5,10], the fluid plasma response modeled by the single fluid ideal MHD equations, ignoring the plasma flow, can be written in terms of various perturbed quantities, such as plasma displacement ξ , perturbed velocity v, perturbed plasma current j , perturbed magnetic field b, and perturbed (fluid) pressure p. The formulation in the full toroidal geometry is described by the following linearized ideal MHD equations in the Eulerian frame: In order to perturb the plasma, the external magnetic fields, which are generated by a source current j coil flowing in the external coils located in the vacuum region, satisfy: In the above formulation, ω = if , f is the rotating frequency of the external fields generated by the coils. The equilibrium quantities of plasma density, magnetic fields, current, and plasma pressure are defined as ρ, B, J ,P respectively. A conventional unit system is assumed with vacuum permeability µ 0 = 1. The ratio of specific heats is Γ = 5/3. The external magnetic perturbations from the coils are assumed to have an e inφ dependence along the toroidal angle φ. Equations (1) and (2) are solved by the MARS-F code [22] to simulate the fluid plasma response.

Kinetic plasma response formulation
The kinetic plasma response requires that the drift kinetic effects are coupled with global MHD behavior, where the formulation of the drift kinetic effects is derived from the perturbed drift kinetic theory and associated with distorted particle orbits by 3D fields [23][24][25]. As indicated in [20,27], by extending the MHD equation, the following hybrid MHDdrift kinetic formulation, implemented in the MARS-K code [4,26], needs to be solved with a self-consistent approach in order to enable interaction between the plasma response and the drift kinetic effects: where p is the anisotropic perturbed pressure tensor, b = B/B is the unit vector of the equilibrium magnetic field, and B is the strength of the equilibrium field. The parallel and perpendicular components of the kinetic pressure p and p ⊥ are defined by (R, φ, Z) which represent the cylindrical coordinate system for the torus. φ is the unit vector in the toroidal direction. Z is the unit vector in the vertical direction in the poloidal plane. V 0 = RΩφ is the plasma equilibrium flow, where Ω is the angular frequency of plasma rotation in the toroidal direction. On the right-hand side of equation (3b), the fourth and fifth terms represent the Coriolis force and the centrifugal force. p and p ⊥ , derived from the drift kinetic theory [24,25], are defined as In equations (4)-(6), the parallel and perpendicular velocity components of the particles are denoted by v and v ⊥ . σ = sign(v ). M α is the particle mass of ions or electrons. The perturbed particle distribution function f 1 L is derived by solving the bounce-averaged perturbed drift kinetic equation [25] for each specie of particles. The summation in equations (4) and (5) can be over the components provided by thermal electrons, thermal ions and energetic particles. The integral is carried out over the velocity space of trapped and passing particles, where the core term, f 1 L , is written as Here, f 0 ε is the energy derivative of the particle equilibrium distribution function. φ (t) = φ(t)− <φ > t, where < · > denotes the average over the particle bounce period. m corresponds to the Fourier harmonic number along the poloidal angle. X m and H ml , defined in [4], are related to the perturbed particle Lagrangian H L where κ = ( b · ∇) b is the magnetic curvature, and µ = Mv 2 ⊥ /2B is the particle magnetic moment. The factor λ l represents the mode-particle resonance operator, where ω * N and ω * T are the diamagnetic drift frequencies associated with the plasma density and temperature gradients, respectively. Details of the hybrid MHD-kinetic modeling can be found in [4]. MARS-K and MARS-F codes share the same module in generating the external magnetic perturbation. While solving the fluid and kinetic plasma response, the vacuum field equations outside the plasma, the thin resistive wall equation, and the coil equations [2] are solved together with the aforementioned MHD equations in MARS-F/K. It is noted that equation (1) can physically recover the model of ideal perturbed equilibrium [1] while ω → 0.

Configuration of NSTX equilibrium and coils
The modeling of ideal and kinetic plasma response introduced in section 2 is applied to understand the observation of the NSTX plasma response experiment presented in [10]. In this experiment, the RWM-EF coils, located at the middle plane of the low field side (LFS), are used to apply the rotating external n = 1 magnetic perturbations with f coil = ±30 Hz field rotating frequencies, where '+' and '−' denote the cocurrent and counter-current direction of the plasma. A pair of upper magnetic sensors, located at LFS, is used to measure the radial perturbed fields of n = 1 plasma response in this study, where the locations of two sensors have a 180 degree difference in the toroidal direction. The 2D schematic of the coil and sensor geometry is illustrated in figure 1 and implemented in the MARS-F/K code. The effect of the resistive wall is also included in the simulation to effectively consider both vacuum vessel and passive plates. There are various ways to define the plasma response based on the magnetic sensor measurements [10,14,20], where these definitions of plasma response can be adapted to each other. Here, the normalized perturbed quantity, δ B tot /δ B vac , the same as in [10], is used to represent the plasma response, and will be compared between the experiment and MARS-F/K simulations, where δ B tot is the total perturbed magnetic field measured by the sensor in the plasma response experiments, δ B vac is the measured perturbations in the presence of the RWM-EF coils only in vacuum. Figure 2 shows the time evolution of the plasma current, q 95 , β N and the measured δ B tot in the NSTX discharge 124 801. Since β N increases continuously in the NSTX experiment, the plasma response has to be extracted in a short time interval to ensure β N changes little. The details of the experimental method to measure the β N dependence of the plasma response are described in [10]. Moreover, the strong plasma rotation and the resistive wall in the NSTX experiments can shield the applied ±30 Hz rotating fields. These factors make the extracted plasma response noisy in the signal analysis. Therefore, using δ B tot to describe the plasma response is preferred in this study since δ B tot is more tolerant of the measurement error than the pure plasma In the NSTX experiments, it is also found that the upper radial sensors, located at LFS, measure the magnetic response with a lower signal noise than other sensors. Using this sensor array to analyze the plasma response can further reduce the error bar of measurement, showing a clear β dependence of the plasma response, δ B tot , for a quantitative comparison with the MARS-F/K simulation while studying the NSTX plasma response. The employed equilibria in the simulation are generated from NSTX discharges (124 808 at 522 ms, 124 801 at 502 ms, 560 ms, and 610 ms), using the LRDFIT code, where the code described in [26] reconstructs the NSTX equilibrium, including plasma rotation, kinetic profiles, etc. The β N values of the reconstructed equilibria vary from 3.64-4.98. MARS-F, based on the ideal MHD simulation, predicts the no-wall beta limit β no−wall N = 4.75. The safety factor, pressure, and current profiles for these cases are plotted in figure 3.

Comparison of fluid/kinetic plasma response with NSTX experiments
To verify the ideal and kinetic plasma response models, the experimental plasma response measured by the upper radial sensor array at LFS is compared with the MARS-F/K simulation. A comparison of the amplitude and toroidal phase of the magnetic response δB tot /δB vac is plotted as a function of figure 4. Two 'fluid' cases in terms of fluid plasma responses with/without the resistive wall effect are presented here. The 'thermal' case corresponding to the kinetic plasma response is also reported in the figure, where the adiabatic contributions from both thermal particles (TPs) and energetic particles (EPs) are included, but the non-adiabatic term includes only the TP contribution in this case. As shown in the figure, each 'fluid' case has two more points than the 'thermal' case. The corresponding equilibria of the two points are obtained by scaling down the pressure of the equilibrium at β N = 0.98β no−wall N , since the 'fluid' cases do not require the kinetic profiles in the simulation. Figure 4 shows that the experimental plasma response has a linear increment and no phase jump while β N increases above β no−wall N . In general, the 'fluid' cases, with/without the resistive wall, are in poor agreement with the experimental measurements for both the amplitude and phase of response, where the no-wall 'fluid' case shows a large amplification and a significant phase change near β no−wall N . The 'fluid' case with the wall clearly shows a significant difference in response amplitude when f = −30 Hz. In contrast, the kinetic 'thermal' case is in reasonably good agreement. It reproduces a linear increment and a smooth phase variation in the experimental measurements. In this comparison, many channels can lead to noise in the plasma response, as discussed in section 3.1. The continuously varied plasma pressure, the shielding effect of the strong plasma rotation and resistive wall can lead to a  large error bar of response in the experiments. The pedestal of plasma pressure in the equilibrium reconstruction is also not well resolved in figure 3(b). Nonetheless, a quantitative agreement between the simulated kinetic response and the experimental measurement is still obtained. The dominant physics of the kinetic thermal particles will be discussed in detail later. When β N approaches and exceeds the limit of β no−wall N , the amplitude of the fluid response without a resistive wall has a peaked amplification due to the vanishing perturbed potential energy δW = δW p + δW vac → 0, where δW p is the plasma potential energy, δW vac is the vacuum energy. Therefore, the plasma is free to be perturbed by external fields. Moreover, this no-wall 'fluid' case has approximately a 180 degree phase change since the sign of δW switches from positive to negative in figures 4(b) and (d), as indicated in [20]. On the other hand, we notice that the amplitude of the plasma response is suppressed in the 'fluid' case with a resistive wall. Since the applied field frequency is |f | = 30 Hz and the resistive wall   time of the NSTX device is τ w ∼ 3.5 ms, the studied NSTX experiments have 2π|f | ∼ 1/τ w , which means that the coil frequency is close to the field penetration frequency of the wall. The applied rotating fields cause an eddy current in the resistive wall which can strongly suppress the plasma response. Due to the dissipative effect of the eddy current, the phase of the fluid response is also adjusted. In the no-wall fluid case, the amplitude of the plasma response is symmetric for f = ±30 Hz. Compared with the no-wall fluid response, the existence of the resistive wall causes the amplitude of the fluid response to become non-symmetric when f = ±30 Hz. But the symmetric property of the linearized MHD equation still remains. This wall effect is discussed further in section 3.3. Finally, as shown in 4, only the 'thermal' kinetic response agrees with the experimental measurements in both f = ±30 Hz cases. In particular, the kinetic MHD simulation, consistent with the experimental observation, predicts that the plasma is stabilized by the kinetic effects when β N > β no−wall N . However, ideal MHD predicts that the plasma, even with a resistive wall, is unstable. This stability prediction is shown in the later Nyquist analysis. Note that, though the steady-state ideal response loses physical meaning due to RWM instability, MARS-F can still compute such a response by inverting the system matrix directly.
MARS-K simulation further reveals that the thermal particles contribute to the predominant kinetic effects in modifying the fluid plasma response and reproducing the experimental measurements. Figure 5 compares four cases of simulated plasma response which all include a resistive wall: 'fluid', 'thermal', 'thermal+fast' and 'full thermal', where the 'full thermal' case assumes that all equilibrium pressure comes from TPs. The kinetic response of the 'thermal+fast' case means that the non-adiabatic contributions from the EPs are added on top of the 'thermal' case. The EPs are modeled with an isotropic slowing down distribution, with fast ion pressure and density computed by the TRANSP code. The 'thermal+fast' case in figure 5 shows similar behavior to other thermal cases, which indicates that the kinetic effects, caused by TPs, play a dominant role in modifying the fluid response leading to the kinetic response. The 'full thermal' case has a slightly larger amplification than the other two kinetic response cases, since one surface term in the adiabatic EP pressure has a damping effect. This term is due to the integration by parts in the particle phase space for the slowing down EPs with finite birth energy [28].
As for the thermal kinetic effects, various kinetic effects can be produced by both thermal ions and thermal electrons from the non-adiabatic kinetic pressure tensor. For instance, both trapped ions and electrons contribute to the kinetic energy through toroidal precession resonance and bounce resonance. The passing ions and electrons provide the kinetic energy through transit resonance. To better understand which kinetic effect is important to the NSTX kinetic response, figure 6 presents the real and imaginary parts of δW K in terms of different kinetic resonances. These energy components are analyzed with the coil frequency f = 30 Hz for the discharge 124 801 at 560 ms which is near the marginal stability β N /β no−wall N = 0.98. The kinetic energy δW K is normalized by the volume integrated inertia δK = ρ(ξ · ∇s) 2 dV . Clearly, thermal ions contribute much more kinetic energy than thermal electrons, since thermal electrons have much higher collision, bounce and transit frequencies than the thermal ions. On the other hand, the imaginary parts of δW K are mainly due to the ion bounce and transit resonances. The reason is that ω E is on the same order of < ω i b > and < ω i t > as shown in figure 7. It means the EXB rotation can match with lω i b and (m − nq + l)ω i t in the phase space, resulting in a strong kinetic resonance and a large imaginary part of δW K . In contrast, the kinetic energy terms studied in the DIII-D case [20] have larger real parts than imaginary parts. This can be understood with the resonant operator (9) at the limit of f → 0, which gives Strong resonance can occur and maintain the finite imaginary part, and the real part of equation (10) approaches to 0, which is the case for NSTX as presented here because ω E is large enough to make b → 0. Therefore, the imaginary part of δW K can be dominant. We have discussed the importance of E × B rotation in determining the kinetic response through kinetic energy. Moreover, the rotational effect can also affect the plasma response through the plasma rotation frequency Ω in the MHD equation (3). To investigate the impact of plasma rotation on plasma response, the 'thermal+fast' case, including plasma rotation, is compared with the same case with zero plasma rotation by assuming Ω = 0 in the MHD part. Figure 8 shows that the kinetic response with plasma rotation (dashed line) has a stronger amplification than the case without flow when β approaches and exceeds β no−wall N . It implies that the plasma rotation in the MHD part plays a role in destabilizing plasma. An energy analysis is performed to further understand what physics causes the destabilizing effect with plasma rotation. In the momentum equation (3b), the plasma rotation enters the left-hand side term and the last two terms on the right-hand side. Therefore, plasma rotation contributes three energy components: plasma inertia, Coriolis force and centrifugal force respectively. Figure 9 compares the three energy components in terms of the 'Thermal+Fast' case near β no−wall N . It clearly shows that the energy provided by the centrifugal force is much smaller than the other two terms. The energy of the Coriolis force has a positive real part which can stabilize the plasma. The imaginary part of this energy, providing the dissipation effect, has an opposite sign to the imaginary part of δW K . It indicates a cancellation of dissipative effects between the Coriolis force and δW K . The plasma inertial energy, as the predominant term, has a negative value which eventually drives the plasma toward marginal stability and results in an amplification of the kinetic response, as shown in figure 8.
In general, the energy analysis indicates that the net contribution, among all the kinetic energy components of the thermal ions and rotation-related energy terms, leads to a modification of the plasma response consistent with the experimental measurements.
With the above analysis, the kinetic plasma response explains the linear increment of the NSTX plasma response while approaching and exceeding β no−wall N . The kinetic plasma response modeling is also applied to understand the frequency dependence of n = 1 plasma response in the NSTX experiments, which is presented in figure 2 of [29]. The measurements show that the maximum amplification of n = 1 response appears between 20-30 Hz while scanning the coil frequency in the experiments. The results imply that a non-negligible dissipative effect of plasma should exist to cause the frequency shift. Here, the frequency response is simulated by using the equilibrium (discharge 124 801 at 560 ms) studied in  this paper. Figure 10 shows that the fluid response has a peak almost at f = 0 Hz. On the other hand, the simulation of the kinetic response qualitatively reproduces the frequency shift where the largest amplification appears at 20 Hz, as shown in [29]. It indicates that kinetic dissipation can be an important candidate in the frequency response. The results also show the importance of hybrid MHD-kinetic modeling in predicting reliable plasma behavior.

Nyquist analysis of NSTX plasmas
In order to understand the change in the plasma stability by the kinetic effects and the plasma flow while β N > β no−wall N in the concerned NSTX experiment, a sequence of frequency scans with RWM coils is performed by MARS-K/F to carry out the Nyquist analysis. The way to form the Nyquist contour is to vary the field rotating frequency, f, from −∞ to +∞. The real and imaginary parts of the radial magnetic perturbations measured by a simulated upper sensor can be plotted in the complex plane to form the Nyquist contour with the coil frequency dependence which is shown in figure 11 where the equilibrium with β N /β no−wall N = 1.05 is used. The points in the figure marked by 'o' and ' ' correspond to the response with f = −2 Hz and f = 2 Hz. The figure clearly shows that the kinetic effects not only change the shape of the Nyquist contour for the fluid response case, but also change the direction of the contour from counter clockwise to clockwise. Here, the 'full thermal' case of the kinetic plasma response is considered in the kinetic response simulation. Comparing the two kinetic Nyquist contours, the plasma flow in the MHD part also has an impact on the change in the Nyquist contour. In the other words, the differences in the Nyquist contours imply a change in plasma stability. Since the plasma response in linear theory is the result of the linear combination of several eigenmodes of plasma responding to an external field, the Padé approximation, based on the linear transfer function (11), can approximate the Nyquist contour to extract the dominant eigenmode [22]: Here, each term in the transfer function corresponds to one eigenmode. γ j is the eigenvalue of the jth eigenmode where the real part of γ j represents the damping/growth rate of the eigenmode, R j is the residual (coefficient) of the eigenmode. The fitted transfer function can be used to infer the plasma stability. Applying equation (11) to fit the three contours in figure 11, the transfer functions extracted from the contour of the 'fluid' case find one unstable eigenmode with Re(γ)= 13.66 Hz >0, which means the plasma in the fluid MHD approach, even with the resistive wall, is unstable since β N > β no−wall N . The least stable modes extracted from the kinetic contours without and with rotation, give Re(γ)= −34. 3 Hz and −13.22 Hz respectively, which means the kinetic effects significantly stabilize the plasma. Comparing the two kinetic cases, it also indicates that the plasma rotation in the MHD equations slightly destabilizes the plasma, which is consistent with the energy analysis of the rotational effect in figure 9.
In figure 4, the fluid response without the resistive wall is symmetric when f = ±30 Hz, since the linearized MHD equation (1) is hermitian. Moreover, though the fluid response with the resistive wall has a different amplitude and phase at the two frequencies, the linear plasma response still retains a symmetric property. Figure 12 assumes two sensors which possess an up and down symmetry at LFS. The Nyquist contour of the plasma response corresponding to each sensor is still symmetric. A similar symmetric property can be found in [31]. In figure 4, a multi-mode response, as a potential candidate, can cause the non-symmetric frequency dependence of the fluid response with a wall. The existence of a multi-mode  fluid response [14] is also indicated in figure 12, because the single mode response should result in a pure circular Nyquist contour. Here, the multi-mode response may be induced by coupling between the plasma response and resistive wall which both respond to the rotating perturbed fields generated by RWM-EF coils. Considering the measured magnetic perturbations on the assumed upper and lower sensors, our study finds that the fluid plasma response with no wall is complex conjugate on the two simulated sensors, and the resistive wall generates an up and down symmetric response. The coupling of the two patterns could produce and enhance the multi-mode response. This numerical simulation eventually points out the complexity of the plasma response while including the resistive wall. In the future, the resistive wall effect should be further studied to help us better understand plasma response in the presence of a rotating field. However, this work indicates that the kinetic effects are dominant in determining the plasma response to reproduce the experimental measurements. Hybrid MHD-kinetic modelling is essential in order to predict reliable plasma behavior in high NSTX plasmas.

Summary
The simulated fluid and kinetic plasma response by MARS-F/K is scrutinized in order to understand the physics of the measured n = 1 plasma response in the NSTX experiments [10]. Similar to the DIII-D experiments [5], the fluid plasma modeling fails to predict the experimental magnetic response while β approaches or exceeds the ideal no wall β limit. As indicated in [20], the kinetic effects can reconcile the disagreement between the predication of plasma response and the NSTX experimental observation. On the other hand, this work reveals that the dissipative term, predicted by [10] in the NSTX experiments, is produced by the kinetic effects, par ticularly from thermal particles. The simulation also finds many unique physical features of plasma response in the NSTX experiments which are different from the DIII-D experiment. Since the high EXB frequency can well match the bounce frequency of trapped ions and the transient frequency of passing ions in the NSTX experiments, the imaginary part of δW K can be dominant and is caused by the bounce resonance and transient resonance of thermal ions, whilst the real part of δW K is dominant in the DIII-D experiments [20]. Based on the equivalence between δW K and NTV torque [18,30], this result implies that the NTV torque produced by these two resonant effects could be important. We also find that the eddy current effect largely suppresses the NSTX plasma response since the ±30 Hz fields' rotating frequency is comparable with the field penetration frequency of the NSTX resistive wall. Moreover, the fast plasma flow in the MHD part plays a destabilizing role to amplify the plasma response. Finally, the aforementioned physical effects are all essential ingredients in determining the NSTX plasma response. In the Nyquist stability analysis, the hybrid MHD-kinetic modeling predicts that the plasma is stable when β exceeds the ideal no wall limit. This is consistent with the NSTX experimental observations. Complementing the validation of the kinetic plasma response in the DIII-D experiments [20], the study of the NSTX plasma response further confirms the validity of the kinetic plasma response modeling. The results demonstrate the importance of solving the hybrid MHD-kinetic equations in order to reliably predict plasma response behavior and plasma stabilities in high performance tokamaks.