3.1 Numerical issues

The poloidal field at the X-point in a divertor configuration vanishes, leaving a saddle point in the poloidal plane and a purely toroidal field. No special behaviour appears in the vicinity of the X-point in strictly two dimensions, though with three-dimensional (3-D) perturbations, the singular behaviour is manifested as Hamiltonian chaos with a complex web of homoclinic tangles (Evans et al. Reference Evans, Moyer, Burrell, Fenstermacher, Joseph, Leonard, Osborne, Porter, Schaffer and Snyder2006). Nevertheless, the X-point is a coordinate singularity and needs to be handled properly. Even though there seems to be no new (non-ideal MHD) physics associated with the singularity, it does present numerical difficulties for quantities such as $q$ that involve surface integrals of $B_{pol}^{-1}$ . In particular, for inverse equilibrium codes based on flux coordinates, such as the CHEASE code (Lutjens et al. Reference Lutjens, Bondeson and Sauter1996), the Jacobian is usually singular (for example if an equal arclength poloidal angle is used) and the coordinate system can become highly distorted in the neighbourhood of the X-point (for straight field line coordinates). For direct equilibrium codes such as EFIT (Lao et al. Reference Lao, Ferron, Groebner, Howl, St. John, Strait and Taylor1990), the singularity is less serious but still requires care.

For stability codes based on a flux coordinate system, the singularity needs to be handled carefully. Two schemes are usually employed – both essentially remove the singular point. In the first approach, the plasma boundary is moved slightly inward so that the X-point lies outside the region considered as plasma. This cutoff procedure is required in numerical calculations that use a Fourier representation in the poloidal angle, but has traditionally been deployed in different ways in different codes – the differences typically being whether the current density and pressure in the removed region are shifted back into the plasma to keep $I$ and ${\it\beta}$ fixed, or are removed completely, and whether force balance is reimposed on the modified plasma profiles.

The MARS (Bondeson et al. Reference Bondeson, Vlad and Lutjens1992) and DCON (Glasser Reference Glasser1997) stability calculations use variations of this cutoff procedure. The equilibria for MARS are recomputed with fixed current density, pressure and boundary shape using the inverse equilibrium code CHEASE (Lutjens et al. Reference Lutjens, Bondeson and Sauter1996), whereas DCON can either use the CHEASE equilibrium or employ a simple cutoff of the edge current and pressure from the original reconstructed direct equilibrium. The resulting $q_{edge}$ is finite but the numerical value depends clearly on how much plasma is cut off. For sufficiently small cutoffs, the different procedures usually result in negligible changes in the results and a consistent answer can be obtained. For larger cutoffs, where significant current density and pressure need to be redistributed, or where a low-order rational surface is removed, the results are sensitive to the cutoff.

The alternative procedure, used in the ERATO (Gruber et al. Reference Gruber, Troyon, Rousset, Kerner and Bernard1981), GATO (Bernard et al. Reference Bernard, Helton and Moore1981) and KINX (Degtyarev et al. Reference Degtyarev, Martynov, Medvedev, Troyon, Villard and Gruber1997) codes, and also effectively in the M3D-C $^{1}$ code (Ferraro, Jardin & Snyder Reference Ferraro, Jardin and Snyder2010), utilizes a finite element representation in the poloidal angle. With a suitable choice of Jacobian (e.g. with an equal arclength poloidal angle), the distortion of the poloidal coordinate lines in the region around the X-point is minimized. Although the Jacobian is still singular at the X-point, the logarithmic singularity is weak and is manifested in surface integrands of the form $\langle F\rangle =\oint F(l)\,\text{d}l/B_{pol}$ , where $l$ is the arclength. Discrete numerical integration then yields finite numbers for these divergent integrals if the singular point where $B_{pol}\rightarrow 0$ is not a grid point. For example, in the case of the edge safety factor, $q=(f/2{\it\pi})\oint R\,\text{d}l/B_{pol},(f=RB_{{\it\varphi}})$ the numerical integration in the GATO code yields $q_{edge}\sim 5$ for the EFIT equilibrium reconstructed from discharge no. 150 513 at 2340 ms. The finite value is, in fact, set largely by how close the 2-D interpolation procedure can approximate the zero gradient in flux at the X-point rather than the fineness of the integration grid. The M3D-C $^{1}$ code also uses essentially the same method; it solves for the perturbed fields in the entire plasma plus vacuum region but avoids calculation of singular surface integrals.

It should be emphasized that the explicit cutoff and the finite (approximate) numerical integration procedures are different. The former explicitly either removes or redistributes pressure and current but the latter does not. The effect of the explicit cutoff has been investigated systematically with KINX in a comparison of the two approaches (Medvedev et al. Reference Medvedev, Hender, Sauter and Villard2001). The conclusion from that work was that the results with no cutoff can be qualitatively different from those that employed an explicit cutoff. A reasonable summary of the situation, consistent with the present and other studies, is that the cutoff matters each time it removes a significant rational surface from the edge. While there are, of course, infinitely many rational surfaces in principle, only the very low order surfaces produce a significant instability. For the low order surfaces, the difference between including the surface in the plasma and removing it is simply that the former is stable but the latter unstable to the corresponding $m/n$ external kink.

In the present case, it means simply that the $q=2$ surface cannot be removed without an expectation that it will affect the stability. From sensitivity studies with respect to variations within experimental uncertainties in the edge profiles, no unstable ideal external kink instabilities with significant $m=2$ could be found using the ideal GATO code for discharge no. 150 513 using reconstructed equilibria. From studies of different cutoffs, it was found that the $q=2$ surface must be removed in order to find the unstable ideal external $2/1$ kink. However, this would require 5 % of the plasma to be cut off. Different cutoffs, provided they were less than a tenth of a percent of the poloidal flux, yielded consistent results except in the narrow ideally unstable $q_{edge}$ ranges just below rational values. As will be discussed in § 4.4, these instabilities cannot easily explain the observed disruption. In contrast, the MARS and DCON codes utilizing a cutoff of 0.3 % provided by the CHEASE code find an edge $q$ just below 3. This might be expected to drive an edge 3/1 external kink or peeling mode but MARS and DCON, as well as GATO also find the equilibrium to be stable. Thus, the disagreement with the experimentally observed instability remains to be resolved.

