PHYSICAL REVIEW E 70, 066409 (2004)

Statistical characterization of the interchange-instability spectrum of a separable
ideal-magnetohydrodynamic model system
1

R. L. Dewar,1,2,* T. Tatsuno,3,2 Z. Yoshida,2 C. Nührenberg,4 and B. F. McMillan1

Department of Theoretical Physics and Plasma Research Laboratory, Research School of Physical Sciences and Engineering,
The Australian National University, ACT 0200, Australia
2
Graduate School of Frontier Sciences, University of Tokyo, 5-1-5 Kashiwanoha, Kashiwa-shi, Chiba 277-8651, Japan
3
Institute for Research in Electronics and Applied Physics, University of Maryland, College Park, Maryland 20742-3511, USA
4
Max-Planck-Institut für Plasma Physik, Teilinstitut Greifswald, D-17491, Germany(Received 3 May 2004; published 27 December 2004)
A Suydam-unstable circular cylinder of plasma with periodic boundary conditions in the axial direction is
studied within the approximation of linearized ideal magnetohydrodynamics (MHD). The normal mode equations are completely separable, so both the toroidal Fourier harmonic index n and the poloidal index m aregood quantum numbers. The full spectrum of eigenvalues in the range 1 艋 m 艋 mmax is analyzed quantitatively,
using asymptotics for large m, numerics for all m, and graphics for qualitative understanding. The density ofeigenvalues scales like m2max as mmax → ⬁. Because finite-m corrections scale as 1 / m2max, their inclusion isessential in order to obtain the correct statistics for the distribution of eigenvalues. Near the largest growth rate,
only a single radial eigenmode contributes to the spectrum, so the eigenvalues there depend only on m and n
as in a two-dimensional system. However, unlike the generic separable two-dimensional system, the statisticsof the ideal-MHD spectrum departs somewhat from the Poisson distribution, even for arbitrarily large mmax.
This departure from Poissonian statistics may be understood qualitatively from the nature of the distribution ofrational numbers in the rotational transform profile.
DOI: 10.1103/PhysRevE.70.066409

PACS number(s): 52.35.Bj, 05.45.Mt

I. INTRODUCTION

complicated geometries in the predesign phase.
For this purpose, a substantial investment in effort has
been expended on developing numerical matrix eigenvalue
programs, such as the three-dimensional TERPSICHORE [8]and CAS3D [9] codes. These solve the MHD wave equations
for perturbations about static equilibria, so that the eigenvalue ␻2 is real due to the Hermiticity (self-adjointness [10])of the linearized force and kinetic energy operators. They use
finite-element or finite-difference methods to convert theinfinite-dimensional PDE eigenvalue problem to an approximating finite-dimensional matrix problem. An alternative approach is to use local analysis using the ballooning representation and to attempt semiclassical quantization to estimate
the global spectrum [11–13].
In order properly to verify the convergence of these codesin three-dimensional geometry, it is essential to understandthe nature of the spectrum—if it is quantum-chaotic, thenconvergence of individual eigenvalues cannot be expected
and a statistical description must be used [14–17].
This is perhaps of most importance in understanding thespectrum in three-dimensional magnetic confinement geometries, in particular the various stellarator experiments currently running or under construction. These devices arecalled three-dimensional because they possess no continuous
geometrical symmetries, and thus there is no separation ofvariables to reduce the dimensionality of the eigenvalue
problem. It has been shown [18] that the semiclassical limit(a Hamiltonian ray tracing problem) for ballooning instabilities in such geometries may be strongly chaotic because
there are no ignorable coordinates in the ray Hamiltonian.
However, the present paper discusses the opposite limit, asystem with a sufficient number of symmetries to make the

The general aim of this paper is to compare and contrast
the spectrum of eigenvalues in typical integrable wave systems (e.g., waves in a rectangular cavity [1]) with the spectrum of instabilities in a cylindrical plasma within the ideal
magnetohydrodynamics (MHD) approximation. This is a
first step in understanding the spectral problem in the complex three-dimensional geometry of the class of magnetic
confinement fusion experiments known as stellarators [2].
In ideal MHD the spectrum of the frequencies, ␻, of normal modes of displacements about a toroidal equilibrium isdifficult to characterize mathematically because the linearized force operator is not compact [3]. In addition to a point(discrete) spectrum of unstable modes 共␻2 ⬍ 0兲, there are the
Alfvén and slow-magnetosonic continuous spectra on thestable side of the origin 共␻2 ⬎ 0兲 and the possibility of dense
sets of accumulation points on the unstable side. In mathematical spectral theory, the stable continua and unstable accumulation “continua” [4] are characterized [5] as belonging
to the essential spectrum. (For a self-adjoint operator L, theessential spectrum is the set of ␭ values for which the rangeof L − ␭ is not a closed set and/or the dimensionality of thenull space of L − ␭ is infinite.)
There is experimental evidence that ideal MHD is relevant in interpreting experimental results [6,7], but perhaps
the greatest virtue of ideal MHD in fusion plasma physics is
its mathematical tractability as a first-cut model for assessing
the stability of proposed fusion-relevant experiments with

*Electronic address: robert.dewar@anu.edu.au
1539-3755/2004/70(6)/066409(11)/$22.50

066409-1©2004 The American Physical Society

PHYSICAL REVIEW E 70, 066409 (2004)

DEWAR et al.

