#  The Relationship between Flux Coordinates and Equilibrium-based Frames of
Reference in Fusion Theory

S.E. Kruger  Tech-X Corporation<EMAIL_ADDRESS>J.M. Greene  Deceased. The
original idea of the $JK$ reference frame and the development of the new
annihilation operator is due strictly to Dr. J.M. Greene.

(September 2, 2024)

###### Abstract

The properties of two local reference frames based on the magnetic field and
the current density are investigated for magnetized plasmas in toroidal
geometry with symmetric angle. The magnetic field-based local frame of
reference has been well-studied for example by Dewar and colleagues [Phys.
Fluids 27, 1723 (1984)] An analogous frame based on the current density vector
is possible because it is also divergence free and perpendicular to the
gradient of the poloidal flux. The concept of straightness of a vector is
introduced and used to elucidate the Boozer and Hamada coordinate systems. The
relationship of these local frames to the more well-known Frenet frame of
reference, which specifies a curve in terms of curvature and torsion, is
given. As an example of the usefulness of the these formal relationships, we
briefly review the ideal MHD theory and their use. We also present a new
annihilation operator, useful for eliminating shorter time scales than the
time scale of interest, for deriving the inner layer equations of Glasser,
Greene, and Johnson [Phys. Fluids 18, 875 (1975)]. Compared to the original
derivation that is based on the local frame of reference in terms of the
magnetic field, the new annihilation operator that is based on the local frame
of reference in terms of the current density simplifies the derivation.

###### pacs:

52.30.-q, 52.65.-y, 52.35.-g, 52.55.Tn, 52.40.Hf, 02.60.Lj, 52.35.Vd

††preprint: DRAFT

##  I Introduction

High density, magnetized plasmas, such as tokamak or stellarator plasmas, are
characterized by a large number of temporal and spatial scales. Laboratory
magnetized plasmas have had great success in creating plasmas with confinement
times much longer than the shortest time scales. Understanding plasma behavior
on these longer times requires analytic, and computational, techniques for
stepping over the shortest time scales to study the time scales of interest. A
key issue for magnetized plasmas is resolving the anisotropy of the magnetic
field. For example, examination of the Braginskii transport equations reveals
different transport time scales for the directions parallel to the magnetic
field, perpendicular to the magnetic flux surfaces, and the binormal direction
that is perpendicular to both directions Braginskii (1966). This anisotropy
arises consistently in fluid and kinetic theories.

A common method for analytically handling this anisotropy is to use a local
frame of reference. By frame of reference, we mean the development of a local
coordinate system at a point in space. For example, the three directions above
enable the development of a set of basis vectors for a local coordinate
system. This local frame of reference has long been used in fluid theory
Greene and Johnson (1962); Coppi et al. (1966) in combination with the use of
a flux coordinates Greene and Johnson (1962) appropriate for toroidal systems.
Dewar and colleagues Dewar et al. (1984), hereafter referred to as DMS,
greatly improved the understanding of this local reference frame by
formalizing many of the relationships between flux coordinates and the local
frame that are commonly used in analytic derivation. This frame of reference
is also implicit in much of the work of Stix Stix (1992) with the “Stix frame”
Wright (1998) being a common tool in RF theory. That frame, like the more
well-known Frenet frame for arbitrary parameterized curves, is expressed in
terms of dimensionless vectors. As we shall see, for most of the theory
calculations, using dimensional local basis vectors is preferred.

Flux coordinate systems make use of two properties of the equilibrium magnetic
field. The first is that the magnetic field is divergence-free. The second is
the equilibrium relationship:

$\mathbold{J}\times\mathbold{B}=\mathbold{\nabla}P.$ (1)

which is the zeroth-order force balance in a magnetized plasma Hazeltine and
Meiss (1992). Dotting both sides with $\mathbold{B}$ shows that
$\mathbold{B}\cdot\mathbold{\nabla}P=0$. This combined with symmetry allows
the definition of a flux function such that
$\mathbold{B}\cdot\mathbold{\nabla}\psi=0$. Quasineutrality enables the charge
continuity equation to be written as $\mathbold{\nabla}\cdot\mathbold{J}=0$.
The equilibrium relationship also shows that
$\mathbold{J}\cdot\mathbold{\nabla}P=0$ or that
$\mathbold{J}\cdot\mathbold{\nabla}\psi=0$ as well. Thus, it is obvious that
another local frame of reference, completely analogous to the magnetic field-
based frame is possible using the current density vector. To our knowledge,
this has never been explored before. In this paper, we formalize the
relationships between the two reference frames and flux coordinate systems
similar in spirit to the work of the DMS paper. We also relate this frame of
reference to the Frenet frame from differential geometry of curves, which
expresses a curve in terms of curvature and torsion. Fusion theory has favored
magnetic shear over torsion, and these relationships are reviewed.

As an example of the usefulness of formalizing the relationships between flux
coordinates and local frames of reference, their usefulness in ideal MHD and
resistive MHD instability theory is reviewed. Resistive instabilities, which
have a time scale that is a hybrid of the Alfvèn and resistive diffusion time
scales, is much slower than the Alfvèn time scales. Analytically, this
requires going to second order in the ordering, which involves significant
algebra, especially in studying instabilities in toroidal geometry with
complicated flux surface shapes. The rest of the paper reviews the derivation
of the inner layer equations in toroidal geometry. The first discussion of how
to derive these equations is contained in the appendix of Johnson and Greene
Johnson and Greene (1967) in 1967, but they were not fully derived and
published until Glasser, Greene, and Johnson (GGJ) in 1967. Unfortunately,
that paper only focused on the analysis of these equations and only presented
them. Part of the point of this paper is that the derivation of these
equations contains pedagogical value, and a new method of deriving these
equations is presented.

The rest of the paper proceeds as follows. First, we introduce the concept of
“equilibrium-based frame of references” and introduce our new frame. Next,
following the work of DMS, we provide the mapping of the equilibrium-based
frames of references to the flux coordinate systems. We then briefly discuss
the relationship of equilibrium reference frames to the Frenet frame and
common geometric quantities. With the mathematical discussion complete, we
then discuss the ideal MHD formulas, and derivation of the inner layer
equations in toroidal geometry. Key to the derivation is the use of
annihilation operators Coppi et al. (1966); Glasser et al. (1975). For the
derivation of the inner layer equations in toroidal geometry, we find that a
new annihilation operator, based on the new equilibrium-based coordinate
system, is more useful in toroidal geometry. Finally, we conclude by
discussing the implications of the derivation, and the relationship to
computational approaches.

##  II Equilibrium-Based Frames of Reference

Frames of references are local coordinate systems that vary throughout space.
An example would be the rotating frame of reference describing the forces felt
by a person on a merry-go-round commonly used to explain the difference
between centripetal and centrifugal forces. Another well-known example from
mathematical physics is the Frenet frame of reference which constructs a local
set of basis vectors along a curve in terms of the tangent of the curve, the
curvature vector which is normal to the tangent, and the binormal vector
created from the cross-product of both the tangent and curvature. This creates
a local orthonormal set of basis vectors at a given point, but it is not a
general coordinate system because it cannot be used to label an arbitrary
point in space.

As discussed in the introduction, the most common frame of reference has the
magnetic field, the gradient of a flux function, and the binormal direction as
the local basis vectors. We term this frame the $BC$ frame and define it as:

$\mathbold{\nabla}\psi,\ \ \mathbold{B},\ \ \mathbold{C}$ (2)

where $\psi$ is the poloidal magnetic flux divided by $2\pi$, and
$\mathbold{C}$ is constructed to make the coordinate system orthogonal:

$\mathbold{C}\equiv\frac{\mathbold{\nabla}\psi\times\mathbold{B}}{\mid\mathbold{\nabla}\psi\mid^{2}}$
(3)

Unlike the Stix frame which normalizes the basis vectors so that they are
unit-less (and orders the parallel direction last for convenience) Wright
(1998), the normalizations chosen here seem to be the most useful as will be
shown.

In this paper, we define a new equilibrium-based frame of reference, the $JK$
frame

$\mathbold{\nabla}\psi,\ \ \mathbold{J},\ \ \mathbold{K}$ (4)

where $\mathbold{K}$ is likewise defined as

$\mathbold{K}\equiv\frac{\mathbold{\nabla}\psi\times\mathbold{J}}{P^{\prime}\mid\mathbold{\nabla}\psi\mid^{2}}.$
(5)

The $P^{\prime}$ denotes the derivative of the pressure with $\psi$. Our use
of this normalization is to make $\mathbold{K}\cdot\mathbold{B}=1$ which is
useful. Although using the other vector in the equilibrium relation,
$\mathbold{J}$, as the basis of an equilibrium-based frame seems obvious, to
our knowledge it has never been studied. In this paper, we show how for long
wavelength instabilities at least, it is useful to explicitly treat it as
such.

Using vector identities, it is easily seen that the other vector components
within each system may be expressed by

$\displaystyle\mathbold{B}=\mathbold{C}\times\mathbold{\nabla}\psi,\ \ \ \ \ \
\ \ \mathbold{\nabla}\psi=\frac{\mathbold{B}\times\mathbold{C}}{C^{2}}$ (6)
$\displaystyle\mathbold{J}=P^{\prime}\mathbold{K}\times\mathbold{\nabla}\psi,\
\ \ \ \ \ \ \
\mathbold{\nabla}\psi=\frac{\mathbold{J}\times\mathbold{K}}{P^{\prime}K^{2}}$
(7)

To express an arbitrary vector in terms of these local basis vectors, we
borrow notation from the curvilinear flux coordinate systems for simplicity.
If one considers these two vector systems as “contravariant-like” basis
vectors, $\mathbold{e}^{i},i\in 1,2,3$, then the concomitant “covariant-like”
basis vectors, $\mathbold{e}_{i}$ are given by
$\mathbold{e}_{1}=\epsilon^{123}{\mathcal{J}}\mathbold{e}^{2}\times\mathbold{e}^{3}$
where
${\mathcal{J}}=(\mathbold{e}^{1}\cdot\mathbold{e}^{2}\times\mathbold{e}^{3})^{-1}$
is the Jacobian of the system, and $\epsilon^{123}$ is the Levi-Cevita symbol.
The Jacobian is used as a normalization so that
$\mathbold{e}^{i}\cdot\mathbold{e}_{j}=\delta^{i}_{j}$ where $\delta^{i}_{j}$
is the Kronecker delta. Any vector is represented in terms of these basis
vectors by the following representation
$\mathbold{A}=\sum_{i}A^{i}\mathbold{e}_{i}=\sum_{i}A_{i}\mathbold{e}^{i},\
i\in\\{1,2,3\\}$. The ${A^{i}}$’s are the contravariant components and the
$A_{i}$’s are the covariant components. By the above relations, the components
may clearly be seen to be given by $A^{i}=\mathbold{A}\cdot\mathbold{e}^{i}$
and $A_{i}=\mathbold{A}\cdot\mathbold{e}_{i}$.

Using this notation for the basis vectors above, the inverse Jacobians of the
two systems are

$\displaystyle
B^{2}=C^{2}\mid\mathbold{\nabla}\psi\mid^{2}=\mathbold{\nabla}\psi\cdot\mathbold{B}\times\mathbold{C}$
(8)
$\displaystyle\frac{J^{2}}{P^{\prime}}=K^{2}\mid\mathbold{\nabla}\psi\mid^{2}=\mathbold{\nabla}\psi\cdot\mathbold{J}\times\mathbold{K}$
(9)

so that the covariant-like basis vectors are:

$\displaystyle\frac{\mathbold{\nabla}\psi}{\mid\mathbold{\nabla}\psi\mid^{2}},\
\ \ \ \ \frac{\mathbold{B}}{B^{2}},\ \ \ \ \frac{\mathbold{C}}{C^{2}}$ (10)
$\displaystyle\frac{\mathbold{\nabla}\psi}{\mid\mathbold{\nabla}\psi\mid^{2}},\
\ \ \ \frac{\mathbold{J}}{J^{2}},\ \ \ \ \frac{\mathbold{K}}{K^{2}}.$ (11)

The use of $\mid\mathbold{\nabla}\psi\mid^{2}$ in the definition of the
binormal vectors $\mathbold{C}$ and $\mathbold{K}$ is to make these covariant-
like basis vectors have this consistent form. An arbitrary vector
$\mathbold{f}$ then can be written as

$\displaystyle\mathbold{f}$
$\displaystyle=f^{\psi}\frac{\mathbold{\nabla}\psi}{\mid\mathbold{\nabla}\psi\mid^{2}}+f^{B}\frac{\mathbold{B}}{B^{2}}+f^{C}\frac{\mathbold{C}}{C^{2}}$
(12)
$\displaystyle=f^{\psi}\frac{\mathbold{\nabla}\psi}{\mid\mathbold{\nabla}\psi\mid^{2}}+f^{J}\frac{\mathbold{J}}{J^{2}}+f^{K}\frac{\mathbold{K}}{K^{2}}.$