The second numerical issue is only a technical complication, though also with possible consequences in practice. It applies in both the limiter and divertor configurations. As discussed in § 2, the equilibria are reconstructed so that $q_{min}>1$ is maintained. The classic internal kink with $q_{min}<1$ (Shafranov Reference Shafranov1970; Bussac et al. Reference Bussac, Edery, Pellat and Soule1976) is then generally avoided. The safety factor profile is then usually slightly non-monotonic and this tends to destabilize the ideal quasi-interchange mode (Wesson Reference Wesson1986) if $q_{min}\leqslant 1+{\it\varepsilon}\;({\it\varepsilon}\sim a/R)$ . In some of the equilibria used in the study, this internal mode is unstable. This is an independent purely $1/1$ mode but can also couple to the external kink to produce coupled $1/1$ and $2/1$ modes at finite ${\it\beta}$ . The external kink component of interest can become dominated by the internal $1/1$ displacement since the $1/1$ mode is easily excited to large amplitudes by a small energy imbalance. The calculations have been repeated for a number of different reconstructions with varying $q_{0}$ and $q_{min}$ , and the external kink modes have varying contributions from the $1/1$ component. The issue then is just interpretation; can the instability be characterized as a predominantly $2/1$ mode rather than a $1/1$ mode. The quasi-interchange can usually be eliminated by raising $q_{min}$ further. Normally, the judicious selection of an appropriate equilibrium reconstruction or time slice with sufficiently high $q_{min}$ is enough to ensure the $2/1$ component is unambiguously dominant.

The stability calculations are discussed separately for the limiter and divertor configurations in the following two sections. The ideal calculations were generally performed using the GATO code, employing no cutoff, but a numerically discrete and finite integration on the diverted surface. The ideal results were confirmed using MARS, DCON, and M3D-C $^{1}$ . Resistive MHD calculations were performed using the MARS code. Here, for the diverted case, an equilibrium recomputed from the CHEASE code was used with the boundary taken at the 99.7 % poloidal flux surface, and the value $q_{edge}\sim 2.85$ , well above 2.0, was obtained. The CHEASE equilibrium is shown overlaid in figure 4 with the original reconstructed equilibrium. Although the cutoff did not destabilize the $m/n=3/1$ ideal external kink, a significant and often dominant $m=3$ , as well as higher $m$ harmonics, were always present in the resistive MHD calculations. This will be discussed. No cutoff was required in the limiter case except to artificially reduce $q_{edge}$ slightly to a value below 2.

Several comments should be made regarding the stability calculations. The ideal calculations with the GATO code generally assumed a wall at infinity since the aim is to show that no ideal mode is unstable in the diverted case. On the other hand, the resistive MARS calculations assumed a resistive wall at the location of the DIII-D vacuum vessel in order to reproduce the instability under experimental conditions. The wall time was taken to be $2.78\times 10^{4}$ plasma Alfvén times, which is the standard value used for DIII-D. This yields an actual resistive wall time of around 9 ms. Recent careful studies, however, have found this to be a factor two longer than the value that provides the closest match to the vacuum field decay time for $n=1$ fields (Hanson et al. Reference Hanson, Bialek, Turco, King, Navratil, Strait and Turnbull2016) (the actual wall times depend weakly on the mode in question). While the actual growth rates are unimportant here, some calculations were nevertheless repeated with no wall to verify that the growth rates are at ideal levels. Finally, in all cases, the density was assumed constant for simplicity. A profile in density changes only the growth rates, with a small effect on the mode structure. Again, the magnitude of the linear growth rates is relatively unimportant here. The effect on mode structure could make a difference in the experimental diagnostic comparisons, but, again, the effect should be small.

3.2 Limiter configuration

Ideal stability calculations using the GATO code for the reconstructed limiter discharge no. 154 907 at 2425 ms with $q_{edge}=2.08$ found complete stability even with no stabilizing wall. This was confirmed with DCON and MARS. Variations in the detailed edge conditions also found stability unless the edge $q_{edge}$ was reduced below 2.0 from the best-fit value of $q_{edge}=2.08$ . This result is completely consistent with the ideal theory.

The ideal $n=1$ stability calculations were performed using GATO for the same discharge no. 154 907 at a slightly extrapolated time, where $q_{edge}=1.95$ . The extrapolation was performed by projecting the profiles to slightly higher current, while keeping the remaining data to be fitted, including the plasma boundary, unchanged. With a perfectly conducting wall on the plasma the equilibrium is unstable to an internal $m=1$ quasi-interchange mode due to the slightly hollow $q$ profile. With no wall, a pair of unstable ideal free boundary modes was found. Figure 5 shows the two modes. Plotted is the Fourier decomposition of the normal displacement $X\equiv {\it\xi}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}{\it\psi}$ with respect to the PEST straight field line angle ${\it\chi}_{P}$ , as a function of $\sqrt{V}$ , with $V$ the normalized volume contained within a given flux surface. The real and imaginary parts are different since the equilibrium is vertically asymmetric but both show $m=1$ coupled to a significant $m=2$ component with different phases. Note also the different scales between (a) and (b) and between (b) and (c); the overall amplitude for each independent eigenmode is arbitrary, however. The respective no wall growth rates are ${\it\gamma}{\it\tau}_{A}=0.00752$ and ${\it\gamma}{\it\tau}_{A}=0.001575$ ( ${\it\tau}_{A}\equiv \sqrt{B_{0}^{2}/{\it\mu}_{0}{\it\rho}R_{0}^{2}}$ is the Alfvén time, with ${\it\rho}$ the mass density). The $2/1$ external kink is destabilized when the wall is removed and couples to the quasi-interchange to produce two separate independent linear combinations of the external and internal modes. The relative contributions from the coupled $1/1$ and $2/1$ system depend on the reconstruction details but the $m/n=2/1$ kink component is the driving instability of the external mode.

Figure 5.

Ideal eigenmodes computed using GATO showing the Fourier decomposition of the normal displacement $X\equiv {\it\xi}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}{\it\Psi}$ with respect to the PEST straight field line poloidal angle ${\it\chi}_{P}$ , as a function of $\sqrt{V}$ , for the two unstable ideal modes found in the limited discharge no. 154 907 near 2425 ms. The equilibrium used is the extrapolated time shown in figure 2. Shown are (a) and (c) the real and (b) and (d) imaginary parts of the two modes. Note the normalization chosen with $\max _{{\it\Psi},m}|Re(X({\it\Psi},m))|^{2}\equiv 1$ results in different scales between the real and imaginary components for each separate eigenmode.

The presence of multiple harmonics and the inclusion of a sizeable $m=n=1$ component are inessential details for the purpose of identifying the disruption cause. Slightly different extrapolation procedures also yield different mixes of poloidal harmonics but all have an $m/n=2/1$ external kink coupled to an $m/n=1/1$ internal kink. The key point is that the predicted instability is consistent with the observed disruption to within a small margin of error.

The MARS code, running under ideal MHD conditions, reproduced the external $m/n=2/1$ kink with $q_{edge}<2$ , utilizing an equilibrium from the CHEASE code (Lutjens et al. Reference Lutjens, Bondeson and Sauter1996) that was obtained from the initial EFIT reconstruction at 2425 ms but with a slightly increased total current chosen to reduce $q_{edge}$ to $q_{edge}=1.995$ . The calculations again found a linearly unstable $2/1$ mode. In this equilibrium, the internal $m/n=1/1$ quasi-interchange was stable and only one unstable external mode was found. Figure 6(a) shows the measured signals from the external magnetic probes located at the wall at 2456 ms. A synthetic diagnostic prediction of these signals was made (Hanson et al. Reference Hanson, Bialek, Baruzzo, Bolzonella, Hyatt, Jackson, King, La Haye, Lanctot and Marrelli2014) using the result from the MARS code (Bondeson et al. Reference Bondeson, Vlad and Lutjens1992) with no wall and is replicated in figure 6(b). The observed $2/1$ structure of figure 6(a) is well reproduced.

