March 26, 1999

Revised October 16, 2003

Fusion plasmas magnetically confined by tokamak devices often exhibit
relaxation oscillations in hydrogen-alpha (or deuterium-alpha) emission
and other quantities in the outer edge of the plasma when in the high-confinement
mode (H-mode) at sufficiently high plasma pressure (beta). These
oscillations are called edge localized modes (ELMs). Although the
ELMs degrade the confinement, they are potentially beneficial for removing
helium ash from a fusion plasma. The ELMs are approximately periodic,
but the period and amplitude generally exhibit variations which are presumably
chaotic. Many theoretical models for the H-mode and ELMs have been
proposed. A particularly simple model was recently advanced by
H. Sugama and W. Horton in Plasma Phys. Control. Fusion **37**, 345-362
(1995). This note shows that their model leads to chaotic dynamics
with some choices of parameters, and it characterizes the chaos thus produced.

The Sugama-Horton model describes the temporal evolution of three variables,
the potential energy *u* associated with the pressure profile, the
turbulent kinetic energy *k*, and the kinetic energy *f* associated
with the background shear flow. In normalized form, these variables
evolve according to the equations:

d*u*/d*t* = *q* - *u*^{1/2}*k*

d*k*/d*t* = *u*^{1/2}*k* - *cu*^{-1/2}*fk*
- *d*^{-1}(*u*)*k*^{2}

d*f*/d*t* = *cu*^{-1/2}*fk* - *cm*(*u*)*f*

where *q* is the normalized potential energy input, which is a
controllable parameter, and *c* is a nondimensional numerical constant.
The function *d*(*u*) is the L-mode anomalous diffusivity, which
depends on *u* through its dependence on background temperature.
The function *m*(*u*) is related to the ion collisional viscosity.
Following Sugama and Horton, we take *d*(*u*) = *u* and
*m*(*u*)
= *u*^{-3/2}(0.95 + 0.05*u*^{5/2}).

Sugama and Horton found limit-cycle solutions for *q* = 1.8 and
*c*
= 5.0 but did not observe chaos. Here we explore the dynamical behavior
of the above system over the whole space of 0 < *q* < 15 and
0 < *c* < 20. The equations were solved using a fourth-order
Runge-Kutta integrator for 32000 steps with a step size of 0.02 and initial
conditions of *u*(0) = *k*(0) = *f*(0) = 0.05. The
character of the solution was identified by calculating the largest
Lyapunov exponent. Exponents greater than 0.008 were assumed
to be chaotic, exponents less than -0.005 were assumed to be fixed points,
and other exponents were assumed to be limit cycles (although they could
also be toruses). A solution was assumed to be unbounded if |*u*|
ever exceeds 5000. The calculation was implemented in PowerBASIC
using the source and (DOS) object
code provided. The resulting plot is shown below:

*q *
*c*

It is likely that some of the unbounded solutions above are a result of initial conditions outside the basin of attraction. It would be even more computationally intensive to test for this possibility.

With a slight modification of the code, the Lyapunov exponent in the chaotic region can be displayed as a gray scale on the same axes as above. The source and object code for accomplishing this are provided. The gray scale is linear in the range of Lyapunov exponents from 0 to 0.2 in the plot shown below:

*q*
*c*

The smallest parameter values that put the dynamics comfortably into
the chaotic region are *q* = 9 and *c* = 9. The remainder
of this note will concentrate on that case. The time evolution of
the three variables were calculated over 32000 time steps starting from
a point on the attractor with a step size of 0.002 using the source
and
object code provided, and the chaotic results
are shown below:

time

In a much longer run, the Lyapunov exponents
(base-*e*) converge to about (0.085, 0, -0.938), giving a Kaplan-York
dimension of about 2.090. An animated view of the attractor is shown
below:

In the image above, the *k* axis points upward, and the object
rotates in the *u*-*f* plane.