J. Plasma Fusion Res. SERIES, Vol. 9 (2010)

Eigenvalue Spectrum of MHD Modes in Cylindrical Tokamak Plasmas
with Small Resistivity
Taro MATSUMOTO1 and Shinji TOKUDA2
1

Fusion Research Development Directorate, Japan Atomic Energy Agency
2

Research Organization for Information Science and Technology(Received: 20 November 2009 / Accepted: 5 March 2010)

For the understandings of the magnetohydrodynamic (MHD) characteristics in plasma, the spectrum of theresistive MHD modes are investigated in detail by solving the eigenvalue problem of the reduced MHD equations
in cylindrical tokamak plasmas. The eigenvalues and eigenfunctions of the resistive MHD modes are clarified forsmall toroidal and poloidal mode numbers, and the discrete eigenvalues with imaginary part are categorized into
several types by analyzing the behavior of their eigenfuctions. In particular, the dependence of the eigenvalue
distribution in the real-imaginary space on the resistivity is studied for small resistivity regime by parallel
computation. It is found that the characteristics of eigenvalue distribution and eigenfunction are changed even forsmall reduction of resistivity in small resistivity region. It is also found that the distribution structure of theeigenvalues is deformed with bifurcation of spectral curve.
Keywords: magnetohydrodynamics, stability, resistive MHD, eigenfunction, eigenvalue, tokamak,
reduced MHD equations, eigenvalue problem, shear Alfven wave, parallel computation

clarified the eigenvalue distribution, its dependency on
resistivity, and characteristics of eigenfunctions.
This paper is organized as follows. In Section 2, we
describe our simulation method to solve the MHD problem,
including the reduced MHD equations, the matrix
eigenvalue problem formulation, the computational
approaches and the equilibrium parameters. In Section 3,
the eigenvalue distribution and the eigenfunction profile
for the m/n = 1/1 and 2/1 modes are described in detail.
Here m and n are the poloidal and toroidal mode number,respectively. The dependency of the eigenvalue distribution
on the plasma resistivity is also shown. Finally, a briefsummary is given in Sec.4.

1. Introduction

For the understandings of the tokamak plasma
dynamics, it is important to investigate and clarify thecharacteristics of the magnetohydrodynamic (MHD)
spectrum, which is also particularly important from theviewpoint of the MHD spectroscopy.
The characteristics of the MHD spectrum had been
investigated by solving the eigenvalue problem based on
the ideal and resistive MHD equations in cylindricaltokamak plasmas [1-2].
In the framework of the ideal MHD model, theeigenvalue problem of a MHD mode with the timedependence proportional to eOt has one or some unstablediscrete eigenvalues on the positive real axis andcontinuum spectrum on the imaginary axis. However, inthe presence of finite plasma resistivity, it was found thatthe continuum spectrum is modified into purely damping
spectrum on the negative real axis and discrete brancheswith finite imaginary and negative real components in thecomplex plane.
Due to the limitation of computational resources inthose days, the eigenvalue distribution had been analyzed
for relatively large resistivity, and the eigenfunctions ofsome distinguished eigenvalues was clarified.
Thus, in this paper, we have analyzed the eigenvalue
distribution of cylindrical tokamak plasmas withasymptotically smaller resistivity than the previous works
by using current computational resources. We have also

2. Simulation Method
2.1 Basic equations

In this research, we employ a reduced MHD model
where plasma can be described with the equations which
do not include the fast magnetosonic wave but the shear
Alfven wave. In this model, the dynamics of cylindricalplasma is described by the following reduced MHD
equations for the vorticity U and the magnetic flux \,BwU[ J ,\ ]linear  in 0 J  Q'I ,
wt
R0
B
w\
[\ , I ]linear  in 0 I  KJ ,
wt
R0

(1)
(2)

author’s e-mail: taro_matsumoto@yahoo.co.jp, matsumoto.taro@jaea.go.jp

568©2010 by The Japan Society of Plasma
Science and Nuclear Fusion Research

