J.
6 (2004)
J. Plasma
Plasma Fusion
FusionRes.
Res.SERIES,SERIES,Vol.
Vol.
6 (2004)40–44
000–000

Quantum Chaos Theory and the Spectrum of Ideal-MHD
Instabilities in Toroidal Plasmas
DEWAR Robert L., NÜHRENBERG Carolin1 and TATSUNO Tomoya2
Research School of Physical Sciences and Engineering, The Australian National University, Canberra ACT 0200, Australia
1
Max-Planck-Institut für Plasmaphysik, Teilinstitut Greifswald IPP-EURATOM Association, D-17489, Greifswald, Germany
2
Institute for Research in Electronics and Applied Physics University of Maryland, College Park, MD 20742-3511, USA(Received: 9 December 2003 / Accepted: 30 April 2004)

Abstract
In a fully 3-D system such as a stellarator, the toroidal mode number n ceases to be a good quantum number –all ns within a given mode family being coupled. It is found that the discrete spectrum of unstable ideal MHD(magnetohydrodynamic) instabilities ceases to exist unless MHD is modified (regularized) by introducing a shortperpendicular-wavelength cutoff. Attempts to use ray tracing to estimate the regularized MHD spectrum fail due tothe occurrence of chaotic ray trajectories. In quantum chaos theory, strong chaos in the semiclassical limit leads toeigenvalue statistics the same as those of a suitable ensemble of random matrices. For instance, the probability
distribution function for the separation between neighboring eigenvalues is as derived from random matrix theory andgoes to zero at zero separation. This contrasts with the Poissonian distribution found in separable systems, showingthat a signature of quantum chaos is level repulsion. In order to determine whether eigenvalues of the regularized
MHD problem obey the same statistics as those of the Schrödinger equation in both the separable 1-D case and thechaotic 3-D cases, we have assembled data sets of ideal MHD eigenvalues for a Suydam-unstable cylindrical (1-D)
equilibrium using Mathematica and a Mercier-unstable (3-D) equilibrium using the CAS3D code. In the 1-D case, wefind that the unregularized Suydam-approximation spectrum has an anomalous peak at zero eigenvalue separation. On
the other hand, regularization by restricting the domain of k^ recovers the expected Poissonian distribution. In the 3D case we find strong evidence of level repulsion within mode families, but mixing mode families produces Poissonian
statistics.
Keywords:
quantum chaos, ideal MHD, interchange spectrum, Suydam, finite Larmor radius, eigenvalue spacing,probability distribution

1. Introduction

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 complicated geometries. For this purpose asubstantial investment in effort has been expended on developing numerical matrix eigenvalue programs, such as thethree-dimensional (3-D) TERPSICHORE [4] and CAS3D [5]codes. These solve the MHD wave equations for perturbations
about static equilibria, so that the eigenvalue w 2  –g 2 is real
due to the Hermiticity (self-adjointness [6]) of the linearized
force and kinetic energy operators. They use finite-element
or finite-difference methods to convert the infinite-dimensional Hilbert-space eigenvalue problem to an approximating
finite-dimensional matrix problem.
In order properly to verify the convergence of thesecodes in 3-D geometry it is essential to understand the nature

In ideal MHD the spectrum of the growth rates, g , of
instabilities is difficult to characterize mathematically because
the linearized force operator is not compact [1]. This gives
rise to the possibility of a dense set of accumulation points(descriptively called the “accumulation continuum” by Spies
and Tataronis [2] though more correctly termed [3] theessential spectrum).
The continuous spectrum in quantum mechanics arises
from the unboundedness of configuration space, whereas the
MHD essential spectrum arises from the unboundedness ofFourier space – there is no minimum wavelength in ideal
MHD. This is an unphysical artifact of the ideal MHD model
because, in reality, low-frequency instabilities with | k^| muchgreater than the inverse of the ion Larmor radius, ai, cannot
exist (where k^ is the projection of the local wavevector into
the plane perpendicular to the magnetic field B).

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

Corresponding author’s e-mail: robert.dewar@anu.edu.au
140Dewar R.L. et al., Quantum Chaos Theory and the Spectrum of Ideal-MHD Instabilities in Toroidal Plasmas

