Here we will show that to first order the conditions for strong isodrasticity
are precisely the weak isodrasticity conditions, and give a computational test
for weak isodrasticity in terms of “Melnikov functions” that is simpler than
that of Section 2. Finally, we interpret the Melnikov functions in terms of
flux of phase-space volume making transitions.

###  8.1. Strong isodrasticity to first order

To compute $\Xi^{\pm}$ to first order in $\tilde{\mu}$, we can integrate the
rate of change of fieldline label along the unperturbed trajectory of ZGCM,
compensated by the rate of change of the fieldline label at the base point.

Recall that the equations to first order can be written in scaled time $\tau$
as

$\frac{dX}{d\tau}=ub+\frac{\sqrt{\tilde{\mu}}}{|B|}\left(u^{2}c_{\perp}+b\times\nabla|B|\right),$

with

(37) $u=\pm\sqrt{2(h-|B|)}$

and $c=\mbox{curl}\,b$ (as already mentioned, $c_{\perp}$ can be written as
$b\times\kappa$ with $\kappa=b\cdot\nabla b$).

Let $\xi$ be a coordinate on $\Sigma^{-}$ (e.g. $x$ or $y$ from the previous
section) and extend it along the fieldlines to a fieldline label. Its
derivative $a=d\xi$ (a covector) can be computed by setting $a(0)=d\xi$ on
tangents to $\Sigma^{-}$ and $a(0)b=0$, and integrating the adjoint equation
along the fieldlines:

(38)
$\frac{da_{i}}{ds}=-a_{j}(s)\partial_{i}b^{j}(x(s)),\quad\frac{dx^{i}}{ds}=b^{i}(x(s)),$

where $s$ is arclength. The covector $a$ quantifies the linearised change in
the fieldline label. For trajectories starting on $N^{-}$, take $h$ in (37) to
be $|B|$ on $\Sigma^{-}$ to first order. Then the first-order change
$\Delta\xi$ in fieldline label $\xi$ along $W^{-}$ is given by integrating the
rate of change of $\xi$ by the perturbed vector field along the unperturbed
trajectory, while simultaneously subtracting off its rate of change at the
base point.

Thus

$\Delta\xi=\sqrt{\tilde{\mu}}\int\frac{u}{|B|}a\,c_{\perp}+\frac{1}{u}\left(a\frac{b\times\nabla|B|}{|B|}-C\right)\
ds,$

where $C=a\frac{b\times\nabla|B|}{|B|}$ at the initial point on $\Sigma^{-}$.
The subtraction of $C$ compensates for the denominator $u$ at the start of the
integration, where $u\sim\sqrt{-|B|^{\prime\prime}}\,s$ for small $s$. The
denominator does not give a problem at a generic bounce point because it
behaves like $\sqrt{s_{0}-s}$ there, which gives an integrable singularity; a
tidy way to compute the integral accurately is to switch variable of
integration from $s$ to $u$ near the bounce.

The first-order change along $W^{+}$ is the negative of $\Delta\xi$. So the
displacement in the chosen fieldline label $\xi$ from $W^{+}$ to $W^{-}$ is to
first order

$2\Delta\xi=2\sqrt{\tilde{\mu}}\mathcal{M}_{\xi},$

with

$\mathcal{M}_{\xi}=\int\frac{u}{|B|}a\,c_{\perp}+\frac{1}{u}\left(a\frac{b\times\nabla|B|}{|B|}-C\right)\
ds.$

It is a function of the initial point on $\Sigma^{-}$ and we call it the
Melnikov function for the fieldline label $\xi$, by analogy with formulae for
splitting of separatrices in other contexts.

If one takes two independent fieldline labels $x,y$ then the Melnikov function
has two components $\mathcal{M}_{x},\mathcal{M}_{y}$. This is particularly
important when one takes initial points on or near a critical point of $|B|$
on $\Sigma^{-}$. Compare the discussion of the mirror machine in Section 6.1.

For a periodic orbit on $N^{-}$, however, we are more interested in the
displacement in fieldline label $h=|B|$ where it crosses $\Sigma^{-}$ than the
possible phase shift along the periodic orbit. Thus we specialise to take $h$
as coordinate on $\Sigma^{-}$ and compute the resulting first-order change in
$h$. In this case several things simplify. Firstly, we can take $a(0)=d|B|$
because $h=|B|$ on $\Sigma^{-}$ and $i_{b}d|B|=0$ on $\Sigma^{-}$. Secondly,
$C=0$. Then letting

$k=\mbox{curl}(ub)$

we obtain

$\mathcal{M}_{h}=\int\frac{ak_{\perp}}{|B|}\ ds$

and $ab=0$ so we can drop the $\perp$. Thus we obtain the simple formula

(39) $\mathcal{M}_{h}=\int\frac{ak}{|B|}\,ds.$

Of course this hides the computation of $k$ and also the fact that $k$ is
singular at $s=0$ (but the singularity is annihilated by $a$).

Note that in practice it may be better to integrate with respect to time $T$
for fieldline flow $dx/dT=B(x)$ than arclength $s$. This would replace
$ds/|B|$ by $dT$, and also simplify the computation of (38) as indicated in
Appendix I.

###  8.2. Relation to weak isodrasticity

Next, we relate the Melnikov function $\mathcal{M}_{h}$ to the ratio of
$dh\wedge d\text{\j}$ to $\beta$ on $\Sigma^{-}$.

######  Theorem 5.

$dh\wedge d\text{\j}=\mathcal{M}_{h}\beta.$

Thus the main part of the weak isodrastic condition, that $d\text{\j}$ and
$dh$ be linearly independent on $\Sigma^{-}$, is the first-order condition for
strong isodrasticity, and $\mathcal{M}_{h}$ is the function $\mathcal{M}$ of
Section 2.

###### Proof.

Recall that $\text{\j}=\int u\,ds=\int ub^{\flat}$ along the segment of
fieldline from $\Sigma^{-}$ to the first bounce point. Hence for a tangent
vector $v$ to $\Sigma^{-}$, $d\text{\j}(v)=\int i_{\tilde{v}}d(ub^{\flat})$,
where $\tilde{v}$ is obtained from $v$ by flowing with the derivative of
$\dot{x}=b(x)$ ($\dot{\tilde{v}}^{i}=\tilde{v}^{j}\partial_{j}b^{i}$). But
$d(ub^{\flat})=i_{k}\Omega$ and displacement along the segment is given by
$b\,ds$, so

$d\text{\j}(v)=\int i_{b}i_{\tilde{v}}i_{k}\Omega\,ds.$

Now $\Omega=\beta\wedge b^{\flat}/|B|$, and $i_{b}\beta=0,i_{b}b^{\flat}=1$,
so this reduces to

$d\text{\j}(v)=\int i_{\tilde{v}}i_{k}\beta\frac{ds}{|B|}.$

Near a point of $\Sigma^{-}$ with $dh\neq 0$ we can choose $h$ as one
fieldline label and can take another one $g$ such that $\beta=dh\wedge dg$ (by
Darboux’s theorem; this is a construction of Clebsch coordinates). Extend $h$
and $g$ to fieldline labels along the fieldlines. Let us take $v$ to have
$dh\,v=0$. Then

$d\text{\j}(v)=\int i_{k}dh\,i_{\tilde{v}}dg\frac{ds}{|B|}=\int
i_{k}dh\frac{ds}{|B|}i_{v}dg=\mathcal{M}_{h}\,i_{v}dg,$

as $i_{\tilde{v}}dg$ is constant along the fieldline. Then evaluate $dh\wedge
d\text{\j}(w,v)$ on $\Sigma^{-}$ for an arbitrary $w$. It is
$dh(w)d\text{\j}(v)=dh(w)\mathcal{M}_{h}\,i_{v}dg=\mathcal{M}_{h}\,\beta(w,v)$.
Since $\Sigma^{-}$ is only two-dimensional, it suffices to evaluate 2-forms on
any pair of independent vectors to determine them. Hence $dh\wedge
d\text{\j}=\mathcal{M}_{h}\,\beta$, i.e. $\mathcal{M}_{h}$ is the ratio of
$dh\wedge d\text{\j}$ to $\beta$ on $\Sigma^{-}$.

At a point of $\Sigma^{-}$ where $dh=0$ then both $\mathcal{M}_{h}$ and
$dh\wedge d\text{\j}$ are zero, so the same result holds. ∎

Computation of the Melnikov function via (39) (rather than by numerical
differentiation of $h$ and ȷ) for some examples will be reported in a future
paper.

To complete the discussion of strong isodrasticity at first order, we look at
$\Sigma^{0}$. The strong isodrasticity condition that
$H=\tfrac{1}{2}u^{2}+|B|$ is constant along components of $N^{0}$, evaluated
to first order is equivalent to $h$ being constant on components of
$\Sigma^{0}$. Linear dependence of $dh$ and $d\text{\j}$ on $N^{-}$ extends to
the boundary, so from $h$ constant on boundary components we also deduce that
ȷ is constant on them.

Hence strong isodrasticity at first order is equivalent to weak isodrasticity,
as claimed.

###  8.3. Interpretation as flux

To conclude this section, we interpret the Melnikov function in terms of the
first-order flux of phase-space volume making the given transition. This is
analogous to the interpretation of $dh\wedge d\j$ as transition flux for
reduced dynamics in Section 2.

For a 2 DoF Hamiltonian system with Hamiltonian $H$, symplectic form $\omega$
and vector field $V$ defined by $i_{V}\omega=dH$, the phase-space volume form
is $\tfrac{1}{2}\omega\wedge\omega$ and the energy-surface volume form
$\varepsilon$ on $H^{-1}(E)$ is defined so that
$\tfrac{1}{2}\omega\wedge\omega=dH\wedge\varepsilon$. A standard calculation
(e.g. [MM, M90]) shows that the energy-surface volume flux form
$\phi=i_{V}\varepsilon$ is just $\omega$.

Specialising to the mirror machine example for illustration, it follows that
the flux of energy-surface volume making the transition from free to bouncing
is the integral of $\omega$ over the lobe between $W^{\pm}$ on the bounce
surface $B=\\{v=0\\}$ in Figure 19. In the scaling we are using,
$\omega=\beta/\sqrt{\tilde{\mu}}+d(ub^{\flat})$, so is just
$\beta/\sqrt{\tilde{\mu}}$ on $\\{v=0\\}$. Choose a local coordinate $g$ on
$\Sigma^{-}$ such that $\beta=dh\wedge dg$ (possible where $dh\neq 0$ by
Darboux’s theorem) and flow the functions $h$ and $g$ along the field. Then
the flux is $\int dh\wedge dg/\sqrt{\tilde{\mu}}$ over the lobe. But the
change in $h$ from $W^{+}$ to $W^{-}$ is to first order
$2\sqrt{\tilde{\mu}}\mathcal{M}_{h}$ so the flux is $\int 2\mathcal{M}_{h}dg$
along the arc of the level set of $h$ between zeroes of $\mathcal{M}_{h}$.

