Analytic theory of ideal-MHD vertical displacements in tokamak plasmas

An analytic derivation of the relevant dispersion relation for vertical displacements in shaped tokamak plasmas is presented, valid for arbitrary values of the ellipticity parameter. The theory is developed within the framework of the reduced ideal-MHD model. A nearby, perfectly conducting wall can provide passive feedback stabilization of vertical displacements on the ideal-MHD timescale. The mechanism for passive stabilization relies on image currents induced on the metallic wall. However, if the plasma extends to the magnetic separatrix, where magnetic X-points are located, as in the case of a divertor tokamak configuration, perturbed axisymmetric currents carried by the plasma in the vicinity of the X-points are triggered. It is shown that these X-point currents can provide passive feedback stabilization, even in the absence of a nearby wall. X-point currents are excited due to the resonant nature of magnetic X-points with respect to toroidal axisymmetric perturbations. An intermediate case, where the plasma boundary is located just inside the magnetic separatrix, is also analyzed, providing additional insight into the stabilization mechanism.


Introduction
Present-day tokamak experiments adopt plasma shaping and magnetic divertor configurations in order to optimize fusion performance, and also to reduce the adverse effects of plasma-wall interactions, as experimental and theoretical studies have confirmed [1][2][3][4][5][6][7][8][9]. Plasma shaping, and in particular cross-section elongation, requires careful design of the * Authors to whom any correspondence should be addressed.
Original Content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. plasma-facing components and feedback control by currents flowing in external coils in order to avoid vertical displacement events, which may endanger the machine's integrity and safe operation [10][11][12][13][14][15][16]. Hence, plasma vertical stability has been the subject of several analytic and numerical investigations, of which [17][18][19][20][21][22][23][24][25][26][27][28][29][30] are a sample.
Typically, in these investigations, the plasma is modeled by ideal-MHD, with boundary conditions representing nearby ideal or resistive wall and external feedback currents. The basic idea is that eddy currents flowing along the surface of the wall, and/or along ad hoc metallic conductors facing the plasma, can lead to passive stabilization of the ideal-MHD vertical instability. Passive stabilization is essential, because if somehow vertical displacements were allowed to grow on the Alfvénic timescale, magnetic fluxes generated by external currents placed outside the wall of the vacuum chamber for active feedback control would not have the time to penetrate the wall and counter the fast-growing ideal-MHD instability. Once passive stabilization of the ideal-MHD vertical instability is achieved, vertical displacements may still grow on the resistive wall timescale. For this relatively slow growth, active stabilization by means of feedback currents flowing in external coils becomes feasible.
Magnetic divertor configurations are also an important feature of present-day tokamak experiments [1,2]. For these configurations, the last closed flux surface is a magnetic separatrix with, typically, either one (single-null) or two (double-null) Xpoint(s) (more complex configurations with multiple X-points have also been explored [31]). These are points on the plasma poloidal cross-section, or magnetic field lines when viewed in three dimensions, where the poloidal magnetic field vanishes and the confining magnetic field is purely toroidal. Vertical displacements, on the other hand, are axisymmetric perturbations with toroidal mode number n = 0. Therefore, a vertical displacement perturbation is resonant at the X-point(s) of a magnetic divertor configuration, in the sense that this perturbation, regardless of its poloidal modulation, is constant along the toroidal field line going through the X-point. Mathematically, the resonance condition is represented by the criterion B eq · ∇χ = 0, where χ is a generic axisymmetric perturbation and B eq the equilibrium magnetic field. This criterion is trivially fulfilled at the X-points, since the toroidal mode number is n = 0, and the poloidal component of the equilibrium magnetic field, B p , vanishes at such points.
As a consequence of the X-point resonance, axisymmetric current sheets localized along the separatrix are likely to form. This process, which leads to current sheet formation in the proximity of magnetic X-points, has been studied in the context of astrophysical plasmas (see, e.g. [32,33]), as well as in connection with fundamental laboratory plasma experiments (see, e.g. [34]), and is well known to researchers working in magnetic fusion (e.g. [28,[35][36][37]). However, as far as we are aware, [38,39] are the first articles where the resonant interaction between n = 0 modes and the X-point(s) of a magnetic divertor configuration was addressed analytically within the framework of a tokamak plasma. Experimental evidence of n = 0 perturbations and current sheets at X-points have been observed in tokamak experiments such as the Joint European Torus (JET) [40][41][42][43]. Current sheets have also been observed in numerical simulations of the vertical instability in tokamaks, with advanced numerical codes that are able to treat correctly the X-point geometry, such as M3D-C 1 , NIMROD and JOREK [10,30,44]. However, analytic understanding of why these current sheets form, and more importantly, the impact they have on the stability of vertical displacements, was not clarified in those numerical works. A caveat of our analysis is that it is restricted to ideal-MHD, while real plasmas are resistive, and the cited numerical works took plasma resistivity into account. Thus, while current sheet formation near magnetic Xpoints is a robust prediction that can be made within the framework of ideal-MHD, our work also suggests that plasma electrical resistivity in a narrow boundary layer along the magnetic separatrix, in addition to wall resistivity, may have a profound impact on the stability of n = 0 vertical displacements.
The main result of [39] is that the current sheets that form along the tokamak divertor separatrix are capable of suppressing the vertical instability on the ideal-MHD timescale, thus providing effective passive stabilization of these modes, even in the absence of a nearby wall and/or ad hoc plasma-facing components. The purpose of this article is to extend the analysis of [39], by providing a more detailed analytic derivation, and to shed additional light on the mechanism leading to the ideal-MHD stabilization of vertical modes as a consequence of the X-point resonance.
This article is organized as follows. In section 2, an heuristic derivation of ideal-MHD vertical stability theory is presented, assuming a rigid vertical displacement of the plasma and in the absence of magnetic X-point effects. Section 3 outlines the model equations and the initial equilibrium configuration. In section 4, we consider the situation where the plasma terminates on an elliptical boundary well within the magnetic separatrix containing the X-points. This configuration, which corresponds to the one first treated in the pioneering work of [17], leads to the conclusion that vertical displacements in elongated tokamak plasmas are ideal-MHD unstable, unless a conducting wall is located relatively close to the elliptical plasma boundary. Therefore, this situation is more appropriately referred to as the limiter tokamak scenario. In section 5, the plasma is assumed to extend to the divertor magnetic separatrix and the impact of the X-point resonance is analyzed. It is found that vertical displacements are stable on the ideal-MHD timescale thanks to current sheets that form along the magnetic separatrix, even in the absence of a wall. Section 6 discusses the intermediate case, where the wall is absent and the plasma terminates close to the magnetic separatrix, however leaving the X-points in the vacuum region. This scenario is studied in order to understand, from a mathematical point of view, how the condition for ideal-MHD marginal stability for the no-wall case can be obtained. Conclusions are presented in section 7. In order to improve the readability of the main text, many mathematical details are provided in the appendixes. Units adopted in this article are c.g.s.

Heuristics
In this section we derive, in an heuristic way, the relevant dispersion relation for a vertical displacement perturbation, assuming the latter to correspond to a rigid shift of the currentcarrying plasma, and on the basis of a simplified toy model. The impact of the X-points is not considered in this section, but it will be discussed in section 5.
Let us consider the equilibrium of three rectilinear current wires, as shown in figure 1(a). At equilibrium, the plasma wire carrying the current I P is located at y = 0, while two equal external wires, each carrying a current I Ext , are located at y = ±l. Note that we take the y-axis as the vertical direction. The currents flow along the z-direction, which mimics the toroidal direction of a tokamak configuration. All three currents are taken to be positive. The external wires are fixed in space, while the plasma wire is free to move along the ydirection. Figure 1(a) also sketches the magnetic flux surfaces produced by the three current wires. A vacuum is assumed to surround the three wires, and therefore, no current is allowed to flow through the magnetic X-points. This configuration is taken to mimic, at a very essential level, the magnetic doublenull equilibrium of a straight tokamak bounded by a magnetic divertor separatrix. Clearly, as we have already pointed out, the more realistic case where the plasma extends to the magnetic separatrix and currents are allowed to flow through the X-points obviously cannot be treated by this toy model, and will be discussed in section 5.
The equation of motion for the plasma wire is (in c.g.s. units): where µ is the linear mass density, c is the speed of light, and an over-dot signifies the time derivative. We neglect self and mutual induction currents. Thus, I P and I Ext remain constant as the plasma wire is displaced. For small y ≪ l, the solution of equation (1) is: where y 0 is an initial displacement, and is an inertial growth rate. At this stage, no particular relation exists between the currents I P and I Ext and the distance l. However, if instead of a plasma wire we consider a diffused plasma column with a uniform current density extending, on the Oxy cross-section, up to an elliptical magnetic surface with minor semi-axis a and major semi-axis b contained within the magnetic separatrix, then, according to the equilibrium solution in [45], a relationship is established: When equation (4) is used in equation (3), the exponential growth rate of the vertical displacement becomes: Therefore, the relevant growth rate depends only on the plasma current I P and on the parameters a and b, but does not depend on the distance l, nor on the value of the external currents. We can rewrite γ H in terms of more familiar plasma parameters. For a cylindrical plasma column with an elliptical cross-section, the linear mass density µ should be replaced by µ → π a b ϱ m , where ϱ m is the volume mass density. We introduce the ellipticity parameter: If we assume e 0 to be small, then a ≈ b, 2 I P ≈ c a B p (a), and the growth rate takes the form: where v A = B p (a)/(4π ϱ m ) 1/2 is the familiar Alfvén velocity based on the poloidal magnetic field. The inverse of γ H is of the order of the characteristic Alfvén time. This is indeed a very fast growth time. For instance, if we take a hydrogen plasma with e 0 = 0.2, a = 1 m, B p (a) = 1 T, and number density n = 10 20 m −3 , we find γ −1 H ≈ 1 µs. The ideal vertical instability is suppressed in the presence of a perfectly conducting wall, provided this wall is not too far away from the plasma. We can understand the stabilization mechanism in the following way. When the plasma current is displaced from its equilibrium position, image currents are induced at the wall, as shown schematically in figure 1(b). The sign of these currents is such that the corresponding forces oppose the motion of the plasma wire. From an heuristic point of view, this effect can be conveniently mimicked by two currents of opposite sign, ±δI, carried by the two external wires, and added to the external currents I Ext (see figure 1(c)). The currents ±δI can be thought of as driven by an induced e.m.f. proportional toẏ. In the perfectly conducting limit, where D is a dimensionless proportionality constant. In the actual tokamak problem, the parameter D depends on the wall geometry (and on other passive stabilization components eventually placed inside the vacuum chamber). Including the effect of the perturbed external currents ±δI, the equation of motion for the plasma wire becomes: where we have taken y ≪ l in the second of equations (9). Stability is obtained for: In this limit, an initial displacement of the plasma wire from its equilibrium position results in oscillatory motion with a characteristic frequency: In [17] (and in section 4 of this article), the geometrical factor D is determined analytically for the case of confocal elliptical plasma and wall cross-sections, and a uniform equilibrium current density distribution up to the elliptical plasma boundary. The following expression is found: where a, b, a w and b w are the minor and major semi-axes of the confocal plasma and wall elliptical cross-sections, respectively, with b 2 − a 2 = b 2 w − a 2 w . According to the equilibrium adopted in [17] and in this article, the marginal stability condition, D = 1, corresponds to the situation where the wall elliptical cross-section intercepts the magnetic X-points; values of D > 1 that are necessary for passive wall stabilization are such that the magnetic X-points lie outside the vacuum chamber.
We also present an alternative derivation of the dispersion equation (9) that relaxes the assumption that the plasma current is simply a filament, but instead makes full use of the equilibrium solution for the equilibrium flux function ψ eq (x, y) derived in [45] (see also the next section). We assume once again that the relevant perturbation corresponds to a rigid shift of the entire plasma column. Suppose that the plasma displaces rigidly in the y-direction by a distance ξ. In the plasma frame, this is equivalent to displacing the two external current filaments in the negative y-direction by the same distance. Suppose that the current in the filament at y = l − ξ is I Ext − δI, whereas the current in the filament at −y = l − ξ is I Ext + δI. The net force per unit length acting on the filaments due to the plasma is: where ψ P (x, y) is the magnetic flux generated by the currents flowing in the plasma. Following the derivation in appendix A, we obtain an equation forξ = ξ/l, where we have normalized δÎ = δI/I Ext andt = γ H t, with γ H as given in equation (7). Thus, in the absence of feedback, δÎ = 0, the result in equation (2) is recovered. With feedback, if we set δÎ = Dξ as in equation (8), the stabilization criterion D > 1 is also recovered. As we shall see, the value of γ H obtained heuristically in this section, as well as the passive stabilization criterion (10), with the parameter D given by equation (12), agree perfectly well with those obtained by the more detailed normal mode analysis of section 4, where the plasma density drops to zero outside the elliptical boundary. This result is no longer valid if the plasma density remains finite around the magnetic Xpoints, as shown in section 5 and first discussed in [38,39].

Model equations, plasma equilibrium, linearized model
The model adopted in this work is the standard reduced ideal-MHD model (see, e.g. [46]), where the reduction is based on the low-β tokamak ordering (β = plasma pressure/magnetic pressure). The magnetic field is represented as B = e z × ∇ψ + B z e z , where e z is the unit vector along the z-direction, and B z is assumed to be nearly constant. We also assume that all physical quantities are independent of the z coordinate. The plasma flow is represented as v = e z × ∇φ.
The reduced ideal-MHD model equations in dimensionless form, for the magnetic flux function ψ, the stream function φ, and non-constant mass density ϱ, are: where the brackets are defined as [χ, η] = e z · ∇χ × ∇η, with χ and η two generic scalar fields, J = ∇ 2 ψ is the normalized current density, and U = ∇ 2 φ is the normalized flow vorticity. Space and time are normalized asr = r/r 0 , where is a convenient equilibrium scale length, andt = t/τ A , where τ A = (4 πρ m ) 1/2 /B ′ p is the relevant Alfvén time and B ′ p the radial derivative of the poloidal magnetic field at the magnetic axis. The dimensionless fields are normalized asψ = ψ/(B ′ p r 2 0 ),φ = (τ A /r 0 2 )φ; the normalized plasma density isρ = ϱ m /ϱ m0 , with ϱ m0 the density on the magnetic axis. In order to simplify the notation, over-hats are dropped in equations (15) and (16), and in the following.
At equilibrium, the fields are independent of time, and equilibrium plasma flows are assumed to be absent. Thus, (15) and (16) reduce to: where the subscript 'eq' denotes equilibrium quantities. The general solution of (18) satisfies: As discussed in section 2, the occurrence of the vertical instability does not depend on details of the equilibrium current density profile, J eq , but on finite ellipticity, e 0 , and on the total current carried by plasma, I P . Therefore, in order to allow for analytic work, a flat equilibrium current density profile can be chosen, as in [17,38], where the total plasma current is uniformly distributed within a region delimited by a convenient elliptical flux surface with minor and major semiaxes a and b, respectively. The area of this region is S = πab, hence I P = πabJ eq . We introduce elliptical coordinates (µ, θ), where: In section 4, the elliptical surface µ = µ b will correspond to the actual plasma boundary. However, in sections 5 and 6, the plasma will extend beyond this convenient flux surface. Thus, consistent with the normalization introduced above, we take is the Heaviside unit step function. Since the ellipse µ = µ b is a magnetic flux surface (this is not true for all other µ = const ellipses), it can also be stated that the equilibrium current density is a function of the equilibrium flux, ψ eq ; hence, (19) is satisfied. In other words, we can also take J eq = 2H(ψ b − ψ eq ), where ψ b is the value of the equilibrium flux at the plasma boundary.
Based on the analysis of [45], it can be shown that the solution, ψ eq , of the equilibrium problem, in the limit where the parameter ε = [(a 2 + b 2 )/l 2 ]e 0 is small, reduces to Gajewski's solution [47], represented as follows. For µ < µ b , the solution of ∇ 2 ψ − eq = 2 that satisfies the regularity condition at the magnetic axis, and that reduces to a constant at the flux surface µ = µ b , in Cartesian coordinates is given by: Here, and in the remainder of this article, we take the subscript '−' to indicate the representation of scalar fields in the region inside the convenient elliptical surface. In this region, the equilibrium magnetic flux surfaces are nested ellipses with constant elongation, κ = b/a. As we stated above, confocal ellipses with µ = const are not flux surfaces, except for the special one corresponding to For µ > µ b , the solution of the vacuum equation ∇ 2 ψ + eq = 0, subject to the boundary conditions that the equilibrium magnetic flux ψ eq (µ, θ) and its normal derivative ∂ψ eq /∂n = n · ∇ψ eq be continuous across the surface µ = µ b , is best obtained in elliptical coordinates: with α 2 = ab/r 2 0 = (1 − e 2 0 ) −1/2 and e 0 the ellipticity parameter defined in equation (6). The subscript '+' indicates the representation of scalar fields in the region outside the convenient elliptical surface. The special flux surface defined by Gajewski's solution approximates well the complete solution given in [45] up to values of |x| and |y| much smaller than l, which, for ε ≪ 1, include the separatrix region.
In the following sections, we shall seek analytic solutions of equations (23) and (24).

Normal mode analysis: limiter tokamak scenario
In this section, we consider the convenient elliptical surface µ = µ b to actually correspond to the plasma boundary. Therefore, the X-points are in the vacuum region, where µ > µ b . We assume a constant plasma density profile up to the elliptical boundary at µ = µ b , and zero outside that boundary: ϱ eq = H(µ b − µ), with H(x) the unit step function. As shown in figure 2, we denote the plasma region by the symbol Ω, while vacuum occupies region ∆ and the V-regions (the meanings of regions V ± and V ′ ± will be discussed in the next section). Let us assume for a moment that, in region Ω (with the exclusion of its boundary), the inertial term in equation (24) is negligible, either because the growth rate γ 2 is small or, more appropriately for the calculation that follows, because the vorticity ∇ 2φ vanishes. Therefore, the perturbed flux satisfies [ψ eq ,J] = 0. There are two possibilities: eitherJ =J(ψ eq ), orJ = 0. We discard the first possibility, as it gives rise to a spurious additional part of the equilibrium current density. Therefore, using elliptical coordinates, we conclude that inside region Ω, is the metric element for elliptical coordinates, and we have used the short-hand notation ∂ µ = ∂/∂µ (and similar).
The general solution forψ − can be written as a summation, over integer elliptical mode numbers m, of products of exponential functions of mµ and sinusoidal functions of mθ. However, as discussed in [17], on the basis of the ideal-MHD energy principle, the most unstable solution corresponds to elliptical mode number m = 1; furthermore, it is odd in θ and even in µ. Thus, the solution of interest is: where ψ 1 is a constant amplitude. The flux-freezing equation (23), in Cartesian coordinates, is: Consistent with the solution forψ − , we find: where we have introduced the constant displacement amplitude ξ, and ψ 1 = −ξ/b 2 . As we can see, solution (27) corresponds to a rigid vertical displacement of constant amplitude ξ, and the perturbed flux can also be written as: Also note that the vorticityŨ = ∇ 2φ = 0, as we anticipated. An interesting feature of the solution we have obtained is that bothψ andφ involve a single m = 1 harmonic in the elliptic angle θ, while, if polar coordinates were used instead, a large number of poloidal harmonics would be coupled together.
In the vacuum region, where µ > µ b , the perturbed flux satisfies ∇ 2ψ+ = 0. Assuming, at first, that the ideal wall is placed at infinity, the relevant solution that decays to zero at infinity, and such thatψ + =ψ − at the elliptical boundary µ = µ b , is:ψ The perturbed flux is continuous across the elliptical boundary, but its derivative is discontinuous, giving rise to a current sheet, where δ(x) is the Dirac delta function. A short calculation provides the modulation of the current sheet with respect to the elliptical angle: where we have used h 2 (µ, θ) = A 2 (cosh 2µ + cos 2θ)/2. The final step is to obtain the dispersion relation for γ 2 . Let us consider the perturbed equation of motion, equation (24). Each of the three terms in this equation is proportional to a delta function centered at µ = µ b . We point out that the bracket terms appearing on the r.h.s. can each be expressed in terms of a divergence; for instance, [ψ, J eq ] = ∇ · [(e z × ∇ψ)J eq ]. An expedient way to obtain the value of γ 2 is to multiply equation (24) by the scale factor h 2 , and integrate across the elliptical boundary over a narrow layer of infinitesimal width, The three terms in (32) are evaluated separately as follows: Balancing the three terms: Note that each of the three terms is proportional to cos θ. Substituting the explicit expressions ofψ andφ, we obtain the dispersion relation, where, for the sake of clarity, physical dimensions for the growth rate, γ, and for the plasma semi-axes, a and b, have been reintroduced. In terms of the ellipticity parameter, e 0 , equation (37) can also be written as: We point out that this result is valid for arbitrary values of e 0 in the interval 0 ⩽ e 0 ⩽ 1. In the limit e 0 = 0, where the plasma boundary becomes circular, the growth rate vanishes, as it should. Neutral stability is obtained also for e 0 = 1, corresponding to the case where the plasma cross-section has an infinite elongation, but also a vanishing area, for a fixed value of a or b. Instability is found for any positive value of e 0 < 1. In the limit of small ellipticity, the growth rate reduces to γ ≈ e 1/2 0 τ −1 A , in full agreement with the heuristic result obtained in equation (7).
If a perfectly conducting wall is present, passive feedback stabilization of the vertical mode can be obtained, depending on the distance between the plasma boundary and the wall. As was done in [17], we assume that the wall is also elliptical, with coordinate µ = µ w ⩾ µ b (if µ w = µ b , the wall coincides with the plasma boundary). Thus, the wall and plasma boundary are confocal ellipses, and therefore b 2 − a 2 = b 2 w − a 2 w , where b w and a w are, respectively, the major and minor semi-axes of the elliptical wall cross-section.
The dispersion relation with passive feedback stabilization can be obtained as follows. The solutions (27) and (28) in the Ω-region are unaffected by the presence of the wall, except, as we shall see, for a redefinition of the displacement amplitude, ξ. On the other hand, the vacuum solution forψ is modified, as the perturbed flux must vanish on the perfectly conducting wall. In the vacuum region, the solution of interest isψ with c 0 and c 1 integration constants to be determined by the conditions that the perturbed flux is continuous at µ = µ b and vanishes at µ = µ w . In a physically more meaningful way, we prefer to represent the vacuum solution as: where ξ ∞ is the amplitude of the rigid vertical displacement in the limit where the wall is moved to infinity, and the term proportional to ξ ext represents the contribution to the perturbed flux due to the image currents that form on the wall when this is at a finite distance from the plasma boundary. Continuity of flux at the plasma boundary requires that ξ = ξ ∞ − ξ ext , and so the actual vertical displacement ξ is reduced, as compared with the no-wall case, by the amount ξ ext . The conditioñ ψ + (µ w , θ) = 0 implies that: Let us now replace the wall-modified vacuum solution in the expression (30) for the current sheet. Instead of equation (30), we obtain: Finally, using equations (39)- (41) in equation (36), we obtain the dispersion relation modified by the presence of the ideal wall: where and γ ∞ is the growth rate (38) found in the limit D → 0, where the wall is at infinity. A short calculation (see appendix B) shows that expression (43) for the parameter D is exactly equivalent to that in equation (12). Therefore, apart from an extra term in the denominator of equation (42) proportional to e 0 , which can be neglected in the limit of small ellipticity, the dispersion relation (42) agrees perfectly well with equation (11) obtained heuristically in section 2 when ellipticity is small. Equation (42) indicates that passive feedback stabilization is obtained for values of D > 1. The maximum D value is found when the wall coincides with the plasma boundary. As a numerical example, take a w = a, b w = b and b/a = 1.4. In this case, D = 5.3. Note that in this limit, where µ w → µ b , the vertical mode is purely oscillatory, with a frequency ω = ±(D max − 1) 1/2 γ ∞ , but also, the amplitude of the displacement, ξ = ξ ∞ − ξ ext , goes to zero, as ξ ext → ξ ∞ . As the wall is placed further away from the plasma boundary, the value of D decreases monotonically and the vertical mode remains purely oscillatory for as long as D remains larger than unity. The marginal stability value, D = 1, is obtained when µ w = 2µ b , which corresponds to the case where the wall intercepts the X-points, as shown in figure 3 and in appendix B. Values of D < 1, for which no passive feedback stabilization is possible, are found when the X-points lie inside the volume bounded by the wall, i.e. inside the tokamak vacuum chamber, while the plasma boundary is well inside the magnetic separatrix. We point out that this result is in perfect agreement with that obtained by Laval et al [17] on the basis of the ideal-MHD energy principle. We note that the conclusion, that passive wall stabilization is possible only if the X-points are located outside the vacuum chamber, is valid for the case where the plasma boundary and the wall are confocal ellipses, but may not be true when the shape of the wall follows more closely the shape of the plasma boundary. On the other hand, the dispersion relation for ideal-MHD vertical displacements given by equation (42), valid in the 'limiter tokamak scenario' for arbitrary values of the ellipticity parameter e 0 , together with its detailed normal mode derivation, are shown in this article for the first time.

Normal mode analysis: divertor tokamak scenario
Next, we consider the case where the plasma is bounded by the magnetic divertor separatrix, which therefore represents the last closed flux surface. In this case, the plasma density extends to the separatrix, and axisymmetric modes are resonant at the X-points. As a consequence, perturbed current sheets are driven in the vicinity of the X-points and along the magnetic separatrix [39]. In this article, we extend the analysis of [39] by providing a more detailed mathematical derivation and additional information.
We remark that a current sheet always forms at the plasma boundary. Indeed, this was the case for the limiter tokamak scenario treated in the previous section. However, if the plasma extends to the magnetic separatrix, the effect of the current sheet at the plasma boundary changes from destabilizing, as for the limiter case, to stabilizing, as was shown in [39] and detailed in the analysis that follows.
We assume a uniform density profile, up to the magnetic separatrix, where the density drops rapidly to zero: ϱ eq (ψ eq ) = H(ψ X − ψ eq ). The equilibrium current density, on the other hand, is the same as that adopted in the previous section, i.e. J eq = 2H(ψ b − ψ eq ), where ψ eq = ψ b = 1/2 corresponds to the convenient elliptical surface, µ = µ b , which now should not be confused with the plasma boundary at ψ eq = ψ X .
In order to obtain the relevant dispersion relation, we follow a procedure similar to the one adopted in the previous section. In particular, the rigid-shift solution for the stream function,φ, given by equation (27), is now valid all the way to the separatrix. Consequently, the perturbed flux in region Ω of figure 2 is still given by equation (28), while, in region ∆, it takes a different form, since the expressions for the equilibrium flux, ψ eq , are different in the two regions. In order to obtainψ ∆ , we use the flux-freezing condition (23), expressed in elliptical coordinates: After a straightforward calculation (see appendix C for details), we obtain: Thus, the perturbed flux in region ∆ is a solution of the Laplace equation,J ∆ = ∇ 2ψ ∆ = 0, corresponding to a current-free perturbation. It still involves a single m = 1 harmonic in the elliptical angle θ. The important aspect is that the solution we have obtained satisfies the ideal-MHD constraint at the Xpoints,ψ(µ X , θ X ) = 0. Moreover,ψ is continuous at the elliptical surface, but its first derivative is discontinuous, and so a perturbed current sheet is still present there: the term [ψ eq ,J] in equation (24) generates a delta function, δ(µ − µ b ), at the elliptical surface, which, however, is exactly canceled by the other delta-function term arising from [ψ, J eq ]. Indeed, In the latter equation, the semi-axes a and b are normalized to r 0 , therefore, (a 2 + b 2 )/(a 2 b 2 ) = 2.
In the vacuum regions, denoted by the letter V in figure 2, the perturbed fluxψ V also satisfies the Laplace equation, whose general solution with the appropriate parity in the angle θ, and decaying to zero at infinity (we consider here only the no-wall case), is: In principle, the coefficients c m are determined by continuity with the perturbed fluxψ ∆ at the magnetic separatrix. However, this turns out to be an analytically difficult task (see also the next section): the plasma boundary, ψ eq = ψ X , is not a constant-µ coordinate surface, and a simple analytic expression for µ as a function of θ along the separatrix cannot be obtained. As a consequence, an infinite number of harmonics in the elliptical angle, θ, are coupled on that surface and into the vacuum solution. This difficulty is resolved ifψ in regions ∆ and V ± is expressed in flux coordinates: where ψ + eq is given by equation (22). Since ∇ 2 u = 0 and ∂ θ u = −∂ µ v, coordinates (u, v) are harmonic and orthogonal in regions ∆ and V ± (but not in region Ω, and so we will not use them in that region). The separatrix corresponds to u = 0 and, at the X-points, v = 0, ±π. On the elliptical boundary, u = u b = −µ b + (e 0 /2) sinh(2µ b ) and v = θ − π/2 + (e 0 /2) sin(2 θ). With reference to figure 2, u is negative in regions ∆, V − and V ′ − , and positive in regions V + and V ′ + ; v ranges from −∞ to +∞.
The summation involves odd integers only. Note that, at the X-points,ψ ∆ (u = 0, v = 0, π) ∝ α m = 0. The coefficients α m and β m can be determined as follows. We exploit the fact that the convenient elliptical surface between regions Ω and ∆ is special, in that both µ = µ b and u = u b are constant on that surface. Then, we setψ where use has been made of (∂ θ u)| µ b = 0. When made explicit, these relations provide a set of two equations for the determination of the coefficients in (50): and J ν (x) are Bessel functions of integer order (since m is odd). One of the essential mathematical aspects of this solution is that the perturbed magnetic flux develops a singularity at the X-points. Indeed, let us consider the perturbed flux on the separatrix, u = 0, as a function of v: The leading asymptotic behavior of the Fourier coefficients for large values of m (see, e.g. [48]) is: where is a positive constant. Thus, the Fourier spectrum for the perturbed flux along the separatrix never decays exponentially for m → ∞, which is indicative of singular behavior. Indeed, consistent with α m ∼ p/m 3/2 for large m, it follows that (see appendix E): This behavior can also be found by direct expansion ofψ ∆ for u = 0 and small v, using equations (48) and (49). As a result, more and more Fourier harmonics are needed in order to properly represent the perturbed flux in the vicinity of the Xpoints. This is an important difference with respect to the case we treated in the previous section, where only one Fourier harmonic in the elliptical angle θ was sufficient. From a mathematical point of view, we may argue that this explains the change in behavior, from unstable to stable (as we shall see shortly), of the rigid shift vertical displacement, when the plasma extends to the magnetic separatrix (and the no-wall case is considered). Next, we are going to construct the vacuum solution in regions V. In V + , the flux coordinate, u, is positive and not bounded, The vacuum solutionψ V+ is no longer periodic in v, and its relevant symmetry is odd with respect to reflections π/2 − v ↔ π/2 + v. The general solution to the Laplace equation that is consistent with the above conditions and is regular at infinity is:ψ Likewise, in V − , u is negative and not bounded, −∞ < u < 0, and v is not cyclic, −∞ < v < +∞. The vacuum solution, ψ V− , is not periodic in v, and its relevant symmetry is even with respect to reflections v ↔ −v. The general solution to the Laplace equation that is consistent with the conditions considered above and is regular at infinity is given as: As for the solutions in V ′ + and V ′ − , symmetry consideration lead to: where the ± sign inψ V ′ − are taken on the positive and negative sides of the branch cut along the negative direction of y-axis.
Imposing flux continuity at the separatrix, from the V + side, yields: and from the V ′ + side, Clearly, the two conditions are equivalent, as they should be, under the transformation Applying an inverse sine Fourier transform leads to the following relation between α(t), β(t) and α m (see appendix F): In order to determine the functions α(t) and β(t), we observe that, because of continuity along the magnetic separatrix, the perturbed flux approaching the X-points from region V + and from region ∆ must have the same behavior as function of the v coordinate. Investigation of (67) indicates that consistency with the continuity condition requires (see appendix G): Therefore, the asymptotic behavior of the perturbed flux along the separatrix, approaching the X-points from regions V + and V − , can be expressed as: ∂ψ Similarly, Another mathematical aspect that is central to our analysis, and that leads to the determination of the constant parameter q in equations (68) and (69), is that the perturbed flux along the separatrix,ψ V+ (0, v), in the vicinity of the upper X-point at v = 0, can be written as the sum of even and odd functions of v,ψ With reference to equation (69), we point out that the even function ψ even (v) arises from the term proportional to sin πt/2, while the odd function ψ odd (v) arises from the term proportional to cos πt/2. The two functions can be considered as independent solutions, and so they can be discussed separately. Clearly, for the even mode, q = p/2, while for the odd mode, q = −p/2. Analogous conclusions can be obtained if the lower X-point at v = π is considered instead. We have obtained asymptotic relations for the perturbed flux in all four regions surrounding the X-point. Thus, the perturbed current density along the separatrix, in the vicinity of the X-points, can be determined fromJ(u, v) = |∇u| 2 Obviously, the perturbed current vanishes everywhere except on the separatrix. Therefore, we can setJ(u, v) = j X (v)δ(u). It follows from (59), (71) and (73), that, for both even and odd solutions, where use has been made of the fact that, near the upper Xpoint, |∇u| 2 (0, v) ∼ (2e 0 /ab)v, which can be easily obtained by direct expansion in the vicinity of X-point (see appendix H). Clearly, for the solution for the perturbed flux that is odd as a function of v near the upper X-point at v = 0, q = −p/2 and the current sheet vanishes altogether. A current sheet is found only for the even-ψ solution, for which p/2 + q = p. The asymptotic behavior of the stream function,φ(0, v), along the magnetic separatrix and in the vicinity of the X-point can also be easily worked out. We find: Note that also the stream function presents a similar type of square-root singularity as a function of v approaching the upper X-point at v = 0, as it should, as this guarantees that the inertial term in the perturbed equation of motion (24) can balance the force density term [ψ eq ,J]. Finally, the dispersion relation can be obtained from the plasma equation of motion, following a procedure that is similar to the one used in section 4, equations (32)- (36). Equation (24) can be written as γ∇ · (ϱ ∇φ) = ∇ · (J e z × ∇ψ eq ). We integrate with respect to u over a narrow interval across the separatrix between regions ∆ and V + , in the vicinity of the upper X-point, where v is positive and small: Note that, as we anticipated below equation (77), a consistent behavior is found, as both terms in equations (78) and (79) are proportional to v −1/2 . Balancing the two terms, we obtain: where ω A = τ −1 A and dimensions have been reintroduced for reasons of clarity.
The following considerations can be drawn at this stage. For the odd-parity solution, q = −p/2, the n = 0 mode is neutrally stable with γ = 0 and no current sheet develops at the magnetic separatrix. The stream function also vanishes in this limit, as φ ∝ γ. This solution can be considered merely as a redefinition of the equilibrium, with the current-carrying plasma shifted vertically by a distance ξ and the equilibrium current density modified by current sheets located at the elliptical flux surface For the even-parity solution, q = +p/2, a current sheet develops at the separatrix. Evidently, this is sufficient to neutrally stabilize the n = 0 mode, which in this case oscillates with a real frequency: where equation (57) has been used for the parameter p ∝ e 1/2 0 , and κ = b/a. For typical elongations of present-day tokamak experiments, for instance, κ = 1.5, which corresponds to e 0 = 0.4, one finds ω ≈ 0.6ω A . Assuming a hydrogen plasma with a = 1 m, B p ′ = 1 T m −1 , and number density n = 10 20 m −3 , we find ω ≃ 200 kHz. Since ω A is the Alfvén frequency based on the poloidal magnetic field and e 0 < 1, the mode frequency falls below the Alfvén continuum spectrum and therefore is unaffected by continuum damping. The even-parity mode can be destabilized by the resonant interaction with fast particle orbits [38,49] and becomes a possible candidate for the interpretation of finite-amplitude n = 0 fluctuations recently observed in JET experiments [42,43].

Searching for ideal-MHD marginal stability
As we have seen, current sheets always form at the plasma boundary, where the density drops to zero. When the plasma boundary coincides with the convenient elliptical surface at µ = µ b (or, equivalently, u = u b = −µ b + (e 0 /2) sinh 2µ b ), vertical displacements were found in section 4 to be ideal-MHD unstable, unless a nearby metallic wall providing passive stabilization was present. By contrast, when the boundary coincides with the magnetic separatrix at u = 0, vertical displacements were found in section 5 to be ideal-MHD stable even in the absence of a nearby wall. Evidently, the nature of the current sheet changes significantly as the plasma boundary approaches the magnetic separatrix. We conjecture that a flux surface u = u marg < 0 exist, such that, when the plasma boundary is at u = u marg , marginal ideal-MHD stability against vertical perturbations is obtained. Finding the value of u marg requires numerical work, which is beyond the scope of the present article. However, in the following, we provide further physical insight into the properties of the dispersion relation as the plasma boundary approaches the magnetic separatrix.
Let us consider the situation where the plasma boundary is at a magnetic flux surface, u = u c < 0, located between the convenient elliptical surface and the magnetic separatrix: 0 > u c > u b . A generic flux surface u = u c is shown in figure 4. Region ∆ is now the region between the elliptical surface (the contour of region Ω) and the plasma boundary at u = u c , while region V ′ denotes the vacuum region between the plasma boundary and the magnetic separatrix. In regions Ω and ∆ we can keep on using the solutions for the perturbed flux,ψ, and for the perturbed stream function,φ, that were found in section 5. Thus, the solution forψ ∆ in region ∆ expressed in flux coordinates (u, v) is given by equation (50), with Fourier coefficients α m and β m fully determined by equations (53) and (54). On the plasma boundary, the perturbed flux has a modulation in v given by: where λ m = α m cosh(mu c ) + β m sinh(mu c ). Figure 5 shows the Fourier coefficients λ m as functions of m for a fixed value of the ellipticity parameter, e 0 = 0.05, corresponding to u b = −1.34, and two negative values of u: u c1 = −1.0 and u c2 = −0.1. As we can see, the Fourier spectrum extends up to m = m max ∼ |u c | −1 when |u c | < 1. This statement can also be confirmed by analytic asymptotic expansion of the coefficients λ m , which shows that λ m → α m ∼ (p/m 3/2 ) exp (mu c ), and so the Fourier spectrum extends to infinity when u c → 0. It is important to note that λ 1 is negative, while all other λ m with m > 1 are positive.
The following conclusion can be drawn at this stage. For values of |u c | ⩾ 1, the perturbed flux and stream function are dominated by the single m = 1 harmonic in the flux coordinate angle v, which, within the region bounded by the flux surface u = u c , can be approximated by the elliptical angle θ. In this case, it can be shown analytically that vertical displacements are ideal-MHD unstable in the absence of a nearby wall. The situation changes when |u c | < 1. In this case, more and more Fourier harmonics in the Fourier series forψ ∆ must be included in the analysis and the nature of the dispersion relation changes significantly. Since λ 1 is negative, while all other Fourier coefficients λ m with m > 1 are positive, we expect that additional Fourier harmonics will give a stabilizing contribution to the relevant dispersion relation. In particular, we expect vertical displacements to become ideal-MHD stable, even in the absence of a nearby wall, for negative values of u c larger than u marg (i.e. approaching the magnetic separatrix at u = 0), with the value of |u marg | < 1 depending in turn on the ellipticity parameter e 0 .
The perturbed flux in the vacuum regions V ± and V ′ satisfies ∇ 2ψ = 0. The solution in elliptical coordinates that decays to zero for µ → ∞ (no wall) is given by equation (47), where the summation involves only odd integers. Now, we are confronted with the task of determining the coefficients c m . We proceed as follows. The vacuum solutions in region V ′ can also be expressed in (u, v) coordinates as: This expression is similar to the one in equation (50); however, the coefficients γ m and η m differ from α m and β m . In order to determine these coefficients, as well as the c m coefficients, we note thatψ V ′ (µ, θ) andψ V ′ (u, v) have equal values on points in region V ′ . There is nothing wrong with extending these functions to region ∆ and all the way to the elliptical surface µ = µ b , which is also a constant u = u b surface. Therefore, we can use the relationsψ A third relation we can use is continuity of the perturbed flux on the plasma boundary u = u c : After relatively straightforward algebra, we arrive at the set of three coupled equations: In these equations, and the coefficients λ m (u c ), which depend on α m and β m , are defined below equation (82). Combining these equations, we finally arrive at: Equation (90) is an infinite set of linear equations that, in principle, can be used to determine the coefficients c n , and hence γ m and η m . In practice, the procedure works if we can truncate the series up to n max = m max , where the value of m max is discussed below equation (82). Then, equation (90) reduces to a set of (m max + 1)/2 equations for the (n max + 1)/2 unknown quantities c 1 , c 3 , . . ., c nmax . The dispersion relation for γ 2 can be obtained from the plasma equation of motion, following a procedure that is similar to the one used in section 4. After straightforward algebra, we obtain: The stream function in region ∆ corresponds to the rigid shift, which, in elliptical coordinates, is given by equation (27). In (u, v) coordinates, g m e −mu + k m e mu sin(mv), where the coefficients g m and k m can be determined following the same procedure that led to the determination of the coefficients α m and β m ; see equations (51) and (52). After straightforward algebra, we obtain: Putting all these pieces together, we can obtain a dispersion relation for γ 2 as a function of u c . However, the determination of the Fourier coefficients appearing in the dispersion relation requires numerical work and will be the subject of a future publication. What we can conclude at this stage is that, for as long as |u c | ∼ 1 and a single harmonic in the angle θ (or equivalently in the angle v) can be used for a good approximation of the fieldsφ andψ, analytic work is relatively simple, and this will result in unstable vertical displacements when the confining metallic wall is moved beyond the X-points. The change in sign of γ 2 , leading to ideal-MHD stable vertical displacements, is mathematically associated with the fact that more and more Fourier harmonics in θ (and in v) become important as |u c | becomes smaller and smaller. Analytic work is also possible when u = u c = 0; this case was treated in section 5 and resulted in stable vertical displacements even in the absence of a nearby wall.

Conclusion
This article contains three main results. First, we have shown that a rigid-shift displacement is indeed the analytic solution for the normal mode analysis of axisymmetric n = 0 modes. This result has been obtained on the basis of the linearized, reduced ideal-MHD model, starting from a relatively simple 'straight tokamak' equilibrium. Our result for the growth rate of ideal-MHD unstable vertical displacements, valid for arbitrary values of the ellipticity parameter e 0 , is given in equation (38). In the limit of small e 0 , the growth rate scales as γ ∼ e 1/2 where τ A is the relevant Alfvén time. With feedback stabilization, the n = 0 mode becomes oscillatory in nature, with a frequency ω ∼ e A . This is what we have called in section 4 the 'limiter tokamak scenario' for n = 0 vertical modes. We have also shown that the small ellipticity limit of the obtained dispersion relation is in perfect agreement with that obtained by a simplified heuristic model that treats the plasma current, as well as the external currents, as three parallel current filaments. These results, however, are correct only as long as the plasma does not extend to the magnetic separatrix, where magnetic Xpoints are located, i.e. as long as the X-point resonance can be neglected.
Secondly, for the more relevant case where the plasma does extend to the magnetic separatrix (the 'divertor tokamak scenario' discussed in section 5), the resonant nature of the magnetic X-points with respect to axisymmetric perturbations is analyzed in detail. It is found that vertical displacements are stable, at least on ideal-MHD time scales, without any need for passive stabilization elements (which could easily be included in our model by changing the boundary conditions for the perturbed flux in the vacuum regions V ± ). The stabilization mechanism is a direct consequence of the ideal MHD fluxfreezing constraint on the X-points. Indeed, as also discussed in [39], the theory discussed in this article presents analogies with the physics of current sheet formation during the evolution of internal kink modes [35,36], and with the magnetic island coalescence problem [37]. In these works, the ideal-MHD constraint causes magnetic flux to pile up near the X-points, leading to perturbed localized currents and a stabilizing effect in the ideal-MHD limit. For the island coalescence problem, it was found that a chain of magnetic islands becomes ideal-MHD unstable when the island width exceeds a critical threshold. In any case, flux pile-up prevents any further nonlinear evolution for both unstable internal kinks and island coalescence, unless the ideal-MHD constraint is relaxed, e.g. by resistivity. In our problem, vertical displacements are found to be linearly stable in the ideal-MHD limit, when the mode resonance at the equilibrium X-points of the divertor separatrix is properly taken into account. As a consequence, axisymmetric perturbed currents are excited in the vicinity of the magnetic X-points, and these currents carried by the plasma are capable of providing passive feedback stabilization of vertical displacements on the ideal-MHD timescale, even in the absence of a nearby metallic wall. Results for the 'divertor tokamak scenario' were presented in [39]; however, mathematical details are provided in this article for the first time.
Finally, we have discussed the scenario where the plasma boundary is located near the magnetic separatrix, but does not extend to include the magnetic X-points. The analysis shows that a rigid-shift displacement is still the analytic solution for the normal mode analysis of axisymmetric n = 0 modes. Noting that current sheets always form at the plasma boundary, where the plasma density drops to zero, it is shown that the nature of the boundary current sheet changes significantly as the plasma boundary approaches the magnetic separatrix. It is argued that a magnetic surface u = u marg < 0 (with u = 0 corresponding to the magnetic separatrix) must exist, such that, when the plasma boundary coincides with this magnetic surface, vertical displacements are marginally stable according to ideal-MHD and in the absence of a nearby wall. Finding the actual value of u marg requires numerical work and will be the subject of a future publication.
The resonant behavior of n = 0 modes at the X-points of a tokamak magnetic separatrix is an ideal-MHD property. Therefore, it is reasonable that it must be studied first according to the ideal-MHD model. Future work will consider extended-MHD effects, but this article, together with its shorter version in [39], represents the starting point for future developments. The most important, and perhaps surprising, result of our work is that, when the plasma density extends to the magnetic separatrix, vertical displacements are found to be stable, on ideal-MHD time scales, without any need for a passive stabilization wall. The stabilization mechanism is a direct consequence of the ideal-MHD flux-freezing constraint on the X-points, which leads to current sheets localized along the magnetic separatrix, exerting a force capable of pushing back the plasma in its vertical motion. This also suggests that plasma electrical resistivity in a narrow boundary layer along the magnetic separatrix, in addition to wall resistivity, may have a profound impact on the stability of n = 0 vertical displacements.

Data availability statement
All data that support the findings of this study are included within the article (and any supplementary files).

Appendix B. Geometrical factor D for ideal elliptical wall stabilization
For confocal ellipses, the following relations hold: It follows that: Making use of (B.3), (B.4), (6) and (40), relation (43) becomes: In the neutral stability case, where the wall intercepts the Xpoints, we have µ w = 2µ b . Substituting into the expressions above and with the help of (B.1) and (B.3), we can easily show that D = 1 in this case.