cylindrical 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, q andz, so the components of x are completely factorizable into
products of functions of the independent variables separately.
In particular, we write the r-component asof the spectrum – if it is quantum-chaotic then convergenceof individual eigenvalues cannot be expected and a statisticaldescription must be used.
It is the thesis of this paper that the language of quantum
chaos [7] theory indeed provides such a statistical framework
for characterizing MHD spectra in that it seeks to classify
spectra statistically by determining whether, and to what
degree, they belong to various universality classes.
In the cylindrical case the eigenvalue problem is separable into three one-dimensional (1-D) eigenvalue problems,with radial, poloidal, and toroidal (axial) quantum numbers l,m, and n, respectively. It is thus to be expected a priori thatthe spectrum will fall within the standard quantum chaostheory universality class for integrable, non-chaotic systems[7]. In particular, it is to be expected that the probability distribution function for the separation of neighboring eigenvalues
is a Poisson distribution. However, the nature of the MHD
spectrum is quite different from that of the typical quantum,
microwave and acoustic systems normally dealt with inquantum chaos theory and it is necessary to test this conjecture by explicit calculation. In fact we find that the result
depends on the method of regularization.
We first present the eigenvalue equation for a reduced
MHD model of a large-aspect-ratio (effectively cylindrical)
stellarator. We study a plasma in which the Suydam criterion[8] for the stability of interchange modes is violated, so thenumber of unstable modes tends to infinity as the smallwavelength cutoff tends to zero. To compute large-m
eigenvalues we transform to a Schrödinger-like form of theradial eigenvalue equation [9], which has essentially the same
form in configuration (r) space as in Fourier (kr) space, thus
allowing easy regularization by restricting the kr domain. To
simplify even further we approximate the effective potential
by a parabola, thus yielding the quantum harmonic oscillator
equation, solvable in parabolic cylinder functions [10].
Real, finite-aspect-ratio stellarators are fully 3-D andtheir ideal-MHD spectra may be expected a priori to fallwithin the universality class appropriate to time-reversible
quantum chaotic systems, where the spectral statistics arefound to be the same as for a Gaussian orthogonal ensembleof random matrices [7] in regions where ray tracing reveals
chaotic dynamics [11]. At the end of this paper we give a
brief report of 3-D calculations peformed with the CAS3D
code on a Mercier-unstable, high-mirror-ratio, high-iota equilibrium representing a Wendelstein 7-X (W7-X) stellarator
variant [12].

rxr = exp(imq ) exp( - inz )j (r ) ,(1)

where the periodic boundary conditions quantize m and n tointegers and we choose to work with the stream function j (r) rxr (r).
Since the primary motivation of this paper is stellarator
physics, we use the reduced MHD ordering for large-aspect
stellarators [13,14], averaging over helical ripple to reduce toan equivalent cylindrical problem [15,16]. 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.
Defining l  w 2 we seek the spectrum of l-values
satisfying the scalar equation

Lj = lMj(2)

under the boundary conditions j (0) = 0 at the magnetic axis
and j (1) = 0, appropriate to a perfectly conducting wall at
the plasma edge (using units such that r = 1 there).
2
The operator M = –—^ and L is given by(3)
where G is a Suydam stability parameter (> 1/4 for instability[8]), proportional to the pressure gradient p¢(r) and theaverage field line curvature [14].
·In this paper we use the notation f  rf ¢(r) for an
arbitrary function f, so j  rdi/dr is a measure of the magnetic
shear and i measures the variation of the shear with radius.
We observe some differences between eq. (2) and thestandard quantum mechanical eigenvalue problem Hy = Ey.
One is of course the physical interpretation of the eigenvalue– in quantum mechanics the eigenvalue E  h- w is linear in
the frequency because the Schrödinger equation is first orderin time, whereas our eigenvalue l is quadratic in the frequency because it derives from a classical equation of motion.
Another difference is that eq. (2) 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.

