2.1.1. Extended MHD equations

The 2-D CASTOR code is based upon the resistive single-fluid MHD equations (Kerner et al. Reference Kerner, Goedbloed, Huysmans, Poedts and Schwarz1998). We extended these equations step by step including the viscous force density $-{ \boldsymbol{\nabla }} \boldsymbol{\cdot }\boldsymbol{\varPi }$ , the heat flux term $(\varGamma -1)({m_i}/{\rho }) \boldsymbol{\nabla } \boldsymbol{\cdot }\boldsymbol{q}$ , the electron pressure term $({\alpha }/{\rho }) {\boldsymbol{\nabla }} p_e$ and the Hall term $({\alpha }/{\rho }) \boldsymbol{j} \times {\boldsymbol B}$ . The extended MHD equations read (Schnack et al. Reference Schnack, Barnes, Brennan, Hegna, Held, Kim, Kruger, Pankin and Sovinec2006)

(2.1) \begin{align} \frac {\partial \rho }{\partial t} = -\boldsymbol{V} \boldsymbol{\cdot }\boldsymbol{\nabla }\rho - \rho \boldsymbol{\nabla } \boldsymbol{\cdot }{\boldsymbol V}, \end{align}

(2.2) \begin{align} \rho \frac {\partial {\boldsymbol V}}{\partial t} = -\rho ({\boldsymbol V} \boldsymbol{\cdot }{ \boldsymbol{\nabla }}) {\boldsymbol V} -{\boldsymbol{\nabla }} p - {\boldsymbol{\nabla }} \boldsymbol{\cdot }\boldsymbol{\varPi } + {\boldsymbol j} \times {\boldsymbol B}, \end{align}

(2.3) \begin{align} \frac {\partial T}{\partial t} = -{\boldsymbol V} \boldsymbol{\cdot }{\boldsymbol{\nabla }} T -(\varGamma -1) T {\boldsymbol{\nabla }} \boldsymbol{\cdot }{\boldsymbol V} - (\varGamma -1) \frac {m_i}{\rho } \boldsymbol{\nabla } \boldsymbol{\cdot }{\boldsymbol q}, \end{align}

(2.4) \begin{align} {\boldsymbol E} = -{\boldsymbol V} \times {\boldsymbol B} -\frac {\alpha }{\rho } {\boldsymbol{\nabla }} p_e + \frac {\alpha }{\rho } {\boldsymbol j} \times {\boldsymbol B}+\eta {\boldsymbol j} \end{align}

with $\alpha$ = $\alpha _{\mathrm{s}}( {m_i}/{e})$ , ion mass $m_i$ , elementary charge $e$ and scaling parameter $\alpha _{\mathrm{s}}$ .Footnote ⁴ These equations are derived from the two-fluid MHD equations (Braginskii Reference Braginskii1965) making several assumptions and approximations, e.g. assuming quasineutrality $n=n_i=n_e$ ( $n_i$ = ion density, $n_e$ = electron density), neglecting the electron mass $m_e$ ( $m_i \gg m_e$ , mass density $\rho = m_i n + m_e n \approx m_i n$ ), using the pressure $p=p_i+p_e=n(T_i+T_e)=nT$ ( $T_i$ = ion temperature, $T_e$ = electron temperature, $T=T_i+T_e$ ) and making the velocity approximation $\boldsymbol{V}=(({m_i \boldsymbol{V}_i + m_e \boldsymbol{V}_e})/({m_i+m_e})\approx \boldsymbol{V}_i)$ with $\boldsymbol{V}_i$ being the ion velocity. The set of equations is closed with the Maxwell equations

(2.5) \begin{align} { \boldsymbol{\nabla }} \boldsymbol{\cdot }{\boldsymbol E} = 0,\,\,\,\,\, {\boldsymbol {\nabla }} \times {\boldsymbol E} = -\frac {\partial {\boldsymbol B}}{\partial t},\,\,\,\,\, { \boldsymbol{\nabla }} \times {\boldsymbol B} = \mu _0 {\boldsymbol j},\,\,\,\,\,{\boldsymbol{\nabla }} \boldsymbol{\cdot }{\boldsymbol B}=0 \end{align}

( $\boldsymbol E$ = electric field, $\boldsymbol B$ = magnetic field, $\boldsymbol j$ = current density). Here, the displacement current $\mu _0 \epsilon _0 ({\partial {\boldsymbol E}}/{\partial t})$ is neglected.

In CASTOR3D, the ion velocity

(2.6) \begin{align} {\boldsymbol V}_i = \underbrace {\frac {{\boldsymbol E} \times {\boldsymbol B}}{B^2}}_{={\boldsymbol V}_{E\textrm {x}B}} + \underbrace {\frac {\alpha }{\rho } \frac {{\boldsymbol B} \times {\boldsymbol{\nabla }} p_{i}}{B^2}}_{={\boldsymbol V}_i^*} + {\boldsymbol V}^{ext} \end{align}

is implemented with ExB-velocity ${\boldsymbol V}_{Ex B}$ , ion diamagnetic drift velocity ${\boldsymbol V}_i^*$ and an externally driven velocity ${\boldsymbol V}^{ext}$ . Drives of the latter are, for example, neutral beam injection, ion cyclotron range of frequencies waves, etc. (Rice Reference Rice2016).

The viscous stress tensor $\boldsymbol{\varPi}$ splits into $\boldsymbol\varPi =\boldsymbol\varPi _{\|}+\boldsymbol\varPi _{\wedge }+\boldsymbol\varPi _{\perp }$ with parallel $\boldsymbol{\varPi }_{\|}$ , cross (gyroviscous) $\boldsymbol{\varPi }_{\wedge }$ and perpendicular $\boldsymbol{\varPi }_{\perp }$ parts. The parallel viscous force density is approximated by

(2.7) \begin{align} -{\boldsymbol{\nabla }} \boldsymbol{\cdot }\boldsymbol{\varPi }_{\|}\approx \mu _{\|} \Delta _{\|} {\boldsymbol V} \end{align}

with $\mu _{\|}$ being the dynamic viscosity coefficient. We neglect the perpendicular viscous force density $-\boldsymbol{\nabla } \boldsymbol{\cdot }\boldsymbol{\varPi }_{\perp }$ and use

(2.8) \begin{align} - \boldsymbol{\nabla } \boldsymbol{\cdot }\boldsymbol{\varPi }_{\wedge }\approx \rho ({\boldsymbol V}^*_i \boldsymbol{\cdot }{ \boldsymbol{\nabla }}) {\boldsymbol V} \end{align}

for the gyroviscous force density (Chang & Callen Reference Chang and Callen1992a ). This expression is generally known as gyroviscous cancellation. Using (2.7) and (2.8), the momentum equation (2.2) reads

(2.9) \begin{align} \rho \frac {\partial {\boldsymbol V}}{\partial t} =- \rho (({\boldsymbol V} - {\boldsymbol V}_i^*) \boldsymbol{\cdot }{\boldsymbol{\nabla }}) {\boldsymbol V} - {\boldsymbol{\nabla }} p + \mu _{\|} {\boldsymbol{\nabla} }_{\|} {\boldsymbol V} + {\boldsymbol j} \times {\boldsymbol B}. \end{align}

Presently, only the parallel part (Braginskii Reference Braginskii1965)

(2.10) \begin{align} \boldsymbol{q}_{\|} = \boldsymbol{q}_{\|e}+\boldsymbol{q}_{\|i}= -\kappa _{\|e} \boldsymbol{\nabla} _{\|} T_e -\kappa _{\|i}\boldsymbol{\nabla} _{\|}T_i \end{align}

of the heat flux $\boldsymbol q={\boldsymbol q}_{\|}+{\boldsymbol q}_{\wedge }+{\boldsymbol q}_{\perp }$ is included in the energy equation (2.3) with $\kappa _{\|e}$ and $\kappa _{\|i}$ being the electron and ion thermal conductivity coefficients and ${\boldsymbol{\nabla} }_{\|} =({1}/{B^2}){\boldsymbol B} ({\boldsymbol B} \boldsymbol{\cdot }{\boldsymbol{\nabla }} )$ . Besides, $\varGamma =5/3$ is used as adiabatic index, $\eta$ denotes the plasma resistivity and $\mu _0$ is the vacuum permeability.

The extended Ohm’s law (2.4) is simplified by using the drift approximation of the Hall term (Schnack et al. Reference Schnack, Barnes, Brennan, Hegna, Held, Kim, Kruger, Pankin and Sovinec2006)

(2.11) \begin{align} {\boldsymbol{\nabla} }_{\perp } p \approx {\boldsymbol j} \times {\boldsymbol B} \end{align}

which leads to

(2.12) \begin{align} {\boldsymbol E} = -{\boldsymbol V} \times {\boldsymbol B} - \frac {\alpha }{\rho } { \boldsymbol{\nabla} }_{\|}p_e +\frac {\alpha }{\rho }{\boldsymbol{\nabla} }_{\perp }p_i+ \eta {\boldsymbol j}. \end{align}

Single-fluid codes which include two-fluid effects often use the assumption $p_i=p_e=0.5p$ , that is, $T_i=T_e$ . Since measured ion and electron temperatures may be different, we use

(2.13) \begin{align} p_i=\alpha _T p =n\alpha _T T\,\,\,\,\,\textrm {and}\,\,\,\,\,\,p_e=(1-\alpha _T)p=n(1-\alpha _T) T \end{align}

with $0\leqslant \alpha _T \leqslant 1$ . The parameter $\alpha _T$ provides an additional degree of freedom that allows to choose a more realistic ratio of ion and electron pressures with respect to the measured temperatures, so that a more realistic treatment of the ion and electron diamagnetic drift effects and the heat flux is possible. Application of the relations (2.13) to the ion drift velocity (2.6), the heat flux (2.10) and the pressure terms of Ohm’s law (2.12) yields

(2.14) \begin{align} {\boldsymbol V}_i^* = \frac {\alpha \alpha _T }{\rho } \frac {{\boldsymbol B} \times { \boldsymbol{\nabla }} p} {B^2}, \end{align}

(2.15) \begin{align} \boldsymbol{q}_{\|}=-(\kappa _{\|e}(1-\alpha _T)+\kappa _{\|i}\alpha _T)\boldsymbol{\nabla} _{\|}T \end{align}

and

(2.16) \begin{align} {\boldsymbol E} = -{\boldsymbol V} \times {\boldsymbol B} - \frac {\alpha }{\rho }{\boldsymbol{\nabla} }_{\|}p +\frac {\alpha \alpha _T}{\rho }{ \boldsymbol{\nabla }} p + \eta {\boldsymbol j}. \end{align}

2.1.2. Linearised extended MHD equations

We use $\boldsymbol B=\boldsymbol{\nabla }\times \boldsymbol A$ ( $\boldsymbol A$ = magnetic vector potential) to fulfil $\boldsymbol{\nabla } \boldsymbol{\cdot }\boldsymbol B=0$ , and the Weyl gauge (electric scalar potential $\varPhi _E =0$ ) yields $\boldsymbol E = - ( {\partial \boldsymbol A}/{\partial t})$ . The time dependence of the perturbed quantities is assumed to be $\sim \text{e}^{\lambda t}$ with eigenvalue

(2.17) \begin{align} \lambda =\gamma + \text{i}\omega . \end{align}

Here, $\gamma$ is the growth rate and $\omega$ the oscillation frequency of the mode. Using the ansatz $f({\boldsymbol r},t)=f_0({\boldsymbol r})+ \text{e}^{\lambda t}f_1({\boldsymbol r})$ for $\rho ,\boldsymbol{V},T$ and $\boldsymbol A$ , the linearised forms of the extended fluid equations (2.1, 2.9, 2.3, 2.16) read

(2.18) \begin{align} \lambda \rho _1 = -{\boldsymbol V}_0 \boldsymbol{\cdot }{\boldsymbol{\nabla }} \rho _1- {\boldsymbol V}_1 \boldsymbol{\cdot }{ \boldsymbol{\nabla }} \rho _0 -\rho _1 {\boldsymbol{\nabla }} \boldsymbol{\cdot }{\boldsymbol V}_0 - \rho _0 {\boldsymbol{\nabla }} \boldsymbol{\cdot }{\boldsymbol V}_1, \end{align}

(2.19) \begin{align} \rho _0 \lambda {\boldsymbol V}_1 &= - \rho _1({\boldsymbol V}_0 \boldsymbol{\cdot }{\boldsymbol{\nabla }}) {\boldsymbol V}_0 - \rho _0({\boldsymbol V}_1 \boldsymbol{\cdot }{ \boldsymbol{\nabla }}) {\boldsymbol V}_0 - \rho _0({\boldsymbol V}_0 \boldsymbol{\cdot }{ \boldsymbol{\nabla }}) {\boldsymbol V}_1 \nonumber \\[3.7pt] &\quad + \rho _1({\boldsymbol V}_{i0}^* \boldsymbol{\cdot }{ \boldsymbol{\nabla }}) {\boldsymbol V}_0 + \rho _0({\boldsymbol V}_{i1}^* \boldsymbol{\cdot }{\boldsymbol{\nabla }}) {\boldsymbol V}_0 + \rho _0({\boldsymbol V}_{i0}^* \boldsymbol{\cdot }{\boldsymbol{\nabla }}) {\boldsymbol V}_1 \nonumber \\[3.7pt] &\quad - \frac {1}{m_i}{ \boldsymbol{\nabla }}(\rho _0 T_1) - \frac {1}{m_i}{ \boldsymbol{\nabla }}(T_0 \rho _1) +\mu _{||} \Delta _{||} {\boldsymbol V}_1 \nonumber \\[3.7pt] &\quad + \frac {1}{\mu _0} ({\boldsymbol{\nabla }} \times ({ \boldsymbol{\nabla }} \times {\boldsymbol A}_1) ) \times ({\boldsymbol{\nabla }} \times {\boldsymbol A}_0) \nonumber \\[3.7pt] &\quad + \frac {1}{\mu _0} ({\boldsymbol{\nabla }} \times ({\boldsymbol{\nabla }} \times {\boldsymbol A}_0) ) \times ({\boldsymbol{\nabla }} \times {\boldsymbol A}_1), \end{align}

(2.20) \begin{align} \lambda T_1 &= - {\boldsymbol V}_1 \boldsymbol{\cdot }{\boldsymbol{\nabla }} T_0 - {\boldsymbol V}_0 \boldsymbol{\cdot }{\boldsymbol{\nabla }} T_1 -(\varGamma -1) T_1 {\boldsymbol{\nabla }} \boldsymbol{\cdot }{\boldsymbol V}_0 \nonumber \\[3.5pt] &\quad -(\varGamma -1) T_0 { \boldsymbol{\nabla }} \boldsymbol{\cdot }{\boldsymbol V}_1 -(\varGamma -1)\frac {m_i}{\rho _0}\boldsymbol{\nabla }\boldsymbol{\cdot }{\boldsymbol q}_{\|1}, \end{align}

(2.21) \begin{align} \lambda \boldsymbol A_1 &= {\boldsymbol V}_1 \times {\boldsymbol B}_0 + {\boldsymbol V_0} \times ({\boldsymbol{\nabla }} \times {\boldsymbol A}_1) +\frac {\alpha }{\rho _0}\frac {1}{B^2_0} {\boldsymbol B}_0 \bigg(({\boldsymbol{\nabla }} \times {\boldsymbol A}_1) \boldsymbol{\cdot }\frac {1}{m_i} {\boldsymbol{\nabla }} (\rho _0 T_0)\bigg) \nonumber \\[3.5pt] &\quad +\frac {\alpha }{\rho _0}\frac {1}{B^2_0} {\boldsymbol B}_0 \bigg({\boldsymbol B}_0 \boldsymbol{\cdot }\frac {1}{m_i} {\boldsymbol{\nabla }}(\rho _1 T_0 + \rho _0 T_1)\bigg) -\alpha \alpha _T \frac {1}{\rho _0}\frac {1}{m_i} { \boldsymbol{\nabla }}(\rho _1 T_0 +\rho _0 T_1) \nonumber \\[3.5pt] &\quad + \alpha \alpha _T \frac {\rho _1}{\rho _0^2} \frac {1}{m_i} {\boldsymbol{\nabla }}(\rho _0 T_0) -\frac {\eta }{\mu _0} ({ \boldsymbol{\nabla }} \times ({ \boldsymbol{\nabla }} \times {\boldsymbol A}_1)). \end{align}