ray Hamiltonian integrable and the eigenvalue problem separable. The geometry is the circular cylinder, periodic in the z
direction to make it topologically toroidal—we shall refer tothe z direction as the toroidal direction and the azimuthal, ␪direction as the poloidal direction. The study of this separable system will provide a baseline for comparison with thethree-dimensional toroidal case in future work. The overall
goal of the paper is to determine if the ideal-MHD spectrum
falls within the same universality class as that of typicalwaves in separable geometries or, if not, what might cause it
to differ.
Berry and Tabor [19] show that the distribution function
P共s兲 for the spacing of adjacent energy levels (suitably
scaled) in a generic separable quantum system with morethan one degree of freedom is exp共−s兲, as for a Poissonprocess with levels distributed at random. They also showthat the spectrum of uncoupled quantum oscillators is nongeneric even when the frequency ratios are not commensurate, in which case P共s兲 peaks about a nonzero value of s (asalso occurs in nonintegrable, chaotic systems—the “level repulsion” effect). A more surprising departure from the Poisson distribution was found by Casati et al. [1] for waves in arectangular box with irrational aspect ratio, but the departure
was very small. Level spacing statistics are discussed also inthe standard monographs on quantum chaos [14–17].
In contrast with quantum mechanics, where the continuous spectrum arises from the unboundedness of configuration
space, the ideal-MHD essential spectrum arises from the unboundedness of Fourier space—there is no minimum wavelength. This is an unphysical artifact of the ideal MHD
model because, in reality, low-frequency instabilities with兩k⬜兩 much greater than the ion Larmor radius, ai, cannot exist(where k⬜ is the projection of the local wave vector into theplane perpendicular to the magnetic field B). Indeed, ideal
MHD breaks down in various ways at large 兩k⬜兩, with dissipative and drift effects coming into play.
In this paper, we do not attempt to model finite-Larmorradius stabilization, but instead simply restrict the poloidal
mode spectrum to m 艋 mmax and study the scaling of thespectrum at large mmax. The nature of the dispersion relation
is such that the toroidal mode numbers n relevant to thespectrum are also restricted. In a matrix eigenvalue codesuch as CAS3D or TERPSICHORE, our procedure corresponds tousing an arbitrarily fine radial mesh but truncating the toroidal and poloidal basis set.
The eigenvalue equation for a reduced MHD model ofa stellarator is presented in Sec. II. We study a plasma inwhich the Suydam criterion [20] for the stability of interchange modes is violated, so the number of unstable modes
is infinite.
Section III is devoted to developing an understanding of
the dependence (the dispersion relation) of the eigenvalues
on the radial, poloidal, and toroidal mode numbers, l, m, andn, respectively. As m and n approach infinity, keeping ␮⬅ n / m fixed, the growth-rate eigenvalues asymptote to a
constant, the Suydam growth rate, depending only on ␮ andthe radial mode number l. We use a combination of perturbation expansion in 1 / m and numerical solution of the eigenvalue equation using a new transformation to Schrödinger
form that is applicable over the whole range of m, from O共1兲to ⬁. This generalizes the approach of Cheremhykh andRevenchuk [21], which was limited to the m = ⬁ Suydam
eigenvalue problem. We compare some of the asymptotic
results in [21] with our numerical solutions. Our perturbation
expansion shows that the correction to the Suydam limit goesas 1 / m2. Contrary to usual experience [22], our numericalsolutions show that the growth rates do not always approach
the Suydam values from below as m → ⬁.
In Sec. IV, we examine the part of the spectrum involving the most unstable modes, which is essentially twodimensional because only the lowest-order radial mode,l = 0, contributes. We relate the considerable amount of structure observed in the spectrum to the Farey sequences of rational values of the rotational transform (winding number) of
the equilibrium magnetic field. Low-order rationals have associated eigenvalue sequences giving a regular distributionof eigenvalues locally more like the spectrum of a onedimensional system than a two-dimensional one.
In Sec. V, we derive the analog of the Weyl formula forthe average density of states, including an asymptotic analysis of the large-l limit. In Sec. VI, we show level spacing
distributions P共s兲. Since we are interested in large m, we firsttry approximating the eigenvalues by their corresponding
asymptotic Suydam limit. This gives a very singular distribution with a ␦-function-like spike at the origin [23] due tothe extremely degenerate nature of the spectrum in this approximation. By contrast, the distribution for the exact spectrum has no spike at the origin, showing that the small 1 / m2
corrections break the degeneracy sufficiently to completelychange the statistics.
We examine the statistics for the l = 0 and l = 1 spectra,
both individually and combined (in the low-growth-rate region where they overlap). We have examined sufficiently
large data sets to show convincingly that the statistical distributions are not Poissonian, though that of the combinedl = 0 and l = 1 spectrum is closest. We also split the l = 0 spectrum into two halves to remove overlap of spectra arisingfrom different parts of the plasma. These split spectra exhibit
a much more dramatic departure from Poisson statistics,
showing that the ideal-MHD interchange spectrum is indeed
nongeneric in the sense of Berry and Tabor [19].
II. CHOICE OF MODEL EIGENVALUE EQUATION

The grand context of this paper is the three-dimensional
linearized ideal MHD problem—to solve, under appropriate
boundary conditions, the equation of motion

␳⳵2t ␰ = F · ␰

共1兲

for small displacements ␰共r , t兲 of the MHD fluid about astatic equilibrium state, where ␳共r兲 is the equilibrium massdensity, r is position, t is time, and F is a Hermitian linearized force operator [10] under the inner product 兰d3x␰ * · F · ␰and suitable boundary conditions. (Superscript * denotes
complex conjugation—we can take ␰ to be complex because
all the coefficients in F are real, so the real and imaginaryparts of ␰ obey the same equation.)
Most modern magnetic confinement fusion experiments,
in particular tokamaks and stellarators, are toroidal. Though

066409-2

STATISTICAL CHARACTERIZATION OF THE …

PHYSICAL REVIEW E 70, 066409 (2004)

not guaranteed for arbitrary three-dimensional systems, theequilibrium magnetic field B共r兲 is normally assumed to be
integrable in the sense that all field lines lie on invariant tori
(magnetic surfaces) nested about a single closed field line(the magnetic axis). Within each toroidal magnetic surface, anatural angular coordinate system is set up, with the poloidal
angle ␪ increasing by 2␲ for each circuit around the short
way and the toroidal angle ␨ increasing by 2␲ for each circuit the long way. Each surface is characterized by a magnetic winding number, the rotational transform –␫, being the
average poloidal rotation of a field line per toroidal circuit,
具d␪ / d␨典, over an infinite number of circuits. (In tokamak
physics the inverse, q ⬅ 1 /–␫, is normally used as the rotation
number.)In this paper, we study an effectively circular-cylindrical
MHD equilibrium, using cylindrical coordinates such thatthe magnetic axis coincides with the z axis, made topologically toroidal by periodic boundary conditions. Thus z andthe toroidal angle ␨ are related through ␨ ⬅ z / R0, where R0 is
the major radius of the toroidal plasma being modeled bythis cylinder. The poloidal angle ␪ is the usual geometriccylindrical angle and the distance r from the magnetic axis
labels the magnetic surfaces (the equilibrium field being
trivially integrable in this case). The plasma edge is at r = a.
In the cylinder there are two ignorable coordinates, ␪ and␨, so the components of ␰ are completely factorizable into
products of functions of the independent variables separately.
In particular, we write the r component asr␰r = exp共im␪兲exp共− in␨兲␸共r兲,

共2兲

where the periodic boundary conditions quantize m and n tointegers and we choose to work with the stream function
␸共r兲 ⬅ r␰r共r兲.
Since the primary motivation of this paper is stellarator
physics, we use the reduced MHD ordering for large-aspect
stellarators [2,24], averaging over helical ripple to reduce toan equivalent cylindrical problem [25,26]. The universality
class should be insensitive to the precise choice of model aslong as it exhibits the behavior typical of MHD instabilities
in a cylindrical plasma, specifically the existence of interchange instabilities and the occurrence of accumulation
points at finite growth rates.
We nondimensionalize by measuring the radius r in unitsof the minor radius of the plasma column, a, and the time t inunits of the poloidal Alfvén time ␶A = R0冑␮0␳ / B0, where B0
is the toroidal magnetic field and ␮0 is the permeability of−1. Defining ␭ ⬅ ␻2, wefree space. Thus ␻ is in units of ␶A
seek the spectrum of ␭ values satisfying the scalar equationL␸ = ␭M ␸

共3兲

under the boundary conditions ␸共0兲 = 0 at the magnetic axisand ␸共1兲 = 0, appropriate to a perfectly conducting wall at theplasma edge. The operators L and M given below are Hermitian under the inner product defined, for arbitrary functions f and g satisfying the boundary conditions, by

具f,g典 ⬅

冕1

f * 共r兲g共r兲rdr.

共4兲

0

The weight factor r in the inner product is a Jacobian factor
coming from d3x = rdrd␪dz.
The operator L is given by
L⬅−

冋

册

1 d
d m2
–␫¨
共n − m–␫兲2r + 2 共n − m–␫兲2 − DS + 共n − m–␫兲 ,
dr r
r dr
m
共5兲

where the Suydam stability parameter DS is
DS ⬅ −