T. Matsumoto and S. Tokuda, Eigenvalue Spectrum of MHD Modes in Cylindrical Tokamak Plasmas with Small Resistivity

U 'I ,
J '\ ,(3)(4)

Safety factor profile
for m=1 mode simulation

(a)
2.0
㼝

Safety Factor : q

which are derived from full MHD equations on theassumption of a large aspect ratio, where t is the time, I is
the electrostatic potential, J is the current density, K is theplasma resistivity and Q is the viscosity. Here, acylindrical coordinates (r,T,z) is used. The brackets in
Eqs.(1) and (2) denote the linear part of the Poissonbracket,
(5)

3.5

2

(6)

Fig.1

A\\

A\I

(7)

where AXY denotes the regular square matrix of the linear
interaction, O is a eigenvalue, and Xmn is a eigenfuction ofphysical quantitysuch as the vorticity U, the magnetic
flux \, the electrostatic potential I, and the current
density J.
By substituting Eqs.(3) and (4) into Eq.(7), an
simplified set of equations,§U · § A

A

'1

UI
O ¨¨ mn ¸¸ ¨¨ UU
1
A
A

\
\I '
© mn ¹ © \U

c
§ AUU
¨
¨ Ac
© \U

AU\  AUJ ' ·§U mn ·¸¨¸A\\  A\J ' ¸¹¨©\ mn ¸¹

AUc \ ·§U mn ·¸¨¸,c ¸¹¨©\ mn ¸¹
A\\

0.8

1.0

0.8

1.0

Safety factor profile
for m=2 mode simulation
q

2.5
2.0
1.5
1.0

0.2

0.4

0.6

Radius : r

Safety factor profile for (a) m = 1 mode and(b) m = 2 mode.

For the numerical calculation of Eq.(8), a eigenvalue
problem code has been developed. Here, at first, EISPACK
is chosen as the science library for a non-parallel
computation in order to solve a two-dimensional matrix.
The validation of this eigenvalue problem code wasconfirmed by comparing with the initial value code which
was independently developed, from the viewpoint of theeigenvalue (the growth rate), the eigenfunction, and grid
number dependency for the m/n = 1/1 and 2/1 unstablemodes.
The number of the radial grid points, imax, is changed
from 512 to 16385. For the large radial grid points morethan 1024, the eigenvalue problem code is parallelized, andScaLAPACK scientific library is used to solve twodimensional block cyclic distributed matrix at eachprocessor element.
As the computational resources, we have used aworkstation and SGI Altix3900 supercomputer with up to
64 processor elements. The computational time for thesame eigenvalue problem is significantly reduced due tothe effect of changing processor element and scientificlibrary, and parallelization. As a result, the computationaltime has successfully reduced by about 4000 times.
From the viewpoint of numerical accuracy, it is foundthat in order to obtain precise eigenvalue distribution, goodresolution should be kept over the whole radius. Thus, theuniform radial grid is used for the physical computation.

The formulation of matrix eigenvalue problem isobtained from Eqs.(1) and (2) as§U mn · §U mn ·
AUJ ·¨\ mn ¸ ¨\ mn ¸¸¨ ¸ O ¨ ¸
A\J ¸
¹¨I mn ¸ ¨I mn ¸¨ ¸ ¨ ¸©J mn ¹ ©J mn ¹

0.6

Radius : r

3.0

0.0
0.0

2.2 Formulation of matrix eigenvalue problem

AUI

0.4

0.5

B0 is the toroidal magnetic field and 2SR0 is thelongitudinal length in z direction.

AU\

0.2

4.0

Safety Factor : qw .
1w § w · 1 w¨r ¸ 
r wr © wr ¹ r 2 wT 2 wz 2

§AUU
¨
¨A
© \U

0.5

0.0
0.0

The delta, ', in Eqs.(3) and (4) denotes the Laplacian in
the cylindrical coordinates, and
'

1.0

(b)

1 §wa wb wa wb ·[a, b]¨¸.
r ©wr wT wT wr ¹

