On the relationship between the multi-region relaxed variational principle and resistive inner-layer theory

Resistive tearing instabilities are of considerable fundamental as well as practical importance in the study of sawtooth oscillations. These oscillations arise from current-induced instabilities that can be caused either internally or externally, such as by electron cyclotron resonance heating [1]. Coppi–Greene–Johnson (CGJ) [2] and Glasser–Greene–Johnson (GGJ) [3] demonstrated that the finite pressure has a stabilizing effect on tearing instabilities. In their magnetohydrodynamic (MHD) theory, the resistive modes were analyzed using the boundary layer theory when the boundary layers occur near the rational surfaces. Basically, the plasma can be analyzed in two regions: an 'outer region', where the plasma is ideal, and an 'inner layer', where non-ideal dissipation is significant. The resistive inner-layer equations in cylindrical geometry were derived by CGJ, and the extensions of these equations in general toroidal geometry were derived by GGJ. In a pressureless plasma [4], the conditions for instability can be understood as $\Delta^{{\prime}} \gt 0$. In a plasma with pressure, an instability can be possible only if $\Delta^{{\prime}} \gt \Delta_\mathrm{crit}$, where $\Delta^{{\prime}}$ is the stability parameter measuring the free energy available for the mode, and $\Delta_\mathrm{crit}$ is a positive threshold value. As a result, several numerical studies have explored the implications of the stabilizing effect of finite pressure on the resistive stability of tokamaks, resulting in a well-understood dependence of tearing mode stability on the pressure gradient $\nabla p$ at the rational surface [5].

In the absence of $\nabla p$ at the resonant surface, two approaches to solving this problem have been reviewed in Ham et al [6]. Firstly, one can deduce $\Delta^{{\prime}}$ from the tearing mode growth rate ω, calculated by a resistive MHD code using the known dispersion relation for the resistive MHD model. Alternatively, one can use a resistive MHD code to obtain a set of basis functions from which $\Delta^{{\prime}}$ can be constructed. Ham et al [6] described an artificial pressure flattening function at the rational surfaces, and its relationship has been established with the calculation of $\Delta^{{\prime}}$. Alongside, Bishop et al [7] discussed a localized pressure flattening perturbation at the resonant surface in order to assess the sensitivity of $\Delta^{{\prime}}$ to such effects. Indeed, the modelling of this in linear tearing theory by including anisotropic thermal transport in the governing equations can significantly modify the Glasser dispersion relation [8, 9], yielding the form of the natural diffusion length scale as,

$$\begin{align} {w}_{d} =2\sqrt{2}\left(\frac{\chi_{\perp}}{\chi_{\Vert}}\right)^{1/4} \left(\frac{r_{s}R_{0}}{n\,s}\right)^{1/2}. \end{align} \tag{ 1 }$$

Here, R₀ is the major radius, r_s is the radial location of a rational surface, n is the toroidal mode number, $s = r q^{{\prime}}/q$ is the positive equilibrium magnetic shear and q is the safety factor, and $\chi_{\perp}$ and $\chi_{\Vert}$ are the perpendicular and parallel transport thermal diffusivities, respectively. Ham et al [6] showed that $\delta_{L}\lt w_{d}$, where δ_L is the characteristic length scale within the resistive singular layer. For anisotropic plasmas [10] ($\chi_{\perp}\lt \chi_{\Vert}$), one could see that $\delta_{L}\lt w_{d}\lt 2\sqrt{2}(r_{s}R_{0}/n \,s)^{1/2}$. Thus, the upper bound of δ_L has no dependence on plasma resistivity.

In this article, we focus on determining pressure-induced tearing instabilities using the variational energy principle of multi-region relaxed MHD (MRxMHD). The MRxMHD variational principle is a generalization of the global Taylor's relaxation conjecture to the partial relaxation, to understand the physical mechanism of magnetic reconnections.

In MRxMHD, the whole plasma Ω is partitioned into the discrete number of Beltrami relaxed plasma regions, $\Omega_{l}$'s, such that the pressure, $p_{l} = const.$ and $\nabla p_{l} = 0$ within these regions. Each relaxed region is then bounded by freely variable toroidal interfaces, on the outer edge by the interface I_l and on the inner edge by $I_{l-1}$ that are assumed invariant ideal surfaces during the minimization of MRxMHD energy. The multi-region relaxed-MHD energy principle [11] minimizes the total potential and thermal energy under invariant topological constraints, known as magnetic helicity,

$$\begin{equation} K_l=\int_{\Omega_l} d^{3}\tau\,\textbf{A} \cdot( \nabla \times \textbf{A}), \end{equation} \tag{ 2 }$$

which takes a variational form [12]

$$\begin{align} \mu_0F = \sum_{l=1}^{N_v} \left[\int_{\Omega_l} \left(\frac{p}{\gamma-1}+\frac{B^2}{2}\right) d^{3}\tau - \frac{\mu_l}{2} ( K_l-K_{l,0}) \right], \end{align} \tag{ 3 }$$

where µ_l is defined as the Lagrange multiplier. Here, the magnetic helicity identifies as the volume-preserving invariant quantity under the gauge transformation of the vector potential, $\textbf{A}\rightarrow \textbf{A}+\nabla \Xi,$ where Ξ is a single-valued gauge potential, alongside the toroidal and poloidal magnetic fluxes. In each $\Omega_l$, the mass and entropy constraints yield an isentropic ideal-gas constraint, $p_{l} V_{l}^{\gamma} = c_{l}$, where V_l is the volume of $\Omega_l$ and c_l is a constant. The volume $\Omega_l$ enclosed by 'ideal interfaces', is constrained to have helicity $K_{l,0}$, the poloidal flux $\Delta\psi_{p,l}$, and the toroidal flux $\Delta\psi_{t,l}$. This theory unifies the ideal MHD energy principle and Taylor's relaxation conjecture [13] by allowing a less-restrictive class of variations in comparison to ideal MHD. These variations allow magnetic reconnection to form islands and chaotic fields. MRxMHD shows no explicit dependence on non-ideal dissipation parameters.

