Identification of multiple eigenmode growth rates towards real time detection in DIII-D and KSTAR tokamak plasmas

The successful application of three-dimensional (3D) magnetohydrodynamic (MHD) spectroscopy enables us to identify the multi-mode eigenvalues in DIII-D and KSTAR tokamak experiments with stable plasmas. The temporal evolution of the multi-modes’ stabilities have been detected. The new method is numerically efficient allowing the real time detection of MHD modes’ stabilities during the discharge. The method performs active detection of the plasma stability by utilizing the upper and lower rows of internal non-axisymmetric coils to apply a wide variety of 3D fields. Multi-mode eigenvalues are extracted using subspace system identification of the plasma response measured by 3D-field magnetic sensors distributed at different poloidal locations. The equivalence of this new method with the one introduced by Wang (2019 Nucl. Fusion 59 024001) has been numerically corroborated. The more robust and efficient calculation developed here will enable real time monitoring of the plasma stability based on the extracted eigenvalues of stable modes.

* Author to whom any correspondence should be addressed.
Original content from this work may be used under the terms of the Creative Commons Attribution 3.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. dangerous to advanced tokamak operation, such as ITER-like tokamaks and future high power fusion reactors [4,5], where the unstable global MHD modes can disrupt the plasma and lead to severe damage of tokamak devices. Therefore, it is essential to develop a technique dedicated to the real time stability detection of these MHD modes with quantitative stability measurements to track and avoid instabilities that lead to disruptions.
Tokamak plasmas can be sensitive to a low level of nonaxisymmetric magnetic perturbations applied by external coils [6], and the response of plasmas contains information related to the stability of various MHD modes. For instance, the resonant field amplification of plasma response, indicating the β limit of ideal kink instability, has been extensively investigated in various experiments [7][8][9][10]. Such responses are often modelled using linear MHD theory [11][12][13][14][15][16], since the magnitude of magnetic perturbation (δB) is often several orders smaller than the equilibrium magnetic field (B). Many previous works [17][18][19] assume one single dominant stable eigenmode in the plasma response. However, the plasma response can contain multiple significant stable MHD eigenmodes, which is not only indicated by theoretical studies [20][21][22][23][24][25][26] but also observed in the EAST [27] and DIII-D experiments [28][29][30]. Therefore, the plasma response is the linear combination of multiple stable eigenmodes in the linear system [29,30].
Theoretically, the active detection of dominant eigenmodes in plasma can be achieved by varying the external perturbations, which was first developed in [6] called MHD spectroscopy, to identify the damping rate and the mode rotating frequency of the stable resistive wall mode. In the method, the single-mode assumption is adopted as done in many previous works [17][18][19]. Since the multi-mode plasma response has been detected in DIII-D experiments [28], it motivates us to develop multi-mode three-dimensional (3D) MHD spectroscopy based on the multi-pole transfer function for more reliable detection of plasma stability [31]. Since the transfer function needs to be extracted by fitting the measured magnetic response and perturbed coil currents in the frequency domain, this method is referred to as the frequency domain method (FDM) in this work. The extracted transfer function not only reveals the contribution of each dominant eigenmode in the plasma response, but also quantitatively indicates the stability of each dominant mode through the extracted eigenvalues which can provide a critical assessment of MHD stability of plasma. Although the FDM has been applied in experiments, the eigenmode extraction was performed offline in post-analysis due to the complicated frequency analysis of signal and nonlinear fitting of transfer function. The FDM method is not efficient for experimental real time detection.
A new time domain method (TDM), working in time domain directly, is developed by employing the subspace system identification (SSI) theory [32], which has previously been utilized to understand MHD stabilities in reversed field pinch and tokamak plasmas through offline analysis, respectively [33,34]. The TDM analysis does not require complicated frequency analysis and nonlinear fitting in the FDM analysis. In the TDM, the linear least square fitting, without requiring the initial guess, not only guarantees the numerical convergence of extracted eigenmodes, but also achieves impressive numerical efficiency. Previous real time spectroscopy [35] used a single-mode model approximation. In this work, it shows for the first time the feasibility of quantitatively detecting the multiple MHD modes' stabilities in real time, based on the active detection and the TDM analysis, which is so-called real time 3D MHD spectroscopy. The technique can be important to monitor the evolution of plasma stability and to avoid the severe plasma disruption in the advanced operation of a fusion reactor.

