Journal of the Korean Physical Society, Vol. 50, No. 1, January 2007, pp. 112∼117

Review Articles

Quantum Chaos? Genericity and Nongenericity in the MHD
Spectrum of Nonaxisymmetric Toroidal Plasmas
R. L. Dewar∗ and B. G. Kenny
Department of Theoretical Physics, Research School of Physical Sciences and Engineering,
The Australian National University, ACT 0200, Australia

C. Nührenberg
Max-Planck-Institut für Plasma Physik, Teilinstitut Greifswald, D-17491 Germany

T. Tatsuno
Center for Scientific Computation and Mathematical Modeling,
University of Maryland, College Park, MD 20742-3289, U.S.A.

B. F. McMillan
Centre de Recherche en Physique des Plasmas, Ecole Polytechnique Fédérale de Lausanne, CH-1015 Lausanne, Switzerland(Received 15 August 2006)
The eigenmode spectrum is a fundamental starting point for the analysis of plasma stability andthe onset of turbulence, but the characterization of the spectrum even for the simplest plasma
model, ideal magnetohydrodynamics (MHD), is not fully understood. This is especially true inconfigurations with no continuous geometric symmetry, such as in a real tokamak when the discrete
nature of the external magnetic field coils is taken into account, or the alternative fusion concept,
the stellarator, where axisymmetry is deliberately broken to provide a nonzero winding number(rotational transform) on each invariant torus of the magnetic field line dynamics (assumed forpresent purposes to be an integrable Hamiltonian system). Quantum (wave) chaos theory provides
tools for characterizing the spectrum statistically, from the regular spectrum of the separable case(integrable semiclassical dynamics) to that where the semiclassical ray dynamics is so chaotic thatno simple classification of the individual eigenvalues is possible (quantum chaos). The MHD spectrum exhibits certain nongeneric properties, which we show, using a toy model, to be understablefrom the number-theoretic properties of the asymptotic spectrum in the limit of large toroidal andpoloidal mode (quantum) numbers when only a single radial mode number is retained. Much more
realistically, using the ideal MHD code CAS3D, we have constructed a data set of several hundred
growth-rate eigenvalues for an interchange-unstable three-dimensional stellarator equilibrium with arather flat, nonmonotonic rotational transform profile. A statistical analysis of eigenvalue spacingsshows evidence of generic quantum chaos, which we attribute to the mixing effect of having a large
number of radial mode numbers.
PACS numbers: 52.35.Bj, 05.45.Mt
Keywords: Fusion plasma, Stellarator, Interchange instability, Suydam, Mercier, Essential spectrum, Quantum chaos, Farey sequence

the other hand, is inherently nonaxisymmetric, and thelack of a continuous symmetry means there are no “goodquantum numbers” to characterize the spectrum.
This makes the numerical computation of the spectrum a challenging task, but numerical matrix eigenvalue programs, such as the three-dimensional TERPSICHORE [1] and CAS3D [2] codes, are routinely used
to assess the ideal magnetohydrodynamic (MHD) stability of proposed fusion-relevant experiments. An example
is the design of the 5-fold-symmetric Wendelstein 7-X

I. INTRODUCTION
The tokamak and the stellarator fusion concepts bothseek to contain a plasma in a toroidal magnetic field. To
a good approximation, the tokamak field is axisymmetric, so the study of the spectrum of normal modes ofsmall oscillations about equilibrium is simplified by theexistence of an ignorable coordinate. A stellarator, on∗ E-mail: robert.dewar@anu.edu.au-112-

Quantum Chaos? Genericity and Nongenericity in the MHD· · · – R. L. Dewar et al.

-113-

Fig. 1. Left: plasma boundary of a 5-periodic toroidal equilibrium geometrically related to W7-X configurations (at hβi =0.05, VMEC calculation). Right three frames: vertical cross-sections of the configuration at three different toroidal angles withina field period, with Z the vertical coordinate and R the distance from the Z axis. The plasma boundary has been scaled suchthat the minor radius of the torus a ≈ 1 and the major radius equals the aspect ratio R/a = A.

