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

We show that the variational energy principle of the multi-region relaxed magnetohydrodynamic (MRxMHD) model can be used to predict finite-pressure linear tearing instabilities. In this model, the plasma volume is sliced into sub-volumes separated by ‘ideal interfaces’, and in each volume the magnetic field relaxes to a Taylor state, where the pressure gradient ∇p=0 . The MRxMHD model is implemented in the Stepped-Pressure Equilibrium Code (SPEC) so that the equilibrium solution in each region is computed while preserving the force balance across the interfaces. As SPEC computes the Hessian matrix (a discretized stability matrix), the stability of an MRxMHD equilibrium can also be computed with SPEC. In this article, using SPEC, we investigate the effect of local pressure gradients and the ∇p=0 in the vicinity of the resonant surface of a tearing mode. For low-beta plasma, we have been able to illustrate a relationship between the resistive singular-layer theory (Coppi et al 1966 Nucl. Fusion 6 101; Glasser et al 1975 Phys. Fluids 18 875–88) and the MRxMHD model. Within the singular layer, the volume-averaged magnetic helicity and the flux-averaged toroidal flux are shown to be the invariants for the linear tearing modes in SPEC simulations. Our technique to compute MRxMHD stability is first tested numerically in a cylindrical tokamak and its application in toroidal geometry is demonstrated. We demonstrate an agreement between the stability boundary obtained with SPEC simulation and the resistive inner-layer theories.


Theoretical review
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 ∆ ′ > 0. In a plasma with pressure, an instability can be possible only if ∆ ′ > ∆ crit , where ∆ ′ is the stability parameter measuring the free energy available for the mode, and ∆ 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 wellunderstood dependence of tearing mode stability on the pressure gradient ∇p at the rational surface [5].
In the absence of ∇p at the resonant surface, two approaches to solving this problem have been reviewed in Ham et al [6]. Firstly, one can deduce ∆ ′ 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 ∆ ′ 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 ∆ ′ . Alongside, Bishop et al [7] discussed a localized pressure flattening perturbation at the resonant surface in order to assess the sensitivity of ∆ ′ 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, Here, R 0 is the major radius, r s is the radial location of a rational surface, n is the toroidal mode number, s = rq ′ /q is the positive equilibrium magnetic shear and q is the safety factor, and χ ⊥ and χ ∥ are the perpendicular and parallel transport thermal diffusivities, respectively. Ham et al [6] showed that δ L < w d , where δ L is the characteristic length scale within the resistive singular layer. For anisotropic plasmas [10] (χ ⊥ < χ ∥ ), one could see that δ L < w d < 2 √ 2(r s R 0 /n s) 1/2 . Thus, the upper bound of δ L has no dependence on plasma resistivity.

Multi-region relaxed MHD model
In this article, we focus on determining pressure-induced tearing instabilities using the variational energy principle of multiregion 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, Ω l 's, such that the pressure, p l = const. and ∇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, which takes a variational form [12] 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, A → A + ∇Ξ, where Ξ is a single-valued gauge potential, alongside the toroidal and poloidal magnetic fluxes. In each Ω l , the mass and entropy constraints yield an isentropic ideal-gas constraint, p l V γ l = c l , where V l is the volume of Ω l and c l is a constant. The volume Ω l enclosed by 'ideal interfaces', is constrained to have helicity K l,0 , the poloidal flux ∆ψ p,l , and the toroidal flux ∆ψ 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 fixedboundary simulations, SPEC requires as inputs the plasma boundary and the N v number of Taylor relaxed volumes, the enclosed poloidal ∆ψ p,l and toroidal flux ∆ψ t,l , and magnetic helicity K l,0 in each volume Ω l , i.e. {p, ∆ψ p,l , ∆ψ t,l , K l,0 }. Alternatively, if the helicity multiplier µ l or parallel current is given, then the equilibrium can be described by {p, ∆ψ p,l , ∆ψ p,l , µ 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][16][17][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, ∆ ′ 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 ∆ ′ , forδ/L < 0.2, whereδ 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 γ = 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.