Figure 6.

(a) Measured and (b) predicted synthetic diagnostic magnetic signals at the wall for the limited DIII-D discharge no. 154 907. The prediction using the MARS code, reproduced from (Hanson et al. Reference Hanson, Bialek, Baruzzo, Bolzonella, Hyatt, Jackson, King, La Haye, Lanctot and Marrelli2014), is from the equilibrium with $I$ increased to reduce $q_{edge}$ to $q_{edge}=1.995$ .

The conclusion for the limiter discharge, no. 154 907 is, as expected, that ideal MHD stability predicts the onset of the $2/1$ external kink mode at $q_{edge}=2.0$ , close to the observed onset at $q_{edge}=2.08$ . Analysis of similar discharges in the same series finds essentially identical behaviour. The predicted ideal kink mode structure also agrees well with the external magnetic signals for this discharge and the growth rate is consistent with an expected $2/1$ ideal MHD mode growing on the resistive wall time scale. One apparently minor point is the estimated onset actually appears to be slightly above $q_{edge}=2.0$ , outside the estimated experimental uncertainty (Hanson et al. Reference Hanson, Bialek, Baruzzo, Bolzonella, Hyatt, Jackson, King, La Haye, Lanctot and Marrelli2014; Martin et al. Reference Martin, Baruzzo, Bialek, Bolzonella, Bonfiglio, Cavazzana, Chacon, Ferraro, Garofalo and Hanson2015). While unaccounted for systematic errors could be present, it does hint that some non-ideal effects might also be important in the limiter cases. Nonetheless, the ideal kink is predicted to be unstable very close to the time when the mode is detected experimentally.

3.3 Divertor configuration

Ideal MHD stability calculations for the diverted discharge no. 150 513 at 2340 ms found complete stability. Different extrapolation and cutoff procedures yielded different values for $q_{edge}$ . The GATO calculations used the EFIT equilibrium with no cutoff but with a numerically calculated edge $q$ of $q_{edge}\sim 5$ . The DCON and MARS codes used the equilibrium recomputed as an inverse equilibrium from CHEASE with the cutoff of 0.3 % in poloidal flux, yielding $q_{edge}=2.85$ , but confirmed the complete stability to external kink modes. An ideal quasi-interchange, however, is unstable that, although irrelevant, does result in complications that will be discussed. The M3D-C $^{1}$ code also used no cutoff and extended the calculation smoothly through the vacuum, with a sharp resistivity transition at the separatrix. This also found complete linear stability. Overall, no external kink instabilities were found with moderate variations in the edge profiles as long as $q_{edge}>2$ and $q_{edge}$ was kept from the narrow peeling instability ranges associated with the higher rational surfaces. Obtaining $q_{edge}\leqslant 2$ requires an unjustifiable 5 % cutoff.

Ideal $n=1$ instabilities do appear in small ranges of $q_{edge}$ near rational values, $(m-1)<q_{edge}^{crit}<q_{edge}<5$ , as predicted (Shafranov Reference Shafranov1970). However, the instabilities for $m\geqslant 3$ consist of a single almost pure $m$ , are highly localized, and have very low growth rates. An example is shown in figure 7 for the $m=3$ peeling instability obtained directly from extrapolating the profiles to 2350 ms and applying a small cutoff, for which $q_{edge}=2.95$ . The almost pure $m=3$ mode is localized strongly to the last few surfaces. It is also only very weakly unstable with ${\it\gamma}{\it\tau}_{A}=8\times 10^{-4}$ with no wall. As will be discussed later in § 4.4, the narrow ranges of plasma current where the mode is linearly unstable are difficult to reconcile with the observed instability time development.

Figure 7.

Ideal peeling eigenmode computed using GATO showing the Fourier decomposition of the normal displacement $X\equiv {\it\xi}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}{\it\Psi}$ with respect to ${\it\chi}_{P}$ , as a function of $\sqrt{V}$ for the diverted discharge no. 150 513 near 2340 ms, using the equilibrium obtained from extrapolation to higher current with a small cutoff with $q_{edge}=2.97$ . Shown are (a) and (b) the real and imaginary parts of the mode, normalized with $\max _{{\it\Psi},m}|Re(X({\it\Psi},m))|^{2}\equiv 1$ .

The simplest explanation is to invoke non-ideal effects. The issue here is that non-ideal effects are generally small perturbations on ideal MHD and this seems inconsistent with the observed large growth rates of the instabilities. Finite resistivity offers the most promising option. To obtain sufficiently large growth rates, the resistivity cannot be a small perturbation as it usually is in the core of the plasma. However, the edge does have large resistivity in a narrow region inside the last closed flux surface. The hypothesis then is that the resistivity only needs to be sufficiently large in the edge region near $q=2$ in order to produce a fast growing, edge-peaked external kink. The MARS code offers an opportunity to explore this.

The resistivity profile is not generally measured in tokamaks. Instead, the resistivity is usually either taken from a model such as the Spitzer formula (Spitzer Reference Spitzer1962) or later improvements to include neoclassical effects (Sauter, Angioni & Lin-Liu Reference Sauter, Angioni and Lin-Liu1999), or is inferred from transport simulations (St. John et al. Reference St. John, Taylor, Lin-Liu and Turnbull1994). The Spitzer model is implemented in the MARS code as

where ${\it\eta}_{0}=S^{-1}$ is an input constant, ( $S\equiv {\it\tau}_{R}/{\it\tau}_{A}$ , the ratio of the resistive diffusion time to the Alfvén time, is the magnetic Reynolds, or Lundquist number), $T_{e}$ is the local electron temperature and $T_{0}\equiv T_{e}(0)$ . $T_{e}$ is measured from Thomson scattering and the profile is shown in figure 8(a). The data in the core from this time are poorly calibrated and ignored but the edge temperature data are reliable; a detailed discussion of the data quality is provided in § 4.4.1 in the context of identifying possible resistive tearing modes in the plasma. The corresponding normalized resistivity profile is shown in figure 8(b). Also overlaid in figure 8(b) are the Sauter and ‘effective’ profiles to be discussed later.

Figure 8.

(a) Measured electron temperature profile for the diverted discharge no. 150 513 at 2340 ms from Thomson scattering. The radial variable is the normalized square root of the toroidal flux. (b) The corresponding Spitzer resistivity profile from (3.1) (solid curve). Also shown in (b) are the model resistivity profiles obtained from the Sauter formula (dashed curve) in (3.3) and the ‘effective’ resistivity (dotted curve) obtained from transport modelling (3.4). These are shown on both linear and log scales against $\sqrt{V}$ .

