Early internal detection of magnetic tearing and implications for tokamak magnetohydrodynamic stability

Tearing modes, both classical and neoclassical, pose a significant challenge to tokamak success, including that of ITER.^1,2 The emergence of tearing modes can limit achievable plasma density, temperature, or confinement time, each of which contributes to the $n T τ$ product that must be sufficiently high for ignition in a nuclear fusion reactor.³ Uncontrolled growth of tearing modes can even cause disruptions, which would be unacceptable for ITER or a fusion-based power plant. One way to achieve fusion conditions is operating a tokamak at a high toroidal plasma current, $I p$⁠, resulting in a low safety factor, $q ∝ 1 / I p$⁠, which can give rise to several tearing/kink mode instabilities.⁴ At lower plasma current and high-q, some of the most confinement limiting modes can be avoided, but this requires operating for sustained periods at large normalized plasma beta for an acceptable reactor scenario.⁵ A successful remediation of the tearing mode problem can be achieved either via complete avoidance by operating a tokamak in safe regimes or with sufficiently early detection/identification to allow mode suppression.^6–8 The detection of instability at earlier times, along with the ability to use models to predict its occurrence, will contribute toward this goal.

The onset of tearing modes often occurs when the normalized plasma beta, $β N$⁠, approaches the ideal wall limit. The normalized plasma beta, $β N = β × a B 0 / I p$⁠, where the plasma beta, $β = 2 μ 0 ⟨ p ⟩ / B 2$⁠, is the ratio of the volume-averaged plasma pressure, $⟨ p ⟩$⁠, and the vacuum toroidal field at the magnetic axis, B₀, in Tesla; a is the minor radius in meters, and $I p$ is the plasma current in mega-amperes. Modeling results have shown that when the ideal-wall $β N$ limit is approached from below, the linear resistive stability parameter, $Δ ′$⁠, increases to a large value, and hence the proximity to the ideal-wall kink limit is taken as a proxy for resistive tearing stability.^9–11 However, in experiments when tearing modes are detected with only the external magnetic sensors, tearing mode onset is not consistently predictable.¹² Some key questions motivating the present work are as follows: Does the correlation between the ideal-wall kink limit and the tearing mode occurrence hold true for all toroidal mode numbers? Could earlier identification of the time of the first occurrence of a mode help us improve the understanding of the magnetohydrodynamic (MHD) stability of a tokamak equilibrium and validate various models and theory that seek to predict the onset of tearing modes? Since the optimal strategies to manage classical and neoclassical tearing modes may be different, can we better identify a responsible instability?

In this Letter, we present detection of a tearing mode with internal measurements of the perturbed radial magnetic field using Faraday-effect polarimetry via the Radial Interferometer Polarimeter (RIP)¹³ diagnostic, and we explore the connection between the appearance of the mode with linear MHD in DIII-D high- $q min$ plasmas.¹⁴ The high- $q min$ scenario is being developed to achieve fully non-inductive operation with high $β N$ for simultaneously high fusion power and high bootstrap current fraction. Plasma equilibria consistent with these goals generally have broadly off-axis total current density and $q min$ > 2.^15,16 The operating point of a typical high- $q min$ discharge is above the ideal no-wall limit but below the ideal with-wall limit for an n = 1 kink mode.^17–20 While such plasmas are mostly resilient to disruptions caused by an n = 1 mode, the emergence of n > 1 modes can lead to reduced confinement and $β N$ and impact the scenario, for example, via the loss of non-inductive sustainment, undesired evolution of the current-profile, and reduced stability.^11,12,17,21 At the very least, scenario design needs to account for n > 1 modes and manage their impacts. With internal radial magnetic field fluctuation measurements, an n = 3 tearing mode is detected earlier than in the measurements with edge sensing coils—the use of such coils mounted on the vessel-wall is a common means of detecting MHD instability in a tokamak.^22,23 The linear stability analysis shows that, at the time of detection of the n = 3 mode, which is the first among the $n = 1 – 3$ modes to be detected, it is the closest to being unstable among them. Also, the equilibrium $β N$ is about 70% of the ideal-wall $β N$-limit for this n = 3 mode. By the time that this mode is both detectable and identifiable with the coils, $β N$ $∼ 100 %$ of the kink limit. Further, it is shown that the $Δ ′$ for this mode, calculated from the linear resistive stability analysis, remains near zero for a long duration after the mode onset. Therefore, while the appearance of the n = 3 mode is consistent with its $β N$ limit being smallest, the resistive stability analysis points toward linear stability. This is the first time to our knowledge that direct, internal measurement of magnetic fluctuations has been applied in this manner.

