This study explores the effects of plasma shaping on magnetohydrodynamic mode stability and rotational stabilization in a tokamak, including both plasma and wall resistivity. Depending upon the plasma shape, safety factor, and distance from the wall, the *β*-limit for rotational stabilization is given by either the resistive-plasma ideal-wall (tearing mode) limit or the ideal-plasma resistive-wall (resistive wall mode) limit. In order to explore this broad parameter space, a sharp-boundary model is developed with a realistic geometry, resonant tearing surfaces, and a resistive wall. The *β*-limit achievable in the presence of stabilization by rigid plasma rotation, or by an equivalent feedback control with imaginary normal-field gain, is shown to peak at specific values of elongation and triangularity. It is shown that the optimal shaping with rotation typically coincides with transitions between tearing-dominated and wall-dominated mode behavior.

## I. INTRODUCTION

Increasing *β*, the volume-averaged ratio of the plasma pressure to the magnetic pressure, increases the fusion power in a tokamak but also drives resistive and ideal magnetohydrodynamic (MHD) instabilities that can destroy the plasma confinement. Focusing on resistive wall modes (RWMs), stabilization methods include but are not limited to plasma rotation with respect to a resistive wall^{1–11} and feedback control.^{12–23} The aim of this paper is to explore the extent to which rotation can raise the MHD *β*-limit of a resistive plasma surrounded by a resistive wall in shaped toroidal geometry.

The study examines the linear onset of MHD instabilities with a broad poloidal harmonic spectrum and fixed toroidal harmonic n = 1, in a plasma that is stable at zero *β* and destabilized at finite *β*. Following Brennan and Finn,^{23} four *β*-limits calculated *without* rotation or feedback control are used to evaluate the extent to which rotation or feedback control can raise the *least stable β limit*—the first of the four limits to go unstable as *β* increases—in a resistive plasma surrounded by a resistive wall. Starting with an ideal-plasma ideal-wall (ip-iw) system, raising *β* produces a kink mode that goes unstable at a relatively high limit, denoted $\beta ip\u2212iw$, with a fast growth rate characterized by the Alfvén timescale *τ _{A}*. Wall resistivity allows the perturbed magnetic flux to penetrate the wall on a resistive timescale

*τ*, introducing a slower growing instability known as the resistive wall mode (RWM), which goes unstable at the

_{w}*ideal-plasma resistive-wall*(ip-rw) limit $\beta ip\u2212rw<\beta ip\u2212iw$. Similarly, plasma resistivity in a system with an ideal wall introduces yet another non-ideal instability known as the tearing mode (TM), which grows on a tearing timescale

*τ*and is also destabilized at a

_{t}*resistive-plasma ideal-wall*(rp-iw) limit $\beta rp\u2212iw<\beta ip\u2212iw$. Lastly, a toroidal system containing both wall resistivity and plasma resistivity become unstable at a

*resistive-plasma resistive-wall*(rp-rw) limit $\beta rp\u2212rw$, below the other three limits. The least stable mode (without rotation) appears when

*β*crosses $\beta rp\u2212rw$, and the mode grows on a timescale that depends upon both

*τ*and

_{w}*τ*, coupling the resistive-wall and the tearing processes by their comparable timescales as well as their similar physics and their mutually inductive perturbed currents. Rather than referring to separate modes, the present unified approach suggests that common designations such as kink, RWM, and TM should be thought of as referring to the dominant MHD mode behavior in different domains of the parameter space. This unified approach was pioneered in a cylindrical study by Finn

_{t}^{24}and a toroidal circular cross-section study by Betti,

^{9}who independently developed models to study MHD stability over a range of

*β*, encompassing both TM and RWM behavior. The approach of keeping both plasma and wall resistivity allows the dominant mode behavior to transition smoothly at each of the limits described above; $\beta rp\u2212rw,\u2009\beta rp\u2212iw,\u2009\beta ip\u2212rw$, and $\beta ip\u2212iw$.

For a circular cross-section tokamak, a typical *β*-limit ordering in this 4-*β* analysis was found to be $\beta rp\u2212rw<\beta rp\u2212iw<\beta ip\u2212rw<\beta ip\u2212iw$ in both the cylindrical model by Finn^{24} and the high aspect-ratio toroidal model by Betti.^{9} Richardson, Finn, and Delzanno^{25} found that the same ordering also applies in a typical reversed field pinch, with the current density parameter $\lambda 0\u2261J\xb7\beta /B2(r=0)$ in place of *β*. Building upon the 4-*β* approach, Brennan and Finn^{23} constructed a cylindrical tokamak model with feedback control, to show that the plasma response to rotation and/or feedback control is characterized by the four *β* domains associated with the four above-mentioned *β*-limits calculated *without* rotation or feedback. With cylindrical geometry that typically satisfies the *β*-limit ordering mentioned above, the study by Brennan and Finn indicates that the limit $\beta rp\u2212rw$ can be raised by rotation up to the next limit $\beta rp\u2212iw$. It was shown^{22} that feedback proportional to the radial field with complex gain is exactly equivalent to rotation of the resistive wall in a cylinder, and therefore to rigid plasma rotation in the opposite direction. The equivalence was shown to be approximate but still useful in the toroidal geometry. Earlier work in Ref. 16 had shown that with two resistive walls, rotating sufficiently fast relative to each other, resistive wall stability can be achieved because the flux cannot penetrate the two walls. A method called the fake rotating shell approach of Ref. 17 uses a feedback array to obtain similar stabilization. In Ref. 22, it was discussed that the equivalence of imaginary gain and wall rotation shows that imaginary gain with two resistive walls can stabilize completely.

While feedback with complex tangential and normal field gains can further stabilize the mode, the present study focuses on the parameter space that permits mode stabilization by rotation, or by an equivalent feedback control with imaginary normal-field gain. By varying tokamak parameters such as the wall radius (or closeness of the wall to the plasma boundary), safety factor, elongation, triangularity, and vertical asymmetry, this study examines the maximum extent to which rigid plasma rotation can raise the *β*-limit of the least stable mode. The parameter space where a non-rotating unstable plasma can be stabilized by rigid rotation is referred to as the rotationally *stabilizable* domain. Stabilizing rotation rates are typically on the order of *τ _{w}* or

*τ*. A general rule, first demonstrated in the present study, is that the

_{t}*β*-domain of rotational stabilization is bounded by the lowest ideal

*β*-limit, which can be given by either $\beta rp\u2212iw$ or $\beta ip\u2212rw$. This rule was examined for the specific case of varying wall distance in studies by Finn

^{24}and Betti.

^{9}This finding has implications for the type of ensuing linear mode behavior. In the common terminology, $\beta rp\u2212iw<\beta ip\u2212rw$ implies rotational stabilization up to the linear onset of a TM-dominated instability, whereas the reversed case $\beta ip\u2212rw<\beta rp\u2212iw$ implies rotational stabilization up to the linear onset of a RWM-dominated instability. The two different stability limit orderings, $\beta rp\u2212iw<\beta ip\u2212rw$ and $\beta ip\u2212rw<\beta rp\u2212iw$, play important roles in the present study and will be discussed at length later on. We suggest that the nature of the linear

*β*-limit ($\beta rp\u2212iw$ or $\beta ip\u2212rw$) determines the dominant type of non-linear behavior (TM or RWM) observed in experiments.

Rotational stabilization of the least stable mode (occurring when *β* crosses $\beta rp\u2212rw$) can be understood in the context of a coupled-mode picture. Namely, a TM with finite flux at the tearing layer and zero flux at the wall and a RWM with finite flux at the wall and zero flux at the tearing layer.^{23} Any system of two coupled modes exhibits a mode interaction which depends upon the proximity of the roots in the complex plane. Plasma rotation in the present plasma-wall system produces a relative phase-shift of the complex roots which changes the coupling of the two modes. This effect tends to raise the $\beta rp\u2212rw$ when the tearing layer and the wall have comparable timescales, with similar rotation timescales. (There is an important exception of destabilization by low rotation when tearing layers have real frequencies.^{26,27}) Rotational stabilization has been verified experimentally^{1–3} and explained theoretically using a number of dissipation mechanisms including sound wave damping,^{4} resistivity,^{5,6} and viscosity,^{7,8} as well as kinetic effects such as the resonance between mode rotation and the precession drift frequency of trapped particles.^{10,11}

