Dynamics of an $n=1$ explosive instability and its role in high-$\beta$ disruptions

TL;DR

This paper investigates the nonlinear behavior of an $n=1$ explosive instability in high-$eta$ plasmas, revealing bifurcations between benign saturated modes and explosive disruptions influenced by background conditions.

Contribution

It identifies the bifurcation mechanisms and nonlinear paths of an $n=1$ mode, highlighting the conditions leading to either stable saturation or explosive disruption in high-$eta$ scenarios.

Findings

Existence of a long-lived saturated mode with minimal confinement loss.

Transition from exponential to explosive growth in ballooning fingers.

Metastability of the saturated mode and its potential for explosive destabilization.

Abstract

Some low-$n$ kink-ballooning modes not far from marginal stability are shown to exhibit a bifurcation between two very distinct nonlinear paths that depends sensitively on the background transport levels and linear perturbation amplitudes. The particular instability studied in this work is an $n=1$ mode dominated by an $m/n=2/1$ component. It is driven by a large pressure gradient in weak magnetic shear and can appear in various high-$\beta,$ hybrid/advanced scenarios. Here it is investigated in reversed shear equilibria where the region around the safety-factor minimum provides favorable conditions. For a certain range of parameters, a relatively benign path results in a saturated "long-lived mode" (LLM) that causes little confinement degradation. At the other extreme, the quadrupole geometry of the $2/1$ perturbed pressure field evolves into a ballooning finger that subsequently transitions from exponential to explosive growth. The finger eventually leads to a fast disruption with precursors too short for any mitigation effort. Interestingly, the saturated LLM state is found to be metastable, it also can be driven explosively unstable by finite-amplitude perturbations. Similarities to some high-$\beta$ disruptions in reversed-shear discharges are discussed.

Figures

Equations

$$\delta W_p$$

$$Q_{\perp }^2\simeq {\xi }_{\perp }^2(m-nq{)}^2\simeq {\xi }_{\perp }^2(n{q}_s^{\prime }{)}^2(r-r_s{)}^2.$$

$$W_K(t)=W_K(t_i)[(t_f-t_i)/(t_f-t){]}^{\nu },$$

$$\ddot{\xi }=-\frac{\partial \varphi }{\partial \xi }-\mu \dot{\xi },\xi (0)=0,\dot{\xi }(0)=v_0\ge 0,$$

Full text

1 Introduction

In tokamaks, crossing certain operational boundaries in plasma density ($n_{e}$)[1], current ($I_{p}$)[2, 3], or pressure ($\beta$)[4] can lead to a disruption, a sudden and uncontrolled loss of thermal and magnetic energy in the plasma. Among these, high-$\beta$ disruptions are particularly challenging, not only because of the high thermal energy content of the plasma, but also because of their extremely fast time-scales in some cases.

Most high-$\beta$ disruptions are mediated by a (neoclassical) tearing mode (NTM), typically with the mode numbers $m=2,n=1,$ that for various reasons lock to the wall, grow in size and eventually cause a loss of confinement[5, 6]. Even when they do not lead to disruptions, NTM’s tend to degrade confinement significantly so that their avoidance or stabilization in the ITER ELMy H-mode baseline scenario has been a high-priority research item (see for example [7, 8]). When the plasma $\beta$ is pushed higher beyond the no-wall limit in “hybrid/advanced tokamak” regimes, more dangerous $n=1$ kink modes can become unstable. In the presence of a close-fitting wall, these are generally transformed into slow-growing resistive wall modes (RWM’s)[9]. Again, if they are allowed to lock to the wall, RWM’s can lead to disruptions. Fortunately, plasma rotation[10, 11, 12], kinetic effects[13, 14, 15], coupled with feedback-control methods[5, 16, 17], can stabilize RWM’s well above the no-wall $\beta$ limit.

Since both NTM’s and RWM’s grow on a slow, resistive time scale, disruptions caused by these modes are easily identified by their long precursors on various diagnostics. In fact, because of their relatively slow time scale, these are precisely the type of disruptions that are targeted by various disruption mitigation schemes, which require at least a few 10’s of milliseconds of warning time[18, 19]. As stated earlier, however, tokamaks disrupt for a wide variety of reasons (see for example [6]), and not all disruptions follow this slow path where their arrival is well-advertised in advance; some in fact occur with little warning. Unfortunately, their very fast time scale apparently makes detailed studies difficult, and it is likely that their rare appearance in the literature does not accurately reflect their actual frequency in the experiments.