Real time 3D MHD spectroscopy
Real time 3D MHD spectroscopy is approached by involving SSI theory [32,34]. Compared with the FDM, which requires the nonlinear fitting to extract the transfer function, SSI improves the numerical efficiency dramatically while extracting the mode stabilities with linear least square fitting. SSI solves the following linear state-space response model (1), which can represent the plasma response modelled by linear  MHD theory in the presence of external 3D perturbations.
( 1 a ) Here, x k represents the system state of the process at discrete time instant k and contains the numerical values of N states. δB k is the system output denoted as the total magnetic response including the contributions of plasma, vacuum vessel and other components in the system. δB k is measured by 3D magnetic sensor arrays at different poloidal locations. The vector δJ k , with dimension M, is the system input of 3D coil current, where M is the number of employed 3D coils. Equation (1a) represents the coupling between the coil current and eigenmodes, and equation (1b) represents the measurements of the plasma response. A, B, C, and D are system matrices to be determined. The N by N matrix A describes the dynamics of the system, where N denotes the number of dominant eigenmodes in the system. In our analysis, the dominant eigenmodes mean the eigenmodes play a dominant role and make major contributions to the plasma response. Here, the eigenmodes are less stable individual modes [31] and can be resonant with the external perturbation due to the small eigenvalues [17]. Matrix B, representing the impact of the applied 3D coil current δJ k , has N by M dimensions. The l by N matrix C, with l being the number of employed sensor arrays, describes the contribution of eigenmodes' response to the measurements δB k . Matrix D, having l by M dimensions, represents the response of vacuum vessel and other components to the applied 3D perturbation. Since this model can also be directly derived from MHD equations, the eigenmodes extracted from this model should be the same as those extracted from the multi-mode transfer function [31]. In other words, these two methods are theoretically equivalent. Evidence of the equivalence between these two methods will be shown later in this work.
Although this work will not focus on the algorithm itself, here the process of the eigenmodes extraction is briefly introduced as follows. The input and output data can be obtained by conducting a dedicated set of system identification experiments. Different from FDM, TDM estimates the system state x k sequence by making a series of mathematical operations on input and output data. Therefore, one can identify the system matrices A, B, C, and D by linear least square fitting, which is the major advantage of this TDM to greatly improve the fitting efficiency and convergence. Finally, the eigenvalues of the eigenmodes, as the stability index of mode, can be computed by the eigenvalue decomposition of the matrix A in TDM analysis. Note that a finite number of x k is used to determine matrix A, which evolves on a longer time scale than the k step.