The RWM-TM interaction is further modified by the geometric mode-coupling induced by the shaped toroidal geometry. In the past, geometric mode-coupling effects of cross-sectional shaping in toroidal geometry have been studied separately for ideal-plasma (kink) modes, and for resistive-plasma (tearing) modes. Early shaping studies focused on the ideal external kink instability, applying a sharp-boundary model with no rational surfaces and no wall. Using a sharp-boundary model with a high aspect-ratio and elliptical cross-section, Freidberg and Haas^{28} found maximum stability for an elongation (height to width ratio) of 2.2. An extension by Freidberg and Grossmann^{29} to a more general shape showed triangularity to be destabilizing in the absence of elongation.

Including internal rational surfaces in a diffuse plasma profile - but still neglecting resistive wall effects - a number of ideal MHD numerical shaping studies of the DIII-D experiment were undertaken by Lazarus *et al.*,^{30} Turnbull *et al.*^{31} and Kessel *et al.*,^{32} Ferron *et al.*^{33} and Holcomb *et al.*^{34} Similar studies were conducted by Menard *et al.*,^{35} Miller *et al.*^{36} and Turnbull *et al.*,^{37} to test for stability in a low aspect-ratio tokamak such as NSTX. These numerical studies all suggest that cross-sectional shaping - most notably a combination of elongation and triangularity - can help raise the stability limit, but generally included relatively few data points to discern an optimal shape or to analyze the physics trends from shaping. In contrast to the domain dominated by the external kink resonance, analytic studies by Bondeson and Bussac^{38} and by Lutjens, Bondeson, and Vlad^{39} showed that the *internal* kink mode, with toroidal number *n* = 1 and low poloidal *m*-numbers (typically *m* = 2), is *destabilized* by vertical elongation. More extensive models by Eriksson and Wahlberg^{40} and by Martynov, Graves, and Sauter^{41} showed that triangularity, on the other hand, stabilizes the internal kink.

Shaping studies focused on *resistive* plasma behavior showed that the TM is generally stabilized by both elongation and triangularity. A semi-analytic model known as the T7 code was developed by Fitzpatrick *et al.*^{42} and recently revisited by Ham *et al.*,^{43} who demonstrated a stabilizing effect of shaping on tearing modes surrounded by an ideal wall. The present paper presents the first shaping study of coupled RWM and TM phenomena spanning a wide range of *β* values.

To explore the effects of plasma shaping on the intrinsic stability limit as well as the limit of rotational stabilization, we have developed a shaped sharp-boundary model in toroidal geometry including resistive resonant surfaces and a resistive wall. The sharp-boundary approach makes it convenient to scan stability over an individual parameter (such as *β*, safety factor, wall distance, elongation, triangularity or vertical asymmetry) while keeping the rest fixed. In this manner, a space of $\u223c103$ stability eigenvalues—including different geometries—can be generated in approximately 1 h. The resulting qualitative study of a broad parameter space is designed to supplement and guide investigation by quantitative stability codes such as MARS^{44} or DCON.^{45}

The new model is based on a study by Fitzpatrick^{46} which incorporated tearing surfaces into a sharp-boundary formulation to examine the effects of shaping on error-field response. In the present model, the geometry in Ref. 46 is generalized to include vertical asymmetry, in order to capture the shape of a single-null diverted plasma. Additionally, the model incorporates toroidal curvature corrections based on an ideal stability model by Freidberg and Grossmann,^{29} retained up to the first order in an expansion in the inverse aspect-ratio *ϵ*. The development of a resistive wall boundary condition facilitates the generalization of the Brennan-Finn^{23} 4-*β* analysis to the shaped toroidal geometry. Using the 4-*β* framework, the new sharp-boundary model reveals qualitative trends of stability and rotational stabilization over a broad range of *β*, safety factor, wall location, elongation, triangularity, and vertical asymmetry.

A main result of this study is that plasma shaping can cause an interchange of the rp-iw and ip-rw *β*-limits, modifying the *β*-limit ordering from $\beta rp\u2212rw<\beta rp\u2212iw<\beta ip\u2212rw<\beta ip\u2212iw$ to $\beta rp\u2212rw<\beta ip\u2212rw<\beta rp\u2212iw<\beta ip\u2212iw$. While $\beta rp\u2212rw$ is always the lowest (least stable) limit and $\beta ip\u2212iw$ sets the upper bound, there is no constraint on the order of the two middle limits $\beta rp\u2212iw$ and $\beta ip\u2212rw$. The discovery of alternate *β*-limit orderings introduced by shaping becomes important in the context of previous observations that the lower of these two ideal limits, $\beta rp\u2212iw$ and $\beta ip\u2212rw$, sets the upper bound for rotational stabilization,^{9,23,24} with small exceptions discussed in Secs. IV B and V.

For a circular cross-section, an interchange of $\beta rp\u2212iw$ and $\beta ip\u2212rw$ was observed by varying the wall radius in stability studies by Finn^{24} and Betti.^{9} The present study reaffirms the mode interchange induced by varying the wall radius, and goes on to demonstrate likewise interchanges induced by varying the safety factor, elongation, and triangularity. This discovery shows that different domains of the tokamak parameter space exhibit rotational stabilization bounded by either TM or RWM type behavior. An optimal shape for stabilization by rotation—or an equivalent feedback with imaginary normal-field gain—was found to typically reside in a window around the transition from TM-limited to RWM-limited domains, which we identify by an interchange of *β*-limits.

The remainder of the paper is structured as follows: Sec. II outlines the new features of the sharp-boundary model, adapted from the formulation of Fitzpatrick.^{46} Section III introduces a new resistive-plasma resistive-wall dispersion relation including effects of shaped toroidal geometry. Section IV presents the resulting growth rate and stability limit calculations for the case of a circular cross-section toroidal plasma. Section V shows how the stability limits are affected by cross-sectional shaping, including first observations of *β*-limit interchanges as the shape is varied. Section VI illuminates a key distinction between ideal *β*-limits in the sharp-boundary model with and without resonant surfaces, and how they relate to ideal *β*-limits in present day simulations. Section VII details a preliminary verification of the qualitative stability features presented in this paper via numerical simulations with NIMROD^{47} and DCON.^{45} Section VIII summarizes the results of the paper. Appendixes A and B provide additional mathematical details of the model geometry and tearing layer response, respectively.

## II. MODIFICATIONS TO THE FITZPATRICK SHARP-BOUNDARY MODEL

The present formulation builds on a resistive-MHD sharp-boundary model developed by Fitzpatrick^{46} to study the effect of plasma shaping on error-field response. We generalize the geometry in Ref. 46 with vertical asymmetry to emulate a single-null diverted plasma (Sec. II A), as well as $O(\u03f5)$ toroidal curvature based on an ideal-MHD model by Freidberg and Grossmann^{29} (Sec. II B). Here, $\u03f5=a/R$ is the usual toroidal inverse aspect ratio. Lastly, the addition of a resistive wall (Sec. II C) facilitates the formulation of a new resistive-plasma resistive-wall dispersion relation in the shaped toroidal geometry (Sec. III).

The sharp-boundary model makes it efficient to scan individual equilibrium parameters—defined at the plasma boundary—without having to adjust the entire plasma profile. This is achieved by approximating the tokamak equilibrium current to consist of a skin-current at the plasma boundary, which results in a discontinuity in the tangential magnetic field. The jump in magnetic pressure across the plasma boundary is balanced by a jump in the fluid pressure, which follows a step function profile. In contrast with the analytic model in Ref. 23 which employed a reduced-MHD scaling (i.e., dominant constant *B _{z}* with small $\delta Bz$) that allowed for a step-function profile in both the pressure and current density, we relax the reduced-MHD scaling assumption so that a pressure step must be balanced by a step in $B2/2$, thus requiring a delta-function profile for the current density.