To numerically access the extremizing states of MRxMHD plasmas, the Stepped-Pressure Equilibrium code (SPEC) [12, 14] was developed. SPEC uses the pseudo Galerkin method with a Fourier–Galerkin discretization, and operates in slab, cylindrical and toroidal geometry. For fixed-boundary simulations, SPEC requires as inputs the plasma boundary and the N_v number of Taylor relaxed volumes, the enclosed poloidal $\Delta \psi_{p,l}$ and toroidal flux $\Delta \psi_{t,l}$, and magnetic helicity $K_{l,0}$ in each volume $\Omega_l$, i.e. $\{p,\Delta \psi_{p,l},\Delta \psi_{t,l},K_{l,0}\}$. Alternatively, if the helicity multiplier µ_l or parallel current is given, then the equilibrium can be described by $\{p,\Delta \psi_{p,l},\Delta \psi_{p,l},\mu_l\}$. Then, as a part of the energy minimization process, the geometries of the ideal interfaces are varied to ensure that force balance is achieved across each barrier.

The linear stability of MRxMHD equilibria in SPEC can also be assessed to analyze the magneto-hydrostatic instabilities. It is shown that the linear stability analysis of MRxMHD can reproduce both ideal and resistive MHD stability results [15–18]. General derivations of the SPEC-stability matrix in toroidal geometry are discussed in Hennenberg et al [19] and Kumar et al [18].

For MRxMHD, a particular effort has been made to clarify the relationship with the outer resistive boundary layer stability conditions, that is, $\Delta^{{\prime}}$ for the case of pressureless slab and cylindrical tokamak plasma. Using SPEC, in slab geometry, Loizu and Hudson [16] found that the variational principle of MRxMHD and the corresponding stability boundary is in exact agreement with the linear tearing mode condition $\Delta^{{\prime}}$, for $\tilde{\delta}/L\lt 0.2$, where $\tilde{\delta}$ is the arbitrarily small thickness (or width) of the resistive current sheet layer and L is the length of the current sheet along the y-direction. For finite pressure cylindrical plasma, CGJ showed that a pressure gradient within the resistive layer can drive a tearing-type instability, and GGJ later showed that in toroidal geometry, where the average curvature is favourable, tearing modes can be strongly stabilized by pressure gradient effects within the resistive layer. To our knowledge, there have not been any studies conducted on the pressure-induced tearing modes for MRxMHD plasmas. The question addressed in this article is: what happens to the MRxMHD resistive layer as a result of tearing instability with finite pressure? We investigate this mode stability using our compressible MRxMHD stability model, and clarify its applicability regimes. This goal is achieved as follows: by investigating a mode stability in cylindrical and tokamak geometry, where the role of the resistive volume layer in MRxMHD can be easily quantified.

This article is structured as follows. Section 2.1 outlines the stability conditions between the CGJ resistive layer theory and MRxMHD model, in cylindrical geometry. We extend our generalized expression for the Hessian matrix of Kumar et al [17] to account for finite compressibility $\gamma = 5/3$. Section 2.2 investigates the role of the resistive inner layer and compares the stability boundary obtained with the SPEC–Hessian and the CGJ model. Section 3.1 extends the work to toroidal geometry and describes the stability conditions between the GGJ resistive layer theory and MRxMHD model. Section 3.2 provides a comparison study of marginal stability prediction of SPEC stability with the GGJ model, and highlights the significance of the MRxMHD energy principle to predict the modified tearing mode. Finally, section 4 discusses the conclusion and identifies future work.

In the CGJ model, all the dynamics of the tearing mode are contained in the linearized set of resistive MHD equations, which can be written as,

$$\begin{align} \rho\omega^{2}\boldsymbol{\xi} = (\nabla \,\times\, \textbf{b}) \times\, \textbf{B} + \textbf{J} \, \times\, \textbf{b} + \nabla( \gamma p \nabla\cdot\boldsymbol{\xi} + \boldsymbol{\xi} \cdot \nabla p), \end{align} \tag{ 4 }$$

$$\begin{align} \textbf{b} - \eta\nabla^{2} \textbf{b} = \nabla \, \times \,(\boldsymbol{\xi} \,\times\, \textbf{B}), \end{align} \tag{ 5 }$$

where ξ and b denote the perturbed displacement and magnetic field, respectively, and η is the plasma resistivity. Fruth et al [22] showed that the approximate balance between the curvature force driving the interchange mode and restoring magnetic forces within the resistive singular layer requires that $qr^{2} \sim \eta$.

Consider the coordinates ($r,\theta,\phi$), such that the equilibria depends only on the radius r. We non-dimensionalize all quantities: scaling length to the plasma-wall boundary (such that a = 1) and the magnetic field to its axis r = 0, such that $B_{z}(0) = 1$. In cylindrical tokamak ordering [23], the CGJ described the resistive inner layer equations as,

$$\begin{align} b_{r}^{{\prime} ^{\prime}}&= Q\,({b}_{r} - r \,{\xi}_{r}), \end{align} \tag{ 6 }$$

$$\begin{align} Q^2 {\xi}_{r}^{{\prime\prime}} &= Q\,r^{2}{\xi}_{r} - D_{s}{Y} - Q_{r} \,{b}_{r}, \end{align} \tag{ 7 }$$

$$\begin{align} {Y}^{{\prime\prime}}&= \left(Q+\frac{Q}{\beta} +\frac{r^{\,2}}{Q}\right){Y} - \left(Q+\frac{Q}{\beta} - \frac{Q\,S}{D_s}\right){\xi}_{r} - \frac{r}{Q}{b_{r}}, \end{align} \tag{ 8 }$$

