Interplay of plasma resistivity and rotation on β limits in free boundary diverted tokamaks

A plasma resistivity-β driving mechanism aimed at explaining the appearance of long wavelength global instabilities in free boundary high-β tokamak plasmas with a divertor is presented. These perturbations resemble very closely the resistive wall mode phenomenon. Performing a proper toroidal analysis, we show that the magnetohydrodynamic stability is worsened by the interplay of plasma β and resistivity. By modelling the effect of a magnetic separatrix through a careful positioning of the resonant surfaces, we find that in an ideal plasma wall effects are effectively screened, so that the ideal β limit becomes independent of the wall position/physics. A lower wall dependent critical β is found if plasma resistivity is allowed. We find that global stability can be improved with a toroidal flow, small enough not to induce equilibrium modification. The rotation stabilisation effectiveness depends upon the proximity of the plasma equilibrium parameters to the resistive marginal boundary.


Introduction
One of the main goals of current tokamak research aims at maximising β, the ratio of plasma pressure over magnetic pressure.Reaching higher β value allows a larger fraction of bootstrap current and higher fusion power yield, both of which are extremely valuable for an economically viable reactor.However, global macroscopic magnetohydrodynamics (MHD) instabilities often limit the the maximum achievable β [1].The accurate evaluation of this critical β is thus of crucial importance for the design of long and steady pulses.
Usually, the maximum β value is extracted from MHD stability analyses performed with close fitting ideal wall boundary conditions.This is typically the highest β that can be reached.Indeed, the stability improvement due to an ideal wall is a well known effect, and this is particularly evident on external kink (XK) modes [2].Nevertheless, experimental evidence often shows the onset of long wavelength MHD activity when β exceeds a threshold which is smaller than the one predicted by the aforementioned MHD analyses with ideal wall boundary conditions [1].This macroscopic pressure driven instability exhibits external features and it is observed to grow on time scales of the order of tens of milliseconds [3] with a rotation frequency much smaller than the one of the bulk plasma [1].Furthermore, complete stabilisation can occur with a sufficiently fast plasma rotation [1,4].
Because of the external features of this MHD instability, the triggering mechanism is believed to be related to the excitation of an XK enhanced by β effects.Since the experimental critical β is lower compared to the one obtained numerically with ideal wall boundary conditions, it is believed that wall resistivity plays a role in determining the instability onset by allowing magnetic flux diffusion.This is known as the resistive wall mode (RWM) whose growth in present machines is of the order of several milliseconds, compatible with the typical resistive wall diffusion times, and much slower than the no-wall XK [1].
In such a configuration, the XK component of the RWM, i.e. the contribution to the available potential energy from the vacuum region, plays a significant role in determining the mode stability.It is well known [2] that for an XK to develop, the mode resonant surface must occur in the plasma-vacuum gap.This is not an issue in limited toroidal plasmas, but the situation is radically different for a diverted geometry in which the safety factor profile q diverges at the separatrix [12].This divergence constrains any mode resonance m/n with m and n the toroidal and poloidal mode numbers respectively such that m/n > q min to occur within the plasma.This in turn imposes strict boundary conditions at the resonance which are rather different to those normally applied in limited plasmas.
Hence, in this work we propose an alternative model aimed at explaining the appearance of these high-β long wavelength instabilities properly accounting for a diverted geometry.We employ the infernal framework [13,14] with a core region of large pressure gradients and low magnetic shear, and a vacuum region between plasma and wall.In order to model the effect of a magnetic separatrix, no modes are allowed to have a vacuum resonance.This framework naturally embodies the toroidicity of the problem, and identifies plasma pressure as the key drive.We find that wall effects, regardless of the fact of whether is ideal or resistive, are screened in an ideal plasma.On the other hand, the inclusion of plasma resistivity is shown to worsen the stability, effectively identifying a smaller critical β compared to the one obtained in the ideal plasma approximation.By allowing for plasma resistivity, we retrieve both the external and internal kink-like features observed in the experiments [15,16].The inclusion of favourable averaged curvature effects, improves stability and permits to tune the mode spatial structure from tearing to kink-like.Furthermore, complete stabilisation can be achieved by including a weak sheared toroidal flow even with modest values depending on the proximity of the plasma parameters to the resistive plasma marginal boundaries.
Hence, the paper is organised as follows: in section 2 a description of the plasma physics model and the equilibrium profiles is given.Sections 3-5 are devoted to the derivation of the eigenmode equations in the sheared, vacuum and core plasma regions respectively.In section 6 the dispersion relation obtained by applying the appropriate matching of the eigensolutions across the different interfaces.A thorough discussion of the physical implications of the dispersion relation is given in section 7. Concluding remarks and future outlook are given in section 8.