Stability calculations in the new model are outlined as follows: the perturbed field response is solved by numerically integrating the vacuum-like magnetic potential, subject to boundary conditions at three surfaces: (i) the plasma boundary obeys the standard ideal MHD conditions of the sharp-boundary theory,^{28,29,46,48–52} enforcing ideal Ohm's law on either side of the perturbed boundary. (ii) The tearing layers, following Fitzpatrick,^{46} are set just outside of the sharp plasma boundary where the presence of a poloidal field produces a finite safety factor *q*. The external region can be thought of as a cold plasma with no equilibrium current and a vacuum-like *q*-profile, as portrayed by Finn.^{6} Following Fitzptatrick,^{46} the resonant layers are compacted near the plasma boundary for numerical convenience without qualitatively modifying the physics of the TM response. While the model is constructed to accept any tearing layer regime, present calculations apply a constant-*ψ* visco-resistive^{53,54} (VR) layer condition. (iii) A resistive thin-wall boundary condition, conformal to a flux surface, is constructed similar to the VR tearing boundary condition but incorporating all poloidal harmonics rather than just a resonant component.

Using the solution of the linear perturbed field problem, the perturbed energy is formulated as a multi-harmonic resistive-plasma resistive-wall dispersion relation. The present calculations use 45 poloidal harmonics ($\u221222<m<22$), comparable with the number used in Refs. 46 and 29, and a fixed toroidal harmonic *n* = 1. The growth rate and mode structure of the dominant mode are found by numerically solving the resulting non-linear eigenvalue problem (i.e., contains different powers of *γ* multiplying the system matrices, as discussed later on). A similar method was implemented by Betti^{9} for the case of step-function current and pressure profiles in a circular cross-section torus. By including only the necessary ingredients for studying tearing and resistive wall physics in the shaped toroidal geometry, the sharp-boundary model facilitates broad qualitative studies, intended to gain physical insights and to guide the investigation using larger quantitative codes.

### A. Vertical asymmetry in cross-section

The new model generalizes the geometry of Fitzpatrick^{46} with a vertically asymmetric cross-section, allowing for shapes that emulate a single-null diverted plasma. As in Ref. 46, the model applies a set of right-handed cylindrical polar coordinates $(R,\varphi ,Z)$, with length scales normalized by *a*, the minor radius in the case of a circular cross-section. To describe a shaped cross-section, *R* and *Z* are related to the curvilinear coordinates *r* and *θ*, representing a radial-like coordinate relative to the magnetic axis *r* = 0 and a poloidal angle-like coordinate *θ* in the cross-sectional plane. The plasma boundary *r* = 1 is parametrized by

independent of the toroidal angle of symmetry $\varphi $. Throughout the paper, subscript “*a*” denotes a function evaluated at the plasma boundary. Figure 1 depicts these coordinates for a typical shaped cross-section. Appendix A explains how the radial extension of this parametrization ($R(r,\theta ),Z(r,\theta )$) is obtained from a high aspect-ratio expansion of the Grad-Shafranov equation, following Connor and Hastie.^{55}

The coordinate basis is defined by the normalized contravariant derivatives, $e\u0302r\u2261\u2207r/|\u2207r|,\u2009e\u0302\theta \u2261\u2207\theta /|\u2207\theta |,\u2009e\u0302\varphi \u2261\u2207\varphi /|\u2207\varphi |$. The elongation parameter *κ* follows the standard definition, being the height to width ratio of the plasma cross-section. The other cross-sectional shaping parameters *δ _{x}* and

*δ*, can be analytically related to top triangularity

_{y}*δ*and bottom triangularity

_{t}*δ*as defined later in Eq. (26) via Taylor expansion of the plasma surface near the points

_{b}*θ*= 0, $\theta =\pi /2$, and $\theta =\u2212\pi /2$, respectively

These “model” definitions of up and down triangularity are used throughout the remainder of the paper, and agree with the standard geometric definitions in Eq. (26), in the limit $(\delta t/\kappa )2\u223c(\delta b/\kappa )2\u226a1$.

To approximate reasonable equilibrium bounds of the shaping parameters, the cross-section is required to maintain a positive curvature. For the vertically symmetric case where $\delta y=0$, this requirement is expressed as $\u2202Ra\u2202\theta |\theta \u2192\pi \u2212<0$. The resulting bound for triangularity is found to be $\delta \u0302<0.54$ (or $\delta x<1/4$), compatible with the typical range of ITER baseline scenarios^{56} with elongation $1.7<\kappa <2.0$, triangularity $0.3<\delta <0.5$, and vertical symmetry.

### B. Toroidal curvature corrections

In addition to vertical asymmetry, the present geometry extends the applicability of Fitzpatrick's model^{46} with toroidal curvature corrections, following Freidberg and Grossmann.^{29} The toroidal field, both inside and outside of the plasma boundary, falls off as $R0/R=1/(1+\u03f5g(r,\theta ))$. The geometric function $g(r,\theta )$ (see Appendix A) appears throughout the present formulation, containing the $O(\u03f5)$ effects of the shaped toroidal geometry on the equilibrium, as well as the perturbed field calculations described in Sec. II C below.

All magnetic field magnitudes are scaled by the vacuum on-axis toroidal field strength *B*_{0}. The internal domain (*r* < 1) contains only a toroidal field since the bulk of the equilibrium current is confined to the plasma boundary (*r* = 1). A constant *B _{i}* defines the relative magnitude of the internal magnetic field along the magnetic axis ($\theta =\pi /2$). With a step-function pressure terminating at

*r*= 1, the boundary poloidal field is determined by the equilibrium pressure balance $B\u03022(r=1+,\theta )=2\mu 0p+B2(r=1\u2212,\theta )$, with $Br(r=1)=0$ and the pressure

*p*constant on the flux surface

*r*= 1. This equilibrium condition, based on the free parameters

*p*and

*B*, as well as the inverse aspect ratio and cross-sectional shape parameters, is used to determine the poloidal field distribution just outside of the plasma boundary

_{i}Here, *α _{a}* is a geometric constant resulting for a flux surface average. The constant

*B*is computed by fixing a value of the edge safety factor

_{i}and numerically solving this transcendental equation. Thus, the equilibrium is characterized by *β* and *q _{a}*, as well as the shape parameters

*ϵ*,

*κ*,

*δ*, and

_{x}*δ*.

_{y}### C. Perturbed field response including a resistive wall

The third major extension to the Fitzpatrick^{46} sharp-boundary model is the addition of a resistive-wall, designed to be conformal to a flux surface *r *=* r _{w}*, geometrically similar to an ideal-wall formulation by Goedbloed.

^{51}

Similar to Ref. 46, the perturbed magnetic field both inside and outside of the plasma boundary is characterized by a scalar potential $\delta B=i\u2207V$. The sharp-boundary current profile results in a curl-free perturbed field satisfying Laplace's equation everywhere except for the layers of surface current. These current sheets, which in Ref. 46 include the plasma boundary and the resonant tearing layers, are now joined by a thin resistive-wall boundary condition

The condition in Eq. (6), expressed in terms of the wall dissipation time *τ _{w}*, is derived analogously to Fitzpatrick's tearing condition except with no rotation-induced Doppler shift since the wall is stationary. The resistive wall is induced with perturbed currents in response to all poloidal harmonics, in contrast with the tearing layers which each respond to only a single

*resonant*harmonic. Shape dependence of the geometric coefficients $g(r,\theta )$ and $h(r,\theta )$, defined in Appendix A, introduces coupling of the poloidal harmonics. The second condition, Eq. (7), enforces continuity of $\delta Br$.

Expanded in the standard Fourier basis, the wall boundary condition in Eq. (6) becomes

The shaped wall coupling matrices

are expressed in terms of a Fourier coupling operator

In the vertically symmetric case $\delta \u0302t=\delta \u0302b=\delta \u0302$, it suffices to calculate Fourier integrals over $0<\theta <\pi $ and take only the real (cosine) component, while the general asymmetric geometry requires evaluation over the entire circumference and inclusion of the imaginary (sine) component. An important feature of the coupling matrix $U(wall)$ is the presence of the geometric terms $gw(\theta )\u2261g(rw,\theta )$ and $hw(\theta )=h(rw,\theta )$, which couple wall currents of different harmonic structures. In addition, the wall boundary condition contains a geometric time factor $\tau \u0302(wall)$ associated with the wall time *τ _{w}*. Both the toroidicity matrix $U(wall)$ and the geometric time factor matrix $\tau \u0302(wall)$ become diagonal in the cylindrical limit. Geometric coupling in the wall plays a key role in this sharp-boundary model, where the total mode structure is strongly affected by the interdependence of the tearing layer and the wall. The harmonic structure of the wall currents is now captured by the jump of the tangential field

written in terms of the resistive wall response matrix