2

1.5

(8)

is obtained, where '-1 is the inversion of the ' of the delta.
Thus, the eigenvalue problem of a restive MHD mode ina cylindrical tokamak plasma is reduced to finding the
eigenvalues and eigenfunctions of Eq.(8).

2.4 Equilibrium parameters

2.3 Computational approach

In this research, the safety factor profiles are assumed

569

T. Matsumoto and S. Tokuda, Eigenvalue Spectrum of MHD Modes in Cylindrical Tokamak Plasmas with Small Resistivity(a)1.0 10-1

5.0 10-2㻱㼕㼓㼑㼚㼢㼍㼘㼡㼑㻌㻰㼕㼟㼠㼞㼕㼎㼡㼠㼕㼛㼚㻌㼛㼒㻌㼙㻩㻝㻌㻹㼛㼐㼑㼑㼠㼍㻜㻩㻝㼑㻙㻡㻘㻌㼢㼕㼟㻜㻩㻜㻚㻜㼑㻙㻜㻘㻌㼕㼙㼍㼤㻩㻡㻝㻞
Eigenvalue

0.08

%
%
%

Imaginary

0.06
0.0 100

0.04'-5.0 10-2-1.0 10 1-8.0 10 0-6.0 10 0Real-4.0 10 0-2.0 10 0

0.0 10 0

0.00)

1.0 10 -1

㻱㼕㼓㼑㼚㼢㼍㼘㼡㼑㻌㻰㼕㼟㼠㼞㼕㼎㼡㼠㼕㼛㼚㻌㼛㼒㻌㼙㻩㻝㻌㻹㼛㼐㼑㼑㼠㼍㻜㻩㻝㼑㻙㻡㻘㻌㼢㼕㼟㻜㻩㻜㻚㻜㼑㻙㻜㻘㻌㼕㼙㼍㼤㻩㻡㻝㻞Imaginary

5.0 10 -2-5.0 10 -2

Fig.2

0.0 10 0Real

5.0 10 -2

-0.03

$

-0.02

'
'

-0.01

'

$

0.00

0.01

which situate off the real axis and have finite imaginary
component and negative real component. Theseeigenvalues have a structure which is imaginarily
symmetrical, as shown in Fig.2(b), where the scale of thereal axis is almost the same to that of the imaginary axis.
This structure has five branches as follows: one branchextending toward the origin, one branch extending toward
the real axis, two branches towards the imaginary axis,
and one branch connecting the above four branches. It is
confirmed that these discrete eigenvalues do not change
as the radial grid number increases while fixing physicalparameters.
There is also one unstable mode on the real axis in
its positive region, that is, the m = 1 kink mode, which
has the eigenvalue of approximately 4 x 10-2.

0.0 10 0-5.0 10 -2

-0.04

)

Fig.3 Detail figure of eigenvalue distribution of m
= 1 mode. Here, K= 1.0 x 10-5, Q= 0, and imax
= 512.

Eigenvalue-1.0 10 -1-1.0 10 -1

%

$

-0.05

(b)

'

%
0.02

-1.0 10 -1-1.2 10 1'

1.0 10 -1

Eigenvalue distribution of m =1 mode in (a)macroscopic scale and (b) microscopic scale.
Here, K= 1.0 x 10-5, Q= 0, and imax = 512.

as shown in Fig.1(a) and (b) for m = 1 (m/n = 1/1) modeand m=2 (m/n = 2/1) mode, respectively. The resonant
magnetic surface locates at almost half plasma radius for m= 1 mode, and at approximately 0.7 of the plasma radius
for m = 2 mode. For simplicity, the toroidal mode number
is chosen to be 1.
The plasma resistivity K is changed from 10-3 to 10-7.
The viscosity v is mostly assumed to be 0. However, it was
confirmed that the eigenvalue distribution and theeigenfunction do not change significantly even for theviscosity equal to the plasma resistivity.

3.2 Eigenfunction of m = 1 mode