Physical model
The geometry of interest is the one of a large aspect ratio tokamak with circular shifted magnetic surfaces of major and minor radii R 0 and a respectively with ε = a/R 0 ≪ 1.We assume a vacuum region separating the plasma from a metallic wall.We adopt the standard low-β ordering, that is β = 2µ 0 p/B 2 0 ∼ ε 2 , where p is the pressure and B 0 is the on axis magnetic field strength.The analysis is conveniently carried out in a right handed straight field line coordinate system (r, ϑ, ϕ) where r is a flux label with the dimensions of length, and ϑ (counter-clockwise in the poloidal plane) and ϕ the poloidal-like and toroidal angles respectively.The equilibrium magnetic field in the plasma is where ψ is the poloidal flux and F = F(r) measures the toroidal magnetic field strength.The plasma is modelled by the resistive MHD equations [17] ρ where v is the plasma MHD velocity, ρ and J = ∇ × B are the mass and current densities (having normalised µ 0 = 1), η the resistivity, and Γ = 5 3 is the adiabatic index.The absence of currents in the vacuum region implies that Allowing for the presence of a resistive wall at radial position b > a, the effect of the wall physics are modelled by means of the thin wall approximation [18].The wall position can vary from b/a = 1 (wall directly interfaced with the plasma) to b/a → ∞ (wall at infinity).
The safety factor profile is constant with value q = q 0 = m/n − δq in the core region extending for 0 < r < r 0 , whereas it increases parabolically for r > r 0 with q = [m/n − δq](r/r 0 ) 2 .We focus on m > n = 1 modes and let 0 < δq < 1.A brief discussion of cases with δq < 0 will be given.In the cylindrical limit, the associated toroidal current is constant for r < r 0 and vanishing for r 0 .The radial position r 0 , which is referred to as the current channel, separates the low-shear region where dq/dr = 0 from the sheared region characterised by rdq/dr ∼ 1 extending from r 0 to a.
In order to model the effect of a magnetic separatrix which induces a logarithmic divergence of q near the plasma edge, we impose that perturbations of helicity m/n and (m + 1)/n resonate within the plasma at radial positions r m = r 0 m/nq 0 and r m+1 = r 0 (m + 1)/nq 0 respectively with r m < r m+1 < a.One finds that the maximum allowed width of the current channel is Knowledge of the fine structure of q near the edge is not required as long as r m+1 is well localised inside the plasma.
We follow [19] in extending the safety factor in x-point geometry beyond the plasma boundary.We shall point out that geometrical implications on the metric tensor coefficients of the presence of an x-point are not considered.We allow for a sheared equilibrium toroidal MHD flow u = Ω(r)∇ϕ which is assumed to be sufficiently weak so that no centrifugal corrections to the equilibrium pressure and mass density profiles [20] are induced.We may take a piece-wise continuous pressure profile of the form employed in [21] as shown in figure 1(b).We consider accordingly a constant mass density profile with value ρ with a monotonically decreasing temperature profile with T i = T e = T such that T(r m ) > T(r m+1 ).However, it will be clear that pressure gradient effects will not play a role in the stability equations in the sheared region apart from inducing a different response at the harmonics resonances due to the temperature variations.Therefore, to simplify the algebra and yet retaining the key effect of core mode coupling we let p 1 , p 2 → 0 and parametrise both pressure and rotation profiles with a Heaviside step-function H, that is with 0 < r p < r 0 and r 1 < r Ω < r 2 , while retaining the effect of a different response at the mode resonances due to a diffuse temperature profile.These equilibrium profiles will be used in the next sections where the equations governing the plasma dynamics in the various plasma regions will be thoroughly discussed.