The wall jump condition in Eq. (12) along with the $\delta Br$ continuity condition in Eq. (7) enter the $\gamma \tau w$-dependent resistive-wall vacuum response matrix

with the wall contribution given by

The first term in Eq. (14) represents the no-wall vacuum response, recovered in the limit $rw\u2192\u221e$. The ideal wall limit $\tau w\u2192\u221e$ forces the condition $\delta Br(rw)=0$ at the wall. In this case, the second term in Eq. (15) vanishes. The uncoupled ideal-wall response of Freidberg and Haas^{48} is recovered in the present formulation by taking the cylindrical limit, whereby

Since the magnitude of the wall factor above is larger than unity, the vacuum gains a positive contribution to the perturbed energy and the effect of the ideal wall is always stabilizing.

Similar plasma boundary and tearing surface conditions are derived by taking different limits of the induction equation. The total perturbed field response is solved by numerically integrating the vacuum-like magnetic potential, subject to boundary conditions at three surfaces: (i) the ideal plasma boundary as in Ref. 29, (ii) the tearing layers as in Ref. 46 (see Appendix B for more details), and (iii) the new resistive thin-wall boundary condition. The response on either side of the plasma boundary is expressed as a multi-harmonic relation between the perturbed tangential and normal fields

The geometric terms define the Fourier-expanded Laplace equation which is expressed by geometric coupling matrices, calculated by FFT (see Appendix A). In the limit of a cylindrical plasma with no wall and no resonant surfaces, internal $r|m|$ and external $r\u2212|m|$ solutions lead to the diagonal response matrices $Cmk(plas)=|m|\u22121\delta mk$ and $Cmk(vac)=\u2212|m|\u22121\delta mk$, respectively. These simple solutions define a basis for the radially launched solutions for each harmonic in order to span the solution space. Broken into a set of coupled first-order ODEs, the system is integrated for each initial condition by a fourth order Runge-Kutta algorithm.

Note that the tearing layer, based on the formulation in an error-field response model by Fitzpatrick,^{46} contributes to the field response of the vacuum region at the plasma edge, which can be thought of as a resistive plasma with negligible equilibrium current as portrayed by Finn.^{6} Additional details of the tearing layer formulation are found in Appendix B.

## III. RESISTIVE-PLASMA RESISTIVE-WALL DISPERSION RELATION

With the formulation of a resistive-wall and shaped toroidal geometry, we are able to formulate the stability problem as a multi-harmonic resistive-plasma resistive-wall dispersion relation, similar to a circular cross-section toroidal formulation by Betti.^{9} Using the field response described in Sec. II C, the total perturbed force matrix takes the form

which contains contributions of the plasma, vacuum (including the resistive wall) and resonant tearing response. The resistive wall and tearing terms are *γ*-dependent and therefore make the force matrix non-self-adjoint. As discussed later in this section, this *γ*-dependent force matrix leads to a non-linear eigenvalue problem. Following Freidberg and Grossmann,^{29} the three matrices *H*, *G*, and $G\u0302$ associated with the surface, internal, and external solutions, respectively, are given by

again using the notation of the Fourier coupling operator in Eq. (11). The energy matrices above make use of the radial dependence of the poloidal field

calculated by Ampere's law, as well as the radial dependence of the metric functions

As in Ref. 29, this form of the force matrix utilizes a perturbed plasma displacement expansion $\xi (\theta ,\varphi )=(1+\u03f5ga(\theta ))\u2211m\xi mei(m\theta \u2212n\varphi )$, with the $(1+\u03f5ga)$ factor included to preserve symmetry of the energy matrices. In terms of the non-ideal perturbed force matrix, the perturbed force balance can be expressed as a non-linear eigenvalue problem. In contrast with the ideal model of Freidberg and Grossmann,^{29} this new non-ideal dispersion relation produces mode behavior on three different timescales, *τ _{A}*,

*τ*, and

_{w}*τ*, as well as the transitions between them.

_{t}The large aspect ratio approximation taken in the model, together with plasma incompressibility (the most unstable perturbations are incompressible) lead to a constraint equation that eliminates the *m* = 0 component of the perturbed normal displacement at the plasma boundary.^{28,46} Following Freidberg and Haas,^{28} a projection operator $P=I\u2212gg\u2020$ is constructed to eliminate the subspace of compressible solutions, with $Pg=0$, where $g\u2261G0m/|G0m|$, and G is defined in Eq. (21). The constraint above is automatically satisfied by with a transformed force matrix $PFP$, and we find that the non-linear eigenvalue problem becomes

where $\gamma D=\gamma +i\Omega $ and Ω is the Doppler shift due to the rigid plasma rotation. This equation is a generalization of the eigenvalue formulation of Freidberg and Haas^{28} to a non-linear eigenvalue problem with a *γ*-dependent force matrix.

The solution is obtained by a direct approach of minimizing the absolute value of the system determinant with respect to *γ*. In the absence of plasma rotation, the dominant mode *γ* is purely real.^{8} Solutions with complex frequencies that appear due to additional tearing layer physics^{27,57} are beyond the scope of the present study. The mode with the largest *γ* is found by initiating a 1-D Newton solver at large *γ* and converging to the solution from above. In order to find the dominant mode, the initial calculation is carried out at large *β* where the dominant growth rate lies well above zero, dominated by ideal-MHD physics. This first calculated growth rate is then used as an initial guess for further calculations as *β* is ramped down toward marginal stability (*γ* = 0). The same method is applied with other parameters as well (e.g., *q _{a}* in the inset of Fig. 2). Calculating the entire growth rate curve, rather than just the marginal stability points, has the benefit of showing how the mode transitions between the ideal, wall, and tearing timescales. The addition of plasma rotation Ω introduces a frequency component in the dominant mode and thus requires a 2-D search over complex

*γ*. Rotation is varied starting with the same initial calculation at high

*β*and zero rotation. Any point on this branch of rotation values (at fixed high

*β*) can be then used to initiate a

*β*scan in the presence of fixed rotation. A reliable 2-D minimizer for this application was determined to be the Nelder-Mead algorithm

^{58}found in the Scientific Python (

*SciPy*

^{59}) library.

## IV. STABILITY OF CIRCULAR CROSS-SECTION TOROIDAL PLASMAS

The purpose of the study in this section is to determine the effects of toroidal geometry on MHD mode stability and rotational stabilization in the presence of tearing layers and a resistive wall. As discussed in the introduction, the least stable mode in a system with both plasma resistivity and wall resistivity goes unstable at a resistive-plasma resistive-wall (rp-rw) *β*-limit with zero plasma rotation. We demonstrate in this work how this least stable limit, denoted $\beta rp\u2212rw$, can be increased by rotation (comparable with the tearing-time or the wall-time) up to the first (lower) ideal limit; either the resistive-plasma ideal-wall (rp-iw) limit $\beta rp\u2212iw$ or the ideal-plasma resistive-wall (ip-rw) limit $\beta ip\u2212rw$. An ideal-plasma ideal-wall (ip-iw) limit $\beta ip\u2212iw$ sets the upper bound for the other limits. The space of rotational stabilization is thus evaluated based upon the order of the four *β*-limits ($\beta rp\u2212rw,\beta rp\u2212iw,\u2009\beta ip\u2212rw,\u2009\beta ip\u2212iw$) over a broad range of plasma parameters. The 4-*β* formalism provides information regarding the nature of the rotationally stabilized MHD mode, which can be dominated by either the tearing mode (TM) or resistive wall mode (RWM) physics.

### A. Ideal and resistive plasma growth rates versus *β*

As an initial qualitative comparison with the growth rate trends obtained by Betti,^{9} Fig. 2 displays the growth rate curves of the four branches (rp-rw, rp-iw, ip-rw, and ip-iw) versus *β*, along with a close-up inset of the four *β*-limits (marginal stability points) plotted logarithmically. The geometry in Fig. 2 is characterized by $\u03f5=0.25,\u2009rw=1.33$, and $qa=2.1$, and the dissipation times are $\tau t=5\xd7104$ and $\tau w=103$ in units of the Alfvén time. The ideal-plasma and ideal-wall limits are obtained by scaling up the dissipation times by a factor of 10^{7}. All present calculations include a single *m* = 2 surface by taking an edge safety factor in the range $2<qa<3$, with $q0>1$.