␤0p⬘共r兲⍀⬘共r兲,
2⑀2

共6兲

with ⑀ ⬅ a / R0 Ⰶ 1 the inverse aspect ratio, p共r兲 the plasmapressure normalized to unity at r = 0, ␤0 ⬅ 2␮0 p0 / B20 the ratioof plasma pressure to magnetic pressure at the magnetic axis,and ⍀⬘ the average field line curvature. Here

冉 冕 冊⍀ ⬅ ⑀2N r2–␫ + 2

r–␫dr ,

共7兲

where the rotational transform is produced by helical current
windings making N Ⰷ 1 turns as ␨ goes from 0 to 2␲, ⍀⬘共r兲
giving the averaged field-line curvature. (Note that ⑀ cancels
out in DS.) We use the notation ḟ ⬅ rf ⬘共r兲 for an arbitrary
function f, so –␫˙ ⬅ rd–␫ / dr is a measure of the magnetic shear
and –␫¨ measures the variation of the shear with radius. The
term ⍀ is a measure of the “magnetic hill” [2] that allows
pressure energy to be released by interchanging field lines,
thus driving the interchange instability.
The operator arising from the inertial term in Eq. (1),2=−M ⬅ − ⵜ⬜

1 d d m2r +,r dr dr r2

共8兲

is easily seen to be positive definite under the inner productEq. (4).
We observe some differences between Eq. (3) and thestandard quantum-mechanical eigenvalue problem H␺ = E␺.
One is of course the physical interpretation of theeigenvalue—in quantum mechanics the eigenvalue E ⬅ ប␻ islinear in the frequency because the Schrödinger equation is
first order in time, whereas our eigenvalue ␭ is quadratic in
the frequency because it derives from a classical equation of
motion.
Another difference is that Eq. (3) is a generalized eigenvalue equation because M is not the identity operator. This is
one reason why it is necessary to treat the MHD spectrum
explicitly rather than simply assume it is in the same universality class as standard quantum-mechanical systems.
Just as in ordinary eigenvalue problems, the eigenvalue
spectrum for the generalized eigenvalue problem is real,
and the eigenfunctions ␸i have a generalized orthogonality
property

066409-3

PHYSICAL REVIEW E 70, 066409 (2004)

DEWAR et al.

FIG. 1. (a) The rotational transform –␫ ⬅ 1 / q, defined by Eq. (11),as a function of r in units of the minor radius; (b) the radial profileof the magnetic shear parameter –␫˙ ⬅ r–␫⬘共r兲. In (a), all distinct rational magnetic surfaces –␫ = ␮ ⬅ n / m are shown for m up to 10.

具␸i,M ␸ j典 = ␦i,j ,共9兲

where the normalization has been chosen to make the coefficient of the Kronecker ␦ unity. Here i and j denote members of the set 兵l , m , n其, where l is the radial node number andthe poloidal and toroidal mode numbers m and n, respectively, are defined in Eq. (2). The negative part of the spectrum, ␭ = −␥2 ⬍ 0, corresponds to instabilities growing exponentially with growth rate ␥.
Equation (3) is very similar to the normal mode equation
analyzed in the early work on the interchange growth rate instellarators by Kulsrud [25]. However, unlike this and mostother MHD studies, we are concerned not with finding thehighest growth rate, but in characterizing the complete set of
unstable eigenvalues.
III. INTERCHANGE SPECTRUM

In this section, we discuss the standard unregularized
ideal MHD spectrum. It is well known that for ␭ ⬎ 0, thespectrum consists of the Alfvén continuum (the slowmagnetosonic continuum being removed in reduced MHD[27]). On the unstable side of the spectrum, ␭ ⬍ 0, it is also
known that there is an infinity of eigenvalues provided theSuydam interchange instability criterion [20]
G⬅

DS 1
⬎
4
–␫2

共10兲

is satisfied over some range of r in the interval (0, 1), but thedetails of the spectrum do not appear to have been published
before.

FIG. 2. (a) The nondimensional pressure profile p共r兲, Eq. (13),used in this paper and (b) the Suydam criterion parameter G共r兲,defined in Eq. (10) (solid line), and the instability threshold 1 / 4(dashed line), showing nearly all the plasma is interchange unstable.

–␫0 = 0.45, –␫0 = 0.2.
Given a rational fraction ␮ = n␮ / m␮ in the interval
关–␫共0兲 ,–␫共a兲兴 (where n␮ and m␮ are mutually prime), there is aunique radius r␮ such that –␫共r␮兲 = ␮. Any pair of integers
共m , n兲␮,␯ ⬅ 共␯m␮ , ␯n␮兲, ␯ = 1 , 2 , 3 , . . . satisfies the resonanceconditionn␮,␯ − m␮,␯–␫共r␮兲 = 0.

共12兲For example, the set of rationals with 1 ⬍ m 艋 10 in
the interval of –␫ shown in Fig. 1(a) is 兵␮其= 兵1 / 2 , 5 / 9 , 4 / 7 , 3 / 5 , 5 / 8其, as shown in the figure.
To understand the global spectrum, we also need to know
something about the pressure profile. In this paper, we use abroad pressure profile that is sufficiently flat near the magnetic axis that the Suydam instability parameter G defined in
Eq. (10) goes to zero at the magnetic axis, and for which p⬘vanishes at the plasma edge
p共r兲 = 1 − 6r5 + 5r6 .

共13兲(Recall from Sec. II that p is in units of the pressure at themagnetic axis.) This profile is shown in Fig. 2(a) and theresulting G profile in Fig. 2(b).
B. High m and n

In this subsection, we choose a particular rational surfacer␮ and restrict attention to pairs 共m , n兲 from the set
兵共m , n兲␮,␯ 兩 ␯ = 1 , 2 , 3 , . . . 其 satisfying the condition Eq. (12).
Defining a scaled radial variable x ⬅ m共r − r␮兲 / r␮, we expand all quantities in inverse powers of m,
L⬅

A. Profiles

m2 共0兲
共L + m−1L共1兲 + m−2L共2兲 + ¯ 兲,r␮2

m2 共0兲
共M + m−1M 共1兲 + m−2M 共2兲 + ¯ 兲.
r␮2Interchange instabilities occur only for values of m and nsuch that n − m–␫ vanishes (or at least can be made very small[26]) and therefore it is important to know something about
the function –␫共r兲. The typical profile of –␫共r兲 in a stellarator ismonotonically increasing in the interval 关0 , a兴 and we shallassume this to be the case here (though it is not always truein modern stellarators). For the numerical work in this paper,
we use a parabolic profile

Also, ␭ = ␭共0兲 + m−1␭共1兲 + m−2␭共2兲, and similarly for ␸. The detailed expressions are given in Appendix A.
We then solve Eq. (3) by equating the LHS to zero order
by order. At O共m0兲, as found by Kulsrud [25], we have thegeneralized eigenvalue equation–␫ = –␫0 + –␫2r2

L共0兲␸共0兲 = 0,

M⬅

共11兲

as illustrated in Fig. 1(a). In this and all subsequent plots,

where

066409-4

共14兲

共15兲

STATISTICAL CHARACTERIZATION OF THE …

L共0兲 ⬅ L共0兲 − ␭共0兲M 共0兲 = −PHYSICAL REVIEW E 70, 066409 (2004)

d 2 2
d共–␫ x − ␭共0兲兲 + –␫2x2 − ␭共0兲 − DS
dx
dx
共16兲