Figure 1 shows experimental parameters for a high- $q min$ discharge with an $I p$ flattop of about 0.8 MA. The plasma current is mostly non-inductively driven by a combination of neutral beam injection (NBI), electron cyclotron heating (ECH), and bootstrap current. The safety factor profile is flat in the core, $q 0 ∼ q min$⁠, and has a steady value at the edge, $q 95 ∼ 6.5$⁠. The line-averaged electron density is < 5 × 10¹⁹ m⁻³; the toroidal magnetic field is approximately 1.6 T. The ramp in the NBI power [Fig. 1(c)] from 5 to 8 MW beginning around 2500 ms results in $β N$ increasing from about 2 to 3.3. NBI power reaches a flat-top at 2600 ms, while the same for $β N$ occurs later at about 2800 ms. Here, the safety factor profile is obtained from kinetic EFIT²⁴ reconstructions.

The polarimetry technique provides a direct measurement of magnetic fluctuations in the core of a fusion-grade plasma;²⁵ other diagnostics can sometimes detect internal magnetic fluctuations indirectly by their impact on plasma parameters. RIP measures the line-integrated equilibrium and fluctuating components of density and magnetic field along the three horizontal chords shown in Fig. 1(e). They are located at the mid-plane and 13.5 cm above and below it. In this work, only the mid-plane chord is used. The Faraday rotation effect is observed in a magnetized plasma by the rotation of the polarization plane of a linearly polarized electromagnetic wave propagating along the magnetic field. The rotation angle is given by $ϕ FR ∝ ∫ n e B R d R$⁠, where R is the major radius, $n e$ is the electron density, and $B R$ is the major radial component of the equilibrium magnetic field. The estimated line-averaged fluctuating radial magnetic field, $δ B R$⁠, is given by²⁶  $δ B R ≡ ∫ n e d B R d R / ∫ ( n e d R )$⁠. The edge sensing coils with which RIP-measured tearing mode amplitudes are compared lie on the outboard midplane of the vacuum vessel wall.^22,23

The detection of the n = 3 mode with RIP and the sensing coils is shown in Fig. 2. The spectrogram corresponding to RIP measurements shows the onset of the n = 3 mode with a rotation frequency of 24 kHz beginning at around 2500 ms, then increasing to 60 kHz by 2950 ms. Extracted from the RIP spectrogram, the root mean squared (RMS) value of the line-averaged radial magnetic field fluctuation amplitude for only the n = 3 mode is shown in Fig. 2(b). The lighter gray lines include bursts from edge localized modes (ELMs) or other non-tearing modes, observable in the spectrogram, whereas the green points exclude most of this transient activity to better highlight the n = 3 evolution. The RMS $δ B p$ amplitude of the n = 3 mode measured by the toroidal array of sensing coils is shown in Fig. 2(c)—the mode is robustly growing from around 2800 ms. The n = 3 mode amplitude is shown here if the corresponding mode frequency lies within the 18–65 kHz range, bracketing the measured frequency range of the mode as measured by RIP in the time window shown. Within this range, the only coherent mode detected by the coils has the same frequency as that measured by RIP. The toroidal mode number is determined for the RIP-measured mode from the toroidal array of sensing coils, when the mode amplitude is large enough on the coils to allow for determination of the mode number. Then, the continuity of this mode on RIP measurements allows tracking the mode evolution back to the earliest time of detection with RIP. This RMS signal, computed automatically, during each DIII-D shot, is a commonly used metric for mode onset detection. Disagreement between the coil RMS and the RIP data on the time of first detection led to a further investigation of individual coils in the toroidal array. Not all coils were found to be equally sensitive—the most sensitive coil in the outboard mid-plane toroidal array of edge sensing coils detects this n = 3 mode at 2620 ms while many of them did not detect this mode until much later. When compared to the commonly used RMS metric, the mode appears on RIP 300 ms earlier, and about 120 ms before the most sensitive of the coils.