For the better understanding, the detailed figure of
the eigenvalue distribution of the m=1 mode is describedin Fig.3, where some characteristic eigenvalues are
identified with denotations.
The eigenfunctions of D8 and D9 on the branch tothe origin are shown in Fig.4(a), with A0 of the internalkink mode. The eigenfunction is normalized so that theabsolute maximum value is unity. It is found that thesethree modes are localized in the core plasma region,
namely within the resonant magnetic surface of the m = 1
mode shown in Fig.1(a). However, compared with A0,
D8 and D9 have fluctuation component, in particulararound the resonant surface.
The eigenfunctions of D1, D2, and D3 on the 1st
branch toward the imaginary axis are shown in Fig.4(b).
It is found that these modes also localize in the coreplasma region, much more inward than D8, D9, and A0.
It is also found that, as departing from the imaginary axis,
that is, as decreasing the real component of theeigenvalue, the eigenfunction comes to have

3. Characteristics of Eigenvalue Distribution andEigenfunctions
3.1 Outline of eigenvalue distribution

At first, the macroscopic view of the eigenvalue
distribution of the m = 1 mode is shown in Fig.2(a), whereK= 1.0 x 10-5, Q= 0, and imax = 512. It should be notedthat the scale of the real axis is 100 times larger than thatof the imaginary axis to emphasize overview of theeigenvalue distribution. There are a lot of eigenvalues on
the real axis in its negative region. These discrete
eigenvalues indicate purely damping modes [3].
On the other hand, there are discrete eigenvalues

570

T. Matsumoto and S. Tokuda, Eigenvalue Spectrum of MHD Modes in Cylindrical Tokamak Plasmas with Small Resistivity

(b)

(a)

D1
D2
D3

D8
D9
A0

0.20
B1-B5

0.16
0.12

(d)

(c)

D6
F1
F2
B11

B1
B2
B3

B25

0.08

F2

0.00
-0.20

(e)

F1

B33

0.04

A977

-0.15

D1-D3

D7
D11

D12
A0

A976

-0.10

-0.05

0.00

Fig.5 Detail figure of eigenvalue distribution of m
= 2 mode. Here, K= 1.0 x 10-5, Q= 0, and imax
= 2048.

B14
A230
A231

MHD mode, m = 2 is chosen. The detailed figure of theeigenvalue distribution of the m = 2 mode is described in
Fig.5. There are also discrete eigenvalues with imaginary
component in the negative real region, which are
identified with denotations, as well as Fig.3. The
eigenvalue distribution has an imaginarily symmetricalstructure which comprises five branches and is similar tothat the m=1 mode although the 1st and 2nd branches
toward the imaginary axis overlap near the imaginary
axis.
There is one unstable mode, A0, on the real axis in
its positive region, that is, the m = 2 tearing mode, which
has the resonant magnetic surface at about r ~ 0.7 asshown in Fig.1(b) and has the eigenvalue ofapproximately 5 x 10-3.
The eigenfunctions of D11 and D12 on the branch tothe origin are shown in Fig.6(a), with the A0 of thetearing mode. It is found that these modes havefluctuation component around the resonant surface aswell as Fig.4(a).
The eigenfunctions of D1, D2, and D3 on the 1st
branch toward the imaginary axis are shown in Fig.6(b).
It is found that these modes localize in the edge plasma
region and are similar to B1 ~ B3 of the m=1 modeshown in Fig.4(c). On the other hand, the eigenfunctionsof B1 ~ B5 on the 2nd branch toward the imaginary axis
are shown in Fig.6(c). These modes are found to localizein the core plasma region and are similar to D1 ~D3 ofthe m=1 mode shown in Fig.4(b).
The eigenfunctions of D7, F1, F2, and B25 on theconnecting branch, and those of B33, A976 and A977 areshown in Fig.6(d) and (e), respectively. The
characteristics of D7, F1, F2, and B25, which are
generally symmetrical about the resonant surface, aresimilar to D6, F1, F2, and B11 of the m = 1 mode asshown in Fig.4(d). The characteristics of B33, A976 andA977, which expand uniformly in the radial direction, areFig.4 Eigenfunctions of m = 1 mode in (a)core-localized modes near the origin, (b)core-localized modes toward imaginary axis,(c) edge-localized modes, (d) resonantsurface symmetric modes, and (e)non-localized modes. Here, K= 1.0 x 10-5,Q= 0, and imax = 512.
short-wavelength oscillations.
The eigenfunctions of B1, B2, and B3 on the 2nd
branch toward the imaginary axis are shown in Fig.4(c).
It is found that these modes localize in the edge plasma
region, differently from the modes on the 1st branch
toward the imaginary axis. It is also found that, as
departing from the imaginary axis, the eigenfunction
comes to have short-wavelength oscillations, similarly tothe modes on the 1st branch toward the imaginary axis.
The eigenfunction of D6, F1, F2, and B11 on theconnecting branch are shown in Fig.4(d). It is found thatF1 and F2 are generally symmetrical about the resonantsurface of the m = 1 mode, that is, almost zero at theresonant surface and finite fluctuation in both core andedge plasma regions. D6 has short-wavelength
oscillations in the core plasma region, and on the other
hand, B12 has short-wavelength oscillations in the edge
plasma region.
The eigenfunction of B14 on the branch toward thereal axis is shown in Fig.4(e), with A230 and A231 of thesubstantially damping modes. It is found that these threemodes do not localize but expand uniformly in the radial
direction with short-wavelength oscillations.