where the ${b_{r}}$, ξ_r and Y denote the radial component of the perturbed magnetic field, the electrostatic potential/displacement vector and the perturbed pressure along the equilibrium magnetic field, respectively. Here, $Q = \omega \delta_{L}^{2}/\eta$ and $\delta_{L} = \eta^{1/3}\rho^{1/6}\left(\frac{qr_{s}}{mq^{{\prime}} B_{\theta}}\right)^{1/3}$, defined as the characteristic resistive thickness of the inner layer, where ρ is the mass density, m is the poloidal mode number and B_θ is the azimuthal component of the equilibrium field. The other components, such as magnetic shear S, Suydam's parameter $D_{s},$ and β depend upon plasma equilibrium quantities. The solutions within the resistive layer then match with the inertia-free outer layer stability condition $\Delta^{{\prime}}$ (Furth et al [22]), using the asymptotic matching technique. For finite compressibility, the dispersion relation for the system of equations (6)–(8) is obtained from the matching condition, given as $a\,\Delta^{{\prime}} = \, \Delta_{}(Q)$, with Q being complex. For the particular choice of parameters, an equilibrium is resistive interchange unstable only if $\Delta^{{\prime}}$ exceeds a critical value, $\Delta_\mathrm{crit} \gt 0$ (see equation (9) of Ham et al [6]), in the vicinity of the resonant location. That is, the stability occurs when $\Delta^{{\prime}} \lt \Delta_\mathrm{crit}$.

To examine the stability threshold of resistive interchange mode in the vicinity of the resonant singular layer, we introduce a localized stability parameter $Z(\delta_{L}^\mathrm{CGJ})$, defined as

$$\begin{align} Z(\delta_{L}^\mathrm{CGJ})=\frac{a\Delta^{{\prime}}}{\Delta_\mathrm{crit}}. \end{align} \tag{ 9 }$$

An instability will occur if $Z(\delta_{L}^\mathrm{CGJ}) \gt 1$ for $\Delta_\mathrm{crit} \gt 0$, such that the marginal stability threshold is determined when $Z(\delta_{L}^\mathrm{CGJ}) = 1$. By doing this, we will be able to bring out the relationship between the resistive layer of CGJ and the MRxMHD model.

In this article, we restrict our stability consideration only to the vicinity of the singular surface where the rotational transform $\iota\!\!\textrm{-} = \iota/2\pi = n/m = 1/2$ is rational. Traditionally, the linear and non-linear tearing mode layer theories predict a stabilizing effect arising from local pressure gradients at the resonant surface coupled to favourable average curvature [24, 25]. In MRxMHD theory, the pressure gradient $\nabla p$ is considered to be zero in the vicinity of the resonant rational surface, to circumvent the Pfirsch–Schlüter current, which takes the form of a $1/x$ singularity. Thus, to satisfy this condition, for a given characteristic radial width of a volume $\delta_{v}^\mathrm{SPEC}$ ($\ll$the plasma minor radius, a) the resonant rational surface where $\textbf{k}\cdot \textbf{B} = 0$ must fall within the Taylor relaxed volume. Thus, the $\delta_{v}^\mathrm{SPEC}$ is understood as a user-defined parameter in SPEC, such that $0\lt \delta_{v}^\mathrm{SPEC}\ll a$. In figure 1, we show a schematic sketch of $\delta_{v}^\mathrm{SPEC}$ and $\delta_{L}^\mathrm{CGJ}$ as a function of r in the vicinity of the $q = m/n = 2/1$ rational surface (dashed grey line).

Zoom In Zoom Out Reset image size

The majority of the pressure gradients are localized on the 'ideal interfaces'. As a consequence, the 'ideal interfaces' enclosed adjacent to a resonant volume must have irrational rotational transform $\iota\!\!\textrm{-}$. This condition is also extremely crucial. If an 'ideal interface' of a rotational transform $\iota\!\!\textrm{-} \in \,\mathbf{Q}^{+}$ persists and the pressure jump is non-zero, then that surface can be unstable to localized ideal modes driven by surface currents [26].

Finite pressure jumps are allowed at the irrational surfaces, and can be interpreted as non-resonant Kolmogorov, Arnold and Moser (KAM) surfaces [27]. To obtain an irrational rotational transform on an 'ideal interface' or a surface (which satisfies the condition $\textbf{B}\cdot\textbf{n} = 0$), Greene and Mackay [28–30] provided very insightful precise methods to determine the existence of a given irrational surface (and it only really makes sense to describe an invariant surface by the rotational transform $\iota\!\!\textrm{-}$). The existence of a given irrational surface is also closely related to the stability of nearby periodic orbits. That is, when $\iota\!\!\textrm{-}$ is irrational, a single field line ergodically covers the flux surface, and the surface is referred to as an irrational surface.

We will now explain how to establish an MRxMHD equilibrium and its stability with SPEC. To find an MRxMHD equilibrium, in each $\Omega_{l}$, SPEC requires the pressure p_l , the enclosed poloidal $\Delta \psi_{p,l}$ and toroidal flux $\Delta \psi_{t,l}$, and the magnetic helicity K_l as input parameters, i.e. $\{p,\Delta \psi_{t},\Delta \psi_{p},K\}_{l = 1,2{\ldots}N_{v}}$. Therefore, at first, for the special case of the relaxed volume containing the resonant rational surface, we let the resonant volume have an arbitrary radial width, denoted by $\delta_{v}^\mathrm{SPEC}$. The input parameters K_l , $\Delta \psi_{t,l}$, $\Delta \psi_{p,l}$ within each $\Omega_{l}$, are determined by discretizing the volume-averaged magnetic helicity $\langle K_{}\rangle$ and the flux-averaged toroidal and poloidal flux ($\tilde{\psi}_{tw}$ and $\tilde{\psi}_{pw}$) profiles over the required number of volumes N_v . These averaged quantities are obtained and evaluated from the equilibrium profiles, which are considered for the simulation. Thus, the size of $\delta_{v}^\mathrm{SPEC}$ is parameterized by the K_l and the enclosed fluxes. Note that the $\langle K_{}\rangle$ is normalized to the total helicity and $\tilde{\psi}_{tw}$ varies between 0 and 1. Then, adjacent to the resonant volume, where we determine the KAM surfaces, we represent them by the functional form of a co-ordinate surface. To adapt such surfaces in SPEC cylindrical geometry, their co-ordinate geometry can be constructed as, $\textbf{x}^{c}_{l}(\theta, \zeta) = R_{l}(\theta,\zeta) \cos \theta \hat{\textbf{i}}+ R_{l}(\theta,\zeta) \sin \theta \hat{\textbf{j}} + \zeta \hat{\textbf{k}}$, where the toroidal angle, ζ, is identical to the cylindrical angle, $\zeta \sim \phi$. The $R_{l}(\theta,\zeta)$ is understood to be in the Fourier summation of $R_{l}(\theta, \zeta) = \sum_{l,i} R_{l,i} \cos(m_{i}\theta - n_{i}N_{p}\zeta)$, where m_i and n_i are $i\textrm{th}$ poloidal and toroidal harmonics and N_p is the field periodicity.