Throughout this paper, plots follow the following convention: the dashed versus solid lines distinguish between a resistive and ideal plasma, while the blue versus red distinguish between a resistive and ideal wall.

The ideal “dome,” reminiscent of Shafranov's reduced-MHD cylindrical analysis of non-resonant ideal kink stability,^{60} agrees qualitatively with the results plotted in Figs. 5 and 6 of Ref. 9. Reference 9 obtains a relatively small separation between the resistive and ideal plasma limits, as observed in the present model for calculations at a high aspect-ratio. Including $O(\u03f5)$ geometric terms consistently throughout, the present model produces four distinct tails of the ideal dome. The resulting *β*-limits are qualitatively similar to those observed by Brennan and Finn.^{23}

The additional stable domain observed for the ip-iw curve at high *β* is the well-known second stability regime.^{61,62} In models with a diffuse profile and full toroidal geometry, second stability results from the high magnetic shear on the outboard side produced by the Shafranov shift,^{63,64} whereas in the sharp-boundary model, a similar effect is created by the wall-induced increase in magnetic pressure.^{9}

Next, we consider the effect of varying the edge safety factor *q _{a}*, still keeping a circular cross-section. Although the

*q*-domain near external kink resonance is generally avoided, it is illuminating to examine this undesirable region in the context of the 4-

_{a}*β*framework. Recall that the ideal dome presented in Fig. 2, for $qa=2.1$, is dominated by

*m*= 2. Keeping $rw=1.33$, Fig. 3 shows how as $qa\u21923$, a lower ideal dome appears corresponding to the

*m*= 3-dominated ideal external kink. This is clearly distinct from the remaining larger dome at higher

*β*which, as in Fig. 2, represents an ideal internal kink dominated by

*m*= 2. The transition is captured in the figure inset by plotting the growth rate versus

*q*at fixed $\beta =0.15$, as in Ref. 60. Observe the analogy between the

_{a}*q*-dome and the

_{a}*β*-dome, corresponding to the current and pressure drive components of the instability. This ability to vary

*β*with

*q*fixed illustrates an advantage of the sharp boundary model.

_{a}Mode transitions, such as observed in the growth rate plot in Fig. 3, can result in a *β*-limit interchange between $\beta rp\u2212iw$ and $\beta ip\u2212rw$ as equilibrium parameters are varied.

### B. *β*-limits for rotational stabilization

This section demonstrates how the first two ideal limits $\beta rp\u2212iw$ and $\beta ip\u2212rw$ can be interchanged by varying the wall distance and edge safety factor in a torus with a circular cross-section. Shaped cross-sections are then considered in Sec. V, with elongation and triangularity also found to similarly modify the *β*-limit ordering.

To first illustrate the effect of rotation, Fig. 4 shows the rotational stabilization of the rp-rw mode for different values of the wall distance *r _{w}* (normalized by the minor radius in the case of a circular cross-section). As before the tearing time and wall time are given, respectively, by $\tau t=5\xd7104$ and $\tau w=103$. For zero rotation, the rp-rw limit is independent of the wall location, being the same as the resistive-plasma no-wall limit. Rotation is observed to raise the rp-rw critical

*β*curve and is thus stabilizing. A sharp transition is observed between the close-wall RWM-limited domain of

*r*and the far-wall TM-limited domain. Consistent with existing studies,

_{w}^{4–8,65}scanning

*r*with fixed

_{w}*β*reveals a finite window for rotational stabilization (marked by the red arrows). The

*r*-window is observed near the transition between the RWM-dominated and TM-dominated domain. The top curve, $\Omega \tau w=40$, defines the saturated bound (same as $\Omega \u2192\u221e$) for rotational stabilization.

_{w}In order to understand the non-monotonic behavior observed in Fig. 4, it is illuminating to examine the higher *β*-limits shown in Fig. 5. The domain of stability for a non-rotating mode is shaded in dark-gray, while the light-gray region shows where the mode is *stabilizable* by rotation or feedback control with imaginary $\delta Br$-gain. A sufficiently close wall is observed to raise $\beta rp\u2212iw$ up to arbitrarily high *β*, which clearly does not provide a complete picture of the rotational stabilization cutoff observed in Fig. 4. Since $\beta rp\u2212iw$ moves up as the wall moves inwards while keeping the $\beta ip\u2212rw$ limit (equivalent for Ω = 0 to the ideal-plasma *no*-wall limit) fixed, a switching of the two limits $\beta rp\u2212iw$ and $\beta ip\u2212rw$ is observed around $rw=1.18$. Similar plots are found in studies by Finn^{24} and Betti^{9} for different models. Below the transition point, which appears as a sharp knee in the light-gray shaded region in Fig. 5, the limit for rotational stabilization corresponds to the $\beta ip\u2212rw$ limit, associated with the onset of a RWM-dominated instability. Here, the non-rotating mode *β*-limit is equivalent to the ideal-plasma no-wall limit and is thus independent of *r _{w}*. For $rw>1.18$, the

*β*-domain where rotational stabilization is possible gradually shrinks with the decrease of the $\beta rp\u2212iw$ limit associated with the onset of a TM-dominated instability. As the wall moves outward (

*r*increases), the $\beta rp\u2212iw$ limit approaches the $\beta rp\u2212rw$ limit. These results demonstrate how a unified treatment of MHD modes over a wide range of

_{w}*β*exhibits both RWM-dominated and TM-dominated behavior. In the interest of extending the 4-

*β*study of Brennan and Finn

^{23}to the shaped toroidal geometry, we begin the shaping studies in the domain that is TM-limited for a circular cross-section, proceeding with $rw=1.33$ for the remaining calculations.

Careful comparison of Figs. 4 and 5 reveals a small exception to the general rule that rotational stabilization of the $\beta rp\u2212rw$ limit is bounded by the first ideal limit. Here, the maximum rotationally stabilized rp-rw limit ($\Omega \tau w=40$ curve) peaks slightly above the ip-rw limit, near the interchange of the ip-rw and rp-iw curves. This small stability window above the first ideal limit, first observed by Finn,^{24} can be explained as a result of complex mode resonances as in the work of Finn and Gerwin.^{26} Bondeson *et al.*^{66} examined this exception in the toroidal geometry and found an even smaller window of stability. A broader investigation of the parameter space should be considered before discounting this higher window of stability. However, for the purposes of this paper, we will ignore this small region of higher rotational stabilization. The appearance of the lower *m* = 3-dominated dome in Fig. 3 is associated with a sharp transition in the *β* limits, portrayed in Fig. 6. In this figure, the edge safety factor varies keeping the central axis safety factor fixed with $1<q0<2$. For $q0<qa<2$, there are no rational surfaces and so the resistive-plasma and ideal-plasma limits coincide. As *q _{a}* is increased above

*q*= 2, the $m/n=2/1$ rational surface is introduced, whereby the ideal plasma limits jump to a high value of

_{a}*β*, representing a mode dominated by an internal $m/n=2/1$ kink. The resistive-plasma modes are continuous at

_{crit}*q*= 2, being equivalent (at Ω = 0) to the non-resonant limits. As $qa\u21923$, and the smaller

_{a}*m*= 3-dominated ideal mode arises as shown in Fig. 3, both ideal-plasma limits are observed to drop off sharply toward their corresponding resistive-plasma limits. This effect produces a finite $qa\u2212\beta $ window for rotational stabilization, lying approximately in $2.0<qa<2.27$, limited by the rp-iw (tearing) mode. This window contains much higher ideal-plasma limits, corresponding to an ideal internal kink. For the shaping calculations that follow in Sec. V, the safety factor is fixed inside the window at $qa=2.1$.

Although not addressed by Brennan and Finn,^{23} the interchange of the two middle *β*-limits due to variation of *r _{w}* or

*q*can be observed even in their cylindrical reduced-MHD model. Next, we consider how cross-sectional shaping in the toroidal geometry can similarly influence the

_{a}*β*-limit ordering and thus, the maximum achievable

*β*by rotational stabilization. The qualitative similarities between the results of the sharp boundary model of this paper with a circular cross section and the cylindrical MHD numerical simulations of Ref. 23 (as well as with the step-function current and pressure profile results of Ref. 23) indicate the validity of the sharp boundary model approach.

## V. *β*-LIMITS FOR ROTATIONAL STABILIZATION IN SHAPED TOROIDAL GEOMETRY