We can treat a range of energies simultaneously. The flux-form for phase space
volume is $i_{V}(\tfrac{1}{2}\omega\wedge\omega)=dH\wedge\omega$. So, given an
area $A$ on $\Sigma^{-}$ corresponding to a transition, the flux of phase-
space volume making the transition is the integral of $dH\wedge\omega$ across
the corresponding region on $\\{v=0\\}$. To first order, $dH=dh$, because we
are treating trajectories that graze the maximum in $|B|$. So the flux is
$\int_{A}2\mathcal{M}_{h}dh\wedge dg=\int_{A}2\mathcal{M}_{h}\beta$. The
factor 2 comes from the convention that $L$ is computed from one zero of $v$
to the other, whereas for a full period it would be twice as much.

So far, we have been treating the problem in scaled variables. To turn this
from scaled variables to the original variables, we use that the real
symplectic form is $\sqrt{m\mu}$ times the scaled one (29), so real phase-
space volume is $m\mu$ times the scaled one, and real time is $\sqrt{m/\mu}$
times scaled time. So real phase-space flux across $A$ is:

$2m^{1/2}\mu^{3/2}\int_{A}\mathcal{M}_{h}\beta.$

A point with given $h$ on $\Sigma^{-}$ corresponds to real energy $\mu h$.

To convert these results to the flux of particles, we need to introduce a
distribution function $\rho$, giving their density with respect to phase-space
volume (this use of $\rho$ is distinct from that for the gyroradius). It is
best to include the dependence of the density on $\mu$, rather than treating a
fixed $\mu$. The phase-space volume form $dq_{1}\wedge dq_{2}\wedge
dq_{3}\wedge dp_{1}\wedge dp_{2}\wedge dp_{3}$ for the 3DoF problem converts
in gyro-coordinates to $m^{2}\tilde{B}_{\parallel}\,\Omega\wedge
dv_{\parallel}\wedge d\mu\wedge d\phi$, with $\Omega$ being ordinary volume
for guiding-centre position and $\phi$ being gyrophase. Integrating over
gyrophase we obtain the gyro-averaged volume-form $2\pi
m^{2}\tilde{B}_{\parallel}\,\Omega\wedge dv_{\parallel}\wedge d\mu$. This is
preserved by the Hamiltonian guiding-centre flow (though not in general by
(9), even if $\tilde{B}_{\parallel}$ is replaced by $|B|$, another reason to
prefer the Hamiltonian version). So the number of particles of given type in a
volume $W$ of gyro-averaged phase-space can be written as
$\int_{W}\rho(X,v_{\parallel},\mu)\,2\pi
m^{2}\tilde{B}_{\parallel}(X)\,\Omega\wedge dv_{\parallel}\wedge d\mu,$ where
$\rho$ denotes the usual scalar guiding-center distribution function.

Note that the gyro-averaged volume-form is the wedge of $2\pi\tfrac{m}{e}d\mu$
with the guiding-centre volume-form. The prefactor is because the standard
convention for action variables was broken by many plasma physicists (though
tends to be respected by A.J. Brizard and high-energy particle physicists). To
see this factorisation of the gyro-averaged volume-form, reall from (20) that
the guiding-centre volume-form $\Lambda=em\tilde{B}_{\parallel}\,\Omega\wedge
dv_{\parallel}$.

Thus to obtain the flux of particles from the flux of guiding-centre volume,
one has to multiply by $2\pi\tfrac{m}{e}\rho\,d\mu$ and integrate over $\mu$.
To leading order in $\mu$ we can replace $\tilde{B}_{\parallel}$ by $|B|$. So
the flux of particles corresponding to an area $A$ on $\Sigma^{-}$, to leading
order in $\mu$, is

$\int_{A}\int_{\mathbb{R}_{+}}\tfrac{4\pi}{e}{(m\mu)^{3/2}}\rho(X,0,\mu)\,d\mu\,\mathcal{M}_{h}\beta.$

A point to note is that the picture to first order depends on $E$ and $\mu$
through only the combination $E/\mu$. The rate of transition has a prefactor
$\sqrt{\mu}$ (or other powers depending on scaling), but is otherwise
independent of $\mu$ given $E/\mu$. As pointed out by Roscoe White (personal
communication), this gives hope that fusion $\alpha$-particles might remain in
their initial class despite slowing down, because it seems for fusion
$\alpha$-particles $E$ and $\mu$ decrease in such a way that $E/\mu$ remains
roughly constant.

##  9\. Discussion

We have generalised the concept of omnigenous magnetic fields and their
analysis by [CS], to remove the requirement of a flux function (and even with
a flux function, our condition is much weaker), obtaining the concept of
weakly isodrastic field. They are magnetic fields for which there are no
transitions between classes of guiding-centre motion, under the assumption
that the longitudinal invariant is conserved. Furthermore, we have provided a
quantification of deviation from the ideal case, namely Melnikov functions,
thereby making available objective functions for optimisation of design.

We have extended our theory, by removing the assumption of conservation of the
longitudinal invariant, to a notion of strong isodrasticity, which provides an
exact prevention of transitions between classes. We have proved that to first-
order in the magnetic moment, strong isodrasticity is weak isodrasticity, thus
justifying weak isodrasticity as a first-order approximation. We have
illustrated how isodrasticity is lost for general perturbations from
axisymmetric fields for toy mirror machines and a toy tokamak. The exact
theory has the advantage that it can in principle be applied to higher order
guiding-centre approximations, of potential relevance to the fusion alpha-
particles.

We have shown how to construct many weakly isodrastic mirror fields, in
particular that are not omnigenous. The main question that remains is whether
isodrasticity (weak or strong) is possible outside axisymmetry for a
stellarator. Quasisymmetry implies isodrasticity. Although perfect
quasisymmetry perhaps does not exist outside axisymmetry, close to
quasisymmetric fields can be made (for recent examples, see [LP]). Close to
quasisymmetric fields should be close to isodrastic, but perhaps there is a
larger class of isodrastic fields than quasisymmetric. Omnigenous fields are
weak isodrastic. Perhaps truly omnigenous fields have to be quasisymmetric,
but close to omnigenous fields can be made, such as in Wendelstein 7-X. As
weak isodrasticity is weaker than omnigenity, there is hope that truly
isodrastic fields can be made. One option is to make all the marginal cases
heteroclinic, but we have shown (not included in this paper) that this reduces
to omnigenity. So we wish to make examples where most of the marginal cases
are homoclinic, as discussed for typical perturbations of a tokamak in Section
3.2, but there will be codimension-one cases of heteroclinic, leading to the
double transition scenario of Appendix G. We are looking for examples. For
general Hamiltonian systems it is certainly true that there is a larger class
of systems with some perfect separatrices than the integrable class (see
Appendix L). Also, the concept of isodrasticity motivates some clear objective
functions that could be fed to an optimisation to automate a search. For
example, one could use the maximum of the Melnikov function or the integral of
its positive part. Even if one does not make exact isodrasticity, these
objective functions could be weighed against others in the design of
stellarators.

The theory also suggests that one might be able to make controlled transitions
between classes by slightly breaking isodrasticity via trim coils. The
Melnikov function for a normally hyperbolic submanifold specifies the flux of
energy-surface volume making a given transition, so if one could learn the
effects of trim coils on the Melnikov function then one could control the
flux.

Finally, the theory sheds light on the effects of imperfections on tokamaks
and quasi-symmetric stellarators. The marginal bouncing trajectories of the
ideal case give rise to double homoclinic manifolds to hyperbolic periodic
orbits for guiding-centre motion. Breaking the symmetry in general breaks both
these separatrices, leading to a stochastic layer in which trajectories
transition between bouncing and co- and counter-circulating. Furthermore, the
guiding-centre dynamics near the magnetic axis is in general significantly
disturbed from the axisymmetric case, leading to mixing in the core (though
this is not necessarily a bad thing, as discussed by [Bo]).

## Acknowledgements

This work was supported by the Simons Foundation (601970, RSM) under the
“Hidden symmetries and fusion energy” collaboration. We also acknowledge
funding support from the US DOE Office of Advanced Scientific Computing
Research (ASCR), DOE-FOA-2493 “Data-intensive scientific machine learning and
analysis”. We are grateful to Nikos Kallinikos for the etymology and for
making some of the figures, to Elizabeth Paul for trying out our ideas in
various realistic nearly quasi-symmetric fields, and to others in the
collaboration for their comments, notably Matt Landreman, Per Helander, Roscoe
White, Eduardo Rodriguez and Gabriel Plunk. Also we thank Anatole Neishtadt
and John Cary for useful pointers to literature.

##  Appendix A Electrostatic and gravitational fields and relativity

###  A.1. Electrostatic and gravitational fields

To add the effects of an electrostatic or gravitational field on guiding
centre motion, one takes

$H=\tfrac{1}{2}mv_{\parallel}^{2}+\mu|B|+e\Phi+mV,$

where $\Phi$ and $V$ are electrostatic and gravitational potentials,
respectively. The symplectic form is unchanged. The resulting equations of
motion are

$\displaystyle\dot{X}$
$\displaystyle=\left({v_{\parallel}}\widetilde{B}+\frac{\mu}{e}b\times\nabla|B|+b\times\nabla\Phi+\frac{m}{e}b\times\nabla
V\right)/\widetilde{B}_{\parallel}$ $\displaystyle\dot{v}_{\parallel}$
$\displaystyle=-\frac{\widetilde{B}}{\widetilde{B}_{\parallel}}\cdot\left(\frac{\mu}{m}\nabla|B|+\frac{e}{m}\nabla\Phi+\nabla
V\right).$

ZGCM is just 1 DoF motion of unit mass in the potential
$\frac{\mu}{m}|B|+\frac{e}{m}\Phi+V$. The longitudinal invariant is modified
to

$L=\int_{s_{1}}^{s_{2}}\sqrt{2m(E-\mu|B|-{e}\Phi-mV)}\,ds.$

The surface $\Sigma$ is modified to be the set of points at which the first
derivative $\mu|B|^{\prime}+e\Phi^{\prime}+mV^{\prime}=0$, where ′ again
denotes derivative along the magnetic field, and decomposes into
$\Sigma^{\pm}$ and $\Sigma^{0}$ according to the sign of the second
derivative. But $\Sigma$ now depends on the ratio $\mu:e:m$. Weak
isodrasticity becomes that $dE\wedge dL=0$ on $\Sigma^{-}$, where $L$ is for
segments starting on $\Sigma^{-}$. The exact FGCM dynamics is still 2DoF and
has NHS with contracting submanifolds, whose intersections can be analysed the
same way.

###  A.2. Relativity

As explained in an appendix to [BKM], relativity can be incorporated in either
laboratory time or proper time. We treat the former here. The first adiabatic
invariant becomes $\mu=\frac{p_{\perp}^{2}}{2m|B|}$ with $p=\gamma mv$,
$\gamma=(1-|v|^{2}/c^{2})^{-1/2}$. The guiding-centre Hamiltonian becomes

$H=c\sqrt{m^{2}c^{2}+p_{\parallel}^{2}+2m\mu|B|}+e\Phi+mV.$

The second adiabatic invariant becomes