In cylindrical geometry SPEC, the stability of any MRxMHD equilibrium [17] can be assessed by considering the infinitesimal variation of the interface force balance term, $f_{l} = -[[p+B^{2}/2\mu_{0}]]_{l}\textbf{n}_{l}$, with respect to the perturbation of interface geometry $\textbf{x}^{c}_{l}$. Here, n_l denotes the unit normal to the interface $\textbf{x}^{c}_{l}$. This form of change is numerically interpreted as the Hessian matrix, which can be written as

$$\begin{align} {\mathbf H}_{j,k,l,{l^{{\prime}}}}&=\frac{\delta}{\delta \textbf{x}^{c}_{l^{{\prime}},k}} \left(\delta F/\delta \textbf{x}^{c}_{l,j}\right), \end{align} \tag{ 10 }$$

where j and k are defined as dummy variables for the Fourier harmonics for the reader's clarity, where $N_{m,n}$ is the total number of Fourier modes, and l and l^ʹ represent the different interface labels. The above equation is expanded as the Fourier summation $\partial {F}_l/\partial \textbf{x}^{c}_{l_{}} = \sum_i\partial {F}_l/\partial \textbf{x}^{c}_{l_{},i}\cos(m_i\theta-n_iN_{p}\zeta)$ and $\textbf{x}^{c}_l = \sum_i \textbf{x}^{c}_{l,i}\cos(m_i\theta-n_iN_{p}\zeta)$. When the matrix H is evaluated at fixed magnetic helicity and enclosed fluxes (toroidal and poloidal), its eigenvalues provide information about the stability corresponding to each Fourier mode harmonics $m_{i},n_{i}$'s. To take account of finite compressibility ($\gamma = \frac{5}{3}$) in SPEC, the pressure variations can be computed adiabatically, as the change in the corresponding volume V_l of relaxed plasma volumes $\Omega_{l}$, that is

$$\begin{align} \frac{\delta p_{l}}{p_{l}} = -\gamma\frac{\delta V_{l}}{V_{l}}, \end{align} \tag{ 11 }$$

where the $\delta V_{l} = (\partial V_{l}/\partial R_{l})\delta R_{l}$. An expression to compute the volume V_l , which is enclosed by the $l\textrm{th}$ and $({l-1})\textrm{th}$ interface, can be obtained by an integral form of

$$\begin{align} V_{l} &= \int_{{\Omega_{l}}} d^{3}\tau = \frac{1}{2}\int_{{\Omega_{l}}} \; \nabla \cdot \textbf{x}^{c}_{l}\,\, d^{3}\tau = \frac{1}{2} \int_{{\partial \Omega_l}} \; \textbf{x}^{c}_{l} \cdot d\textbf{S}, \end{align} \tag{ 12 }$$

$$\begin{align} &= \frac{1}{2} \int_{0}^{2\pi} d{\theta} \int_{0}^{2\pi/N_{p}} d{\zeta} (\textbf{x}^{c}_{l} \cdot {\textbf{e}_{\theta}} \times {\textbf{e}_{\zeta}}),\quad\quad\quad\, \end{align} \tag{ 13 }$$

$$\begin{align} &= \frac{1}{2} \int_{0}^{2\pi}\!\!\!d{\theta} \int_{0}^{2\pi/N_{p}} d{\zeta} \; R_{l}^{2},\quad\quad\quad\quad\quad\quad\quad\,\,\,\, \end{align} \tag{ 14 }$$

where we have used $\nabla \cdot \textbf{x}^{c}_{l} = 2$ (because it is 2D), and have assumed that the domain is periodic in the angles. The above equation is understood as a summation of the Fourier harmonics as

$$\begin{align} V_{l} = \frac{1}{2} \; \frac{4\pi^{2}}{N_{p}} \; \sum_{i}\sum_{j}\,R_{l,i}R_{l,j}\,\oint\!\!\!\oint \!\! d\theta d\zeta \;\cos\alpha_{i}\,\cos\alpha_{j}, \end{align} \tag{ 15 }$$

where $i\textrm{th}$ and $j\textrm{th}$ are the Fourier harmonics of R_l , $\alpha_{i} = m_{i}\theta - n_{i}\zeta$ and $\alpha_{j} = m_{j}\theta - n_{j}\zeta$. The required partial derivative $\frac{\partial V_{l}}{\partial R_{l,i}}$, with their trigonometrical quantities, can be obtained as

$$\begin{align} \frac{\partial V_{l}}{\partial R_{l,i}} = \frac{1}{2} \; \frac{4\pi^{2}}{N_{p}} \; \sum_{i}\sum_{j}\,R_{l,j}\,[2\cos(\alpha_{i}-\alpha_{j})+\,2\cos(\alpha_{i}+\alpha_{j})] .\end{align} \tag{ 16 }$$

The symmetry of H, means that all its eigenvalues are real numbers [31]. Using the principle axis theorem [32], the quadratic form $\delta \textbf{{x}}^T \cdot \textbf{H} \cdot \delta \textbf{{x}}$ can be condensed as

$$\begin{align} \delta {\textbf{x}}^T \cdot \textbf{H} \cdot \delta {\textbf{x}}= \sum_{j}\lambda_{j}v_{j}^{2}, \end{align} \tag{ 17 }$$