Sheared region
We start by analysing the high-shear region.Fixing the toroidal mode number n = 1, let us assume that the radial fluid perturbation ξ can be decomposed into a dominant harmonic with mode numbers (m, n) and two subdominant ones with Note that the (m − 1, n) harmonic does not resonate since q > (m − 1)/n by hypothesis.Thanks to the choice of a stepped pressure profile (see (7)), harmonic coupling is prevented because of the absence of pressure gradients [14].Hence, far from resonances r m and r m+1 , the fluid perturbation obeys equation ( ′ ≡ d/dr) where ℓ = m, m ± 1 and ι = 1/q.Note that this is essentially the small aspect ratio tearing equation, expressed in terms of the radial displacement which is linked to the perturbed poloidal flux through the relation It is immediate to verify that to leading order r Br ℓ = ψℓ .It is worth noting that even in the presence of small pressure gradients, but with a magnetic shear of the order of unity, coupling between the m ± 1 and m harmonics is small.This is because at leading order ξ m obeys (9) for r > r m .The fulfilment of the boundary condition ξ m = 0 either at the ideal wall or at infinity forces ξ m (r > r m ) = 0. Therefore, in this region, ξ m±1 is independent of ξ m .
The general solution of (9) with the shape of the safety factor discussed in the earlier section is formally written as where A ±ℓ and C ±ℓ are some constants which depend on the mode number ℓ, and the sign ± stands for r ≷ r ℓ .A detailed knowledge of the behaviour of the eigenmode solutions near the resonance is needed.As (11) approaches its associated resonance located at x = (r − r ℓ )/r ℓ we have which has been obtained by expanding the safety factor to second order in x.This expression can be conveniently written as where ∆ r,ℓ and ∆ l,ℓ contain plasma inertia contributions, i.e. the growth rate and resistive effects.Note that (11) deviates from the solution of ( 9) obtained with more generic q profiles.However, the near-resonance asymptotic behaviour is still expected to be of the form (12), apart from some logarithmic terms which nevertheless do not play a role in the asymptotic matching with the layer solution [22].The parameters ∆ r,ℓ and ∆ l,ℓ may take a different form when q has a different shape, but these quantities are not expected to be altered too dramatically, so that the results presented here should remain qualitatively the same.Thus, by comparing ( 12) and ( 13) we have The quantity ∆ l,ℓ is determined by matching with the solution in the inertial-resistive layer, whereas ∆ r,ℓ , which includes the wall physics, is obtained by joining smoothly with the plasma and vacuum solutions.This is worked out in the next subsection.

The vacuum-wall solution
At the plasma-vacuum interface we must impose the appropriate boundary condition to the solution of (9).Since this equation is essentially derived in the cylindrical limit, we may employ in the vacuum the cylindrical approximation for the Laplacian operator (cf equation ( 5)), so that the magnetic perturbation obeys equation [2] r r Br Using the thin wall approximation [18], the solution of the equation above can be cast as with where τ w = bd/(2η w ) is a characteristic wall diffusion time with d the wall thickness and η w the wall resistivity (recall that we use the normalisation µ 0 = 1).A more rigorous approach would give Br written in terms of Bessel functions [23].
Noting that the relation between ξ r ℓ and (r Br ) ℓ in the plasma can be written as (cf (10)) We can introduce, for mathematical ease, a fictitious vacuum displacement ξ which relates to the magnetic fluctuation through the relation above.Using this trick, the vacuum perturbation can be written in a form similar to (11) with a near-resonance asymptotic behaviour given by ( 13).This means that the vacuum perturbation obeys equation (9).Combining ( 16) with ( 12) and ( 13) finally yields It is worth stressing that ∆ r,ℓ is a bounded quantity as long as r ℓ < b, and we can approximate 1, which usually holds for the innermost resonating mode or if the wall is at very large distance from the plasma boundary.We notice that the expression of ∆ r,ℓ with b → ∞ coincides with the one obtained by letting γ → 0 with τ w finite.This means that the no-wall marginal boundaries are the identical to the ones obtained by including a resistive-wall.Finally, we point out that, as long as the mode resonance is sufficiently far from the boundary, one can integrate (9) across a even if q diverges in the narrow edge region.This is because 1/q → 0 in (9).We point out that in employing this procedure to obtain the boundary conditions at the plasma-vacuum interface, we ignored the geometry of the flux surfaces and we only considered the cylindrical contribution.However, in doing so, we retained the essential information of the divergence of q at the plasma edge.Note that this is similar to what has been done in [24], in which the stability analysis of a simplified cylindrical equation was performed on top of profiles obtained from the solution of the full toroidal equilibrium.Hence, this leads us to infer that ξ ℓ and ξ ′ ℓ are continuous at the plasmavacuum interface, meaning that (19) should hold with a separatrix as well.It is nevertheless important to stress that this a strong approximation, and a proper analysis would be better accomplished by numerical tools designed to deal with geometries with x-points [19,25].