Experimental analysis in DIII-D tokamak
The TDM has been used to detect the time evolution of the growth rate for multiple stable MHD modes in DIII-D plasmas. In these experiments, the upper and lower rows of internal 3D coils [9] are applied to generate n = 1 3D magnetic perturbations. The number of system input, the coil current δJ k , is thus M = 2. Both the coil current phasing Δφ and rotation frequency f are scanned to measure the variation of the plasma response. This provides more information of system output δB k for better fitting equation (1). Here, the coil phasing Δφ = φ up − φ low is the difference between the upper coil current phase, φ up , and the lower coil current phase, φ low . In addition, both the square wave and travelling wave coil oscillations were applied to assess any difference in the reliability between the two which can help to optimize the coil waveform in the future development of real time 3D MHD spectroscopy. In this paper, the plasma response is denoted as the total magnetic perturbations measured by the toroidally distributed 3D magnetic sensor arrays located at the middle plane on both the low field side (LFS) and high field side (HFS) to observe more aspects of plasma response. The configuration of magnetic probes are described in [25,36]. coils is 2 kA. Since the ratio of the magnetic response to the equilibrium toroidal field is less than 0.1%, figure 1(a) indicates the equilibrium with little variation is hardly perturbed by the applied 3D fields. Therefore, the plasma response is in the linear regime in the DIII-D experiments, which has already been validated in [31]. In the stable discharge 178 583, the coil current, illustrated in figure 1(b), is applied to generate n = 1 perturbed fields in the flattop phase, where figure 1(a) shows plasma current (I p ), safety factor (q 95 ) and normalized beta (β N ) have little variation. Figures 1(c) and (d) present the time evolution of n = 1 magnetic response measured by the radial sensors located at the middle plane of LFS and HFS, respectively.
With the input of coil current (δJ k ) and the output of measured 3D magnetic response (δB k ), TDM analysis is made by fitting equation (1) with SSI method [32,34]. Figures 1(c) and (d) report a good agreement between the fitting results and DIII-D experimental data for both HFS and LFS magnetic response. By fitting the time interval from 1.1 s to 2 s, TDM finds two dominant stable eigenmodes, where the least stable mode has the eigenvalue γ 1 = −9.54 + 1.88i Hz and the secondary mode has γ 2 = −73.24 + 0.85i Hz. It is noted that the real part of the eigenvalue, Re(γ), is identified as the damping rate of mode, where the more negative damping rate indicates the mode is more stable. The imaginary part of eigenvalue represents the natural mode rotating frequency. Though there are only stable modes found in this stable plasma, it is important to note that the SSI method, employed by TDM analysis, is theoretically capable for the detection of unstable modes.
For the purpose of developing real time 3D MHD spectroscopy to detect the eigenmodes' stabilities in real time, equation (1) needs to be fit using a time sequence of experimental data within a short time window. In figure 2, equation (1) is fit in a short time window ∼100 ms at both the beginning of the travelling and square waves covered by the shade respectively. The predictions, made by the travelling wave fit model and the square wave fit model, are represented by the dash-dotted curves in yellow and dashed curves in red in figures 2(a) and (b). The results show that both fit models can provide good predictions agreeing with the experimental measurements. Here, the prediction errors with y p being predicted data and y o being experiment data, have been evaluated as the fitting tolerance. Defining a good prediction as E p < 5%, it is found that the minimum time window is about 40 ms in the presence of the square wave perturbations. The minimum time window is 60 ms when the travelling wave is applied.
Here, the feasibility of extracting the dominant eigenmodes with a short time window is verified. The streaming analysis can be carried out to trace the evolution of eigenmodes' stabilities. The temporal evolution of the extracted stable eigenmodes is shown in figure 3, where the slipping time window size is Δt = 150 ms, and the window moving step, which is the actual updating time of eigenvalues, is δ = 10 ms. Note that, in the experiment, δ depends on the performance of the fitting process in the plasma control system (PCS). A DIII-D PCS test indicates that TDM analysis can achieve high efficiency and update the eigenvalues every δ = 4 ms, which is sufficient for real time detection. It is noted that, due to the existence of the resistive wall, the high frequency magnetic signal could be shielded, which implies the eigenmodes with very stable eigenvalues and high nature frequencies can be hardly resonant with the applied 3D fields. Thus, real time 3D MHD spectroscopy mainly focuses on the detection of low frequency global MHD modes. In figure 3, zero is the marginal stability boundary of each eigenmode. Re(γ) > 0 indicates the eigenmode is unstable. During the discharge, the two dominant eigenmodes have negative Re(γ) fluctuating with little change, since the equilibrium parameters in the flattop remain almost constant.
The MARS-F code [37,38] is employed to simulate the ideal plasma response for better understanding the extracted eigenmodes in the experiment. Figure 1 shows the equilibrium has little variation in the flattop phase. The equilibrium, reconstructed at t = 2200 ms of discharge 178 583 for the MARS-F simulation, is preferred since there is a quiet plasma and no The TDM analysis reports γ 1 = −9.96 + 0.16i Hz and γ 2 = −47.73 − 0.040i Hz. The similar eigenvalues, comparing the FDM and TDM analysis, further confirm the reliability of TDM analysis. In particular, the simulated eigenvalue of the least stable mode is close to the experimental measurement. The eigenvalue of the secondary mode is slightly different from experimental one, which could be due to the sensitivity of equilibrium reconstruction and the impact of other physics such as plasma rotation and resistivity which are not included in the simulation. By decomposing the contribution of each dominant mode in the plasma response based on FDM analysis [31], figure 4 shows the least stable mode is dominant at the HFS and the second mode is dominant at LFS. This result illustrates how necessary it is to measure the plasma response at multiple locations, which greatly help to better extract the multiple eigenmodes, as shown in [28][29][30]. Moreover, the MARS-F code can also solve eigenvalue problem of MHD equations with driving terms turned off in the simulation [39,40]. Note that it is challenging to solve the individual MHD mode when the mode is stable due to the existence of Alfvén continuum spectra. Alternatively, using the eigenvalue, extracted by real time 3D MHD spectroscopy, as the initial guess greatly helps to find the dominant eigenmodes. Figure 5 shows the different structures of two stable dominant n = 1 eigenmodes calculated by MARS-F. The least stable mode shows the core dominant structure in the radial displacement. The mode, with a kink like structure, shows a single dominant poloidal mode triggered in plasma regions with strong gradients of the plasma current. Although the m = 1, −1 poloidal Fourier harmonics of the radial displacement have big amplitudes near the magnetic axis, the large magnetic perturbation is observed at the plasma boundary in radial magnetic perturbations in figure 5(a), since B r ∼ (m − nq)ξ r is amplified while increasing the safety factor q. The secondary mode has more edge perturbations, and the toroidal coupling between poloidal harmonics indicates this eigenmode could be the seed of a peeling or a ballooning mode.