(W7-X) stellarator, where CAS3D was used [3] to study
a number of different cases.
The configurations foreseen to be studied experimentally in W7-X are MHD stable, but in the present paper
we are concerned with an unstable case from this study, ahigh-mirror-ratio, high-rotational-transform equilibrium(Fig. 1). Due to its less pronounced shaping, this caseis quite unstable, which contrasts with the properties ofgenuine W7-X configurations. The three-dimensional nature of the equilibrium breaks all continuous symmetries,
coupling both poloidal (m) and toroidal (n) Fourier harmonics and thus precluding separation of variables andsimple classification of the eigenvalues.
These eigenvalues, ω 2 ≡ −γ 2 , are real due to theself-adjointness [4] of the force and kinetic energy operators in ideal MHD, linearized about a static equilibrium. This is analogous to the Hermitian nature ofquantum mechanics, so we might a priori expect to beable to take over the mathematical techniques used inquantum mechanics. Thus, we study the W7-X Mercier(interchange)-unstable case mentioned above by using
statistical techniques from the theory of quantum chaos[5, eg].
This is of practical importance for numerical analysisof the convergence of eigenvalue codes because, if the system is quantum-chaotic, convergence of individual eigenvalues cannot be expected, and a statistical description
must be used. However, there is a fundamental question
as to whether the ideal MHD spectrum lies in the same
universality classes as the spectra of typical quantum
mechanical systems.
This question has been addressed recently [6, 7] bystudying the interchange unstable spectrum in an effectively cylindrical model of a stellarator. In the cylindricalcase the eigenvalue problem is separable into three onedimensional eigenvalue problems, with radial, poloidal,
and toroidal (axial) quantum numbers l, m, and n, respectively. If the spectrum falls within the generic quantum chaos theory universality class for integrable, nonchaotic systems [8], then the probability distributionfunction for the separation of neighboring eigenvalues isa Poisson distribution. However, the work of Ref. 6 indicates that the universality class depends on the methodof regularization (ie truncation of the countably infinite
set of ideal-MHD interchange growth rates): a smooth,physically-motivated finite-Larmor-radius roll-off in the
spectrum appears to give the generic Poisson statistics
for separable systems, but a sharp truncation in m andn gives highly non-generic statistics. The latter case isless physical, but corresponds closely to the practice in
MHD eigenvalue studies of using a restricted m, n basis
set, but a relatively fine mesh in the radial direction.
A careful analysis of the spectrum of ideal-MHD
interchange modes in a separable cylindrical approximation [7] revealed non-generic behavior of the spectral statistics-a bimodal probability distribution function (PDF), rather than the expected Poisson distribution. The non-genericity of this separable case indicatesthat caution must be applied in applying conventional
quantum chaos theory in non-separable geometries. The
study [7] indicated that the non-generic behavior of idealMHD interchange modes is due to the peculiar featureof the dispersion relation for these modes that the eigenvalues in the short-wavelength limit depend only on the
direction of the wave vector, not on its magnitude. (Thisis unusual behavior, but it is shared with internal gravity
waves in geophysical fluid dynamics.) It was suggested
in Ref. 7 that the detailed features of the spectrum could
be understood from the properties of Farey sequences.
In the present paper, we discuss fictitious eigenvalues
from a toy model, used to elucidate the importance of
number-theoretic effects, that illustrate how nongeneric
the MHD eigenvalue spectrum can be if only one radial
eigenmode is used. Then we present the results of aquantum chaos analysis for the W7-X case.

II. TOY EIGENVALUE PROBLEM
To gain insight into the nongeneric behavior found in
the separable case [7] with only one radial eigenmode, we
study the energy spectrum {En,m } for a “toy” quantum

-114-

Journal of the Korean Physical Society, Vol. 50, No. 1, January 2007

The list G(mmax ) ≡ {En,m }, sorted into a non-decreasing
sequence {Ei |i = 1, 2, . . . , Nmax }, is very similar to theFig. 2. Number N (E) of renormalized eigenvalues Em,nbelow a given value E for mmax = 10: (a) for the Hamiltonian
H = pφ /pθ ; (b) for the Farey sequence.