Figure 3.

Cylindrical ( $R,\varphi ,Z$ ) and magnetic flux coordinates ( $s,v,u$ ). The flux surfaces described by $s= \text{const.}$ lines are shown for a toroidal cross-section of a tokamak plasma.

Here, $\rho _0,{\boldsymbol V}_0, T_0,{\boldsymbol A}_0$ are the equilibrium quantities of density, ion velocity, temperature ( $T_0=T_{0e}+T_{0i}$ ) and magnetic vector potential, while $\rho _1,{\boldsymbol V}_1, T_1, {\boldsymbol A}_1$ are their perturbations. The perturbation of the ion diamagnetic drift velocity (2.14) and the parallel heat flux term $\boldsymbol{\nabla }\boldsymbol{\cdot }{\boldsymbol q}_{\|}$ read

(2.22) \begin{align} {\boldsymbol V}_{i1}^*&=-\alpha \alpha _T \frac {\rho _1}{\rho _0^2}\frac {{\boldsymbol B}_0 \times {\boldsymbol{\nabla }} p_0}{B_0^2}+\alpha \alpha _T \frac {1}{\rho _0} \frac {{\boldsymbol B}_0 \times (T_0 { \boldsymbol{\nabla }} \rho _1 + T_1 {\boldsymbol{\nabla }} \rho _0 +\rho _0 {\boldsymbol{\nabla }} T_1 +\rho _1 { \boldsymbol{\nabla }} T_0)}{B_0^2} \nonumber \\[5pt] &\quad +\alpha \alpha _T\frac {1}{\rho _0} \frac {({\boldsymbol{\nabla }} \times {\boldsymbol A_1}) \times {\boldsymbol{\nabla }} p_0}{B_0^2} -2\alpha \alpha _T\frac {1}{\rho _0} \frac {{\boldsymbol B}_0 \boldsymbol{\cdot }({ \boldsymbol{\nabla }} \times {\boldsymbol A_1})}{B_0^4}({\boldsymbol B}_0 \times { \boldsymbol{\nabla }} p_0) \end{align}

and

(2.23) \begin{align} \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{q}_{\|1}&=-(\kappa _{\|e}(1-\alpha _T) +\kappa _{\|i}\alpha _T)\boldsymbol{\cdot }\nonumber \\[5pt] &\quad \left [(\boldsymbol{B}_0\boldsymbol{\cdot }\boldsymbol{\nabla }) \frac {((\boldsymbol{\nabla }\times \boldsymbol{A}_1)\boldsymbol{\cdot } \boldsymbol{\nabla })T_0}{B_0^2} + \boldsymbol{\nabla }\boldsymbol{\cdot }\frac {\boldsymbol{B}_0 (\boldsymbol{B}_0\boldsymbol{\cdot }\boldsymbol{\nabla })T_1} {B_0^2}\right ]. \end{align}

2.1.3. Weak formulation of the MHD equations

The accuracy and the numerical effort of the solution of the linearised MHD equations (2.18–2.21) depends strongly on the chosen coordinates. Since the behaviour of the plasma is extremely anisotropic with respect to the parallel and the perpendicular direction to the magnetic field $\boldsymbol{B}_0$ , magnetic flux coordinates are the preferential choice. Figure 3 shows flux surfaces of a tokamak plasma at a toroidal cross-section $\varphi = \text{const.}$ and the directions of right-handed cylindrical ( $R,\varphi ,Z$ ) and magnetic flux coordinates ( $s,v,u$ ). The toroidal $v$ - and poloidal $u$ -coordinates are angle-like variables with $0 \leqslant v \lt 1$ and $0 \leqslant u \lt 1$ . Coordinate surfaces with constant radial coordinate $s$ correspond to flux surfaces with $s=0$ being the magnetic axis and $s=1$ representing the plasma boundary. Using the right-handed magnetic flux coordinates and the so-called weak formulation, an eigenvalue problem is derived from the linearised MHD equations. For this purpose, we use a procedure analogous to the one described in Huysmans et al. (Reference Huysmans, Goedbloed and Kerner1993) and Kerner et al. (Reference Kerner, Goedbloed, Huysmans, Poedts and Schwarz1998).

In a first step, a complex ansatz is made for the perturbed quantities

(2.24) \begin{align} \rho _1 = \frac {1}{s} \hat \rho , \end{align}

(2.25) \begin{align} {\boldsymbol v_1} =\frac {R^2}{\sqrt {g}}{\hat v}^s {\boldsymbol r},_s + \frac {R^2}{\text{i}\sqrt {g}}{\hat v}^u {\boldsymbol r},_u + \frac {R^2 }{\text{i} \iota {\hat f}} {\hat v}^v {\boldsymbol B_0}, \end{align}

(2.26) \begin{align} T_1 = \frac {1}{s} \hat T \end{align}

and

(2.27) \begin{align} {\boldsymbol A_1} = -\text{i} {\hat A}_s \boldsymbol{\nabla }s + {\hat A}_u \boldsymbol{\nabla }u + {\hat A}_v \boldsymbol{\nabla }v \end{align}

with Jacobian $\sqrt g$ , $\hat f$ = $({\partial \varPhi _0}/{\partial s})$ and $\varPhi _0$ being the toroidal flux. The ansatz is inserted into the linearised equations (2.18–2.21) which are summarised in an eigenvalue problem

(2.28) \begin{align} \lambda {{\unicode{x1D647}}} \boldsymbol{w} = {{\unicode{x1D64D}}} \boldsymbol{w}. \end{align}

The symbolic vector

(2.29) \begin{align} \boldsymbol{w} = (\hat \rho , \hat v^s, \hat v^u, \hat v^v, \hat T, \hat A_s, \hat A_u, \hat A_v) \end{align}

comprises the eight variables.

For the numerical calculations, each component of $\boldsymbol w$ is approximated by a Fourier seriesFootnote ⁵ in the $u$ - and $v$ -directions and a finite element interpolation in the $s$ -direction:

(2.30) \begin{align} \hat w_{f}(s,v,u) &= \sum _{{m,n,j,q}}\big ( a^{{\text{s},m,n}}_{{f,j,q}} h^{j}_{{f,q}}(s) \textrm {sin}(2\pi ({\mathit m} u +{\mathit n} v)) + a^{{\text{c},m,n}}_{{f,j,q}} h^{ j}_{{f,q}}(s) \textrm {cos}(2\pi ({m} u + {\mathit n} v))\big) \nonumber \\ &= \sum _{{m,n,j,q}} c^{{m,n}}_{{f,j,q}} h^{j}_{{f,q}}(s) \text{e}^{2\pi \text{i} ({m} u+{n} v)} + \!\sum _{{m,n,j,q}}\check {c}^{{m,n}}_{{f,j,q}} h^{j}_{{f,q}}(s) \text{e}^{-2\pi \text{i}({m} u+{n} v)} =\hat w_{{f},c}+\hat w_{{f},\check {c}} \end{align}

with $f$ labelling $\hat \rho , \hat v^s, \hat v^u, \hat v^v, \hat T, \hat A_s, \hat A_u, \hat A_v$ . The Hermite polynomials $h^{j}_{{f,q}}(s)$ are either quadratic ( ${ f}= \hat \rho , \hat v^u, \hat v^v, \hat T, \hat A_s$ )Footnote ⁶ or cubic ( ${\mathit f}= \hat v^s, \hat A_u, \hat A_v$ ). The complex Fourier coefficients $\{c^{{m,n}}_{{f,j,q}}, \check c^{{m,n}}_{{f,j,q}}\}$ are related with the complex coefficients $\{a^{\text{s},{m,n}}_{{f,j,q}},a^{\text{c},{m,n}}_{{f,j,q}}\}$ via

(2.31) \begin{align} c^{{m,n}}_{{f,j,q}} =\frac {1}{2}a^{{\text{c},m,n}}_{{f,j,q}}+ \frac {1}{2\text{i}} a^{{\text{s},m,n}}_{{f,j,q}} \,\,\,\,\,\textrm {and}\,\,\,\,\, \check c^{{m,n}}_{{f,j,q}} =\frac {1}{2}a^{{\text{c},m,n}}_{{f,j,q}}- \frac {1}{2\text{i}} a^{{\text{s},m,n}}_{{f,j,q}}. \end{align}

The indices m and n denote the poloidal and toroidal Fourier harmonics of the perturbation spectrumFootnote ⁷ and $q=1,2$ labels the two quadratic and the two cubic Hermite polynomials at a radial grid point $j$ (Kerner et al. Reference Kerner, Goedbloed, Huysmans, Poedts and Schwarz1998).

The integral

(2.32) \begin{align} (\boldsymbol{z}, \boldsymbol{w}) = \int \text{d}\tau \boldsymbol{z}^* \boldsymbol{\cdot }\boldsymbol{w} \end{align}

denotes an inner product in the appropriate Hilbert space. Taking the inner product of (2.28) by using a weighting function $\boldsymbol z$ and by integrating over the plasma volume $\tau$ ,

(2.33) \begin{align} \lambda (\boldsymbol{z}, {{\unicode{x1D647}}} \boldsymbol{w}) = (\boldsymbol{z}, {{\unicode{x1D64D}}}\boldsymbol{w}), \end{align}

we get the matrix eigenvalue problem

(2.34) \begin{align} \lambda {{\unicode{x1D63D}}} \boldsymbol{x} = {{\unicode{x1D63C}}} \boldsymbol{x} \end{align}

with the eigenvector

(2.35) \begin{align} \boldsymbol{x}= \big(\{c^{{m,n}}_{{f,j,q}}, \check c^{{m,n}}_{{f,j,q}}\}\big) \end{align}

and the complex eigenvalue $\lambda$ . The components of the weighting vector $\boldsymbol{z}^*$ read

(2.36) \begin{align} \hat z^{*{{\mathrm{s},\tilde{{m}},\tilde{{n}}}}}_{{{\tilde{{f}},\tilde{{j}},\tilde{{q}}}}}(s,v,u) &= \hat t_{\tilde{{f}}} h^{\tilde{{j}}}_{{\tilde{{f}},\tilde{{q}}}}(s)\textrm {sin}(2\pi (\tilde{m} u +\tilde{n} v)) \nonumber \\ &\quad =\hat t_{\tilde{{f}}} h^{\tilde{{j}}}_{{\tilde{{f}},\tilde{{q}}}}(s)\frac {1}{2\text{i}} \big(\text{e}^{2\pi \text{i}(\tilde{{m}} u + \tilde{{n}} v)} -\text{e}^{-2 \pi \text{i}(\tilde{{m}} u + \tilde{{n}} v)}\big) \end{align}

and

(2.37) \begin{align} \hat z^{*{{\mathrm{c},\tilde{{m}},\tilde{{n}}}}}_{{{\tilde{{f}},\tilde{{j}},\tilde{{q}}}}}(s,v,u) &= \hat t_{\tilde{{f}}} h^{\tilde{{j}}}_{{\tilde{{f}},\tilde{{q}}}}(s) \textrm {cos}(2\pi (\tilde{m} u +\tilde{n} v)) \nonumber \\ &=\hat t_{{\tilde{f}}} h^{\tilde{{j}}}_{{\tilde{{f}},\tilde{{q}}}}(s)\frac {1}{2} \left (\text{e}^{2\pi \text{i}(\tilde{{m}} u + \tilde{{n}} v)} +\text{e}^{-2 \pi \text{i}(\tilde{{m}} u + \tilde{{n}} v)}\right ) \end{align}

with

(2.38) \begin{align} \hat t_{\hat \rho }= \frac {1}{s}, \end{align}

(2.39) \begin{align} \hat t_{\hat v_s}= \frac {R^2}{\sqrt {g}} \boldsymbol{r},_s,\,\,\, \hat t_{\hat v_u}= \frac {\text{i} R^2}{\sqrt {g}} \boldsymbol{r},_u,\,\,\, \hat t_{\hat v_v}= \frac {\text{i} R^2}{\sqrt {g}} \boldsymbol{r},_u - \frac {\text{i} J_0}{4 \pi ^2 \iota \hat f} \boldsymbol{r},_v \end{align}

(2.40) \begin{align} \hat t_{\hat T}= \frac {1}{s}, \end{align}

(2.41) \begin{align} \hat t_{\hat A_s} = \text{i} \boldsymbol{\nabla }s,\,\,\,\hat t_{\hat A_u}=\boldsymbol{\nabla }u,\,\,\, \hat t_{\hat A_v} = \boldsymbol{\nabla }v \end{align}

and $J_0$ being the equilibrium poloidal current.

More details of the weak formulation concerning the Hermite polynomials and the computation of the matrix elements are described in the Appendices A.1 and A.2. Eventually, the structures of the matrices $\unicode{x1D63C}$ and $\unicode{x1D63D}$ , their dimensions and possibilities to reduce them are discussed in § 3.5.

2.1.4. Boundary conditions

At the magnetic axis, the boundary conditions

(2.42) \begin{align} v^s_1(s=0)=0,\,\,\,A_{1u}(s=0)=0,\,\,\,A_{1v}(s=0)=0 \end{align}

result from the regularity conditions (Huysmans Reference Huysmans1991).

If the plasma is enclosed by a perfectly conducting wall, the boundary conditions read

(2.43) \begin{align} \boldsymbol{n}_{\mathrm{p}} \boldsymbol{\cdot }\boldsymbol{v}_1=0,\,\,\,\,\, \boldsymbol{n}_{\mathrm{p}} \boldsymbol{\cdot }\boldsymbol{B}_1=0,\,\,\,\,\,\, \boldsymbol{n}_{\mathrm{p}}\times \boldsymbol{E}_1 =0 \end{align}

with $\boldsymbol{n}_{\mathrm{p}}$ being the exterior unit normal vector at the plasma boundary. That is,

(2.44) \begin{align} v^s_1(s=1)=0,\,\,\,A_{1u}(s=1)=0,\,\,\,A_{1v}(s=1)=0. \end{align}

Both boundary conditions (2.42) and (2.44) are essential boundary conditions. They must be applied to the variables and the test functions (Kerner et al. Reference Kerner, Goedbloed, Huysmans, Poedts and Schwarz1998). As in the CASTOR code, these conditions are implemented in the CASTOR3D code by removing the rows and columns of the corresponding matrix elements on the magnetic axis and the plasma boundary.

If the plasma is surrounded by a vacuum region, the corresponding boundary conditions implemented in CASTOR and CASTOR3D are restricted to plasma equilibria without surface currents. That is, the equilibrium pressure has to fall to zero at the plasma boundary, $p_0 (s=1)=0$ , so that the equilibrium magnetic field is continuous across the plasma boundary, $\boldsymbol{B}_0=\boldsymbol{B}_0^{vac}$ . Then, in the case of an ideal plasma, the normal component of the magnetic field perturbation and the total pressure perturbation

(2.45) \begin{align} \boldsymbol{B}_1\boldsymbol{\cdot }\boldsymbol{n}_{\mathrm{p}} = \boldsymbol{B}_1^{vac}\! \boldsymbol{\cdot }\boldsymbol{n}_{\mathrm{p}}\,\,\,\,\,\textrm {and}\,\,\,\,\,p_1+ \boldsymbol{B}_0 \boldsymbol{\cdot }\boldsymbol{B}_1 = \boldsymbol{B}_0^{vac} \boldsymbol{\cdot }\boldsymbol{B}_1^{vac} \end{align}

are continuous at plasma boundary. It is $p_1=\rho _0 T_1 + T_0 \rho _1$ .

In the case of resistive stability studies, surface currents are no longer allowed in the plasma perturbation. This leads to the additional condition

(2.46) \begin{align} \boldsymbol{n}_{\mathrm{p}} \times \boldsymbol{B}_1 = \boldsymbol{n}_{\mathrm{p}} \times \boldsymbol{B}^{vac}_1. \end{align}

Then, the tangential components of the magnetic field perturbations are also continuous,