A cap, ${\it\eta}_{max}$ , to the edge value, traditionally set at 100 times the core value, is applied in the code. This can crudely model some neoclassical corrections (Sauter et al. Reference Sauter, Angioni and Lin-Liu1999) and avoids numerically localized singular spikes right on the last surface. These are unphysical and are generated by the several decade variation in ${\it\eta}$ over the last few grid points as ${\it\eta}\rightarrow \infty \;(T_{e}\rightarrow 0)$ . For our purposes, the default cap is often too low and so was varied over a range, determined only by when physically real modes were obtained.

Figure 9.

Resistive kink eigenmodes computed using MARS for the diverted DIII-D discharge no. 150 513 showing the Fourier decomposition of the normal displacement $X\equiv {\it\xi}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}{\it\Psi}$ with respect to ${\it\chi}_{P}$ , as a function of $\sqrt{V}$ , for the case with the Spitzer profile with (a) and (b) ${\it\eta}_{0}=10^{-4}$ and (c) and (d) ${\it\eta}_{0}=10^{-5}$ . The equilibrium used is the CHEASE equilibrium in figure 4. Shown are (a) and (c) the real and (b) and (d) imaginary parts for the two cases, normalized with $\max _{{\it\Psi},m}|Re(X({\it\Psi},m))|^{2}\equiv 1$ .

In addition, numerical issues in MARS require that ${\it\eta}\rightarrow 0$ at the magnetic axis. This has no effect on the results here since the aim is to find an instability from the large edge resistivity. In fact, with later modifications to the profile to be discussed subsequently, the resistivity was sometimes forced to vanish faster in the core than in (3.1). Finally, additional options exist in the code to crudely model the effect of trapped particles, as contained in the Sauter model. However, except for the cap, these were ignored and the Sauter model was obtained independently from ONETWO transport code modelling (Pfeiffer et al. Reference Pfeiffer, Davidson, Miller and Waltz1980; St. John et al. Reference St. John, Taylor, Lin-Liu and Turnbull1994).

The stability calculations performed for discharge no. 150 513 for the ideal case with the MARS code were repeated using the Spitzer resistivity profile in (3.1) with varying ${\it\eta}_{0}$ . For ${\it\eta}_{0}<10^{-5}$ , no instability was found. For ${\it\eta}_{0}=10^{-4}$ and ${\it\eta}_{0}=10^{-5}$ , external kink-like solutions were obtained and are shown in figure 9. These were obtained with caps of ${\it\eta}_{max}=20{\it\eta}_{0}$ for ${\it\eta}_{0}=10^{-4}$ and ${\it\eta}_{max}=100{\it\eta}_{0}$ for ${\it\eta}_{0}=10^{-5}$ . The real and imaginary parts of the poloidal Fourier decomposition of each of the modes are shown. In both cases the $m=2$ component is large and is maximized at the $q=2$ surface for both the real and imaginary parts. There is also a significant $m=3$ component peaked at the edge and some $m=1$ is also present in the core. In figure 9(c,d), the singular higher $m$ are also present at the edge. These are eliminated with a lower cap $({\it\eta}_{max}<50{\it\eta}_{0})$ and a corresponding small reduction in growth rate. The resistive kink modes have a finite real frequency due to the resistive wall and resistive plasma. For ${\it\eta}_{0}=10^{-4}$ the complex growth rate, normalized to an inverse Alfvén time, ${\it\tau}_{A}^{-1}$ , is $({\it\gamma}+\text{i}{\it\omega}){\it\tau}_{A}=0.0844-0.00928\text{i}$ , and for ${\it\eta}_{0}=10^{-5}$ , $({\it\gamma}+\text{i}{\it\omega}){\it\tau}_{A}=0.0568+0.00203\text{i}$ . With no wall, these become respectively $({\it\gamma}+\text{i}{\it\omega}){\it\tau}_{A}=0.1371-0.00429\text{i}$ , and $({\it\gamma}+\text{i}{\it\omega}){\it\tau}_{A}=0.1342-0.00360\text{i}$ . The mode structures are largely unchanged.

These are ideal kink level growth rates and the instabilities are kink-like with no sudden phase changes, characteristic of tearing, at the rational $q=2$ surface. Note also that the resistive wall does slow the modes but the growth rates remain at ideal levels. For now, they are referred to as resistive kink modes, since, even though the scaling with ${\it\eta}$ is not determined at this point, they clearly exist only because of the finite resistivity. The scaling with key characteristic resistivity values will be considered later in § 4.

One point of possible contention should be addressed before proceeding. The modes in figure 9 have an $m=2$ component that is not peaked at the edge, but in fact turns over at the $q=2$ surface. This is in contrast to the ideal external kink mode calculated for the limiter discharge no. 154 907 where the $m=2$ amplitude increases fairly monotonically toward the edge (figure 5). The $m=3$ harmonic is peaked at the edge in figure 9 and, according to the usual experimental categorization of selecting the largest poloidal harmonic in the diagnostic signals, one might be tempted to classify the instabilities as $m=3$ rather than $m=2$ external kinks. However, the $m=2$ harmonic must peak near the $q=2$ surface. This is the signature of kink or interchange parity. It only peaks at the edge when the $q=2$ surface lies in the vacuum. In the diverted case there are always corresponding higher $m$ components up to $m=n([q_{edge}]+1)$ , where $[q_{edge}]$ is the integer part of $q_{edge}$ . Variations in the profile parameters ${\it\eta}_{q=2}\equiv {\it\eta}\;(q=2)$ and ${\it\eta}_{max}$ find considerable sensitivity of the ratio of the $m=2$ to $m=3$ peaks. Note that the $m=2$ peak is dominant at least for the case with ${\it\eta}_{0}=10^{-4}$ . The identification of the mode with a given single $m$ can make qualitative sense, however, when comparing the signals measured on the external diagnostics. This will be considered later in this section.

Efforts to find resistive external kink instabilities for this equilibrium at lower core resistivity, ${\it\eta}_{0}<10^{-6}$ , with the Spitzer profile were unsuccessful. Keeping the cap ratio ${\it\eta}_{max}/{\it\eta}_{0}$ fixed while ${\it\eta}_{0}$ is reduced tended to result in only stable predictions, whereas lower ${\it\eta}_{0}$ , with a fixed cap value ${\it\eta}_{max}$ , tended to yield only internal $m=1$ quasi-interchange-like modes. For sufficiently low ${\it\eta}_{0}$ , increasing ${\it\eta}_{max}$ tends to produce the spurious singular edge ‘peeling-like’ modes, usually $m\geqslant 3$ , and localized on the last one or two surfaces at the edge. This is known behaviour for the MARS code and is one reason the cap is imposed; these spurious $m\geqslant 3$ are even more singular at the edge than the weakly unstable ideal peeling mode in figure 7.