mechanical Hamiltonian H = pφ /pθ , where the configuration space is the 2-torus θ ∈ [0, 2π), φ ∈ [0, 2π) withperiodic boundary conditions. In the semiclassical approximation, we see that H depends only on the direction of p, not its magnitude, as for MHD interchange
modes and internal gravity waves.
The eigenvalue problem is the time-independent
Schrödinger equation, Hψ = Eψ, the eigenfunctions being exp[i(mθ + nφ)]/4π 2 , where m and n are integers.
The eigenvalues are given by En,m = n/m (m 6= 0).
Note the singular nature of the spectrum-it is discrete,
yet infinitely dense, the rationals being dense on the realline. Also, the spectrum is infinitely degenerate as eigenvalues are repeated whenever m and n have a common
factor. Mathematically, such a spectrum, neither point
nor continuous, belongs to the essential spectrum [9].
In order to analyze this spectrum by using standard quantum chaos techniques, we first regularize itby bounding the region of the m, n lattice studied, andthen allowing the bound to increase indefinitely. Fortunately, the PDF, P (s), is independent of the precise
shape of the bounding line when we follow standard practice [5] in renormalizing (unfolding) the energy levels to
make the average spacing unity. Thus, we adopt thesimplest choice, taking the bounded region to be the triangle 0 ≤ n ≤ m ≤ mmax . As the points (n, m) forma lattice in the plane with mean areal density of 1, we
can estimate the asymptotic, large-mmax behavior of thenumber of levels, Nmax , from the area of the bounding
triangle: the m axis, the line n = m and the line m = 1,
which gives the “Weyl formula” [5] Nmax ∼ m2max /2.
The Farey sequence F(Q) is the set of all rational numbers p/q between 0 and 1, 1 < q ≤ Q, arranged in orderof increasing arithmetic size and keeping only mutuallyprime pairs of integers p and q. Farey sequences are important in number theory [10] and have application invarious dynamical systems problems, such as the theory of mode locking in circle maps [11]. They even havea connection with the famous Riemann hypothesis [12].
They are conveniently generated using the Mathematica
algorithm [13]Farey[Q_] := Union[{0,1},
Flatten[Table[p/q,{q, Q},{p, q-1}]]].

Fig. 3. Separation statistics: (a) for the model Hamiltonian; (b) for the Farey sequence. The solid curves are fromthe Farey spacing measure, Eq. (2). In (b) the short-dashedcurve is for the Poisson Process of the generic integrable
problem and the long-dashed curve is that for the Gaussianorthogonal ensemble of random matrices (quantum chaoticcase).

Farey sequence F(mmax ) except for the high degeneracy(multiplicity) of numerically identical levels, especiallywhen n/m is a low-order rational.
Define the renormalized (unfolded) energy as En,m ≡Nmax En,m . The normalization by Nmmax ensures that
ENmax = Nmax , so the mean slope of the Devil’s staircase shown in Fig. 2(a) is unity. The large vertical steps
visible in Fig. 2(a) are due to the high degeneracy at
low-order rational values of n/m. This high degeneracy
is also the cause of the delta function spike at the origin visible in the level-separation probability distribution
plot shown in Fig. 3(a), which is very similar to Fig. 9(a)in [7].
The tail in Fig. 3(a) is due to the non-degenerate component of the spectrum, obtained by reducing all fractions p/q to lowest terms and deleting duplications. Thus
the eigenvalues in this set are NQ times the terms of the
Farey sequence F(Q).
To study the statistics of this non-degenerate component, one naturally defines the Farey spectrum {EiF } as
N F (Q) times the terms of the Farey sequence F(Q),where N F (Q) is the number of terms in F(Q). The
asymptotic behavior of N F (Q) in the large-n limit isgiven [14, p. 391] by N F (Q) ∼ 3Q2 /π 2 +O(Q ln Q). The
staircase plot and the separation distribution P F (s) for
the Farey spectrum are given in Fig. 2(b) and Fig. 3(b),respectively.
It is a standard result in the theory of Farey sequences[10, p.301] that the smallest and the largest nearestneighbor spacings in F(Q) are given, respectively, by
1
Q(Q − 1)

and

1
,
Q

(1)

so that the support of the tail component of P (s) in
Fig. 3(a) becomes [1/2, ∞) in the limit Q → ∞, whilethat of P F (s) in Fig. 3(b) is [3/π 2 , ∞). Augustin et al.
[12], Eq. (1.9), derived the spacing density for the Farey
sequence asQuantum Chaos? Genericity and Nongenericity in the MHD· · · – R. L. Dewar et al.

g1 (t) ≡


0,


 6

π 2 t2 ln





12
π 2 t2 ln

 2 
π t

h 23 
π t
6

,
1−

q

1 − π122 t

i

for
for

0 ≤ t ≤ π32 ,312
π2 ≤ t ≤ π2 ,for12
π 2 ≤ t.

-115-

(2)

Fig. 4. Left: the rotational transform ι- = ι-(s) profile versus normalized toroidal flux s (s ∝ r2 near the magnetic axis).
Right: a measure of the Mercier stability versus normalized
toroidal flux for the 3-dimensional W7-X-like case of Fig. 1.
A negative value indicates instability.

The solid curve in Fig. 3(b) is obtained by setting
P F (s) = g1 (s) and is seen to agree well with the numerical results. The solid curve in Fig. 3(a) is obtained by
setting P (s) = [N F (Q)/N (Q)]2 g1 (N F s/N ) and agrees
well with the tail of the histogram. The ratio of thearea of the tail in Fig. 3(a) to the strength of the delta
function in Fig. 3(a) is N F (Q)/[N (Q) − N F (Q)] ≈ 1.55.
We have verified that the probability distributions remain unchanged if subranges of the spectra are used, inagreement with the result included in Theorem 1.1 ofRef. 12 that the convergence to a probability measure isindependent of the interval chosen.