There do exist some documented high-$\beta$ disruptions with precursors of the order of a millisecond or less. For example: $\beta$-limit disruptions in TFTR due to toroidally localized ballooning modes in the presence of $n=1$ magnetohydrodynamic (MHD) activity [20], localized resistive interchange modes that couple to a global $n=1$ mode and lead to a disruption in negative central shear (NCS) discharges in DIII-D[21, 22], and disruptions following an internal transport barrier (ITB) collapse in JET[23]. In these discharges, some of the important details were clearly different: in TFTR, at least initially, the $q=1$ surface was involved, whereas DIII-D and JET presumably both had $q_{min}\simeq 2.$ But generally, a large pressure gradient in regions of weak magnetic shear is believed to have played an essential role. Thus, one of the things we will do in this work will be a short review of the resistive and ideal stability of such configurations.

However, linear stability analysis alone cannot explain the fast time scale of these disruptions. The mode that is involved has to be growing near Alfvénic rates to account for the time scale, but it is not clear how a discharge evolving on the slow transport time scale can generate an unstable mode with a near-Alfvénic growth rate without producing a long series of precursor oscillations during its sub-Alfvénic period.

The point raised above is in fact part of a more generic problem: If an event (e.g., a sawtooth crash, edge-localized mode (ELM) crash, disruption, etc.) that makes macroscopic changes in the state of a discharge in a time scale $\tau_{e}$ is attributed to some global instability, the instability growth rate $\gamma_{e}$ at the time the event is observed has to be commensurate with that time scale, i.e., we need to have $\gamma_{e}\sim{\cal O}(1/\tau_{e}).$ We can safely assume $1/\gamma_{e}\sim\tau_{e}\ll\tau_{t},$ where $\tau_{t}$ is a characteristic transport time scale; otherwise the “event” cannot be distinguished from ordinary transport. If we assume that the changes in the mode growth rate occur entirely because of modifications to the background equilibrium by transport, then there will necessarily be a long period during which $1/\tau_{t}<\gamma(t)<\gamma_{e}$, i.e., a time when the mode is growing faster than the transport rate but does not yet have the eventual “crash” rate. But then we are faced with two related questions: (i) During this period, can the mode grow without being detected? The short answer is probably “no,” since it is hard to imagine how such a mode could avoid generating a long series of precursors during the period mentioned. (ii) Since it is now growing faster than the transport time scale, would it not “self-stabilize” and saturate without causing the “event” by modifying the equilibrium faster than transport processes? Here a general answer is again difficult, but a mode that depends very sensitively on local conditions for stability can probably “self-stabilize” and saturate more easily than a global mode like an $m=1$ resistive kink or one that is responsible for a major disruption.

Thus, disruptions or other events that occur without long precursors seem to require a different evolution scenario than the one proposed above. Instead of the mode growth rate $\gamma(t)$ slowly evolving with the background equilibrium, we have to consider mechanisms that can make changes in $\gamma(t)$ at a rate much faster than expected from transport alone.

A mechanism proposed by Hastie[24] for fast sawtooth crashes and by Callen[25] for the DIII-D disruption mentioned earlier assumes that, in response to a linearly increasing plasma pressure driven by auxiliary heating, the growth rate of a pressure-driven mode grows as $\gamma(t)\propto\beta^{1/2}\propto t^{1/2},$ where the plasma $\beta$ is defined as $\beta=2\mu_{0}\left\langle p\right\rangle/B^{2}$ and $\left\langle.\right\rangle$ denotes a volume average. Then for large times, the plasma displacement, assumed to satisfy ${\ddot{\xi}}=\gamma^{2}\xi$, can be shown to evolve as $\xi(t)\propto\xi_{0}\exp{(t/\tau)^{3/2}}$, where $\tau\simeq(\tau_{MHD}^{2}\tau_{h})^{1/3}$. Here $\tau_{MHD}$ is the fast MHD time scale, $\tau_{MHD}=L/v_{A},$ where $v_{A}=B/\sqrt{\mu_{0}\rho_{m}}$ is the Alfvén velocity, $L$ is a global length scale, and $\tau_{h}$ is the slow “heating time,” comparable to the “transport time” mentioned earlier, $\tau_{h}\sim\tau_{t}$. Thus, this mechanism seems to lead to a faster-than-exponential growth and possibly explain the near-absence of precursors before some sawtooth crashes or disruptions. However, Cowley[26] has shown that, because of the large separation in the MHD and transport time scales, this path to super-exponential growth requires an unrealistically small initial perturbation and can be ruled out.