2. One-dimensional model eigenvalue
equationIn this paper we study an effectively circular-cylindrical
MHD equilibrium, using cylindrical coordinates such that themagnetic axis coincides with the z-axis, made topologicallytoroidal by periodic boundary conditions. Thus z and thetoroidal angle z are related through z  z/R0, where R0 is the
major radius of the toroidal plasma being modeled by thiscylinder. The poloidal angle q is the usual geometric2
41

Dewar R.L. et al., Quantum Chaos Theory and the Spectrum of Ideal-MHD Instabilities in Toroidal Plasmas

In this paper we work only to lowest order in 1/m, theSuydam approximation. As found by Kulsrud [15], we havethe generalized eigenvalue equation

Equation (2) is very similar to the normal mode equation
analyzed in the early work on the interchange growth rate instellarators by Kulsrud [15]. 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.
Suydam instabilities occur only for values of m and nsuch that n – mk vanishes. For the 1-D numerical work in thispaper we use a parabolic transform profile k = k 0 + k 2r2 as
illustrated in Fig. 1.
Given a rational fraction m = nm/mm in the interval [k(0),k(a)] (where nm and mm are mutually prime) there is a uniqueradius rm such that k (rm) = m. Any pair of integers (m, n)m,n (n mm, n nm), n = 1, 2, 3, … satisfies the resonance condition

nm ,n - mm ,n k (rm ) = 0 .

(0) (0)L j

rm2 Ê ( 0 ) 2 Ë L - l( 0 ) M ( 0 ) ˆ¯ j ( 0 ) = 0 ,m(5)where, more explicitly,(6)with G 2  –l(0)/ j 2 and j and G evaluated at rm. Under theboundary conditions as j (0) Æ 0 as g Æ ±•, eq. (5) can besolved to give a square-integrable eigenfunction, with growthrate g = j G, provided l(0) < 0 is one of the eigenvalues lm ,l.
The radial mode number l = 0, 1, 2, … denotes the numberof nodes of the eigenfunction j (0)  j m ,l (r). Note that lm ,ldepends only on m = n/m and is otherwise independent of themagnitude of m and n.
Restricting attention to unstable modes, so that g 1/2(–l) is real, we transform eq. (5) to the Schrödinger form
[9]

(4)

We use a broad pressure profile that is sufficiently flatnear the magnetic axis that the Suydam instability parameter
G goes to zero at the magnetic axis, and for which p¢ vanishes
at the plasma edge. The resulting G-profile is shown in Fig.
2.
Defining a scaled radial variable x = m(r – rm)/rm, we
can find the large-m spectrum of eq. (2) by expanding all
quantities in inverse powers of m, and equating the LHS tozero order by order.

d 2y+ Q(h)y = 0 ,
dh 2

(7)

Q  G - 14 - 14 sech 2 h - G 2 cosh 2 h ,

(8)

where

Fig. 1

with h defined through x  g sinhh/ j (rm), and y  (coshh)1/2
j (x).
From, e.g., eq. (4.7) of [9] we see that, provided theSuydam criterion G > 1/4 is satisfied, there is an infinity of geigenvalues accumulating exponentially toward the origin
from above (so the l-values accumulate from below) in thelimit 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. (2) at eachlm, l as m Æ • with l fixed. (Although l(0) is infinitely
degenerate, we can break this degeneracy by proceeding
further with the expansion in 1/m, thus showing that lm, l isan accumulation point.) Since the rationals m are dense on
the real line, there is an “accumulation continuum” [2]between g = 0 and the maximum growth rate, g = gmax.

The rotational transform k(r )  1/q (r ) with k0 = 0.45, k2 =0.2. All distinct rational magnetic surfaces m = n/m areshown for m up to 10.

3. Regularization

Fig. 2

The accumulation points of the ideal MHD spectrum
found above are mathematically interesting but exist only asa singular limit of equations containing more physics,
including finite-Larmor-radius (FLR) effects and dissipation,that regularize the spectrum.
In order to proceed further we need to be explicit about
the nature of this singular limit. As we are primarily concerned with the universality class question, we seek only a

The Suydam criterion parameter G (r ) (solid line), andthe instability threshold 1/4 (dashed line), showing
nearly all the plasma is Suydam unstable.

342

Dewar R.L. et al., Quantum Chaos Theory and the Spectrum of Ideal-MHD Instabilities in Toroidal Plasmas