The coil-measured $δ B p$ for the n = 3 mode [Fig. 2(c)] grows barely above the measurement threshold around 2700 ms but does not persist as it does on RIP measurements. This small momentary increase in $δ B p$ is likely caused by an uptick in both the mode amplitude and the frequency. Because the edge coils measure a time-derivative of the poloidal field, $d / d t ( δ B p )$⁠, which is weighted by the mode frequency, a faster rotating mode in the laboratory frame results in a measurement larger than the noise floor. The RIP-data spectrogram shows that the n = 3 mode frequency is increasing up to around 2750 ms, then drops a bit, likely reducing the edge coil signal below the detection limit.

Starting at about 2800 ms, the amplitude of the poloidal field fluctuations [Fig. 2(c)] grows linearly in time—such linear growth is a signature of a neoclassical tearing mode.^27,28 The evolution measured by RIP [Fig. 2(b)] during this time is also roughly linear, but the precise evolution is more challenging to determine from RIP data given, for example, the residual impact of ELMs. It is also likely affected by the method with which the mode amplitude is determined from RIP data. A peak in the power spectrum is first determined for each time-window in a series of consecutive time-windows, and then for each time-window, the spectral power density is integrated over a 2 kHz range around this peak. Whereas for the edge coils, the mode amplitude is determined after fitting a toroidal mode structure over a finite temporal window. Because RIP measurements are scaled by the density profile, peaking and evolution of the density profile in the vicinity of the mode resonant surface may also contribute. The differing slopes between the two diagnostics during the linear growth phase is expected given the fundamental difference in measured mode amplitudes, as described further below.

The poloidal mode number of a mode can, in principle, be deduced from the poloidal array of sensing coils, but it is often not possible in the case of large m (m > 4) either due to strong plasma shaping or due to the limited number of coils in the poloidal array. Therefore, the poloidal mode number is determined here using the measured mode rotation frequency, the carbon-ion toroidal rotation profile,²⁹ and the EFIT-reconstructed q-profile. Within our analysis window indicated in Fig. 2, the n = 3 mode is rotating at frequencies consistent with the toroidal plasma rotation rate in the region $ρ ∼ 0.3 – 0.4$⁠. Here we consider only the toroidal plasma rotation, $f m , n = n ν tor / 2 π R m , n$⁠, where $f m , n$ is the mode frequency and ν_tor is the plasma rotation velocity at the major radial location, $R m , n$⁠, corresponding to the mode.³⁰ Further, the safety factor profile from the best available equilibrium reconstruction indicates that the closest n = 3 surface is the $q = 6 / 3 = 2$⁠.

The line-integrated polarimetry measurement can be more sensitive than that of the sensing coils for perturbations with high-m that originate in the plasma core. Estimated by $B r ̃ ∝ ( r / r s ) − ( m + 1 )$⁠, the eigenfunction of a high-m mode falls off more rapidly from the resonant surface than that of a low-m mode.³¹ The perturbed magnetic field at the location of a rational surface, $q = m / n$⁠, is much larger than the line-averaged value from RIP measurements. On the other hand, the perturbed field at the radius of the sensing coils is much smaller for high-m modes.³¹ Such modes can, therefore, go undetected by the coils, especially when they are of smaller intrinsic amplitude. The linear ideal MHD equation in cylindrical geometry has been used elsewhere to estimate the ratio of the radial field at the resonant surface and poloidal fluctuations measured at the edge by the sensing coils, $δ B R$/ $δ B p$⁠.^31,32 A similar calculation for comparison with line-integrated RIP measurements for this case of an m = 6 mode located at $ρ ∼ 0.35$⁠, assuming a flat density profile, shows that the ratio $δ B R$/ $δ B p$ is within an order of magnitude of the measurement. However, a precise ratio would depend on more precise knowledge of the mode resonant surface location and the perturbed eigenfunction—analysis in that direction will be presented in a future publication.