The layer eigenfunction
Let us consider the harmonic with poloidal mode number ℓ. Close to its associated resonant surface r ℓ , i.e.where ℓ/q − n = 0, we allow for plasma inertia and Glasser-Greene-Johnson (GGJ), namely curvature, effects, the latter denoted by the symbol ν.The quantity ν is not specified as it is used as a free parameter to tune the strength of the GGJ stabilisation, although it is assumed to be β dependent.Note that, as mentioned before, although the pressure gradient is vanishing in our model at the resonance r ℓ , GGJ effects are nevertheless expected to be present in a realistic situation, hence we retain ν corrections even if they are strictly speaking dependent upon p ′ 0 .A non vanishing plasma resistivity η is also taken into account.The layer analysis is more easily performed in Fourier k-space via the transformation [26] ξ where x is the variable defined in (13).Assuming the pressure to be weak enough not to induce coupling between neighbouring harmonics, the associated eigenfunction for the mode ℓ reads [26,27] d dy where y = k/ℓ, G = ℓ 2 η a 2 γ , and which has to be evaluated at r ℓ with ω A = B 0 /(R 0 √ ρ) the Alfvén frequency and s the magnetic shear.Here ν ≪ 1 plays the role of the resistive interchange parameter in [28].
The equation above has been solved in previous works [22,29], and the far from resonance real space asymptotic behaviour of the even (ξ e ℓ ) and odd (ξ o ℓ ) solutions is [22,30] where Γ is the Gamma function [31] and Thus, in order to match (21) with (13) we augment ( 9) by a GGJ-like term ν ≪ 1 yielding having defined k || = ℓι − n.Expressing (23) in terms of the perturbed poloidal flux ψℓ through (10), we perform a WKB expansion with small parameter 1/ℓ [32] which at leading order gives where in the last passage we exploited the fact that ν is small and we are not too close to r ℓ .Hence, far from the resonance, the term proportional to ν is dropped and one can approximate When r ℓ is approached, ν/k 2 || becomes larger, hence expanding about the resonance while keeping ν small we get Close to r ℓ in the limit ℓ 2 ≫ 1, we reduce (23) to where here ξ ℓ is to be intended as a function of x.Perturbing in ν, ξ ℓ can be formally expanded as ξ ℓ = ξ (0) + ξ (1) with ξ (1) /ξ (0) ∼ ν and solved to the first two orders in ν.This gives at leading order in ν where A and B are some constants, the upper/lower sign is for x ≷ 0 and Ei is the exponential integral [31].Since ξ ℓ ∼ ψℓ /x to leading order, the asymptotic behaviour of ( 26) is recovered by the expression above for large x.In the opposite limit x → 0, by applying the transformation Ei(−x) = −E 1 (x) for x > 0 and carrying out the appropriate expansions for ν ≪ 1 which are summarised in appendix A one has which has the same asymptotic behaviour of (21).Thus, by combining the equation above with (26) we can match the layer solution with (12) in the limit ℓ ≫ 1.Therefore, by letting ν ≪ 1 one finally obtains The unknown of ( 30) is ∆ l,ℓ with ℓ = m, m + 1.While ∆ r,ℓ has already been evaluated and the result given by ( 19), ∆ R,ℓ is computed from the solution of the layer equations.
Assuming that the temperature decreases with the radius, we allow for an ideal response at the inertial layer of the dominant mode m, whereas resistivity is allowed at r m+1 .We neglect Small pressure corrections in the ideal inertial layer at r m , and letting M → ∞ in (30) gives where ω A = B 0 /R 0 √ ρ and s m the magnetic shear at r m .For sufficiently small γ and Ω 0 , we let ∆ R,m ≫ ∆ r,m , so that from (30) one has Conversely, for the mode m + 1, if ν is small enough we have [27,29,33] where c 0 = 1 + 2(m + 1) 2 /n 2 /(ns m+1 ) with s m+1 the magnetic shear at r m+1 , S = a 2 ω A /η the Lundquist number, and We point out that with the choice of the safety factor above, as long as q 0 > m/n the tearing stability index of the mode m + 1 is never positive, that is no classical tearing modes can develop when ν → 0 [34] (the expression for the tearing stability index will be given later).It then follows that with ∆ R,m+1 given by ( 33).This equation is the equivalent of ( 30) but for the m + 1 instead of the m harmonic.Notice that ∆ R,m+1 → ∞ in the ideal (S → ∞) limit.