Both the Stix frame and the Frenet frame use basis unit vectors and decompose
vectors into components that have the same units. Using these “contravariant-
like” components is perhaps un-intuitive, but they have long been used in
stability analyses Coppi et al. (1966). As we show in Section VI, this is
because they identify the correct components that equilibrate quickly due to
compressional Alfvèn waves.

Using force balance $\mathbold{J}\times\mathbold{B}=\mathbold{\nabla}P$, and
defining a variable $\sigma$ for the parallel current, the relationships
between the $BC$ and $JK$ frames are:

$\displaystyle\mathbold{J}=\sigma\mathbold{B}-P^{\prime}\frac{\mathbold{C}}{C^{2}},$
$\displaystyle\mathbold{K}=\frac{\sigma}{P^{\prime}}\mathbold{C}+\frac{\mathbold{B}}{B^{2}}$
(13) $\displaystyle\mathbold{B}=\sigma
B^{2}\frac{\mathbold{J}}{J^{2}}+\frac{\mathbold{K}}{K^{2}},$
$\displaystyle\mathbold{C}=\frac{\sigma
C^{2}}{P^{\prime}}\frac{\mathbold{K}}{K^{2}}-P^{\prime}\frac{\mathbold{J}}{J^{2}}$
(14)

The metrics of this transformation are

$\displaystyle\mathbold{J}\cdot\mathbold{B}=\sigma B^{2},$
$\displaystyle\mathbold{J}\cdot\mathbold{C}=-P^{\prime}$ (15)
$\displaystyle\mathbold{K}\cdot\mathbold{B}=1,$
$\displaystyle\mathbold{K}\cdot\mathbold{C}=\frac{\sigma C^{2}}{P^{\prime}}.$

The angle between the coordinate systems then is related to
$\sigma/P^{\prime}$, which is a flux function in cylindrical geometry but not
in toroidal geometry. To restore symmetry to the coordinate systems, one could
also use $P^{\prime}$ in the definition of $\mathbold{C}$, or equivalently,
use $\mathbold{\nabla}P$ as the radial vector, which has a certain elegance
since all three vectors of the equilibrium relation would then be used as the
basis vectors. The choice of normalizations used here seem to be the easiest
to use in practice however. The contravariant vector components are related by

$\displaystyle f^{C}=\frac{P^{\prime}}{J^{2}}(\sigma
B^{2}f^{K}-P^{\prime}f^{J}),$ $\displaystyle f^{B}=\frac{1}{K^{2}}(\sigma
C^{2}f^{J}+f^{K})$ (16) $\displaystyle f^{K}=\frac{1}{P^{\prime}B^{2}}(\sigma
B^{2}f^{C}+P^{\prime}f^{B}),$ $\displaystyle f^{J}=\frac{1}{C^{2}}(\sigma
C^{2}f^{B}-P^{\prime}f^{C}).$ (17)

The metrics presented here appear in the Mercier criterion; thus, one sees
even here that stability can be viewed as a complicated function of the angle
between these two frames.

We summarize the expressions in this section in Table 1. While these frames of
reference are useful, it is necessary to make use of geometry-based coordinate
systems to study the physics in the appropriate magnetic configuration. The
rest of this paper focuses on the doubly-periodic systems relevant to tokamaks
and symmetric stellarators.

##  III Flux Coordinates and the Relationship to Equilibrium-based Frames

While the local frames of reference are useful for getting at the physics of
the problem, they do not take advantage of the periodicity of the toroidal
systems that we wish to study. Following the work of DMS, the relationship
between the equilibrium-based frames of reference and the geometry-based
coordinate systems is explicitly examined. We start by briefly reviewing flux
coordinates Hazeltine and Meiss (1992), labeled in our notations as the
$\psi,\Theta,\zeta$ system. In this paper, only systems with at least one
degree of symmetry are considered with $\zeta$ being the ignorable coordinate.
This includes helically symmetric systems.

We assume nested flux surfaces and will consider various averages related to
those surfaces. A volume average of a quantity $Q$ at a surface $\psi$ is
denoted as

$\left<Q\right>_{V}\equiv\frac{1}{V(\psi)}\int Q\
dV=\frac{1}{V(\psi)}\int_{0}^{2\pi}\int_{0}^{2\pi}\int_{\psi_{min}}^{\psi}Q\
{\mathcal{J}}d\psi^{\prime}d\Theta d\zeta$ (18)

where $V(\psi)$ is the volume enclosed by a flux surface,
$V=\int\int\int{\mathcal{J}}d\psi^{\prime}d\Theta d\zeta$. The flux-surface
average of a quantity $Q$ at a surface $\psi$ is denoted as

$\left<Q\right>=\frac{1}{V^{\prime}(\psi)}\int_{0}^{2\pi}\int_{0}^{2\pi}Q\
{\mathcal{J}}d\Theta d\zeta$ (19)

where $V^{\prime}(\psi)$ is the derivative of $V$ with respect to $\psi$,
$V^{\prime}=\oint\oint{\mathcal{J}}d\Theta d\zeta$. Also, we define the theta
average of a quantity $Q$ as

$\left<Q\right>_{\Theta}=\frac{1}{2\pi}\int_{0}^{2\pi}Q\ d\Theta d\zeta.$ (20)

For axisymmetric quantities, this average is closely related to the flux-
surface average; i.e., if $Q$ is independent of $\zeta$, then
$\left<{\mathcal{J}}Q\right>_{\Theta}=V^{\prime}/(4\pi^{2})\left<Q\right>$.

###  III.1 Properties of divergence-free fields perpendicular to flux surfaces

Flux coordinates are usually derived using the properties of the divergence-
free nature of the magnetic field along with the definition of flux surface as
a surface perpendicular to the magnetic field. Our equilibrium-based local
frames of reference include two other divergence-free vectors, $\mathbold{J}$
and $\mid\mathbold{\nabla}\psi\mid^{2}vC$, so it is useful to first understand
general properties of vector-free fields perpendicular to flux surfaces. The
divergence-free property of $\mid\mathbold{\nabla}\psi\mid^{2}vC$ can se seen
using $\mathbold{J}\cdot\mathbold{\nabla}\psi=0$:

$\mathbold{\nabla}\psi\cdot\mathbold{J}=\mathbold{\nabla}\psi\cdot\mathbold{\nabla}\times\mathbold{B}=\mathbold{\nabla}\cdot\mathbold{B}\times\mathbold{\nabla}\psi=\mathbold{\nabla}\cdot\mid\mathbold{\nabla}\psi\mid^{2}\mathbold{C}=0.$
(21)

No similar relation can be found for $\mathbold{K}$. This makes this vector
fundamentally different from the other basis vectors.

Consider the general properties of a vector $\mathbold{V}$ with
$\mathbold{\nabla}\cdot\mathbold{V}=0$ and
$\mathbold{V}\cdot\mathbold{\nabla}\psi=0$. These properties allow
$\mathbold{V}$ to be written in a Clebsch representation:

$\mathbold{V}=\mathbold{\nabla}\beta_{V}\times\mathbold{\nabla}\psi.$ (22)

The most general form for $\beta_{V}$ in toroidal geometry is Boozer (1981):

$\beta_{V}=T_{V}\left(\psi\right)\zeta-
Z_{V}\left(\psi\right)\Theta+\hat{\beta}_{V}\left(\psi,\Theta,\zeta\right)$
(23)

where $\hat{\beta}_{V}$ is a periodic function of $\Theta$ and $\zeta$. This
allows $\mathbold{V}$ to be written as

$\displaystyle\mathbold{V}$
$\displaystyle=T_{V}\mathbold{\nabla}\zeta\times\mathbold{\nabla}\psi+Z_{V}\mathbold{\nabla}\Theta\times\mathbold{\nabla}\psi+\mathbold{\nabla}\hat{\beta}_{V}\times\mathbold{\nabla}\psi$
$\displaystyle={\mathcal{J}}^{-1}\left[\left(T_{V}+\frac{\partial\hat{\beta}_{V}}{\partial\zeta}\right)\mathbold{e}_{\Theta}-\left(Z_{V}-\frac{\partial\hat{\beta}_{V}}{\partial\Theta}\right)\mathbold{e}_{\zeta}\right]$

which is the covariant form for $\mathbold{V}$.

Using the expressions for the contravariant components of $\mathbold{V}$ in
Eq. (III.1), the poloidal flux of $\mathbold{V}$ may be easily calculated:

$\displaystyle\Gamma_{V}^{\Theta}\left(\psi\right)$
$\displaystyle\equiv\frac{1}{(2\pi)^{2}}\int\oint\mathbold{V}\cdot\mathbold{\nabla}\Theta\
{\mathcal{J}}d\zeta d\psi$ $\displaystyle=\frac{1}{2\pi}\int
T_{V}d\psi+\frac{1}{(2\pi)^{2}}\int\oint\frac{\partial\hat{\beta}_{V}}{\partial\zeta}d\zeta
d\psi.$

The last term is zero due to periodicity. Following a similar procedure for
the toroidal flux, $\Gamma_{V}^{\zeta}\left(\psi\right)\equiv
1/(2\pi)^{2}\int\oint\mathbold{V}\cdot\mathbold{\nabla}\zeta\
{\mathcal{J}}d\Theta d\psi$, one obtains the relationships

$T_{V}=2\pi{\Gamma_{V}^{\Theta}}^{\prime};\ \ \ \ \ \
Z_{V}=2\pi{\Gamma_{V}^{\zeta}}^{\prime}.$ (24)

which relates the contravariant components to the fluxes.

The fact that $\mathbold{V}\cdot\mathbold{\nabla}\psi=V^{\psi}=0$ allows the
concept of straightness to be defined. A vector $\mathbold{V}$ with zero
$\psi$ contravariant component is said to be straight if for each flux surface
defined by $\psi$, the ratio of the remaining contravariant components is
constant; i.e., $\mathbold{V}$ is straight if
$V^{\zeta}/V^{\Theta}=f\left(\psi\right)$. Looking at the contravariant
components for $\mathbold{V}$ in Eq. (III.1), it is apparent that this is
true, if and only if $\hat{\beta}_{V}=0$. In this case,

$\frac{V^{\zeta}}{V^{\Theta}}=\frac{Z_{V}}{T_{V}}=\frac{{\Gamma_{V}^{\zeta}}^{\prime}}{{\Gamma_{V}^{\Theta}}^{\prime}}=f\left(\psi\right).$
(25)

Calculations involving a straight vector are obviously greatly simplified
because of the simplification of the contravariant components to flux
functions.

To convert between $\mathbold{B}$ and $\mathbold{C}$, and $\mathbold{J}$ and
$\mathbold{K}$, the following identity for the concomitant vector,
$\mathbold{\nabla}\psi\times\mathbold{V}/\mid\mathbold{\nabla}\psi\mid^{2}$,
is useful:

$\displaystyle\frac{\mathbold{\nabla}\psi\times\mathbold{V}}{\mid\mathbold{\nabla}\psi\mid^{2}}$
$\displaystyle=\mathbold{\nabla}\psi\times\left(\mathbold{\nabla}\beta_{V}\times\mathbold{\nabla}\psi\right)$
$\displaystyle=\mathbold{\nabla}\beta_{V}-\frac{\mathbold{\nabla}\beta_{V}\cdot\mathbold{\nabla}\psi}{\mid\mathbold{\nabla}\psi\mid^{2}}\mathbold{\nabla}\psi$
$\displaystyle\equiv\mathbold{\nabla}_{t}\beta_{V}$

where in the last line a new operator is introduced for convenience. The
$\mathbold{\nabla}_{t}$ operator with respect to $\mathbold{\nabla}\psi$ is
analogous to the more familiar $\mathbold{\nabla}_{\perp}$ operator Kruger et
al. (1999) with respect to $\mathbold{B}$. The $t$ notation is used because
$\mathbold{\nabla}_{t}f$ is tangential to the constant $\psi$ surfaces. Using
this notation and the equation for $\beta_{V}$, Eq. (23), the concomitant
vector may be written as

$\frac{\mathbold{\nabla}\psi\times\mathbold{V}}{\mid\mathbold{\nabla}\psi\mid^{2}}=\left(T_{V}+\frac{\partial\hat{\beta}_{V}}{\partial\zeta}\right)\mathbold{\nabla}_{t}\zeta-\left(Z_{V}-\frac{\partial\hat{\beta}_{V}}{\partial\Theta}\right)\mathbold{\nabla}_{t}\Theta.$
(26)

Since $\mathbold{V}$ has a convenient contravariant form, the concomitant
vector will have a convenient covariant form, especially if $\mathbold{V}$ is
straight.