Given that no external modes were found for ${\it\eta}_{0}<10^{-5}$ , in the context of the hypothesis that the resistivity in the edge region destabilizes the mode, the resistivity at $q=2$ is expected to largely determine the mode onset condition. From the calculations, this onset value is found to be of the order of ${\it\eta}^{crit}(q=2)\equiv {\it\eta}_{q=2}^{crit}\sim 10^{-4}$ in dimensionless inverse Lundquist number units. For the Spitzer profile, ${\it\eta}_{q=2}\sim 50{\it\eta}_{0}$ . To reach values of ${\it\eta}_{q=2}>{\it\eta}_{q=2}^{crit}$ using the Spitzer profile, the ${\it\eta}_{0}$ values required ${\it\eta}_{0}\sim {\it\eta}_{q=2}/50>0.2\times 10^{-5}$ , are then much larger than the experimentally expected core resistivity values.

The program to explain the observed instability as a resistive kink, requires ${\it\eta}_{0}$ to be reduced to more realistic values of order ${\it\eta}_{0}\sim 10^{-6}$ , to ${\it\eta}_{0}\sim 10^{-8}$ or so, while maintaining sufficiently high ${\it\eta}$ near $q=2$ . Two strategies were tried. First, the profile in ${\it\eta}$ was modified to drop faster at the axis, making it essentially zero in most of the core. Second, a bump was introduced in the profile at $q=2$ to maintain ${\it\eta}_{q=2}>{\it\eta}_{q=2}^{crit}$ as ${\it\eta}_{0}$ is reduced. In addition, as mentioned, the edge resistivity cap ratio ${\it\eta}_{max}/{\it\eta}_{0}$ was increased as ${\it\eta}_{0}$ was reduced to maintain ${\it\eta}_{max}>{\it\eta}_{q=2}^{crit}$ , while avoiding spurious singular modes right at the edge. Several justifications will be given for this procedure in § 4.

The resistivity profile was accordingly modified from the Spitzer form in (3.1) by including two form factors, $s^{{\it\alpha}}$ and $f(s)$ , so that

with

Choosing ${\it\alpha}\geqslant 1$ forces the profile ${\it\eta}(s)$ to vanish quickly in the core. Here, $s=\sqrt{{\it\psi}}$ , with ${\it\psi}$ the poloidal flux, is the radial coordinate employed in the MARS code. The Gaussian bump $f(s)$ , on the other hand, allows an enhancement to ${\it\eta}(s)$ near $u(s)=u_{\ast }\equiv u(s_{\ast })$ by a factor $1+A_{0}$ with width $w$ . The parameter $u(s)$ is simply a stretching coordinate, $u(s)=as^{2}+(1-a)s$ , with $a$ chosen to squeeze the Gaussian more rationally within the steep edge region. Then $s_{\ast }$ is usually taken in the vicinity of $q=2$ or $q_{95}$ .

Figure 10 shows the modified profile (dashed blue curve) for ${\it\eta}_{0}=10^{-5}$ with ${\it\alpha}=0$ , $A_{0}=1.5$ , $w=0.25$ , $s_{\ast }=0.975$ (i.e. at $q=q_{95}$ ), and ${\it\eta}_{max}=100{\it\eta}_{0}$ . This is shown in comparison to the unmodified Spitzer profile with ${\it\eta}_{0}=10^{-5}$ (solid blue curve). The major effect of the Gaussian bump is to move the steep edge region inward; a quite large change in the local magnitude of ${\it\eta}(s)$ corresponds to only a small radial shift inward.

Figure 10.

Spitzer profiles (solid curves) for the cases with ${\it\eta}_{0}=10^{-5}$ and ${\it\eta}_{max}=100{\it\eta}_{0}$ (blue curves; left scale) and ${\it\eta}_{0}=10^{-6}$ , and ${\it\eta}_{max}=200{\it\eta}_{0}$ (red curves; right scale). Modified resistivity profiles (dashed curves) for the same two cases with modification parameters ${\it\alpha}=0,A_{0}=1.5$ , $w=0.25$ and $s_{\ast }=0.975$ , and ${\it\alpha}=1$ , $A_{0}=5.0$ , $w=0.50$ and $s_{\ast }=0.985$ , respectively.

Figure 11.

Resistive kink eigenmodes computed using MARS for the diverted DIII-D discharge no. 150 513 showing the Fourier decomposition of the normal displacement $X\equiv {\it\xi}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}{\it\Psi}$ with respect to ${\it\chi}_{P}$ , as a function of $\sqrt{V}$ , for the case with the two modified Spitzer profiles shown in figure 10. The equilibrium used is the equilibrium recomputed from CHEASE in figure 4. Shown are (a) and (c) the real and (b) and (d) the imaginary parts for the two cases, normalized with $\max _{{\it\Psi},m}|Re(X({\it\Psi},m))|^{2}\equiv 1$ .

The results of a calculation with ${\it\eta}_{0}=10^{-5}$ but including this moderate enhancement at $q=2$ are shown in figure 11(a,b). Comparing with figure 9(c,d), the Gaussian bump has significantly increased the relative amplitude of the $m=2$ peak. The growth rate is also increased by approximately 50 % to $({\it\gamma}+\text{i}{\it\omega}){\it\tau}_{A}=0.0298-0.00065\text{i}$ . Figure 11(c,d) shows a similar successful calculation for ${\it\eta}_{0}=10^{-6}$ , in this case with ${\it\alpha}=1$ , $A_{0}=5.0$ , $w=0.50$ , $s_{\ast }=0.985$ and ${\it\eta}_{max}=200{\it\eta}_{0}$ . With a resistive wall, the growth rate is $({\it\gamma}+\text{i}{\it\omega}){\it\tau}_{A}=0.0101-0.00161\text{i}$ . Here, there is a large $m=1$ in the core in addition to the $m=2$ peak and the $m=3$ ‘peeling’ component. The parameters for the Gaussian modification were selected to minimize this $m=1$ component and maximize the $m=2$ harmonic with a moderate value for the amplitude $A_{0}$ , and the modified resistivity profile is shown as the dotted red curve in figure 10. For ${\it\eta}_{0}\leqslant 10^{-6}$ , the $m=1$ component was present in most of the calculations that revealed an unstable numerically resolved instability. Taking ${\it\alpha}=1$ tended to reduce but not eliminate the $m=1$ component.

The calculations were repeated using the Sauter model for the resistivity profile and for an ‘effective’ resistivity profile obtained from transport simulations for this discharge. The Sauter formula (Sauter et al. Reference Sauter, Angioni and Lin-Liu1999) is given by

with the effective trapped particle fraction given in terms of the charge state $Z_{eff}$ and the actual trapped particle fraction $f_{t}$ by

The ‘effective’ resistivity is calculated using the ONETWO $1(1/2)D$ transport code (Pfeiffer et al. Reference Pfeiffer, Davidson, Miller and Waltz1980; St. John et al. Reference St. John, Taylor, Lin-Liu and Turnbull1994) from the Ohm’s law relation between the parallel electric field $E_{\Vert }$ and current density $j_{\Vert }$ ,