$L=\int p_{\parallel}\,ds=\int\sqrt{(\tfrac{E}{c}-e\Phi-
mV)^{2}-m^{2}c^{2}-2m\mu|B|}\,ds$

at $H=E$. If $\Phi$ and $V$ are constant we again have $\Sigma$ the set of
points where $|B|^{\prime}=0$ and the same analysis of weak isodrasticity.
With $\Phi$ or $V$, however, the shape of the argument of the square root now
depends on $E$ so $\Sigma$ becomes $E$-dependent and the nice picture breaks
down.

Nonetheless the exact FGCM still has a NHS continuing the non-relativistic
$\Sigma$ to the relativistic case, at least for $\mu$ not too large. The same
picture of its contracting manifolds holds. For magnetically confined fusion
devices, probably adiabatic invariance of $\mu$ fails before relativistic
effects come in (the alpha particles produced by DT fusion have speed 4.3% of
the speed of light and in a 1T field have gyroradius 0.266 metres times the
sine of their pitch angle). But perhaps in some astrophysical contexts,
relativistic guiding-centre motion is relevant.

##  Appendix B Fields with a flux function

We take the opportunity to review the theory of magnetic fields possessing a
flux function, including some results that we have not found in the
literature.

Existence of a flux function is automatic for non-degenerate
magnetohydrostatic (MHS) fields, i.e. those satisfying $i_{J}i_{B}\Omega=-dp$
($J\times B=\nabla p$) for some function $p$ with $dp\neq 0$ ($\nabla p\neq
0$) almost everywhere, where $i_{J}\Omega=dB^{\flat}$ ($J=\mbox{curl}\ B$):
one can take $\psi=p$ (though it is generally preferred to take the toroidal
flux enclosed by a level set of $p$). It is also automatic for most
axisymmetric fields: let $u$ be the vector field $\partial_{\phi}$ in
cylindrical coordinates, then $di_{u}i_{B}\Omega=0$ ($\mbox{curl}(B\times
u)=0$) so $i_{u}i_{B}\Omega=d\psi$ ($B\times u=\nabla\psi$) for some function
$\psi$ locally, and generally there is no cohomological obstruction, so then
$\psi$ is a global function, $i_{B}d\psi=0$ and $d\psi\neq 0$ except where $B$
is parallel to $u$.

The components of level sets of $\psi$ are called flux surfaces. They are
called regular if $d\psi\neq 0$ ($\nabla\psi\neq 0$) on them. The bounded
boundaryless regular flux surfaces are co-oriented by $\nabla\psi$ and hence
oriented. Because they carry a nowhere-zero vector field $B$ their Euler
characteristic is zero, so they are tori. Furthermore, they carry an area-form
$\mathcal{A}$ satisfying $\mathcal{A}\wedge d\psi=\Omega$, e.g.

$\mathcal{A}=i_{n}\Omega$

($\mathcal{A}(\xi,\eta)=n\cdot(\xi\times\eta)$) with
$n=\nabla\psi/|\nabla\psi|^{2}$. The fieldline flow preserves $\mathcal{A}$ on
the flux surfaces, because applying $di_{B}$ to $\mathcal{A}\wedge
d\psi=\Omega$ produces $di_{B}\mathcal{A}\wedge d\psi=0$ and so
$di_{B}\mathcal{A}=0$ on tangents to the flux surface. So there is no
asymptotic convergence of fieldlines in either direction of time. This forces
the flow to have a cross-section (a closed curve transverse to the vector
field such that the trajectory of every point crosses it in forward and
backward time) because the only other option for a nowhere-zero vector field
on a 2-torus has a “Reeb component” (an annulus bounded by periodic orbits in
opposite directions) and area can not be preserved. See sec.4 of [BGKM] for a
summary of the theory of vector fields on a 2-torus.

Furthermore, integrating $i_{B}\mathcal{A}$ from a reference point produces a
local coordinate on the flux surface that is preserved by the flow. It follows
that the return map to a cross-section is a rigid translation (in this
coordinate). Thus, the fieldline flow on each torus is equivalent to that of a
(non-zero) constant vector field on a standard torus
$\mathbb{R}^{2}/\mathbb{Z}^{2}$ up to a possible time-change (and examples can
be made for which that is necessary) (equivalently, conjugate to a vector
field with a constant direction). In particular, the fieldlines wind around
the torus with a well-defined winding ratio (limit of ratio of number of turns
in two angles), called its rotational transform $\iota$, and in the rational
case all the fieldlines are closed.

In the case of Diophantine winding ratio, i.e. $|k\iota-m|\leq Ck^{-\sigma}$
for some $C>0$, $\sigma\geq 1$ for all integers $k,m$ with $k>0$, if the flow
is smooth enough then it is conjugate to a constant vector field without any
need for a time-change factor, but the conjugacy is in general less smooth
than the flow (by about $\sigma$ derivatives, depending whether one works in
$C^{r}$ or Sobolev spaces). The $C^{\infty}$ case is treated by [KH] in Prop
2.9.5, where $C^{\infty}$ conjugacy results, but the method of proof there
gives results for less smoothness.

##  Appendix C Omnigenity for passing trajectories

In [LC], omnigenity is shown to imply that $L$ is constant on flux surfaces
for bouncing trajectories of given energy and class (actually, the extension
to more than one class comes in [PCHL]). As mentioned at the end of [BKM],
however, omnigenity also implies that $L$ is constant for circulating
trajectories of given energy on each rational flux surface, and to deduce
omnigenity both conditions must be satisfied. In the meantime, we realised
that constancy of $L$ for circulating trajectories on a rational surface is
implied by constancy for bouncing trajectories.

In this appendix we first prove that omnigenity requires constancy of $L$ for
all circulating trajectories on a rational flux surface. Then we prove that
this is implied by the same statement for bouncing trajectories. After this,
we give an equivalent proof to [He] that the time-average of
$V\cdot\nabla\psi$ is zero for circulating trajectories on irrational tori.

It follows that constancy of $L$ for bouncing trajectories of given energy and
class on each flux surface implies omnigenity, a result that we feel had not
been established correctly before. The continuity argument of [He] from
irrational to rational surfaces shows only that the flux-surface average of
$i_{V}d\psi$ for circulating particles is zero on each rational surface. But
the time-averages along individual fieldlines on a rational surface need not
agree with the flux-surface average.

###  C.1. Omnigenity for circulating trajectories on rational flux surfaces

The longitudinal adiabatic invariant for general periodic orbits of zeroth-
order guiding-centre motion (ZGCM) (not just bouncing ones) is

$L=\int_{\gamma}eA^{\flat}+mv_{\parallel}b^{\flat},$

where $\gamma$ is the segment of fieldline covered. The term in $A^{\flat}$
gives zero when integrated along a bouncing trajectory because it backtracks
exactly, which is why it is usually left out. Also it is a constant for given
flux surface when $\gamma$ is a closed loop restricted to the flux surface, so
plays no role in the present discussion. Thus we take
$L=\int_{\gamma}mv_{\parallel}b^{\flat}$ for both.

For circulating periodic orbits $\gamma$ of ZGCM for a field with a flux
function $\psi$, the first-order drift in $\psi$ averaged over one period $T$
is

$\langle\dot{\psi}\rangle=\frac{1}{T}\int_{\gamma}V\cdot\nabla\psi\
dt=\frac{1}{T}\int_{\gamma}\frac{i_{V}d\psi}{v_{\parallel}}b^{\flat},$

with

(40) $\tfrac{1}{2}mv_{\parallel}^{2}=E-\mu|B|.$

The definition of omnigenity uses the non-Hamiltonian velocity

$V=v_{\parallel}b+\frac{mv_{\parallel}^{2}}{e|B|}c_{\perp}+\frac{\mu
v_{\parallel}}{e|B|}b\times\nabla|B|$

from (9) rather than the Hamiltonian one of (6-8). We drop the
$\dot{v}_{\parallel}$ component because we need only $V\cdot\nabla\psi$. Then
using $B\cdot\nabla\psi=0$, we obtain

$\langle\dot{\psi}\rangle=\frac{1}{T}\int_{\gamma}\frac{1}{e|B|}(mv_{\parallel}c+\frac{\mu}{v_{\parallel}}b\times\nabla|B|)\cdot\nabla\psi\
b^{\flat}=\frac{1}{T}\int_{\gamma}(\tfrac{mv_{\parallel}}{e|B|}c\cdot\nabla\psi-\tfrac{\mu}{ev_{\parallel}}\xi\cdot\nabla|B|)b^{\flat},$

where

$\xi=(b\times\nabla\psi)/|B|.$

To test whether $L$ is constant on a flux surface (for given $E$), it is
enough to test whether $dL\ \xi=0$, because tangents to a flux surface are
linear combinations of the vector fields $B$ and $\xi$, and the change of $L$
along $B$ is zero. Using the Lie derivative $L_{\xi}=i_{\xi}d+di_{\xi}$
(distinguish the function $L$ and the Lie derivative $L_{\xi}$ along $\xi$) on
differential forms,

$dL\
\xi=\int_{\gamma}L_{\xi}(eA^{\flat}+mv_{\parallel}b^{\flat})=\int_{\gamma}ei_{\xi}dA^{\flat}+mi_{\xi}d(v_{\parallel}b^{\flat}),$

the total derivative term integrating to zero because $\gamma$ is closed. Now
$dA^{\flat}=i_{B}\Omega$ so along fieldline $\gamma$ it gives zero. The second
term gives

$dL\
\xi=\int_{\gamma}m(i_{\xi}dv_{\parallel}\,b^{\flat}-v_{\parallel}i_{\xi}db^{\flat}).$

Differentiating (40) we obtain $mv_{\parallel}dv_{\parallel}=-\mu d|B|$; and
using $db^{\flat}=i_{c}\Omega$, we obtain
$i_{\xi}db^{\flat}=-i_{c}i_{\xi}\Omega=i_{c}(b^{\flat}\wedge d\psi)/|B|$. Thus

$dL\ \xi=\int_{\gamma}-\frac{\mu}{v_{\parallel}}i_{\xi}d|B|\
b^{\flat}-\frac{mv_{\parallel}}{|B|}i_{c}b^{\flat}d\psi+\frac{mv_{\parallel}}{|B|}b^{\flat}i_{c}d\psi.$

The middle term is zero because $d\psi$ applied to a tangent to $\gamma$ is
zero. We are left with

$dL\ \xi=\frac{e}{T}\langle\dot{\psi}\rangle.$

This proves that for passing particles on a rational flux surface,
$\langle\dot{\psi}\rangle=0$ iff $L$ is constant on it for given energy.

###  C.2. Omnigenity for bouncing trajectories implies omnigenity for
circulating ones on rational surfaces

We scale out the magnetic moment $\mu$ by writing $E=h\mu$ and
$L=\sqrt{m\mu}\,j$ with

$j(h)=\int_{\gamma}\sqrt{2(h-|B|)}\,ds,$

for segment $\gamma$ of fieldline. As remarked in section 2, the function $j$
is the Abel transform of the length function $\ell$ for a fieldline. $\ell(v)$
is the length of the subsegments of $\gamma$ for which $|B|<v$. Its Abel
transform (with a factor of $\sqrt{2}$) is