After the development of this formalism, the calculations for the specific
cases of
$\mathbold{V}\in\left\\{\mathbold{B},\mathbold{J},\mid\mathbold{\nabla}\psi\mid^{2}\mathbold{C}\right\\}$
are greatly simplified. In particular, our goal is to find expressions for
$\left\\{T,Z,\hat{\beta}\right\\}_{\left\\{B,C,J\right\\}}.$

###  III.2 Flux coordinates

In this paper, only flux coordinate systems where $\mathbold{B}$ is straight
($\hat{\beta}_{B}=0$) are considered. Because the radial coordinate has been
chosen to be $\psi\equiv 2\pi\Gamma_{B}^{\Theta}$, the components of
$\mathbold{B}$ have the nice form

$T_{B}=1;\ \ \ \ \ \ Z_{B}=q\left(\psi\right)$ (27)

where

$q\left(\psi\right)\equiv\frac{B^{\zeta}}{B^{\Theta}}=\frac{{\Gamma_{B}^{\zeta}}^{\prime}}{{\Gamma_{B}^{\Theta}}^{\prime}}.$
(28)

is the safety factor. This simple form for $\mathbold{B}$ allows
$\mathbold{C}$ to be written in the contravariant form

$\mathbold{C}=\mathbold{\nabla}_{t}\zeta-q\mathbold{\nabla}_{t}\Theta$ (29)

using Eq. (26).

Considering the covariant form of
$\mid\mathbold{\nabla}\psi\mid^{2}\mathbold{C}$ from Eq. (III.1) and using Eq.
(26) to find the contravariant form of $\mathbold{B}$ (Eq. (6)) yields

$\mathbold{B}=-\left(T_{C}+\frac{\partial\hat{\beta}_{C}}{\partial\zeta}\right)\mathbold{\nabla}_{t}\zeta+\left(Z_{C}-\frac{\partial\hat{\beta}_{C}}{\partial\Theta}\right)\mathbold{\nabla}_{t}\Theta.$
(30)

At this point, the ignorability of $\zeta$ is used to simplify the $\zeta$
covariant component of $\mathbold{B}$. We define

$F\left(\psi\right)\equiv-
T_{C}=-{\mathcal{J}}C^{\Theta}=-2\pi{\Gamma_{C}^{\Theta}}^{\prime}$ (31)

In axisymmetric systems, $F$ is known as the toroidal flux function
$F=RB_{toroidal}$. To find a useful expression for $Z_{C}$, the covariant form
of $\mathbold{B}$ (Eq. (30)) is dotted with the contravariant form of
$\mathbold{B}$ to give

$B^{2}={\mathcal{J}}^{-1}\left(Z_{C}-\frac{\partial\hat{\beta}_{C}}{\partial\Theta}\right)+{\mathcal{J}}^{-1}qF.$
(32)

For convenience, we define a new variable, $g$,

$g\equiv
Z_{C}-\frac{\partial\hat{\beta}_{C}}{\partial\Theta}={\mathcal{J}}B^{2}-qF.$
(33)

Using the $\Theta$ average allows us to annihilate the second term and solve:

$\displaystyle Z_{C}=\left<g\right>_{\Theta}\equiv G\left(\psi\right)$
$\displaystyle=\left<{\mathcal{J}}B^{2}\right>_{\Theta}-qF=\frac{V^{\prime}}{4\pi^{2}}\left<B^{2}\right>-qF;$
$\displaystyle\frac{\partial\hat{\beta}_{C}}{\partial\Theta}$
$\displaystyle=\left<{\mathcal{J}}B^{2}\right>_{\Theta}-{\mathcal{J}}B^{2}=\frac{V^{\prime}}{4\pi^{2}}\left<B^{2}\right>-{\mathcal{J}}B^{2}.$

Evidently, to make $\mathbold{C}$ straight, ${\mathcal{J}}B^{2}$ needs to be a
flux function. A coordinate system with this property is known as the Boozer
flux coordinate system Boozer (1981, 1982). Using these definitions, the
covariant form of $\mathbold{B}$ can be written as
$\mathbold{B}=F\mathbold{\nabla}_{t}\zeta+g\mathbold{\nabla}_{t}\Theta$.

Taking the curl of $\mathbold{B}$ to give the current density gives

$\mathbold{J}=-F^{\prime}\mathbold{\nabla}\zeta\times\mathbold{\nabla}\psi-\left[\frac{\partial
g}{\partial\psi}-\frac{\partial}{\partial\Theta}\left(F\frac{g^{\psi\zeta}}{g^{\psi\psi}}+g\frac{g^{\psi\Theta}}{g^{\psi\psi}}\right)\right]\mathbold{\nabla}\Theta\times\mathbold{\nabla}\psi.$
(34)

The expression for $T_{J}$ is simply $T_{J}=-F^{\prime}$; thus, the function
$F$ introduced above is proportional to the poloidal current enclosed by a
flux surface as seen from Eq. (24) ($F=-2\pi\Gamma_{J}^{\Theta}$). Similarly,
one can show
$G\left(\psi\right)=\left<g\right>_{\Theta}=2\pi\Gamma_{J}^{\zeta}$ is
proportional to the toroidal current enclosed by a flux surface.

The above expression for $\mathbold{J}$ was derived using Ampere’s law. The
expression for $Z_{J}$ can be simplified by deriving another covariant form of
$\mathbold{J}$ using Eq. (13) and the covariant forms for $\mathbold{B}$ and
$\mathbold{C}$ already derived. Using the covariant forms yields

$\mathbold{J}=\left(\sigma+\frac{P^{\prime}F}{B^{2}}\right)\mathbold{\nabla}\zeta\times\mathbold{\nabla}\psi-\left(q\sigma-\frac{P^{\prime}g}{B^{2}}\right)\mathbold{\nabla}\Theta\times\mathbold{\nabla}\psi.$
(35)

Equating the $\Theta$ contravariant components yields an expression for the
parallel current:

$\sigma=-F^{\prime}-\frac{P^{\prime}F}{B^{2}}.$ (36)

This relation allows us to rewrite the $\zeta$ contravariant component of
$\mathbold{J}$ in Eq. (35) as

$h\equiv
Z_{J}-\frac{\partial\hat{\beta}_{J}}{\partial\Theta}=q\sigma-\frac{P^{\prime}g}{B^{2}}=-\left(P^{\prime}{\mathcal{J}}+qF^{\prime}\right)$
(37)

where $h$ is defined analogously to $g$. Using the $\Theta$ average again
allows us to annihilate the second term and solve for each term:

$\displaystyle Z_{J}=\left<h\right>_{\Theta}=G^{\prime}\left(\psi\right)$
$\displaystyle=-\left(P^{\prime}\left<{\mathcal{J}}\right>_{\Theta}+qF^{\prime}\right)=-\left(\frac{V^{\prime}}{4\pi^{2}}P^{\prime}+qF^{\prime}\right)=2\pi\Gamma_{J}^{\zeta};$
$\displaystyle\frac{\partial\hat{\beta}_{J}}{\partial\Theta}$
$\displaystyle=P^{\prime}\left({\mathcal{J}}-\left<{\mathcal{J}}\right>_{\Theta}\right)=P^{\prime}\left({\mathcal{J}}-\frac{V^{\prime}}{4\pi^{2}}\right).$

Evidently, to make $\mathbold{J}$ straight, the Jacobian needs to be a flux
function. A coordinate system with this property is known as the Hamada flux
coordinate system.

Equating the $\zeta$ contravariant components yields

$\frac{\partial
g}{\partial\psi}+\frac{\partial}{\partial\Theta}\left(F\frac{g^{\psi\zeta}}{g^{\psi\psi}}+g\frac{g^{\psi\Theta}}{g^{\psi\psi}}\right)=h$
(38)

which is a generalized Grad-Shafranov equation.

The contravariant form of $\mathbold{J}$ can be derived using Eq. (13) and the
covariant forms for $\mathbold{B}$ and $\mathbold{C}$ already derived. Using
the contravariant forms yields

$\mathbold{J}=\left(\sigma
F-\frac{P^{\prime}}{C^{2}}\right)\mathbold{\nabla}_{t}\zeta+\left(\sigma
g+\frac{qP^{\prime}}{C^{2}}\right)\mathbold{\nabla}_{t}\Theta.$ (39)

Using the contravariant and covariant forms of $\mathbold{J}$, expressions for
$\mathbold{K}$ may be easily derived using Eq. (5). The results are

$\mathbold{K}=\frac{-1}{P^{\prime}\mid\mathbold{\nabla}\psi\mid^{2}}\left[\left(\sigma
F-\frac{P^{\prime}}{C^{2}}\right)\mathbold{\nabla}\zeta\times\mathbold{\nabla}\psi\left(\sigma
g+q\frac{P^{\prime}}{C^{2}}\right)\mathbold{\nabla}\Theta\times\mathbold{\nabla}\psi\right]$
(40)

and

$\mathbold{K}=\frac{-1}{P^{\prime}}\left[F^{\prime}\mathbold{\nabla}_{t}\zeta+h\mathbold{\nabla}_{t}\Theta\right]$
(41)

We summarize the expressions which were derived above for all four of our
basis vectors in Table 1.

###  III.3 Helical angle coordinates

In addition to the $\psi,\Theta,\zeta$ coordinate system, another useful
coordinate system is the $\psi,\Theta,\alpha$ coordinate system where $\alpha$
is the Clebsch angle for $\mathbold{B}$

$\alpha\equiv\beta_{B}=\zeta-q\Theta.$ (42)

such that $\mathbold{B}=\mathbold{\nabla}\alpha\times\mathbold{\nabla}\psi$.
If $\zeta$ is an ignorable coordinate, then $\alpha$ is an ignorable
coordinate also. The primary advantage of this coordinate system is that the
tangential gradient of $\alpha$ is equal to $\mathbold{C}$ as can easily be
seen by Eq. (29):

$\mathbold{C}=\mathbold{\nabla}_{t}\alpha=\mathbold{\nabla}\alpha-\frac{g^{\psi\alpha}}{g^{\psi\psi}}\mathbold{\nabla}\psi.$
(43)

Expressions for the basis vectors using this coordinate system can be easily
evaluated. We summarize the expressions which were derived above for all four
of our basis vectors in Table 1.

Examination of the equation for $\mathbold{\nabla}_{t}\alpha$, Eq. (43), in
more detail, shows that the last term may be written as

$\displaystyle\frac{g^{\psi\alpha}}{g^{\psi\psi}}$
$\displaystyle=\frac{1}{\mid\mathbold{\nabla}\psi\mid^{2}}\left(\mathbold{\nabla}\zeta\cdot\mathbold{\nabla}\psi-\mathbold{\nabla}(q\Theta)\cdot\mathbold{\nabla}\psi\right)=-\left(q^{\prime}\Theta+\frac{qg^{\psi\Theta}-g^{\psi\zeta}}{g^{\psi\psi}}\right)$
$\displaystyle=-\left(q^{\prime}\Theta+\mathcal{R}(\psi,\Theta)\right).$

The $\mathcal{R}$ term is known as the residual shear and will be discussed in
the next section. As evidenced by the $q^{\prime}\Theta$ term in the previous
equation, $\alpha$ is not a periodic coordinate. While this coordinate system
is still useful for short wavelength modes, it is often more useful to have a
coordinate system which is periodic, yet closely related to $\alpha$. We
define the coordinate $u$ as

$u\equiv\zeta-q_{s}\Theta=\alpha+(q-q_{s})\Theta$ (44)

where $q_{s}$ is the value of $q$ at a particular surface. For a vector
$\mathbold{F}\perp\mathbold{\nabla}\psi$,
$\mathbold{F}\cdot\mathbold{\nabla}u=\mathbold{F}\cdot\mathbold{\nabla}\alpha+(q-q_{0})\mathbold{F}\cdot\mathbold{\nabla}\Theta$
which may be used to find the $u$ contravariant components and following a
similar derivation as above, the covariant components as well. The summary of
these relationships is shown in vectors in Table 1 as well.

###  III.4 Special cases of flux coordinate systems

Until this point, $\Theta$ and $\zeta$ are only specified such that the
magnetic field lines are straight:

$\frac{\mathbold{B}\cdot\mathbold{\nabla}\zeta}{\mathbold{B}\cdot\mathbold{\nabla}\Theta}=q(\psi).$
(45)

This is one constraint on two quantities. To completely specify a straight-
field line coordinate system, another requirement is needed; generally the
Jacobian is specified.

Previously we mentioned two specific coordinates, the Boozer and Hamada
coordinate systems, and here we describe their Jacobians and properties.
Recall that Hamada coordinates specify the constraint that $\mathbold{J}$ is
straight,