where the λ_j is the eigenvalue of H and v_j is the corresponding eigenvector for $j = (1,2,3,\ldots,N_{mn}$). The stability of an equilibrium can be predicted from the sign of the eigenvalue λ_j : that is, if a j exists such that $\lambda_j \lt 0$, then an equilibrium is said to be unstable and, if all the $\lambda_j \gt 0$, then an equilibrium is said to be stable. These eigenvalues hereafter regarded as $\lambda^{m,n}_\mathrm{SPEC}$, from H, are evaluated numerically using the SPEC–Hessian calculation.

To compute the stability conditions of CGJ with SPEC, we conform to their notation and express the smallest negative eigenvalue (normalized to its maximum value), referred to as ${\textrm{min}}\{\lambda^{m,n}_\mathrm{SPEC}\}$ in terms of the stability parameter $Z(\delta_{v}^\mathrm{SPEC})$, defined by,

$$\begin{align} Z(\delta_{v}^\mathrm{SPEC})= 1- {\textrm{min}}\{\lambda^{m,n}_\mathrm{SPEC}(\delta_{v}^\mathrm{SPEC})\}. \end{align} \tag{ 18 }$$

Therefore, the stability conditions of $Z(\delta_{v}^\mathrm{SPEC})$ can be interpreted as, if $\lambda^{m,n}_\mathrm{SPEC}(\delta_{v}^\mathrm{SPEC}) \lt 0$ exists, then an instability occurs when $Z(\delta_{v}^\mathrm{SPEC}) \gt 1$, and the stability is determined for $Z(\delta_{v}^\mathrm{SPEC}) \lt 1$.

To clarify the concepts described in the previous section, we discuss the resistive interchange instability of a cylindrical equilibrium considered in Izzo et al [33]. The equilibrium of interest is described by the pressure profile

$$\begin{align} p(r) = p_{0}(0.001+0.028\,r^2 -0.059\, r^4 + 0.03\, r^6), \end{align} \tag{ 19 }$$

and the safety factor profile as

$$\begin{align} q(r) = q_{0}(1+(r/0.8)^2), \end{align} \tag{ 20 }$$

where $q_{0} = 1.6$, the aspect ratio A = 5 and $0 \leqslant p_{0} \leqslant1$. For the mode perturbation, m = 2, n = 1 this equilibrium is always tearing unstable. The physical motivation for using the equilibrium with the pressure gradient reversed at the resonant layer is to simulate the effects of good average curvature. Izzo et al [33] investigated what happens to the resistive tearing mode at a fixed resistivity as the pressure parameter p₀ is gradually increased.

We investigate the role of the resistive singular layer in both models by considering the equilibrium configuration in the scenario of fixed $p_{0} = 0.25$ (see figure 2). In figure 3(a) we have plotted $Z(\delta_{L}^\mathrm{CGJ})$ as a function of $\hat{\delta}$. We see that as $\delta^\mathrm{CGJ}_{L}$ increases, $Z(\delta^\mathrm{CGJ}_{L})$ crosses the stability threshold line (Z = 1) and predicts instability for the mode perturbation $m/n = 2/1$. Here, the stability threshold in terms of $\delta_{L}^\mathrm{CGJ}$ is approximated as $1.77\times 10^{-3}$. Now, the eigenvalues $\lambda_\mathrm{SPEC}^{2,1}$ from equation (10) are evaluated numerically using the SPEC for different values of $\delta_{v}^\mathrm{SPEC}$, which is understood in terms of equation (18). We observed that, for the sufficiently small value of $\delta_{v}^\mathrm{SPEC}$, $Z(\delta_{v}^\mathrm{SPEC})$ crosses its the stability threshold, and coincides with $Z(\delta_{L}^\mathrm{CGJ}) \gt 1$ (see figure 3(a)). Since the $m/n = 2/1$ mode is unstable, this indicates that the effect of the width of the singular layer in both models is comparable, close to the marginal stability locus. Finally, figure 3(b) shows the radial structure of the SPEC eigenfunction $\boldsymbol{\xi}\cdot\nabla s$ for the $m/n = 2/1$ unstable equilibrium case with $\delta_{v}^\mathrm{SPEC}\sim1.87 \times 10^{-3}$. In the vicinity of the q = 2 rational surface, typical spatial behaviour of this tearing eigenfunction can be observed. Here, the $\boldsymbol{\xi}\cdot\nabla s$ is defined as the radial perturbed component of the interface displacement. These results show that the marginal stability threshold of CGJ and MRxMHD theory coincide, when the $\delta_{v}^\mathrm{SPEC}$ is proportional to the $\delta^\mathrm{CGJ}_{L}$, i.e. $\delta_{v}^\mathrm{SPEC}\sim \delta_{L}^\mathrm{GGJ}$.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Moreover, we postulate that as $\delta_{v}^\mathrm{SPEC}$ decreases, the 'ideal interfaces' surrounding the $\iota\!\!\textrm{-} = 1/2$ resonant surface came sufficiently close to rational surfaces, and induce shielding currents. That is, in the limit of the vanishing width of the resonant volume $\delta_{v}^\mathrm{SPEC}\rightarrow{} 0$, the parallel current density becomes infinite, such that the parallel current within the volume region becomes finite and non-zero. In accordance, the emergence of shielding currents from the 'ideal interfaces' could also be a potential reason for stabilization of this resistive mode. It is the current sheet that allows the small solutions on either side of a singular surface to be disconnected in ideal MHD, screening one side from the other. In ideal MHD, field-line reconnection is forbidden by the frozen-in flux conditions; therefore, the current sheet must form to prevent the tearing mode island (of arbitrarily small amplitude in the linearized approximation) that forms in relaxed MHD. This phenomenon has also been observed in the formation of current sheets in magnetic reconnection [34, 35].