(41) $j(h)=\int_{h_{\min}}^{h}\sqrt{2(h-v)}\,d\ell(v),$

where $h_{\min}$ is the minimum of $|B|$ along it (one can replace $h_{\min}$
by $-\infty$ because $d\ell(v)=0$ for $v<h_{\min}$). This applies equally well
to a closed fieldline with $h>h_{\max}$, the maximum of $|B|$ on it.

The condition that $j$ be the same for all bouncing trajectories of given
class on a flux surface with given $h$, for all values of $h$ for this class,
implies that the length function $\ell$ is the same for each segment of
fieldline for the class on the flux surface, by the Abel inversion formula,
e.g. [Ke]:

$\ell(v)=\frac{2}{\pi}\int_{h_{\min}}^{v}\frac{dj(h)}{\sqrt{2(v-h)}}.$

This is a reformulation of a result of [SS] (see also [CS]). One could write
in more detail about the treatment of multiple classes of bouncing trajectory
[PCHL], but hopefully the above is clear enough. Note that it includes that
$h_{\min}$ and $h_{\max}$ are the same for all fieldlines on the same flux
surface.

Our point is that it follows from the length function being the same for each
bouncing class and the Abel representation (41) that $j(h)$ is the same for
all circulating trajectories on the rational surface with given $h>h_{\max}$.

###  C.3. Omnigenity for circulating trajectories on irrational surfaces

ZGCM at given energy for circulating particles on an irrational flux surface
in given direction is uniquely ergodic (i.e. it has a unique ergodic
probability measure). Up to normalisation, the unique ergodic probability
measure is given by the area-form $\frac{|B|}{v_{\parallel}}\mathcal{A}$
(where $\mathcal{A}$ was defined in Appendix B). This is because $\mathcal{A}$
is preserved by $B$ but the time spent anywhere by ZGCM is a factor
$|B|/v_{\parallel}$ longer than for fieldline flow. Thus, the time-average of
any continuous function along ZGCM at energy $E$ is equal to its surface-
average with respect to $\frac{|B|}{v_{\parallel}}\mathcal{A}$. On a flux
surface, $V\cdot\nabla\psi\ \mathcal{A}=i_{V}\Omega$ because
$\Omega=\mathcal{A}\wedge d\psi$ and $d\psi$ on tangents to a flux surface is
zero. Now using the non-Hamiltonian form for $V$, (40) and other relations as
above, we obtain
$\frac{|B|}{v_{\parallel}}i_{V}\Omega=\frac{|B|}{\tilde{B}_{\parallel}}i_{B}\Omega+\frac{m}{e}d(v_{\parallel}b^{\flat})$.
The first term is zero on tangents to the flux surface and the integral of the
second term over the flux surface is zero.

One might worry that the time required for convergence to the time-average
might be long for some irrationals, but actually under the omnigenity
condition we believe a uniform estimate is possible.

It is curious that the results of this subsection and subsection C.1 are exact
for the non-Hamiltonian form for $V$, but would be true only to first order if
we used the Hamiltonian form (6).

##  Appendix D Relation to pseudo-symmetry

A weaker notion than omnigenity was introduced by Mikhailov, e.g. [M+], called
pseudo-symmetry. There are various formulations, e.g. [Sk], but the way we
choose is that a magnetic field is said to be pseudo-symmetric if it has a
flux function and on each flux surface the contours of $|B|$ are nowhere
tangent to the magnetic field.

A pseudo-symmetric field is not necessarily isodrastic. Although the above
formulation constrains the local maxima of $|B|$ for a pseudo-symmetric field
to form (non-contractible) closed curves on each flux surface, it does not
force the separatrix area $j$ to be the same for each homoclinic orbit coming
from a local maximum with the given value of $|B|$.

Conversely, an isodrastic field is not necessarily pseudo-symmetric. Firstly,
an isodrastic field need not have a flux function, though if the marginal
trajectories are all heteroclinic then it does ($h=|B|$ on $\Sigma^{-}$
extended along the fieldlines is a flux function). Secondly, even if an
isodrastic field has a flux function, the contours of $|B|$ on it could have
an island chain instead of a curve of minima, as long as the fieldlines see
only local minima on crossing the island chain (else new $\Sigma^{-}$ is
created).

##  Appendix E $C^{4}$-generic $\Sigma^{0}$

A property is generic in a topological space if it happens on a countable
intersection of open dense subsets.

######  Theorem 6.

Restricting attention to bounded subsets, $C^{4}$-generically $\Sigma$ is a
$C^{3}$ surface and $\Sigma^{0}$ is a $C^{2}$-curve on $\Sigma$, separating it
into $\Sigma^{\pm}$.

###### Proof.

First recall that $\Sigma$ is defined by $|B|^{\prime}=0$ so if $B$ is $C^{r}$
with $r\geq 2$ then by the submersion theorem [La] the subset where
$d|B|^{\prime}\neq 0$ is a $C^{r-1}$ surface. The joint condition that
$|B|^{\prime}=0$ and $d|B|^{\prime}=0$ is four conditions ($\mbox{div}\ B=0$
does not make any restriction) on three variables, so generically does not
happen. Thus $C^{r}$-generically, $\Sigma$ is a $C^{r-1}$-surface. Note that
the condition $d|B|^{\prime}\neq 0$ is met in particular on $\Sigma^{\pm}$,
which are the subsets where $|B|^{\prime\prime}>0,<0$ respectively.
Furthermore, it follows that $\Sigma^{\pm}$ are transverse to $B$.

Secondly, the set $\Sigma^{0}$ is defined by
$|B|^{\prime}=0,|B|^{\prime\prime}=0$, so by the submersion theorem the subset
where the derivative $D(|B|^{\prime},|B|^{\prime\prime})$ has rank 2 is a
$C^{r-2}$ curve. The joint condition that
$|B|^{\prime}=0,|B|^{\prime\prime}=0$ and $D(|B|^{\prime},|B|^{\prime\prime})$
does not have rank 2 is 4 conditions on 3 variables, so generically does not
happen. Thus $C^{r}$-generically, $\Sigma^{0}$ is a $C^{r-2}$ curve.
$\Sigma^{0}$ lies on $\Sigma$ because $|B|^{\prime}=0$ on it. It separates
$\Sigma$ into $\Sigma^{\pm}$ because it is the subset with
$|B|^{\prime\prime}=0$. ∎

Next we study the generic behaviour of $|B|$ near $\Sigma^{0}$. Firstly, near
a generic point of $\Sigma^{0}$, there is a choice of fieldline labels $x,y$
and a coordinate $t$ along fieldlines, near to arc length $s$, such that

(42) $|B|=f(x,y)+yt+xt^{2}+kt^{3}$

for some function $f$ and constant $k\neq 0$. Note that there is no remainder
in this expression; this can be achieved by normal form principles of
singularity theory [Mi]. The restriction $\mbox{div}\ B=0$ plays no role in
this, only in the mapping between $(x,y,t)$ and physical space (compare the
proof of realisability in Section 4.2). Then
$|B|^{\prime}=(y+2xt+3kt^{2})t^{\prime}$, where $t^{\prime}$ denotes
$\frac{\partial t}{\partial s}$ holding $x,y$ constant, so $\Sigma$ is locally
the surface

(43) $y=-2xt-3kt^{2}.$

Thus we can use $(x,t)$ as coordinates on $\Sigma$. Furthermore

$|B|^{\prime\prime}=(2x+6kt)t^{\prime
2}+(y+2xt+3kt^{2})t^{\prime\prime}=(2x+6kt)t^{\prime
2}+|B|^{\prime}t^{\prime\prime}/t^{\prime}.$

Now $|B|^{\prime}=0$ on $\Sigma$, and $t^{\prime}$ is near $1$, so
$\Sigma^{0}$ is locally the curve

$x=-3kt,\ y=3kt^{2}.$

Its projection to fieldline labels $(x,y)$ is the curve $3ky=x^{2}$. It is a
standard fold in $\Sigma$. $\Sigma^{\pm}$ are the parts of $\Sigma$ with
$x+3kt>0,<0$ respectively. The result is shown in Figure 24. We see that there
are short bouncing segments on some fieldlines, namely those for which the
cubic has a well, i.e. for which $y<-2xt-3kt^{2}$.

Figure 24. Generic neighbourhood of a point $(0,0,0)$ on $\Sigma^{0}$. $x,y$
are fieldline labels and $s$ is arclength along the field. It shows that
$\Sigma^{0}$ forms a smooth (blue) curve on $\Sigma$, separating $\Sigma$ into
$\Sigma^{\pm}$ on the right and left, respectively. Three fieldlines are
indicated, with $x=0$, $y=-1,0,+1$ and the shapes of $|B|$ along them are
sketched.

Secondly, we have to take care of exceptional points $p$ of $\Sigma^{0}$ at
which $k=0$. Then (42) is modified generically to

$|B|=f(x,y)+yt+xt^{2}+at^{4}$

for some $a\neq 0$. Repeating the analysis, we obtain that $\Sigma$ is locally
$y=-2xt-4at^{3}$, and $\Sigma^{0}$ is locally $x=-6at^{2},\ y=8at^{3}$. Its
projection to fieldline labels $(x,y)$ is a standard semi-cubic cusp:
$27ay^{2}+16x^{3}=0$.

We can combine the treatment of generic points of $\Sigma^{0}$ and those
having $k=0$ by including both the $kt^{3}$ and $at^{4}$ terms in the
expression for $|B|$, with $k,a$ not simultaneously zero. Then on $\Sigma$,

$|B|=f(x,-2xt-3kt^{2}-4at^{3})-xt^{2}-2kt^{3}-3at^{4}.$

In particular, along $\Sigma^{0}$,

$|B|=f(-3kt-6at^{2},3kt^{2}+8at^{3})+kt^{3}+3at^{4},$

so using $t$ now as a coordinate along $\Sigma^{0}$,
$\frac{d|B|}{dt}=-3kf_{,x}$ at $t=0$, where $f_{,x}$ denotes the partial
derivative of $f$ with respect to $x$ at $(0,0)$. At points of $\Sigma^{0}$
where $kf_{,x}\neq 0$ the field is not weakly isodrastic, because $|B|$ is not
constant along $\Sigma^{0}$ (a little more analysis shows that the level sets
of $|B|$ cross $\Sigma^{0}$ with cubic tangency). A closed component of
$\Sigma^{0}$ has at least one minimum and one maximum of $|B|$, so there are
points on it where $f_{,x}=0$ or $k=0$. A question is whether both types of
point have to occur, but generically they do not happen simultaneously. We
analyse the local picture for the $|B|$ levels on $\Sigma$ near such points.

At a point of $\Sigma^{0}$ where $f_{,x}=0$ one has generically a non-
degenerate local maximum or minimum of $|B|$ along $\Sigma^{0}$. The level
sets of $|B|$ on $\Sigma$ look locally like Figure 25.