$\frac{\mathbold{J}\cdot\mathbold{\nabla}\zeta}{\mathbold{J}\cdot\mathbold{\nabla}\Theta}=f(\psi),$
(46)

and Boozer coordinates specify that the magnetic binormal vector is straight
111A useful mnemonic: Allen Boozer made C straight as simple as ABC. Hamada Is
J straight as simple as HIJ.,

$\frac{\mathbold{C}\cdot\mathbold{\nabla}\zeta}{\mathbold{C}\cdot\mathbold{\nabla}\Theta}=f(\psi).$
(47)

In addition to these two coordinate systems, we also discuss the “symmetry
coordinates”. Practical aspects related to calculating the coordinate systems
numerically may be found in DMS Dewar et al. (1984).

####  III.4.1 Hamada Coordinate System

Hamada coordinates are elegant because it makes the current density field
lines straight as well, thus continuing the analogous relationships between
$\mathbold{B}$ and $\mathbold{J}$. From Eq. (37), we see that

$\frac{\mathbold{J}\cdot\mathbold{\nabla}\zeta}{\mathbold{J}\cdot\mathbold{\nabla}\Theta}=\frac{h}{F^{\prime}}=q+{\mathcal{J}}\frac{P^{\prime}}{F^{\prime}}$
(48)

To make these straight, we need to make the Jacobian a flux-function. This
flux function may be found from the relationship:

$\left<{\mathcal{J}}^{-1}\right>=\frac{4\pi^{2}}{V^{\prime}}.$ (49)

which implies that the Jacobian for the Hamada coordinate system may be
written as

${\mathcal{J}}_{H}=\frac{V^{\prime}}{4\pi^{2}}.$ (50)
$\mathbold{J}=-F^{\prime}(\psi)\mathbold{\nabla}\zeta\times\mathbold{\nabla}\psi-G^{\prime}(\psi)\mathbold{\nabla}\Theta\times\mathbold{\nabla}\psi.$
(51)

with $G^{\prime}/F^{\prime}=q+4\pi^{2}/(F^{\prime}V^{\prime})$.

$-qF^{\prime}-G^{\prime}=\frac{V^{\prime}}{4\pi^{2}}P^{\prime}.$ (52)

which is the Grad-Shafranov equation in Hamada coordinates (using $\psi$ as
the dependent variable). This relationship is used extensively in Greene and
Johnson Greene and Johnson (1962), and in the original derivation of the inner
layer equations of GGJ.

####  III.4.2 Boozer Coordinate System

Boozer coordinates have nice properties because $\mathbold{C}$ is straight as
well as $\mathbold{B}$. From Eq. (33), we see that

$\frac{\mathbold{C}\cdot\mathbold{\nabla}\zeta}{\mathbold{C}\cdot\mathbold{\nabla}\Theta}=\frac{-g}{F}=q-\frac{{\mathcal{J}}B^{2}}{F}$
(53)

To make these straight, we need to have

${\mathcal{J}}_{B}B^{2}=G(\psi)$ (54)

To find this flux function, we have

$\left<\frac{{\mathcal{J}}^{-1}_{B}}{B^{2}}B^{2}\right>=\frac{{\mathcal{J}}^{-1}_{B}}{B^{2}}\left<B^{2}\right>=\frac{4\pi^{2}}{V^{\prime}}$
(55)

or

${\mathcal{J}}_{B}B^{2}=G(\psi)=\frac{V^{\prime}\left<B^{2}\right>}{4\pi^{2}}.$
(56)

####  III.4.3 Symmetry Coordinates

In an axisymmetric system, the Boozer and Hamada surfaces of constant $\zeta$
will not correspond to surfaces of constant $\phi$ where $\phi$ is the
“symmetry” angle. That is, given a cylindrical $R,\phi,Z$ coordinate system,
one is interested in the case of $\zeta=-\phi$ 222The minus sign is to enable
the right handed-ness to remain in the $(R,\phi,Z)$ and $(\psi,\Theta,\zeta)$
coordinate systems.

In this coordinate system,
$\mathbold{\nabla}\zeta\cdot\mathbold{\nabla}\zeta=1/R^{2}$,
$\mathbold{\nabla}\zeta\cdot\mathbold{\nabla}\psi=0$, and
$\mathbold{\nabla}\zeta\cdot\mathbold{\nabla}\Theta=0$. This means that the
covariant and contravariant components are parallel:

$\mathbold{\nabla}\Theta\times\mathbold{\nabla}\psi=-{\mathcal{J}}^{-1}\mathbold{e}_{\zeta}={\mathcal{J}}^{-1}R^{2}\mathbold{\nabla}\zeta$
(57)

which allows us to write the magnetic field as

$\mathbold{B}={\mathcal{J}}^{-1}qR^{2}\mathbold{\nabla}\zeta+\mathbold{\nabla}\zeta\times\mathbold{\nabla}\psi.$
(58)

Because of the orthogonality of the $\zeta$ surfaces and $\psi$ surfaces in
this coordinate system, the tangential derivative is the same such that we can
write

$\mathbold{B}=F\mathbold{\nabla}\zeta+g\mathbold{\nabla}_{t}\Theta.$ (59)

Equating the two expressions for the covariant components allows us to write
the familiar form for the symmetric angle Freidberg (1987); Hazeltine and
Meiss (1992),

$\mathbold{B}=F\mathbold{\nabla}\zeta+\mathbold{\nabla}\zeta\times\mathbold{\nabla}\psi,$
(60)

with the Jacobian in this system given as

${\mathcal{J}}=\frac{qR^{2}}{F}.$ (61)

This is sometimes referred to as the PEST angle because the first paper
describing it was in a paper by Grimm, Greene, and Johnson Grimm et al. (1976)
describing the PEST ideal MHD code.

##  IV Geometric Quantities and Frenet’s Local Frame of Reference

Expressions using the quantities of the first section such as $\sigma$ or
$P^{\prime}$, or the metric elements of the flux coordinate system, can be
difficult to use for gaining physical intuition. In this section, common
geometric quantities are introduced and related to the reference frames of the
previous sections.

Considerations of a particle tracing out a curve in three-dimensional
Euclidean space lead to studies in the mid-19th century of the geometry of
curves. This led to the development of the Frenet frame of reference, which
uses as its local basis vectors the tangent of the curve, the curvature vector
which is normal to the tangent, and the binormal vector. The fundamental
theorem of space curves states that in three-dimensional space every curve
with non-zero curvature has its shape completely determined by curvature and
torsion. Because of the importance of the $BC$ frame of reference, the
relationship of this frame to the curvature and torsion has been considered
for some time Dewar et al. (1984); Hegna (2000). That is, while the curvature
and torsion of $\mathbold{J}$ might also be of interest, we will only consider
the curvature and torsion of the magnetic field. Our goal is to express the
Frenet quantities of torsion and curvature in terms of our quantities from the
$BC$ and $JK$ reference frames.

Hegna Hegna (2000) related the quantities of torsion to the normalized version
of our $BC$ frame (also known as the Stix frame). The relationships are:

$\displaystyle\left(\hat{\mathbold{b}}\cdot\mathbold{\nabla}\right)\hat{\mathbold{n}}$
$\displaystyle=-\kappa_{n}\hat{\mathbold{b}}-\tau_{n}\hat{\mathbold{n}}\times\hat{\mathbold{b}}$
$\displaystyle\left(\hat{\mathbold{b}}\cdot\mathbold{\nabla}\right)\hat{\mathbold{b}}$
$\displaystyle=\kappa_{n}\hat{\mathbold{n}}-\kappa_{g}\hat{\mathbold{n}}\times\hat{\mathbold{b}}$
$\displaystyle\left(\hat{\mathbold{b}}\cdot\mathbold{\nabla}\right)\hat{\mathbold{n}}\times\hat{\mathbold{b}}$
$\displaystyle=\tau_{n}\hat{\mathbold{n}}+\kappa_{g}\hat{\mathbold{b}}$

The curvature, rather than acting as the normal component in the Frenet frame,
is now decomposed into the normal curvature, $\kappa_{n}$ that is
perpendicular to a flux surface, and the geodesic curvature, $\kappa_{g}$,
that is parallel to the magnetic binormal direction. In these formulas, the
only component of torsion that matters is the normal torsion, $\tau_{n}$.
Here, we briefly review the relationship of these quantities to our frames,
and derive other important geometric relationships.

###  IV.1 Magnetic Shear

Although a single curve can be described by curvature and torsion, it does not
convey how neighboring curves behave. In magnetically confined systems, how a
magnetic field line moves relative to a neighboring field line is important
for modes that have finite width. The most useful measure of this variation is
the magnetic shear, $\mathcal{S}$, which is defined as Greene and Chance
(1981); Johnson and Greene (1967):

$\mathcal{S}\equiv-\frac{1}{C^{2}}\mathbold{C}\cdot\mathbold{\nabla}\times\mathbold{C}.$
(62)

Definition of local shear often differs from this definition by either a minus
sign Hegna (2000) or the factor of $C^{2}$ Greene and Chance (1981); Johnson
and Greene (1967).

Although the $\mathbold{C}$ component is generally the most useful component
of $\mathbold{\nabla}\times\mathbold{C}$, we will consider the other
components as well. The normal component can be shown to be zero using the
divergence-free nature of $\mathbold{B}$:

$\mathbold{\nabla}\cdot\mathbold{B}=\mathbold{\nabla}\cdot\left(\mathbold{C}\times\mathbold{\nabla}\psi\right)=\mathbold{\nabla}\psi\cdot\mathbold{\nabla}\times\mathbold{C}=0.$
(63)

Using the relation between $\mathbold{C}$ and $\alpha$ (Eq. 43), the curl of
$\mathbold{C}$ may be written as

$\mathbold{\nabla}\times\mathbold{C}=\left(\mid\mathbold{\nabla}\psi\mid^{2}\mathbold{C}\cdot\mathbold{\nabla}\frac{g^{\psi\alpha}}{g^{\psi\psi}}\right)\frac{\mathbold{B}}{B^{2}}-\left(\mathbold{B}\cdot\mathbold{\nabla}\frac{g^{\psi\alpha}}{g^{\psi\psi}}\right)\frac{\mathbold{C}}{C^{2}}$
(64)

Using Eq. (III.3), the local magnetic shear then is given by

$C^{2}\mathcal{S}=-\mathbold{B}\cdot\mathbold{\nabla}\frac{g^{\psi\alpha}}{g^{\psi\psi}}={\mathcal{J}}^{-1}\left(q^{\prime}+\frac{\partial\mathcal{R}}{\partial\Theta}\right)$
(65)

Thus $q^{\prime}$ is termed the global shear, and $\mathcal{R}$ is termed the
integrated residual shear because its derivative give the variation of the
local shear within a flux surface.

Performing a similar calculation for the parallel component of the curl of
$\mathbold{C}$ allows us to write the total curl as

$\mathbold{\nabla}\times\mathbold{C}=-C^{2}\left[F\mathcal{S}\frac{\mathbold{B}}{B^{2}}+\mathcal{S}\frac{\mathbold{C}}{C^{2}}\right]$
(66)

###  IV.2 Torsion

The normal component of the torsion is given by

$\tau_{n}=\hat{\mathbold{b}}\cdot\mathbold{\nabla}\frac{\mathbold{C}}{\mid\mathbold{C}\mid}\cdot\frac{\mathbold{\nabla}\psi}{\mid\mathbold{\nabla}\psi\mid}=\frac{1}{B^{2}}\mathbold{B}\cdot\mathbold{\nabla}\mathbold{C}\cdot\mathbold{\nabla}\psi.$
(67)

Using the identity

$\displaystyle
2\mathbold{a}\cdot\mathbold{\nabla}\mathbold{b}\cdot\mathbold{c}$
$\displaystyle=\mathbold{a}\cdot\mathbold{\nabla}\left(\mathbold{b}\cdot\mathbold{c}\right)-\mathbold{b}\cdot\mathbold{\nabla}\left(\mathbold{c}\cdot\mathbold{a}\right)+\mathbold{c}\cdot\mathbold{\nabla}\left(\mathbold{a}\cdot\mathbold{b}\right)$
$\displaystyle-\mathbold{a}\times\mathbold{b}\cdot\mathbold{\nabla}\times\mathbold{c}+\mathbold{b}\times\mathbold{c}\cdot\mathbold{\nabla}\times\mathbold{a}-\mathbold{c}\times\mathbold{a}\cdot\mathbold{\nabla}\times\mathbold{b}$

the normal torsion can be related to other geometric quantities

$\displaystyle
2\mathbold{B}\cdot\mathbold{\nabla}\mathbold{C}\cdot\mathbold{\nabla}\psi$
$\displaystyle=\mathbold{C}\times\mathbold{\nabla}\psi\cdot\mathbold{\nabla}\times\mathbold{B}-\mathbold{\nabla}\psi\times\mathbold{B}\cdot\mathbold{\nabla}\times\mathbold{C}$
$\displaystyle=\mathbold{B}\cdot\mathbold{J}-\mid\mathbold{\nabla}\psi\mid^{2}\mathbold{C}\cdot\mathbold{\nabla}\times\mathbold{C}$