where $E_{\Vert }$ is obtained from Faraday’s law using successive equilibrium reconstructions at two different times

Here, $R$ is the major radius, $V_{loop}$ is the loop voltage and $E_{{\it\phi}}=(B_{{\it\phi}}/B)E_{\Vert }$ is the toroidal electric field ( $B$ and $B_{{\it\phi}}$ are the total magnetic field and toroidal component, respectively). The parallel current $j_{\Vert }$ is then extracted from the equilibrium reconstruction for the total current density at an intermediate time, after subtracting the non-Ohmic contributions.

The unmodified profiles for the three models are compared in figure 8(b). These are shown on both linear and log scales against the square root of the volume, which is close to the toroidal flux variable in figure 8(a). The Sauter and ‘effective’ models include a pedestal at the inside foot of the steep gradient region that is not present in the simpler Spitzer model. The ‘effective’ profile is a little broader but all three are similar. While the ‘effective’ profile is subject to fairly large uncertainties since it is derived from a time difference of two full equilibrium reconstructions, the profile is nevertheless consistent with the others.

MARS calculations for the divertor discharge no. 150 513 using the unmodified Sauter and ‘effective’ profiles from (3.3) and (3.4) yield weakly unstable and poorly resolved resistive kink instabilities. Only a small enhancement following (3.2) is needed to obtain clear, well-resolved $m=2$ kink-like unstable modes in this case. Figure 12 shows the resulting instability for one of the modified Sauter profiles, chosen because it revealed the clearest $m=2$ solution with the least enhancement, namely with ${\it\alpha}=0$ , $A_{0}=2.25$ , $s_{\ast }=0.98$ , $w=0.30$ and ${\it\eta}_{max}=180{\it\eta}_{0}$ . For the Sauter profile, ${\it\eta}_{0}=1.1\times 10^{-6}$ . The complex growth rate with the wall in this case is $({\it\gamma}+\text{i}{\it\omega}){\it\tau}_{A}=0.0631-0.00385\text{i}$ and with no wall $({\it\gamma}+\text{i}{\it\omega}){\it\tau}_{A}=0.1309-0.00462\text{i}$ . The $m=2$ structure is clear, again with no sign change associated with tearing at the $q=2$ surface. The inset to figure 12 shows the resistivity profile used, overlaid with the original Sauter profile. The modification is visible in this case as an actual bump just outside $q=2$ but well inside the edge. The main effect, as in the Spitzer case, however, is to simply move the steep gradient inward.

Figure 12.

Resistive kink eigenmodes computed using MARS for the diverted DIII-D discharge no. 150 513 showing the Fourier decomposition of the normal displacement $X\equiv {\it\xi}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}{\it\Psi}$ with respect to ${\it\chi}_{P}$ , as a function of $\sqrt{V}$ , for the case with the Sauter profile modified with an additional enhancement near $q=2$ using ${\it\alpha}=0,A_{0}=2.25,w=0.30$ , $s_{\ast }=0.98$ , ${\it\eta}_{0}=1.1\times 10^{-6}$ and ${\it\eta}_{max}=180{\it\eta}_{0}$ . The profile is shown in the inset overlaid with the original Sauter profile. The equilibrium used is the equilibrium recomputed from CHEASE in figure 4. Shown are (a) the real and (b) the imaginary parts, normalized with $\max _{{\it\Psi},m}|Re(X({\it\Psi},m))|^{2}\equiv 1$ .

The calculations were repeated for two of the modified ‘effective’ resistivity profiles. Again, the enhancements needed were relatively small, at least in terms of the shift in the steep gradient region. Figure 13 shows the results. In figure 13(a,b), the modifications to the effective profile were taken as ${\it\alpha}=0$ , $A=10$ , $s_{\ast }=0.975$ , $w=0.30$ and ${\it\eta}_{max}=10^{5}{\it\eta}_{0}$ . From the transport modelling, ${\it\eta}_{0}=3.4\times 10^{-7}$ ; the resistivity profile is shown in the inset to figure 13(a). This was again selected as having the clearest $m=2$ structure and in this case, $({\it\gamma}+\text{i}{\it\omega}){\it\tau}_{A}=0.0796-0.00422\text{i}$ . With no wall, $({\it\gamma}+\text{i}{\it\omega}){\it\tau}_{A}=0.1602-0.00428\text{i}$ . A second case, shown in figure 13(c,d), with ${\it\alpha}=0$ , $A=2$ , $s_{\ast }=0.98$ , $w=0.30$ and ${\it\eta}_{max}=1.2\times 10^{3}{\it\eta}_{0}$ , was chosen as having the smallest modification but a clear $2/1$ resistive kink instability with $({\it\gamma}+\text{i}{\it\omega}){\it\tau}_{A}=0.0513-0.00392\text{i}$ for the resistive wall calculation and $({\it\gamma}+\text{i}{\it\omega}){\it\tau}_{A}=0.1041-0.00485\text{i}$ for the no wall case. In both, the structure is well resolved and quite similar to the Sauter case in figure 12.

Figure 13.

Resistive kink eigenmodes computed using MARS for the diverted DIII-D discharge no. 150 513 showing the Fourier decomposition of the normal displacement $X\equiv {\it\xi}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}{\it\Psi}$ with respect to ${\it\chi}_{P}$ , as a function of $\sqrt{V}$ , for the case with the ‘effective’ profile modified with additional enhancements near $q=2$ using (a) and (b) ${\it\alpha}=0,A_{0}=10$ , $w=0.30,s_{\ast }=0.975,{\it\eta}_{0}=3.4\times 10^{-7}$ , and ${\it\eta}_{max}=10^{5}{\it\eta}_{0}$ and (c) and (d) ${\it\alpha}=0,A_{0}=2$ , $w=0.30,s_{\ast }=0.980$ , ${\it\eta}_{0}=3.4\times 10^{-7}$ and ${\it\eta}_{max}=1.2\times 10^{3}{\it\eta}_{0}$ . The profiles are shown in the inset overlaid with the original ‘effective’ resistivity profile. The equilibrium used is the equilibrium recomputed from CHEASE shown in figure 4. Shown are (a) and (c) the real and (b) and (d) the imaginary parts, normalized with $\max _{{\it\Psi},m}|Re(X({\it\Psi},m))|^{2}\equiv 1$ .