There are, however, nonlinear processes in plasmas that can generate explosive (faster-than-exponential) growth while the underlying mode is still not far from marginal stability. In a numerical study of the semi-collisional/collisionless $m=1$ mode using a reduced two-fluid model, nonlinearities involving the parallel pressure gradient were shown to give a near-exponential increase in the growth rate of the mode[27], providing a possible explanation for precursor-less, fast sawtooth crashes. Similarly, Cowley and colleagues[26, 28, 29, 30] have shown that the nonlinear evolution of pressure-driven modes can generate finite-time singularities, again demonstrating how a long period of precursors can be avoided during a fast disruptive event.

In this work we extend our study of a specific example[31], a pressure-driven $n=1$ kink-ballooning mode that can continue to grow exponentially well into its nonlinear regime and become explosive with an apparent finite-time singularity at the end. We show that it can actually exhibit two very different types of nonlinear behavior depending on small differences in the assumed transport levels and linear perturbation amplitudes. In addition to the explosive behavior, it can also display a more benign evolution and saturate in a “long-lived mode” (LLM) with only minor confinement degradation[32]. The LLM itself is shown to be a meta-stable state; it can be pushed into the explosive regime with small changes in the transport coefficients or with a finite-size perturbation.

The experimental context for this computational study is KSTAR discharges with $q_{95}\simeq 7,~{}q_{0}\simeq 2$, and a low inductive current fraction, similar to some hybrid/advanced scenarios[33]. With on axis electron cyclotron resonance heating (ECRH), the pressure profile peaks and drives an ideal $m/n=2/1$ mode that saturates at a small amplitude. The resulting long-lived mode (LLM) survives many tens of seconds (as long as ECRH is maintained), with only a small effect on confinement[32]. Although there are no documented examples for the explosive version of this mode in KSTAR, in the absence of any detailed study of KSTAR disruptions, their existence cannot be ruled out. The computational tool used here is the CTD code, which solves the nonlinear MHD equations in toroidal geometry (see [34] and the references therein).

Before moving on to a discussion of the nonlinear results, in the next section we briefly review the salient features of the linear stability of pressure-driven modes.

2 Linear stability

A general understanding of the stability of pressure-driven modes can be obtained from a cursory examination of the ideal MHD energy integral, written here in its “intuitive form” (plasma contribution only)[9, 35]:

[TABLE]

The largest stabilizing contribution tends to be the $|\bm{Q}_{\perp}|^{2}$ term in the first integral representing the line-bending energy, where $\bm{Q}=\nabla\times{\bm{\xi}}\times{\bm{B}}_{0}$ is the perturbed field. The destabilizing pressure-gradient and parallel current terms are grouped together in the second integral. The pressure gradient makes a destabilizing contribution to $\delta W_{p}$ only in those regions where the field line curvature ${\bm{\kappa}}\equiv{\bm{b}}\cdot\nabla{\bm{b}}=(\mu_{0}/B^{2})\nabla_{\perp}(p+B^{2}/2\mu_{0})$ is “unfavorable,” i.e., where ${\bm{\kappa}}\cdot\nabla p>0.$ By having the displacement vanish where the curvature is favorable, ${\bm{\kappa}}\cdot\nabla p<0$, the net destabilizing contribution from the pressure forces can be maximized. This is the path the ballooning modes take, but they pay a price in excess line-bending energy since the perturbation is not constant along the field lines.

Another path for pressure-driven instabilities opens up if the magnetic shear is weak in a region of finite width. Simplifying and expanding $\bm{Q}$ around a rational surface, we have $\bm{Q}_{\perp}\simeq{\bm{\xi}}_{\perp}(i{\bm{k}}\cdot{\bm{B}}_{0})$, and

[TABLE]

Thus, if the global shear $s\equiv rq^{\prime}/q$ is weak enough in regions with strong pressure gradients, interchange-like modes become possible even in “Mercier-stable” equilibria with $q^{2}>1,$ first recognized by Zakharov[36]. In fact, a rational surface is not necessary for instability. With $q\simeq(m+\epsilon)/n,~{}0<\epsilon\ll 1,$ and $s\ll 1,$ the line-bending energy can be minimized again since $Q_{\perp}^{2}\simeq\xi_{\perp}^{2}\epsilon^{2},$ which can be overcome by a strong-enough pressure gradient.

This simple piece of physics, strong pressure drive coupled with weak shear, is behind the quasi-interchange mode[37, 38, 39, 40] (for $q_{min}\simeq 1$) and the “infernal” modes[41, 42, 43] (for $q_{min}>1$), both pressure-driven modes in low-shear equilibria. The former was studied in the context of fast sawtooth crashes caused by an internal $m/n=1/1$ mode. The latter are particularly dangerous global modes that can be unstable much below the $n\rightarrow\infty$ ballooning limit and lead to major disruptions. They are typically thought of as “low-$n$” modes, but of course the same physics can also make the $n=1$ mode unstable, which will be the focus of this work.