or, relating the normal component of the torsion to the above quantity:

$2\tau_{n}=\sigma+\mathcal{S}.$ (68)

This relationship is discussed in both Ref. Greene (1996) and Hegna (2000) and
can be viewed as a type of Grad-Shafranov equation.

###  IV.3 Curvature

The curvature of the magnetic field line is defined as

$\mathbold{\kappa}\equiv\left(\hat{\mathbold{b}}\cdot\mathbold{\nabla}\right)\hat{\mathbold{b}}.$
(69)

Using the equilibrium force balance
($\mathbold{J}\times\mathbold{B}=\mathbold{\nabla}P$) and Ampere’s law
($\mathbold{J}=\mathbold{\nabla}\times\mathbold{B}$), the magnetic curvature
can be expressed as

$\mathbold{\kappa}=\frac{1}{B^{2}}\mathbold{\nabla}_{\perp}\left(P+\frac{B^{2}}{2}\right).$
(70)

Because the curvature vector is perpendicular to the magnetic field, we can
write the curvature vector in the contravariant form:
$\mathbold{\kappa}=\kappa_{\psi}\mathbold{\nabla}\psi+\kappa_{C}\mathbold{C}$.
The covariant components are trivially related to the normal and geodesic
curvatures by normalization constants;
$\kappa_{n}=\kappa_{\psi}\mid\mathbold{\nabla}\psi\mid$,
$\kappa_{g}=\kappa_{C}\mid\mathbold{C}\mid$. Using Eq. 70, expressions for the
covariant components of the curvature may be easily calculated:

$\kappa_{\psi}=\frac{1}{B^{2}}\frac{\partial}{\partial\psi}\left(P+\frac{B^{2}}{2}\right)+\frac{g^{\psi\Theta}}{g^{\psi\psi}}\frac{1}{B^{2}}\frac{\partial}{\partial\Theta}\left(\frac{B^{2}}{2}\right)$
(71)
$\kappa_{C}=\frac{1}{C^{2}B^{2}}\mathbold{C}\cdot\mathbold{\nabla}\Theta\frac{\partial}{\partial\Theta}\left(\frac{B^{2}}{2}\right)=-F{\mathcal{J}}^{-1}\frac{1}{B^{4}}\frac{\partial}{\partial\Theta}\left(\frac{B^{2}}{2}\right).$
(72)

This will be used in the next section.

An important relation for the geodesic curvature can be derived using
$\mathbold{\nabla}\cdot\mathbold{J}=0$ and Eqs. (13), 21, 72:

$\displaystyle\mathbold{\nabla}\cdot\mathbold{J}$
$\displaystyle=\mathbold{\nabla}\cdot\left(\sigma\mathbold{B}-P^{\prime}\frac{\mathbold{C}}{C^{2}}\right)=\mathbold{B}\cdot\mathbold{\nabla}\sigma-\mathbold{\nabla}\cdot\frac{P^{\prime}}{B^{2}}\mid\mathbold{\nabla}\psi\mid^{2}\mathbold{C}=\mathbold{B}\cdot\mathbold{\nabla}\sigma+P^{\prime}2\kappa_{C}$
(73)

or,

$2\kappa_{C}=-\mathbold{B}\cdot\mathbold{\nabla}\frac{\sigma}{P^{\prime}}.$
(74)

Thus we can see that the angle between the $BC$ coordinate system and the $JK$
coordinate system which was shown to be related to $\sigma/P^{\prime}$ is
related to the geodesic curvature.

###  IV.4 Curl($\mathbold{K}$)

Due to its relationship to the local magnetic shear, the curl of
$\mathbold{C}$ is an important quantity. The curl of $\mathbold{K}$ is also an
important quantity and here we examine its components. Similar to the curl of
$\mathbold{C}$, the curl of $\mathbold{K}$ can be shown to be perpendicular to
$\mathbold{\nabla}\psi$ using $\mathbold{\nabla}\cdot\mathbold{J}=0$ which
allows us to decompose the curl using $J$ and $K$ components, or $B$ and $C$
components. The $K$ component is the definition of the shear of the current
density field. Experience shows that this is not a very useful quantity, so we
opt for examining the $BC$ components of the curl. Calculating the curl of
$\mathbold{K}$:

$\displaystyle\mathbold{\nabla}\times\mathbold{K}$
$\displaystyle=\mathbold{\nabla}\times\frac{\mathbold{B}}{B^{2}}+\mathbold{\nabla}\times\frac{\sigma}{P^{\prime}}\mathbold{C}=\sigma\frac{\mathbold{B}}{B^{2}}+\frac{\mathbold{B}\times\mathbold{\nabla}(P+B^{2})}{B^{4}}-\mathbold{C}\times\mathbold{\nabla}\frac{\sigma}{P^{\prime}}+\frac{\sigma}{P^{\prime}}\mathbold{\nabla}\times\mathbold{C}$
(75)

The parallel component is

$\mathbold{B}\cdot\mathbold{\nabla}\times\mathbold{K}=\sigma-B^{2}\frac{\partial}{\partial\psi}(\frac{\sigma}{P^{\prime}})-\frac{\sigma
C^{2}}{P^{\prime}}F\mathcal{S}$ (76)

where we have used Eq. (66). The perpendicular component is

$\mathbold{C}\cdot\mathbold{\nabla}\times\mathbold{K}=-2\kappa_{\psi}+\frac{P^{\prime}}{B^{2}}-\frac{\sigma
C^{2}}{P^{\prime}}\mathcal{S}$ (77)

Using the generalized Grad-Shafranov equation (Eq. (38)) and equations for the
normal curvature (Eq. (71)) and shear (Eqs. (65) and (III.3)), one can derive

$\mathbold{C}\cdot\mathbold{\nabla}\times\mathbold{K}={\mathcal{J}}^{-1}\frac{\partial}{\partial\psi}{\mathcal{J}}+{\mathcal{J}}^{-1}\frac{\partial}{\partial\Theta}\left({\mathcal{J}}\frac{g^{\psi\Theta}}{g^{\psi\psi}}\right)+\frac{F^{\prime}C^{2}}{P^{\prime}}\mathcal{S}$
(78)

The importance of this component is discussed in the derivation of the inner
layer equations.

A partial listing of the expressions in this section are summarized in Table
1.

Summary of equilibrium-based local frames of reference
---
$\mathbold{\nabla}\psi=\frac{\mathbold{B}\times\mathbold{C}}{C^{2}}$ | $\mathbold{\nabla}\psi=\frac{\mathbold{J}\times\mathbold{K}}{P^{\prime}K^{2}}$
$\mathbold{B}=\mathbold{C}\times\mathbold{\nabla}\psi$ | $\mathbold{J}=P^{\prime}\mathbold{K}\times\mathbold{\nabla}\psi$
$\mathbold{C}\equiv\frac{\mathbold{\nabla}\psi\times\mathbold{B}}{\mid\mathbold{\nabla}\psi\mid^{2}}$ | $\mathbold{K}\equiv\frac{\mathbold{\nabla}\psi\times\mathbold{J}}{P^{\prime}\mid\mathbold{\nabla}\psi\mid^{2}}$
$C^{2}=\frac{B^{2}}{\mid\mathbold{\nabla}\psi\mid^{2}}$ | $K^{2}=\frac{J^{2}}{P^{\prime 2}\mid\mathbold{\nabla}\psi\mid^{2}}\rule[-5.16663pt]{0.0pt}{0.0pt}$
$\mathbold{B}=\sigma B^{2}\frac{\mathbold{J}}{J^{2}}+\frac{\mathbold{K}}{K^{2}}$ | $\mathbold{J}=\sigma\mathbold{B}-P^{\prime}\frac{\mathbold{C}}{C^{2}},\rule{0.0pt}{19.80551pt}$
$\mathbold{C}=\frac{\sigma C^{2}}{P^{\prime}}\frac{\mathbold{K}}{K^{2}}-P^{\prime}\frac{\mathbold{J}}{J^{2}}$ | $\mathbold{K}=\frac{\sigma}{P^{\prime}}\mathbold{C}+\frac{\mathbold{B}}{B^{2}}\rule[-5.16663pt]{0.0pt}{0.0pt}$
$\mathbold{f}=f^{\psi}\frac{\mathbold{\nabla}\psi}{\mid\mathbold{\nabla}\psi\mid^{2}}+f^{B}\frac{\mathbold{B}}{B^{2}}+f^{C}\frac{\mathbold{C}}{C^{2}}$ | $=f^{\psi}\frac{\mathbold{\nabla}\psi}{\mid\mathbold{\nabla}\psi\mid^{2}}+f^{J}\frac{\mathbold{J}}{J^{2}}+f^{K}\frac{\mathbold{K}}{K^{2}}.\rule{0.0pt}{19.80551pt}$
$f^{B}=\frac{1}{K^{2}}(\sigma C^{2}f^{J}+f^{K})$ | $f^{J}=\frac{1}{C^{2}}(\sigma C^{2}f^{B}-P^{\prime}f^{C})$
$f^{C}=\frac{P^{\prime}}{J^{2}}(\sigma B^{2}f^{K}-P^{\prime}f^{J})$ | $f^{K}=\frac{1}{P^{\prime}B^{2}}(\sigma B^{2}f^{C}+P^{\prime}f^{B})$

Summary of flux coordinate’s relationship to equilibrium-based frames
---
$\mathbold{B}$ |  $=\mathbold{\nabla}\zeta\times\mathbold{\nabla}\psi-q\mathbold{\nabla}\Theta\times\mathbold{\nabla}\psi$ | flux coordinates: |
|  $=F\mathbold{\nabla}_{t}\zeta+g\mathbold{\nabla}_{t}\Theta$ | $q=q(\psi)$ |  $\Rightarrow\mathbold{B}$ is straight
$\mathbold{C}$ |  $=\frac{-1}{\mid\mathbold{\nabla}\psi\mid^{2}}\left[F\mathbold{\nabla}\zeta\times\mathbold{\nabla}\psi+g\mathbold{\nabla}\Theta\times\mathbold{\nabla}\psi\right]$ |  $g\equiv{\mathcal{J}}B^{2}-qF$; |  $G\left(\psi\right)=\left<g\right>_{\Theta}=\frac{V^{\prime}}{4\pi^{2}}\left<B^{2}\right>-qF$
|  $=\mathbold{\nabla}_{t}\zeta-q\mathbold{\nabla}_{t}\Theta$ | Boozer: $g=G$ |  $\Rightarrow\mid\mathbold{\nabla}\psi\mid^{2}\mathbold{C}$ is straight
$\mathbold{J}$ |  $=-F^{\prime}\mathbold{\nabla}\zeta\times\mathbold{\nabla}\psi-h\mathbold{\nabla}\Theta\times\mathbold{\nabla}\psi$ |  $h\equiv-{\mathcal{J}}P^{\prime}-qF^{\prime}$; |  $G^{\prime}\left(\psi\right)=\left<h\right>_{\Theta}=-\frac{V^{\prime}}{4\pi^{2}}P^{\prime}-qF^{\prime}$
|  $=\left(\sigma F-\frac{P^{\prime}}{C^{2}}\right)\mathbold{\nabla}_{t}\zeta+\left(\sigma g+\frac{qP^{\prime}}{C^{2}}\right)\mathbold{\nabla}_{t}\Theta$ | Hamada: $h=G^{\prime}$ |  $\Rightarrow\mathbold{J}$ is straight
$\mathbold{K}$ |  $=\frac{-1}{P^{\prime}\mid\mathbold{\nabla}\psi\mid^{2}}\left[\left(\sigma F-\frac{P^{\prime}}{C^{2}}\right)\mathbold{\nabla}\zeta\times\mathbold{\nabla}\psi\right.$ |  |
|  $\ \ \ \ \left.+\left(\sigma g+q\frac{P^{\prime}}{C^{2}}\right)\mathbold{\nabla}\Theta\times\mathbold{\nabla}\psi\right]$ |  |
|  $=\frac{-1}{P^{\prime}}\left[F^{\prime}\mathbold{\nabla}_{t}\zeta+h\mathbold{\nabla}_{t}\Theta\right]$ |  |
Summary of helical angle coordinate’s relationship to equilibrium-based frames
$\mathbold{B}$ |  $=\mathbold{\nabla}\alpha\times\mathbold{\nabla}\psi$ |  |
|  $=F\mathbold{\nabla}_{t}\alpha+d\mathbold{\nabla}_{t}\Theta$ |  |
$\mathbold{C}$ |  $=\frac{-1}{\mid\mathbold{\nabla}\psi\mid^{2}}\left[F\mathbold{\nabla}\alpha\times\mathbold{\nabla}\psi-d\mathbold{\nabla}\Theta\times\mathbold{\nabla}\psi\right]$ |  $d\equiv{\mathcal{J}}B^{2}$; |  $D\left(\psi\right)=\left<d\right>_{\Theta}=\frac{V^{\prime}}{4\pi^{2}}\left<B^{2}\right>$
|  $=\mathbold{\nabla}_{t}\alpha$ | Boozer: $d=D$ |  $\Rightarrow\mid\mathbold{\nabla}\psi\mid^{2}\mathbold{C}$ is straight
$\mathbold{J}$ |  $=-F^{\prime}\mathbold{\nabla}\alpha\times\mathbold{\nabla}\psi+P^{\prime}{\mathcal{J}}\mathbold{\nabla}\Theta\times\mathbold{\nabla}\psi$ |  |
|  $=\left(\sigma F-\frac{P^{\prime}}{C^{2}}\right)\mathbold{\nabla}_{t}\alpha+\sigma d\mathbold{\nabla}_{t}\Theta$ | Hamada: ${\mathcal{J}}={\mathcal{J}}(\psi)$ |  $\Rightarrow\mathbold{J}$ is straight
$\mathbold{K}$ |  $=\frac{-1}{P^{\prime}\mid\mathbold{\nabla}\psi\mid^{2}}\left[\left(\sigma F-\frac{P^{\prime}}{C^{2}}\right)\mathbold{\nabla}\alpha\times\mathbold{\nabla}\psi\right.$ |  |
|  $\ \ \ \ \left.+\left(\sigma d\right)\mathbold{\nabla}\Theta\times\mathbold{\nabla}\psi\right]$ |  |
|  $=\frac{-1}{P^{\prime}}\left[F^{\prime}\mathbold{\nabla}_{t}\alpha-P^{\prime}{\mathcal{J}}\mathbold{\nabla}_{t}\Theta\right]$ |  |