In the vicinity of a resonant volume, SPEC allows a transition from partial Taylor relaxation to ideal MHD as $\delta_{v}^\mathrm{SPEC}\rightarrow{}0$. Here, it should be noted that the 'ideal interfaces' will not overlap, even when SPEC computes an equilibrium solution as $\delta_{v}^\mathrm{SPEC}$ decreases in the vicinity of the resonant surface. Overlapping of the 'ideal interfaces' is not allowed on both conceptual and computational grounds.

In the GGJ model [3], while an equilibrium configuration is ideally stable for $q_{0}\gt1$ and satisfies the resistive interchange stability condition $D_{R}\lt 0$, it is unstable to a special case of tearing mode known as the modified tearing mode if $\Delta^{{\prime}} \gt \Delta_\mathrm{crit}$ in the vicinity of the resonant location. Here, we define

$$\begin{align} D_{R} &= \frac{2q^{4} p^{{\prime}}}{rB_{0}^{2} q^{{\prime} 2}} \left(\!\!1-\frac{1}{q^2} +\frac{qq^{{\prime}}}{r^3}\right) \!\!\int_{0}^{r} \!\!dt\!\! \left(\frac{t^{3}}{q^2} - \frac{2 R_{0}t^{2}p^{{\prime}}}{B_{0}^{2}}\right)\!, \end{align} \tag{ 21 }$$

$$\begin{align} \Delta_\mathrm{crit}&\approx 1.54 \left(\frac{V_{S}}{X_{0}}\right)|D_{R}|^{5/6}, \end{align} \tag{ 22 }$$

and $\Delta^{{\prime}}$ is the generalized jump in the logarithmic derivative of the perturbed magnetic flux across the resistive layer. Here, $V_{s}/X_{0}$ is defined as the ratio of the macroscopic resistive scale length to a MHD scale that varies as $\eta^{1/3}$.

To examine the stability of a modified tearing mode in the vicinity of a singular layer, where the magnetic winding number is rational, we introduce a localized stability parameter $Z(\delta^\mathrm{GGJ}_{L},\upsilon)$, which is a function of both $\delta^\mathrm{GGJ}_{L}$ and a magnetic shear parameter $\upsilon \in \mathbf{Z}^{+}$. We define the $Z(\delta^\mathrm{GGJ}_{L},\upsilon)$ as,

$$\begin{align} Z(\delta^\mathrm{GGJ}_{L},\upsilon) = \frac{\Delta^{{\prime}}}{\Delta_\mathrm{crit}}, \end{align} \tag{ 23 }$$

where $\Delta_\mathrm{crit}$ is rewritten in terms of the resistive width of a singular layer $\delta^\mathrm{GGJ}_{L}$ as,

$$\begin{equation} \Delta_\mathrm{crit}\approx \frac{1.54\,|D_{R}|^{5/6}|Q|^{1/4}}{\delta^\mathrm{GGJ}_{L}}, \end{equation} \tag{ 24 }$$

with Q being a complex and dimensionless variable. It is obtained by solving a dispersion relation when the boundary layer solutions of the inner layer are matched with the outer layer. This dispersion relation is written as

$$\begin{equation} \Delta(Q)=\Delta^{{\prime}}, \end{equation} \tag{ 25 }$$

where

$$\begin{align} \Delta^{{\prime}}(Q)&=\frac{\pi D_{R}(2V_{s}/X_{0})^{1-2H} \Gamma(1/4)\Gamma^{2}(1-{H}/{4})}{(1-2H)\Gamma(1-{H}/{2})}\\ &\quad \times \frac{\Gamma(3/4-H/2)Q^{{(2H-1)}/{4}}}{\left(\cos({\pi H}/{2})\Gamma({1+H}/{2})\Gamma(1-H)\right)^{2}}\nonumber\\ &\quad \times \left( \frac{Q^{3/2}}{D_{R}}-\frac{\,\Gamma(3/4)\Gamma^{2}(1/2-H/4)\Gamma(1/4-H/2)}{4\Gamma(1/4)\Gamma^{2}(1-H/4)\Gamma(3/4-H/2)}\right),\nonumber \end{align} \tag{ 26 }$$

where Γ is the gamma function and the analytical expressions for H, $V_{s}/X_{0}$ and $\delta^\mathrm{GGJ}_{L}(r)$ are written as

$$\begin{align} H &= - \frac{2q^{5} p^{{\prime}}}{r^{\,4}B_{0}^{2} q^{{\prime} }} \!\!\int_{0}^{r} \!\!dt\! \left(\frac{t^{3}}{q^2} - \frac{2 R_{0}t^{2}p^{{\prime}}}{B_{0}^{2}}\right), \end{align} \tag{ 27 }$$

$$\begin{align} \frac{V_{S}}{X_{0}}&=\left(\frac{n B_{0}}{\eta R}\frac{q^{{\prime}}(r)}{q(r)} \right)^{1/3} (1+2q(r)^{2})^{-1/6}\rho^{-1/6}, \end{align} \tag{ 28 }$$

$$\begin{align} \delta^\mathrm{GGJ}_{L} &= \left(\left(\frac{\eta R}{n B_{0 }}\frac{q(r)}{q^{{\prime}}(r)}\right)^{2}\! (1+2q(r)^{2})\right)^{1/6}\!\!|Q|^{1/4}\rho^{1/6}. \end{align} \tag{ 29 }$$

The above equations (27)–(29) are available from Glasser et al [3] (see equations (A21), (A31) and (A32)).

An instability occurs (either pressure-induced or modified tearing mode), if a solution of equation (25) with $Re(Q) \gt 0$ exists, for a given value of $\Delta^{{\prime}}$. The modified tearing mode is unstable if the $Z(\delta^\mathrm{GGJ}_{L},\upsilon)\gt1$ condition is satisfied, otherwise stable below the critical condition of $Z(\delta^\mathrm{GGJ}_{L},\upsilon) = 1$, which denotes the marginal stability conditions for this kind of tearing mode.