Low-shear region
In region 0 < r < r 0 , we allow for toroidicity driven mode coupling between a main mode ξ m and its satellite harmonics ξ m±1 [14,35].The standard infernal ordering is adopted with higher order harmonics ignored.We assume that the safety factor profile drops significantly below m/n so that inertial contributions can be neglected.Nevertheless, we assume δq = m/n − q 0 to remain sufficiently small in order to consider it as a small parameter ordered as (δq/q) 2 ∼ εα.This is somehow similar to the analysis of the m = 1 internal kink mode [36].We point out that dropping inertia is consistent within the framework of slowly growing modes.Hence, the model equations employed for the perturbation analysis for r < r 0 are [14,37] where α = −(2R 0 p ′ 0 q 2 )/B 2 0 and Note that we allowed δq to be small enough in approximating DM .Equation ( 37) can be integrated once yielding where L ± are two constants of integration.The regularity of the lower m − 1 mode on the magnetic axis implies that L − = 0, so that ξ m effectively couples with the upper harmonic ξ m+1 only.Therefore, we find that the equation for ξ m is [14] Because of the weak flow assumption, modifications to Q and DM due to toroidal rotation [20] are neglected.Using the standard techniques [20,35,37], the constant L + is obtained by solving for ξ m+1 in the region r p < r < r 0 and imposing smooth matching across r 0 .This yields [14,35] where c = r 0 ξ ′ m+1 (r 0 )/ξ m+1 (r 0 ).It is worth to point out that from (37) a discontinuity of ξ ′ m+1 across r 0 is expected if pressure gradients are allowed at this point.This is because of the sudden and sharp variation of the main harmonic at r 0 [36].

Dispersion relation
The dispersion relation is obtained by integrating (38) across r p giving where and ε p = r p /R 0 .The toroidal β obtained with the profiles given by ( 7) reads β = β * (r p /a) 2 /q 2 0 .Plasma inertia, namely the growth rate, and wall effects are contained in the first term on the left hand side of (40) and in the coefficient c.These two contributions are worked out as follows.
Exploiting the step-like behaviour of the pressure, we can solve (38) for ξ m on either side of r p giving where continuity at r p has been imposed.Note that ( 41) is valid for r < r 0 .The ratio a 0 /b 0 is determined by imposing the continuity of the logarithmic jump of ( 41) with the sheared region solution, i.e. equation (11), at r 0 .This is equivalent to the requirement that ξ m is smooth and continuous at r 0 .This gives where for the sake of simplicity we write Information about inertia at r m is contained in the unknown A m /C m .The smooth matching of the solution valid in region r 0 < r < r m with the one in the inertial layer is immediately obtained through equations ( 12), ( 13) and (32).One then sees that the ratio A m /C m evaluated from ( 14) in the limit ∆ R,m ≫ 1, i.e. small growth rate and toroidal rotation (cf equation ( 32)), yields Using ( 41) one gets to the first two leading orders in (44) where the coefficients c ± and d ± are given by By means of the formulae above, and expanding to leading order in δq, a little algebra shows that where the last approximation holds for r p /r 0 sufficiently small and δq not too close to unity.It is interesting to note that in case of no m/n surface with r 0 /b small enough and small growth rates, we must substitute (44) with [38] rξ This is briefly discussed in appendix B. This shows that with a flat core safety factor and no resonance field line bending stabilisation is weaker, and hence allowing the plasma to be more unstable [15].For the simple case of an ideal plasma, the mode structure in the low-shear region is obtained from equations ( 41) and (37).Few integrations across r p of the latter provide jump conditions for the sidebands to be combined with the requirement of regularity at the magnetic axis and smoothness at r 0 [37].The resulting structure of eigenfunction of the main mode compares favourably with the one observed in experiments and numerical simulations (see figure 2), when both the m/n resonance occurs in the plasma or in the vacuum region [4,15] (the shape of the upper sideband will be discussed later).
Using (11) we can evaluate the constant c giving where similarly to the notation introduced above, we let (A −,m+1 , C −,m+1 ) → (A m+1 , C m+1 ).By means of ( 14) one finally obtains where Z = (r m+1 /r 0 ) 2m+2 .Thus collating these results together allows to recast the dispersion relation equation (40) as We focus our attention on n = 1 modes, which are the ones of most interest.The expression above can be rearranged in a more compact form as follows where the meaning of each term appearing in ( 54) is now discussed: the magnitude of the ideal mode growth rate is measured by where we point out that the first term on the right-hand-side is always positive for 0 < δq < 1.Both A and B are positive for m > 1, and they read having neglected the δq dependence in q 2 0 B, which proves to be accurate enough if r p /r 0 is not too close to unity.It is found numerically that A ∼ B ∼ O (1).Finally, measures the stability of the system against tearing-like perturbations at the resonant surface of the m + 1 mode.Since ∆ r,m+1 is bounded, the quantity above can be taken to be of the order of unity.For δqZ sufficiently large and the wall at infinity (b → ∞) the expression above reduces to One sees that this expression of ∆ ′ T conforms to the classical free boundary tearing stability index of a perturbation with poloidal mode number m + 1 fulfilling the dispersion relation ∆ R,m+1 = ∆ ′ T .Note that the most tearing unstable case is obtained for b → ∞ (cf (19)) and δq = 0, which gives ∆ ′ T = 0.A thorough analysis of the dispersion relation equation ( 54) is carried out in the next section.