Figure 25. Level sets of $|B|$ on $\Sigma$ for $f_{,x}=0$ in the case that
$f_{,y}>0$ and (a) $f_{,xx}>0$, (b) $f_{,xx}<0$ where $H$ and $L$ denote high
and low regions of $|B|$.

So the level curves of $|B|$ cross $\Sigma^{0}$ except in the hyperbolic
sector of the first case.

For points of $\Sigma^{0}$ where $k=0$, note that on $\Sigma^{0}$, $k=0$ iff
$B$ is tangent to $\Sigma^{0}$. So one has generically one of the cases of
Figure 26.

Figure 26. The disposition of $\Sigma$ and $\Sigma^{0}$ relative to the
magnetic field lines near a point $x=y=0$ with $k=0$; (a) quartic minimum, (b)
quartic maximum. In each case, three fieldlines are shown (denoted by red,
green, and blue lines) and a bouncing segment on each of them is denoted by
bar arrows. The $|B|$ profiles along the fieldlines are shown in the bottom
panel, together with the chosen level of $E/\mu$ for the bouncing segment.

In the case of $|B|$ coming to a quartic minimum, all level curves of $|B|$
collide with $\Sigma^{0}$. The case of $|B|$ coming to a quartic maximum
produces a curve of switch of lowest hill. This “Maxwell curve” is $y=0$,
$x=2t^{2}$. All $|B|$ levels on $\Sigma^{+}$ cross $\Sigma^{0}$.

If we want weak isodrasticity and there is a curve of $\Sigma^{0}$ then we
must have $kf_{,x}=0$ everywhere on it.

To analyse the dynamics more thoroughly near $\Sigma^{0}$ one needs to compute
$j$. In particular, to test isodrasticity near but not on $\Sigma^{0}$ one
needs ȷ on $\Sigma^{-}$. We can compute this to a good approximation for the
short homoclinics that occur near the generic points of $\Sigma^{0}$ in the
well of the cubic. From a point of $\Sigma^{-}$ labelled by $(x,t)$, we use
$t^{\prime}(s)\approx 1$ to get the approximation

$\text{\j}=\int_{s_{0}}^{s_{1}}\sqrt{2(C-ys-xs^{2}-ks^{3})}\,ds,$

where $C$ is the value at the local maximum $s_{0}$ of $ys+xs^{2}+ks^{3}$ and
$s_{1}$ is the other point at the same height. Write $s=s_{0}+v/k$ and
$\alpha=\sqrt{-|B|^{\prime\prime}/2}$ at the point of $\Sigma^{-}$ and recall
that $|B|^{\prime\prime}=2x+6ks_{0}$. Then

$\text{\j}=\frac{\sqrt{2}}{k^{2}}\int_{0}^{\alpha^{2}}v\sqrt{\alpha^{2}-v}\,dv.$

Put $u^{2}=\alpha^{2}-v$, to obtain

$\text{\j}=\frac{2\sqrt{2}}{k^{2}}\int_{0}^{\alpha}(\alpha^{2}-u^{2})u^{2}\,du=\frac{(-|B|^{\prime\prime})^{5/2}}{15k^{2}}.$

Using $|B|^{\prime\prime}=2x+6kt$ on $\Sigma$, we obtain

$d\text{\j}=\frac{(-|B|^{\prime\prime})^{3/2}}{3k^{2}}(dx+3k\,dt).$

Recalling the expression for $h=|B|$ on $\Sigma$,

$dh=f_{,x}dx+f_{,y}(-2tdt-(2x+6kt)dt)-t^{2}dx-(2xt+6kt^{2})dt.$

Thus

$dh\wedge
d\text{\j}=\frac{(-|B|^{\prime\prime})^{3/2}}{3k^{2}}(-2xf_{,y}-3kf_{,x}+y)\,dx\wedge
dt,$

using $y=-2xt-3kt^{2}$ (43). It is zero iff $y=2xf_{,y}+3kf_{,x}$. This PDE
can be solved for $f$, resulting in $h$ being an arbitrary function of
$x^{2}-3ky$ (the deviation from $\Sigma^{0}$), but the main conclusion can be
read off immediately by putting $(x,y)=(0,0)$: $kf_{,x}=0$ at $(0,0)$. This
implies for $k\neq 0$ that the derivative of $h$ along $\Sigma^{0}$ is zero,
confirming one of our necessary conditions for weak isodrasticity.

It is useful to compare $dh\wedge d\text{\j}$ to $\beta$. On $\Sigma$, locally
$\beta=c(x,y)\,dx\wedge dy$ for some non-zero smooth function $c$, because $x$
and $y$ are fieldline labels. Using $y=-2xt-3kt^{2}$ and the expression for
$|B|^{\prime\prime}$, we obtain $dx\wedge dy=-c|B|^{\prime\prime}dx\wedge dt$.
So the Melnikov function

$\mathcal{M}=\frac{(-|B|^{\prime\prime})^{1/2}}{3ck^{2}}(y-2xf_{,y}-3kf_{,x}).$

Finally, for isodrasticity one also needs the corresponding condition for
segments leaving $\Sigma^{-}$ in the opposite direction, which is not
obtainable by local analysis. A typical picture for the phase space $F_{j}$
around a fieldline with a cubic critical point is sketched in Figure 27, using
the normal form (42) and a generic assumption that the separatrix-area of the
cubic critical point varies at non-zero rate along the fold curve $\Sigma^{0}$
(but note that this assumption is incompatible with weak isodrasticity,
because we proved that $\j$ is constant along $\Sigma^{0}$ for weak isodrastic
fields.)

Figure 27. Sketch of the phase space $F_{j}$ near a fieldline with a cubic
critical point with generic unfolding. $\Sigma^{0}$ is the curve with cubic
critical points: $3ky=x^{2}$ in fieldline labels $(x,y)$. The cusped wedge is
excluded.

The cubic critical point unfolds in the downward vertical direction, to a
local maximum and a local minimum. Moving horizontally to the right, the
separatrix-area for the cubic critical point decreases, thus to maintain
constant $j$, the energy-level changes as indicated. The boundary of $F_{j}$
is a cusped curve with width $\Delta x$ asymptotically proportional to
$(-y)^{5/4}$. This follows from the above computation of $\j$ for the little
well in the cubic. The lefthand curve is determined by making the separatrix
area for the main well be $j$; on the righthand curve the sum of the
separatrix areas for the main well and the little well is $j$.

In conclusion, generic $\Sigma^{0}$ is incompatible with isodrasticity except
with special design. The typical ways isodrasticity fails near $\Sigma^{0}$
reveals a sensitivity of axisymmetric and quasisymmetric designs to
imperfections.

##  Appendix F Derivative of reduced Hamiltonian

Here we give an alternative proof of Proposition 1 that illustrates its
general connection to the theory of nearly-periodic systems.

###### Proof.

In the bouncing region of the guiding center phase space the guiding center
vector field $V_{\epsilon}$ defines a Hamiltonian nearly-periodic system with
limiting roto-rate $R_{0}=\frac{T}{2\pi}\,V_{0}$ and exact
$\epsilon$-dependent symplectic form $\omega_{\epsilon}$. Here Hamiltonian
means $\iota_{V_{\epsilon}}\omega_{\epsilon}=dH_{\epsilon}$, for some
Hamiltonian $H_{\epsilon}$, and $T$ denotes the true period for bounce motion.
As for all Hamiltonian nearly-periodic systems [BS], there exists an all-
orders roto-rate $R_{\epsilon}$ such that
$\iota_{R_{\epsilon}}\omega_{\epsilon}=dJ_{\epsilon}$, where $J_{\epsilon}$
denotes the all-orders bounce adiabatic invariant. By [BS] the formula for the
first-order term in the roto-rate is

$\displaystyle R_{1}=\mathcal{L}_{R_{0}}I_{0}\widetilde{V}_{1},$

where $I_{0}$ denotes the inverse of $\mathcal{L}_{V_{0}}$ restricted to the
subspace of vector fields with vanishing $U(1)$-average, and
$\widetilde{V}_{1}=V_{1}-\langle V_{1}\rangle$ is the first-order guiding
center vector field $V_{1}$ less its $U(1)$-average. It is easy to show that
$H_{0}=J_{0}=0$, and that $H_{1},J_{1}$ correspond to the usual leading-order
expressions for the guiding center energy and bounce invariant. Since both
$R_{\epsilon}$ and $V_{\epsilon}$ are Hamiltonian vector fields we have

$\displaystyle\iota_{V_{1}}\omega_{0}+\iota_{V_{0}}\omega_{1}$
$\displaystyle=dH_{1}$
$\displaystyle\iota_{R_{0}}\omega_{0}+\iota_{R_{0}}\omega_{1}$
$\displaystyle=dJ_{1}.$

In particular
$dH_{1}=\iota_{V_{1}}\omega_{0}+\frac{2\pi}{T}(dJ_{1}-\iota_{R_{1}}\omega_{0})$.
This formula is useful because it simplifies computing the derivative of
$H_{1}$ along $U(1)$-invariant vector fields
$\mathcal{X}=\bm{u}\cdot\partial_{\bm{X}}+a\,\partial_{v_{\parallel}}$ that
preserve $J_{1}$, i.e. $dJ_{1}(\mathcal{X})=0$. In particular,

$\displaystyle dH_{1}(\mathcal{X})$
$\displaystyle=\omega_{0}(V_{1},\mathcal{X})-\frac{2\pi}{T}\omega_{0}(R_{1},\mathcal{X})$
$\displaystyle=\omega_{0}(V_{1},\mathcal{X})-\frac{2\pi}{T}\omega_{0}(\mathcal{L}_{R_{0}}I_{0}\widetilde{V}_{1},\mathcal{X})$
$\displaystyle=\omega_{0}(V_{1},\mathcal{X})-\frac{2\pi}{T}\mathcal{L}_{R_{0}}I_{0}\bigg{[}\omega_{0}(\widetilde{V}_{1},\mathcal{X})\bigg{]}$
$\displaystyle=\omega_{0}(V_{1},\mathcal{X})-\omega_{0}(\widetilde{V}_{1},\mathcal{X})$
$\displaystyle=\langle\omega_{0}(V_{1},\mathcal{X})\rangle$
$\displaystyle=\frac{1}{T}\int_{0}^{T}\bigg{(}\mu\bm{u}\cdot\nabla_{\perp}|B|-mv_{\parallel}^{2}\,\bm{b}\cdot\bm{c}_{\perp}\times\bm{u}\bigg{)}\,dt.$

This establishes the Proposition after using $dt=b^{\flat}/v_{\parallel}$ and
energy conservation. ∎

##  Appendix G Double transitions

As a special event, a fieldline segment can approach double transition. Some
ways this can occur are indicated in Figure 28. They involve formation of a
heteroclinic cycle. There are others involving formation of a degenerate
critical point too, but they are more special.

Figure 28. Some ways to double transition.

In the 2D space of fieldlines, a heteroclinic cycle requires only one
condition, namely that $|B|$ have two local maxima at the same height, so it
happens generically along curves in the space of fieldlines.