From the perspective of the MRxMHD model, we now proceed to characterize the axisymmetric equilibria and its stability. In accordance with the previous section 2.1, we first let the resonant volume have an arbitrary radial width, denoted by $\delta_{v}^\mathrm{SPEC}$. Then, the input parameters, such as K_l , $\Delta \psi_{t,l}$, $\Delta \psi_{p,l}$, are determined by discretizing the volume-averaged magnetic helicity $\langle K_{} \rangle$ and the flux-averaged toroidal and poloidal flux ($\tilde{\psi}_{tw}$ and $\tilde{\psi}_{pw}$) profiles over the required number of volumes N_v . These averaged quantities are evaluated from the given equilibrium profiles, which are considered for the simulations. Thus, the size of the $\delta_{v}^\mathrm{SPEC}$ is also parameterized by the K_l and the enclosed fluxes. Here, $\langle K_{}\rangle$ is also normalized to the total helicity and the $\tilde{\psi}_{w}$ varies between 0 and 1. Then, adjacent to a resonant volume, where we determine the KAM surfaces, SPEC's toroidal co-ordinate is constructed as $\textbf{x}_{l}(\theta,\zeta) = R_{l}(\theta,\zeta)\hat{\textbf{e}}_{R}+Z_{l}(\theta,\zeta)\hat{\textbf{e}}_{Z}$. Here, $\hat{\textbf{e}}_{R} = \cos\phi \hat{\textbf{i}} +\sin\phi \hat{\textbf{j}}$ for the toroidal angle $\zeta\sim\phi$, and the $R_{l}(\theta,\zeta)$ and $Z_l(\theta,\zeta)$ are an even and odd function of $(\theta,\zeta)$, respectively. The symmetric and non-symmetric variables are discretized in the Fourier basis function as $R_l(\theta,\zeta) = \sum_{i} R_{l,i} \cos(m_{i}\theta-n_{i} N_{p}\zeta)$ and $Z_l(\theta,\zeta) = \sum_{i} Z_{l,i} \sin(m_{i}\theta-n_{i} N_{p}\zeta)$.

In toroidal geometry, the stability of an MRxMHD equilibrium can be assessed by considering the infinitesimal variation of the interface force balance term, $f_{l} = -[[p+B^{2}/2\mu_{0}]]_{l}\textbf{n}_{l}$, with respect to the poloidal and toroidal perturbation of interface geometry x_l [18]. Similar to the previous cylindrical stability implementation in SPEC, this form of change is also numerically interpreted as the Hessian matrix, which can be written as [18]

$$\begin{align} {\mathbf H}_{j,k,l,{l^{{\prime}}}}&=\frac{\delta}{\delta \textbf{x}_{l^{{\prime}},k}} \left(\delta F/\delta \textbf{x}_{l,j}\right), \end{align} \tag{ 30 }$$

where j and k are defined as dummy variables for the Fourier harmonics for clarity, where $N_{m,n}$ is the total number of Fourier modes, and l and l^ʹ represent the different interface labels. To include the effects of finite compressibility in toroidal geometry, the pressure variation can be computed using equation (11), where $\delta V_{l} = (\partial V_{l}/\partial R_{l})\delta R_{l}+(\partial V_{l}/\partial Z_{l})\delta Z_{l}$. The expression for V_l can be obtained by the integral

$$\begin{align} V_{l} &= \int_{{\Omega_l}} d^{3}\tau = \frac{1}{3}\int_{{\Omega_l}} \; \nabla \cdot \textbf{x}_{l}\,\, d^{3}\tau = \frac{1}{3} \int_{{\delta \Omega_l}} \; \textbf{x}_{l} \cdot d\textbf{S}, \end{align} \tag{ 31 }$$

$$\begin{align} &= \frac{1}{3} \int_{0}^{2\pi} d{\theta} \int_{0}^{2\pi/N} d{\zeta} (\textbf{x}_{l} \cdot {\textbf{e}_{\theta}} \times {\textbf{e}_{\zeta}}),\quad\quad\quad\;\;\;\! \end{align} \tag{ 32 }$$

$$\begin{align} &= \frac{1}{3} \int_{0}^{2\pi}\!\!\!d{\theta} \int_{0}^{2\pi/N} d{\zeta} \; R_{l} \left( Z_{l} R_{l,\theta} - R_{l} Z_{l,\theta} \right),\quad\, \end{align} \tag{ 33 }$$

where we have considered $\nabla \cdot \textbf{x}_l = 3$. Upon expansion of equation (33) as a summation of the Fourier harmonics, we have

$$\begin{align} \begin{aligned} V_{l} & = \frac{1}{3} \; \sum_i \sum_j \sum_k R_{l,i} \left(Z_{l,j} R_{l,k} - R_{l,j} Z_{l,k} \right) (+m_k)\\ & \quad \times \oint \!\!\! \oint \!\! \; \cos\alpha_i \cos\alpha_j \cos\alpha_k \,d\theta d\zeta, \end{aligned} \end{align} \tag{ 34 }$$

where $i\textrm{th}$, $j\textrm{th}$ and $k\textrm{th}$ are the Fourier harmonics of $R_{l},Z_{l}$. Then, the partial derivatives $\frac{\partial V_{l}}{\partial R_{l,i}}$ and $\frac{\partial V_{l}}{\partial Z_{l,i}}$ are obtained as

$$\begin{align} \frac{\partial V_{l}}{\partial R_{l,i}} & = \left(Z_{l,j} R_{l,k} m_k - R_{l,j} Z_{l,k} m_k - R_{l,j} Z_{l,k} m_k \right)\\ &\quad \times \oint\!\!\!\oint \!\! d\theta d\zeta \; \cos\alpha_i \cos\alpha_j \cos\alpha_k \nonumber\\ & \quad + \left( - Z_{l,j} R_{l,k} m_k + R_{l,j} Z_{l,k} m_k + R_{l,j} Z_{l,k} m_k \right)\nonumber \\ & \quad \times \oint \!\!\! \oint \!\! d\theta d\zeta \; \cos\alpha_i \sin\alpha_j \sin\alpha_k, \nonumber \end{align} \tag{ 35 }$$

and