Analysis of the dispersion relation
We identify two limiting cases: the ideal and resistive plasma approximations.In the former case, the inertial response at both resonances r m and r m+1 is assumed to be ideal.Let us assume that the equilibrium toroidal flow is small and that the analysis is carried out in a neighbourhood of the marginal boundary (Re(γ) → 0).In such a case ∆ R,m+1 ≫ ∆ ′ T , so that (54) can be easily recast as where ω m+1 is computed from (31) with the obvious substitutions.One sees that the stability boundary is identified by the relation λ H = 0, and that the effect of rotation is only to produce a frequency Doppler shift.Since wall effect do not appear in (59) we can infer the following: because of the ideal response at each resonance, wall effects are completely screened, regardless of whether the wall is ideally conducting or resistive.This is a direct consequence of the strict boundary conditions to be applied to the structure of the eigenfunction at the r m and r m+1 resonances.This leads us to infer that in an ideal plasma the stability properties, namely the critical β beyond which MHD activity is triggered, are completely independent of the wall physics.
We shall now analyse resistive stability.In doing so, we discriminate between cases with ν = 0 and no rotation and cases for which ν ̸ = 0 with rotation.These are discussed in the following subsections.

Vanishing GGJ effects
us assume that ν = 0 with resistivity only taken into account at the r m+1 surface, and vanishing toroidal rotation.In such a case, (33) yields the classical tearing response ∆ R,m+1 ∼ γ 5/4 .Inspecting (54), it is immediate to see that in case of vanishing pressure, and thus A → 0, the dispersion relation reduces to the one for tearing modes, that is This is because λ H ̸ = 0 for β * → 0 if q 0 < m/n and 1/∆ R,m ∼ γ which can be made small enough in a neighbourhood of the marginal stability boundary.Recall that β * is linearly proportional to β.
Let us now assume that β * ̸ = 0.The marginal boundaries are identified by ∆ R,m+1 → 0 and ∆ R,m → ∞, which plugged into (54) give This is the resistive marginal boundary.In order to extend the applicability of (54), the quantity ∆ ′ T , negative in our model, may be regarded as a free parameter which can vary from −∞ to 0. Thus, it follows from (61) that the marginal λ H is negative with the stability increasingly worsened as ∆ ′ T becomes less negative.If ∆ ′ T → 0, one sees that λ H /A → −∞ which indicates that the system is always unstable.An example of the critical β as function of δq for different values of ∆ ′ T is shown in figure 3.
Contrary to the ideal case, wall effects enter through the quantity ∆ r,m+1 contained in ∆ ′ T .We note that in the case of an ideal wall (τ w → ∞), the closer r m+1 to b the stronger the wall effects (see figure 4).
Figure 4 also shows that moving the current channel away from the pressure step improves stability as higher values of T is computed consistently with the equilibrium parameters.
β at constant r p are reached (this is particularly evident in the b → ∞ case).The effect of the pressure peaking on the marginal q on the axis is shown in figure 5.
Our choice of the safety factor gives ∆ ′ T < 0 of (m − 1)/n < q 0 < m/n.Cases with ∆ ′ T > 0 obtained with more general shapes of q can be tackled by deploying a slight different approach for the solution of the sideband.According to [39], after performing a WKB expansion in the high-shear region, we may substitute where B0 < 0 and s 0 is the magnetic shear at r 0 + ϵ with ϵan infinitesimally small positive quantity.Using this result in (40) gives The marginal boundary is then obtained by requiring 1/∆ R,m = ∆ R,m+1 = 0.In such a case, the first term in the expression above is negative (cf equation ( 44)) and so is the last if ∆ ′ T > 0 indicating that no marginal threshold is present, i.e. the mode is always unstable.