where n = l in the unregularized case, kmax = hmax = •. In theeven-l, regularized case n may be found by solving for a zeroof M (n/2, 1/2, (4G 2 – 1)1/2 h2max / 2), where M is Kummer’sfunction. For l = 0, n becomes exponentially small as hmax Æ•, which allows an approximate regularization formula to be
derived.
We study the spectrum between the maximum l = 1
growth rate, gmax (l = 1)  maxm gm,1, and the maximum overall
growth rate, gmax = maxm gm,0. Only the l = 0 modes exist in
this range of g , which corresponds to the range in m between
mmin ª 0.522 and mmax ª 0.628. Throughout this range G is >1/2, so that Q has a single minimum [9] and the quadratic
approximation of this section is appropriate. In this range
there are only four low-order rationals n/m with m < 10.
Taking r* = 0.001, all pairs of integer values m, n in the
fan-shaped region 1 £ m £ mmax (n/m), mmin £ n/m £ mmax were
evaluated, giving an initial dataset of over 32,000 points (m,n). The corresponding set of unregularized eigenvalues wascalculated by solving eq. (13) with n = 0 and the eigenvalues
were sorted and numbered from the top to give the integrated
density of states “staircase” function N (g ).
The curve 0.3523 – 9.5733 ¥ 10–11 N 2 – 1.1625 ¥ 10–20 N 4
was found to give a good fit to the smoothed behavior of thisfunction. Inverting this function gives the smoothed function
Ñ(g ) which is used to “unfold” [7] spectra by defining a new“energy eigenvalue” E  Ñ (g ), such that N (E) increases
linearly on average.
This means that the average separation of eigenvalues is
now unity, making comparison with spectra from otherphysical systems meaningful and allowing universal behavior
to become apparent if present. However, Fig. 3 shows thatthe probability distribution of eigenvalue spacings s is far
from universal for the unregularized Suydam spectrum,
exhibiting a delta-function-like spike at s = 0. This ispresumably because, although we have truncated the spectrum
in m, we have not removed the degeneracies arising for loworder rationals m in the range mmin < m < mmax.
Fig. 4 on the other hand shows that when a similar
procedure is applied to the regularized spectrum (retaining
only regularized eigenvalues above gmax (l = 1), the universal

minimal modification of eq. (2) that has some physical basis
but makes as little change to ideal MHD as possible. To preserve the Hermitian nature of ideal MHD we cannot use thedrift correction used for estimating FLR stabilization ofinterchange modes by Kulsrud [15]. However it is possible
to effect a pseudo-FLR regularization of ideal MHD by
restricting k^ to a disk of radius less than the inverse ion
Larmor radius. In our nondimensionalized, large-aspect ratio
model this implies( kq2 + kr2 )1/ 2 r* £ 1 ,(9)where kr and kq are the radial and poloidal components of the
wavevector, respectively, r* and is the ion Larmor radius (ata typical energy) in units of the minor radius.
To apply eq. (9) precisely we need to relate kr and kq tothe eigenvalue problem discussed above. From eq. (1) we seethat kq = m/r. We define kr as the Fourier variable conjugate
to r. Fourier transformation of eq. (2) is only practical in thelarge-m limit, when modes are localized near the resonantsurfaces r = rm, which is why we have restricted the discussion
to leading order in the 1/m expansion.
Using the stretched radial coordinate x  m(r – rm)/rm wedefine kr  mk/rm, where k is the Fourier-space independentik,
variable conjugate to x. With the substitutions d/dx
id/dk, and using the fact that kd/dk and (d/dk)k  1 +x
kd/dk commute, eq. (5) transforms toÈ- d (1 + k 2 ) d + G 2 (1 + k 2 ) - G  j = 0 .
ÍÎ d k
 k
dk

(10)

The transformation k = sinhh then leads back to Eq. (7), withl now to be interpreted as a distorted Fourier-space
independent variable, rather than as a real-space coordinate !
Equation (9) implies that Eq. (10) is to be solved on thedomain –kmax £ k £ kmax whereÈÊ rm ˆ 2 
k max ( m )  ÍÁ˜ - 1
ÍÎË mr* ¯


1/ 2

.

(11)

This exists provided |m| < mmax, where

mmax ( m )  rm / r* .

(12)

Analogously to quantum mechanical box-quantization we use
Dirichlet boundary conditions at ±kmax.

4. Spectral statistics in the 1-D case
As only the qualitative nature of the spectrum is important, we approximate the function Q by Q(0) + –12 Q (0) h2,
so eq. (7) can be solved in parabolic cylinder functions [10].
We find the dispersion relation
Fig. 32

n + 12 =- 12
1/ 2 ,

G-G( 4G 2 - 1)

(13)

4
43

The histogram shows an estimate, based on a data setof about 32,000 unregularized eigenvalues, of theprobability distribution function for the eigenvalue
separation s. The plot is dominated by the spike at s =0.

Dewar R.L. et al., Quantum Chaos Theory and the Spectrum of Ideal-MHD Instabilities in Toroidal Plasmas

statistics are analyzed within the three mode families the
eigenvalue spacing distribution function is closer to the
Wigner conjecture form found for generic chaotic systems [7]than to the Poisson distribution for separable systems, as
might be expected from [11]. However, when the spectra fromthe three uncoupled mode familes are combined, there areenough accidental degeneracies that the spacing distribution
becomes close to Poissonian.

Poisson distribution expected from a separable system isobtained to a good approximation, thus leading to the expectation that generic quantum chaos theory is applicable once
any physically reasonable regularization is performed.
Further support for this hypothesis is obtained from a
CAS3D study of a W7-X variant equilibrium with a nonmonotonic, low-shear transform profile (kaxis = 1.1066, kmin =1.0491, k edge = 1.0754). As seen from Fig. 5, when the

Acknowledgments:
Part of this work was performed at the University ofTokyo, Graduate School of Frontier Sciences (RLD, TT). We
thank Professor Zensho Yoshida for his hospitality, supportand useful discussions. RLD was partially supported by theAustralian Research Council.

References

Fig. 4

Eigenvalue spacing distribution for a data set of about
12,000 regularized eigenvalues. The exponential curveshows the Poisson distribution.

Fig. 5

Eigenvalue spacing distributions from mode family
datasets N = 0 (137 values), N = 1 (214 values) and N =2 (178 values) from a W7-X-like equilibrium, and the
distribution for the combined spectrum, N = 0, 1 and 2.

[1] A.E. Lifschitz, Magnetohydrodynamics and Spectral
Theory (Kluwer, Dordrecht, The Netherlands, 1989), p.
416.
[2] G.O. Spies et al., Phys. Plasmas 10, 413 (2003).
[3] E. Hameiri, Commun. Pure Appl. Math. 38, 43 (1985).
[4] D.V. Anderson et al., Int. J. Supercomp. Appl. 4, 34(1990).
[5] C. Schwab, Phys. Fluids B5, 3195 (1993).
[6] I.B. Bernstein et al., Proc. R. Soc. London Ser. A 244,
17 (1958).
[7] F. Haake, Quantum Signatures of Chaos, 2nd ed.
(Springer-Verlag, Berlin, 2001).
[8] B.R. Suydam, Proc. Second Int. Conf. on the PeacefulUses of Atomic Energy (United Nations, Geneva, 1958),
Vol. 31, p. 157.
[9] O.K. Cheremhykh et al., Plasma Phys. Control. Fusion
34, 55 (1992).
[10] Handbook of Mathematical Functions, edited by M.
Abramowitz et al. (Dover, New York, 1965).
[11] R.L. Dewar et al., Phys. Rev. Lett. 86, 2321 (2001).
[12] C. Nührenberg, Phys. Plasmas 3, 2401 (1996).
[13] H.R. Strauss, Plasma Phys. 22, 733 (1980).
[14] M. Wakatani, Stellarator and Heliotron Devices, No. 95
in The international series of monographs on physics(Oxford University Press, New York, 1998).
[15] R.M. Kulsrud, Phys. Fluids 6, 904 (1963).
[16] T. Tatsuno et al., Nucl. Fusion 39, 1391 (1999).

544