For computational economy, our earlier work[31] focused on the nonlinear evolution of pressure-driven $n=1$ modes in circular geometry, in both monotonic and weakly-reversed $q$ profiles. Here we will extend it to non-circular geometry and provide a brief review of the linear properties of the relevant modes. Partly because of KSTAR’s recent interest in advanced scenarios with internal transport barriers (ITB’s), we will mainly consider reversed-shear equilibria with $q_{min}>2.$

Typical equilibrium profiles used in the linear and some of the nonlinear calculations are shown in Fig. 1. The shaped geometry has $\kappa=1.5$ (elongation) and $\delta=0.6$ (triangularity) within a perfectly conducting boundary; these geometric parameters are held fixed, except when we revisit circular-geometry. The CTD code uses a conformal transform from the poloidal plane to a unit circle in $(\rho,\omega)$ coordinates to deal with weakly-shaped equilibria[44]. The coordinate axis is shifted to approximately align the $\rho=const.$ surfaces with flux surfaces, but $\rho$ is not a flux coordinate. For this reason, the plots as in Fig. 1 (b,c) show both the $\omega=0$ (outboard) and $\omega=\pi$ (inboard) sections of the mid-plane. Note that a simple pressure profile without an internal or edge transport barrier is used to simplify the discussion.

As expected, resistivity enlarges the instability domain for the $n=1$ mode ($n>1$ stability is not considered in this work) so that an unstable mode is observed well below the ideal MHD stability limits. However, we find that the nature of the unstable resistive mode can be confusing. The “infernal” mode theory predicts a mode with a tearing scaling, $\gamma\tau_{A}\propto S^{-3/5},$ at low $\beta$. Close to the ideal stability boundary, the resistivity scaling is weaker, $\gamma\tau_{A}\propto S^{-3/13},$ becoming independent of $S$ beyond the ideal limit[42]. Here we define $\tau_{A}$ as the shear-Alfvén time, $\tau_{A}=R_{0}/v_{A},~{}v_{A}=B/\sqrt{\mu_{0}\rho_{m}}$. Then the magnetic Reynolds (Lundquist) number is given by $S=\tau_{R}/\tau_{A},$ where $\tau_{R}=\mu_{0}a^{2}/\eta$ is the resistive diffusion time, and $a,R_{0}$ are the minor and major radii of the torus, respectively. Throughout this work, $S$ is defined in terms of the value of resistivity at the coordinate axis, but the resistivity itself is in general a function of the poloidal coordinates such that $\eta(\rho,\omega)J_{\zeta 0}(\rho,\omega)=E_{0}=const.$ The electric field $E_{0}$ is associated with the “loop voltage.” Numerically it is used to prevent the Ohmic diffusion of the equilibrium current during long nonlinear calculations.

The well-known interchange theory also predicts a resistive mode with the usual interchange scaling, $\gamma\tau_{A}\propto S^{-1/3}$[45], but only for reversed-shear equilibria. If we briefly recall the relevant theory, the Mercier (ideal interchange) modes are unstable for $D_{I}\equiv D_{M}-1/4>0,$ where $D_{M}=-(2\mu_{0}rp^{\prime}/B_{\zeta}^{2}s^{2})(1-q^{2})$ in circular geometry[46]. Although rare, Mercier modes have actually been observed experimentally[47]. Resistive instabilities require $D_{R}\equiv D_{I}+(H-1/2)^{2}>0,$ where $H\propto-p^{\prime}/q^{\prime}.$ We see that $D_{R}>0$ is possible, even when Mercier stable ($D_{I}<0$), with weakly-reversed shear ($H<0$) at high enough $\beta.$ At lower $\beta$ this mode also reverts to the tearing scaling with $S^{-3/5}$.

Stability of the equilibrium in Fig. 1 is summarized in Fig. 2, where we plot the growth rate of the $n=1$ mode as a function of the magnetic Reynolds number $S$ for various values of the normalized $\beta,~{}\beta_{N}=\beta({\bar{a}}B)/I_{p}.$ Here ${\bar{a}}=a[(1+\kappa^{2})/2]^{1/2}$ is an equivalent minor radius defined for an equilibrium with elongation $\kappa$, and $I_{p}$ is the plasma current. The $S$-scans are performed at a constant magnetic Prandtl number, $P_{M}=\mu/\eta=10,$ where $\mu$ is the normalized viscosity coefficient. Although normalized viscosity tends to be higher than resistivity in fusion plasmas, this value of $P_{M}$ is chosen entirely for numerical reasons.