Experimental analysis in KSTAR tokamak
TDM analysis is also applied in KSTAR experiments for crossdevice validation. In this experiment, an n = 1 travelling wave is generated by the upper, middle and lower rows of internal 3D coils. The magnetic perturbations are measured by the toroidally distributed magnetic pick-up sensor arrays located at the middle plane of the LFS, where the configuration of magnetic probes are described in [41]. Despite the smaller number of sensors and a lower signal to noise ratio in the KSTAR data, a reliable TDM analysis is still possible by using a larger time window of 250 ms. Figures 6(c) and (d) present the applied  figure 6(a). The Re(γ) curves have higher variance (with prediction error E p ∼ 15%) than the case in figure 3, despite the steady KSTAR equilibrium parameters. A main reason is the higher noise level of magnetic measurements, since the sensor arrays are not dedicated to the 3D magnetic measurement in KSTAR. However, the general tendency of modes' stabilities is almost constant. One may eliminate the variance by averaging the eigenvalue within a short time window. To improve quantitative accuracy of the extracted eigenvalues, more sensors and/or better wave-filtering techniques are needed for the high quality of signals. Nevertheless, the strong anti-noise feature of the TDM guarantees the qualitative correctness of stability as long as the noise is incomparable to the plasma response.

Conclusion
In this paper, 3D MHD spectroscopy towards real time detection of multi-mode plasma stabilities is successfully developed and applied in both DIII-D and K-STAR stable plasmas. A new TDM has been applied to experimentally extracted the eigenvalues of multiple dominant eigenmodes by fitting the magnetic response measured on 3D sensors in a short time window (about 50 ms). In this work, the post-processing TDM analysis proves the feasibility of tracing the evolution of multimode instability. As with FDM analysis [31], the TDM is capable of quantifying the number of dominant modes and the corresponding eigenvalues. Meanwhile, the TDM shows great numerical efficiency towards the real time detection. Successful cross-machine application of TDM in both DIII-D and KSTAR tokamak indicates the technique can be transferred to multiple machines equipped with 3D coils and magnetic sensors. In comparison with the DIII-D case, the extracted multiple eigenvalues from KSTAR experimental data having higher variance shows how clean 3D magnetic measurements at multiple poloidal locations can greatly reduce the fitting time window of TDM and improve the reliability of eigenmode extraction. A short time window is critical for real time detection in the future to facilitate the stability monitoring to avoid plasma disruption. The updating time window relies on the system hardware, which in ITER would not have much difference with in DIII-D (less than 4 ms). The fitting time window is generally determined by the lowest frequency of applied perturbation, such as 10 Hz in DIII-D experiment. Provided the quality of response signal in ITER should be no worse than that in DIII-D, it is expected that the fitting time window should not be quite different from DIII-D one. The efficiency of the TDM will be investigated for ITER in the future.
Besides the great convergence and efficiency for real time detection, the TDM, as real time 3D MHD spectroscopy, can be combined with FDM to understand the contribution of the each dominant mode in the 3D plasma response by using the TDM extracted eigenmodes as the initial guess to reduce the nonlinear fitting uncertainty in the FDM analysis. Real time 3D MHD spectroscopy can also provide the reliable initial guess for the simulation of stable individual MHD modes and reveal the structure of modes. For instance, MARS-F simulation indicates that the secondary mode in figure 5 has more edge perturbation. Here, one may consider an application to find the optimal frequency and phasing to amplify this mode, which may help with the suppression of the edge localized mode. One key feature of TDM analysis is to identify the MHD stability quantitatively through the real parts of extracted eigenvalues. This work suggests using the square wave in the real time detection, since the shortest fitting time window of the square wave is shorter than that of the travelling wave. But, a more comprehensive comparison between the two waveform should be further investigated. TDM has been implemented in DIII-D PCS for testing the real time detection of plasma stability in live experiments. This form of real time 3D MHD spectroscopy detection of plasma stability could be significant for the prediction of MHD stabilities and therefore the avoidance of the severe disruptions in ITER and future fusion reactors.

Disclaimer
This report was prepared as an account of work sponsored by an agency of the United States Government. Neither the United States Government nor any agency thereof, nor any of their employees, makes any warranty, express or implied, or assumes any legal liability or responsibility for the accuracy, completeness, or usefulness of any information, apparatus, product, or process disclosed, or represents that its use would not infringe privately owned rights. Reference herein to any specific commercial product, process, or service by trade name, trademark, manufacturer, or otherwise does not necessarily constitute or imply its endorsement, recommendation, or favouring by the United States Government or any agency thereof. The views and opinions of authors expressed herein do not necessarily state or reflect those of the United States Government or any agency thereof.