GGJ and sheared rotation stabilisation
Let us now assume that ν ̸ = 0 in (cf equation (33)).Furthermore, we take λ H < 0, that is we analyse an ideally stable configuration, and assume that we are sufficiently far from the ideal Boundary.By means of (31), the dispersion relation equation ( 54) is then recast as having allowed for a sheared toroidal rotation of the form given by (7).Since the rotation is vanishing at r m+1 , this suggests the presence of two branches, one with rotation frequency Ω 0 and another which is not rotating in the lab-frame.We refer to the former as the fast-frequency root, whereas the second branch is called the low-frequency root.
For the fast-frequency root we write γ = iΩ 0 + δ where δ ≪ Ω 0 (for the sake of simplicity we take Ω 0 > 0).Plugging this form of the eigenvalue γ into (64) gives From (33), we see that ∆ R,m+1 ∼ S 3/4 γ 5/4 where γ ∼ iΩ 0 .Thus, for sufficiently large S, the second term on the righthand-side of the equation above can be made much smaller than |λ H |, so that this root is stable.
For the low-frequency root, we approximate γ − iΩ 0 ≈ −iΩ 0 in ∆ R,m [10].With this approximation and far from the ideal boundary (i.e.λ H sufficiently negative) equation ( 54) becomes By defining Ω = BΩ 0 /ω m , this equation can be written as where Ĉ0 = 2Γ(5/4)/Γ(3/4).We assume that the sum of the first two terms on the right-hand-side of the equation above is non-vanishing and greatly exceeds the absolute value of the third one, i.e. we assume that we are sufficiently far away from the resistive marginal boundary.Since M ∼ γ 3/2 , in the limit of ν/M small we write and define a characteristic growth rate We notice that γ T scales as S −3/5 , and corresponds to growth rates of the order of tens of milliseconds with ω A of the order of megahertz and S ∼ 10 6 − 10 7 .Hence, within the abovementioned assumptions we obtain By assuming a sufficiently small rotation, far from the resistive marginal boundary we can treat the second term in the bracket as a small correction, so that we exploit once more the smallness of the term proportional to ν giving where s m+1 measures the magnetic shear at r m+1 (in our model s m+1 = 2).Here we have approximated 1 + 2m 2 /n 2 ≈ √ 2m with n = 1.Thanks to finite GGJ effects, a threshold γ T , or equivalently |λ H | for fixed Ω, is introduced [29,33,40], and the limit Ω → 0 is identified by the relation We call this marginal boundary, identified by the dotted line in figure 6, the GGJ resistive boundary, and is seen to approach the ideal (λ H = 0) as ν is increased (this can be achieved by e.g. an increase of S).As discussed previously, we shall point out that (70) should be used when |λ H | is not too small, i.e. far from the ideal boundary.Noticing that the rotation frequency appears in the denominator of (68), we can invert (70) to obtain a critical marginal rotation which reads Although s m = 2 with the current profile of section 2, for modelling realistic smooth profiles one may approximate s m ≈ δq/m if δq is sufficiently small.Moreover, since ν is generally Proportional to the pressure gradient [28], we may approximate ν ∼ Cβ 5/6 where C ∼ 10-100 for S ∼ 10 6 -10 8 with m = 2 and s m+1 = 2.A plot of the critical rotation values for a model equilibrium is shown in figure 6.We note that stabilisation can be achieved for rather modest values of the rotation frequency as long as the plasma state is close to the GGJ resistive boundary.Furthermore, one can expect from (64) that as ν increases the local structure of the eigenfunction about r m+1 transitions from tearing-like to kink-like parity, approaching the one obtained for an ideal plasma (S → ∞) [4].This is because in absence of toroidal flows, with |λ H | → 0 implies ∆ R,m+1 → ∞ (see equation ( 64)), thus giving a response similar to that of an ideal mode at the sideband resonance.Finally, letting ∆ ′ T → 0 and Ω ≪ |λ H |, we see from ( 68) and (69) that Because of the S −3/5 dependence of γ T , it follows that Im(γ) ≪ Ω, i.e. the mode has a rotation frequency which is much smaller compared to the one of the bulk plasma.