3.3 Eigenfunction of m = 2 mode

In order to understand the characteristics ofeigenvalue and eigenfunction of a different resistive571

T. Matsumoto and S. Tokuda, Eigenvalue Spectrum of MHD Modes in Cylindrical Tokamak Plasmas with Small Resistivity

(a)

(b)

(a)

D1
D2
D3

D11
D12
A0

1.0 10 -1

Eigenvalue Distribution of m=1 Modeeta0=1e-4, vis0=0.0e-0, imax=512

K = 1.0e-4(c)Imaginary

8.0 10 -2

(d)
B1
B2
B3
B4
B5

D7
F1
F2
B25

6.0 10 -2

4.0 10 -2

2.0 10 -2

0.0 10 0-5.0 10 -2 -4.0 10 -2 -3.0 10 -2 -2.0 10 -2 -1.0 10 -2 0.0 10 0 1.0 10 -2

Real

(e)

(b)

B33
A976
A977

1.0 10 -1

Eigenvalue Distribution of m=1 Modeeta0=2e-5, vis0=0.0e-0, imax=512

K = 2.0e-5

Imaginary

8.0 10 -2

6.0 10 -2

Fig.6 Eigenfunctions of m = 2 mode in (a)core-localized modes near the origin, (b)edge-localized modes toward imaginary axis,(c) core-localized modes, (d) resonant surface
symmetric modes, and (e) non-localized
modes. Here, K= 1.0 x 10-5, Q= 0, and imax =2048.

4.0 10 -2

2.0 10 -2

0.0 10 0-5.0 10 -2 -4.0 10 -2 -3.0 10 -2 -2.0 10 -2 -1.0 10 -2 0.0 10 0 1.0 10 -2

Real

(c)

1.0 10 -1

Eigenvalue Distribution of m=1 Modeeta0=1.5e-6, vis0=0.0e-0, imax=512

K = 1.5e-6

similar to B33, A976 and A977 of the m = 1 mode asshown in Fig.4(e)Imaginary

8.0 10 -2

3.4 Dependency on plasma resistivity.