It leads to corners in the reduced spaces $F_{j}$, where boundaries
corresponding to two different transitions meet. See Figure 29.

Figure 29. An example of a corner in the reduced space $F_{j}$.

One could make explicit examples, for example with $|B|$ a quartic in
arclength whose coefficients depend on two fieldline labels, as for the
formulae leading to Figure 26(a). Although evaluating $j$ requires elliptic
integrals, the boundary cases can be computed explicitly (generalising the
cubic case of Section 4 and Appendix E).

Near points of double transition, transitions between several classes can
occur. Isodrasticity requires that $\ker dH_{j}$ contain the tangent to the
boundary at single transition points. Taking the limit to a transverse corner,
this implies that $dh=0$ there. The consequences of this are left to a future
publication.

Note that the times spent by a periodic orbit near the saddles of a
heteroclinic cycle of a Hamiltonian system are asymptotically proportional to
the Lyapunov times of the saddles (i.e. the inverses of their positive
Lyapunov exponents). So using (15), along the curve for equal height, $dH_{j}$
is asymptotic to the convex combination of $dh$ at the two ends, weighted by
their Lyapunov times.

Note also that $d\j$ goes to infinity at generic corners because the time
spent near the second saddle grows logarithmically as it is approached along a
level curve of $\j$, and the derivative of $\j$ is related to this time.

##  Appendix H Dipole field

As a simple illustration of $\Sigma$ and the reduced Hamiltonian $H_{j}$, we
consider GCM in a dipole field

$B=\frac{3\cos\theta\ \hat{r}-\hat{z}}{r^{3}}$

in spherical polar coordinates (radius $r$, colatitude $\theta$, longitude
$\phi$), illustrated in Figure 30. The dipole strength has been scaled to
$4\pi$.

Figure 30. A meridional section through a dipole field, showing some
fieldlines and the field strength in colour (blue to red indicates weak to
strong and white indicates $|B|$ exceeds a threshold). $\Sigma^{+}$ is the
equatorial plane.

The fieldlines are $r=r_{e}\sin^{2}\theta$, $\phi=$ constant, where $r_{e}$ is
the radius at which they cross the equatorial plane. Along the fieldlines,

$|B|=\frac{\sqrt{3\cos^{2}\theta+1}}{r_{e}^{3}\sin^{6}\theta}.$

The guiding centres bounce across the equatorial plane between the regions of
stronger field near the poles. Then $\Sigma$ is the equatorial plane and it
consists of only $\Sigma^{+}$ (see Figure 2(a)).

For a bouncing segment, $j$ can be written as $j={r_{e}^{-1/2}}F(hr_{e}^{3})$,
where $h$ is $|B|$ at the bounce points and

$F(k)=\int_{\theta_{0}}^{\theta_{1}}\sqrt{2\left(k-\frac{\sqrt{3\cos^{2}\theta+1}}{\sin^{6}\theta}\right)}\sqrt{2\cos^{2}\theta+1}\,\sin\theta\,d\theta,$

with $\theta_{i}(k)$ being the zeroes of the first square root (symmetric
about $\pi/2$). $F$ is defined for argument greater than or equal to $1$, has
$F(1)=0$, positive derivative (including at $k=1$ where
$F^{\prime}(1)=\tfrac{\pi}{3}$), and goes to infinity as $k\to\infty$. FGCM is
described in the adiabatic approximation with given value of $j$ (and scaled
time $\tau$) by the Hamiltonian $H_{j}$ on $\Sigma^{+}$, defined by
$H_{j}=r^{-3}F^{-1}(\sqrt{r}j)$ (we write $r_{e}=r$ on $\Sigma^{+}$), and the
symplectic form that is just the magnetic flux-form $r^{-2}dr\wedge d\phi$ on
$\Sigma^{+}$. Now

$\frac{dH_{j}}{dr}=-\frac{3}{r^{4}}F^{-1}(\sqrt{r}j)+\frac{j}{2r^{7/2}F^{\prime}(k)}=r^{-4}\left(-3k+\frac{F}{2F^{\prime}}\right),$

where $k=H_{j}r^{3}$. We didn’t check, but presumably $F/F^{\prime}<{6k}$
(this is certainly true for $k=1$ where $F/F^{\prime}$ is zero, and at
$k=\infty$ where $F/F^{\prime}\sim{2k}$) so $H_{j}$ is a strictly decreasing
function of $r$ for given $j$. In any case, the level sets of $H_{j}$ are
concentric circles, and the bouncing segments precess around the dipole axis
at constant rate $\frac{d\phi}{d\tau}=r^{2}\tfrac{dH_{j}}{dr}$.

If axisymmetry is broken but not too much in $C^{1}$, as perhaps for the
earth’s magnetosphere, and we ignore changes far away (such as due to the
solar wind), then $\Sigma^{+}$ deforms into a nearby surface and $H_{j}$ into
a nearby function. Its level sets deform to closed curves near the original
circles. The segments continue to precess, but in general no longer at
constant rate. The precession period can be calculated from (18).

As there is no $\Sigma^{-0}$, there are no transitions and the field is
automatically isodrastic. The only thing to check for confinement is which set
of precessing segments to populate. For example, in the context of the earth’s
magnetosphere, the desired region is the outside of the earth. Then in the
axisymmetric case, scaling the earth’s radius to 1 and the field strength to 1
at the earth’s equator to fit with the above, a calculation shows that the
segments that do not hit the earth are those for which $h<\sqrt{4-3/r_{e}}$,
where $h$ is the ratio of the energy to $\mu$.

The earth’s magnetosphere is perturbed not only by deviations from a dipole of
magnetic generation in the earth but also by interaction with the magnetic
field of the solar system and solar wind, which change the arrangement of the
fieldlines, notably introducing magnetic nulls. We do not pursue those effects
here.

##  Appendix I Computational practicalities

Firstly, we show how square roots can be avoided in computation of many of the
quantities required to find $\Sigma$ and the Melnikov function (though not
all). Then we give suggestions for the computation of ȷ on $\Sigma^{-}$.

###  I.1. Eliminating square roots

We have written expressions like $|B|^{\prime}=b\cdot\nabla|B|$ but if $B$ is
given in components, $|B|$ involves taking a square root, so $b=B/|B|$
involves dividing by a square root. So does differentiating $|B|$. It can be
better to write such quantities in terms of $|B|^{2}$. For example,

$|B|^{\prime}=\tfrac{1}{2}B\cdot\nabla|B|^{2}/|B|^{2}.$

Here, we collect various such formulae.

The grad-B drift involves $b\times\nabla|B|$, which can be written as
$\frac{1}{2}B\times\nabla|B|^{2}/|B|^{2}$.

The curvature drift involves $(\mbox{curl}\,b)_{\perp}/|B|$. Now

$\mbox{curl}\,b=\mbox{curl}\,\frac{B}{|B|}=\frac{J}{|B|}-\frac{\nabla|B|}{|B|^{2}}\times
B=\frac{J}{|B|}+\tfrac{1}{2}B\times\frac{\nabla|B|^{2}}{|B|^{3}}.$

So

$(\mbox{curl}\
b)_{\perp}/|B|=J_{\perp}/|B|^{2}+\tfrac{1}{2}B\times\nabla|B|^{2}/|B|^{4}.$

We can also remove square roots from the computation of the covector $a$ in
(38) by switching to fieldline-flow time $T$. We have
$da_{i}/dT=-|B|a_{j}\partial_{i}b^{j}$ and
$\partial_{i}b^{j}=\frac{1}{|B|}\partial_{i}B^{j}-\frac{B^{j}}{2|B|^{3}}\partial_{i}|B|^{2}$.
So

$\frac{da_{i}}{dT}=-a_{j}\left(\partial_{i}B^{j}-\frac{B^{j}}{2|B|^{2}}\partial_{i}|B|^{2}\right)=-a_{j}\partial_{i}B^{j},$

because $aB=0$.

###  I.2. Integrating along a fieldline

Next, ȷ is defined at a point $X_{0}\in\Sigma^{-}$ by taking $h=|B(X_{0})|$
and then letting $\text{\j}=\int\sqrt{2(h-|B|)}\,ds$ with respect to arclength
$s$ along the fieldline through $X_{0}$ to the first bounce, i.e. where
$|B|=h$ again.

Suppose we are integrating in the positive direction along $B$ (the obvious
changes apply for the other direction). It is slightly more convenient to
integrate with respect to fieldline flow time $T$ than arclength $s$, i.e.
start at $T=0$ with $j=0$ and $X=X_{0}$ and integrate

$\frac{dX}{dT}=B(X),\quad\frac{dj}{dt}=v|B|,$

with $v=\sqrt{2(h-|B|)},$ until the first bounce, thereby banishing square
roots to only the second of the two equations.

More importantly, it is better to switch to integration with respect to $w=-v$
when approaching the first bounce, i.e. by eliminating time from ZGCM,

$\frac{dX}{dw}=-\frac{wb}{|B|^{\prime}},\quad\frac{dj}{dw}=\frac{w^{2}}{|B|^{\prime}},$

starting from $w=-\sqrt{2(h-|B|)}$ at the switch point and stopping at $w=0$.
Then put ȷ equal to the final value of $j$. The switch can be made at any
point between the last local minimum of $|B|$ along the fieldline and the
bounce point.

Nonetheless, in the codimension-1 case that the first place where $|B|=h$
again is also a zero of $|B|^{\prime}$, this switch to integration with
respect to $w$ is not appropriate. Instead one can just integrate with respect
to $T$ with termination condition $|B|^{\prime}=0$ when near $B=h$ again.

These tips are also relevant for computing the Melnikov function.

##  Appendix J Perturbed toy tokamak

Here we give the expressions used in the calculations for the weak and strong
form of isodrasticity for the toy tokamak.

The perturbation to break the axisymmetric field is given by

$\displaystyle\nabla\times
A_{\varepsilon}=\left(-{\varepsilon\cos\phi},0,\dfrac{\varepsilon
z\cos\phi}{R}\right),$

where $\varepsilon$ is the perturbation parameter (and not the GC parameter in
Sec. 1).

The non-axisymmetric field magnitude,

$\displaystyle|B|=\dfrac{\sqrt{C^{2}+(\varepsilon
R\cos\phi+z)^{2}+(R-R_{0}+\varepsilon z\cos\phi)^{2}}}{R},$

and

$\displaystyle\dfrac{d|B|}{ds}$
$\displaystyle=|B|^{\prime}=\frac{f_{1}+f_{2}+f_{3}}{R^{2}\left(C^{2}+\left(R\varepsilon\cos\phi+z\right)^{2}+\left(R-R_{0}+\varepsilon
z\cos\phi\right)^{2}\right)},$ $\displaystyle\text{where}\;f_{1}$
$\displaystyle=-C\varepsilon\left(R\left(R\varepsilon\cos\phi+z\right)+z\left(R-R_{0}+\varepsilon
z\cos\phi\right)\right)\sin\phi$ $\displaystyle f_{2}$
$\displaystyle=R\left(R-R_{0}+\varepsilon
z\cos\phi\right)\left(R\varepsilon\cos\phi+\varepsilon\left(R-R_{0}+\varepsilon
z\cos\phi\right)\cos\phi+z\right)$ $\displaystyle f_{3}$
$\displaystyle=\left(R\varepsilon\cos\phi+z\right)\left(C^{2}-R_{0}(R-R_{0})+2(R-R_{0})\varepsilon
z\cos\phi+\varepsilon^{2}z^{2}\cos^{2}{\left(\phi\right)}+z^{2}\right).$