with –␫˙ and DS evaluated at r␮. For ␭共0兲 ⬍ 0, Eq. (15) can besolved to give a square-integrable eigenfunction under theboundary conditions ␸共0兲 → 0 as r → ± ⬁ when ␭共0兲 is one of
the eigenvalues ␭␮,l, l = 0 , 1 , 2 , . . ., denoting the number ofradial nodes of the eigenfunction ␸共0兲 = ␸␮,l. Note that ␭␮,ldepends only on ␮ = n / m and is otherwise independent of themagnitude of m and n. We assume that the ␸␮,l共r兲, whencombined with the continuum generalized eigenfunctions for␭共0兲 ⬎ 0, form a complete set.
1. Suydam approximation

The leading term in the expansion of the eigenvalue in
1 / m gives the growth rate in the limit m → ⬁, known as theSuydam approximation. Restricting attention to unstablemodes, so that ␥ ⬅ 共−␭兲1/2 is real, we transform Eq. (15) to
the Schrödinger form [21]d 2␺+ Q共␩兲␺ = 0,
d␩2

共17兲

where
Q = Q0共␩兩␥, ␮兲 ⬅ G − 41 − 41 sech2 ␩ − ⌫2 cosh2 ␩ , 共18兲with G ⬅ G共r␮兲 defined as in Eq. (10), ⌫ ⬅ ␥ /–␫˙共r␮兲, ␩ defined
through x ⬅ ␥ sinh ␩ /–˙␫共r 兲, and ␺ ⬅ 共cosh ␩兲1/2␸共x兲. [In Ref.
␮[21], Eq. (18) is derived from the Fourier transform of Eq.
(15), but we can also use the real-space version as the equation shares with the quantum oscillator the remarkable property of having the same general form in both Fourier spaceand real space.]
Cheremhykh and Revenchuk [21] (CR) have made an extensive study of the eigenvalues of Eq. (17) using the semiclassical quantization condition

冖

Q0共␩兲1/2d␩ = 共2l + 1兲␲ ,

共19兲

␺
which
follows
from
the
WKB
ansatz= A共␩兲exp± i 兰 Q1/2d␩.
CRderiveseveralapproximations,0
useful in appropriate limits, improving on the earlier result ofKulsrud [25]. In this paper, we use two of their results tocompare with numerical solutions of Eq. (17). The first is Eq.
(4.5) of [21],

⌫⬇

冤

4␴
exp −
e

冉 冊

冥

1
␲
2
1
− 2 ,
2␴
4␴

l+

FIG. 3. (a) m = ⬁ eigenfunctions for the l = 0 (solid line), l = 1(short dashes), and l = 2 (short and long dashes) modes at ␮ = 1 / 2,
arbitrary normalization. (b) Growth rates ␥ (in units of the inverse
poloidal Alfvén time) vs resonant –␫ ⬅ ␮. Dashed lines show approximations Eq. (21) (for l = 0) and Eq. (20) (for l = 1 and 2).

As is seen from Fig. 3, Eq. (21) gives a remarkably good
approximation to the growth rate of the most unstable radial
eigenmode, l = 0, and Eq. (20) gives a good approximation
for the higher-l modes (the semiclassical quantization being
strictly justifiable only for large l). The growth-rate maxima
for each l occur close to the maximum of G (and hence ⌫),but not exactly owing to the –␫˙ factor in the definition ⌫⬅ ␥ /–␫˙共r␮兲.
From Eq. (21) we see that, provided the Suydam criterionG ⬎ 1 / 4 is satisfied, there is an infinity of growth rate eigenvalues accumulating exponentially toward the origin fromabove (so the ␭ values accumulate from below) in the limit
l → ⬁.
Perhaps less widely appreciated (because m and n arenormally taken to be fixed) is the fact that there is also apoint of accumulation of the eigenvalues of Eq. (3) at each␭␮,l as mmax → ⬁ with l fixed. To break the degeneracy of␭共0兲, we must proceed further with the expansion in 1 / m.
2. 1 / m2 corrections

Proceeding with the expansion Eq. (14), the calculation
goes through much as in standard time-independent quantum
perturbation theory [28], e.g..
The lowest-order eigenvalues and eigenfunctions are, as
found in Sec. III B, ␭共0兲 = ␭␮,l and ␸共0兲 = ␸␮,l共x兲, respectively.
The O共1 / m兲 correction, ␭共1兲, vanishes identically from parity
considerations—␸␮,l共x兲 is either an even or odd function so
its contribution to the matrix elements of L共1兲 and M 共1兲 between ␸共0兲 and ␸共0兲 is even. On the other hand, L共1兲 and M 共1兲
are odd, so ␭共1兲 ⬅ 0. (This contrasts with the finite-aspectratio toroidal case where toroidal coupling of Fourier harmonics of different m to form ballooning modes leads to anonvanishing 1 / n correction [29,30].)
The first nonvanishing correction term is thus

共20兲␭共2兲 = 具␮,l兩L共2兲兩␮,l典 −which [combining the criteria given in CR’s Eqs. (4.4) and(4.12)] is applicable when ␴ ⬅ 共G − 1 / 4兲1/2 Ⰷ 1 / 2. The second CR result we use is their Eq. (4.7),⌫2 ⬇

G − 共2l + 1兲G1/2,
1 + 共4G兲−1applicable when G ⲏ ⌫2 Ⰷ 1.

共21兲

兺

l⬘⫽l

具␮,l兩L共1兲兩␮,l⬘典具␮,l⬘兩L共1兲兩␮,l典
␭␮,l⬘ − ␭␮,l,共22兲
where the sum over l⬘ is taken to include an integration overthe continuum. The operators L共i兲 ⬅ L共i兲 − ␭␮,lM 共i兲 are thehigher-order generalizations of L共0兲, defined by Eq. (16). Them = ⬁ matrix elements of any operator F are defined by

066409-5

PHYSICAL REVIEW E 70, 066409 (2004)

DEWAR et al.

FIG. 4. Nondimensional growth rates ␥l,m,n ⬅ 共−␭l,m,n兲1/2 vs m−2
for l = 0 and 1, found by numerical solution of Eq. (3) (Sec. III C):(a) n / m = 1 / 2 共m = 2 , 4 , 6 , . . . 兲 and (b) n / m = 11/ 19 共m= 19, 38, 57, . . . 兲. At high m, the dependence becomes linear, in
qualitative agreement with Sec. III B 2.

具␮,l⬘兩F兩␮,l⬙典 ⬅

冕

⬁

−⬁

*

␸␮,l⬘共x兲F␸␮,l⬙共x兲dx,

共23兲