At $\beta_{N}=1.82$ (Fig. 2 (a)), there is a weakly unstable resistive mode; both the infernal mode and the resistive interchange theory seem to predict here a tearing-like scaling with $\gamma\tau_{A}\propto S^{-3/5}$ (the blue dashed line), but we find that a stronger dependence with $\gamma\tau_{A}\propto S^{-3/4}$ (the red line) is a better fit to the numerical data. Neither one of these theories takes into account viscous effects. The classical viscous-tearing mode theory that assumes $P_{M}<1$ predicts a mode with the $S^{-2/3}$ scaling[48], which is somewhat weaker than the $S^{-3/4}$ scaling we observe. It is possible that for $P_{M}\gg 1$, the $S^{-2/3}$ scaling changes to $S^{-3/4}$, but that possibility has not been investigated.

At $\beta_{N}=2.62$ (Fig. 2 (b)), the mode is still resistive and has a clear resistive interchange scaling, $\gamma\tau_{A}\propto S^{-1/3}$. Here it is possible that there is a resistive infernal mode (with the $S^{-3/13}$ scaling[42]) that is in competition with the interchange, but it is not observed numerically. The weak reversed shear (not considered in the infernal mode theory) may be making the resistive interchange the dominant mode in this particular parameter regime.

At an even higher $\beta$ ($\beta_{N}=3.35$, panel (c)), the mode is ideally unstable (with wall). During this study, $\beta_{N}$ was increased in large steps so that the locations of the transition points between various regimes are not known, a task left for a future work.

Since there is no $q=2$ rational surface in the plasma for the series of equilibria considered here ($q_{min}>2$), and because of the wide region of weak magnetic shear around $q_{min}$ (see Fig. 1), the eigenfunctions do not exhibit a distinctive “singular” behavior there. In fact, they have the features of a global kink mode, as seen in Fig. 3. For $\beta_{N}=1.82$ (panels (a,b)), there is a strong coupling to an $m/n=3/1$ mode near the edge, a $\Delta^{\prime}$-stable tearing mode. At $\beta_{N}=2.62$ (panels (c,d)), this coupling seems to become weaker, probably explaining the change in scaling from “tearing”-like to $S^{-1/3}$ mentioned earlier. Note that although only two modes are shown, the linear calculations for $n=1$ included all poloidal mode numbers in the range $m\in[-7,57].$

Above we discussed in some detail the linear stability for $q_{0}=2.15,~{}q_{min}=2.02,$ mainly to place our nonlinear calculations below in some context. Summarizing our other linear results, for a more deeply-reversed equilibrium with $q_{0}=2.57,~{}q_{min}=2.02$ we find no ideal instability for $\beta_{N}\leq 3.08,$ the limit of our numerical explorations for this $q$-profile. On the other hand, for an equilibrium with $q_{0}=2.03,~{}q_{min}=2.02$ (much weaker central shear), we find that $\beta_{N}\geq 2.72$ is ideally unstable, although the actual stability boundary has not been explored and is probably lower (but higher than $\beta_{N}=1.90$, where we find a resistive mode).

3 The explosive instability and disruptions

Nonlinearly the pressure-driven $n=1$ mode can turn into an explosive instability, as it was first demonstrated in circular geometry[31]. Here we discuss the nonlinear evolution of the mode in shaped geometries using the linear results of the previous section as starting points.

An ideally unstable $n=1$ mode (e.g., one from Fig. 2 (c)) with its large growth rate will naturally lead to a fast disruption. As discussed at some length in the Introduction, however, MHD modes are not “born” in this robustly unstable state. They tend to come into existence as weak resistive instabilities as the equilibrium slowly (on transport time scale) passes through some marginal stability point due to the evolving discharge conditions. Hence our goal in this section is to demonstrate how a weak resistive instability can evolve into a robust mode that will result in a fast disruption with only a brief period of precursors.

Thus we start with a weakly unstable equilibrium similar to that of Fig. 1; some of the relevant parameters are: $q_{0}=2.145,q_{min}=2.023,q_{l}=9.424,\beta_{N}=1.67$, which results in an even weaker instability than that of Fig. 2(a). The nonlinear calculations are performed at $S=10^{6},~{}P_{M}=\mu/\eta=10,$ using 21 toroidal Fourier modes. The poloidal Fourier expansions have $m\in[0,64]$ for $n=0$ and $m\in[-5,64]$ for $n\in[1,20].$ The finite difference scheme in the radial direction uses $192$ grid points. Some of the algorithmic details of the CTD code used here can be found in [34] and the references therein.