Cross-sectional shaping is shown in this section to also influence the *β*-limit ordering. Starting with a realistic inverse aspect ratio of $\u03f5=0.3$, a wall located at $rw=1.33$ and a safety factor of $qa=2.1$, Fig. 7 shows how the critical *β* curves of the four branches vary with elongation *κ*. As before, the central axis safety factor is set to $1<q0<2$, so that with $2<qa<3$, there is only an $m/n=2/1$ surface included in the calculations. For a circular cross-section and $qa=2.1$, this yields an ideal mode dominated by *m* = 2 as discussed in Sec. IV A.

Figure 7 shows how, for these equilibrium parameters, the ideal plasma modes are predominantly *destabilized* by elongation beyond $\kappa =1.1$ for $\beta ip\u2212rw$ and beyond $\kappa =1.3$ for $\beta ip\u2212iw$. This behavior is consistent with shaping models of the ideal internal kink.^{38,39} The resistive plasma modes, on the other hand, are stabilized up to an elongation in the neighborhood of $\kappa =1.9$, in qualitative agreement with the non-resonant result of Freidberg and Haas.^{28} The *β*-limit in the presence of plasma rotation is observed to peak at $\kappa =1.86$, the interchange point of the rp-iw and ip-rw curves (where $\beta rp\u2212iw=\beta ip\u2212rw$), beyond which increasing elongation severely reduces the ip-rw limit and the resulting extent of rotational stabilization. The resulting optimal elongation exhibits a good agreement with the range of optimal values calculated for DIII-D by Kessel *et al.*,^{32} as well as typical values in the ITER design.^{56}

In order to highlight how vertical elongation can open a window of higher *β* in the presence of rotational stabilization, analogous to the *r _{w}*-window observed in Fig. 4, Fig. 8 plots the growth rates of the least stable (rp-rw) mode versus

*κ*. The plotted rotation values, normalized by the wall-time, are $\Omega \tau w=0,\u2009\Omega \tau w=2$, and $\Omega \tau w=6$. The fixed $\beta =0.11$ corresponds to a horizontal cut across the plot in Fig. 7, above the peak

*β*value of non-rotating stability $\beta crit=0.08$, but below the peak

*β*value in the presence of rotational stabilization, occurring at $\beta crit=0.13$. Starting from $\Omega \tau w=0$, where the mode is unstable for any elongation, an increase of rotation to $\Omega \tau w=2$ is observed to saturate the extent of rotational stabilization for any elongation above $\kappa >2.0$. For these highly elongated shapes, the ip-rw limit (associated with the linear onset of a RWM) is observed to block the mode from being completely stabilized. A slight exception is found near the transition point from the TM-dominated (left) branch to the RWM-dominated (right) branch of the $\Omega \tau w=2$ curve. Similar to the case discussed in Sec. IV B for

*β*versus

_{crit}*r*, this is another example of slight stabilization above the ip-rw (RWM) limit as discovered by Finn.

_{w}^{24}

We observe that the shaped case in Fig. 8 hits the transition point at $\Omega \tau w=2$ rather than $\Omega \tau w=40$ as observed in Fig. 4 for a circular cross-section. At lower elongation values, increasing rotation continues to stabilize the mode up to the rp-iw (TM) limit, saturating slightly above $\Omega \tau w=6$. The $\Omega \tau w=6$ also exhibits a slight crossing of the ip-rw limit near the transition point between the ip-rw and rp-iw curves, in this case for stable values of *γ*. Increasing the fixed *β* would shift all of the curves upwards, until the marginal stability ($\gamma \tau w=0$) line coincides with the bottom of the high rotation curve $\Omega \tau w=6$, which occurs at $\beta =0.13$. The point of maximal *β* in the presence of rotational stabilization is observed to coincide with the crossing of the rp-iw and ip-rw curves. This result suggests that, similar to the window of *r _{w}* for rotational stabilization (Fig. 4), the optimal window of elongation

*κ*is found around the transition from the TM-dominated to the RWM-dominated domain.

Next, Fig. 9 introduces triangularity $\delta \u0302$ (obtained from either $\delta \u0302t$ or $\delta \u0302b$ in (3) with $\delta y=0$), fixing elongation at the locally maximal value of $\kappa =1.86$. Starting at the mode interchange where the two middle *β* limits coincide, increasing triangularity is observed to create a slight separation and then another interchange at $\delta \u0302=0.29$ ($\delta x=0.14$). Beyond this critical triangularity, the ip-rw limit decreases rapidly and diminishes the *β*-limit achievable with rotational stabilization, similar to the high elongation effect observed in Fig. 7. While this value provides only a local optimum around a fixed aspect ratio and elongation, as well as wall position and safety factor, it provides a proof of concept for parameter optimization based on the maximal achievable *β*-limit in the presence of stabilization by plasma rotation or an equivalent feedback control with imaginary $\delta Br$-gain.

As with elongation, the peak rotational stabilization is found nearby the peak of intrinsic (non-rotating) stability. Further investigation is required to determine if this is a trend or mere coincidence. In contrast with the relatively high range of optimal triangularity, $0.5<\delta <0.8$ calculated for ideal external kink modes in DIII-D by Kessel *et al.*,^{32} present results suggest that such high triangularity may be detrimental in a real system with plasma dissipation. The triangularity range in the ITER design,^{56} $0.3<\delta <0.5$, is in the slightly more conservative range but still possibly beyond optimal stabilization by rotation or the equivalent feedback control with imaginary $\delta Br$-gain.

Lastly, starting with the (locally) optimal shaping parameters $\kappa =1.86$ and $\delta \u0302=0.29$ (or $\delta x=0.14$), Fig. 10 plots the four *β _{crit}* curves with increasing vertical asymmetry, measured by the difference between the approximate top and bottom triangularity given in Eq. (3) (see the diagram in Fig. 1). While the vertical asymmetry is found to destabilize the least stable (rp-rw) mode, it is not seen to substantially modify the relative height of the first ideal

*β*-limit. No

*β*-limit reordering is observed. Thus, although destabilizing, vertical asymmetry is not predicted to severely reduce the effectiveness of rotational stabilization or the equivalent feedback control with imaginary $\delta Br$-gain. One more observation is the increase in the ip-iw limit as the asymmetry is increased. This opposite behavior of the rp-rw and ip-iw limits should serve as a caveat for ideal-plasma ideal-wall models, which may falsely conclude that vertical asymmetry is generally stabilizing.

## VI. IDEAL *β*-LIMITS IN THE SHARP-BOUNDARY MODEL WITH AND WITHOUT RESONANT SURFACES

An important subtlety of the analysis presented in this paper is the distinction between ideal *β*-limits obtained in the present sharp-boundary model *with* resonant surfaces and the ideal *β*-limits found in previous sharp-boundary models *without* resonant surfaces.^{28,29,48} Resonant-ideal boundary conditions—as used in the present model—enforce $\delta Br=0$ for each resonant harmonic on its respective surface, whereas the previous models^{28,29,48} have no internal resonant layers. In the present model where the rational surfaces lie at the plasma boundary (just outside the equilibrium current layer), the resonant ideal limit $\tau t\u2192\u221e$ shields resonant perturbations from reaching the wall. The non-resonant case (i.e., the sharp-boundary model without rational surfaces) is recovered in this model by taking the limit $\tau t\u21920$. Figure 11 plots the four resonant branches along with the non-resonant resistive-wall and non-resonant ideal-wall branches. Here, the cross-section is taken to be circular, with $qa=2.1,\u2009\u03f5=0.25,\u2009rw=1.33a$ and Ω = 0. The finite dissipation timescales are set in units of Alfvén time to $\tau w=103$ and $\tau t=5\xd7104$. In the limit $\tau t\u21920$, the new model presented in this paper recovers the classic non-resonant ideal kink *β*-limits of Freidberg and Haas^{48} (both with and without a perfectly conducting wall); see Fig. 11. Here, the cross-section is taken to be circular, with $qa=2.1,\u2009\u03f5=0.25,\u2009rw=1.33a$, and Ω = 0. The finite dissipation timescales are set in units of Alfvén time to $\tau w=103$ and $\tau t=5\xd7104$. In the opposite limit $\tau t\u2192\u221e$, the resonant boundary condition becomes similar to the ideal boundary conditions of modern stability codes such as DCON.^{45} The resistive-plasma (finite *τ _{t}*) branches, while connecting to the resonant-ideal branches at high $\gamma \tau A$, are observed to coincide with the non-resonant ($\tau t=0$) limits, independent of the value of