(2.47) \begin{align} \boldsymbol{B}_0 \boldsymbol{\cdot }\boldsymbol{B}_1 = \boldsymbol{B}_0^{vac} \boldsymbol{\cdot }\boldsymbol{B}^{vac}_1, \end{align}

and it has to be $p_1=0$ at the plasma boundary. In order to fulfil this condition, only equilibria with $p_0(s=1)=0$ and $( {\partial p_0}/{\partial s})\big |_{s=1}=0$ are allowed for free-boundary resistive stability studies because then this condition is trivially satisfied (Huysmans et al. Reference Huysmans, Goedbloed and Kerner1993).

The weak formulation of the linearised MHD equations (2.18)–(2.21) yields the following boundary terms which result from partial integrationFootnote ⁸ of the corresponding terms in the momentum equation (2.19):

(2.48) \begin{align} W_{\hat v_s,\tilde{{q}}}^{{\mathrm{s},\tilde{{m}},\tilde{{n}}}} = - \frac {1}{\mu _0} \int _0^1 \text{d}v\int _0^1 \text{d}u\, H_{\tilde{{q}}} \sin (2\pi (\tilde{{m}}u+\tilde{{n}}v)) R^2 \boldsymbol{B}_0 \boldsymbol{\cdot }\boldsymbol{B}^{vac}_1)|_{s=1}, \end{align}

(2.49) \begin{align} W_{\hat v_s,\tilde{{q}}}^{{\mathrm{c},\tilde{{m}},\tilde{{n}}}} = -\frac {1}{\mu _0} \int _0^1 \text{d}v \int _0^1 \text{d}u\,H_{\tilde{{q}}} \cos (2\pi (\tilde{{m}}u+\tilde{{n}}v)) R^2 \boldsymbol{B}_0 \boldsymbol{\cdot }\boldsymbol{B}^{vac}_1)|_{s=1} \end{align}

and Ohm’s law (2.21),

(2.50) \begin{align} W_{\hat A_u,\tilde{{q}}}^{{\mathrm{s},\tilde{{m}},\tilde{{n}}}} &= \frac {1}{\mu _0} \int _0^1\text{d}v \int _0^1\text{d}u\, H_{\tilde{{q}}} \sin (2\pi (\tilde{{m}}u+\tilde{{n}}v))(-\eta ) B_{1v}^{vac}|_{s=1},\\[-10pt]\nonumber\end{align}

(2.51) \begin{align} W_{\hat A_u,\tilde{{q}}}^{{\mathrm{c},\tilde{{m}},\tilde{{n}}}} & = \frac {1}{\mu _0}\int _0^1\text{d}v \int _0^1\text{d}u\, H_{\tilde{{q}}} \cos (2\pi (\tilde{{m}}u+\tilde{{n}}v))(-\eta ) B_{1v}^{vac}|_{s=1},\\[-10pt]\nonumber\end{align}

(2.52) \begin{align} W_{\hat A_v,\tilde{{q}}}^{{\mathrm{s},\tilde{{m}},\tilde{{n}}}} & = \frac {1}{\mu _0}\int _0^1\text{d}v \int _0^1\text{d}u\, H_{\tilde{{q}}} \sin (2\pi (\tilde{{m}}u+\tilde{{n}}v)) \eta B_{1u}^{vac}|_{s=1},\\[-10pt]\nonumber\end{align}

(2.53) \begin{align} W_{\hat A_v,\tilde{{q}}}^{{\mathrm{c},\tilde{{m}},\tilde{{n}}}} & = \frac {1}{\mu _0}\int _0^1\text{d}v \int _0^1\text{d}u\, H_{\tilde{{q}}} \cos (2\pi (\tilde{{m}}u+\tilde{{n}}v)) \eta B_{1u}^{vac}|_{s=1}. \end{align}

The cubic Hermite polynomials are $H_1(s=1)=1$ and $H_2(s=1)=0$ . That is, only terms with index $\tilde{q}=1$ do not vanish. Furthermore, the boundary terms of Ohm’s law completely disappear in the case of ideal stability studies.

2.2.1. Energy functional and boundary conditions

In order to describe the basic elements, we concentrate on a configuration with one smooth closed wall, a complex multiply connected structure and ${N}_c$ coils. Using the indices ${i}=1$ or $\text{p}$ for the plasma boundary, ${i}=2$ or $\textrm {w}_{\textrm f}$ for the closed wall, ${i}=3$ or $\textrm {w}_\mathrm{t}$ for the multiply connected wall and ${j}$ = 1, … , ${N}_c$ for the coils, we introduce the Lagrangian

