J. Plasma
J.
Plasma Fusion
Fusion Res.
Res. SERIES,SERIES,Vol.
Vol.6 6(2004)(2004)207–209
000–000

Theory of the Newcomb Equation and Applications to MHD
Stability Analysis of a Tokamak
TOKUDA Shinji and AIBA Nobuyuki1
Naka Fusion Research Establishment, Japan Atomic Energy Research Institute, Ibaraki, 311-0193, Japan
1
The Graduate University for Advanced Studies, Toki, 509-5292, Japan(Received: 9 December 2003 / Accepted: 25 April 2004)

Abstract
Recent progress in the theory of the Newcomb equation is reported. Emphasis is put on the analysis of externalmodes including peeling modes (high n kink modes), where n is the toroidal mode number. A theory for low n
external modes is developed so that it is also useful for the analysis of resistive wall modes.
Keywords:
Newcomb equation, eigenvalue problem, ideal MHD stability, external mode, edge stability, tokamak

1. Introduction

Here, R0 is the position of the magnetic axis (for simplicity,
the mirror symmetry of the equilibrium is assumed:ψ (R, –Z ) = ψ (R, Z )); q(r) is the safety factor, as usual. We
define the poloidal angle θ so that the magnetic field lines
are straight; the Jacobian of the coordinate system (r, θ, φ) is→⎯√ g = rR2 /R0. Let ζ be an infinitesimal plasma displacementwith the toroidal mode number n that is incompressible→(∇ · ζ = 0), and let

It is well known that the Newcomb equation, the inertia
free linear ideal magnetohydrodynamic (MHD) equation [1],plays fundamental roles in the MHD stability theory. A code
MARG2D has been developed which solves numerically the2-dimensional Newcomb equation and the associatedeigenvalue problem [2]. In this paper, recent research on theNewcomb equation is reported. The main focus is to develop
tools for the analysis of low n and high n external modes,
where n is the toroidal mode number. The high n kink modes,
called peeling modes [3], recently get attention in the study
on MHD stability of tokamak edge plasmas [4]. For the lown external modes, we develop a theory that expresses thechange of potential energy due to the plasma displacement
by a quadratic form with respect to the values of thedisplacement at the plasma surface. This formulation is usefulfor the analysis of resistive wall modes. For the analysis ofpeeling modes we extend the MARG2D formulation into thevacuum region by expressing the perturbation of magnetic
fields in vacuum by a suitable vector potential.

→X (r ,θ ) := ζ ⋅∇r ,→1→V (r ,θ ) := r ζ ⋅∇θ − ζ ⋅∇φ ,q(Wp = π ∫ L drd θ ,ψqdψ .
F(5)

L = a | Dθ ( X ) |2 + c | ∂ r (rX ) + ∂ θV | 2+ b | inV +1∂ (rX ) + hX + r βr θ Dθ ( X ) | 2
q r+ e | X |2 .

(6)

Here, the operator Dθ (X) is defined by(1)Dθ ( X ) :=where the cylindrical coordinate system (R, Z, φ) is employed;ψ (R, Z) and F(ψ ) are, respectively, the poloidal flux function
and toroidal field function. We define the radial coordinateby0(4)and the Lagrangian density function L reads

In an axisymmetric toroidal system such as a tokamak,equilibrium magnetic fields are expressed asr 2 := 2 R0 ∫)then the change of the plasma potential energy Wp due to the→displacement ζ reads [5]2. Newcomb equation

Beq = ∇φ × ∇ψ + F ∇φ ,(3)1∂θ X − inX ,q(7)and the other coefficients are given in ref. [2].
By minimizing Wp with respect to V(r, θ), we obtain thereduced energy integral(2)

Corresponding author’s e-mail: tokudas@fusion.naka.jaeri.go.jp1207©2004 by The Japan Society of Plasma
Science and Nuclear Fusion Research

Tokuda S. et al., Theory of the Newcomb Equation and Applications to MHD Stability Analysis of a TokamakaW p [ X , X ] = 2π 2 ∫ L[ X , X ]dr .

(8)