*τ*. Just as the no-wall and resistive-wall stability limits coincide,

_{t}^{4–8}so do the non-resonant ($\tau t=0$) and resonant resistive (finite

*τ*) limits. This implies that, at zero plasma rotation, the resistive-plasma values of

_{t}*β*in this paper can be compared with previous results of non-resonant ideal models.

_{crit}## VII. COMPARISON WITH NIMROD AND DCON RESULTS

A preliminary verification of the qualitative stability features presented in this paper has been obtained via numerical simulations using the model DIII-D shaped equilibrium shown in Figs. 12 and 13 using NIMROD,^{47} together with computations using DCON.^{45} Using the geometric definitions of top and bottom triangularity found in both codes, i.e.,

where $a\xaf=(Rmax\u2212Rmin)/2$ and $R\xaf=(Rmax+Rmin)/2$, the DIII-D equilibrium shown in Fig. 12, the standard triangularity calculated by this method on the computed equilibrium in Fig. 12 is $\delta t,b=0.77$. However, the standard calculation does not take into account the *δ _{x}* factor used in this paper. Also, the analytic form detailing the shape in Eqs. (1) and (2) cannot exactly match the shape of this numerical equilibrium which is based on a DIII-D experimental reconstruction. For a more meaningful comparison, the shape of the equilibrium in Fig. 12 can be approximated using the standard definitions for the geometric parameters $R0=1.67$,

*a*= 0.645, and $\kappa =1.86$, and using $\delta x\u223c0.17$ to give an effective $\delta \u223c0.35$. The profiles of the equilibrium are parametric and chosen to represent a realistic (L-mode) case which is stable to resistive MHD without a wall at low

*β*. The pressure in this case $P=P0\u2009exp\u2009(\u22122\psi )+Pe$, where

*P*is specified to reduce the edge pressure to zero. At low

_{e}*β*, the current profile is specified with a parametric function to be flat in the core and rapidly reduces to zero at the boundary, giving a monotonically increasing

*q*profile. The total current is then specified to set the safety factor on axis $q0=1.25$ with a toroidal field of 1

*T*at

*R*

_{0}. Given this stable equilibrium at low

*β*, the pressure is increased multiplicatively while holding the

*q*profile fixed, finding a new current profile to solve the equilibrium. Resultant profiles for a higher

*β*equilibrium are shown in Fig. 13, where a moderate variation in the current profile shape can be seen.

For the NIMROD simulations, the DIII-D wall structure as shown in Fig. 12 is used, the Lundquist number is $S=107$ with Spitzer resistivity $\eta \u223cT\u22123/2$, *q _{min}* = 1.25, and there is no equilibrium rotation. The NIMROD results also have the viscosity kept fixed, with the Prandtl number approximately equal to 50, and the resistive wall time is set to $2\xd710\u22124$ s. This wall time is set to keep the ratio of tearing to wall time (the relevant timescale for the RPRW mode) for the dominant mode surface ($m/n=2/1$) of the simulation comparable to the ratio found in DIII-D plasmas with $S\u223c108$ and $\tau w\u223c2.4\xd710\u22123$ s. The results are shown in Fig. 14, where the growth rate of the least stable mode (as

*β*is increased) is shown using NIMROD for either a resistive or ideal DIII-D like wall. DCON was used to obtain the ideal plasma

*β*-limits in the figure, termed the resonant-ideal

*β*-limits in Sec. VI, which effectively sets the resonant $B\u22a5$ to zero on each respective rational surface. The DCON cases were run with and without an ideal wall (the ideal plasma no wall limit is the IPRW

*β*-limit in the case of zero rotation). The DCON simulations use the VACUUM code

^{67}for modeling the region outside the plasma, set with a conformal wall outside the last closed flux surface at a normalized minor radius of $r/a=1.4$. This wall distance is typical for matching the ideal limits with the actual DIII-D wall shape for equilibria like the ones used here, although it can vary with equilibrium details. For the equilibrium used in these simulations, the middle two

*β*-limits are very close (see Fig. 14), with the IPRW limit being slightly lower than the RPIW limit. Recall in Secs. IV B and V, that rotation stabilizes the RPRW mode up to the next ideal limit as

*β*is increased; this limit is either $\beta rp\u2212iw$, associated with the linear onset of a TM-dominated instability, or is $\beta ip\u2212rw$, associated with the linear onset of a RWM-dominated instability. The fact that the two limits IPRW and RPIW are so close suggests that the equilibrium for numerical simulation is near a point in parameter space where the limits would interchange, and thus, the mode behavior under rotational stabilization would change. This result is consistent with the effective triangularity $\delta \u223c0.35$ being close to and somewhat above the transition in Fig. 9 at $\delta \u223c0.29$. As the analytic form in Eqs. (1) and (2) cannot accurately capture this numerical boundary shape, the specific value of the effective triangularity is less important that the fact that it is clearly somewhat above 0.29, the transition point in the order of the IPRW and RPIW limits. Nevertheless, the qualitative behavior of stability, including the changing of the order between the

*β*limits RPIW and IPRW, is well captured by the sharp boundary model.

## VIII. CONCLUSIONS

A new sharp-boundary model has been developed to study the impact of shaping on the stability of MHD modes in a tokamak, including both plasma and wall resistivity. The model adapts an error field response model by Fitzpatrick^{46} to include (i) vertical asymmetry, (ii) toroidal curvature corrections, and (iii) a resistive wall surrounding the vacuum region outside the plasma. By scanning a broad parameter space, the model was used to examine the maximum extent to which rigid plasma rotation can raise the *β*-limit of the least stable mode, denoted $\beta rp\u2212rw$ for its combined resistive-plasma resistive-wall properties.

These modifications of Ref. 46 facilitate the derivation of a new resistive-plasma resistive-wall dispersion relation, for exploring the growth rates and stability limits over a variety of timescales; ideal, resistive-wall, and tearing. The sharp-boundary approach makes it efficient to scan stability over an individual parameter (such as *β*, safety factor, wall distance, elongation, triangularity or vertical asymmetry) while keeping the rest fixed. The qualitative study of a broad parameter space is designed to supplement and guide investigation by quantitative stability codes such as MARS,^{44} DCON,^{45} or NIMROD.^{47} By locating the sub-domains of interest, the model may be used to guide the investigation of larger codes in determining optimal design parameters for future tokamak devices.

The presented calculations demonstrate a new approach to shape optimization, based on maximizing the *β*-limit achievable with stabilization by rigid rotation or an equivalent feedback control with imaginary $\delta Br$-gain.

The analysis in this study is based on the 4-*β* formalism of Brennan and Finn,^{23} dividing the plasma response to rotation or feedback according to the four *β*-limits; resistive-plasma resistive-wall ($\beta rp\u2212rw$), resistive-plasma ideal-wall ($\beta rp\u2212iw$), ideal-plasma resistive-wall ($\beta ip\u2212rw$), and ideal-plasma ideal-wall ($\beta ip\u2212iw$), calculated *without* rotation or feedback. In the absence of rotation, increasing *β* causes the least stable mode to go unstable at the $\beta rp\u2212rw$ limit. To achieve higher *β*, the mode is shown to be stabilized by rotation (comparable with the wall-time or the tearing-time) up to the first ideal limit; this limit is either $\beta rp\u2212iw$, associated with the linear onset of a TM-dominated instability, or is $\beta ip\u2212rw$, associated with the linear onset of a RWM-dominated instability.

Extending existing predictions that the plasma-wall distance can interchange the order of the rp-iw and ip-rw *β*-limits,^{9,24} present results show that the safety factor, elongation, and triangularity can all introduce similar mode transitions, which affect the *β*-limit achievable in the presence of plasma rotation. The shaping window which maximizes the *β*-limit with rotational stabilization is found to reside around the transition point between $\beta ip\u2212rw$ and $\beta rp\u2212iw$, which defines a local optimum in the parameter space. Beyond this point, excessive elongation or triangularity is found to severely reduce the *β*-limit achievable with rotation by strongly reducing the $\beta ip\u2212rw$ limit (associated with RWM-dominated behavior) below the $\beta rp\u2212iw$ limit (associated with TM-dominated behavior).

## ACKNOWLEDGMENTS