Figure 14 shows the measured magnetic signals obtained from the magnetic diagnostics on DIII-D for this discharge at the onset of the $2/1$ instability. Unlike discharge no. 154 907 (figure 6 a), for discharge no. 150 513, the full internal poloidal and toroidal array of poloidal magnetic field ( $B_{pol}$ ) probes was not installed. However, data were available from three poloidal arrays at different toroidal locations but limited to the outboard side, as well as a single array at one location with full poloidal coverage and a single toroidal array around the outboard midplane. From the toroidal array, which shows almost entirely an $n=1$ structure, and using differences to remove the $n=0$ component, the signal can be reconstructed over the entire vessel assuming the poloidal arrays also the measure the same $n=1$ mode. In addition, outside the conducting vessel, two full poloidal coverage saddle-loop arrays at different toroidal locations and a separate toroidal array on the outboard midplane were available. These measure the normal field (usually referred to as the radial field, $B_{rad}$ ) penetrating the vessel. Figure 14(a) shows the contours of the perturbed poloidal field at the wall location in poloidal–toroidal angle space reconstructed from the internal arrays. With respect to the geometrical angle subtended from the vessel centre, this shows a mixed $m=2$ and $m=3$ structure (figure 14 b); the Fourier decomposition using this data indicates that the $m=2$ is dominant; the $m=2$ comprises 54 %, whereas the $m=3$ comprises 23 % with most of the rest $m=1$ and $m=4$ . Except for the $90^{\circ }$ phase shift between $B_{pol}$ and $B_{radial}$ , the radial field from the saddle loops outside the vessel is similar, although there, the $m=3$ is slightly more dominant; the $m=2$ comprises 31 %, whereas the $m=3$ comprises 38 % with most of the rest $m=1$ .

Figure 14.

Measured (a) poloidal magnetic field signals at the wall for diverted DIII-D discharge no. 150 513 reconstructed from a single full coverage array at the $320^{\circ }$ toroidal angle and the outboard poloidal arrays, and assuming an $n=1$ mode. The probe locations are indicated. (b) Poloidal Fourier decomposition of loop signals.

Synthetic diagnostic predictions for these signals are shown in figure 15(a–c) for the Sauter case of figure 12 and the two ‘effective’ resistivity cases of figure 13. The predictions from all three cases are similar and match well with the measured data in figure 14. Note that the Fourier content in figure 14(b) is for the poloidal field and is not taken with a straight field line angle so is not directly comparable to the displacement mode structures in figures 9 and 11–13. Nevertheless, the overall $m=2,3$ structure is readily identifiable. The relative contributions of the computed mode structures within the plasma from $m=1,2,3$ and higher $m$ are sensitive to significant changes in the resistivity profile. Generally, it is found that, consistent with naïve expectations, increasing ${\it\eta}_{q=2}$ tends to increase the $m=2$ harmonic, increasing ${\it\eta}_{max}$ tends to increase the higher $m$ harmonics and increasing ${\it\eta}_{0}$ tends to increase the $m=1$ core component. For sufficiently high ${\it\eta}_{0}$ and low ${\it\eta}_{q=2}$ , the computed instabilities are internal, resistive quasi-interchange modes and are discarded from consideration. At very low ${\it\eta}_{0}$ , the irrelevant ideal quasi-interchange instability is found with very low growth rate. This independent mode is usually masked by the faster growing resistive kink modes when the latter are unstable and is also ignored.

Figure 15.

Synthetic diagnostic prediction for the magnetic signals at the wall for the diverted DIII-D discharge no. 150 513 from the three resistive kink modes in figures 12 and 13.

For very large and steep profiles, the physical resistive kink is lost and a numerical artefact only is present, consisting of $m>3$ harmonics localized almost entirely on the last surface. These modes are also discarded as unphysical. In the intermediate range, some modes consist of a well-resolved $m\leqslant 3$ part and a poorly resolved edge. While not always accurately resolved, these cases are included in the scaling studies discussed in the next section. In all other cases, the modes display the general features shown in figures 9, 11 and 12, with some $m=1$ , significant $m=2$ peaking at $q=2$ (i.e. with kink-like parity) but remaining finite at the edge, and with an $m=3$ peaked at the edge.

The resistive kink is then a good candidate to explain the observed instabilities in a diverted cross-section. While it has not been shown that the resistive kink is necessarily unstable, it has been shown that the mode can be found within the considerable experimental uncertainties in the resistivity profile and that the computed instabilities reproduce the observed magnetic signals. The hypothesis also makes sense in that at some point, for sufficiently large resistivity at $q=2$ , the plasma edge region should be indistinguishable from the vacuum on the other side of the separatrix.

5.1 Future directions

The results open new avenues that can be pursued in order to resolve the remaining concerns. Most important, the analysis in § 3 was performed for $L$ -mode plasmas since experiments with good data were available for those cases and the comparison could be made directly with the limited $L$ -mode discharge. The edge conditions in $H$ -mode discharges, with a large pedestal pressure and associated bootstrap current, are significantly different from those in $L$ -mode. In particular, the temperatures near the 95 % flux surface tend to be higher and the resistivity correspondingly lower. At some point close enough to the edge, the resistive kink should still be destabilized and it needs to be shown that this corresponds to $q_{95}=2$ .

Some confidence in this respect can be found from the fact that one sequence shown in figure 19 used a Spitzer resistivity profile taken from a previously analysed $H$ -mode discharge. For this sequence, the width of the bump was generally larger but the enhancement factors were similar to those using the Sauter and ‘effective’ profiles. Despite the different conditions, it is significant that these data are consistent with the rest of the data, which correspond to an $L$ -mode resistive edge. It can then be expected then that the quantitative results in $H$ -mode would be similar to that reported here for the $L$ -mode case since the same physics is applicable. Detailed analysis of the $H$ -mode discharges in the same sequence as the $L$ -mode cases analysed here should confirm this.

The growth rate scaling studies have revealed some additional surprises. First, at both low and high ${\it\eta}$ , the diverted and limited cases appear to scale somewhat differently. The reason for this is not clear; it may be an artefact of the finite data set or due to the influence of the nearby X-point in the diverted case. In general, the transitions between different fractional scalings are unexpected since there seems to be no definite transition in mode structure in either discharge type. This is in contrast to the transitions between ${\it\nu}=3/5$ and ${\it\nu}=1/3$ for the internal resistive modes that coincide with a distinct change from tearing to interchange parity (Hastie et al. Reference Hastie, Hender, Carreras, Charlton and Holmes1987; Charlton et al. Reference Charlton, Carreras, Holmes and Lynch1988; Huysmans et al. Reference Huysmans, Goedbloed and Kerner1993). For the diverted case the fractional ${\it\nu}=1/2$ scaling is unusual. Given the scatter in the data, a ${\it\nu}=3/5$ scaling fits almost as well for the lowest ${\it\eta}$ values (figure 19 dashed green line) but the modes are interchange or kink parity.

The smooth transitions in scaling exponent ${\it\nu}$ for the limiter case are especially clear from the continuous sequence in figure 20 with increasing bump amplitude. It is interesting to speculate that the resistive external kink modes might generally have a scaling that transitions smoothly through various fractional values, ${\it\nu}=2/3$ , ${\it\nu}=3/5$ , ${\it\nu}=1/2$ and ${\it\nu}=1/3$ , to the ideal scaling with ${\it\nu}=0$ as ${\it\eta}$ increases from low to vacuum-like values. While the linear regime may not always be present (for example in diverted equilibria), the ${\it\nu}=0$ scaling transition may be present in the divertor case. Such a continuous transition, however, would question the validity of fractional power-law scalings for the external resistive kink modes over more than a limited range of values. An analytic model that includes the sharp increase in edge resistivity is probably needed to resolve the scaling issues.