with the eigenfunctions ␸␮,l共x兲 being normalized so that具␮ , l兩M 共0兲兩␮ , l典 = 1. Note that, with the operators L and M
defined as in Eqs. (5) and (8), L共i兲 is Hermitian under the
inner product used in Eq. (23) only for i = 0. However, it canbe made Hermitian at arbitrary order by the redefinitions
L 哫 rL and M 哫 rM, which puts the eigenvalue equation into
Sturm-Liouville form.
As in quantum mechanics [28] e.g., if L共1兲 is Hermitian
the contribution of the second term on the right-hand side ofEq. (22) is always negative for the lowest eigenvalue, ␭共0兲0 ,共0兲⬎ 0. However, in ideal MHD a positive
because ␭l − ␭共0兲l⬘contribution from the first term usually dominates and theinfinite-m mode is most unstable [22]. As seen in Fig. 4, this共2兲
is not always the case: ␭l=0
is negative for ␮ = 1 / 2 but positive for ␮ = 11/ 19= 0.578 947. . . .
The latter value of ␮ is very close to the value giving the
global maximum Suydam growth rate (see Fig. 3). Thus, inthe special case studied here and in accordance with conventional wisdom, the global maximum interchange growthrate occurs at m = ⬁. Both these results are intuitively
reasonable—the eigenfunctions become increasingly localized as m → ⬁, so the highest growth rate is obtained by
localizing in the “most unstable” region of the plasma, where–␫ ⬇ 11/ 19. On the other hand, modes which localize in “lessunstable” regions as m → ⬁ can achieve a higher growth rate
at finite values of m because their more extended finite-m
eigenfunctions overlap the more unstable region and tap into
the free energy from the pressure gradient in this region.
Since the eigenvalues approach ␭␮,l as 1 / m2, there is an
infinity of modes in the neighborhood of each ␭␮,l in thelimit mmax → ⬁. That is, they are finite-growth-rate accumulation points of the complete spectrum. Because the rationals
␮ are dense on the interval (–␫共r1兲 ,–␫共r2兲), where 共r1 , r2兲 is theregion in which the Suydam instability criterion is satisfied,and because ␭␮,l in general depends continuously on ␮, the2, 0兲 densely.
accumulation points ␭␮,l fill the interval 共−␥maxThis is the part of the unstable spectrum called the “accumulation continuum” by Spies and Tataronis [4], though “accumulation essential spectrum” might be better terminology
mathematically.
C. Finite m and n

In order to calculate arbitrarily high- or low-m eigenfunctions, we generalize the transformation in Sec. III B 1 by the

FIG. 5. Nondimensional growth rates vs radial node number l:(a) found by numerical solution of Eq. (3) (Sec. III C) for 共m , n兲
= 共2 , 1兲 (diamonds) and (8, 4) (triangles), normalized to the ␮= 1 / 2, infinite-m results (filled boxes); (b) ␮ = 1 / 2, infinite-m results(points) and asymptotic result Eq. (20) (dashed line).

change of variable from r to a new independent variable ␩such thatm–␫共r兲 − n ⬅ ␥ sinh ␩ ,共24兲

and a new dependent variable ␺共␩兲 such that␸ ⬅ 共–␫˙ cosh ␩兲−1/2␺共␩兲,共25兲

so that Eq. (3) becomes the Schrödinger equation (17), butwith Q0 replaced by
Q ⬅ G共␩兲 −+1 sech2 ␩ ␥2−−
cosh2 ␩
4
4
–␫˙2

冉 冊

tanh ␩ d–␫˙ 1 d2–␫˙
1 d–␫˙−2 +2–␫˙ d␩ 2–␫˙ d␩
4–␫˙2 d␩

2

,

共26兲

where –␫˙ ⬅ rd–␫ / dr is as defined in previous sections, but expressed in terms of ␩.
Differentiating Eq. (24), we finddr ␥ cosh ␩ r
.
=
m
d␩
–˙␫

共27兲

Thus, in the large-m limit, equilibrium parameters such as G
and –␫˙ are slowly varying functions of ␩, e.g., d–␫˙ / d␩
= O共1 / m兲 and d2–␫˙ / d␩2 = O共1 / m2兲. Comparing Eq. (26) with
Eq. (18), we see that, to leading order in 1 / m, Q = Q0 but
with –␫˙ now a slow variable rather than a strict constant.
With the simple form for –␫, Eq. (11), assumed in thispaper, Eq. (24) is easily inverted to give r共␩兲, and also acancellation occurs between the terms –␫˙⬘共␩兲tanh ␩ / 2–␫˙ and–␫˙⬙共␩兲 / 2–␫˙, so that the exact Q is not much more complicatedthan Q0. The eigenvalues in Figs. 4 and 5 were computed byintegrating Eq. (17) with Q0 replaced by the exact Q andwith the appropriate finite boundary conditions. Low-m results were checked against those from an untransformed
shooting code. The dashed lines represent the results of scansthrough unquantized, noninteger values of m to show thesmooth, but not necessarily monotone, functional dependence of ␥ on m

066409-6

STATISTICAL CHARACTERIZATION OF THE …

PHYSICAL REVIEW E 70, 066409 (2004)

FIG. 6. (a) Lattice of quantum numbers on which the part of thespectrum between the “ground state” and the threshold for the entryof the l = 1 mode is defined, and the unbounded contours of constanteigenvalue. (b) The same, in ␮ ⬅ n / m, m space.
IV. l = 0 SPECTRUM

The most unstable modes are those with radial node number l = 0. Thus we first consider the set S0 ⬅ 兵␭0,m,n 兩 1 艋 m
艋 mmax , m␮min ⬍ n ⬍ m␮max其, where m and n are integers and␮min and ␮max are chosen to give the desired range of ␥. As
we shall be rescaling the eigenvalues prior to statisticalanalysis, it makes no difference whether we work with thespectrum of growth rates ␥ or the eigenvalues ␭ ⬅ −␥2. However, the latter choice makes it clearer that the analog of thequantum-mechanical ground state is the most rapidly growing mode—denoting the maximum growth rate of the l = 02= −␥20.
mode by ␥0, the minimum ␭ is ␭0 = −␥max
The spectrum is defined on the fanlike subset of the twodimensional quantum-number lattice depicted in Fig. 6(a).
Also shown are contours of constant ␥ (or ␭), regarded as acontinuous function of m and n, which are seen more clearlyin Fig. 6(b). Here we see a striking contrast with more generic systems [19], where the constant-eigenvalue contours
are segments of topological circles enclosing the origin. In
the ideal-MHD case, the contours are topologically hyperbolic, with asymptotes radiating from the origin toward infinity.
An interesting representation of the l = 0 spectrum isshown in Fig. 7. A great deal of structure can be discerned,
determined by the number-theoretic properties of the intervalof ␮ depicted. For instance, focusing on the low-order rational number 4 / 7, we define spectral subsets S0共N / M 兩 4 / 7兲⬅ 兵␭0,m,n 兩 m = M + 7k , n = N + 4k , k = 0 , 1 , 2 , . . . , 关共mmax − N兲 / 7兴其,
where 关x兴 denotes the largest integer 艋x.

FIG. 7. Nondimensional growth-rate eigenvalues of l = 0 modes
near the maximum growth rate vs ␮ ⬅ n / m. The ensemble shown isfor mmax = 100.