Conclusions
In this work we analysed the impact of plasma resistivity on the marginal stability boundaries of long wavelength global instabilities in free boundary toroidal geometry.The effect of a magnetic separatrix has been taken into account through a careful choice of the positioning of the relevant resonances, which have been constrained to occur within the plasma region not in the vacuum.By adopting stepped profiles, which allow for a full analytic treatment performed within the infernal mode framework, it has been found that with an ideal plasma wall effects are completely screened by the resonance ideal responses.This implies that the threshold β above which this global instability is triggered carries no information on the wall physics.This rules out current driven XKs a possible trigger for the instability.However, if plasma resistivity is allowed at the outermost resonance, a lower marginal β is found.This depends upon the wall physics through a modification of the perturbation magnetic tail extending in the vacuum region.It has been found that the closer the outer resonance to the plasma boundary, the stronger the wall effects.Favourable curvature effects improve stability, and can tune the response of the outermost resonance making it transition from tearing to kink-like.We found that the structure of the associated eigenfunction compares favourably with the one observed both in experiments and numerical simulations [4,15].
By including a sheared toroidal flow, we identified a slowgrowing root with growth time-scales proportional to S −3/5 (of the order of the order of several milliseconds for experimental tokamak parameters) whose rotation frequency is significantly smaller than the one of the bulk plasma.This root is stabilised by the shearing of the rotation, and the magnitude of the flow needed to achieve stabilisation is found to be of few percent of the Alfvén speed depending on the proximity to the resistive marginal boundary obtained with the inclusion of curvature effects but no rotation.
Several features observed experimentally (such as critical β values, mode structure, growth time-scales, rotation stabilisation) are retrieved.Further work is nevertheless required to assess more thoroughly the perturbation ballooning nature exhibited near the edge, where the accumulation of several resonances occur, and the dependence of mode stability on more realistic (smooth) current and pressure profiles.Possible comparisons with MHD stability codes may be addressed with a careful tailoring of the mode spectrum in order to avoid harmonics resonating in the vacuum, or with numerical frameworks which include intrinsically separatrix effects in the equilibrium geometry.

Figure 1 .
Figure 1.Model safety factor (a), and pressure and toroidal rotation profiles (b).In panel (b), the units of the y-axis are arbitrary.

Figure 4 .
Figure 4. Critical β for the dominant m = 2, n = 1 mode as function of the ideal wall distance for three different values of the current channel extension r 0 /a.Here δq = 0.3, rp/a = 0.4 and ε = 1/3.The radial position of the resonance of the m + 1 mode is indicated.∆ ′T is computed consistently with the equilibrium parameters.

Figure 5 .
Figure 5. Critical δq for the dominant m = 2, n = 1 mode as function of the pressure peaking pax/ < p >= (a/rp) 2 for three fixed values of the β in an ideal (a) and resistive (b) plasma.Here r 0 /a = 0.6, b/a = 1.05 and ε = 1/3.With these parameters, the wall position has a weak influence on the resistive boundaries.The region of instability lies below each curve.The angular brackets indicate volume average, that is ⟨•⟩ = ´V r(•)dV/ ´V rdV where dV is the volume element.

Figure 6 .
Figure 6.Critical rotation values for the resistive mode stabilisation for dominant mode (m = 2, n = 1) with r 0 /a = 0.6, rp/a = 0.4, ε = 1/3 and b → ∞.Here we let ∆ ′ to be a free parameter and we set ∆ ′ T = −0.1.The ideal and resistive (no GGJ with Ω 0 = 0) marginal boundaries are identified by the solid and dashed curves respectively.Ideal instability occurs in region A, whereas region B is made stable by GGJ effects with C = 10.By using this value of C to compute ν, the marginal points obtained from (73) are indicated by the dotted line.Equation (74) holds below the dot-dashed line which identifies the A/|λ H | + ∆ ′ T − ν ≈ 1 level.