Resistive interchange mode in a cylindrical tokamak
In this section, we examine the variational energy principle of MRxMHD for resistive interchange modes in a cylindrical tokamak, and compare it with the CGJ compressible resistive layer model. Much of the physical picture underlying this instability has been well known [20,21]. The major result of this section lies with the establishment of a clear relationship between the resistive inner layer of the CGJ and MRxMHD model, in the cylindrical geometry.

Stability conditions between CGJ and the MRxMHD model
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, 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 ∼ η. Consider the coordinates (r, θ, ϕ), 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, 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, , 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 ∆ ′ (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 ∆ ′ = ∆(Q), with Q being complex. For the particular choice of parameters, an equilibrium is resistive interchange unstable only if ∆ ′ exceeds a critical value, ∆ crit > 0 (see equation (9) of Ham et al [6]), in the vicinity of the resonant location. That is, the stability occurs when ∆ ′ < ∆ 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(δ CGJ L ), defined as An instability will occur if Z(δ CGJ L ) > 1 for ∆ crit > 0, such that the marginal stability threshold is determined when Z(δ CGJ L ) = 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 ι-= ι/2π = 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 ∇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 δ SPEC v (≪the plasma minor radius, a) the resonant rational surface where k · B = 0 must fall within the Taylor relaxed volume. Thus, the δ SPEC v is understood as a user-defined parameter in SPEC, such that 0 < δ SPEC v ≪ a. In figure 1, we show a schematic sketch of δ SPEC v and δ CGJ L as a function of r in the vicinity of the q = m/n = 2/1 rational surface (dashed grey line).
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 ι-. This condition is also extremely crucial. If an 'ideal interface' of a rotational transform ι-∈ 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 B · n = 0), Greene and Mackay [28][29][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 ι-). The existence of a given irrational surface is also closely related to the stability of nearby periodic orbits. That is, when ι-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 Ω l , SPEC requires the pressure p l , the enclosed poloidal ∆ψ p,l and toroidal flux ∆ψ t,l , and the magnetic helicity K l as input parameters, i.e. {p, ∆ψ t , ∆ψ p , K} l=1,2...Nv . 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 δ SPEC v . The input parameters K l , ∆ψ t,l , ∆ψ p,l within each Ω l , are determined by discretizing the volume-averaged magnetic helicity ⟨K⟩ and the flux-averaged toroidal and poloidal flux (ψ tw andψ 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 δ SPEC v is parameterized by the K l and the enclosed fluxes. Note that the ⟨K⟩ is normalized to the total helicity andψ 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, where m i and n i are ith 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, n l , with respect to the perturbation of interface geometry x c l . Here, n l denotes the unit normal to the interface x c l . This form of change is numerically interpreted as the Hessian matrix, which can be written as 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 ∂F l /∂x c l = 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 (γ = 5 3 ) in SPEC, the pressure variations can be computed adiabatically, as the change in the corresponding volume V l of relaxed plasma volumes Ω l , that is where the δV l = (∂V l /∂R l )δR l . An expression to compute the volume V l , which is enclosed by the lth and (l − 1)th interface, can be obtained by an integral form of where we have used ∇ · 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 where ith and jth are the Fourier harmonics of R l , α i = m i θ − n i ζ and α j = m j θ − n j ζ. The required partial derivative ∂V l ∂R l,i , with their trigonometrical quantities, can be obtained as The symmetry of H, means that all its eigenvalues are real numbers [31]. Using the principle axis theorem [32], the quadratic form δx T · H · δx can be condensed as where the λ j is the eigenvalue of H and v j is the corresponding eigenvector for j = (1, 2, 3, . . . , 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 λ j < 0, then an equilibrium is said to be unstable and, if all the λ j > 0, then an equilibrium is said to be stable. These eigenvalues hereafter regarded as λ m,n 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 min{λ m,n SPEC } in terms of the stability parameter Z(δ SPEC v ), defined by, Therefore, the stability conditions of Z(δ SPEC v ) can be interpreted as, if λ m,n SPEC (δ SPEC v ) < 0 exists, then an instability occurs when Z(δ SPEC v ) > 1, and the stability is determined for Z(δ SPEC v ) < 1.

Numerical observation of marginal stability threshold
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 p(r) = p 0 (0.001 + 0.028 r 2 − 0.059 r 4 + 0.03 r 6 ), (19) and the safety factor profile as where q 0 = 1.6, the aspect ratio A = 5 and 0 ⩽ p 0 ⩽ 1. 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 0 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(δ CGJ L ) as a function ofδ. We see that as δ CGJ L increases, Z(δ 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 δ CGJ L is approximated as 1.77 × 10 −3 . Now, the eigenvalues λ 2,1 SPEC from equation (10) are evaluated numerically using the SPEC for different values of δ SPEC v , which is understood in terms of equation (18). We observed that, for the sufficiently small value of δ SPEC v , Z(δ SPEC v ) crosses its the stability threshold, and coincides with Z(δ CGJ L ) > 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 ξ · ∇s for the m/n = 2/1 unstable equilibrium case with δ SPEC v ∼ 1.87 × 10 −3 . In the vicinity of the q = 2 rational surface, typical spatial behaviour of this tearing eigenfunction can be observed. Here, the ξ · ∇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 δ SPEC v is proportional to the δ CGJ L , i.e. δ SPEC v ∼ δ GGJ L . Moreover, we postulate that as δ SPEC v decreases, the 'ideal interfaces' surrounding the ι-= 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 δ SPEC v → 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 δ SPEC v → 0. Here, it should be noted that the 'ideal interfaces' will not overlap, even when SPEC computes an equilibrium solution as δ SPEC v decreases in the vicinity of the resonant surface. Overlapping of the 'ideal interfaces' is not allowed on both conceptual and computational grounds.  (20). Here, the number of the volume, Nv = 180, is considered in SPEC.

Modified tearing mode in large aspect ratio axisymmetric plasma
The preceding section showed how the CGJ and MRxMHD model are related in cylindrical geometry. This section compares the GGJ compressible resistive layer model to the variational energy principle of the MRxMHD for the modified tearing modes in a large aspect ratio toroidal geometry.
The modified tearing instability normally occurs during the course of tokamak discharges due to thermal instability in ohmic plasmas [36]. Usually, the stabilizing effect of magnetic shear on ideal interchange instability is eliminated by these modified tearing modes. GGJ were the first to investigate the tearing stability threshold in relation to shear, pressure and toroidicity.
Our primary goal in this section is to investigate the presence of the resistive volume layer width within MRxMHD, in comparison to the classic GGJ model. In GGJ, the resistive width of a singular layer (here, we denoted as δ GGJ L ) depends on the plasma inertia (growth rate of an unstable ideal MHD mode) and the plasma resistivity. It approaches zero as the growth rate or the resistivity reduces to zero simultaneously.

Stability conditions between GGJ and the MRxMHD model
In the GGJ model [3], while an equilibrium configuration is ideally stable for q 0 > 1 and satisfies the resistive interchange stability condition D R < 0, it is unstable to a special case of tearing mode known as the modified tearing mode if ∆ ′ > ∆ crit in the vicinity of the resonant location. Here, we define ∆ crit ≈ 1.54 and ∆ ′ 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 η 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(δ GGJ L , υ), which is a function of both δ GGJ L and a magnetic shear parameter υ ∈ Z + . We define the Z(δ GGJ L , υ) as, where ∆ crit is rewritten in terms of the resistive width of a singular layer δ GGJ L as, 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 where where Γ is the gamma function and the analytical expressions for H, V s /X 0 and δ GGJ L (r) are written as 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) > 0 exists, for a given value of ∆ ′ . The modified tearing mode is unstable if the Z(δ GGJ L , υ) > 1 condition is satisfied, otherwise stable below the critical condition of Z(δ GGJ L , υ) = 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 δ SPEC v . Then, the input parameters, such as K l , ∆ψ t,l , ∆ψ p,l , are determined by discretizing the volume-averaged magnetic helicity ⟨K⟩ and the flux-averaged toroidal and poloidal flux (ψ tw andψ 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 δ SPEC v is also parameterized by the K l and the enclosed fluxes. Here, ⟨K⟩ is also normalized to the total helicity and theψ 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 x l (θ, ζ) = R l (θ, ζ)ê R + Z l (θ, ζ)ê Z . Here,ê R = cos ϕî + sin ϕĵ for the toroidal angle ζ ∼ ϕ, and the R l (θ, ζ) and Z l (θ, ζ) are an even and odd function of (θ, ζ), respectively. The symmetric and non-symmetric variables are discretized in the Fourier basis function as R l (θ, ζ) = ∑ i R l,i cos(m i θ − n i N p ζ) and Z l (θ, ζ) = ∑ i Z l,i sin(m i θ − n i N p ζ).
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µ 0 ]] l 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] 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 δV l = (∂V l /∂R l )δR l + (∂V l /∂Z l )δZ l . The expression for V l can be obtained by the integral where we have considered ∇ · x l = 3. Upon expansion of equation (33) as a summation of the Fourier harmonics, we have where ith, jth and kth are the Fourier harmonics of R l , Z l . Then, the partial derivatives ∂V l ∂R l,i and ∂V l ∂Z l,i are obtained as ×˛˛dθdζ cos α i cos α j cos α k ×˛˛dθdζ cos α i sin α j sin α k , and ∂V l ∂Z l,i = (−R l,k R l,j m i )˛˛dθdζ cos α i cos α j cos α k (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, λ m,n 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 min{λ m,n SPEC } in terms of the stability parameter Z(δ SPEC v , υ) as, such that the stability conditions of Z(δ GGJ L , υ) can be interpreted the same as described before.

Numerical observation of marginal stability threshold
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 where υ ∈ Z + is the shear parameter, J 0 = 2B 0 /q 0 R 0 and p 0 = β p (B 0 a/q 0 R 0 (1 + υ)) 2 . The analytical expression for poloidal plasma beta is obtained as ) 2´a 0 rp(r) dr. We investigate this equilibrium model in the scenario with q 0 = 1.1, β 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. , which is understood in terms of equation (37). For υ = 2, we observe that the Z(δ SPEC v , 2) predicts instability as it crosses the marginal stability threshold line (Z = 1), for sufficiently small value of δ SPEC v , which is approximated as 2.8 × 10 −3 . In addition, the value of δ GGJ L at which the Z(δ GGJ L , 2) crosses its stability threshold coincides with that at which Z(δ SPEC v , 2) is greater than 1. This confirms the potential relationship δ SPEC v ∼ δ GGJ L . Now, for υ = 3, it is observed that Z(δ SPEC v , 3) predicts instability for the smaller value of δ SPEC v than for υ = 2. This is because as υ increases, the current channel shrinks in the vicinity of the rational surface, and both δ SPEC v and δ GGJ L reduce. For υ = 3, the threshold δ SPEC v is approximated as 1.72 × 10 −3 , and similar threshold behaviour is found for Z(δ GGJ L , 3). Thus, the MRxMHD stability boundary is in agreement with the linear modified tearing mode theory. However, if δ SPEC v becomes sufficiently large compared to δ 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.
Finally, figure 4(b) shows the spatial structure of the SPEC eigenfunction ξ · ∇s as a function of the effective radius r eff. . The corresponding m/n = 2/1 unstable equilibrium for ν = 2 and 3 are considered with δ SPEC v equal to 3.5 × 10 −3 and To obtain Z(δ GGJ L , υ), we consider the mass density ρ = 1, µ 0 = 1 and the Alfven speed v A = ∥ ⃗ Beq∥/ √ ρ 0 µ 0 in SI units, where ⃗ Beq is the equilibrium magnetic field. It follows that, in our units system, v A , τ A and r are unity, and thus η = S −1 , where S is the Lundquist number. The solid black line indicates the marginal stability threshold conditions Z = 1; (b) On the left axis: the SPEC computed perturbed surface displacement (ξ · ∇s) vs r eff ∼ ψt/ψ edge for unstable m/n = 2/1 mode; On the right axis: the SPEC computed q profile vs r eff ∼ ψt/ψ edge . Here, the number of the volume, Nv = 180, is considered in SPEC.
2.5 × 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.