Fig. 5. Choice of basis set of toroidal and poloidal Fourier
harmonics.

III. W7-X RESULTS
The W7-X variant equilibrium studied was generated
with the VMEC [15] code, which assumes the magnetic
field to be integrable, so that all magnetic field lines lie on
nested toroidal flux surfaces, which we label by s, the enclosed toroidal magnetic flux divided by the toroidal flux
enclosed by the plasma boundary. The magnetic field ischaracterized on each flux surface by its winding number
ι-(s). (In tokamaks its inverse, q ≡ 1/ι-, is more commonly
used.) As seen in Fig. 4, the rotational transform profile is nonmonotonic and has low shear (ι-axis = 1.1066,
ι-min = 1.0491, ι-edge = 1.0754), so it is close to, but
greater than, unity over the whole plasma. Figure 4 also
shows the equilibrium to be interchange unstable, because the Mercier stability criterion is violated over thewhole plasma.
The CAS3D code expands the eigenfunctions in a finite
Fourier basis in the toroidal and the poloidal angles, selected so as to include all |n| ≤ nmax and all m such thatn/m lies in a band including the range of ι-. The Fourier
tableau is depicted graphically in Fig. 5. In this code,

Fig. 6. Unfolded eigenvalue spacing distributions frommode-family datasets N = 0 (137 values), N = 1 (214 values) and N = 2 (178 values) calculated by CAS3D for ourW7-X-like equilibrium, and the distribution for the combined
spectrum, N = 0, 1 and 2.

the radial dependence of the perturbation functions istreated by using a hybrid Finite-Element approach, witha linear interpolation for the normal displacement andpiecewise constant interpolations for the scalar components that describe the MHD displacement within themagnetic surfaces. In the calculations discussed here,
301 radial grid points were used. The kinetic energy wasused as the normalization; therefore, the unstable eigen-

-116-

Journal of the Korean Physical Society, Vol. 50, No. 1, January 2007

Fig. 7. Dyson-Mehta spectral rigidity as a function ofsubinterval length L. N = 0 family: black short-dashed;
N = 1: gray long-dashed; N = 2: gray solid; combined data
set: gray short-dashed; Poisson-process: black long-dashed.

form found for generic chaotic systems [5] than to the
Poisson distribution for separable systems, as might be
expected from Ref. 18. However, when the spectra fromthe three uncoupled mode-families are combined, there
are enough accidental degeneracies for the spacing distribution to become close to Poissonian.
To test the sensitivity to the precise method of unfolding chosen, we did the statistics using two different
methods: the Gaussian unfolding [17] and a fit with exponentials. The results, shown in Fig. 8, indicate little
sensitivity to the unfolding method.

IV. CONCLUSION

Fig. 8. Unfolded eigenvalue spacing distributions formode-family N = 1, calculated using two different unfolding methods. The results are seen to be consistent to within
statistical error.

values λ may be converted
p to a nondimensional growthrate γ viz. γτA = R0 (0) |λ|/B0 (0). Here, R0 (0) is the
major radius, and B0 (0) the equilibrium magnetic fieldmeasured on the magnetic axis.
Because of the 5-fold symmetry of the equilibrium, any
toroidal Fourier harmonic n in an eigenfunction is coupled to toroidal harmonics n ± 5. With the poloidal harmonics chosen to be positive, m ≥ 0, there are just three
uncoupled mode-families N = 0, 1, 2 (compare [2]).
We characterize the statistics of the ensembles of
eigenvalues within the three mode-families using two
standard measures from quantum chaos theory [5, 16,17], first renormalizing (“unfolding”) the eigenvalues sotheir average separation is unity. The first measure,shown in Fig. 6, is the probability distribution function,
P (x), for the eigenvalue separation x (the s of Sec. II).
The other, shown in Fig. 7, is the Dyson-Mehta rigidity∆3 (L), where L is the subrange of unfolded eigenvalues
used.
As seen from Fig. 6, when the statistics are analyzed
within the three mode-families the eigenvalue spacing
distribution function is closer to the Wigner conjecture