These spectral sequences all accumulate toward the same
Suydam eigenvalue ␭4/7,0 as mmax → ⬁ independently of the
choice of M and N. However, the rapidity of this approach issensitive to the choice of M / N. For instance, we see in Fig.
7 the most rapidly converging sequence, S0共4 / 7 兩 4 / 7兲, as aset of points accumulating vertically from below toward theSuydam eigenvalue. Other sequences on either side of
S0共4 / 7 兩 4 / 7兲 approach the accumulation point obliquely andmuch more slowly—for mmax = 100 they visibly have somedistance to go. The sequence immediately to the left of
S0共4 / 7 兩 4 / 7兲 is S0共1 / 2 兩 4 / 7兲, while that to the right is
S0共3 / 5 兩 4 / 7兲, 1 / 2 and 3 / 5 being the immediate neighbors of
4 / 7 in the Farey sequence ([31], p. 300) of order 7 (the firstorder at which 4 / 7 appears), with the ␮ values corresponding to S0共1 / 2 兩 4 / 7兲 and S0共3 / 5 兩 4 / 7兲 providing the immediate neighbors of 4 / 7 in each higher-order Farey sequence.
In discussing the structure of the spectrum, it is useful topartition S0 into two subsets, S−0 and S+0 , as the points are to
the left or right, respectively, of the dashed vertical line
shown in Fig. 7 passing through the point of maximum
growth rate.
The sequences S0共1 / 2 兩 4 / 7兲 and S0共3 / 5 兩 4 / 7兲 accumulate
toward ␭4/7,0, but slower 关O共1 / mmax兲兴 than does S0共4 / 7 兩 4 / 7兲2
兲, from Sec. III B 2]. Thus there is a gap contain[O共1 / mmax
ing ␭4/7,0 within which S0共4 / 7 兩 4 / 7兲 contributes O共mmax兲
points to S−0 , while other sequences contribute at most a set of
O共1兲 points.
Within the gap, the spectrum S−0 is essentially onedimensional, being indexed by the single quantum number k.
In the full spectrum, S0 = S−0 艛 S+0 , O共mmax兲 unrelated eigenvalues from S+0 appear in the gap, making the spectrum locally more random and two-dimensional.
V. WEYL FORMULA

As discussed in Sec. III B 2, the overall maximum growthrate for the l = 0 and 1 modes (and, we assume, for all l)occurs at m = ⬁. Thus the threshold value when a given model first starts contributing to the spectrum is at ␭ = −␥2l , where␥l is the maximum over ␮ of ␥共␮ , l兲. We denote the corresponding value of ␮ by ␮l.
For fixed l and large mmax, the number of eigenvalues
Nl共␮兲 in an interval of n / m between ␮l and ␮ is asymptotically equal to the area in the m, n plane [see Fig. 6(a)] of thetriangle bounded by the lines n = ␮m, n = ␮lm, and m = mmax.
2
That is, Nl共␮兲 ⬃ 21 兩␮ − ␮l兩mmax.
Since contours of constant ␭ (or ␥) asymptote to lines of
constant ␮ as m → ⬁, we can estimate the number of eigenvalues between two values of ␭ (or ␥) by inverting the function ␭␮,l for ␮ and substituting this into the above expression
for Nl共␮兲. The inverse is double-valued: ␮ = ␮+l 共␭兲 ⬎ ␮l and␮−l 共␭兲 ⬍ ␮l. Then the number of eigenvalues between the
ground state ␭␮0,l and ␭ is approximately
1
2.
N̄±l 共␭兲 ⬅ 兩␮±共␭兲 − ␮0兩mmax
2

共28兲

The asymptotic dependence of the total spectrum S⬅ S−0 艛 S+0 艛 S−1 艛 S+1 艛 . . . is thus

066409-7

PHYSICAL REVIEW E 70, 066409 (2004)

DEWAR et al.

冋 冉 冊册

␴共␮兲 = ␴max 1 −␮ − ␮0 2+ O关共␮ − ␮0兲3兴,
⌬␮

共30兲

where 共⌬␮兲2 ⬅ −2␴max / ␴⬙共␮0兲. To leading order all other parameters are evaluated at the maximum point ␮ = ␮0. The
quadratic correction to ␴ need only be retained in the term
involving the expansion parameter l, so, to leading order,
FIG. 8. (a) Eigenvalue sequence number N共␥兲 for the combined
spectral set S0 艛 S1 in the case mmax = 100. The Weyl formula for S0,
N̄−0 共␥兲 + N̄+0 共␥兲, is shown as a dashed line. (b) Closeup of the region
containing eigenvalues associated with ␮ = 3 / 5 and ␮ = 5 / 9.
⬁N̄共␭兲 ⬅

N̄±l 共␭兲.
兺
兺
l=0 ±

冋

␥l ⬇ ␥0 exp −

共31兲

2
兲.
where ␥0 ⬅ 4–␫˙共r␮0兲共␴max / e兲exp共−␲ / 4␴max − 1 / 4␴max
Solving for ␮, we find

冋

共29兲

This is the analog of the Weyl formula [14, p. 258] for the
integral of the smoothed spectral density (“density ofstates”).
Approximating ␭␮,l = ␭l + 21 共⳵2␭␮l,l / ⳵␮2l 兲共␮ − ␮l兲2, we get
␮±l 共␭兲 = ␮l ± 冑2共␭ − ␭l兲1/2 / 共⳵2␭␮l,l / ⳵␮2l 兲1/2. Thus there is asquare-root singularity at each mode threshold.
A comparison between the Weyl formula for S0 and theset of points 兵共␥N , N兲其, where N is the sequence number obtained by sorting the set of l = 0 and l = 1 growth-rate eigenvalues from largest to smallest, is shown in Fig. 8(a), showing excellent agreement above the threshold for S1. The
plotted points may also be regarded as the locations of thesteps in the “staircase plot” of the piecewise-constant integrated density-of-states function N共␥兲, but the scale in thisplot is too coarse to resolve the staircase structure.
A finer-scale plot is shown in Fig. 8(b), in which significant deviations from the Weyl curve are seen in the microstructure. The range shown in Fig. 8(b) is unusual in that it
contains two well-defined accumulation sequences in close
proximity. These are associated with low-order values of ␮
occurring on either side of the growth-rate maximum near
␮ = 11/ 19⬇ 0.579—the sequence associated with ␮ = 5 / 9⬇ 0.556 is in S−0 and the one associated with ␮ = 3 / 5 = 0.6 isin S+0 . There are very few eigenvalues associated with highorder rational values of ␮ in the range shown and the two
low-order sequences present are practically unmixed, eitherwith each other or with eigenvalues associated with unrelated higher-order rational values of ␮. [In fact there is only
one such high-order mode in the region of the accumulation
sequences, ␮ = 51/ 92, the closest approximant to 5 / 9 in theset corresponding to S0共1 / 2 兩 5 / 9兲, which causes the slightjump seen in the ␮ = 3 / 5 sequence.] Also, the wide gap containing no eigenvalues is because the intersection of the gaps
associated with the two low-order rationals is nonempty.
The spectrum near the marginal stability point, ␥ = 0, will
involve the superposition of many branches of radial eigenvalue l. To estimate the asymptotic behavior of N共␥兲 as ␥
→ 0, we use the approximate dispersion relation Eq. (20).
Taking l to be large, we see from Eq. (20) that the Suydam
growth rates ␥l共␮兲 are sharply peaked about the location of
the maximum, ␮0, of G共r␮兲, where ␴共␮兲 ⬅ 共G − 1 / 4兲1/2 isalso a maximum. Thus we can expand ␴共␮兲 about ␮0,

册

␲l
,
2␴共␮兲

␮±l 共␥兲 = ␮0 ± ⌬␮ 1 −llmax共␥兲

册

1/2

,

共32兲

where lmax共␥兲 ⬅ 共2 / ␲兲␴max ln共␥0 / ␥兲. Substituting Eq. (32) in
Eq. (29) and approximating the sum over l by an integral, we
find the leading-order asymptotic behavior of the number of
eigenvalues to beN̄共␥兲 ⬃

4⌬␮
␥02
␴maxmmax
ln ,
3␲
␥

共33兲

which diverges logarithmically as ␥ → 0.
VI. NEAREST-NEIGHBOR STATISTICS

Preparatory to the statistical analysis of eigenvalue spacing, it is standard practice to rescale, or unfold, the eigenvalues so as to make their average separation unity, thus making
possible the comparison of different systems on the samefooting.
We can unfold the spectra by using the Weyl formulas
above, e.g., for ␭i 苸 S±0 we can define rescaled eigenvalues E±i
by
E±i ⬅ N̄±0 共␭i兲.

共34兲

For the set S0 = S+0 艛 S−0 , we can unfold with the combined
Weyl function, 兺±N̄±0 . However, for practical purposes wehave in this section used empirical least-square fits of N共␥兲
to a linear superposition of the basis functions 共␥max − ␥兲1/2,
共␥max − ␥兲, 共␥max − ␥兲3/2, which captures the square-root singularity but avoids having to invert ␥␮,l.
When mmax is large, the great majority of eigenvalues
␭l,m,n are very close to the corresponding m = ⬁ eigenvalue
with the same ␮ ⬅ n / m, ␭␮,l. Thus one might supposethat the statistics of the spectrum are asymptotically the same
⬅ 兵␭n/m,0 兩 1 艋 m 艋 mmax ,as those of an ensemble SSuydam
0m␮min ⬍ n ⬍ m␮max其.
In Fig. 9(a), we show the distribution of nearest-neighbor, and in Fig. 9(b)unfolded eigenvalue spacings for SSuydam
0that for the set S0 with the correct finite-m eigenvalues. It isseen that the two distributions are radically different—eventhough low-order rational values of ␮ are rare and the distribution is coarse-grained, the high-m approximation induces
sufficient extra degeneracy that the Suydam spectrum isdominated by a large, but spurious, ␦-function-like spike at

066409-8

STATISTICAL CHARACTERIZATION OF THE …

PHYSICAL REVIEW E 70, 066409 (2004)

FIG. 9. (a) Nearest-neighbor eigenvalue spacing distribution forusing ␮min = 0.5044, ␮max
the approximate spectral set SSuydam0= 0.6288, mmax = 100 (625 eigenvalues). (b) The same, for the corresponding set of accurate eigenvalues S0.

FIG. 11. (a) Nearest-neighbor eigenvalue spacing distribution
for the spectral set S−0 for mmax = 1000 (37 932 eigenvalues). (b)Nearest-neighbor eigenvalue spacing distribution for the spectral set
S+0 (24 412 eigenvalues).

s = 0. (The range of ␮ used in Fig. 9 corresponds to the rangeof l = 0 growth rates above the maximum l = 1 rate, in which
S0 is the only contributor to the spectrum.)
The reason why finite-m effects are so important, despite
the smallness of the O共1 / m2兲 corrections found in Sec.
III B 2, is seen from the Weyl formula, Eq. (28), whichshows that the average eigenvalue spacing in a set containing all values of n / m within the range of interest scales as−2, which is the same order as the smallest O共1 / m2兲 cormmax
rection within a set containing only n / m = const. Thus in theset of accurate eigenvalues S0 there is a strong interminglingof eigenvalues with different n / m that does not occur in the.
approximate set SSuydam
0This explains why the nearest-neighbor eigenvalue spacing distribution in Fig. 9(b) is much closer to the Poisson
distribution exp共−s兲 obtained for a random distribution of
numbers on the real line, and also predicted for generic separable systems [19], than that in Fig. 9(a). Nevertheless, theset of 625 eigenvalues used in Fig. 9(b) is too small to say
convincingly that the distribution is or is not Poissonian, so
we need to analyze larger data sets to determine how close to
generic the ideal-MHD spectrum is.
A cutoff at mmax = 1000 gives a set S0 containing about
62 254 eigenvalues in the range between the maximum l = 0
growth rate and the maximum l = 1 growth rate. [Note the2approximately mmax
scaling in the size of S0, as predicted by
the Weyl formula, Eq. (28).] In Fig. 10(a), we show thenearest-neighbor distribution for this set. Close examinationof the region near the origin reveals no trace of the spikeseen in Fig. 9(a), not even the tiny spike found by Casati etal. [1] for the spectrum of waves in an incommensurate rectangular box. However, it is clear that the statistics are notexactly Poissonian.

In Fig. 10(b), we show the Dyson-Mehta rigidity parameter ⌬3共L兲 [17, pp. 321-323], defined as the least-squares
deviation of the unfolded eigenvalue staircase N共E兲 from thebest-fitting straight line in an interval of length L. Again, thebehavior is similar to that for the completely random spectrum (Poisson process) in that ⌬3 increases linearly with L,
but the slope is slightly less than the 1 / 15 expected for the
Poisson processs.
In order to understand the departure from Poisson statistics better, we show in Fig. 11 the spacing distribution for thecorresponding sets S−0 and S+0 . The departure from Poisson
statistics is now quite striking.
In Fig. 12(a), we show the spacing distribution for the l= 1 spectrum, which is seen to be very much like the l = 0
spectrum of Fig. 10(a) in its departure from the Poisson distribution. However, we might expect that mixing the l = 0
with the l = 1 spectrum will make the levels appear more
“random” and Fig. 12(b) confirms that the level spacing distribution does indeed become more like the exponential expected for a Poisson process.

FIG. 10. (a) Nearest-neighbor eigenvalue spacing distribution
for the spectral set S0 using ␮min = 0.5044, ␮max = 0.6288, mmax
= 1000 (62,254 eigenvalues). (b) The Dyson-Mehta spectral rigidity
for this set (solid line) compared with that for the Poisson process(dashed line).

FIG. 12. (a) Nearest-neighbor eigenvalue spacing distribution
for the first 72 500 eigenvalues of the l = 1 spectral set S1, mmax
= 1000. (b) Nearest-neighbor eigenvalue spacing distribution for themixed spectral set S0 艛 S1 over the same range of eigenvalues as in(a) (total of 95 000 eigenvalues).

VII. CONCLUSION

We have demonstrated that the statistical nature of theideal-MHD interchange spectrum deviates significantly fromthe random Poisson process of generic separable systems dueto the number-theoretic structure of the eigenvalue distribution. The similarity between the two level-spacing distributions in Fig. 11, which correspond to two different parts of
the rotational transform profile, suggests the possibility that
there may nevertheless be some universality in the statistics.
If so, we have found a new universality class.

066409-9

PHYSICAL REVIEW E 70, 066409 (2004)

DEWAR et al.

The crude regularization used in this paper, simply restricting the poloidal mode numbers to m 艋 mmax, is not very
physical but corresponds closely to what is done in the large
three-dimensional eigenvalue codes CAS3D [9] and TERPSICHORE [8]. Thus, apart from fundamental mathematical interest, the primary motivation of this paper has been the numerical analysis of the three-dimensional ideal-MHD
spectrum as produced by these codes. Preliminary results[23] on an interchange-unstable stellarator test case show
spectra with eigenvalue separation statistics similar to thoseof strongly quantum chaotic systems. However, the results ofthe present paper indicate that some caution should be taken
in interpreting ideal-MHD spectra in terms of conventional
quantum chaos theory because of the radically different nature of the dispersion relation.
In subsequent work it will be important to examine theeffect of the finite Larmor radius on the spectrum. However,this typically makes the problem non-Hermitian and lesseasy to compare with standard quantum chaos theory.
A preliminary attempt [23] to regularize MHD via a
Hermiticity-preserving inverse-Larmor-radius cutoff in 兩k⬜兩found a much more Poissonian eigenvalue-spacing histogram for the l = 0 spectrum than that found in the present
paper. This suggests that finite-Larmor radius regularization
is not only more physical than simple truncation in m and n,
but also makes the spectrum more generic (provided the eigenvalues remain real).

APPENDIX A: 1 / m2 CORRECTIONS

The coefficients of the expansion Eq. (14) are found by
Taylor expansion of the geometric and equilibrium quantities
in Eq. (5),
L共0兲 = −L共1兲 = xd ˙2 2 d ˙2 2
+ –␫ x − DS ,
–␫ x
dx
dx

d ˙2 2 ddd− –␫˙–␫¨x3 + –␫˙共–␫¨ − 3–␫˙兲x3
–␫ x
dx
dx dx
dx

+ 共2DS − ḊS − –␫˙–␫¨兲x,
L共2兲 = − x2⫻冉dd ˙2 2 ddd+ x –˙␫–␫¨x3 + x4
–␫ x
dx
dx
dx
dx dx

冊

–␫˙2 –␫¨2 –␫˙–␫¨ –␫˙ត–␫ d x2− + −+12 4 2 3 dx 2

⫻共5–␫˙–¨␫ − 2–␫˙ត–␫ − –␫¨2 − 6DS + 5ḊS − D̈S兲 + x4⫻

冉冊

71–˙␫2 –␫¨2 7–␫˙–␫¨ –␫˙ត
–␫
+
+ −
12
4
2
3

共A1兲

and Eq. (8),M 共0兲 = −M 共1兲 = x

d2
+ 1,
dx2

d2
d d
− x − 2x,dx2 dx dx

ACKNOWLEDGMENTS

M 共2兲 = − x2d dd2
x + 3x2 ,
2 +x
dx dx
dx

共A2兲

One of us (R.L.D.) acknowledges the support of the Australian Research Council and useful discussions with H.
Friedrich, R. Mennicken, H. Schomerus, G. Spies, J. Wiersig, and N. Witte on spectral and quantum chaos issues forming the background of this paper.

where, as in the main text, dots denote dimensionless derivatives, –␫˙ ⬅ r␮d␫ / dr, etc., and all equilibrium quantities areevaluated at r = r␮.

[1] G. Casati, B. V. Chirikov, and I. Guarneri, Phys. Rev. Lett. 54,
1350 (1985).
[2] M. Wakatani, Stellarator and Heliotron Devices, No. 95 in The
International Series of Monographs on Physics (Oxford University Press, New York, 1998).
[3] A. E. Lifschitz, Magnetohydrodynamics and Spectral Theory(Kluwer, Dordrecht, The Netherlands, 1989).
[4] G. O. Spies and J. A. Tataronis, Phys. Plasmas 10, 413 (2003).
[5] E. Hameiri, Commun. Pure Appl. Math. 38, 43 (1985).
[6] F. Troyon, R. Gruber, H. Saurenmann, S. Semenzato, and S.
Succi, Plasma Phys. Controlled Fusion 26, 209 (1984).
[7] J. R. Ferron et al., Phys. Plasmas 7, 1976 (2000).
[8] D. V. Anderson, W. A. Cooper, R. Gruber, S. Merazzi, and U.
Schwenn, Int. J. Supercomput. Appl. 4, 34 (1990).
[9] C. Schwab, Phys. Fluids B 5, 3195 (1993).
[10] I. B. Bernstein, E. A. Frieman, M. D. Kruskal, and R. M.
Kulsrud, Proc. R. Soc. London, Ser. A 244, 17 (1958).

[11] R. L. Dewar and A. H. Glasser, Phys. Fluids 26, 3038 (1983).
[12] P. Cuthbert, J. L. V. Lewandowski, H. J. Gardner, M. Persson,
D. B. Singleton, R. L. Dewar, N. Nakajima, and W. A. Cooper,
Phys. Plasmas 5, 2921 (1998).
[13] M. H. Redi, J. L. Johnson, S. Klasky, J. Canik, R. L. Dewar,
and W. A. Cooper, Phys. Plasmas 9, 1990 (2002).
[14] M. C. Gutzwiller, Chaos in Classical and Quantum Mechanics, No. 1 in Interdisciplinary Applied Mathematics Series(Springer-Verlag, New York, 1990).
[15] H. J. Stöckmann, Quantum Chaos: An Introduction (Cambridge University Press, Cambridge, 1999).
[16] F. Haake, Quantum Signatures of Chaos, 2nd ed. (SpringerVerlag, Berlin, 2001).
[17] M. L. Mehta, Random Matrices, 2nd ed. (Academic Press, San
Diego, 1991).
[18] R. L. Dewar, P. Cuthbert, and R. Ball, Phys. Rev. Lett. 86,2321 (2001).

066409-10

STATISTICAL CHARACTERIZATION OF THE …

PHYSICAL REVIEW E 70, 066409 (2004)[19] M. V. Berry and M. Tabor, Proc. R. Soc. London, Ser. A 356,
375 (1977).
[20] B. R. Suydam, in Proceedings of the Second International
Conference on the Peaceful Uses of Atomic Energy (United
Nations, Geneva, 1958), Vol. 31, p. 157.
[21] O. K. Cheremhykh and S. M. Revenchuk, Plasma Phys. Controlled Fusion 34, 55 (1992).
[22] H. Sugama and M. Wakatani, J. Phys. Soc. Jpn. 58, 1128
(1989).
[23] R. L. Dewar, C. Nührenberg, and T. Tatsuno, in Proceedings ofthe 13th International Toki Conference, Toki, Japan, 2003 [J.
Plasma Fusion Res. (to be published)].
[24] H. R. Strauss, Plasma Phys. 22, 733 (1980).
[25] R. M. Kulsrud, Phys. Fluids 6, 904 (1963).
[26] T. Tatsuno, M. Wakatani, and K. Ichiguchi, Nucl. Fusion 39,1391 (1999).

[27] B. F. McMillan, R. L. Dewar, and R. G. Storer, Plasma Phys.
Controlled Fusion 46, 1027 (2004).
[28] L. D. Landau and E. M. Lifshitz, Quantum Mechanics (Nonrelativistic Theory), 3rd ed. No. 3 in Course of Theoretical
Physics (Pergamon, Oxford, 1991).
[29] J. W. Connor, R. J. Hastie, and J. B. Taylor, Proc. R. Soc.
London, Ser. A 365, 1 (1979).
[30] R. L. Dewar, M. S. Chance, A. H. Glasser, J. M. Greene, andE. A. Frieman, Tech. Rep. PPPL-1587, Princeton University
Plasma Physics Laboratory (1979), available from National
Technical Information Service, U.S. Department of Commerce, 5285 Port Royal Road, Springfield, VA 22151.
[31] I. Niven, H. S. Zuckerman, and H. L. Montgomery, An Introduction to the Theory of Numbers, 5th ed. (Wiley, New York,
1991).

066409-11