$$\begin{align} \frac{\partial V_{l}}{\partial Z_{l,i}} & = \left( - R_{l,k} R_{l,j} m_i \right) \oint \!\!\! \oint \!\! d\theta d\zeta \; \cos\alpha_i \cos\alpha_j \cos\alpha_k\\ & \quad + \left( - R_{l,k} R_{l,j} m_{k} \right) \oint \!\!\! \oint \!\! d\theta d\zeta \; \cos\alpha_i \sin\alpha_j \sin\alpha_k.\nonumber \end{align} \tag{ 36 }$$

When this matrix H is evaluated at fixed magnetic helicity and enclosed fluxes, its eigenvalues provide information about the stability corresponding to each Fourier mode harmonic $m,n$. These eigenvalues, $\lambda^{m,n}_\mathrm{SPEC}$ from H, are evaluated numerically using the SPEC–Hessian calculation.

Following equation (18) of section 2.2, we now express the smallest negative eigenvalue (normalized to its maximum value), referred to as ${\textrm{min}}\{\lambda^{m,n}_\mathrm{SPEC}\}$ in terms of the stability parameter $Z(\delta_{v}^\mathrm{SPEC},\upsilon)$ as,

$$\begin{align} Z(\delta_{v}^\mathrm{SPEC},\upsilon)= 1- {\textrm{min}}\{\lambda^{m,n}_\mathrm{SPEC}(\delta_{v}^\mathrm{SPEC},\upsilon)\}, \end{align} \tag{ 37 }$$

such that the stability conditions of $Z(\delta^\mathrm{GGJ}_{L},\upsilon)$ can be interpreted the same as described before.

In this section, we discuss the modified tearing instability of a model circular tokamak (large aspect ratio) equilibrium considered in Glasser et al [3]. The aspect ratio A = 8.4, where the major radius $R_{0} = 8.4$ m and minor radius a = 1 m, such that all the equilibrium scalars are independent of the toroidal angle φ about the axis of symmetry. The equilibrium toroidal current density and the parabolic pressure profiles are described as a function of r, which is

$$\begin{align} J_{\phi}(r)=J_{0}/(1 + \upsilon \,r^{2}/a^{2})^{2},\,\,p(r)=p_{0}(1-r^{2}/a^{2}), \end{align} \tag{ 38 }$$

where $\upsilon \in \mathbf{Z}^{+}$ is the shear parameter, $J_{0} = 2B_{0}/q_{0}R_{0}$ and $p_{0} = \beta_{p}(B_{0}a/q_{0}R_{0}(1+\upsilon))^{2}$. The analytical expression for poloidal plasma beta is obtained as $\beta_{p} = \left(\frac{2Rq(a)}{B_{0}a^{2}}\right)^{2}\int_{0}^{a} rp(r) \,dr$. We investigate this equilibrium model in the scenario with $q_{0} = 1.1$, $\beta_{p} = 0.8$ and the shear parameter values υ = 2 and υ = 3. This shear parameter υ plays a critical role in the destabilizing and stabilizing factor of this plasma configuration.

Figure 4(a) compares the results of $Z(\delta_{v}^\mathrm{SPEC},\upsilon)$ and $Z(\delta^\mathrm{GGJ}_{L},\upsilon)$ as a function of $\delta_{v}^\mathrm{SPEC}$ and $\delta^\mathrm{GGJ}_{L}$, respectively. The eigenvalues $\lambda_\mathrm{SPEC}^{2,1}$ from equation (30) are evaluated numerically using the SPEC for different values of the $\delta_{v}^\mathrm{SPEC}$, which is understood in terms of equation (37). For υ = 2, we observe that the $Z(\delta_{v}^\mathrm{SPEC},2)$ predicts instability as it crosses the marginal stability threshold line (Z = 1), for sufficiently small value of $\delta_{v}^\mathrm{SPEC}$, which is approximated as $2.8 \times 10^{-3}$. In addition, the value of $\delta_{L}^\mathrm{GGJ}$ at which the $Z(\delta^\mathrm{GGJ}_{L},2)$ crosses its stability threshold coincides with that at which $Z(\delta_{v}^\mathrm{SPEC},2)$ is greater than 1. This confirms the potential relationship $\delta_{v}^\mathrm{SPEC} \sim \delta_{L}^\mathrm{GGJ}$. Now, for υ = 3, it is observed that $Z(\delta_{v}^\mathrm{SPEC},3)$ predicts instability for the smaller value of $\delta_{v}^\mathrm{SPEC}$ than for υ = 2. This is because as υ increases, the current channel shrinks in the vicinity of the rational surface, and both $\delta_{v}^\mathrm{SPEC}$ and $\delta_{L}^\mathrm{GGJ}$ reduce. For υ = 3, the threshold $\delta_{v}^\mathrm{SPEC}$ is approximated as $1.72 \times 10^{-3}$, and similar threshold behaviour is found for $Z(\delta_{L}^\mathrm{GGJ},3)$. Thus, the MRxMHD stability boundary is in agreement with the linear modified tearing mode theory. However, if $\delta_{v}^\mathrm{SPEC}$ becomes sufficiently large compared to $\delta^\mathrm{GGJ}_{L}$, it can be conceptualized that the pressure flattening in SPEC can indeed remove the stabilizing effects and considerably affect the stability boundary of the mode. For $p_{l} = const.$ over a larger volume width, the mode can still be strongly destabilized in SPEC and finds a different stability threshold.

Zoom In Zoom Out Reset image size

Finally, figure 4(b) shows the spatial structure of the SPEC eigenfunction $\boldsymbol{\xi}\cdot \nabla s$ as a function of the effective radius $r_\mathrm{eff.}$. The corresponding $m/n = 2/1$ unstable equilibrium for ν = 2 and 3 are considered with $\delta_{v}^\mathrm{SPEC}$ equal to $3.5\times 10^{-3}$ and $2.5 \times 10^{-3}$, respectively. We would like to remark that the SPEC-stability results shown in figure 4 are converged in the sense that the Fourier resolution and the radial basis function are increased. Here, $N_{v} = 180$ is considered.