In $\tilde{B}_{\parallel}=|B|+\sqrt{\tilde{\mu}}u(\mbox{curl}\,b)\cdot b$, we
rewrite

$(\mbox{curl}\,b)\cdot b=\frac{-CR-CR_{0}+C\varepsilon
z\cos\phi-R(R-R_{0})\varepsilon\sin\phi+\varepsilon
z^{2}\sin\phi}{R\left(C^{2}+\varepsilon^{2}\cos^{2}{\left(\phi\right)}(R^{2}+z^{2})+(R-R_{0})^{2}+2\varepsilon
z\cos\phi(2R-R_{0})+z^{2}\right)}.$

##  Appendix K Slow manifold computation

The procedure from [M04] to compute a first-order symplectic slow manifold
$N^{-}$ for a Hamiltonian system from a zeroth-order one
($\Sigma^{-}\times\\{u=0\\}$ in our case) is to compute a symplectically
orthogonal foliation to the latter and to find the unique nearby critical
point of $H$ on each leaf.

The tangents to $\Sigma^{-}\times\\{u=0\\}$ are the vectors $(\delta x,\delta
u)\in\mathbb{R}^{3}\times\mathbb{R}$ satisfying $\delta
x\cdot\nabla|B|^{\prime}=0,\delta u=0$. The symplectic form is
$\tilde{\omega}=\tilde{\mu}^{-1/2}\beta+u\,db^{\flat}+du\wedge b^{\flat}$.
This requires $\mu>0$ but the method works the same. Thus a vector $(\xi,w)$
is symplectically orthogonal to $\Sigma^{-}\times\\{u=0\\}$ iff

$\tilde{\mu}^{-1/2}\beta(\xi,\delta x)+wb\cdot\delta x=0$

for all tangents $\delta x$ to $\Sigma^{-}$. Recall that $\beta=i_{B}\Omega$.
For a generic point of $\Sigma^{-}$ where $b$ is not perpendicular to
$\Sigma^{-}$, the choice $\delta x=b\times\nabla|B|^{\prime}$ shows that
$\xi=\alpha b+\gamma b\times\nabla|B|^{\prime}$ for some $\alpha$ and
$\gamma$. Next, the choice $\delta
x=\nabla|B|^{\prime}\times(b\times\nabla|B|^{\prime})$ yields

(44) $w=-\tilde{\mu}^{-1/2}\frac{\beta(\gamma
b\times\nabla|B|^{\prime},\nabla|B|^{\prime}\times(b\times\nabla|B|^{\prime}))}{b\cdot(\nabla|B|^{\prime}\times(b\times\nabla|B|^{\prime}))}=-\tilde{\mu}^{-1/2}\gamma|B||B|^{\prime\prime},$

after some cancellation.

We take the locally linear foliation given by the above symplectically
orthogonal planes to $\Sigma^{-}\times\\{u=0\\}$. So the leaf associated to a
point $x\in\Sigma^{-}$ is parametrised by displacement $\alpha$ from
$\Sigma^{-}$ along $b$, $\gamma b\times\nabla|B|^{\prime}$ tangent to
$\Sigma^{-}$ and scaled parallel velocity
$w=-\tilde{\mu}^{-1/2}\gamma|B||B|^{\prime\prime}$.

The Hamiltonian $H=\tfrac{1}{2}w^{2}+|B|$ constrained to this leaf is

$H=\tfrac{\gamma^{2}}{2\tilde{\mu}}|B|^{2}|B|^{\prime\prime 2}+|B(x+\alpha
b+\gamma b\times\nabla|B|^{\prime})|.$

Now $|B|^{\prime}=0$ on $\Sigma^{-}$, so neglecting the effect of curvature of
$\Sigma^{-}$, critical points with respect to displacement $\alpha$ along $b$
have $\alpha=0$ (it would be good to estimate the error). Criticality with
respect to $\gamma$ is given to leading order by

$\frac{\gamma}{\tilde{\mu}}|B|^{2}|B|^{\prime\prime
2}+(b\times\nabla|B|^{\prime})\cdot\nabla|B|=0.$

There’s also a term proportional to $\gamma$ from the second derivative of
$|B|$ but the $w$ term dominates it by the factor $1/{\tilde{\mu}}$. So we end
up with

$\gamma=-\tilde{\mu}\frac{b\times\nabla|B|^{\prime}\cdot\nabla|B|}{|B|^{2}|B|^{\prime\prime
2}}.$

In particular, using (44),

$w=\sqrt{\tilde{\mu}}\frac{b\times\nabla|B|^{\prime}\cdot\nabla|B|}{|B||B|^{\prime\prime}}.$

##  Appendix L Systems with perfect separatrices

There exist Hamiltonian systems that have one or more perfect separatrices but
are not integrable.

Firstly, one can make area-preserving twist maps with this property. The idea
is to use the construction by de la Llave in the appendix to [Mat]. Given a
(lift of a) degree-one homeomorphism $g$ of the circle
$\mathbb{R}/\mathbb{Z}$, define

$h(x)=g(x)+g^{-1}(x)-2x,$

and define map $T$ for $(x,y)$ on the cylinder
$\mathbb{R}/\mathbb{Z}\times\mathbb{R}$ by

$y^{\prime}=y+h(x),\quad x^{\prime}=x+y^{\prime}$

It is an area-preserving twist map and the circles $y=x-g(x)$ and
$y=x-g^{-1}(x)$ are invariant, with dynamics $x^{\prime}=g^{-1}(x)$ and
$x^{\prime}=g(x)$ respectively. So we choose $g$ to be a circle homeomorphism
with two fixed points, e.g. $g(x)=x+k\sin 2\pi x$ with $0<k<\tfrac{1}{2\pi}$.
Then we obtain a map $T$ with two period-one islands with perfect
separatrices. But in general it is non-integrable, as illustrated in Figure 31
(we presume that one could prove this if desired).

Figure 31. Some orbits of the map $T$ for $k=0.1$ with the $y=x-g^{-1}(x)$
invariant circle shown as the red and the $y=x-g(x)$ as the blue curve.

One could modify the choice of $g$ to make an example with a single period-one
island but it would require careful matching of derivatives of $g$ to the left
and right of the fixed point if one wants $T$ to be smooth.

Similarly, one can make continuous-time Hamiltonian systems with a perfect
separatrix without imposing integrability. For example, given a function
$S:\mathbb{R}^{2}\to\mathbb{R}$ with constants $A,B$ such that
$S(x+m,y+n)=S(x,y)+mA+nB$, and the vector field $\dot{q}=\nabla S(q)$ induced
on $\mathbb{R}^{2}/\mathbb{Z}^{2}$, then
$H=\tfrac{1}{2}|p|^{2}-\tfrac{1}{2}|\nabla S(q)|^{2}$ has invariant graph
$p=\nabla S(q)$ on which $\dot{q}=\nabla S(q)$. So choose $S$ with a periodic
orbit repelling on one side, attracting on the other. What is not clear to us
is how to make examples with perfect separatrices for all energies
simultaneously.

##  Appendix M Persistence of $\Sigma^{0}$ and $\Sigma^{+}$

In this appendix, we address the question of whether also $\Sigma^{0}$ and
$\Sigma^{+}$ have continuations to invariant submanifolds for small $\mu$.

We begin with $\Sigma^{0}$. If a 2DoF Hamiltonian system has a 2D manifold
$N^{0}$ consisting of elementary saddle-centre periodic orbits (“elementary”
means that some generic conditions are satisfied) then firstly it is part of a
3D manifold $N$ consisting of periodic orbits, which decomposes into
$N^{-},N^{0}$ and $N^{+}$, with the orbits on $N^{-}$ being hyperbolic and
those on $N^{+}$ being elliptic. Secondly, all nearby Hamiltonian systems have
a nearby such $N$ locally. These results follow from [Mey], which is
formulated in the context of area-preserving maps.

This is a useful result once we have $\mu>0$, but unfortunately it does not
apply to $\mu=0$ because firstly the points of $\Sigma^{0}\times\\{0\\}$ are
all equilibria, not periodic orbits. Secondly, the Poisson bracket is
degenerate at $\mu=0$ so we are not starting from a genuine Hamiltonian
system.

We suspect that under the conditions that $\Sigma^{0}$ be a generic curve and
$|B|$ be constant along its components then at least
$\Sigma^{-0}\times\\{0\\}$ persists to an invariant $N^{-0}$. The $|B|$
constant condition is necessary for the reduced dynamics for $\mu>0$ to have
$\Sigma^{0}$ invariant.

One might ask why we do not conjecture that the whole of $\Sigma\times\\{0\\}$
persists to an invariant $N$. The answer is that generically we expect
resonances to break $\Sigma^{+}\times\\{0\\}$. We give some explanation.

$\Sigma^{+}\times\\{0\\}$ is a slow manifold but normally elliptic instead of
hyperbolic. In general, the best one can deduce for a normally elliptic slow
manifold is that the true system has a sequence of submanifolds $N^{+}_{m}$
that are invariant to $n^{th}$ order, but in general the sequence does not
converge [M04]. In this low-dimensional case, however, the reduced dynamics on
$\Sigma^{+}$ consists of periodic orbits that are non-degenerate except at
integer resonances (in this context, integer resonances means the linearised
bounce frequency is an integer multiple of the precession frequency). If we
cut out a neighbourhood of these integer resonances, there is a true invariant
submanifold $N^{+}$ nearby consisting of periodic orbits of the guiding-centre
dynamics, which are elliptic except near half-integer resonances (where they
generically turn inversion hyperbolic). This procedure to construct $N^{+}$
fails at integer resonances (including at $\Sigma^{0}$ where the elliptic
frequency goes to zero).

Integer resonances all correspond to saddle-centre periodic orbits but they
are unlikely to be elementary: that would require in particular that the field
strength happen to have a turning point at them. Thus we expect $\Sigma^{+}$
to break at integer resonances, if there is something for them to resonate
with. For axisymmetric fields the resonances are not excited by turning on
$\mu>0$, but for general non-axisymmetric ones they are excited. The (angular)
bounce frequency is $\omega_{b}=\sqrt{\tfrac{\mu}{m}|B|^{\prime\prime}}$ and
the (angular) precession frequency is
$\omega_{p}=2\pi\tfrac{\mu}{e}\frac{\partial h}{\partial\Phi}$, so the
condition for integer resonance $n$ is

$\frac{\partial\Phi}{\partial h}\sqrt{|B|^{\prime\prime}}=2\pi
n\sqrt{\tilde{\mu}},$