01t
MH := ( M + M ).
2Here, the vector function X(r) is defined as
X(r) := {X– L f (r), . . . , XLf (r)}t ,(17)Now let us make a vector function Y m(r) ∈ S(9)mmtY m(r) = (Y–Lf (r), . . . , YL f (r)) ,by using the poloidal Fourier harmonics Xl (r)for m = 0, ±1, . . . , ± L f , where each poloidal harmonics
Ylm(r) satisfies the condition

LfX (r ,θ ) = ∑ X l (r ) exp(ilθ ),(18)(10)

l =− L fYlm(a) = 0

where Lf is the truncated poloidal mode number. And thereduced Lagrangian density is(l ≠ m),Ymm(a) = 1,

l = 0, ±1, . . . , ± L f .

(19)

〈〉
dXdX
(11)
+〈|M |X〉 +〈X|M |,drdr 〉The set {Y (r)} forms a basis [6]. An external mode can be
expressed by using an arbitrary set {xm} of real numbers as

where L, M, K are matrices; L and K are hermitian, the detailsof which are given in ref. [2], andThe change of the potential energy due to ξ is given by thequadratic form of the vector x,〈 X | K | X 〉 := ∑ X j K jk X k .

W[ξ , ξ ] = 〈x| A|x〉,

L [ X, X ] =mdX
dX| L|+〈X|K|X〉
drdr→

ξ (r ) = ∑ xm Y m(r ).

t→→ →j ,kN X := − d ⎛⎜ L dX ⎞⎟ − d ( M tX )dr ⎝ dr ⎠ drdX
+ K X = 0.
dr

N X = –λR X,(13)where R is a multiplicative and diagonal operator whose
components are Rm,m ∝ (m/q – n)2. The natural boundary
condition for Xm is imposed at the rational surface of rm (m =nq(rm)); the continuous conditions are imposed for otherharmonics X l (l ≠ m). A code MARG2D which solves
eqs. (12) and (13) has been developed by using a finite
element method. The code has been applied to identify stablestates for ideal MHD internal perturbations [2].

→→N X m(r) = –N Zm(r),→(15)→ →→→(25)The vacuum contribution to the change of potential
energy is represented by the matrix MV in eq. (14), which iscomputed by using Green’s function of the Laplace operator[9]. This method deals flexibly with the shape of a conductingwall. However, the method is limited to low n modes since itis difficult to evaluate numerically special functions thatappear in constructing M V when n becomes large. It isconvenient to express magnetic fields by a vector potential

→dη a
1 →
ξ a | L|2dr→

1 dξ a→+| L |η a ,2 dr

X m(a) = 0.

4. Extension of the MARG2D formulation
into the vacuum

and the matrix M V stands for the contribution from the→→vacuum region. Let S = {ξ |N ξ = 0} be a set of functions→→
that satisfy the Newcomb equation. If ξ (r), η (r) ∈ S, thenwe haveWp [ξ ,η ] = 〈 ξ a | M H | η a 〉 +(24)Since eq. (25) can be solved by the MARG2D code, we canconstruct the basis {Y m(r)} and the stability matrix A.
The present formalism has been satisfactory applied tothe stability analysis of low n, typically n = 1 or 2, external
modes, which is reported at this conference [8].

where→(23)where Z (r) given analytically satisfies the inhomogenous
boundary condition, eq. (19). Then, we have an inhomogeneous equation for X m(r) with the homogenous
boundaryn condition(14)ξ a = ξ (a),Mp(l,m) = Wp[Y l, Y m ].

m

The bilinear form associated with eq. (8) is given by→ →(22)Y m(r) = X m(r) + Zm(r),3. Application to the theory of externalmodes
W[ξ , η ] = Wp[ξ , η ] + 〈ξ a | MV |η a〉,

A = Mp + M V ,We call A the stability matrix for external modes. If theminimum eigenvalue of A is negative, then the plasma isunstable against ideal external kink modes. The matrix A also
plays an important role in the stability of resistive wall modes[6].
The basis {Y m(r)} can be constructed by using theresponse formalism [7]. Let us write Y m(r) as(12)The eigenvalue problem associated with eq. (12) is given by→ →(21)where the matrix A, which is real and symmetric, is given by

From L[X, X], we have the 2-D Newcomb equation+M(20)m(16)2208

Tokuda S. et al., Theory of the Newcomb Equation and Applications to MHD Stability Analysis of a Tokamak