Although not presented here, when all unstable eigenmodes (i.e., all l, m, and n) are included, the eigenvalue
spacing statistics for the ideal-MHD interchange eigenvalue spectrum in a separable cylindrical approximation
is close to that of generic separable wave equations [8],despite our earlier finding [6,7] that the spectrum in thesubspace of the most unstable radial eigenmode l = 0
is nongeneric, as explained by the model presented inSec. II. . In this paper we have shown that, in a strongly
three-dimensional stellarator equilibrium related to W7X, the unstable interchange (Mercier) mode spectrum
has, to within statistical uncertainties, similar statistics
to that of generic quantum chaotic systems. That is,
the overwhelming majority of eigenvalues are not “goodquantum numbers” and can, thus, be expected to display sensitivity to small perturbations. This needs tobe borne in mind when doing convergence studies using
stability codes such as CAS3D.
An interesting question for further work is whether
other modes, such as drift waves, are quantum chaotic
in stellarators, or if this is a peculiarity of MHD modes.
There is already evidence that kinetic effects make thesemiclassical (WKB) dynamics closer to integrable [19].
Another question is whether quantum chaos of Mercier
modes occurs in machines with more field periods than
the 5 of W7-X. Earlier work [20] suggested that, in a 10field-period heliotron equilibrium related to the Large
Helical Device (LHD), the spectrum is close to that ofan equivalent axisymmetric torus, and thus not chaotic.
However, a relatively few modes were studied and thespacing statistics were not calculated.

ACKNOWLEDGMENTS
One of us (RLD) acknowledges support by the Australian Research Council and thanks the organizers of
the Fourth Dynamics Days Asia Pacific (DDAP4) conference at which this work was presented.

Quantum Chaos? Genericity and Nongenericity in the MHD· · · – R. L. Dewar et al.

REFERENCES[1] D. V. Anderson, W. A. Cooper, R. Gruber, S. Merazzi
and U. Schwenn, Int. J. Supercomp. Appl. 4, 34 (1990).
[2] C. Schwab, Phys. Fluids B 5, 3195 (1993).
[3] C. Nührenberg, Phys. Plasmas 3, 2401 (1996).
[4] I. B. Bernstein, E. A. Frieman, M. D. Kruskal and R. M.
Kulsrud, Proc. R. Soc. London Ser. A 244, 17 (1958).
[5] F. Haake, Quantum Signatures of Chaos, 2nd ed.
(Springer-Verlag, Berlin, 2001).
[6] R. L. Dewar, C. Nührenberg and T. Tatsuno, J.
Plasma Fusion Res. SERIES 6, 40 (2004), Proceedingsof the 13th International Toki Conference (Toki, 2003,
arXiv:physics/0409070.
[7] R. L. Dewar, T. Tatsuno, Z. Yoshida, C. Nührenberg
and B. F. McMillan, Phys. Rev. E 70, 066409 (2004),
arXiv:physics/0405095.
[8] M. V. Berry and M. Tabor, Proc. R. Soc. Lond. A 356,
375 (1977).
[9] E. Hameiri, Commun. Pure Appl. Math. 38, 43 (1985).
[10] I. Niven, H. S. Zuckerman and H. L. Montgomery, An
Introduction to the Theory of Numbers, 5th ed. (Wiley,
New York, 1991).
[11] R. Artuso, P. Cvitanović and B. G. Kenny, Phys. Rev.

-117-

A 39, 268 (1989).
[12] V. Augustin, F. P. Boca, C. Cobeli and A. Zaharescu,
Math. Proc. Camb. Phil. Soc. 131, 23 (2001).
[13] P. Abbott, The Mathematica Journal 2 Tricks of
the Trade Section (1992), http://www.mathematicajournal.com/issue/v2i2/.
[14] P. Cvitanović, R. Artuso, R. Mainieri, G. Tanner
and G. Vattay, Classical and Quantum Chaos, 10thed. (ChaosBook.org Niels Bohr Institute, Copenhagen,2003), Webbook: http://chaosbook.org/.
[15] S. P. Hirshman and O. Betancourt, J. Comput. Phys.
96, 99 (1991).
[16] M. L. Mehta, Random Matrices, 2nd ed. (Academic
Press, San Diego, 1991).
[17] H. Bruus and J.-C. Anglès d’Auriac, Phys. Rev. B 55,
9142 (1997).
[18] R. L. Dewar, P. Cuthbert and R. Ball, Phys. Rev. Letters86, 2321 (2001), e-Print arXiv:physics/0102065.
[19] B. F. McMillan and R. L. Dewar, Nucl. Fusion 46, 477(2006).
[20] W. A. Cooper, D. B. Singleton and R. L. Dewar, Phys.
Plasmas 3, 275 (1996), Erratum: Phys. Plasmas 3, 3520
(1996).