where $\tilde{\mu}=\tfrac{m}{e^{2}}\mu$. Unfortunately, for typical fields
with the derivative of the lefthand side non-zero this means that there is a
large set of resonances for $\tilde{\mu}$ small. Nonetheless, if the lefthand
side is non-zero they have large $n$ and for smooth enough fields the
resonances can be expected to be very weak. See [M04] for an example where
they are exponentially weak for large $n$. Also, adiabatic invariance of $L$
typically makes the motion near $\Sigma^{+}$ bounded.

We close this appendix by illustrating the problem of integer resonances for
$\Sigma^{+}$ by the tokamak example of 3.2. Then $\Phi=2\pi\psi$ plus a
constant, so

$\frac{1}{2\pi}\frac{\partial\Phi}{\partial
h}\sqrt{|B|^{\prime\prime}}=-\frac{r^{3/2}R(C^{2}+r^{2})^{1/4}}{C^{2}-r}$

on $\Sigma^{+}$. So we get resonance $-n$ at $r\approx
C\tilde{\mu}^{1/3}n^{2/3}$. As $\mu\to 0$ each $-n$ resonance tends along
$\Sigma^{+}$ to the magnetic axis, so Figure 21 would generically show a
sequence of breaks on breaking axisymmetry.

## References

  * [Ar] Arnol’d VI, Instability of dynamical systems with many degrees of freedom, Sov Math Dokl 6 (1964) 581–5.
  * [BGKM] Baesens C, Guckenheimer J, Kim S, MacKay RS, Three coupled oscillators: Mode-locking, global bifurcations and toroidal chaos, Physica D 49 (1991) 387–475.
  * [Ba] Balescu R, Transport Processes in Plasmas, Neoclassical Transport Vol.II (North-Holland, 1988).
  * [B+] Beidler CD, Kolesnichenko YaI, Marchenko VS, Sidorenko IN, Wobig H, Stochastic diffusion of energetic ions in optimized stellarators, Phys Plasma 8 (2001) 2731–8.
  * [Bo] Boozer AH, Why carbon dioxide makes stellarators so important, Nucl Fusion 60 (2020) 065001.
  * [Bo84] Boozer AH, Time-dependent drift Hamiltonian, Phys Fluids 27 (1984) 2441–5.
  * [BSQ] Burby JW, Squire J, Qin H, Automation of the guiding centre expansion, Phys Plasma 20 (2013) 072105.
  * [BE] Burby JW, Ellison CL, Toroidal regularization of the guiding center Lagrangian, Phys Plasma 24 (2017) 110703.
  * [BS] Burby JW, Squire J, General formulas for adiabatic invariants in nearly-periodic Hamiltonian systems, J Plasma Phys 86 (2020) 835860601.
  * [BH] Burby JW, Hirvijoki E, Normal stability of slow manifolds in nearly-periodic Hamiltonian systems, J Math Phys 62 (2021) 093506.
  * [BKM] Burby JW, Kallinikos N, MacKay RS, Some mathematics for quasi-symmetry, J Math Phys 61 (2020) 093503.
  * [BKM2] Burby JW, Kallinikos N, MacKay RS, Approximate symmetries of guiding centre motion, J Phys A 54 (2021) 125202
  * [CS] Cary JR, Shasharina SG, Omnigenity and quasihelicity in helical plasma confinement systems, Phys Plasmas 4 (1997) 3323–33.
  * [CF] Cerfon AJ, Freidberg JP, “One size fits all” analytic solutions to the Grad-Shafranov equation, Phys Plasmas 17 (2010) 032502.
  * [Co] Conley CC, Low energy transit orbits in the restricted three-body problem, SIAM J Appl Math 16 (1968) 732–746.
  * [Du] Dumas SH, The KAM story (World Sci, 2014).
  * [F+] Faustin JM, Cooper WA, Graves JP, Pfefferlé D, Geiger J, Fast particle loss channels in Wendelstein 7-X, Nucl Fusion 56 (2016) 092006.
  * [Fe] Fenichel N, Persistence and smoothness of invariant manifolds for flows, Indiana U Math J 21 (1971) 193–225.
  * [GB] Garren DA, Boozer AH, Magnetic field strength of toroidal plasma equilibria, Phys Fluids B 3 (1991) 2805–21.
  * [G+] Garren AA, Riddell RJ, Smith L, Bing G, Henrich LR, Northrop TG, Roberts JE, Individual particle motion and the effect of scattering in an axially symmetric magnetic field, Univ of California Radiation Laboratory report UCRL-8076 (1958); permalink https://escholarship.org/uc/item/9rb7q30n
  * [GRR] Golab AJ, Robinson JC, Rodrigo JL, On 2D harmonic extensions of vector fields and stellarator coils. arXiv:2108.07643 (2021).
  * [GT] Goldston RJ, Towner HH, Effects of toroidal field ripple on suprathermal ions in tokamak plasmas, J Plasma Phys 26 (1981) 283–307.
  * [HM] Hall LS, McNamara B, Three-dimensional equilibrium of the anisotropic, finite-pressure guiding-center plasma: Theory of the magnetic plasma, Phys Fluid 18 (1975) 552–65.
  * [He] Helander P, Theory of plasma confinement in non-axisymmetric magnetic fields, Rep Prog Phys 77 (2014) 087001.
  * [HM82a] Holmes, P. & Marsden, J. Horseshoes in perturbations of Hamiltonian systems with two degrees of freedom. Commun. Math. Phys., 82, (1982), 523–544.
  * [HM82b] Holmes, P. & Marsden, J. Melnikov’s method and Arnold diffusion for perturbations of integrable Hamiltonian systems. J Math Phys 23,4, (1982), 669–675.
  * [HPS] Hirsch MW, Pugh CC, Shub M, Invariant manifolds, Lect Notes Math 583 (Springer, 1977).
  * [KIVM] Kallinikos N, Isliker H, Vlahos L, Meletlidou E, Integrable perturbed magnetic fields in toroidal geometry: An exact analytical flux surface label for large aspect ratio, Phys Plasma 21 (2014) 064504.
  * [KH] Katok A, Hasselblatt B, Introduction to the modern theory of dynamical systems (Cambridge U Press, 1995).
  * [Ke] Keller JB, Inverse problems, Am Math Monthly 83 (1976) 107–18.
  * [KW] Krajnak V, Waalkens H, The phase space geometry underlying roaming reaction dynamics, J Math Chem 56 (2018) 2341–78.
  * [Kr] Kruskal MD, Asymptotic theory of Hamiltonian and other systems with all solutions nearly periodic, J Math Phys 3 (1962) 806.
  * [Ku] Kuehn C, Multiple time scale dynamics (Springer, 2015).
  * [LC] Landreman M, Catto PJ, Omnigenity as generalized quasisymmetry, Phys Plasma 19 (2012) 056103.
  * [LP] Landreman M, Paul E, Magnetic fields with precise quasisymmetry for plasma confinement, Phys Rev Lett 128 (2022) 035001.
  * [La] Lang S, Introduction to Differentiable Manifolds, 3rd edn (New York: Interscience, 1967)
  * [Li] Littlejohn RG, Variational principles of guiding centre motion, J Plasma Phys 29 (1983) 111–25.
  * [M90] MacKay RS, Flux over a saddle, Phys Lett A 145 (1990) 425–7.
  * [M94] MacKay RS, On the motion of guiding centres, in: Transport, chaos and plasma physics, eds Benkadda S, Doveil F, Elskens Y (World Sci, 1994) 96–101.
  * [M04] MacKay RS, Slow manifolds, in: T Dauxois, A Litvak-Hinenzon, RS MacKay, A Spanoudaki (eds), Energy localisation and transfer (World Sci, 2004) 149–92.
  * [M20] MacKay RS, Differential forms for plasma physics, J Plasma Phys 86 (2020) 925860101
  * [MM] MacKay RS, Meiss JD, Flux and differences of action for continuous-time Hamiltonian systems, J Phys A 19 (1986) L255–9.
  * [MR] Marsden JE, Ratiu TS, Introduction to mechanics and symmetry, 2nd ed (Springer, 1999).
  * [Mat] Mather JN, Non-existence of invariant circles, Ergod Th Dyn Sys 4 (1984) 301–9.
  * [Men] Menyuk CR, Particle motion in the field of a modulated wave, Phys Rev A 31 (1985) 3282–90.
  * [Mey] Meyer KR, Generic bifurcation of periodic points, Trans Am Math Soc 149 (1970) 95–107.
  * [M+] Mikhailov MI, Cooper WA, Isaev MYu, Shafranov VD, Skovoroda AA, Subbotin AA, Improved stellarator systems, in ISPP-18 “Piero Caldirola”, Theory of Fusion plasmas, eds Connor JW, Sindoni E, Vaclavik J (SIF, Bologna, 1999) 185–198.
  * [Mi] Milnor J, Morse theory (Princeton U Press, 1969).
  * [Mo] Moser JK, Stable and random motions in dynamical systems (Princeton U Press, 1973)
  * [MCB] Mynick H, Chu, Boozer A, Class of Model Stellarator Fields with Enhanced Confinement, Phys Rev Lett 48 (1982) 322–6.
  * [Ne] Neishtadt AI, On the change in the adiabatic invariant on crossing a separatrix in systems with two degrees of freedom, PMM USSR 51:5 (1987) 586–92.
  * [PCHL] Parra FI, Calvo I, Helander P, Landreman M, Less constrained omnigeneous stellarators, Nuclear Fusion 55 (2015) 033005.
  * [P+] Paul EJ, Bhattacharjee A, Landreman M, Velasco JL, Energetic particle loss mechanisms in reactor-scale equilibria close to quasi-symmetry, Nuclear Fusion 62 (12), (2022), 126054.
  * [Poi] Poincaré H, Sur le problème des trois corps et les équations de la dynamique, Acta Math 13 (1890) 1–270.
  * [Pos] Post RF, The magnetic mirror approach to fusion, Nucl Fusion 27 (1987) 1579–739.
  * [RR] Renardy M, Rogers RC, An introduction to partial differential equations, 2nd ed (Springer, 2004).
  * [Ro] Robinson RC, Horseshoes for autonomous Hamiltonian systems using the Melnikov integral, Ergod Th Dyn Sys 8* (1988) 395–409.
  * [RBNV] Ross, S., BozorgMagham, A., Naik, S. & Virgin, L. Experimental validation of phase space conduits of transition between potential wells. Phys. Rev. E. 98 (5), (2018), 052214.
  * [SZ] Salamon D, Zehnder E, KAM theory in configuration space, Comment Math Helv 64 (1989) 84–132.
  * [Sk] Skovoroda AA, Pseudosymmetry near a magnetic surface in a plasma confinement system, Plasma Phys Reports 26 (2000) 550–9.
  * [SS] Skovoroda AA, Shafranov VD, Isometric magnetic confinement systems, Plasma Phys Rep 21 (1995) 886–906.
  * [So] Solov’ev LS, The theory of hydromagnetic stability of toroidal plasma configurations, Sov Phys JETP 26 (1968) 400–7.
  * [We] Weinstein A, Connections of Berry and Hannay type for moving Lagrangian submanifolds, Adv Math 82 (1990) 133-59.