for middle n (n = 2, 3, . . . , 10) or high n (n > 10) externalmodes, which was shown in ref. [10]. We adopt the vector
potential method in the MARG2D code for the analysis of
peeling modes. It is shown that the Lagrangian density
function for the change of energy in the vacuum has the same
form as that for the change of plasma potential energy, andthat the MARG2D formulation is easily extended to suchexternal modes.
Let ψ V (R, Z ) be a function defined in the vacuum. Thecoutour ψ V (R, Z ) = ψs coincides with the plasma surfaceψ (R, Z ) = ψs. We also assume that the outermost contour ofψ V (R, Z ) = const. coincides with the cross section of theconducting wall. Next, we introduce a vector field CV in thevacuume byCV = ∇φ × ∇ψ V + TV (R, Z )∇φ .

anda :=L (Y, Y ) =.

(36)

Y (ψ V , θ ) = Σ Y l (ψ V ) exp (il θ ) ,(37)Y(ψV) := {Y–L f (ψV), . . . , YL f (ψV)}t.

(38)

(39)

dY L dY + Y K Ydψ Vdψ V+dY M t Ydψ V+ Y M dY . (40)dψ V

We have reported the recent progress in the theory ofthe Newcomb equation and applications to MHD stability
analysis. A theory of external modes has been developed byusing the Newcomb equation, which is also useful for theanalysis of resistive wall modes. The MARG2D formulation
has been extended in the vacuum region. Such extension
enable us to analyze the stability of the edge region of atokamak, and code developing is now on going.

(27)

and the vector potential A is expressed by→A = ξ V × CV,(29)

Acknowledgements

where ξ V is the unknown vector to be determined. By
introducing the functions→Y(ψV , θ ) := ξ V · ∇ψV ,The authors wish to thank to Dr. H. Ninomiya and Dr.
Y. Kishimoto at JAERI for supporting our research.

(30)→
1 →V (ψ V , θ ) := ξ V ⋅ ∇θ – q ξ V ⋅ ∇φ ,References[1] W.A. Newcomb, Ann. Phys. 10, 232 (1960).
[2] S. Tokuda and T. Watanabe, Phys. Plasmas 6, 3012
(1999).
[3] D. Lortz, Nucl. Fusion 15, 49 (1975).
[4] W. Connor, R.J. Hastie and H.R. Wilson, Phys. Plasmas
5, 2687 (1998).
[5] I.B. Bernstein, E.A. Frieman, M.D. Kruskal and R.M.
Kulsrud, Proc. Roy. Soc. London A 244, 17 (1958).
[6] A.H. Boozer, Phys. Plasmas 5, 3350 (1998).
[7] S. Tokuda and T. Watanabe, J. Plasma Fusion Res. 73,1141 (1997).
[8] N. Aiba, S. Tokuda, T. Ishizawa and M. Okamoto, The
effect of the aspect-ratio on the external kink-ballooning
instability in high beta tokamaks, this conference.
[9] R. Gruber et al., Comput. Phys. Commun. 21, 323(1981).
[10] R. Gruber et al., Comput. Phys. Commun. 24, 363(1981).

(31)

sthe change of energy in the vacuum is given by

WV = π qs ∫ L dψ V dθ ,(32)L = a D θ (Y) + T ∂∂Vθ + ∂∂ψY22V
2∂Y+ b inq s V +β D (Y) .
∂ ψ V ψθ θ(33)1Dθ (Y ) := q ∂θ – in Y ,s(34)

Here,

TV

WV = 2π 2qs ∫ TV (ψ V L (Y, Y) dψ V ,(28)2, b :=

5. Summary

B = ∇ × A,∇ψ V

2The reduced form has the same form in eq. (11):at all points in the vacuum, where qs is the safety factor atthe edge.
The perturbation of magnetic fields is given by∇ψ V · ∇θ

22l(26)

CV ⋅∇φ= q s (= const.),CV ⋅ ∇θβ ψθ := q sR ∇ψ V∇ψ VWe can eliminate the function V by the same procedure inSec. 2 by using the poloidal Fourier harmonics

Here we assume TV is independent of φ. It is easy to see that
CV is a solenoidal vector [10], div CV = 0. The poloidal angle
θ and the function TV can be defined such that→TV(35)3209