Brief summary of key geometric quantities
---
$\mathcal{S}\equiv-\frac{1}{C^{2}}\mathbold{C}\cdot\mathbold{\nabla}\times\mathbold{C}$ | $\mathbold{\nabla}\times\mathbold{C}=-C^{2}\left[F\mathcal{S}\frac{\mathbold{B}}{B^{2}}+\mathcal{S}\frac{\mathbold{C}}{C^{2}}\right]$
$\tau_{n}\equiv\frac{1}{B^{2}}\mathbold{B}\cdot\mathbold{\nabla}\mathbold{C}\cdot\mathbold{\nabla}\psi$ | $2\tau_{n}=\sigma+\mathcal{S}.$
$\mathbold{\kappa}\equiv\left(\hat{\mathbold{b}}\cdot\mathbold{\nabla}\right)\hat{\mathbold{b}}$ | $\mathbold{\kappa}=\frac{1}{B^{2}}\mathbold{\nabla}_{\perp}\left(P+\frac{B^{2}}{2}\right).$
$\mathbold{C}\cdot\mathbold{\nabla}\times\mathbold{K}=-2\kappa_{\psi}+\frac{P^{\prime}}{B^{2}}-\frac{\sigma C^{2}}{P^{\prime}}\mathcal{S}$ |

Table 1: Summary of the flux and helical-angle coordinate systems with
demonstration of the special cases.

##  V Brief Example of equilibrium-based coordinate calculations: Ideal MHD

Before discussing the derivation of the inner layer equations, we give a brief
example of how the new formalism simplifies computations that arise in MHD.
The energy principle is given by Bernstein et al. (1958); Freidberg (1987):

$W_{F}=\frac{\mu_{0}}{2}\int
dV\left[\left[\frac{\vec{Q}}{\mu_{0}}+\left(\mathbold{J}\times\hat{\mathbold{n}}\right)\left(\mathbold{\xi}\cdot\hat{\mathbold{n}}\right)\right]^{2}+\Gamma
P\left(\mathbold{\nabla}\cdot\mathbold{\xi}\right)^{2}-2\left[\mathbold{B}\cdot\mathbold{\nabla}\hat{\mathbold{n}}\cdot\mathbold{J}\times\hat{\mathbold{n}}\right]\left(\mathbold{\xi}\cdot\hat{\mathbold{n}}\right)^{2}\right]$
(79)

Expressing this in terms of our equilibrium-based coordinate systems shows
that the vector $\mathbold{K}$, which has not been discussed before, features
prominently:

$W_{F}=\frac{\mu_{0}}{2}\int
dV\left[\left[\frac{\vec{Q}}{\mu_{0}}+\frac{\xi^{\psi}}{P^{\prime}}\mathbold{K}\right]^{2}+\Gamma
P\left(\mathbold{\nabla}\cdot\mathbold{\xi}\right)^{2}+\left[2P^{\prime}B\mid\mathbold{\nabla}\psi\mid\hat{\mathbold{b}}\cdot\mathbold{\nabla}\hat{\mathbold{n}}\cdot\mathbold{K}\right]\left(\xi^{\psi}\right)^{2}\right]$
(80)

The fact that $\mathbold{K}$ appears in every both the non-compressive
stabilizing term, and the de-stabilizing term indicates that this vector is
important in MHD studies.

The factor in the brackets of the second term involving $\mathbold{K}$ can be
simplified using Eqs. (13) and (IV):

$\displaystyle
2P^{\prime}B\mid\mathbold{\nabla}\psi\mid\hat{\mathbold{b}}\cdot\mathbold{\nabla}\hat{\mathbold{n}}\cdot\mathbold{K}$
$\displaystyle=2\mid\mathbold{\nabla}\psi\mid^{2}C\hat{\mathbold{b}}\cdot\mathbold{\nabla}\hat{\mathbold{n}}\cdot\left[\frac{\sigma}{P^{\prime}}\mathbold{C}+\frac{\mathbold{B}}{B^{2}}\right]$
$\displaystyle=-2P^{\prime}B^{2}\left[\frac{\sigma}{P^{\prime}}\tau_{n}+\frac{1}{BC}\kappa_{n}\right]$
$\displaystyle=-2P^{\prime}\kappa_{n}\mid\mathbold{\nabla}\psi\mid-2\tau_{n}\sigma
B^{2}$

which gives a nice expression for the destabilizing term in terms of geometric
quantities. Using Eq. (68), the torsion can be related to the parallel current
and the shear as well. These give two different forms of $W_{F}$:

$\displaystyle W_{F}$ $\displaystyle=\frac{\mu_{0}}{2}\int
dV\left[\left[\frac{\vec{Q}}{\mu_{0}}+\frac{\xi^{\psi}}{P^{\prime}}\mathbold{K}\right]^{2}+\Gamma
P\left(\mathbold{\nabla}\cdot\mathbold{\xi}\right)^{2}-\left[2P^{\prime}\kappa_{n}\mid\mathbold{\nabla}\psi\mid+2\tau_{n}\sigma
B^{2}\right]\right]$ (81) $\displaystyle=\frac{\mu_{0}}{2}\int
dV\left[\left[\frac{\vec{Q}}{\mu_{0}}+\frac{\xi^{\psi}}{P^{\prime}}\mathbold{K}\right]^{2}+\Gamma
P\left(\mathbold{\nabla}\cdot\mathbold{\xi}\right)^{2}-\left[2P^{\prime}\kappa_{n}\mid\mathbold{\nabla}\psi\mid+\sigma
B^{2}\mathcal{S}-\sigma^{2}B^{2}\right]\right]$ (82)

The first form was first given in DMS Dewar et al. (1984). Both forms are
presented in Greene Greene (1996) along with a discussion of their advantages.
One can see that $\mathbold{K}$ arises naturally in MHD although it has not
formally been discussed in that way. By making its use explicit and clarifying
its relationships, derivations involving its use are easier.

##  VI Derivation of Inner Layer Equations

Although fusion plasmas are very nearly ideal, resistivity can still change
the topology near surfaces where the magnetic winding number is rational. Just
as the flow of wind over an airplane wing, which may be analyzed using
boundary layer theory, resistive modes can be analyzed using boundary layer
theory when the boundary layers occur near rational surfaces. In this theory,
the plasma is analyzed in two regions: an “outer region" where the plasma is
ideal, and an “inner layer" where dissipation is important.

In this section, linear equations are found for the inner resistive layer
based on the narrow-layer-width approximation. The inner-layer equations in
cylindrical geometry were originally derived by Coppi, Greene, and Johnson
(CGJ). In toroidal geometry, the derivation is much more difficult because the
$1/R$ dependence of the magnetic field causes the poloidal harmonics to be
coupled. The formalism used in references Johnson and Greene (1967); Glasser
et al. (1975) is used in deriving these equations in toroidal geometry, which
are presented in the Glasser-Greene-Johnson (GGJ) paper Glasser et al. (1975).
Here we show how the knowledge of the coordinate systems allows one to more
easily derive these relations.

###  VI.1 Resistive instability ordering

We will consider a narrow layer width at the $q_{s}=M/N$ rational surface
where we adopt the usual resistive layer ordering Johnson and Greene (1967);
Glasser et al. (1975)

$x\equiv\psi-\psi_{s}\sim\gamma\sim\eta^{1/3}\sim\epsilon,\ \ \
\frac{\partial}{\partial\psi}\sim\frac{1}{\epsilon}.$ (83)

Equilibrium quantities will be considered to be approximately constant across
the layer.

We begin with the linearized version of the MHD equations. Our notation is
such that equilibrium quantities are shown with capital letters, and perturbed
quantities are denoted with lower case letters. That is, for an arbitrary
quantity $Q$, we write

$Q=Q+\epsilon q_{\left(1\right)}+\epsilon^{2}q_{\left(2\right)}+...$ (84)

The linearized MHD equations are then:

$\displaystyle\gamma^{2}\rho\xi$
$\displaystyle=(\mathbold{\nabla}\times\mathbold{b})\times\mathbold{B}+\mathbold{J}\times\mathbold{b}-\mathbold{\nabla}p$
(85) $\displaystyle\mathbold{b}$
$\displaystyle=\frac{\eta}{\gamma}\mathbold{\nabla}\times\left(\mathbold{\nabla}\times\mathbold{b}\right)+\mathbold{\nabla}\times\left(\xi\times\mathbold{b}\right)$
(86) $\displaystyle p$
$\displaystyle=\xi\cdot\mathbold{\nabla}P+\gamma_{h}P\left(\mathbold{\nabla}\cdot\xi\right)$
(87)

where $\gamma_{h}$ is the ratio of specific heats, $\xi$ is the velocity
divided by complex growth rate $\gamma$.

Because we are considering perturbed quantities near the rational surface, the
$\psi,\Theta,u$ coordinate system is the most useful because it allows easy
identification of resonant perturbations when Fourier expanding the perturbed
quantities:

$\displaystyle f$
$\displaystyle=\sum_{m,n}f_{m,n}\left(\psi\right)e^{i\left(m\Theta-n\zeta\right)}e^{i\gamma
t}$ $\displaystyle=\underbrace{f_{M,N}\left(\psi\right)e^{-iNu}e^{i\gamma
t}}_{resonant\
perturbation}+\underbrace{\sum_{m,n\not=M,N}f_{m,n}\left(\psi\right)e^{-inu}e^{i\gamma
t}e^{i\left(m-nq_{s}\right)\Theta}}_{non-resonant\ perturbation}$
$\displaystyle=\sum_{n}f_{n}\left(\psi,\Theta\right)e^{i\gamma t-iNu}.$

where again $q_{s}=M/N$ is the rational surface of interest. The last line is
the procedure we will adopt for simplicity. If after transforming, the Fourier
transform of a perturbed quantity, $f_{n}$, is independent of $\Theta$, then
only the resonant harmonic is present.

The helical coordinate also facilitates the representation of the parallel
gradient operator:

$\displaystyle\mathbold{B}\cdot\mathbold{\nabla}f$
$\displaystyle={\mathcal{J}}^{-1}\left(\frac{\partial
f}{\partial\Theta}+(q-q_{s})\frac{\partial f}{\partial u}\right)$
$\displaystyle\approx{\mathcal{J}}^{-1}\frac{\partial
f}{\partial\Theta}-\epsilon iN\lambda xf$

where $\lambda={\mathcal{J}}^{-1}q^{\prime}$ and $x=\psi-\psi_{s}$. In the
last line, $q-q_{s}$ was Taylor expanded about the rational surface after
Fourier expanding to give the approximation $q-q_{s}\approx\epsilon
q^{\prime}x$. Other directional derivatives may be easily computed in a
similar manner