Conclusion and future work
In this article, we have investigated the impact of the variational energy principle of the MRxMHD model to predict the finite-pressure linear tearing stability of tearing modes.
For low-pressure plasma, we have investigated a technique with which we have been able to establish a relationship between the resistive singular-layer theories of CGJ, GGJ and the MRxMHD model. Our analyses show that the SPEC shows stabilizing effects as the width of the resistive volume layer is decreased. Indeed, if δ SPEC v ∼ δ CGJ L and ∼ δ GGJ L , that is, the effects of finite resistivity and the pressure gradient roughly compensate, and the overall marginal stability of the mode is the same for SPEC, CGJ and GGJ. Physical insights into the spatial structure of the eigenfunction of the pressure-driven tearing modes computed from SPEC have clarified the applicability regime of the MRxMHD model. Our results indicate the possibility to couple the MATCH code [37], which solves the resistive inner-layer equations in toroidal geometry, with SPEC, not only to predict the stability of MRxMHD plasma but also to approximate the growth rates (quantitatively). This may also be a numerical treatment for determining the stability of coupled tearing modes, i.e. the modes (m, n) coupled with other poloidal mode harmonics (m + 1, n) and (m − 1, n) due to finite pressure effects and axisymmetric toroidal geometry.
In addition to these studies, we anticipate that it may be possible to establish a relationship between the pressure-flattening model (in the vicinity of resonant rational surfaces) discussed in [6,7] and our model. In fact, as pressure increases, it is commonly observed that the Mercier indices move apart and it becomes difficult to obtain the large and small solutions in the vicinity of the rational surface [6,38]. This restriction can be overcome in both MRxMHD and the Ham et al [6] model, due to the pressure flattening at the rational surface. We therefore aim to address this in our future investigations.
When using MRxMHD to predict non-linear tearing mode saturation, the difference between the potential energy corresponding to the equilibrium and the secondary minimized total energy can be interpreted as the second variation in the nonlinear stability case. Following this, in slab geometry, Loizu et al [39] demonstrated that the non-linear saturation of tearing modes can be predicted directly with SPEC using appropriate constraints, without resolving the complex resistivitydependent dynamics and without free parameters. It should be noted that while linear stability analysis based on the MRxMHD principle retrieves the linear tearing stability theory, the volume-averaged magnetic helicity is not an invariant as tearing modes grow non-linearly. In fact, the non-linear evolution of the resistive tearing modes is too slow for helicity to be well conserved, and instead a better invariant is the flux-surface average of the equilibrium toroidal current profile (see, e.g. Loizu et al [39]). To extend Loizu's work for finite beta cylindrical or toroidal plasma, our technique to compute an initial unstable MRxMHD equilibrium state can be utilized. We intend to investigate this work further in the future.

Data availability statement
The data cannot be made publicly available upon publication because no suitable repository exists for hosting data in this field of study. The data that support the findings of this study are available upon reasonable request from the authors.