The results also have a bearing on the conventional cutoff procedures used in ideal stability calculations. The conventional procedure used in ideal MHD codes, as defined from earlier studies (Medvedev et al. Reference Medvedev, Hender, Sauter and Villard2001), and set out in § 3.1 is largely validated; small cutoffs are valid so long as they do not remove the key rational surfaces and generate localized peeling modes such as in figure 7. However, this can now be refined in light of the results presented here. Given that the actual instability is a resistive kink, an ideal code can still reproduce an instability that exhibits most of the key features of the resistive kink. If the discharge is unstable, an ideal stability code can find an equivalent ideal instability only if the respective rational surface is removed; the cutoff of the rational surface is necessary to ‘reproduce’ the instability as an ideal mode. Presumably in that case ${\it\eta}_{q=2}\geqslant {\it\eta}^{crit}$ . On the other hand, if the discharge is actually stable, then the cutoff should not remove the rational surface since otherwise an (ideal) instability will be found. In this case, the resistive kink is stable in reality so that ${\it\eta}_{q=2}<{\it\eta}^{crit}$ and the cutoff should not exceed the surface where ${\it\eta}={\it\eta}^{crit}$ .

The standard procedure requires that the actual stability be known beforehand so that the appropriate surfaces are removed. The new alternative is that if ${\it\eta}^{crit}$ is known, the cutoff can be imposed there. For high and intermediate $n$ peeling-ballooning stability calculations, where the pure high $m$ peeling modes are considered irrelevant, the standard practice of using only small cutoffs with a minimum distance to the next rational surface (Eldon et al. Reference Eldon, Boivin, Groebner, Osborne, Snyder, Turnbull, Tynan, Boedo, Burrell and Kolemen2015) is probably justified. For the lower $m$ modes like the $2/1$ , the instabilities are highly relevant and the standard procedure should be replaced. Resistive kink calculations can be used to quantify the ${\it\eta}^{crit}$ values for these lower $n$ modes, making the cutoff selection for ideal codes quantitative.

The question then remains to determine ${\it\eta}^{crit}$ . From the calculations, a rough critical value for ${\it\eta}_{q=2}$ appears to exist, above which the $m=2$ resistive kink is unstable and below which it is stable. The value is in the range ${\it\eta}_{q=2}\sim 10^{-4}$ , but also depends on details of the profile in the neighbourhood of $q=2$ . The region just inside $q=2$ appears to be the most crucial. The appropriate value presumably depends on the instability in question.

Additionally, the Glasser stabilization effect (Glasser, Greene & Johnson Reference Glasser, Greene and Johnson1975) probably imposes a dependence of ${\it\eta}^{crit}$ on the plasma ${\it\beta}$ , though for the low pressure equilibria considered here with ${\it\beta}_{p}\sim 0.1$ , this is probably small. Previous studies also found this effect to be much weakened for instabilities closer to the edge. The effect is possibly important in $H$ -mode, however, since there is a significant pressure gradient in the region near $q=2$ .

Finally, in the RFX-mod experiments, the closed-loop feedback was turned off at different times in the evolution and the open-loop growth rate trend (Hanson et al. Reference Hanson, Bialek, Baruzzo, Bolzonella, Hyatt, Jackson, King, La Haye, Lanctot and Marrelli2014; Piovesan et al. Reference Piovesan, Hanson, Martin, Navratil, Turco, Bialek, La Haye, Lanctot, Okabayashi and Paz-Soldan2014) strongly suggested a stable region exists below $q_{edge}\sim 1.5$ . This is consistent with the original picture (Shafranov Reference Shafranov1970; Wesson Reference Wesson1978), which found a region stable to ideal kink modes when $1.0\leqslant q_{edge}\leqslant q_{crit}\sim 1.5$ for sufficiently diffuse current profiles. This appears to be a promising regime, allowing very high current and with no low $n$ rational surfaces in the plasma. The major obstacle to reaching it at present is the lack of sufficient power supplies to maintain active stabilization until the stable region is reached. Given sufficient power supplies in the future, an obvious concern would then be whether the ideally stable region also allows stability to the resistive kink. The answer seems likely to depend on the resistivity profile. Resistive kink stability could be controlled partially by varying the edge resistivity, either through heat addition or removal, influx of neutral particles or by modifying the neoclassical corrections through the current and pressure profiles.

Probably a more serious concern is the internal quasi-interchange, since this is unstable typically when $q_{min}\sim 1+{\it\varepsilon}$ , which would reduce the stable $q$ profile range to between $1+{\it\varepsilon}$ in the core and about 1.4 at the edge. For most scenarios in DIII-D, ${\it\varepsilon}\sim 0.05$ , but the limiting lower value of $1+{\it\varepsilon}$ increases weakly with decreasing $q_{edge}$ . Nevertheless, if it can be created and pushed through the unstable band using active stabilization, and if it can subsequently be maintained by active profile control, the resulting high current, potentially high ${\it\beta}$ configuration would result in a highly attractive Advanced Tokamak scenario.

5.2 Significance of $q_{95}$

The remaining question that originally motivated this study, but has not been answered, is why $q_{edge}$ in limiter discharges is replaced so effectively in major tokamaks by $q_{95}$ when the cross-section is diverted. The glib but partial answer provided here is that $q_{95}$ is simply where the resistivity is sufficiently large that the region outside this is effectively a vacuum as far as the kink mode is concerned. This is not entirely satisfactory, since it begs the question why it is so precisely $q_{95}$ in all major experiments. The only answer that can be given at present is that the width of the steep resistivity region is determined largely by processes, such as neutral particle fluxes across the separatrix, that scale with the physical device size. Most present day diverted tokamaks are roughly of the same size, namely of the order of a metre or two in major radius. Thus $q_{95}$ tends to fall at roughly the same distance (within a factor two) inside the discharge. This happens to coincide with the steep rise in the resistivity profile. In addition, while the limit is consistent with all known cases, it has not been carefully and quantitatively tested in all conditions, particularly at low aspect ratio or high elongations.

While this explanation is not very satisfactory, if it turns out to be true, then it has some serious implications for ITER, which will be an additional factor two larger in physical size than the largest existing experiments – a factor four times DIII-D for example. If the distance to $q_{95}$ scales with machine size as expected, but the steep resistivity location scales differently (e.g. is a fixed distance like a few cm), then $q_{95}$ may no longer be the relevant parameter in ITER. The ITER edge temperatures are also expected to be substantially hotter than conventional experiments as a result of the alpha heating and the resistivity should be correspondingly lower. The scaling of the width of the high resistivity region with machine size and the distance of the 95 % flux surface from the edge, need to be determined.