With weak shear and $q_{min}>2$, an important feature of the linear eigenfunction is the dominance of the $m=2$ poloidal component. This is clearly seen in Fig. 4(a), which shows the quadrupole geometry of the pressure perturbation (some coupling to an $m=3$ on the outside is also visible). This perturbation leads to an elliptical deformation of the flux surfaces in the core plasma that eventually forms a ballooning finger, as seen in Figs. 4(b-c). The finger pushes through the flux surfaces on near-Alvénic time scales and brings the core plasma in contact with the boundary (Fig. 4(d)). Using typical parameters for modern tokamaks ($B=3~{}T,~{}n_{e}=10^{19}~{}cm^{-3},R_{0}=3~{}m$), the state with no visible deformation in Fig. 4(b) and the final state in Fig. 4(d) are separated by less than $1~{}ms.$ Thus an actual disruption following this path would have a very short warning time.

The rapidity of the final disruptive phase is clearly due to a nonlinear increase in the growth rate rather than a slow, transport time scale change. This is seen in Fig. 5 where the kinetic energy in the $n\geq 1$ modes (excluding the $n=0$ equilibrium flows) is plotted. Two points are immediately obvious: (i) The mode continues to exponentiate well into the nonlinear regime with a growth rate $\gamma\tau_{A}=1.70\times 10^{-3}$ (the dashed red line in panel (a)). (ii) Instead of saturation, the late stages are characterized by a superexponential or explosive growth. In fact this phase has the appearance of a a finite-time singularity (panel (b)) where the growth has the form

[TABLE]

where $t_{i}=2773.9,~{}t_{f}=4277.6,$ and the exponent $\nu=2.05.$ In our earlier calculations in circular geometry the explosive phase was even faster, with $\nu=3.37$[31]. The slowdown seen here has both physical and numerical sources: The underlying $n=1$ mode is inherently more stable in shaped geometry. And, because the poloidal spectrum for each toroidal mode is much wider in shaped geometry, fewer ($n_{max}=20$) toroidal modes were used here. In circular geometry we were able to use $n_{max}=30,$ which allowed us to continue the calculations further into the explosive phase.

The ballooning finger of Fig. 4 has an extended structure along the field lines with $q_{s}=m/n=2/1$ helicity. As seen in Fig. 6, this symmetry is preserved as the finger moves outward. Because of magnetic shear, however, it becomes more localized in both parallel (along the field) and perpendicular directions (compare panels (a) and (b)) as it moves through regions of increasing safety factor in order to reduce the line-bending energy. Thus, under experimental conditions when it eventually comes into contact with the wall or limiter, its thermal content will be deposited in a small area and possibly cause serious damage.

There have been observations of ballooning fingers during high-$\beta$ disruptions in TFTR[20, 49]. JET also has reported similar results. During disruptions following an ITB collapse, localized disturbances in the ECE data that propagate from the ITB to the edge at velocities approaching 3km/s are seen[23]. These experimental observations can be associated with the radial propagation of a ballooning finger as shown here. In fact, a comparison of synthetic ECE diagnostics from our calculations with the JET data from Ref. [23] (their Fig. 4) shows good agreement, as seen in Fig. 7.

Although there are several observations in KSTAR that may be associated with fast, high-$\beta$ disruptions without significant precursors, the necessary analysis to link them formally to an explosive instability has not been carried out.

4 Bifurcated states

Generally, away from exact marginal points, we expect small changes in a relevant parameter to result in similarly small changes in the evolution of an unstable mode. Thus, small differences in the resistive dissipation level rarely have a significant impact on the saturation width of a tearing mode. However, there are counter examples where, for instance, an increase in the Prandtl number (ratio of viscous to resistive dissipation) beyond a threshold leads to a qualitatively different nonlinear regime[50].

Here we demonstrate an extreme case where a small change in a transport coefficient leads to a bifurcation between a benign, saturated state (the long-lived mode, LLM) and an explosive instability for the $n=1$ kink-ballooning mode. For computational economy, we expand upon our earlier results[31, 32] while staying in circular geometry. The bifurcation is summarized in Fig. 8 where we follow the nonlinear evolution of the mode starting with the same initial conditions and linear perturbation, but using slightly different transport coefficients. With $S=10^{6}$, thermal conductivity $\kappa_{\perp}=4\times 10^{-6}$ and viscosity $\mu=1\times 10^{-5}$, the mode goes through an exponential growth phase but saturates at a small amplitude (curve (1) in Fig. 8(a)). This regime is identified with the long-lived mode (LLM) observed in KSTAR, where an $m/n=2/1$ perturbation is seen in the electron cyclotron imaging (ECEI) data for many tens of seconds during the current flat-top period. The experimental conditions under which the LLM was observed is described in more detail in [32, 51].