Measurements of the electron temperature profile show that the growing n = 3 mode causes flattening near the mode resonant surface. Figure 3 shows temperature profiles, $T e ( ρ )$⁠, at multiple times, obtained by fitting Thomson scattering measurements.³³ The vertical shaded region, $ρ ∼ 0.3 – 0.4$⁠, indicates the estimated range of the spatial location of the mode resonant surface as described earlier. The T_e profile flattening becomes more pronounced as the mode amplitude increases later in time. Additional measurements using the electron cyclotron emission (ECE)³⁴ show temperature fluctuations at the frequency of the mode starting from 2500 ms. The temperature fluctuations are evident in several adjacent spatial channels whose location is consistent with the flattening in Fig. 3.

To better understand the MHD stability of this equilibrium and the origin of the 6/3 mode, linear ideal and resistive stability analyses were carried out. The Direct Criterion of Newcomb (DCON),³⁵ an ideal MHD stability code, is used to examine the proximity of the experimental equilibrium to kink stability for an ideal wall for toroidal mode numbers $n = 1 – 3$⁠. The results are shown in Figs. 4(a)–4(c), where the vertical lines again represent the times of initial RIP and edge coil detections of the n = 3 mode, as in Fig. 2. The total energy dW determines whether a mode is stable (positive) or unstable (negative) to a perturbation of magnetic field lines. Each toroidal mode, n, in Fig. 4(a) encompasses multiple poloidal (m) harmonics. The $n = 1 – 3$ are all stable at about 2000 ms, but the n = 3 is closest to being unstable. Up to 2500 ms, the time at which the n = 3 mode is first detected by RIP, dW is roughly constant for all three modes, but then it begins to decrease. The stability limit for the n = 4 mode (not shown) is lower than that of n = 3 for t > 2000 ms. No n = 4 mode is detected with the edge coils. From the RIP spectrogram, a mode that may be n = 4 is discerned from around 2500 to 2600 ms by further processing the RIP data to remove ELMs.

Considering cases with an ideal wall, the $β N$-kink limits for $n = 1 – 3$ are compared in Fig. 4(b) with the experimental $β N$ (black solid line). The $β N$-limits are obtained from DCON by keeping the current profile and the total current the same as in the experiment, and scaling the plasma pressure until the equilibrium becomes unstable.¹² The experimental $β N$ is initially lower than the ideal kink limits for all three modes, but as it increases due to the increased NBI power [Fig. 1], it quickly approaches the n = 3 limit. Figure 4(c) shows the experimental $β N$ as a fraction of n = 3 ideal-wall kink limit—the n = 3 mode is detected at about 70% of the limit. In previous analyses¹² of shots like those described here, an n = 2 is the only robustly detected mode with edge coils. Its appearance was thought to be due to the proximity to the ideal-wall beta limit, and it affected the achievable $β N$⁠. As shown in Fig. 2, the n = 2 mode does appear in this shot, but only after the n = 3 mode has persisted for some time.

While the focus of this work is on MHD activity occurring during the mid-shot ramp in $β N$⁠, such activity is observed earlier in the shot as well. However, the shape and elongation of the plasma is changing until around 2000 ms when the final configuration is attained, and no coherent MHD activity is observed from 2000 to 2500 ms. Even though the early observed n = 3 (⁠ $m / n = 6 / 3$⁠) resides on the same surface as would the lower harmonics, $n = 2 , 1$ (4/2, 2/1), neither of these lower harmonics, for which the $β N$ limits are higher, are observed until several hundred milliseconds after the detection of the n = 3 mode. In spite of the above arguments suggesting a mechanism similar to previous observations^9,11,12 where approaching $β N$ limits causes tearing mode onset, it is unlikely that this mode originates as a classical tearing mode based on the resistive stability analysis discussed next.