(2.54) \begin{align} \mathcal{L} &= \frac {\mu _0}{8\pi } \sum _{{i,i'=1}}^{3,3} \int _{S_{i}} \text{d}\kern-2pt f_{{i}} \int _{S_{{i'}}} {\text{d}\kern-2pt f}_{{i'}} \frac {{\boldsymbol{K}}_{{i}}\boldsymbol{\cdot }{\boldsymbol{K}}_{{i'}}} {|{\boldsymbol{r}}_{{i}} - {\boldsymbol{r}}_{{i'}}|} + \frac {\mu _0}{8\pi } \sum _{{j,j'=1}}^{{N}_c, {N}_c} I_{{j}} I_{{j'}} \int _{S_{{j}}} \text{d}\kern-2pt f_{{j}} \int _{S_{{j'}}} \text{d}\kern-2pt f_{{j'}} \frac {{\boldsymbol{k}}_{{j}} \boldsymbol{\cdot }{\boldsymbol{k}}_{{j'}}}{|{\boldsymbol{r}}_{{j}} - {\boldsymbol{r}}_{{j'}}|} \nonumber \\ &\quad +\frac {\mu _0}{4\pi } \sum _{{j,i}=1}^{{N}_c,3} I_{{j}} \int _{S_{{j}}} \text{d}\kern-2pt f_{{j}} \int _{S_{{i}}} \text{d}\kern-2pt f_{{i}} \frac {{\boldsymbol{K}}_{{i}} \boldsymbol{\cdot }{\boldsymbol{k}}_{{j}}} {|{\boldsymbol{r}}_{{i}} - {\boldsymbol{r}}_{{j}}|} - \int _{S_{\text{p}}} \text{d}\kern-2pt f_{\text{p}} \boldsymbol{K}_{\text{p}} \boldsymbol{\cdot }\boldsymbol{A}_1 \nonumber \\ &\quad + \frac {1}{2\lambda } \sum _{{i}=2}^{3}\int _{S_{{i}}} \text{d}\kern-2pt f_{{i}} \frac {1}{\sigma _{{i}}d_{{i}}} {\boldsymbol{K}}_{{i}} \boldsymbol{\cdot }{\boldsymbol{K}}_{{i}} + \frac {1}{2 \lambda } \sum _{{j}=1}^{{N}_c} R_{{j}} I_{{j}}^2 \int _{S_{{j}}} \text{d}\kern-2pt f_{{j}} {\boldsymbol{k}}_{{j}}^2 ,\end{align}

which describes the magnetic energy of the currents in the plasma boundary and in the external wall and coil structures. The latter are considered as passive conductors. The vector quantities $\boldsymbol{K}_{{i}}$ are the surface current densities on the plasma and wall surfaces, and $\boldsymbol{k}_{{j}}$ are the corresponding normalised quantities on the coil surfaces. Here $I_{{j}}$ denotes the magnitude of the current in coil j and $R_{j}$ its resistances, and $\sigma _{{i}}$ and $d_{{i}}$ are the specific conductivity and thickness of wall structure i.Footnote ¹⁰ The surface integrals $\int _{S_{{i}}}\text{d}\kern-2pt f_{{i}}$ and $\int _{S_{{j}}}\text{d}\kern-2pt f_{{j}}$ are performed over the plasma surface and the wall and coil surfaces. The first term of the Lagrangian (2.54) describes the magnetic energy related to the interaction of the virtual current on the plasma boundaryFootnote ¹¹ with itself and with the wall currents and the interactions of the wall currents with themselves. The second and the third terms are related to the interaction of the coil currents with themselves and with the surface current densities on the plasma boundary and the wall structures. The last three terms are boundary terms. The first boundary term depends on the vector potential perturbation $\boldsymbol{A}_1$ of the plasma while the last two terms depend on the eigenvalue $\lambda$ of the plasma perturbation. The boundary terms result from the following boundary conditions.

Applying the condition that the tangential component of the electric field perturbation is continuous at the plasma boundary,

(2.55) \begin{align} \boldsymbol{n}_{\mathrm{p}} \times \boldsymbol{E}^{vac}_1=\boldsymbol{n}_{\mathrm{p}} \times \boldsymbol{E}_1 \end{align}

together with the relation $\boldsymbol{E}_1$ = $-\lambda \boldsymbol{A}_1$ and performing the transformations

(2.56) \begin{align} \boldsymbol{n}_{\mathrm{p}} \boldsymbol{\cdot }(\boldsymbol{K}_{\mathrm{p}} \times (\boldsymbol{n}_{\mathrm{p}} \times \boldsymbol{A}^{vac}_1)) = \boldsymbol{n}_{\mathrm{p}} \boldsymbol{\cdot }(\boldsymbol{K}_{\mathrm{p}} \times (\boldsymbol{n}_{\mathrm{p}} \times \boldsymbol{A}_1)), \end{align}

we get

(2.57) \begin{align} \boldsymbol{K}_{\mathrm{p}} \boldsymbol{\cdot }\boldsymbol{A}^{vac}_1 = \boldsymbol{K}_{\mathrm{p}} \boldsymbol{\cdot }\boldsymbol{A}_1, \end{align}

that is,

(2.58) \begin{align} \int _{S_{\mathrm{p}}} \text{d}\kern-2pt f_{\mathrm{p}} \boldsymbol{K}_{\mathrm{p}} \boldsymbol{\cdot }\boldsymbol{A}^{vac}_{1} = \int _{S_{\mathrm{p}}} \text{d}\kern-2pt f_{\mathrm{p}} \boldsymbol{K}_{\mathrm{p}} \boldsymbol{\cdot }\boldsymbol{A}_1. \end{align}

Similarly, the magnetic vector potentials of the wall and coil surfaces are expressed by the corresponding relations derived from Ohm’s law $\boldsymbol{E}=\sigma {\boldsymbol j}$ with $\boldsymbol{j}$ being the current density. That is,

(2.59) \begin{align} \boldsymbol{A}^{vac}_1=\begin{cases} \dfrac {1}{\lambda \sigma _{{i}} d_{{i}}} {\boldsymbol{K}_{{i}}} &\mathrm{\,\,\,on\,resistive\,wall\,surfaces} \,S_{{i}}, \\[14pt] \dfrac {R_{{j}}}{\lambda }I_{{j}} {\boldsymbol{k}_{{j}}} &\mathrm{\,\,\,on\,resistive\,coil\,surfaces} \,S_{{j}}. \\ \end{cases} \end{align}

2.2.2. Vacuum equations

The divergence-free surface currents $\boldsymbol{K}_{{i}}$ are derived from current potentials $\varPhi ^i$ ,

(2.60) \begin{align} {\boldsymbol{K}}_{{i}}={\boldsymbol{n}}_{{i}} \times \boldsymbol{\nabla }\varPhi ^{{i}} \end{align}

with

(2.61) \begin{align} \boldsymbol{n}_{{i}} = \frac {\boldsymbol{r},_v^{{i}}\times \boldsymbol{r},_u^{{i}}} {|\boldsymbol{r},_v^{{i}}\times \boldsymbol{r},_u^{{i}}|} \end{align}

being the exterior unit normal vector on surface i.

A single-valued Fourier ansatz plus net-poloidal, $J^{{i}}_{pol}$ , and net-toroidal, $I^{{i}}_{tor}$ , currents is made for the current potential

(2.62) \begin{align} \varPhi ^{{i}} (u,v) &= J^{{i}}_{pol} v + I^{{i}}_{tor}u \nonumber \\ &\quad +\sum _{{\substack{{m}_i \geqslant 0,{n}_i \ne 0\\ {m}_i \gt 0,{n}_i=0}}} \big({\hat \varPhi }^{{i,\text{s}}}_{{m_i,n_i}} \sin (2\pi ({m_i} u+ {n_i} v))+ {\hat \varPhi }^{{i,\text{c}}}_{{m_i,n_i}} \cos (2\pi ({m_i} u+ {n_i} v))\big) \end{align}

on the plasma boundary and on the closed wall with smooth geometry. Here $\varPhi ^{{i}}$ is a many-valued potential. Such potentials are necessary to treat the m = 0, n = 0 harmonic of a vertical instability correctly. An m = 0, n = 0 normal displacement which corresponds to the dilation/contraction of the plasma induces a net-poloidal and a net-toroidal flux in the vacuum region. These fluxes are related with net currents in the plasma boundary and an external conductor (Lüst & Martensen Reference Lüst and Martensen1960). Inserting (2.61) and (2.62) in to (2.60) yields

(2.63) \begin{align} \boldsymbol{K}_{{i}} &= \frac {1}{|{\boldsymbol{r}},_v^{{i}} \times {\boldsymbol{r}},_u^{{i}}|} \left (J^{{i}}_{pol} {\boldsymbol{r}},_u^{{i}}-I^{{i}}_{tor} {\boldsymbol{r}},_v^{{i}}\right ) \nonumber \\ &\quad + \boldsymbol{n}_{{i}} \times \boldsymbol{\nabla }\sum _{{\substack{{m}_i \geqslant 0,{n}_i \ne 0\\ {m}_i \gt 0,{n}_i=0}}} \big({\hat \varPhi }^{{i,\text{s}}}_{{m_i,n_i}} \sin (2\pi ({m_i} u+ {n_i} v))+ {\hat \varPhi }^{{i,\text{c}}}_{{m_i,n_i}} \cos (2\pi ({m_i} u+ {{n_i}} v))\big). \end{align}

Figure 4.

Schematic representations of (a) a wall triangle, (b) a coil segment (reproduced with permission from Nucl. Fusion vol. 57, 2017, p. 016032 and EUROfusion; copyright EURATOM 2016) and (c) a wall segment.

Walls with complex, multiply connected structures are discretised into triangles. The position vector of a triangle t (see figure 4 a) of surface i is given by

(2.64) \begin{align} {\boldsymbol{r}}^{{i,t}}= {\boldsymbol{r}}^{{i,t}}_1 + \alpha {\boldsymbol{r}}^{{i,t}}_{21} + \beta {\boldsymbol{r}}^{{i,t}}_{31},\,\,\, 0\lt \alpha +\beta \lt 1,\,\,\,{\boldsymbol{r}}^{{i,t}}_{{lk}}= {\boldsymbol{r}}^{{i,t}}_{{l}}- {\boldsymbol{r}}^{{i,t}}_{{k}},\,\,\,{\mathit{l,k}}=1,2,3 \end{align}

with the vertices numbered anticlockwise. Assuming that the surface current density is constant on a triangle, it can be written as

(2.65) \begin{align} {\boldsymbol{K}}_{{i,t}}= \frac {\phi ^{{i,t}}_1 {\boldsymbol{r}}^{{i,t}}_{23} + \phi ^{{i,t}}_2 {\boldsymbol{r}}^{{i,t}}_{31} + \phi ^{{i,t}}_3 {\boldsymbol{r}}^{{i,t}}_{12}} {\big|{\boldsymbol{r}}^{{i,t}}_{21} \times {\boldsymbol{r}}^{{i,t}}_{32} \big|} \end{align}

or

(2.66) \begin{align} {\boldsymbol{K}}_{{i,t}}= \sum _{{h=1}}^3 \varPhi ^{{i,t}}_{{h}} {\boldsymbol{e}}^{{t}}_{{i,h}}, \,\,\,\,\,{\boldsymbol{e}}^{{t}}_{{i,h}}=\frac {1}{2} \epsilon _{{hkl}}^{{i,t}} \frac {{\boldsymbol{r}}^{{i,t}}_{{kl}}} {\big|{\boldsymbol{r}}^{{i,t}}_{21} \times {\boldsymbol{r}}^{{i,t}}_{31}\big|}. \end{align}

Then, the total surface current density of wall structure ${i}$ reads

(2.67) \begin{align} {\boldsymbol{K}_{{i}}} = \sum _{{t=1}}^{{N_i}} {\boldsymbol{K}}_{{i,t}}= \sum _{{t=1}}^{{N_i}} \sum _{{h=1}}^3 \varPhi ^{{i,t}}_{{h}} {\boldsymbol{e}}^{{t}}_{{i,h}} \end{align}

with ${N}_{i}$ being its total number of triangles.

Figure 4(c) shows a wall segment composed of triangles. Each vertex is characterised by a global node number V with $V=1,\ldots ,N_{{V_i}}$ and $N_{{V_i}}$ being the total number of vertices of the triangulated domain i. The current potential of a multiply connected wall with boundaries and/or holes has to be constant along the edges. That is, the current-potential values have to be the same along an edge. So, the number of independent global nodes $N_{{G_i}}$ is smaller than the number of vertices $N_{{V_i}}$ . The connectivity tensor $\unicode{x1D643}$ dedicates the independent current potential values $\hat \varPhi ^{{i}}_{{g}}$ at the global nodes to potential values $\varPhi ^{{i,t}}_{{h}}$ at the vertices of the individual triangles:

(2.68) \begin{align} \varPhi ^{{i,t}}_{{h}}= \sum _{{g=1}}^{{N_{G_i}}}{\unicode{x1D643}}({h,t,g}) \hat \varPhi ^{{i}}_{{g}}\,\,\,\,\,{\text{with}}\,\,\,\,\, {h=1,2,3}\,\,\,\,\,\text{and}\,\,\,\,\,{t=1,\ldots ,3N_i}. \end{align}

Inserting (2.68) into (2.67) yields the total surface current density of wall structure i as function of the independent current potential values:

(2.69) \begin{align} {\boldsymbol{K}_{{i}}} = \sum _{{t=1}}^{{N_i}} \sum _{{h=1}}^3 \sum _{{g=1}}^{{N_{G_i}}}{\unicode{x1D643}}({h,t,g}) {\boldsymbol{e}}^{{t}}_{{i,h}} \hat \varPhi ^{{i}}_{{g}}. \end{align}

A coil is represented by an infinitely thin band which is discretised into triangles as shown in figure 4(b). The vertices of the triangles are located on the boundaries of the coil band. All vertices of one side are characterised by the same current potential. Usually, one potential is set to zero while the other corresponds to the amount of the total current $I_{{j}}$ flowing through coil j. Here, we define $\varPhi _{\mathrm{A}}=0$ and $\varPhi _{\mathrm{B}}=I_{{j}}$ .

Using again the definition of the position vector (2.64) and assuming that the surface current density is constant on a triangle, $\boldsymbol{K}_{{j,t}}$ is given by

(2.70) \begin{align} \boldsymbol{K}_{{j,t=1}} = \varPhi _{\mathrm{B}}^{{j}} \frac {\boldsymbol{r}^{{j}}_{23}} {\big|\boldsymbol{r}^{{j}}_{12} \times \boldsymbol{r}^{{j}}_{32}\big|} = I_{{j}} \boldsymbol{k}_{{j,t=1}}, \end{align}

(2.71) \begin{align} \boldsymbol{K}_{{j,t=2}} = \varPhi _{\mathrm{B}}^{{j}} \frac {\boldsymbol{r}^{{j}}_{14}} {\big|\boldsymbol{r}^{{j}}_{43} \times \boldsymbol{r}^{{j}}_{31}\big|} = I_{{j}} \boldsymbol{k}_{{j,t=2}}, \end{align}

etc. with $\boldsymbol{k}_{{j,t}}$ being the normalised surface current density of triangle ${t}$ . Then, the total surface current density of coil ${j}$ reads

(2.72) \begin{align} {\boldsymbol{K}_{{j}}} = I_{{j}} \sum _{{t=1}}^{{N_j}}{\boldsymbol{k}}_{{j,t}} \end{align}

with $N_{j}$ being its total number of triangles.

Inserting (2.63), (2.69) and (2.72) into the Lagrangian (2.54) and making the variations with respect to $J^{\mathrm{p}}_{pol},I^{\mathrm{p}}_{tor}, \{\hat \varPhi ^{\mathrm{p,s}}_{{m_\mathrm{p},n_\mathrm{p}}}\}, \{\hat \varPhi ^{\mathrm{p,c}}_{{m_\mathrm{p},n_\mathrm{p}}}\}, J^{\mathrm{w_f}}_{pol},I^{\mathrm{w_f}}_{tor}, \{\hat \varPhi ^{\mathrm{w_f},s}_{{m_\mathrm{w_f}},{n_\mathrm{w_f}}}\}, \{\hat \varPhi ^{\mathrm{w_f,c}}_{{m_\mathrm{w_f}},{n_\mathrm{w_f}}}\}, \{\hat \varPhi ^{\mathrm{w_t}}_{{g}}\}$ and the coil currents { $I_{{j}}$ }, we obtain the following set of equations:

(2.73) \begin{align} {\unicode{x1D648}}_{\mathrm{pp}}\boldsymbol{x}_{\mathrm{p}} + {\unicode{x1D648}}_{\mathrm{pc}}\boldsymbol{x}_{\mathrm{c}} + {\unicode{x1D648}}_{\mathrm{pw}}\boldsymbol{x}_{\mathrm{w}} = {\unicode{x1D64D}}_{\mathrm{pd}}\boldsymbol{x}_{\mathrm{d}}, \nonumber \\ {\unicode{x1D648}}_{\mathrm{cp}} \boldsymbol{x}_{\mathrm{p}} + {\unicode{x1D648}}_{\mathrm{cc}} \boldsymbol{x}_{\mathrm{c}} + {\unicode{x1D648}}_{\mathrm{wc}} \boldsymbol{x}_{\mathrm{w}} = -\frac {1}{\lambda } {\unicode{x1D64D}}_{\mathrm{cc}} \boldsymbol{x}_{\mathrm{c}}, \nonumber \\ {\unicode{x1D648}}_{\mathrm{wp}} \boldsymbol{x}_{\mathrm{p}} + {\unicode{x1D648}}_{\mathrm{wc}} \boldsymbol{x}_{\mathrm{c}} + {\unicode{x1D648}}_{\mathrm{ww}} \boldsymbol{x}_{\mathrm{w}} = -\frac {1}{\lambda } {\unicode{x1D64D}}_{\mathrm{ww}} \boldsymbol{x}_{\mathrm{w}}. \end{align}

The symbolic vectors are

(2.74) \begin{align} \boldsymbol{x}_{\mathrm{p}} = \big(\hat \varPhi ^{\mathrm{p,s}}_1,\ldots , \hat \varPhi ^{\mathrm{p,s}}_{{N_\mathrm{p}}},I^{\mathrm{p}}_{pol}, \hat \varPhi ^{\mathrm{p,c}}_1,\ldots ,\hat \varPhi ^{\mathrm{p,c}}_{{N_\mathrm{p}}}, J^{\mathrm{p}}_{tor}\big) \end{align}

with $N_\mathrm{p}$ (m, n)-Fourier harmonics of the plasma boundary,

(2.75) \begin{align} \boldsymbol{x}_{\mathrm{w}} = \left(\hat \varPhi ^{\mathrm{w_f,s}}_1,\ldots , \hat \varPhi ^{\mathrm{w_f,s}}_{{N_\mathrm{w_f}}}, I^{\mathrm{w_f}}_{pol},\hat \varPhi ^{\mathrm{w_f,c}}_1,\ldots , \hat \varPhi ^{\mathrm{w_f,c}}_{{N_\mathrm{w_f}}},J^{\mathrm{w_f}}_{tor}, \hat \varPhi ^{\mathrm{w_t}}_1,\ldots ,\hat \varPhi ^{\mathrm{w_t}}_{{N_\mathrm{G_{w_t}}}}\right) \end{align}

with ${N}_{\mathrm{w_f}}$ (m, n)-Fourier harmonics of the closed wall and $N_{\mathrm{G_{w_t}}}$ global nodes of the discretised wall structure,

(2.76) \begin{align} \boldsymbol{x}_{\mathrm{c}} = (I_1,\ldots .,I_{{N_\mathrm{c}}}) \end{align}

with $N_\mathrm{c}$ coil currents $I_j$ and $\boldsymbol{x}_{\mathrm{d}}$ containing the corresponding boundary elements of the eigenvector (2.35). The boundary term ${\unicode{x1D64D}}_{\mathrm{pd}}\boldsymbol{x}_{\mathrm{d}}$ represents the connection between vacuum and plasma while ${\unicode{x1D64D}}_{\mathrm{cc}}\boldsymbol{x}_{\mathrm{c}}$ and ${\unicode{x1D64D}}_{\mathrm{ww}}\boldsymbol{x}_{\mathrm{w}}$ are the coil and wall boundary terms. The matrices ${\unicode{x1D648}}_{\mathrm{pp}}$ , ${\unicode{x1D648}}_{\mathrm{pw}}$ , … denote the plasma–plasma ${\unicode{x1D648}}_{\mathrm{pp}}$ , the plasma–wall ${\unicode{x1D648}}_{\mathrm{pw}},{\unicode{x1D648}}_{\mathrm{wp}}$ , the plasma–coil ${\unicode{x1D648}}_{\mathrm{pc}},{\unicode{x1D648}}_{\mathrm{cp}}$ , the wall–coil ${\unicode{x1D648}}_{\mathrm{wc}},{\unicode{x1D648}}_{\mathrm{cw}}$ and the coil–coil ${\unicode{x1D648}}_{\mathrm{cc}}$ terms. For details concerning the computation of the matrix elements see (Merkel & Strumberger Reference Merkel and Strumberger2015).

2.3.1. Fixed-boundary eigenvalue problem

Assuming that the plasma boundary is a perfectly ideal conducting wall, the eigenvalue problem (2.34) reads

(2.77) \begin{align} \lambda \left (\begin{array}{c@{\quad}c} {{\unicode{x1D63D}}}_{\mathrm{ll}} & {{\unicode{x1D63D}}}_{\mathrm{ld}} \\ {{\unicode{x1D63D}}}_{\mathrm{dl}} & {{\unicode{x1D63D}}}_{\mathrm{dd}} \\ \end{array}\right ) \left (\begin{array}{c} {\boldsymbol x}_{\mathrm{l}}\\ {\boldsymbol x}_{\mathrm{d}} \end{array}\right ) = \left (\begin{array}{c@{\quad}c} {{\unicode{x1D63C}}}_{\mathrm{ll}} & {{\unicode{x1D63C}}}_{\mathrm{ld}} \\ {{\unicode{x1D63C}}}_{\mathrm{dl}} & {{\unicode{x1D63C}}}_{\mathrm{dd}} \\ \end{array}\right ) \left (\begin{array}{c} {\boldsymbol{x}}_{\mathrm{l}}\\ {\boldsymbol{x}}_{\mathrm{d}} \end{array}\right ) \end{align}

with the eigenvector (2.35)

(2.78) \begin{align} \boldsymbol{x}=(\boldsymbol{x}_{\mathrm{l}}, \boldsymbol{x}_{\mathrm{d}}). \end{align}

The index l characterises matrix and eigenvector elements of the plasma interior while d marks the boundary components. In accordance with the boundary conditions (2.44), the corresponding matrix elements of submatrix ${{\unicode{x1D63C}}}_{\mathrm{dd}}$ are set to zero.

2.3.2. Free-boundary eigenvalue problem

If the plasma is surrounded by a vacuum region with perfectly conducting wall and coil structures, the boundary terms of the Lagrangian (2.54) that are related to the dissipation of the wall and coil resistivities disappear and with them the $\lambda$ -dependence of the vacuum equations (2.73). They are reduced to

(2.79) \begin{align} \left (\begin{array}{c@{\quad}c@{\quad}c} {\unicode{x1D648}}_{\mathrm{pp}} & {\unicode{x1D648}}_{\mathrm{pw}} & {\unicode{x1D648}}_{\mathrm{pc}} \\[2pt] {\unicode{x1D648}}_{\mathrm{wp}} & {\unicode{x1D648}}_{\mathrm{ww}} & {\unicode{x1D648}}_{\mathrm{wc}} \\[2pt] {\unicode{x1D648}}_{\mathrm{cp}} & {\unicode{x1D648}}_{\mathrm{cw}} & {\unicode{x1D648}}_{\mathrm{dc}} \\ \end{array}\right ) \left (\begin{array}{c} \boldsymbol{x}_{\mathrm{p}} \\ \boldsymbol{x}_{\mathrm{w}} \\ \boldsymbol{x}_{\mathrm{c}} \end{array}\right ) = \left (\begin{array}{c} {\unicode{x1D64D}}_{\mathrm{pd}} \boldsymbol{x}_{\mathrm{d}}\\[2pt] 0 \\ 0 \\ \end{array}\right ). \end{align}

Using the boundary condition (2.58) and inserting the definitions of the current surface density $\boldsymbol{K}_{\mathrm{p}}$ (2.60) and the normal vector $\boldsymbol{n}_{\mathrm{p}}$ (2.61) into the first boundary term of the Lagrangian (2.54), making a partial integration and using $\boldsymbol{B}_1= \boldsymbol{\nabla }\times \boldsymbol{A}_1$ , leads to

(2.80) \begin{align} \int _{\mathrm{S_p}} \text{d}\kern-2pt f_{\mathrm{p}} \boldsymbol{K}_{\mathrm{p}} \boldsymbol{\cdot }\boldsymbol{A}^{vac}_1 &= \int _{\mathrm{S_p}} \text{d}\kern-2pt f_{\mathrm{p}} \left (\frac {\sqrt {g}\boldsymbol{\nabla }s} {|\boldsymbol{r},_{{v}} \times \boldsymbol{r},_{{u}}|} \times \boldsymbol{\nabla }\varPhi ^{\mathrm{p}}\right) \boldsymbol{\cdot }\boldsymbol{A}_1 \nonumber \\[2pt] &=\int _{\mathrm{S_p}} \text{d}u \text{d}v \left ( -\frac {\partial \varPhi ^{\mathrm{p}}} {\partial u}A_{1u} + \frac {\partial \varPhi ^{\mathrm{p}}}{\partial v}A_{1v}\right )\nonumber \\[2pt] &= -\int _{\mathrm{S_p}}\text{d}u \text{d}v \varPhi ^{\mathrm{p}} \left (\frac {\partial A_{1u}} {\partial v} - \frac {\partial A_{1v}}{\partial u} \right ) \nonumber \\[2pt] &=-\int _{\mathrm{S_p}} \text{d}u \text{d}v \varPhi ^{\mathrm{p}} \sqrt {g} \boldsymbol{\nabla }s \boldsymbol{\cdot }(\boldsymbol{\nabla }\times \boldsymbol{A}_1) =-\int _{\mathrm{S_p}} \text{d}u \text{d}v \varPhi ^{\mathrm{p}} \sqrt {g}B^s_1. \end{align}

Inserting $\varPhi ^{\mathrm{p}}$ (2.62) and the ansatz

(2.81) \begin{align} \sqrt {g} B^s_1 = \sum _{{m,n}} d^{\mathrm{s}}_{{m,n}}{\sin }(2\pi ({m} u +{n} v)) + d^{\mathrm{c}}_{{m,n}}{\cos }(2\pi ({{m}} u + {{n}} v)) \end{align}

into (2.80) and making the variation with respect to $\varPhi ^{\mathrm{p}}$ , we obtain

(2.82) \begin{align} -\int _{S_{\mathrm{p}}} \text{d}u \text{d}v \delta \varPhi ^{\mathrm{p}} \sqrt {g} B^s_1 = \boldsymbol{R}_{\mathrm{pd}} \boldsymbol{d} \end{align}

with

(2.83) \begin{align} \boldsymbol{d}= \big(\big\{d^{\mathrm{s}}_{{m,n}}\big\}, \big\{d^{\mathrm{c}}_{{m,n}}\big\}\big), \end{align}

and, eventually, we substitute $\boldsymbol{R}_{\mathrm{pd}}\boldsymbol{x_{\mathrm{d}}}$ by $\boldsymbol{R}_{\mathrm{pd}}\boldsymbol{d}$ in (2.79). Because of the linearity of equation system (2.79), no detailed information of the plasma magnetic field perturbation is necessary for its solution. It is sufficient to replace the vector $\boldsymbol{d}$ by a set of $2N_{{mn}}$ Footnote ¹² vectors describing unit perturbations and solving (2.79) individually for each vector. All components of such a unit vector are zero except one, e.g.

(2.84) \begin{align} \boldsymbol{u}=\big(0,\ldots ,0, u^{\mathrm{s}}_{{m,n}}=1,0,\ldots ,0\big). \end{align}

The solutions of the equation system yield the vacuum magnetic field perturbation

(2.85) \begin{align} \boldsymbol{B}^{vac}_1 &= \sum _{{i}=1}^3 \boldsymbol{B}_{1,{i}}^{vac}+ \sum _{{j}=1}^{{N_\mathrm{c}}}\boldsymbol{B}_{1,{j}}^{vac} \nonumber \\ & =\frac {\mu _0}{4\pi } \sum _{{i}=1}^2 \int _{S_{{i}}} \text{d}\kern-2pt f_{{i}}\boldsymbol{K}_{{i}} \times \frac {(\boldsymbol{r}-\boldsymbol{r}_{{i}})} {|\boldsymbol{r}-\boldsymbol{r}_{{i}}|^3}- \frac {\mu _0}{4\pi }\sum _{{t}=1}^{{N_\mathrm{w_t}}} \boldsymbol{K}_{{t}} \times \boldsymbol{\nabla }\int _{S_t} \text{d}\kern-2pt f_{{t}} \frac {1}{|\boldsymbol{r}-\boldsymbol{r}_{{t}}|} \nonumber \\ &\quad -\frac {\mu _0}{4\pi } \sum _{{j}}^{{N_\mathrm{c}}} I_{{j}} \sum _{{t}}^{{N_j}}\boldsymbol{k}_{{j,t}} \times \boldsymbol{\nabla }\int _{S_{{j,t}}}\text{d}\kern-2pt f_{{j,t}} \frac {1}{|\boldsymbol{r}-\boldsymbol{r}_{{j,t}}|} \end{align}

and, with it, the boundary elements (2.48)–(2.53) of matrix $\unicode{x1D63C}$ (2.34) as functions of the unit perturbations. The set of unit perturbations matches the set of poloidal and toroidal Fourier harmonics taken into account for the plasma perturbations (2.30). So, the vacuum response to the unit perturbations contains all information of the vacuum needed for the free-boundary solution of the eigenvalue problem (2.34).

2.3.3. Extended free-boundary eigenvalue problem

In the case of resistive conducting structures in the vacuum region, the equation system (2.73)

(2.86) \begin{align} \lambda \left (\begin{array}{c@{\quad}c@{\quad}c@{\quad}c} -{\unicode{x1D64D}}_{\mathrm{pd}} & {\unicode{x1D648}}_{\mathrm{pp}} & {\unicode{x1D648}}_{\mathrm{pw}} & {\unicode{x1D648}}_{\mathrm{pc}} \\ \boldsymbol{0}_{\mathrm{wd}} & {\unicode{x1D648}}_{\mathrm{wp}} & {\unicode{x1D648}}_{\mathrm{ww}} & {\unicode{x1D648}}_{\mathrm{wc}} \\ \boldsymbol{0}_{\mathrm{cd}} & {\unicode{x1D648}}_{\mathrm{cp}} & {\unicode{x1D648}}_{\mathrm{cw}} & {\unicode{x1D648}}_{\mathrm{cc}} \end{array}\right ) \boldsymbol{\cdot }\left (\begin{array}{c} \boldsymbol{x}_{\mathrm{d}} \\ \boldsymbol{x}_{\mathrm{p}} \\ \boldsymbol{x}_{\mathrm{w}} \\ \boldsymbol{x}_{\mathrm{c}} \end{array}\right ) = -\left (\begin{array}{c@{\quad}c@{\quad}c} \boldsymbol{0}_{\mathrm{pp}} & \boldsymbol{0}_{\mathrm{pw}} & \boldsymbol{0}_{\mathrm{pc}} \\ \boldsymbol{0}_{\mathrm{wp}} & {\unicode{x1D64D}}_{\mathrm{ww}} & \boldsymbol{0}_{\mathrm{wc}} \\ \boldsymbol{0}_{\mathrm{cp}} & \boldsymbol{0}_{\mathrm{cw}} & {\unicode{x1D64D}}_{\mathrm{cc}} \end{array}\right ) \boldsymbol{\cdot }\left (\begin{array}{c} \boldsymbol{x}_{\mathrm{p}} \\ \boldsymbol{x}_{\mathrm{w}} \\ \boldsymbol{x}_{\mathrm{c}} \end{array}\right ) \end{align}

has to be closed with the weak form of the linearised MHD equations (2.34)

(2.87) \begin{align} \lambda \left (\begin{array}{c@{\quad}c} {\unicode{x1D63D}}_{\mathrm{ll}} & {\unicode{x1D63D}}_{\mathrm{ld}} \\ {\unicode{x1D63D}}_{\mathrm{dl}} & {\unicode{x1D63D}}_{\mathrm{dd}} \\ \end{array}\right ) \left (\begin{array}{c} \boldsymbol{x}_{\mathrm{l}} \\ \boldsymbol{x}_{\mathrm{d}} \\ \end{array}\right ) = \left (\begin{array}{c@{\quad}c@{\quad}c@{\quad}c@{\quad}c} {\unicode{x1D63C}}_{\mathrm{ll}} & {\unicode{x1D63C}}_{\mathrm{ld}} & \boldsymbol{0}_{\mathrm{lp}} & \boldsymbol{0}_{\mathrm{lw}} & \boldsymbol{0}_{\mathrm{lc}} \\ {\unicode{x1D63C}}_{\mathrm{dl}} & {\unicode{x1D63C}}_{\mathrm{dd}} & {\unicode{x1D64D}}_{\mathrm{dp}} & {\unicode{x1D64D}}_{\mathrm{dw}} & {\unicode{x1D64D}}_{\mathrm{dc}} \\ \end{array}\right ) \left (\begin{array}{c} \boldsymbol{x}_{\mathrm{l}} \\ \boldsymbol{x}_{\mathrm{d}} \\ \boldsymbol{x}_{\mathrm{p}} \\ \boldsymbol{x}_{\mathrm{w}} \\ \boldsymbol{x}_{\mathrm{c}} \\ \end{array}\right ). \end{align}

Again, the indices l and d name plasma interior and boundary related elements, respectively. The matrices ${\unicode{x1D64D}}_{\mathrm{dp}}$ , ${\unicode{x1D64D}}_{\mathrm{dw}}$ and ${\unicode{x1D64D}}_{\mathrm{dc}}$ are obtained by inserting the definition of the magnetic field perturbation (2.85) into the boundary terms (2.48)–(2.53). Combining (2.86) and (2.87) yields the extended eigenvalue problem,

(2.88) \begin{align}& \lambda \left (\begin{array}{c@{\quad}c@{\quad}c@{\quad}c@{\quad}c} {\unicode{x1D63D}}_{\mathrm{ll}} & {\unicode{x1D63D}}_{\mathrm{ld}} & \boldsymbol{0}_{\mathrm{lp}} & \boldsymbol{0}_{\mathrm{lw}} & \boldsymbol{0}_{\mathrm{lc}} \\ {\unicode{x1D63D}}_{\mathrm{dl}} & {\unicode{x1D63D}}_{\mathrm{dd}} & \boldsymbol{0}_{\mathrm{dp}} & \boldsymbol{0}_{\mathrm{dw}} & \boldsymbol{0}_{\mathrm{dc}} \\ \boldsymbol{0}_{\mathrm{pl}} & -{\unicode{x1D64D}}_{\mathrm{pd}} & {\unicode{x1D648}}_{\mathrm{pp}} & {\unicode{x1D648}}_{\mathrm{pw}} & {\unicode{x1D648}}_{\mathrm{pc}} \\ \boldsymbol{0}_{\mathrm{wl}} & \boldsymbol{0}_{\mathrm{wd}} & {\unicode{x1D648}}_{\mathrm{wp}} & {\unicode{x1D648}}_{\mathrm{ww}} & {\unicode{x1D648}}_{\mathrm{wc}} \\ \boldsymbol{0}_{\mathrm{cl}} & \boldsymbol{0}_{\mathrm{cd}} & {\unicode{x1D648}}_{\mathrm{cp}} & {\unicode{x1D648}}_{\mathrm{cw}} & {\unicode{x1D648}}_{\mathrm{cc}} \\ \end{array}\right ) \left (\begin{array}{c} \boldsymbol{x}_{\mathrm{l}} \\ \boldsymbol{x}_{\mathrm{d}} \\ \boldsymbol{x}_{\mathrm{p}} \\ \boldsymbol{x}_{\mathrm{w}} \\ \boldsymbol{x}_{\mathrm{c}} \\ \end{array}\right ) \nonumber \\[5pt] &\quad=\left (\begin{array}{c@{\quad}c@{\quad}c@{\quad}c@{\quad}c} {\unicode{x1D63C}}_{\mathrm{ll}} & {\unicode{x1D63C}}_{\mathrm{ld}} & \boldsymbol{0}_{\mathrm{lp}} & \boldsymbol{0}_{\mathrm{lw}} & \boldsymbol{0}_{\mathrm{lc}} \\ {\unicode{x1D63C}}_{\mathrm{dl}} & {\unicode{x1D63C}}_{\mathrm{dd}} & {\unicode{x1D64D}}_{\mathrm{dp}} & {\unicode{x1D64D}}_{\mathrm{dw}} & {\unicode{x1D64D}}_{\mathrm{dc}} \\ \boldsymbol{0}_{\mathrm{pl}} & \boldsymbol{0}_{\mathrm{pd}} & \boldsymbol{0}_{\mathrm{pp}} & \boldsymbol{0}_{\mathrm{pw}} & \boldsymbol{0}_{\mathrm{pc}} \\ \boldsymbol{0}_{\mathrm{wl}} & \boldsymbol{0}_{\mathrm{wd}} & \boldsymbol{0}_{\mathrm{wp}} & -{\unicode{x1D64D}}_{\mathrm{ww}} & \boldsymbol{0}_{\mathrm{wc}} \\ \boldsymbol{0}_{\mathrm{cl}} & \boldsymbol{0}_{\mathrm{cd}} & \boldsymbol{0}_{\mathrm{cp}} & \boldsymbol{0}_{\mathrm{cw}} & -{\unicode{x1D64D}}_{\mathrm{cc}} \end{array}\right ) \left (\begin{array}{c} \boldsymbol{x}_{\mathrm{l}} \\ \boldsymbol{x}_{\mathrm{d}} \\ \boldsymbol{x}_{\mathrm{p}} \\ \boldsymbol{x}_{\mathrm{w}} \\ \boldsymbol{x}_{\mathrm{c}} \end{array}\right ) \end{align}

or in abbreviated form

(2.89) \begin{align} \lambda {\unicode{x1D63D}}^{\,ext}\boldsymbol{x}^{ext}= {\unicode{x1D63C}}^{ext}\boldsymbol{x}^{ext}. \end{align}

The capability to solve this extended eigenvalue problem makes the CASTOR3D code a very flexible and powerful numerical tool for linear stability studies of tokamak and stellarator equilibria. A concise description of CASTOR3D will be given in the following section.

3.2.1. Density and temperature profiles

Providing a density profile $n_0(s)$ , the temperature profile is obtained from $T_0=p_0/n_0$ . As described in § 2.1.1, we introduced the parameter $\alpha _T$ which allows a more realistic ratio of ion $T_i=\alpha _T T$ and electron temperature $T_e=(1-\alpha _T)T$ .

3.2.2. Plasma resistivity

Besides user defined resistivity profiles $\eta _0(s)$ , it is also possible to choose between Spitzer (Wesson Reference Wesson1997), Sauter (Sauter et al. Reference Sauter, Angioni and Lin-Liu1999, Reference Sauter, Angioni and Lin-Liu2002) and Redl (Redl et al. Reference Redl, Angioni, Belli and Sauter2021) resistivity in CASTOR3D. The neoclassical resistivity derived by Sauter et al. takes effects due to trapped particles as well as impurities into account. An improved expression for this neoclassical resistivity is given by the resistivity of Redl, which is much more accurate in the high collisionality regime.

3.2.3. Heat conductivity coefficients

We follow Chen, Xu & Lei (Reference Chen, Xu and Lei2019) in order to determine expressions for the heat conductivity coefficients of the parallel heat flux (2.15). The well-known classical Spitzer–Härm (superscript SH) expressions (Spitzer & Härm Reference Spitzer and Härm1953; Braginskii Reference Braginskii1965) for the parallel thermal conductivity coefficients are rigorous in the short-mean-free-path (collisional) regime and assume that the lowest-order particle distribution function is a Maxwellian,

(3.1) \begin{align} \kappa _{\|e}^{SH}=3.2 n_e v^2_{th,e}\tau _e,\,\,\,\,\,\, \kappa _{\|i}^{SH}=3.9 n_i v^2_{th,i}\tau _i \end{align}

with thermal velocities $v_{th,{j}}=(T_{{j}}/m_{{j}})^{1/2}$ and collision times $\tau _{{j}}$ ( ${j}=e,i$ ). This approximation breaks down and can greatly overestimate the parallel heat flux when the mean-free-path is long compared with the temperature gradient scale length. In a collisionless tokamak plasma, typically, the upper limit to the parallel heat flux is of the order of the free streaming (superscriptFS) flux

(3.2) \begin{align} \boldsymbol{q}_{\|{j}}^{\,FS} \simeq -\kappa _{\|{j}}^{FS}\boldsymbol{\nabla} _{\|} T_{{j}}\,\,\,\,\,\text{with}\,\,\,\,\, \kappa _{\|{j}}^{FS} \simeq n_{{j}} v_{th,{j}}q R_0. \end{align}

Here, $q$ is the local safety factor and $R_0$ the major radius.

In CASTOR3D, the user can choose between Spitzer–Härm expressions, free streaming expressions and combinations of both for the ion and electron thermal conductivitiesFootnote ¹⁵ (Chen et al. Reference Chen, Xu and Lei2019):

(3.3) \begin{align} \kappa _{\|{j}} =\,{\min}\big(\kappa _{\|{j}}^{FS},\kappa _{\|{j}}^{SH}\big) \,\,\,\,\,\textrm {or}\,\,\,\,\, \kappa _{\|{j}}=\left (\frac {1}{\kappa _{\|{j}}^{FS}}+ \frac {1}{\kappa _{\|{j}}^{SH}}\right )^{-1}. \end{align}

3.2.4. Parallel dynamic viscosity coefficient

The parallel dynamic viscosity coefficient $\mu _{\|}$ (2.7) has to be provided by the user. Rough estimates of $\mu _{\|}$ are

(3.4) \begin{align} \mu _{\|}=\begin{cases} 0.96 \dfrac {v_{th,i}^2}{\nu _i}\rho _i &\mathrm{for}\ \nu _{i}\gt k_{\|}v_{th,i} \\[9pt] \kappa _{\|}\dfrac {v_{th,i}}{k_{\|}}\rho _i &\mathrm{for\,}\nu _{i}\lt k_{\|}v_{th,i} \\ \end{cases} \end{align}

for a collisional (Braginskii Reference Braginskii1965) and a collisionless plasma (Waltz Reference Waltz1988), respectively. The latter expression is used to model the ion Landau damping effect. In (3.4), $\nu _i$ is the ion collision frequency, $k_{\|}$ denotes the parallel wave vector with $k_{\|}\sim 1/R_0$ , $\rho _i$ is the ion mass density and $\kappa _{\|}$ is a strength parameter (Chu et al. Reference Chu, Greene, Jensen, Miller, Bondeson, Johnson and Mauel1995).

3.2.5. Ion diamagnetic drift velocity

Boozer magnetic flux coordinates ( $s,v,u$ ) (Boozer Reference Boozer1982) have the advantage that in these coordinates the poloidal and toroidal contravariant components of the diamagnetic drift velocity are constant on flux surfaces (Strumberger et al. Reference Strumberger, Günter, Lackner and Puchmayr2023). It is

(3.5) \begin{align} \boldsymbol{V}_{i0}^*&=\frac {m_i}{e}\frac {1}{\rho _0} \frac {\boldsymbol{B}_0 \times \boldsymbol{\nabla }p_{i0}}{B_0^2} =V^{v*}_{i0}{\boldsymbol r}_{,v} +V^{u*}_{i0}{\boldsymbol r}_{,u} \nonumber \\ &=\frac {m_i}{e}\frac {1}{\rho _0}\frac {\partial p_{i0}}{\partial s} \left (\frac {I_0}{\varPhi '_0J_0+\varPsi '_0I_0}\boldsymbol{r}_{,v} - \frac {J_0}{\varPhi '_0J_0+\varPsi '_0I_0}\boldsymbol{r}_{,u}\right ) = \frac {1}{2\pi } \varOmega _{tor}^{i*} {\boldsymbol r}_{,v} + \frac {1}{2\pi } \varOmega _{pol}^{i*} {\boldsymbol r}_{,u} \end{align}

with the toroidal and poloidal frequencies

(3.6) \begin{align} \varOmega _{tor}^{i*}(s)=2\pi \frac {m_i}{e}\frac {1}{\rho _0} \frac {\partial p_{i0}}{\partial s} \frac {I_0}{\varPhi '_0J_0+\varPsi '_0I_0}\,\,\textrm {and}\,\, \varOmega _{pol}^{i*}(s)=-2\pi \frac {m_i}{e}\frac {1}{\rho _0} \frac {\partial p_{i0}}{\partial s} \frac {J_0}{\varPhi '_0J_0+\varPsi '_0I_0}.\nonumber\\[1pt] \end{align}

Here, $\varPhi '_0$ and $\varPsi '_0$ are the derivatives of the toroidal $\varPhi _0(s)$ and poloidal $\varPsi _0(s)$ fluxes.

3.2.6. The ExB-velocity

Excluding radial flows which would be generated by non-radial components of the electric field $\boldsymbol{E}_0$ , we assume that the electric potential $\varPhi _E$ is constant on flux surfaces, that is,

(3.7) \begin{align} \boldsymbol{E}_0= -\frac {\partial \varPhi _E}{\partial s}\boldsymbol{\nabla }s = E_s \boldsymbol{\nabla }s \end{align}

with $E_s$ being the radial covariant component. Then, analogous to the ion diamagnetic drift velocity (3.5), the ExB-velocity can be expressed in Boozer coordinates by (Strumberger et al. Reference Strumberger, Puchmayr and Hindenlang2025)

(3.8) \begin{align} \boldsymbol{V}_{Ex B}=\frac {\boldsymbol{B}_0 \times \boldsymbol{\nabla }{\varPhi _E(s)}}{{B_0^2}}= \frac {1}{2\pi } \varOmega _{tor}^{ExB}(s)\boldsymbol{r}_{,v} + \frac {1}{2\pi } \varOmega _{pol}^{ExB}(s)\boldsymbol{r}_{,u}. \end{align}

Assuming $T_{i0} \approx T_{e0}$ and $\hat \nu _i /\sqrt {m_i} \approx \hat \nu _e /\sqrt {m_e}$ with $\hat \nu _i$ and $\hat \nu _e$ being the ion and electron collision frequencies normalised to the corresponding particle masses with $m_i \gg m_e$ , the well-known ambipolarity condition (Stroth, Manz & Ramisch Reference Stroth, Manz and Ramisch2011; Viezzer Reference Viezzer2012) simplifies to

(3.9) \begin{align} E_s \approx \frac {m_i}{e}\frac {1}{\rho _0}\frac {\partial p_{i0}}{\partial s}. \end{align}

In this case, the ion diamagnetic drift velocity is balanced by the ExB-velocity.

3.2.7. External velocity $\boldsymbol{V}^{ext}$

The general formulation of the velocity

(3.10) \begin{align} \boldsymbol{V}^{ext}&= \frac {1}{2\pi }\varOmega _{tor}(s)f_{vv}(s,v,u) \boldsymbol{r},_v+ \frac {1}{2\pi }\varOmega _{pol}(s)f_{vu}(s,v,u) \boldsymbol{r},_v\nonumber\\& \quad + \frac {1}{2\pi }\varOmega _{tor}(s)f_{uv}(s,v,u) \boldsymbol{r},_u+ \frac {1}{2\pi }\varOmega _{pol}(s)f_{uu}(s,v,u) \boldsymbol{r},_u \end{align}

was added to the ion velocity (2.6) because it allows us to express a given velocity in different kinds of flux coordinates.

Using Boozer coordinates, the ion drift velocity (3.5) and/or the ExB-velocity (3.8) can also be expressed in the $\boldsymbol{V}^{ext}$ -formulation (3.10) with $f_{vv}=f_{uu}=1$ and $f_{vu}=f_{uv}=0$ . This possibility allows to test CASTOR3D with respect to code errors concerning the implementation of these velocities.

Driven, for example, by neutral beam injection, tokamak plasmas can quickly rotate in toroidal direction. Using NEMEC or 2-D straight field line coordinates, $\boldsymbol{V}^{ext}$ simply reads

(3.11) \begin{align} \boldsymbol{V}^{ext}= \frac {1}{2\pi }\varOmega _{tor}(s) \boldsymbol{r},_v\!. \end{align}

However, if Boozer coordinates are used for a 3-D tokamak equilibrium, this velocity has to be described by toroidal and poloidal components ( $f_{vv}\ne 0$ and $f_{uv}\ne 0$ ) because $\boldsymbol{r},_v$ in Boozer coordinates is not parallel to $\boldsymbol{r},_v$ in NEMEC or 2-D straight field line coordinates. The transformation of the velocity into Boozer coordinates is performed with the COTRANS code.

Analogous to tokamaks, the plasma of a QA stellarator enjoys ‘intrinsic ambipolarity’ and is free to rotate in the direction of symmetry (Helander & Simakov Reference Helander and Simakov2008), which corresponds to the toroidal direction ( $v$ -direction) in Boozer coordinates. So, using Boozer coordinates, the externally driven flow can be expressed by (3.10).

3.2.8. Normalisation

The CASTOR3D code needs its input and yields its output in SI units. However, the solution of the eigenvalue problem is performed in normalised quantities. Thus, in the physics part, all input quantities are normalised to radius $R_a$ , magnetic field strength $B_a$ and density $\rho _a$ at the magnetic axis.

3.6.1. Determination of the mode type

Multiplying the components of the eigenvector (2.35) with the associated Hermite polynomials, we obtain the Fourier spectra

(3.17) \begin{align} \hat w_{{f},c}^{{m,n}}(s)=\sum _{{j,q}} c^{{m,n}}_{{f,j,q}} h^{{j}}_{{f,q}}(s) \,\,\,\,\,\textrm {and}\,\,\,\, \hat w_{{f},\check c}^{{m,n}}(s)=\sum _{{j,q}} \check c^{{m,n}}_{{f,j,q}} h^{{j}}_{{f,q}} (s) \end{align}

of the perturbations. These spectra reveal important information about location and structure of a mode. Usually, we consider the real and/or imaginary parts of the radial velocity perturbation vicariously for all perturbed quantities, that is, $\text{Re}(\hat v^s_{{m,n}}(s))=\text{Re}(\hat w_{{\hat v^s},c}^{{m,n}}(s))$ , $\text{Im}(\hat v^s_{{m,n}}(s))=\text{Im}(\hat w_{{\hat v^s},c}^{{m,n}}(s))$ , $\text{Re}(\hat {\bar v}^s_{{m,n}}(s))= \text{Re}(\hat w_{v^s,\check c}^{{m,n}}(s))$ and $\text{Im}(\hat {\bar v}^s_{{m,n}}(s))= \text{Im}(\hat w_{v^s,\check c}^{{m,n}}(s))$ . While in the case of axisymmetric configurations the modes can be easily characterised by their toroidal mode number ${n}$ , several ${n}$ -harmonics couple together in 3-D plasmas. Therefore, we introduced the dominant toroidal harmonic ${n}^*$ to name the modes (Strumberger & Günter Reference Strumberger and Günter2019). In order to determine $ n^*$ , we compute the following four functions:

(3.18) \begin{align} f_{\text{Re}}( n)=\sum _ m \int \big(\text{Re}\big(\hat v^s_{{m,n}}(s)\big)\big)^2 \text{ d}s, \,\,\,\,\,\,\,\,\, \tilde f_{\text{Re}}( n)=\sum _{{m}} \int \big(\text{Re}\big(\hat {\bar v}^s_{{m,n}}(s)\big)\big)^2 \text{ d}s, \end{align}

(3.19) \begin{align} f_{\text{Im}}( n)=\sum _ m \int \big(\text{Im}\big(\hat v^s_{{m,n}}(s)\big)\big)^2 \text{ d}s, \,\,\,\,\,\,\,\,\, \tilde f_{\text{Im}}( n)=\sum _ m \int \big(\text{Im}\big(\hat {\bar v}^s_{{m,n}}(s)\big)\big)^2 \text{ d}s \end{align}

and look for their maxima.

Figure 10.

Cosine-type external kink modes of a QA equilibrium. (a) Real parts of the eigenfunction Fourier spectrum for the radial velocity perturbation of the ${n}^*=1$ mode. (b) Contributions of the toroidal Fourier harmonics ( $f_{\text{Re}}$ (solid lines), $\tilde f_{\text{Re}}$ (symbols)) to the cosine-type mode for various values of the parallel viscosity coefficient $\mu _{\|}$ (from Strumberger et al. (Reference Strumberger, Puchmayr and Hindenlang2025); licensed under a Creative Commons Attribution licence (CC BY)).

Figure 10(a) shows the Fourier spectrum of the radial velocity perturbation of an external kink mode of a QA equilibrium with a parallel viscosity of $\mu _{\|}=3\text{ kg ms}^{-1}$ . Since, both, the n = 1 (green) and n = 5 (blue) harmonics contribute strongly to the mode, n $^*$ cannot be determined from this figure. In figure 10(b), however, the functions (3.17) indicate that n = 1 is the dominant harmonic in the case of $\mu _{\|}=3\text{ kg ms}^{-1}$ . Since, it is $\hat v^s_{{m,n}}(s)= \hat {\bar v}^s_{{m,n}}(s)$ in figure 10(a), this $n^*$ = 1 mode is a cosine-type mode, while sine-type modes are characterised by $\hat v^s_{{m,n}}(s)= -\hat {\bar v}^s_{{m,n}}(s)$ (see figure 26).

Besides the Fourier spectra, the diagnostic part of the CASTOR3D code also yields the mode structure in real space (e.g. see figures 12 c and 12 d, 15 b–d and 17 a and 17 b).

3.6.2. Analysis of the energetic drives of the modes

The energy functional of ideal MHD was extended to the regime of finite resistivity and transformed into an intuitive form (Puchmayr et al. Reference Puchmayr, Dunne and Strumberger2022) that separates pressure-gradient drive, current-density drive and stabilising phenomena analogously to the ideal version of the functional derived by Greene and Johnson (Greene & Johnson Reference Greene and Johnson1968). The resulting resistive energy functional is not in the form of a solvable eigenvalue problem, but is intended to analyse the energetic composition of eigenfunctions which were calculated in advance by solving the MHD eigenvalue problem (2.34) or the extended eigenvalue problem (2.89). The functional was decomposed into energy contributions related to the physical mechanisms which drive or stabilise the instabilities. The decomposition is based on the linearised resistive MHD equations, excluding equilibrium flows and any extended MHD effects,

(3.20) \begin{align} \rho _0 \lambda ^2 \boldsymbol{\xi } = \boldsymbol{\nabla }(\boldsymbol{\xi }\boldsymbol{\cdot }\boldsymbol{\nabla }p_0) + \boldsymbol{\nabla }(p_0 \varGamma (\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{\xi })) + \frac {1}{\mu _0} (\boldsymbol{\nabla }\times \boldsymbol{B}_0) \times \boldsymbol{B_1} + \frac {1}{\mu _0} (\boldsymbol{\nabla }\times \boldsymbol{B_1}) \times \boldsymbol{B}_0, \end{align}

(3.21) \begin{align} \boldsymbol{B_1} = \boldsymbol{\nabla }\times (\boldsymbol{\xi } \times \boldsymbol{B}_0 ) - \lambda ^{-1} \boldsymbol{\nabla }\times \bigg(\frac {\eta }{\mu _0} \boldsymbol{\nabla }\times \boldsymbol{B_1}\bigg), \end{align}

where $\boldsymbol{\xi }=\lambda ^{-1}\boldsymbol{v_1}$ is the displacement. The initial form of the resulting energetic decomposition is similar to the energy functional of Bernstein et al. (Reference Bernstein, Frieman, Kruskal and Kulsrud1958),

(3.22) \begin{align} - \lambda ^2 \int \rho _0 |\boldsymbol{\xi }|^2 \text{dV} = 2 \delta W = 2 \delta W_V + 2 \delta W_B \end{align}

with the volume energy contributions $\delta W_V$ and boundary contributions $\delta W_B$ . The left-hand side of (3.22) is related to the kinetic energy of the instability, while the right-hand side $2\delta W = 2 \delta W_V + 2 \delta W_B$ describes the change of the potential energy caused by the perturbation of the equilibrium.

The derivation of the intuitive form of the resistive energy functional is described in detail in Puchmayr et al. (Reference Puchmayr, Dunne and Strumberger2022). Here, we only present the intuitive form of the volume energy contribution. It reads for resistive MHD,

(3.23) \begin{align} 2 \delta W_V&= \underbrace {\int \vphantom {\frac {1}{\boldsymbol{B}_0^2}} \frac {| \boldsymbol{B}_{1,\perp } |^2}{\mu _0} \text{dV}}_{\delta W_{\text{SHA}}} + \underbrace {\int \frac {\boldsymbol{B}_0^2}{\mu _0} \Big \vert \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{\xi }_\perp + 2 (\boldsymbol{\xi }_\perp \boldsymbol{\cdot }\boldsymbol{\kappa }) + \frac {1}{B_0^2}\boldsymbol{B}_0 \boldsymbol{\cdot }\boldsymbol{B}_R \Big \vert ^2 \text{dV}}_{\delta W_{\text{CPA}}} \nonumber \\&\quad +\underbrace {\varGamma \int p_0 \vert \boldsymbol{\nabla } \boldsymbol{\cdot }\boldsymbol{\xi }\vert ^2 \text{dV}}_{\delta W_{\text{SND}}} +\underbrace {\frac {1}{\mu _0^2 \lambda ^*} \int \eta | \boldsymbol{\nabla }\times \boldsymbol{B}_1 |^2 \text{dV}}_{\delta W_{\text{RCD}}} -\underbrace {\int \frac {j_{0,\|}} {| \boldsymbol{B}_0 |} (\boldsymbol{\xi }_\perp ^\ast \times \boldsymbol{B}_0) \boldsymbol{\cdot }\boldsymbol{B}_{1,\perp } \text{dV}}_{\delta W_{\text{CUR}}} \nonumber \\ &\quad -\underbrace {\int 2 (\boldsymbol{\xi }_\perp \boldsymbol{\cdot } \boldsymbol{\nabla }p_0) (\boldsymbol{\xi }_\perp ^\ast \boldsymbol{\cdot }\boldsymbol{\kappa }) \text{dV}}_{\delta W_{\text{DP}}} -\underbrace {\int \frac {1}{B_0^2} (\boldsymbol{\xi }_\|^\ast \times \boldsymbol{\nabla }p_0) \boldsymbol{\cdot }(\boldsymbol{B}_0 \times \boldsymbol{B}_R) \text{dV}}_{\delta W_{\text{RD,}\|}} \nonumber \\ &\quad -\underbrace {\int \frac {1}{B_0^2} (\boldsymbol{\xi }_\perp \boldsymbol{\cdot }\boldsymbol{\nabla }p_0) (\boldsymbol{B}_0 \boldsymbol{\cdot }\boldsymbol{B}_R^\ast ) \text{dV}}_{\delta W_{\text{RD,}\perp }} \end{align}

with the self-induced ohmic field $\boldsymbol{B}_R = \lambda ^{-1} \boldsymbol{\nabla } \times (\eta \boldsymbol{j}_1)$ and the equilibrium curvature vector $\boldsymbol{\kappa } = (\boldsymbol{b} \boldsymbol{\cdot } \boldsymbol{\nabla })\boldsymbol{b}$ with $\boldsymbol{b}=\boldsymbol{B}_0/|\boldsymbol{B}_0|$ . There are four stabilising contributions to the volume energy which are related to perpendicular magnetic perturbations (shear Alfvén waves) $\delta W_{\text{SHA}}$ , parallel magnetocompressional perturbations (compressional Alfvén waves) $\delta W_{\text{CPA}}$ , adiabatic compression (sound waves) $\delta W_{\text{SND}}$ and resistive current diffusion $\delta W_{\text{RCD}}$ . The destabilising terms are the parallel current-density drive $\delta W_{\text{CUR}}$ , the field-line-curvature-dependent pressure-gradient $\delta W_{\text{DP}}$ drive and two resistive corrections to the pressure-gradient drive. The latter two account for the effect of resistive diffusion on the work done by the induced magnetic perturbation that is generated by parallel $\delta W_{RD,\|}$ and perpendicular $\delta W_{RD,\perp }$ compression of the plasma, respectively. It should be noted that the resistive correction to the relation between perpendicular compression and its induced magnetic perturbation also appears in the stabilising energy of magnetocompressional perturbations $\delta W_{\text{CPA}}$ . In order to obtain energetic contributions which are independent of the arbitrary linear amplitude of the perturbation, it is natural to normalise the energy terms by

(3.24) \begin{align} \delta W_{\text{norm}} \equiv \int \rho _0 |\boldsymbol{\xi }|^2 \text{dV}. \end{align}

Concerning the boundary contributions $\delta W_B$ , it should be mentioned that no additional surface or vacuum energy contributions have to be considered in the case of finite resistivity.

The terms of the resistive energy functional were implemented in the CASTOR3D code. The energetic decomposition contributes to the physical understanding of the perturbations (see § 4.3) and allows to distinguish current-density and pressure-gradient driven instabilities (see § 4.1.2).

4.1.1. Ideal modes

Figure 11.

Growth rate $\gamma$ and absolute value of the oscillation frequency $|\omega |$ as function of (a) the toroidal mode number n and (b) the scaling parameter $\alpha _s$ (see § 2.1.1). The black dashed line represents the ideal growth rate $\gamma _{ideal}$ obtained without taking the diamagnetic drift effects into account. (a) The blue stars mark the absolute value of half of the ion diamagnetic drift frequency ${|\omega ^*_i|}/{2}$ (4.2). (From Strumberger et al. (Reference Strumberger, Günter, Lackner and Puchmayr2023); licensed under a Creative Commons Attribution licence (CC BY).)

Figure 11. Long description

The image contains two line graphs labeled as Panel A and Panel B. Panel A shows the growth rate and absolute value of the oscillation frequency as a function of the toroidal mode number n. The vertical axis represents the growth rate and oscillation frequency in units of 1 per second, and the horizontal axis represents the toroidal mode number n. The black dashed line represents the ideal growth rate obtained without considering diamagnetic drift effects. The red line with red dots represents the growth rate, and the blue line with blue stars marks the absolute value of half of the ion diamagnetic drift frequency. Panel B shows the growth rate and absolute value of the oscillation frequency as a function of the scaling parameter. The vertical axis represents the growth rate and oscillation frequency in units of 1 per second, and the horizontal axis represents the scaling parameter. The black dashed line represents the ideal growth rate, the red line with red dots represents the growth rate, and the blue line with blue stars marks the absolute value of half of the ion diamagnetic drift frequency. The graphs illustrate how the growth rate and oscillation frequency vary with the toroidal mode number and scaling parameter, highlighting the differences between the ideal growth rate and the growth rate considering diamagnetic drift effects.

The considered axisymmetric, high- $\beta$ equilibrium (volume averaged plasma beta: $\langle \beta \rangle =2.97\,\,\%$ ) was ideal unstable. Figure 11 illustrates the strong stabilising impact and the oscillation frequency mainly caused by the gyroviscosity. Because of the high plasma- $\beta$ , the contribution of the electron diamagnetic drift effect was negligibly small as shown in Strumberger et al. (Reference Strumberger, Günter, Lackner and Puchmayr2023). We used $\alpha _T=0.4$ , that is, $p_i=0.4p$ and $p_e=0.6p$ . In good agreement with the results of the MINERVA-DI code and the dispersion relation

(4.1) \begin{align} \gamma +\text{i}\omega =\text{i} \frac {|\omega ^*_i|}{2} \pm \sqrt {\gamma ^2_{ideal}- \left (\frac {\omega ^*_i}{2}\right )^2} \end{align}

presented by Aiba (Reference Aiba2016), the modes became stabilised when the absolute value of the oscillation frequency $|\omega |$ of the mode exceeded the ideal growth rate, $|\omega |\gtrapprox \gamma _{ideal}$ . The modes were stable for $n\gt 8$ as shown in figure 11. In the case of simple equilibria where mainly one poloidal harmonic m and one toroidal harmonic ${n}$ contribute to the Fourier spectrum of the mode, the ion diamagnetic drift frequency $\omega ^*_i=\boldsymbol V^*_i \boldsymbol{\cdot }\boldsymbol{k}$ ( $\boldsymbol{k}$ = wave number vector) can be approximated by $\boldsymbol{k}=-\text{i}(( { \boldsymbol{\nabla }f})/{f})$ with function $f=f_0(s)\exp (2\pi \text{i}({m}u-{n}v))$ . The latter describes the perturbation caused by the leading $({m,n})$ -Fourier harmonic. Expressing $\boldsymbol{V^*_i}$ in Boozer coordinates (3.5), we obtained

(4.2) \begin{align} \omega ^*_i(s)=2\pi ({\mathrm{m}V^{*u}_i-\mathrm{n}V^{*v}_i)} \end{align}

for the ion diamagnetic drift frequency. In figure 11(a) the blue stars mark $({|\omega ^*_i|}/{2})$ calculated with the analytic expression (4.2) at the rational $q=1$ flux surface. These values agree very well with the oscillation frequencies resulting from the solutions of the eigenvalue problem.

Besides the axisymmetric equilibrium, we also investigated, for the first time, diamagnetic drift effects in an equivalent one-periodic, 3-D equilibrium with helical core (Strumberger et al. Reference Strumberger, Günter, Lackner and Puchmayr2023). This equilibrium was the result of a 3-D fixed-boundary NEMEC equilibrium calculation using an axisymmetric, circular plasma boundary (Strumberger & Günter Reference Strumberger and Günter2017). Neglecting the pressure terms in Ohm’s law (2.16), the results obtained for ideal kink and interchange modes confirmed the growing stabilising diamagnetic drift effect with increasing toroidal mode number (which is the dominant toroidal harmonic ${n}^*$ ). Furthermore, the oscillation frequencies of the modes corresponded in good approximation to $({\omega ^*_i}/{2})$ . However, additionally to these similarities with the axisymmetric equilibrium, 3-D specific effects were found: sufficiently strong diamagnetic drift effects are needed (i) to cause a degeneration of otherwise non-degenerate growth rates and (ii) to let the mode oscillate (see also § 4.3). (iii) Diamagnetic drift effects change the mode structure as shown in figure 12. The figure illustrates the changes of the composition of the Fourier spectrum and the mode location for the ${n}^*=6$ mode. Here, we only show the radial velocity perturbation of the cosine-type function because the sine-type function yielded no new information. The components related with the complex and complex conjugate $\text{e}$ -functions (3.17) of the Fourier series are equal (figure 12 a) without diamagnetic drift effects. However, including diamagnetic drift effects, either the complex or complex conjugate components strongly dominate the Fourier spectrum (figure 12 b). Here, it should be mentioned that not only drift effects lead to such a decoupling, but any kind of plasma flow (Puchmayr et al. Reference Puchmayr2023; Strumberger et al. Reference Strumberger, Puchmayr and Hindenlang2025). Furthermore, besides the reduction of the growth rate and the decoupling of the the complex and complex conjugate parts, the number of contributing toroidal harmonics is reduced. This leads to a more uniformly distributed mode structure in the real space as shown in figures 12(c) and 12(d). Here, the coupling of the poloidal harmonics played only a minor role for the mode structure because of the very flat $q$ -profile and the circular plasma boundary.

Figure 12.

(a,b) Fourier spectra of the real part of the radial velocity perturbation, and (c,d) mode structures of the ${n}^*=6$ mode (a,c) without and (b,d) with ion diamagnetic drift effect. In (a) and (b), the solid and dashed lines denote the Fourier spectra of the radial velocity perturbation related with the complex and complex conjugate $\text{e}$ -functions. The largest contributions are marked by their poloidal and toroidal harmonics ${m/n}$ . In (c) and (d), the vectors characterise the velocity perturbation in the $R$ – $Z$ -plane at toroidal angle $\varphi =0^{\textrm {o}}$ and time $t=0$ . (From Strumberger et al. (Reference Strumberger, Günter, Lackner and Puchmayr2023); licensed under a Creative Commons Attribution licence (CC BY).)

4.1.2. Resistive modes

Figure 13.

Energetic decomposition of the resistive internal perturbations of the considered low- $\beta$ equilibrium: (a) relative current-density drive as function of the toroidal mode number for resistivities of $\eta =10^{-7}$ Ωm (blue) and $\eta =10^{-6}$ Ωm (orange). Normalised energy densities corresponding to (b,d) the current-density and (c) the pressure-gradient drives as functions of the square root of the normalised toroidal flux for resistivities of (b) $\eta =10^{-7}$ Ωm and (c,d) $\eta =10^{-6}$ Ωm. (From Strumberger et al. (Reference Strumberger, Günter, Lackner and Puchmayr2023); licensed under a Creative Commons Attribution licence (CC BY).)

Figure 13. Long description

The image contains four line graphs analyzing the energetic decomposition of resistive internal perturbations in a low-beta equilibrium. Panel A shows the relative current-density drive as a function of the toroidal mode number for resistivities of 1 times 10 to the power of -7 ohm meters (blue) and 1 times 10 to the power of -6 ohm meters (orange). The x-axis represents the toroidal mode number (n) ranging from 1 to 15, and the y-axis represents the relative current-density drive in percent. Both lines show a decreasing trend as the toroidal mode number increases, with the blue line starting higher and converging towards the orange line. Panel B displays the normalized energy density corresponding to the current-density drive as a function of the square root of the normalized toroidal flux for a resistivity of 1 times 10 to the power of -7 ohm meters. The x-axis ranges from 0.28 to 0.40, and the y-axis ranges from -0.5 to 1.5 times 10 to the power of -4. Multiple lines represent different toroidal mode numbers (n=1, 2, 3, 5, 10, 15), showing various peaks and troughs around the square root of the normalized toroidal flux of 0.36. Panel C shows the normalized energy density corresponding to the pressure-gradient drive as a function of the square root of the normalized toroidal flux for a resistivity of 1 times 10 to the power of -6 ohm meters. The x-axis ranges from 0.28 to 0.40, and the y-axis ranges from -1.2 to 0 times 10 to the power of -4. Multiple lines represent different toroidal mode numbers (n=1, 2, 3, 5, 10, 15), showing a significant trough around the square root of the normalized toroidal flux of 0.36. Panel D displays the normalized energy density corresponding to the current-density drive as a function of the square root of the normalized toroidal flux for a resistivity of 1 times 10 to the power of -6 ohm meters. The x-axis ranges from 0.28 to 0.40, and the y-axis ranges from -0.5 to 1.5 times 10 to the power of -4. Multiple lines represent different toroidal mode numbers (n=1, 2, 3, 5, 10, 15), showing various peaks and troughs around the square root of the normalized toroidal flux of 0.36.

Figure 14.

(a) Growth rates and (b) oscillation frequencies as functions of parameter $\alpha _T$ for various n-harmonics. The dotted vertical lines mark the $\alpha _T$ -values for which the considered resistivity value corresponds to the approximated Spitzer resistivity. The stars denote $\omega _i^*$ . (From Strumberger et al. (Reference Strumberger, Günter, Lackner and Puchmayr2023); licensed under a Creative Commons Attribution licence (CC BY).)

In order to study the diamagnetic effects on resistive modes, the same axisymmetric plasma configuration was used as for the ideal modes, but with a much lower plasma- $\beta$ of $\langle \beta \rangle =0.97\,\%$ (Strumberger et al. Reference Strumberger, Günter, Lackner and Puchmayr2023). This equilibrium was unstable with respect to internal resistive modes located at the $q=1$ flux surface. Their slow growth rates approximately scaled with $\eta ^{1/3}$ which is typical for resistive kink and interchange modes. The analysis of the resistive energy functional (§ 3.6.2) prevailed the pressure-gradient driven character of these modes. Figure 13(a) shows the relative current-density drive $\delta W_{CUR}/(\delta W_{CUR}+\delta W_{DP})$ as function of the toroidal mode number. It adds up to 3.5 % for ${n}=1$ and $\eta =1 \times10^{-7}$ Ωm and strongly decreases with increasing n. Although the relative current-density drive is nearly independent of the resistivity, the perturbation and the corresponding energy densities are broadened with increasing resistivity (see figures 13 b and 13 d). For increasing mode numbers, however, the modes become more and more localised at the $\iota =1$ flux surface (see figure 13 d). While the integrated current-density drive $\delta W_{\text{CUR}}$ (3.23) remains small because of the cancellation of stabilising and destabilising contributions in different radial regions of the plasma (see figures 13 c and 13 d), the pressure-gradient drive $\delta W_{\text{DP}}$ (3.23) strongly grows with increasing n.

Figures 14(a) and 14(b) show the growth rates and oscillation frequencies as functions of the parameter $\alpha _T$ for various ${n}$ -modes. They include the limiting cases of cold ions ( $\alpha _T=0$ ) and cold electrons ( $\alpha _T=1$ ). In the case of cold ions, only the electron diamagnetic drift remains that is described by the parallel electron pressure-gradient term in the extended Ohm’s law (2.16), while, in the case of cold electrons, only the ion diamagnetic drift influences the stability via the gyroviscosity term in the momentum equation (2.9) and the ion pressure-gradient term in (2.16). All combinations of ion and electron diamagnetic drift effects have a large stabilising effect (Strumberger et al. Reference Strumberger, Günter, Lackner and Puchmayr2023).

The oscillation frequencies of the resistive modes scale linearly with $\alpha _T$ as shown in figure 14(b). Computing the ion diamagnetic drift frequency (4.2) for the toroidal mode numbers and a chosen $\alpha _T$ -parameter (indicated by the dotted lines in the figures 14 a and 14 b), it turned out that these frequencies agree very well with the oscillation frequencies resulting from the solution of the eigenvalue problem, that is, $\omega \approx \omega _i^*$ . In the case of ideal modes, however, it is $\omega \approx \omega _i^*/2$ as shown in § 4.1.1.

4.4.1. The ELMs in 2-D and 3-D AUG tokamak plasmas

Figure 22.

Toroidal cross-sections (a) for the axisymmetric and (b) magnetically perturbed equilibria of the AUG discharges $\#$ 40180 (solid, black) and $\#$ 40181 (dashed, orange) as well as (c) pressure and safety factor profiles. For visibility, the corrugation in (b) is amplified by a factor of 2.5. (d) Growth rates for a range of equilibria between AUG $\#$ 40180 and AUG $\#$ 40181 with and without ion diamagnetic drift effects. The numbers annotating the data points refer to the toroidal harmonic ${n}^*$ of the most unstable instability. The dotted blue line indicates the empirical stability limit from Suttrop et al. (Reference Suttrop2018) assuming equal ion and electron pressures. (From Puchmayr et al. (Reference Puchmayr, Dunne, Strumberger, Willensdorfer, Zohm and Hindenlang2024); licensed under a Creative Commons Attribution 4.0 licence (CC BY).)

Type I ELMs limit the pressure gradient and lead to significant energy and particle losses in tokamak plasmas. In order to avoid severe damage in future fusion devices (e.g. ITER (Loarte et al. Reference Loarte2003)), they have to be suppressed or, at least, mitigated. The application of non-axisymmetric MP fields in various experiments (Evans et al. Reference Evans2004; Nazikian et al. Reference Nazikian2015; Suttrop et al. Reference Suttrop2018) led to a mitigation or even suppression of ELMs because of an enhanced transport which reduces the pedestal density and pressure gradient (Kirk et al. Reference Kirk2011; Suttrop et al. Reference Suttrop2018). Since the onset of Type I ELMs is very well described by the growth of edge-localised peeling–ballooning modes (Snyder et al. Reference Snyder, Wilson, Ferron, Lao, Leonard, Osborne, Turnbull, Mossessian, Murakami and Xu2002; Wilson et al. Reference Wilson, Cowley, Kirk and Snyder2006), a predictive linear stability analysis of ELMs is possible.

Puchmayr (Reference Puchmayr2023) and Puchmayr et al. (Reference Puchmayr, Dunne, Strumberger, Willensdorfer, Zohm and Hindenlang2024), the impact of symmetry-breaking induced by MP fields on the MHD stability limit was analysed for low-n kink-peeling mode up to high-n peeling–ballooning mode limited regimes considering, for the first time, ion diamagnetic drift effects in realistic non-axisymmetric (3-D) tokamak geometry. The application of MP fields leads to a corrugation of the equilibrium flux surfaces (Chapman et al. Reference Chapman2012; Willensdorfer et al. Reference Willensdorfer2019). Figures 22(a) and 22(b) show the axisymmetric and the corrugated 3-D plasma boundaries of the AUG discharges $\#$ 40180 (ELM suppression) and $\#$ 40181 (ELM mitigation).Footnote ²⁴ The maximum corrugation amounts to 7 mm for the ELM suppression and 13 mm for the ELM mitigation case at the low-field side. Their pressure and safety factor profiles are illustrated in figure 22(c). In order to analyse the transition from the ELM mitigation to the ELM suppression case, a range of equilibria with increasing pedestal top pressure was created inbetween these two cases. Figure 22(d) shows the normalised growth rate as function of the pedestal top pressure for axisymmetric and 3-D equilibria with and without diamagnetic drift effects. The figure illustrates the decrease of the critical pedestal top pressure caused by the symmetry-breaking effect of the applied MP fields. In contrast to the experimental observations, the 2-D and 3-D numerical results predict highly unstable ballooning modes with ${n}^*=15{-}19$ if no ion diamagnetic drift effect is taken into account. Inclusion of this effect leads to a considerable stabilisation and a reduction of the toroidal mode number to ${n}^*=5$ or 6. Here, the peeling–ballooning modes are suppressed below the empirical MP stability limit and the finite growth rates correspond to mitigated peeling–ballooning modes. These results demonstrate the importance of 3-D stability studies and the inclusion of diamagnetic drift effects.

4.4.2. External kink and peeling–ballooning modes in a QA stellarator

Figure 23.

(a) Growth rates as functions of ${n}^*$ without (black circles) and with (red circles) taking parallel viscosity into account. (b) Growth rates and oscillation frequencies as functions of ${n}^*$ taking parallel and gyroviscosity into account. The results with and without additional ion diamagnetic drift and ExB-velocities are marked by red/green and black/blue circles, respectively. (From Strumberger et al. (Reference Strumberger, Puchmayr and Hindenlang2025); licensed under a Creative Commons Attribution licence (CC BY).)

Figure 23. Long description

Two scatter plots are presented. Panel A shows growth rates as functions of n* without and with taking parallel viscosity into account. The x-axis represents n* and the y-axis represents growth rates in units of 1 per second. Black circles represent data without parallel viscosity, while red circles represent data with parallel viscosity. Panel B shows growth rates and oscillation frequencies as functions of n* taking parallel and gyroviscosity into account. The x-axis represents n* and the y-axis represents growth rates in units of 1 per second and oscillation frequencies in units of radians per second. Red and green circles represent results with additional ion diamagnetic drift and ExB-velocities, while black and blue circles represent results without these factors.

In § 4.3, we described low- ${n}^*$ locked modes in a two-periodic QA equilibrium and the effect of an external flow on these modes. Strumberger et al. (Reference Strumberger, Puchmayr and Hindenlang2025), the investigations of the same equilibrium were extended from low- ${n}^*$ external kink modes to high- ${n}^*$ peeling–ballooning modes. In addition to parallel viscosity and an external toroidal flow in symmetry direction, gyroviscosity, ion diamagnetic drift velocity and ExB-velocity were also considered in this work. The growth rates as function of the dominant toroidal mode number ${n}^*$ are illustrated in figure 23(a,b). Figure 23(a) compares the growth rates obtained without and with a parallel viscosity of $\mu _{\|}=4\text{ kg ms}^{-1}$ .Footnote ²⁵ The parallel viscosity has a stabilising effect which increases with ${n}^*$ . The high- ${n}^*$ peeling–ballooning modes are clearly less unstable if parallel viscosity is included. Nevertheless, the high- ${n}^*$ peeling–ballooning modes are still the most unstable modes.

Assuming that there is no external flow drive and that the ambipolarity condition (3.9) holds ( ${\boldsymbol V}_{ExB}=-{\boldsymbol V}_{i0}^*$ ), all terms containing ${\boldsymbol V}_0$ disappear in the linearised MHD equations (2.18–2.21), and the momentum equation (2.19) reads

(4.4) \begin{align} \rho _0 \lambda {\boldsymbol V}_1&= - \frac {1}{m_i}{\boldsymbol{\nabla }}(\rho _0 T_1) - \frac {1}{m_i}{\boldsymbol{\nabla }}(T_0 \rho _1) + \rho _0({\boldsymbol V}_{i0}^* \boldsymbol{\cdot }{ \boldsymbol{\nabla }}) {\boldsymbol V}_1 +\mu _{||} \Delta _{||} {\boldsymbol V}_1 \nonumber \\[5pt] &\quad + \frac {1}{\mu _0} ({\boldsymbol{\nabla }} \times ({\boldsymbol{\nabla }} \times {\boldsymbol A}_1) ) \times ({\boldsymbol{\nabla }} \times {\boldsymbol A}_0) + \frac {1}{\mu _0} ({ \boldsymbol{\nabla }} \times ({\boldsymbol{\nabla }} \times {\boldsymbol A}_0) ) \times ({\boldsymbol{\nabla }} \times {\boldsymbol A}_1). \end{align}

While the parallel viscosity only reduces the growth rate, but does not cause mode oscillations, the gyroviscosity leads to oscillations. Figure 23(b) compares growth rates and oscillation frequencies obtained without and with including ion diamagnetic drift and ExB-velocities additionally to parallel and gyroviscosity. Here, the ion diamagnetic drift velocity and the ExB-velocity increase the oscillation frequency by more than a factor of 2, but they have a negligible effect on the growth rate. However, all effects together stabilise the high- ${n}^*$ peeling–ballooning modes so much that the ${n}^*=2$ external kink mode becomes the most unstable mode. Analogous to the external flow, gyroviscosity, drift and ExB-velocity are able to unlock modes (Strumberger et al. Reference Strumberger, Puchmayr and Hindenlang2025) and can lead to overstable modes as will be described in § 4.5.2.

Interesting is also the combination of external flow, gyroviscosity, ion diamagnetic drift velocity and ExB-velocity. In Strumberger et al. (Reference Strumberger, Puchmayr and Hindenlang2025), we considered a toroidal velocity, $\boldsymbol{V}^{ext}_{tor}=({1}/({2\pi })) \varOmega ^{ext}_{tor}(s)\boldsymbol{r},_v$ , externally driven in the symmetry direction $v$ (Boozer toroidal coordinate). We assumed $\varOmega ^{ext}_{tor}(s)= \varOmega _0^{ext}(1-0.8 s^4)$ with $-25\,000 \leqslant \varOmega _0^{ext} \leqslant 25\,000\text{ rad s}^{-1}$ . Figure 24(a–d) show growth rates and oscillation frequencies of ${n}^*=2$ , 3 and 7 external kink modes and ${n}^*=22$ and 47 peeling–ballooning modes as functions of the externally driven flow frequency $\varOmega _0^{ext}$ . Depending on the sign of $\varOmega ^{ext}_0$ , the additional external drive increases or decreases the damping effects caused by gyroviscosity, ion diamagnetic drift velocity and ExB-velocity and also changes the oscillation frequencies of the modes. The growth rate decreases if the chosen flow direction leads to an increase of the oscillation frequency of the mode and vice versa. This behaviour is in good agreement with the results obtained by Aiba (Reference Aiba2016). Using the MINERVA-DI code, he investigated axisymmetric equilibria with circular cross-section. Yet, an additional effect appears in the case of the considered 3-D equilibria. Appropriate external flow compensates the oscillation frequencies caused by gyroviscosity, diamagnetic drift velocity and ExB-velocity and, therefore, locks the low- ${n}^*$ modes. While the oscillation frequencies become zero, the growth rates split into two values.

Finally, the presented results show that direction and magnitude of the external flow may decide whether a low- ${n}^*$ external kink or a high- ${n}^*$ peeling–ballooning mode will be the most unstable mode. In figure 24(c), the black stars mark the growth rates of the $n^*=2$ mode for $\varOmega ^{ext}_0= -25\,000\text{ rad s}^{-1}$ and $25\,000\text{ rad s}^{-1}$ , respectively. Here, the growth rate of the ${n}^*=22$ mode has almost reached the value of the ${n}^*=2$ mode for $\varOmega ^{ext}_0=-25\,000\text{ rad s}^{-1}$ . That is, a sufficiently fast external flow in negative $v$ -direction is able to compensate the stabilising effects of gyroviscosity, ion diamagnetic velocity and ExB-velocity on the high- ${n}^*$ peeling–ballooning modes and to make them to the most unstable modes of the considered equilibrium, while the growth rates of low- ${n}^*$ modes are only weakly affected by the flows (see figures 23 b and 24 a).

Figure 24.

(a,c) Growth rates and (b,d) oscillation frequencies of (a,b) ${n}^*=2$ , 3 and 7 external kink modes, and (c,d) ${n}^*=22$ and 47 peeling–ballooning modes as functions of the externally driven toroidal flow frequency $\varOmega ^{ext}_{0}$ taking parallel viscosity ( $\mu _{||}=4\text{ kg ms}^{-1}$ ), gyroviscosity, ion diamagnetic drift velocity and ExB-velocity into account. The solutions were obtained solving (a,b) the full eigenvalue problem and (c,d) the reduced eigenvalue problem. (c) The two black stars located in the left- and right-hand corners mark the growth rates of the ${n}^*=2$ mode. (From Strumberger et al. (Reference Strumberger, Puchmayr and Hindenlang2025); licensed under a Creative Commons Attribution licence (CC BY).)

4.5.1. Axisymmetric equilibria

In axisymmetric equilibria, an overstable mode can occur if, at least, two unstable modes of the same n-type exist and their growth rates overlap within certain ranges of plasma parameters. When the growth rates as functions of a single parameter (e.g. pressure gradient (Kerner et al. Reference Kerner, Jakoby and Lerbinger1986), resistivity (Huysmans et al. Reference Huysmans, Goedbloed and Kerner1993), etc.) become the same, an overstable oscillating mode evolves. We will illustrate this process for an example of DTMs. These modes appear in the case of reversed shear $q$ -profiles (e.g. Günter et al. Reference Günter, Schade, Maraschek, Pinches, Strumberger, Wolf and Yu2000). Mostly, the two tearing modes forming around the positions of equal $q$ -values couple asymmetrically, but a symmetric coupling is also possible (Otto & Birk Reference Otto and Birk1992).

Figures 25(a) and 25(b) show the growth rates and oscillation frequencies in dependence of the plasma resistivity for 2/1 DTMs of an AUG-type equilibrium. The two $q=2$ surfaces being located at $s_1=0.36$ and $s_2=0.52$ have a distance of $\Delta s=0.16$ . While the growth rate of the asymmetric DTM becomes smaller with decreasing resistivity, the growth rate of the symmetric DTM increases. When the growth rates of both modes come close, an oscillating overstable mode forms. Then, the eigenvalue problem yields two overstable solutions with $\gamma _1=\gamma _2$ and $\omega _1=-\omega _2$ . Both solutions describe the same mode with the opposite signs of the frequencies indicating a phase shift in the oscillation cycle. Figures 25(c) and 25(d) illustrate the eigenfunction Fourier spectra of the radial velocity perturbations of the asymmetric- and the symmetric-coupled DTM for $\eta = 6\times10^{-7}$ Ωm.

Figure 25.

Double tearing modes of an AUG-type equilibrium. (a) Growth rates and (b) oscillation frequencies as functions of the resistivity. (c,d) Real parts of the eigenfunction Fourier spectra of the radial velocity perturbation for (c) an asymmetric- and (d) a symmetric-coupled DTM for $\eta = 6\times10^{-7}$ Ωm.

Figure 25. Long description

Panel A: A line graph shows the growth rates as a function of resistivity. The x-axis represents resistivity in ohm meters, ranging from 10^-8 to 10^-4. The y-axis represents growth rates in per second, ranging from 10^0 to 10^5. The graph includes three distinct curves labeled asymmetric DTM, overstable DTMs, and symmetric DTM. The asymmetric DTM curve shows an increasing trend, the overstable DTMs curve shows a peak and then a decline, and the symmetric DTM curve shows a decreasing trend. Panel B: Another line graph displays oscillation frequencies as a function of resistivity. The x-axis represents resistivity in ohm meters, ranging from 10^-8 to 10^-4. The y-axis represents oscillation frequencies in radians per second, ranging from -800 to 800. The graph includes three distinct curves labeled asymmetric and symmetric DTMs, and overstable DTMs. The asymmetric and symmetric DTMs curve shows a peak and then stabilizes, while the overstable DTMs curve shows a flat line. Panel C: A line graph shows the real parts of the eigenfunction Fourier spectra of the radial velocity perturbation for an asymmetric DTM. The x-axis represents a dimensionless parameter s, ranging from 0 to 1. The y-axis represents the real part of the eigenfunction in arbitrary units. The graph includes multiple colored lines labeled with different ratios, such as 2/1, 3/1, 4/1, and 5/1. The lines show various peaks and troughs. Panel D: Another line graph shows the real parts of the eigenfunction Fourier spectra of the radial velocity perturbation for a symmetric DTM. The x-axis represents a dimensionless parameter s, ranging from 0 to 1. The y-axis represents the real part of the eigenfunction in arbitrary units. The graph includes multiple colored lines labeled with different ratios, such as 2/1, 3/1, 4/1, and 5/1. The lines show various peaks and troughs.

Figure 26.

Overstable ${n}^*=1\leftrightarrow {n}^*=9$ external kink mode of a QA equilibrium: (a) real and (b) imaginary part of the eigenfunction Fourier spectrum for the radial velocity perturbation of the sine-type overstable mode (parallel viscosity: $\mu _{||}=1.1\text{ kg ms}^{-1}$ ). (From Strumberger et al. (Reference Strumberger, Puchmayr and Hindenlang2025); licensed under a Creative Commons Attribution licence (CC BY).)

4.5.2. The 3-D equilibria

In contrast to axisymmetric equilibria, overstable modes can also evolve if the growth rates of different ${n}^*$ -type modes of the same mode family overlap within a certain parameter range because the toroidal harmonics are coupled in 3-D equilibria. Investigating the growth rates of internal resistive modes (Strumberger et al. Reference Strumberger and Günter2020) and ideal external modes (Strumberger et al. Reference Strumberger, Puchmayr and Hindenlang2025) as functions of resistivity, viscosity and plasma flows, several overstable modes between various ${n}^*$ -types were found. Considering the QA equilibrium discussed in Strumberger & Günter (Reference Strumberger and Günter2019) and Strumberger et al. (Reference Strumberger, Puchmayr and Hindenlang2025), figure 26 illustrates the oscillation of an overstable ${n}^*=1\leftrightarrow {n}^*=9$ external kink mode that appeared within a certain range of the parallel viscosity. Here, the real part of the Fourier spectrum shows an ${n}^*=9$ structure, while the imaginary part of the Fourier spectrum corresponds to an ${n}^*$ = 1 mode. That is, the overstable mode oscillates between these two structures.