In a slightly less dissipative system with $\kappa_{\perp}=3.5\times 10^{-6}$ but $\mu=1\times 10^{-5}$ still, the mode initially follows a similar path and goes through an exponential-growth phase (curve (2) in (a)). However, instead of saturating, it gradually enters a super-exponential regime where the growth rate itself increases rapidly, eventually becoming explosive. A similar effect is seen with reduced viscosity: with $\mu=6\times 10^{-6},~{}\kappa_{\perp}=4\times 10^{-6}$ the mode again transitions into the explosive regime (curve (3)).

This bifurcation can be understood qualitatively if we assume the explosive phase has a finite threshold in the perturbation amplitude. Dissipation affects both the linear growth rate and the nonlinear saturation amplitude of the unstable mode. The higher dissipation level clearly causes the mode to saturate below the apparent threshold. This point is confirmed in the next section where we show that the curve (1) of Fig. 8(a) represents a continuous set of metastable states.

5 Metastability

The $n=1$ mode for the set of equilibria we consider here is linearly unstable, which implies that an infinitesimal perturbation will grow in time exponentially, at least until the mode attains a finite amplitude. The results of the previous section imply that whether it turns into a LLM or becomes explosive is determined by a critical perturbation amplitude that itself is a function of the dissipation coefficients, $\xi_{cr}=\xi_{cr}(\eta,\mu,\kappa_{\perp})$. If it nonlinearly saturates below the threshold, $\xi<\xi_{cr}$, the result is a benign long-lived mode. Above $\xi_{cr}$ it becomes explosive. Thus, the pressure-driven $n=1$ mode is said to exhibit metastability. In this section this theoretically predicted behavior[30] is demonstrated numerically.

To show metastability, a long-lived mode, identified by the dashed curve (1a) in Fig. 9, is perturbed at various points along its trajectory with $n=1$ perturbations of varying amplitude. The relevant transport coefficients for this baseline, metastable state were: $S=10^{6},$ thermal conductivity $\kappa_{\perp}=4.75\times 10^{-6},$ viscosity $\mu=1.0\times 10^{-5},$ and the initial equilibrium was perturbed with the $n=1$ linear eigenfunction using a perturbation amplitude of $\epsilon=1.0\times 10^{-6}$ (The exact meaning of this parameter is not as important as its relative amplitude.). The same equilibrium, when perturbed with $\epsilon=2\times 10^{-6}$, evolves into an explosive mode (curve 1b). Similar numerical experiments are performed at later points in the nonlinear evolution of the baseline LLM: At $t=1012,~{}\epsilon=4.0\times 10^{-6}$ (curve 2a) is stable, however $\epsilon=5.0\times 10^{-6}$ becomes explosively unstable (curve 2b). At a later time, $t=1256,$ the same perturbation ($\epsilon=5\times 10^{-6}$) leads to large-amplitude, damped oscillations but the mode remains stable (curve 3a). A slightly larger perturbation ($\epsilon=6\times 10^{-6}$) again becomes explosively unstable (curve 3b).

Clearly, the quasi-equilibrium states representing the early nonlinear phase of the long-lived mode are metastable with respect to finite (as opposed to infinitesimal) perturbations. Transition to the explosive regime requires a perturbation amplitude above a threshold. The critical amplitude increases in time, which is expected since the background pressure profile relaxes due to dissipation, thus gradually reducing the free-energy source for the mode.

6 One dimensional model

In this section we present a simple system that exhibits linear instability, metastability and explosive behavior, which should make the nonlinear results of the previous sections more intuitive and easier to understand. Of course the results of this one dimensional model are not meant to be a quantitative explanation for the complex nonlinear behavior seen with the full MHD equations.

We start with the $1D$ equation of motion for a particle moving in a potential $\phi$ and experiencing a damping force:

[TABLE]

where $\mu$ is the damping coefficient. We choose the potential to be the quartic polynomial $\phi(\xi)=-(\xi-c_{0})(\xi-c_{1})(\xi-c_{2})(\xi-c_{3}),$ where $c_{0}=-0.48,~{}c_{1}=0.8,~{}c_{2}=2,~{}c_{3}=3$. As seen in Fig. 10 (a), the point $\xi=0$ represents a linearly unstable equilibrium (we will consider only $\xi\geq 0$ here). With any positive velocity perturbation ($v_{0}>0$) the particle will move down the potential hill, exhibiting “linear instability.” But if the damping coefficient is large enough ($\mu>\mu_{cr}(v_{0})$), it will not be able to climb out of the well and eventually settle at the metastable equilibrium point marked with ⓜ at $\xi=1.39$. This behavior is shown in Fig. 10 (b), where, after a period of exponential growth, the displacement $\xi(t)$ exhibits damped oscillations. However, with only slightly lower dissipation, the particle is able to move past the saddle point ⓢ at $\xi=2.60$ and become “explosively unstable,” as seen in panel (c). Of course the same result can be achieved by keeping the dissipation level constant but increasing the initial perturbation (panel (d)).