The resistive stability of the time-evolving equilibrium is examined with the resistive-DCON (RDCON) code.³⁶ The resistive stability parameter, $Δ ′$⁠, for the $m / n =$ 6/3 and 2/1 modes is shown in Figs. 4(d) and 4(e), respectively. Initially, $Δ 6 / 3 ′ ∼ 0$ and increases to a large value later in time, when the equilibrium is ideally unstable. While $Δ ′ > 0$ is often considered to be marginal for linear stability,³⁷ it may exceed a critical value before the equilibrium is actually unstable.^38,39 The $Δ 2 / 1 ′$ is positive and increases rapidly later when the equilibrium is classically unstable at 2800 ms, but a 2/1 mode is not observed in the experiment at this time. The n = 3 mode is clearly detected at 2500 ms when the $Δ 6 / 3 ′$ is at a small saturated value, well before the time when it is rapidly increasing. If the mode were only detected at 2800 ms then it would point to a mechanism where a classical mode grows in proximity to the ideal kink limit accompanied by a rapidly increasing $Δ ′$⁠. Further, the linear growth of the mode after 2800 ms [Fig. 2] shows that the 6/3 mode from this point onward is likely a neoclassical tearing mode triggered by a mixed seeding mechanism.³⁸ Also note that, in general, the results of linear stability analysis do not hold once a mode has grown to sufficient amplitude for an island to form. The flattening of the temperature profile, likely indicating the presence of an island, is most apparent starting around 2800 ms, as the mode amplitude begins to increase.

The DCON computations shown here utilize equilibrium reconstructions obtained by constraining the central safety factor, $q 0 ( t )$⁠, to vary smoothly with time. Equilibria without such a constraint resulted in unphysically large variation in $q 0 ( t )$⁠, as seen in Fig. 1(b), likely because of inadequate constraint on the core current density from available data, further leading to a larger variation of $Δ ′$ for consecutive time slices. Variations in $Δ ′$⁠, not presented here, have been shown to be more sensitive to the equilibrium reconstructions than the standard deviation in its computation.^9,40 In Fig. 4(d), the average temporal trend of $Δ ′$ remaining close to zero for a long time interval is, therefore, taken to indicate the plasma being linearly stable to the n = 3 perturbation.

In summary, we have shown that internal RIP measurements allow direct detection of an n = 3 tearing mode well before the sensing coils. The mode detection precedes temperature profile flattening and likely affects the local pressure and current profile evolution. Though the proximity to $β N$-limits is certainly observed, the mode is not classically driven. Additionally, the subsequent linear growth of the n = 3 mode suggests neoclassical tearing, but the cause for the initial onset of the mode is not clear. Among the possible triggers or drivers of neoclassical tearing modes relevant to the discussed experiment include ELMs,²⁸ and the fast-ions from neutral beam injection.^41,42 Finally, recall that the $m / n = 6 / 3$ mode occurs at the q = 2 surface without its lower harmonics. The 4/2 and 2/1 modes do appear later in time, though establishing a connection between the n = 3 with the emergence of the n = 2, 1 modes is a topic of future work.

The use of internal measurements to validate computational models such as DCON is an important step toward employing them as predictive tools. Detecting the instabilities early in time and/or at lower amplitude may also be important, especially in cases where high-m instabilities arising in the core potentially lead to the emergence of lower-m instabilities. Even for $q min$ > 2 cases, lower-m instabilities have been observed at the q = 3 surface following the appearance of higher-m instabilities.¹⁸ The combination of understandings gained from internal measurements and resulting stability analysis may point to potential new means by which tearing modes can either be prevented or controlled.