$\mathbold{J}\cdot\mathbold{\nabla}f\approx
P^{\prime}iNf-F^{\prime}{\mathcal{J}}^{-1}\frac{\partial
f}{\partial\Theta}-\epsilon F^{\prime}iN\lambda xf$ (88)
$\mid\mathbold{\nabla}\psi\mid^{2}\mathbold{C}\cdot\mathbold{\nabla}f\approx
B^{2}iNf-F{\mathcal{J}}^{-1}\frac{\partial f}{\partial\Theta}-\epsilon
FiN\lambda xf$ (89)

The fact that straight-field line coordinate systems give convenient forms for
$\mathbold{B}\cdot\mathbold{\nabla}$ operators has long been recognized. This
shows that $\mathbold{C}\cdot\mathbold{\nabla}$ and
$\mathbold{J}\cdot\mathbold{\nabla}$ are greatly simplified as well.

As is obvious from the above relations, terms of the form
${\mathcal{J}}^{-1}\partial Q/\partial\Theta$ will often arise. It will be
found useful to eliminate these terms which arise due to non-resonant
harmonics to reduce the problem from a two-dimensional problem in
$\psi,\Theta$ to a one-dimensional problem in $\psi$. To do so, we introduce
the averaging operator:

$\left<Q\right>=\frac{2\pi}{V^{\prime}}\oint Q{\mathcal{J}}d\Theta,$ (90)

If $Q$ is symmetric such that there is no $\zeta$ dependence, this average is
the same as the flux surface average given in Eq. (19).

For the vector quantities, decomposition using the covariant basis vectors is
the most convenient. The perturbed quantities are decomposed and ordered as

$\displaystyle\mathbold{b}$
$\displaystyle=b^{\psi}_{\left(2\right)}\frac{\mathbold{\nabla}\psi}{\mid\mathbold{\nabla}\psi\mid^{2}}+b^{B}_{\left(1\right)}\frac{\mathbold{B}}{B^{2}}+b^{C}_{\left(1\right)}\frac{\mathbold{C}}{C^{2}}$
$\displaystyle=b^{\psi}_{\left(2\right)}\frac{\mathbold{\nabla}\psi}{\mid\mathbold{\nabla}\psi\mid^{2}}+b^{J}_{\left(1\right)}\frac{\mathbold{J}}{J^{2}}+b^{K}_{\left(1\right)}\frac{\mathbold{K}}{K^{2}}$
(91) $\displaystyle\mathbold{\xi}$
$\displaystyle=\xi^{\psi}_{\left(2\right)}\frac{\mathbold{\nabla}\psi}{\mid\mathbold{\nabla}\psi\mid^{2}}+\xi^{B}_{\left(1\right)}\frac{\mathbold{B}}{B^{2}}+\xi^{C}_{\left(1\right)}\frac{\mathbold{C}}{C^{2}}$
$\displaystyle=\xi^{\psi}_{\left(2\right)}\frac{\mathbold{\nabla}\psi}{\mid\mathbold{\nabla}\psi\mid^{2}}+\xi^{J}_{\left(1\right)}\frac{\mathbold{J}}{J^{2}}+\xi^{K}_{\left(1\right)}\frac{\mathbold{K}}{K^{2}}$
(92) $\displaystyle p$ $\displaystyle=p_{\left(1\right)}.$ (93)

The $\psi$ components of the vector quantities are lower order to help satisfy
the divergence criterion as shown below. The ordering also arises from viewing
the vectors as arising from a potential of order $\epsilon^{2}$, and the
non-$\psi$ components are of lower order because of the $d/d\psi$ derivative
which arises in those components Kruger et al. (1999).

Because the perturbed magnetic field is divergence-free
($\mathbold{\nabla}\cdot\mathbold{b}=0$), the components can be related to
each other:

$\displaystyle\mathbold{\nabla}\cdot\mathbold{b}$
$\displaystyle={\mathcal{J}}^{-1}\left[\frac{\partial}{\partial\psi}\left({\mathcal{J}}b^{\psi}\right)+\frac{\partial}{\partial\Theta}\left(\frac{1}{B^{2}}\left(b^{B}-Fb^{C}\right)\right)+\frac{\partial}{\partial
u}\left({\mathcal{J}}b^{C}\right)\right]+\mathcal{O}(\epsilon^{2})$ (94)

Taking the average of the last expression, we have

$\frac{\partial}{\partial\psi}\left<b^{\psi}\right>=iN\left<b^{C}\right>.$
(95)

The orderings shown in Eq. VI.1 are motivated by this balance. Looking at the
pressure equation, Eq. (93), we see that the first two terms are of order
$\epsilon^{2}$ which implies that the the plasma is incompressible to lowest
order, $\mathbold{\nabla}\cdot\mathbold{\xi}\approx 0$. An equation similar to
Eq. 95 for the displacement vector can then be derived.

###  VI.2 First-order equations and magnetosonic waves

Our goal is to obtain equations for averaged perturbed quantities
$\left<b^{\psi}\right>,\left<\xi^{\psi}\right>,\left<b^{B}\right>$ by taking
the appropriate projections of the linearized MHD equations and ordering. The
details of this calculation are presented in the appendix. Here we summarize
the results of the lowest order equations which gives the components whose
resonant harmonics are dominant at the rational surface:

$\displaystyle\mathbold{\nabla}\psi\cdot{\rm(Momentum\ Eq.)}$
$\displaystyle\rightarrow p=-b^{B}$ (96)
$\displaystyle\mathbold{\nabla}\psi\cdot{\rm(Induction\ Eq.)}$
$\displaystyle\rightarrow\left<\xi^{\psi}\right>=\xi^{\psi}$ (97)
$\displaystyle\mathbold{B}\cdot{\rm(Momentum\ Eq.)}$
$\displaystyle\rightarrow\left<b^{B}\right>=b^{B}$ (98)
$\displaystyle\mathbold{J}\cdot{\rm(Momentum\ Eq.)}$
$\displaystyle\rightarrow\left<b^{J}\right>=b^{J}$ (99)
$\displaystyle\mathbold{C}\cdot{\rm(Induction\ Eq.)}$
$\displaystyle\rightarrow\left<\xi^{C}\right>=\xi^{C}$ (100)
$\displaystyle\mathbold{K}\cdot{\rm(Induction\ Eq.)}$
$\displaystyle\rightarrow\left<\xi^{K}\right>=\xi^{K}$ (101)

The first equation is the equilibration of compressional Alfvèn waves to
lowest order. The second, third, and fifth equations are the equilibration of
sound waves to lowest order. The fourth equation is a statement of
$\mathbold{\nabla}\cdot\mathbold{J}=0$ lowest order, and the final equation is
a statement of $\mathbold{\nabla}\cdot\mathbold{V}=0$ to lowest order. Using
Eq. 17 and $\left<b^{J}\right>=b^{J}$, we can find the variation of $b^{C}$
within a surface:

$b^{C}=\frac{C^{2}}{\left<C^{2}\right>}\left<b^{C}\right>+\frac{1}{P^{\prime}}b^{B}\left(\sigma
C^{2}-\left<\sigma C^{2}\right>\frac{C^{2}}{\left<C^{2}\right>}\right).$ (102)

Similarly, using Eq. 17 and $\left<\xi^{K}\right>=\xi^{K}$, we can find the
variation of $\xi^{B}$ within a surface:

$\xi^{B}=\frac{B^{2}}{\left<B^{2}\right>}\left<\xi^{B}\right>-\frac{1}{P^{\prime}}\xi^{C}\left(\sigma
B^{2}-\left<\sigma B^{2}\right>\frac{B^{2}}{\left<B^{2}\right>}\right).$ (103)

One can find similar relations for $b^{K}$ and $\xi^{J}$ using Eqs. (16), 98,
and 100, but they are not needed for this derivation.

###  VI.3 Annihilation operators

We are now ready to look at the next order of the equations. The
$\mathbold{C}$ component of the momentum equation will yield the same
information as the $\mathbold{\nabla}\psi$ component of the magnetic field to
lowest order; i.e., that magnetosonic waves need to be eliminated (Eq.
LABEL:magnetosonic). To obtain additional information from this equation, we
must go to higher order and subtract the two components from each other. A
convenient way to do this is to find an operator that annihilates the lowest
order information, apply it to the full equation before ordering, and then
take the lowest-order terms. This is the concept of annihilation introduced by
Kruskal. CGJ Coppi et al. (1966) and GGJ Glasser et al. (1975) used as the
annihilation operator

$\mathbold{\nabla}\cdot\frac{\mathbold{B}}{B^{2}}\times$ (104)

This operator was so important to MHD theory that the Princeton theory group
dubbed it the Grand Old Operator Coppi . In toroidal geometry and considering
the effects of compressibilityGlasser et al. (1975), this operator leaves a
higher order term, $\xi^{C}_{\left(3\right)}$. To cancel this term, the GGJ
equation is formed by

$\left<\mathbold{\nabla}\cdot\frac{\mathbold{B}}{B^{2}}\times(Momentum\
Eq.)-\frac{\sigma}{P^{\prime}}\frac{\partial}{\partial\psi}\mathbold{B}\cdot(Momentum\
Eq.)\right>$ (105)

Comparing the equation for the complete formation of the annihilated equation
to the relation of $\mathbold{K}$ to $\mathbold{B}$ and $\mathbold{C}$ (Eq.
(13))motivates a new annihilation operator:

$\mathbold{\nabla}\cdot\mathbold{K}\times$ (106)

We study the derivation of the annihilated equation in detail to show its
advantages.

Applying this operator to the momentum equation, we have

$\mathbold{\nabla}\cdot\mathbold{K}\times\left[\gamma^{2}\rho\xi=(\mathbold{\nabla}\times\mathbold{b})\times\mathbold{B}+\mathbold{J}\times\mathbold{b}-\mathbold{\nabla}p\right].$
(107)

The first term is

$\displaystyle\gamma^{2}\rho\mathbold{\nabla}\cdot\mathbold{K}\times\xi$
$\displaystyle=\gamma^{2}\rho\mathbold{\nabla}\cdot\xi^{\psi}\frac{\mathbold{J}}{P^{\prime}\mid\mathbold{\nabla}\psi\mid^{2}}-\xi^{J}\frac{\mathbold{\nabla}\psi}{P^{\prime}\mid\mathbold{\nabla}\psi\mid^{2}}=\frac{-\gamma^{2}\rho}{P^{\prime}}\frac{\partial\xi^{J}}{\partial\psi}+{\mathcal{O}\left(\epsilon^{2}\right)}.$
(108)

The second term is

$\displaystyle\mathbold{\nabla}\cdot\mathbold{K}\times(\mathbold{\nabla}\times\mathbold{b})\times\mathbold{B}$
$\displaystyle=\mathbold{\nabla}\cdot\left(\mathbold{K}\cdot\mathbold{B}(\mathbold{\nabla}\times\mathbold{b})\right)-\mathbold{\nabla}\cdot\left(\mathbold{K}\cdot(\mathbold{\nabla}\times\mathbold{b})\mathbold{B}\right)$
$\displaystyle=\mathbold{B}\cdot\mathbold{\nabla}\left(\mathbold{\nabla}\cdot\mathbold{K}\times\mathbold{b}-\mathbold{b}\cdot\mathbold{\nabla}\times\mathbold{K}\right)$
$\displaystyle=\mathbold{B}\cdot\mathbold{\nabla}\left[\mathbold{\nabla}\cdot\left(b^{\psi}\frac{\mathbold{J}}{P^{\prime}\mid\mathbold{\nabla}\psi\mid^{2}}-b^{J}\frac{\mathbold{\nabla}\psi}{P^{\prime}\mid\mathbold{\nabla}\psi\mid^{2}}\right)-\mathbold{b}\cdot\mathbold{\nabla}\times\mathbold{K}\right]$
$\displaystyle=iN\lambda x\frac{1}{P^{\prime}}\frac{\partial
b^{J}}{\partial\psi}-{\mathcal{J}}^{-1}\frac{\partial}{\partial\Theta}\left(...\right)+{\mathcal{O}\left(\epsilon^{2}\right)}.$

The cancellation used in going from the first step to the second required
$\mathbold{K}\cdot\mathbold{B}=1$. This is the elimination of the magnetosonic
wave required of our annihilation operator.

The third term is

$\displaystyle\mathbold{\nabla}\cdot\mathbold{K}\times\left(\mathbold{J}\times\mathbold{b}\right)$
$\displaystyle=\mathbold{\nabla}\cdot
b^{K}\mathbold{J}=\mathbold{J}\cdot\mathbold{\nabla}b^{K}=-F^{\prime}{\mathcal{J}}^{-1}\frac{\partial
b^{K}}{\partial\Theta}+P^{\prime}iNb^{K}+{\mathcal{O}\left(\epsilon^{2}\right)}.$
(109)

The fourth term is