The damped oscillations of Fig. 10 (b) correspond to the long-lived mode described by curve (1) in Fig. 8 (a). The explosive instability in Fig. 10 (c) that develops at a lower dissipation level corresponds to the curves (2) and (3) of Fig. 8 (a). Finally the explosive instability in Fig. 10 (d) that results from a higher initial velocity corresponds to the explosive curve $1b$ in Fig. 9 that follows a larger perturbation of the initial equilibrium. Curves $2b,3b$ of Fig. 9 have not been simulated with the $1D$ model, but clearly they correspond to large-velocity perturbations of the damped oscillations in Fig. 10 (b).

7 Summary and discussion

In this work we demonstrate that the nonlinear evolution of a pressure-driven $n=1$ kink-ballooning mode can exhibit a bifurcation between a benign final state with little confinement degradation–a long-lived mode (LLM), and an explosive instability that results in a fast disruption with very short precursors. The bifurcation depends sensitively on assumed transport levels and the initial perturbation amplitude. Large diffusive transport or too small a perturbation leads to a saturated $n=1$ LLM. Equivalently, there is a transport-dependent critical perturbation amplitude, $v_{cr}=v_{cr}(\eta,\mu,\kappa_{\perp},\ldots)$, such that $v>v_{cr}$ leads to explosive behavior. The long-lived mode itself is metastable and can be pushed into the explosive regime, again with a finite perturbation above a threshold. Thus it is possible that a LLM can abruptly terminate with a fast disruption.

Since a benign LLM is a possible end state, it is clear that the initial $n=1$ instability has to be weak and not too far from an instability threshold; a robust and ideally unstable $n=1$ is unlikely to saturate without serious deleterious effects on confinement. Thus we have concentrated on weak, resistive modes far from ideal instability boundaries. This choice follows also from the expectation that an MHD mode does not come into existence as a robustly unstable ideal mode with a large growth rate; the resistive thresholds tend to be much lower and the instability generally appears first as a weak, resistive mode. But this feature (weak instability) that makes a LLM possible would at the same time seem to make it difficult to explain a fast disruption with little or no precursors, the other possible end state of the bifurcation.

This difficulty is resolved by the numerical observation, with some theoretical support (e.g., [30]), that a weakly growing $n=1$ mode can become explosive nonlinearly. Thus, a feeble resistive instability can transform into a robust ideal mode in a short period of time. Although a detailed understanding of the nonlinear process responsible for this transformation is lacking at this point, a simple one-dimensional model is presented that mimics essentially all its important features: a linear instability that can either saturate in a metastable state or lead to a nonlinear explosive instability. Of course the most important consequence of this explosive instability is that it obviates the need for a long period of mode evolution on transport time scales where it gradually becomes stronger with the changing background equilibrium. Thus a long series of precursor oscillations are avoided.

One particular feature of the pressure-driven $n=1$ mode that plays an important role in the transition to the explosive regime is the quadrupole geometry of the pressure perturbation due to the dominant $m/n=2/1$ harmonic. This perturbation naturally adds an elliptical deformation to the flux surfaces, which can nonlinearly turn into a ballooning finger. Once formed, the finger rapidly moves outward, pushing through flux surfaces while essentially maintaining its original $2/1$ symmetry. To minimize bending of the field lines as it moves into regions with $q\gg 2$, it becomes localized both in parallel and perpendicular directions. Thus, although it is originally quite extended along the field lines, it turns into a highly local nonlinear structure as it moves to the edge, transporting a significant portion of the energy in the core to a small area on the wall. Although not shown here but discussed elsewhere[31], the rapid ejection of the core is accompanied by stochastization of the field outside the finger, while the finger itself remains well-confined by regular flux surfaces. These features are consistent with jet-like flows in some high-$\beta$ disruptions in TFTR[20] and JET[23]. It is likely that the global $n=1$ mode responsible for the well-documented fast high-$\beta$ disruption in DIII-D[21, 22] also has similar origins.

Finally, ITER disruption mitigation efforts seem to be based on the anticipation that a resistive wall mode or a tearing mode would slowly lock to the wall prior to the disruption, giving at least a few tens of milliseconds of warning time. Without a perfect pressure-profile control system, fast disruptions of the type discussed in this work would make the efficacy of this approach questionable in some of the advanced scenarios planned for ITER.

8 Acknowledgements

This work was supported by MSIP, the Korean Ministry of Science, ICT and Future Planning, through the KSTAR project.