### APPENDIX A: GRAD-SHAFRANOV SOLUTION FOR A SHARP-BOUNDARY MODEL

Let $(R,\varphi ,Z)$ be the standard right-handed cylindrical polar coordinates. The generalized cross-sectional coordinates, $(r,\theta )$, are mutually perpendicular to the axisymmetric coordinate $\varphi $. Following Connor and Hastie,^{55} (*R*, *Z*) are expressed as a Fourier series over $(r,\theta )$. The present formulation extends that of Ref. 55 to include the vertical asymmetry, by including both symmetric and asymmetric shaping parameters, respectively, labeled *S _{n}* and

*C*. The signs preceding the shaping parameters are chosen to guarantee the orthogonality of the external contravariant basis vectors, as in Eq. (A6). The radial shape functions $Sn(r)$ and $Cn(r)$ are determined by an $O(\u03f5)$ correction to the Grad-Shafranov pressure balance. The Grad-Shafranov equation is expanded up to $O(\u03f52)$, following the steps used by Fitzpatrick, Gimblett, and Hastie,

_{n}^{69}resulting in a shape equation for the vertically asymmetric terms $Cn(r)$, found to independently satisfy

This is the same equation as for the vertically symmetric components $Sn(r)$. Here, $f1(r)=\psi \u2032(r)$ describes the poloidal field distribution. The shape functions are resolved in the limit of a *δ*-function current at the plasma boundary. In order to permit a non-trivial Grad-Shafranov solution, a small flat-current region is maintained inside the plasma. The innermost surface *r *=* r*_{0} is kept finite in order to avoid numerical issues near *r* = 0.

A general flat current profile corresponds to a radial poloidal field distribution $f1(r)\u221dr$, according to Ampere's law. In the vacuum region beyond the plasma current, the poloidal field decreases as $f1(r)\u221dr\u22121$. These poloidal field functions result in a pair of shape equations independent of the equilibrium current magnitude and pressure. The regular solutions for the shape functions [according to Eq. (A1)] are given by

In the present model parametrization

we include the vertically symmetric terms representing elongation $e(r)=S2$ and triangularity $t(r)=\u2212S3$, and a vertically asymmetric term $d(r)=\u2212C3$ which introduces a separate top and bottom triangularity. For convenience, the non-constant part of the major radius is denoted $g(r,\theta )$, since a recurring factor of $R/R0=1+\u03f5g$ appears throughout the analysis.

Note that the coordinate parametrization in Eqs. (A4) and (A5) is designed so that the covariant derivatives satisfy the orthogonality relations

which guarantee the orthogonality of the basis vectors, $\u2207r\xb7\u2207\theta =0$. Orthogonality simplifies the external metric coefficient

required for the *θ*-line element $dl\theta =rhd\theta $, the surface area element $ds=rhRd\theta d\varphi $, and the volume element $d3x=rh2Rdrd\theta d\varphi $.

The radial structure required by Eqs. (A2) and (A3), which extends the original Fitzpatrick parametrization^{46} with vertical asymmetry, is given by

The finite innermost surface $r=r0>0$ is designated to prevent numerical issues, and to initialize integration through $r0<r<1$ with decoupled cylindrical-like solutions $Vm(r0)\u221dr0m$. Throughout this paper, we set $r0=0.3$. This set of non-orthogonal internal coordinates and orthogonal external coordinates is depicted in Fig. 15. The constants are related to the boundary shaping parameters in Eqs. (1) and (2) by

The resulting contravariant derivatives are obtained by the dual relations $\u2207r=J\u22121(\u2202Z\u2202\theta R\u0302\u2212\u2202R\u2202\theta Z\u0302)$ and $\u2207\theta =J\u22121(\u2202R\u2202rZ\u0302\u2212\u2202Z\u2202rR\u0302)$. Finally, the toroidal component satisfies $|\u2207\varphi |=R\u22121(r,\theta )$.

The field perturbations throughout the volume must satisfy Laplace's equation, neglecting the axisymmetric $O(\u03f52)$ terms. In the internal region, ($r0<r<1$), this equation retains the form derived by Fitzpatrick in Ref. 46 Eqs. (48)–(56), but with an additional toroidal factor $1+\u03f5g(r,\theta )$ multiplying each of the metric coefficients.

In the external region, $1<r$, including the toroidal curvature term neglected in Ref. 46 modifies the form of Laplace's equation to

in terms of the vacuum region coupling matrix, which can be expressed analytically as

The Kronecker-*δ* matrices $\delta m,k\xb11$ and $\delta m,k\xb12$ show how elongation introduces additional $m\xb11$ coupling, and triangularity and vertical asymmetry introduce $m\xb12$ coupling.

### APPENDIX B: FITZPATRICK MODEL TEARING LAYER RESPONSE

The tearing layer formulation employed by Fitzpatrick in Ref. 46, is adapted to incorporate the effects of toroidal curvature in the outer region. While the model is constructed to accept any tearing layer regime, present calculations apply a constant-*ψ* visco-resistive^{53,54} (VR) layer condition. Each tearing surface *r _{s}* imposes a boundary condition similar to the resistive-wall condition in Eq. (6), except applied separately to each resonant harmonic. In addition, the growth rate in the tearing equation, unlike that appearing in the wall equation, is Doppler shifted due to plasma rotation. Fitzpatrick's straight field-line angle is modified slightly by the toroidal curvature factor $1+\u03f5g(r,\theta )$, and becomes

Here, the subscript *m* denotes evaluation at the resonant *m*/*n*-surface, with fixed *n* = 1. The equilibrium poloidal field $B\u0302m=B\u0302\theta (r=rm)$ can be approximated by the same functional form at *r _{m}* as at $r=1+$, with the resonant-surface poloidal field distribution solved analogously to Eqs. (4) and (5), except for $qm=m/n$ leading to a different free parameter to replace

*B*.

_{i}A modified tearing surface coupling matrix

presently yields the resonant layer response

for each rational surface *r _{m}*. Each has a corresponding Doppler-shifted growth rate $\gamma Dm=\gamma +in\Omega m$, layer time

*τ*, and associated geometric time factor

_{m}analogous to the wall matrix in Eq. (10) corresponding to a multi-harmonic geometric time factor. The formulation of Ref. 46 neglects this tearing layer integral, effectively applying a cylindrical low-*β* approximation. Noting the typical tokamak ordering of the poloidal field $B\u0302m=B\u0302\theta (r=rm)\u223c\u03f5$, the integrand above appears as $O(\u03f53)$ but is actually $O(\u03f5)$.

Following Ref. 46, the layer response is simplified in the limit where the rational surfaces all lie just outside of the plasma edge for numerical convenience in the outer region. In this way, the resonant surfaces are combined by taking $rm\u21921+,\u2009\tau m\u2192\tau t,\u2009\Omega m\u2192\Omega ,\u2009hm(\theta )\u2192ha(\theta )$, and $gm(\theta )\u2192ga(\theta )$. Artificial singularities at the plasma boundary are avoided by calculating the equilibrium field distribution $B\u0302m(\theta )$ at the resonant layer with $qm\u2260qa$. With a unique poloidal field distribution, each layer maintains its own geometric coupling matrix $U(tear)$ and geometric time factor $\tau \u0302(tear)$. After taking the limit and summing over all rational surfaces, $q0<qm<qa$, the total resonant response at $rm\u2192rs\u21921+$ is written as

in terms of the tearing layer response matrix

which combines all of the relevant rational surfaces in the range $nq0<m0\u2264\u2026\u2264m1<nqa$. This tearing formulation agrees with that of Ref. 46 for the case of high aspect-ratio and vertical symmetry. This tearing response matrix enters the model dispersion relation via the perturbed field relation in Eq. (18).

In the language of the tearing theory, the inner solution is given by a constant-*ψ* visco-resistive^{53,54} diagonal matrix $\Delta mk(\gamma D\tau t)=\gamma D\tau t\delta mk$ (whereas the other terms in Eq. (B5) all belong to the outer solution $\Delta \u2032$). Other tearing regimes may be conveniently substituted into this formulation by replacing the linear $\gamma D\tau t$ with different functional forms (for example, $(\gamma D\tau t)5/4$ in the resistive-inertial regime), while keeping the outer geometric coupling terms unchanged. This is yet another advantage of the sharp boundary model construction presented in this paper.