Next, in order to investigate the dependency of
eigenvalue distribution on the plasma resistivity, theeigenvalue problem of the m = 1 mode is solved for
different plasma resistivity. Figure 7 shows thedependency of the eigenvalues distribution of the m = 1
mode on the plasma resistivity. It is confirmed that, for awide regime of the resistivity K from 1.0 x 10-4 to 1.5 x
10-6, the shape and location of these curves (branches) isalmost independent of the resistivity, although the densityof the eigenvalues on these curves increases.
However, for further small resistivity less than 1.5 x
10-6, it is found that the characteristics of eigenvalue
distribution and eigenfunction is changed even for small
reduction of resistivity, as shown in Fig.7(d) where the
distribution structure is deformed with bifurcation ofspectral curve.
The eigenvalue distribution of the m = 1 mode isalso investigated with much larger radial grid number up
to 16384, and this deformation of the eigenvalue
distribution is also confirmed.
Moreover, for the much smaller resistivity, less than
1.0 x 10-6, further deformation with bifurcation isobserved. We are investigating this behavior of the

6.0 10 -2

4.0 10 -2

2.0 10 -2

0.0 10 0-5.0 10 -2 -4.0 10 -2 -3.0 10 -2 -2.0 10 -2 -1.0 10 -2 0.0 10 0 1.0 10 -2

Real

(d)

1.0 10 -1

Eigenvalue Distribution of m=1 Modeeta0=1e-6, vis0=0.0e-0, imax=1024

K = 1.0e-6

Imaginary

8.0 10 -2

6.0 10 -2

4.0 10 -2

2.0 10 -2

0.0 10 0-5.0 10 -2 -4.0 10 -2 -3.0 10 -2 -2.0 10 -2 -1.0 10 -2 0.0 10 0 1.0 10 -2

Real

Fig.7 Eigenvalue distribution of m = 1 mode for (a)K= 1.0 x 10-4, (b) K= 2.0 x 10-5, (c)K= 1.5 x
10-6, and (d)K= 1.0 x 10-6. Here, Q= 0, andimax is changed from 512 to 2048.

572

T. Matsumoto and S. Tokuda, Eigenvalue Spectrum of MHD Modes in Cylindrical Tokamak Plasmas with Small Resistivity

eigenvalue distribution in detail by massively parallel
computation.

4. Summary

For understandings of the MHD spectrum in acylindrical tokamak plasma the characteristics of the m = 1
and m = 2 resistive MHD modes are investigated bysolving an eigenvalue problem based on the reduced MHD
model. The eigenvalue problem code was developed andwas successfully parallelized for much small plasma
resistivity.
The eigenvalues and eigenfunctions are clarified indetail, in particular for the structure of eigenvalue
distribution which discrete eigenvalue comprise. It is
confirmed that eigenvalue distribution is almost
independent of the resistivity larger than 10-6.
However, it is found that eigenvalue distribution is
topologically deformed by reducing the resistivity to lessthan about 10-6. It is also found that the deformation is very
sensitive to small change of the resistivity in this small
resistivity region.
As a future work, the WKBJ method will be applied
to understand the strange phenomenon observed in Fig.
7(d). It will also be confirmed the complex spectrum issensitive to the method used to reduce the equations [4].
The eigenvalue and eigenfunction in the inner resistive
layer will be also analyzed and will be compared withthose obtained by the eigenvalue problem based on theresistive MHD equations.

Acknowledgements

The authors would like to express their gratitude toDrs. M. Kikuchi, H. Ninomiya and Y. Miura (JAEA) for
their useful suggestions and encouragement.
The authors would like to thank Drs. M.Yagi, N.
Miyato (JAEA) and Y. Kishimoto (Kyoto Univ.) for theirhelpful comments.
This study was partly supported by Grant-in-Aids
from the Japan Society for the Promotion of Science(No.19560832).

References
[1]
[2]
[3]
[4]

W. Kerner, K. Lerbinger and K. Riedel, Phys. Fluids, 29,2975 (1986).
B. van der Holst, A.J. Bellen, J.P. Goedbloed, et al., Phys.
Plasmas, 6, 1554 (1999).
R. L. Dewar and B. Davies, J. Plasma Phys. 32, 443 (1984)
B. F. McMillan, R. L. Dewar, and R. G. Storer, Plasma
Phys. Control. Fusion 46, 1027 (2004).

573