$\displaystyle\mathbold{\nabla}\cdot\mathbold{K}\times\mathbold{\nabla}p$
$\displaystyle=\mathbold{\nabla}p\cdot\mathbold{\nabla}\times\mathbold{K}=-\frac{\partial
b^{B}}{\partial
u}\mathbold{\nabla}_{t}u\cdot\mathbold{\nabla}\times\mathbold{K}=\mathbold{C}\cdot\mathbold{\nabla}\times\mathbold{K}\
iNb^{B}+{\mathcal{O}\left(\epsilon^{2}\right)}.$ (110)

Putting all of the terms together and taking the average, we have

$\frac{\gamma^{2}\rho}{P^{\prime}}\frac{\partial}{\partial\psi}\left<\xi^{J}\right>+iNP^{\prime}\left<b^{K}\right>+iN\left<\mathbold{C}\cdot\mathbold{\nabla}\times\mathbold{K}\right>b^{B}+iN\Lambda
x\frac{1}{P^{\prime}}\frac{\partial b^{J}}{\partial\psi}=0.$ (111)

This is the annihilated momentum equation. All that remains is to convert the
above variables to the variables we are solving for
$\left<b^{\psi}\right>,\left<\xi^{\psi}\right>,\left<b^{B}\right>$ using (Eqs.
(16)-17) and our previous relationships.

###  VI.4 The inner layer equations

Similar to the annihilated momentum equation, the parallel induction equation
leaves higher order terms in toroidal geometry. The best annhilation operator
for this term is $\mathbold{K}\cdot$ for this equation. Applying it yields

$\gamma
b^{K}=\eta\frac{\partial^{2}b^{K}}{\partial\psi^{2}}\mid\mathbold{\nabla}\psi\mid^{2}-\mathbold{\nabla}\cdot\mathbold{\xi}+\mathbold{B}\cdot\mathbold{\nabla}\xi^{K}+\mathbold{C}\cdot\mathbold{\nabla}\times\mathbold{K}\xi^{\psi}$
(112)

Taking the average and substituting in the expression for
$\mathbold{\nabla}\cdot\mathbold{\xi}$ from the pressure equation yields

$\left<b^{K}\right>=\frac{\eta}{\gamma}\left<\frac{\partial^{2}}{\partial\psi^{2}}b^{K}\mid\mathbold{\nabla}\psi\mid^{2}\right>-iN\Lambda
x\xi^{K}+\frac{b^{B}}{\gamma\gamma_{h}P}+\frac{\xi^{\psi}}{\gamma}\left(\left<\mathbold{C}\cdot\mathbold{\nabla}\times\mathbold{K}\right>+\frac{P^{\prime}}{\gamma_{h}P}\right)$
(113)

Again, we need to convert the variables in the above equation. After a bit of
algebra (see Appendix), the normalized inner layer equations are:

$\displaystyle{\Psi}_{XX}-H{\Upsilon}_{X}=Q\left({\Psi}-\bar{x}{\Xi}\right),$
(114) $\displaystyle{\Gamma}_{X}=H{\Psi}_{XX}+F{\Upsilon}_{X},$ (115)
$\displaystyle
Q^{2}{\Xi}_{XX}-\bar{x}^{2}Q{\Xi}+E{\Upsilon}+Q\bar{x}{\Psi}+{\Gamma}=0,$
(116)
$\displaystyle\frac{1}{Q^{2}}{\Upsilon}_{XX}-\frac{\bar{x}^{2}}{Q^{2}}{\Upsilon}+G{\Upsilon}+\frac{\bar{x}}{Q^{2}}{\Psi}-K{\Gamma}+(G-KE){\Xi}=0.$
(117)

where

$\displaystyle
E=\frac{-\left<C^{2}\right>P^{\prime}}{\Lambda^{2}}\left[\left<\mathbold{C}\cdot\mathbold{\nabla}\times\mathbold{K}\right>+\frac{\Lambda}{P^{\prime}}\frac{\left<\sigma
B^{2}\right>}{\left<B^{2}\right>}\right]$ (118) $\displaystyle
F=\frac{\left<C^{2}\right>}{\Lambda^{2}}\left[P^{\prime
2}\left<\frac{1}{B^{2}}\right>+\left<\sigma^{2}C^{2}\right>-\frac{\left<\sigma
C^{2}\right>^{2}}{\left<C^{2}\right>}\right]$ $\displaystyle
H=\frac{\left<C^{2}\right>}{\Lambda}\left[\frac{\left<\sigma
B^{2}\right>}{\left<B^{2}\right>}-\frac{\left<\sigma
C^{2}\right>}{\left<C^{2}\right>}\right]$ $\displaystyle
G=\frac{\left<B^{2}\right>}{\gamma_{h}P{M}};\ \ \ \ \ \ \ \ \
K_{GGJ}=\frac{\Lambda^{2}R}{P^{\prime 2}{M}}$ $\displaystyle
M=\frac{\left<C^{2}\right>}{P^{\prime 2}}\left[P^{\prime
2}\left<\frac{1}{C^{2}}\right>+\left<\sigma^{2}B^{2}\right>-\frac{\left<\sigma
B^{2}\right>^{2}}{\left<B^{2}\right>}\right]$

Eqs. (116) and (117) are the normalized annihilated momentum and induction
equations. As can be seen by the complicated factors in $F$ and $H$, many of
the complications appearing in going from Eqs. (111) and (113) to these
equations are due to converting the $J$ and $K$ contravariant components which
arise naturally into the normal and parallel components. Examination of Eqs.
(16) and (17) shows easily how many of the factors arise.

In GGJ, it shown that a plasma is ideal unstable if $D_{I}>0$ where
$D_{I}=E+F+H$ and that it is unstable to resistive interchange modes if
$D_{R}=E+F+H>0$. From the expression for $E$ given above it is not obvious
that our form reduces to forms given previously in the literature. We now show
that it does.

Using the equation for $\mathbold{C}\cdot\mathbold{\nabla}\times\mathbold{K}$
(Eq. (78)) and the above relation, one can write

$\mathbold{C}\cdot\mathbold{\nabla}\times\mathbold{K}={\mathcal{J}}^{-1}\frac{\partial}{\partial\psi}{\mathcal{J}}+\frac{F^{\prime}}{P^{\prime}}C^{2}\mathcal{S}+{\mathcal{J}}^{-1}\frac{\partial}{\partial\Theta}\left({\mathcal{J}}\frac{g^{\psi\Theta}}{g^{\psi\psi}}\right).$
(119)

Taking the average of this equation and substituting into the equation for
$E$, yields

$E=-\frac{\left<C^{2}\right>P^{\prime}}{\Lambda^{2}}\left[\frac{V^{\prime\prime}}{V^{\prime}}+\frac{4\pi^{2}F^{\prime}q^{\prime}}{V^{\prime}P^{\prime}}+\frac{\Lambda}{P^{\prime}}\frac{\left<\sigma
B^{2}\right>}{\left<B^{2}\right>}\right]$ (120)

which gives the more familiar $V^{\prime\prime}$ criterion of GGJ Glasser et
al. (1975); Greene (1997). The direct calculation of $E$ is easier in Hamada
coordinates because in calculating the curl of $\mathbold{K}$, it is easier
when $\mathbold{J}$ is straight and $\mathbold{K}$ has a nice contravariant
form seen in Eq. (41).

##  VII Discussion and Summary

Two equilibrium-based local frame of references have been introduced. The $BC$
reference frame is widely used even if not explicitly stated as such. When
normalized, it is closely related to the explicitly named Stix frame
(differing only in that the parallel direction is last and the binormal
changed accordingly) Wright (1998), and to the Frenet-like frame discussed by
Hegna Hegna (2000). Decomposition of fields into the $BC$ local frame was
explicitly performed by MHD theory in the 60’s Coppi et al. (1966). The
usefulness of having $\mathbold{C}$ straight was pointed out by Boozer in the
early 80’s Boozer (1981, 1982).

It was DMS Dewar et al. (1984) that made the transformation between flux
coordinates and the local $BC$ frame explicit. This systematic investigation
is useful for the algebraic expressions that consistently arise in analytic
derivations. The original motivation of the DMS paper is in discussing
numerical implementations. Because different analytic forms of the linearized
MHD equations can have different convergence properties, being able to
analytically transform easily is important.

An equivalent systematic investigation of the transformation properties of the
$JK$ frame has never been performed, and is new to this paper. The fusion
theory community has developed a battery of techniques, both analytic and
computational, for dealing with the stiff time scales of plasma oscillations,
compressional Alfvèn waves, and (not discussed in this paper) cyclotron
motion. One of the main methods for dealing with compressional Alfvèn waves
analytically, the Grand Old Operator, actually does not quite work in toroidal
geometry. It is seen that $\mathbold{\nabla}\cdot\mathbold{K}\times$ is the
appropriate annihilation operator, and that the dominant term of $D_{I}$
appears naturally. The other terms in $D_{I}$ in this derivation are the
metric elements of the $BC$ and $JK$ frames that arise when converting the
annihilated equation variables to the fundamental variables that we wish to
solve for. By having a clean separation from the two frames, the derivation is
clearer. The magnetic well term, as embodied in the $V^{\prime\prime}$
criterion, is shown to be related to the magnetic binormal direction of the
curl of the current density binormal. The important role of $\mathbold{K}$
perhaps should have been obvious sooner given that the ideal MHD potential
energy term, $W_{F}$, is nicely expressed in terms of $\mathbold{K}$ and the
manipulations thereof.

The presentation here focused on the calculation of $W_{F}$ and the inner
layer equations, but the formulas and techniques presented here should be
useful for any theoretical analysis in toroidal geometry where one wishes to
separate the length scales associated with different directions.

###### Acknowledgements.

This manuscript would not exist without the ideas and insights of Dr. John M.
Greene. Dr. Greene asked the primary author for aid in expanding these ideas
after helping him with his thesis work. It is unfortunate that the fulfillment
of these ideas occurred after illness and subsequent death prevented him from
providing further inspiration. The primary author would also like to thank
Drs. Chris Hegna, James Callen, and Alan Turnbull for useful discussions, and
Dr. Tom Jenkins for reviewing the manuscript. This material is based on work
supported by US Department of Energy, Office of Science, Office of Fusion
Energy Sciences under award numbers DE-SC0019067.

## References

  * Glasser et al. (1975) A. Glasser, J. Greene, and J. Johnson, Phys. Fluids 18, 875 (1975).
  * Braginskii (1966) S. Braginskii, _Transport Processes in a Plasma_ (Consultants Bureau, New York, 1966), classic on derivation of two-fluid equations, etc.
  * Greene and Johnson (1962) J. Greene and J. Johnson, Phys. Fluids 5, 510 (1962), remarkable paper that derive Mercier criterion using Hamada coordinates. It made both acceptable but they didn’t get their name on either.
  * Coppi et al. (1966) B. Coppi, J. Greene, and J. Johnson, Nucl. Fusion 6, 101 (1966).
  * Dewar et al. (1984) R. Dewar, D. Monticello, and W.-C. Sy, Phys. Fluids 27, 1723 (1984).
  * Stix (1992) T. Stix, _Waves in Plasmas_ (American Institute of Physics, New York, NY, 1992).
  * Wright (1998) J. C. Wright, Ph.D. thesis, Princeton University (1998).
  * Hazeltine and Meiss (1992) R. Hazeltine and J. Meiss, _Plasma Confinement_ (Addison-Wesley Publishing Company, Redwood City, CA, 1992).
  * Johnson and Greene (1967) J. Johnson and J. Greene, Plasma Physics 9, 611 (1967).
  * Boozer (1981) A. Boozer, Phys. Fluids 24, 1999 (1981).
  * Kruger et al. (1999) S. Kruger, C. Hegna, and J. Callen, Phys. Plasmas 6 (1999).
  * Boozer (1982) A. Boozer, Phys. Fluids 27, 2441 (1982).
  * Freidberg (1987) J. Freidberg, _Ideal Magnetohydrodynamics_ (Plenum Publishing Company, New York, NY, 1987).
  * Grimm et al. (1976) R. Grimm, J. Greene, and J. Johnson, J. Comp. Phys. p. 253 (1976), the PEST paper. Very good.
  * Hegna (2000) C. Hegna, Phys. Plasmas 7, 3921 (2000).
  * Greene and Chance (1981) J. Greene and M. Chance, Nucl. Fusion 21, 453 (1981).
  * Greene (1996) J. Greene, Phys. Plasmas 3, 8 (1996).
  * Bernstein et al. (1958) I. Bernstein, E. Freiman, M. Kruskal, and R. Kulsrud, Proc. R. Soc. London Ser. A 244, 17 (1958).
  * (19) B. Coppi, Private communication (????).
  * Greene (1997) J. Greene, Comments on Modern Physics: Part E 6, 331 (1997).