Regions without flux surfaces of given class for magnetic fields in toroidal geometry

A converse KAM method for 3D vector fields, establishing regions through which passes no invariant 2-tori transverse to a given direction field, is tested on some helical perturbations of an axisymmetric magnetic field in toroidal geometry. It finds regions corresponding to magnetic islands and chaos for the fieldline flow. The minimization of these regions is proposed as a tool to help in the design of plasma confinement devices of tokamak and stellarator type.


Introduction
KAM theory provides sufficient conditions for the existence of invariant tori in Hamiltonian systems.In particular, many invariant tori persist from generic integrable Hamiltonian systems under smooth and small enough perturbations (for a semi-popular introduction, see [1]).Nevertheless, it is still hard work to prove existence of a realistic fraction of the tori that are suggested to exist by numerical simulation, e.g.[2].
An alternative approach is to determine regions through which no invariant tori of given class pass.Termed Converse KAM theory [3,4], it is much easier to implement than KAM theory and gives close to optimal conclusions without excessive computation.
The present work is an application to magnetic fields of Converse KAM theory, as extended in [5] to allow more general classes of tori than the earlier references and to treat 3D vector fields rather than positive-definite Lagrangian systems.It follows the main points of the implementation presented in [6], which was for the planar circular restricted three-body problem on level sets of the Jacobi constant.
The case of study in this work is magnetic fields in toroidal configurations, in particular the identification of regions through which pass no invariant tori (flux surfaces) of a given class.We define a class of tori by specifying a direction field almost everywhere and asking for tori that are transverse to that direction field.The principal choice of direction field is the gradient of a suitable notion of distance from a closed fieldline (magnetic axis), with respect to a chosen metric.
Our method uses the magnetic flux-form as principal representation of a magnetic field.Given a vector field B preserving a volume-form Ω, its flux-form is β = i B Ω. The integral S β over any surface S (with boundary allowed) represents the magnetic flux through that surface.Appendix A presents a summary of relevant background.For a tutorial about the use of differential forms in plasma physics, see [7].
In Section 2, we introduce the magnetic fields to be studied in this paper.In Section 3, we explain how to apply the Converse KAM method to magnetic fields.In Section 4, we present the results of numerical implementation of the method on the chosen fields.Section 5 discusses the results.Finally, three appendices give pedagogical introductions to some of the mathematics.

Toroidally helical magnetic fields
The magnetic fields that we choose to illustrate the Converse KAM method here are perturbations of a circular tokamak field by helical modes, based on [8].They have the advantages that: 1. there is an explicit magnetic axis and an easily specified class of tori that surround it; 2. with a single helical mode, the field is still integrable, but has a computable island; the invariant tori outside the island all belong to the chosen class and none of those inside the island do, so the method can be tested on its ability to detect the island; 3. with more than one helical mode, the field can be expected to have the typical mix of invariant tori of the original class, islands, and chaos, so the method can be tested on such cases; 4. they show how to handle fields presented in non-trivial coordinate systems, which is the typical case for tokamak and stellarator fields.
Given a coordinate system (x 1 , x 2 , x 3 ), we will express a magnetic field B in terms of its contravariant components B i , rather than its physical ones; they differ by length factors (see Appendix B for a summary about components of vector fields in curvilinear coordinates).An advantage is that the equations of motion for fieldline flow are just ẋi = B i (x); here "time" is to be understood along the magnetic field lines.
Our fields are simplest described and treated in an adapted toroidal coordinate system (ψ, ϑ, φ), which is a variant of the standard toroidal coordinates (r, θ, φ).It is not straightforward to describe the coordinate system, and it needs first considerations on B, in particular its toroidal component B φ , but we do it in the next subsection.It might seem demanding, but working in a non-trivial coordinate system is likely to be part of any application to realistic fields.
First we recall the standard toroidal coordinates (r, θ, φ).They are related to Cartesian coordinates (x, y, z) through for some R 0 > 0 and 0 ≤ r < R 0 .As shown in Figure 1, R 0 is the radius of the magnetic axis and R represents the cylindrical radius relative to the z-axis.In these coordinates, the metric tensor is represented by the matrix diag(1, r 2 , R 2 ).

Construction of adapted coordinates
Following [8], we introduce coordinates (ψ, ϑ) to make the restriction β T of the magnetic flux-form β to a poloidal section (φ = constant) take the form where ∧ denotes the exterior product of differential forms (see [7] for a tutorial).First, we define ψ as the toroidal magnetic flux across the poloidal disk of radius r about a point of the magnetic axis, divided by 2π.Thus, integrating over a poloidal disk of radius r yields ψ.Then ϑ can be constructed by equating (1) and (2), i.e., the condition rRB φ dr ∧ dθ = dψ ∧ dϑ.
Therefore the transformation (r, θ) −→ (ψ, ϑ) basically relies on the toroidal component B φ .We choose our magnetic fields to all have with B 0 > 0. This is a simple but realistic form for B φ , corresponding to external poloidal current 2πR 0 B 0 /µ 0 .Hence we arrive at [8,9] ( Note that ψ ∼ B 0 r 2 /2 as r/R 0 −→ 0, ψ is restricted to non-negative values less than B 0 R 2 0 , and has a coordinate singularity at 0.
The magnetic flux-form β plays a key role in the application of the Converse KAM method, so it is useful to simplify its expression by suitable coordinates.Here, in particular, we focused on the restriction β T of β to poloidal sections, because we will see that for our examples we can implement the method using only β T .For more general implementation though, it is essential to use the full flux-form β.
Remark 2.1 Another way of thinking of (ψ, ϑ) and β T is related to the standard Hamiltonian formulation of magnetic fields.The latter typically uses φ as time along the field lines, assuming B φ = 0. Field line flow can be written then as a time-dependent Hamiltonian system with Hamiltonian function H = − A φ (A being the vector potential, see below) and symplectic form ω = β T .Thus, bringing (2) to the form (1) simply amounts to finding canonical coordinates (ψ, ϑ) for ω.The Hamiltonian treatment though is neither necessary nor simpler for either Converse KAM or the helical fields we use.The flux-form β is key instead.

Magnetic fields studied
To enforce volume-preservation by the fieldline flow, we specify B as the curl of a vector potential A. In terms of the covariant components of A, the contravariant components of B are given by where is the Levi-Civita symbol and |g| is the determinant of the matrix g representing the metric tensor, ds 2 = g ij dx i dx j .In our adapted toroidal coordinates, the volume factor |g| is 1/B φ = R 2 /(B 0 R 0 ).This can be shown without finding g, by computing the volume-form Ω = dx ∧ dy ∧ dz = rR dr ∧ dθ ∧ dφ and using (3)-( 4) for the toroidal flux.We take a vector potential with helical modes introduced in its toroidal component, of the form (in covariant components) where w 1 ∈ R, w 2 = 0, m, n are integers with m ≥ 2, f mn are smooth functions and ζ mn arbitrary phases.The factor ψ m/2 is to make the resulting vector potential smooth at ψ = 0 (this was not done in [8] but its authors were interested there only in the neighbourhoods of islands).The coefficient w 2 produces shear.Extension to examples with change of sign of shear could be achieved by adding a term w 3 ψ 3 ; that would be a good next test case for the method, because in contrast to [4] the method here does not require shear, but we leave it for future work.
The vector potential (6) gives rise to the magnetic field B = (B 0 R 0 /R 2 )V where the components of the auxiliary vector field V are The cylindrical radius R occurring in the conversion from V to B can be expressed in our adapted coordinates via but we can avoid the conversion by applying the Converse KAM method to V rather than B, as will be explained.
Because we take m ≥ 2, the fields all have ψ = 0 as a closed fieldline, as claimed, which we call the magnetic axis.
We define the principal class of tori to be the differentiable tori that are transverse to ∇ψ.For example, with no helical modes the field is integrable with integral ψ and the invariant tori ψ = constant belong to the principal class.So do all C 1 -small deformations of them.Specifying ∇ψ entails a choice of Riemannian metric, but there is no need to use the Euclidean one, especially as in the adapted toroidal coordinates its computation would add extra work.It is preferable to choose a metric so that ∇ψ is in the same direction as ∂ ψ (see Appendix B for the distinction; in particular, this fails for the Euclidean metric: the adapted toroidal coordinates are not orthogonal).Then the principal class of tori consists of the graphs of ψ as a differentiable function of (ϑ, φ).We keep the more general specification ∇ψ, however, for flexibility.
The last ingredient to describe is the full magnetic flux-form β (as opposed to just its restriction β T to poloidal sections).This is defined by β = i B Ω where Ω is the volume form, or equivalently by β = dA , where (indeed, it is better to think of the vector potential A as a 1-form potential A for β).Thus Because V φ = 1, we see that restricted to a poloidal section, β = β T = dψ ∧ dϑ, as claimed earlier.

Single helical mode
A nice feature of our chosen form of field is that with a single helical mode, the field is still integrable [8].Indeed, it has the invariant (i.e., integral of motion) This can be checked directly.Alternatively, it can be derived from the symmetry u = n∂ ϑ + m∂ φ , as follows.The vector field u preserves the components A ψ , A ϑ , A φ in (6), and therefore A , i.e., L u A = 0. Thus, the rate of change of u • A along B is This result holds not only for fields with a single helical mode but also for any field with potential (6) in which A φ is a function of only ψ and a single combination mϑ − nφ of the angle variables.Note that although β is invariant under u, since L u β = dL u A = 0, the magnetic field B itself is not, as u is not volume-preserving.
In general, the integral gives rise to a family of invariant tori of the principal class and (if the signs are appropriate) a family that foliate an island.See Figure 2 for an example on a poloidal section.
To help orient the reader, we plotted this figure first in Cartesian coordinates, but in future we will just plot in "symplectic" coordinates on the poloidal section φ = 0. Near the magnetic axis this is a small distortion (especially for large aspect ratio r/R 0 −→ 0) of the true yz-plane x = 0, but with area equal to toroidal flux and the magnetic axis shifted to the origin.
In particular, for m = 2 and f (ψ) = f 0 + f 1 ψ, completing the square shows that the tori are the components of the sets where with ζ = 2ϑ−nφ+ζ 2n .These are graphs of ψ as a function of (ϑ, φ) if and only if the righthand side is everywhere positive.Supposing w 2 > 0, w 1 < n/2, (n/2 − w 1 )f 1 + 2w 2 f 0 = 0 and ε > 0 is small enough, we obtain the island explicitly as the set where The formula for other combinations of signs can be obtained if desired, but if |(n/2 − w 1 )f 1 + 2w 2 f 0 | < ε there is a more complicated island with four critical points.Note that the tori outside the island might not all be transverse to ∇ψ if the metric is not well chosen and ε is not small, but they are transverse to ∂ ψ .This explains our preference for a metric such that ∇ψ is in the direction of ∂ ψ .

Converse KAM for magnetic fields
The basic idea of Converse KAM methods is to consider how infinitesimal displacements ("tangent vectors" in mathematical terminology) rotate under a flow.If a 3D flow has an invariant orientable surface of a given class then it prevents infinitesimal displacements from rotating from one side of it to the other.So if an infinitesimal displacement from some trajectory rotates incompatibly with this restriction then there is no invariant surface containing that trajectory.
To make this precise, we have to specify a class of surfaces, in particular tori, and make clear what qualifies as rotating from one side to the other for all candidate tori in this class.We achieve these by choosing a "direction field" ξ and a 1-form λ, both to be explained below, and using the magnetic flux-form β.In the special case of fields with "stellarator symmetry" we explain how to streamline the method for symmetric trajectories.
Converse KAM theory in continuous time was first developed for Hamiltonian systems.It is standard knowledge that magnetic fields can be regarded as Hamiltonian systems (Appendix C describes the way we prefer to do this), but it is more straightforward to work directly with the magnetic flux-form.

Direction field
For a vector field B on an oriented 3D space, the Converse KAM method of [5] eliminates regions through which pass no invariant tori of B transverse to a given 1D foliation.A continuous choice of orientation can be assigned to the leaves of the foliation and thus a continuous choice of non-zero vectors ξ tangent to the foliation can be made, indicating the orientation.Because only the direction matters, not the magnitude, we call ξ a direction field (in standard differential-geometric terminology, ξ is the distribution associated to the foliation).
As presented in the previous section, we choose direction field for the principal class of tori in our examples to be ∇ψ with respect to some metric (which need not be the Euclidean one).In practise, we chose the metric to make ξ be in the direction of ∂ ψ , so that we can be sure of the classification of tori for integrable fields with one helical mode.We present the method for general ∇ψ, however, for compatibility with [6] and potential applications to include island tori where we'd replace ψ by Ψ of equation ( 8), and to guiding-centre motion.
An important extension is required, however, to cater for classes of tori around a magnetic axis.Namely, we allow the direction field ξ to have zeroes.This is the case for ∇ψ on the magnetic axis, for example.Note that if a torus is transverse to ξ then a fortiori it does not intersect the zero-set of ξ.

Nonexistence condition
The method of [5] gives a sufficient condition for the non-existence of invariant tori of a 3D vector field through a given point, transverse to a direction field ξ.We describe here its adaptation to magnetic fields B. We will assume B is nowhere zero in the domain of interest; equivalently, the kernel of β is one-dimensional at every point.
Before we start, for any positive function f , the vector field V = B/f has the same invariant tori as B. So it is a good idea to choose a function f to simplify the expression of V .See the previous section, for example.In general, V no longer preserves the same volumeform Ω as B but it preserves the related volume-form f Ω.Also the important relation i B β = 0 is inherited by V : i V β = 0. We will treat B in what follows, but one should bear in mind this possibly useful pre-processing.
Given an initial point s 0 = s(0) in 3D space, take initial tangent vector η s 0 = ξ s 0 .For t positive or negative, compute the propagation of s(t) and η s(t) under the dynamics ṡ = B(s) and the linearised dynamics ηs = DB s η s , respectively.If there is an invariant torus T passing though s 0 that is transverse to ξ, then η s(t) must stay on the same side of T for all t.In particular, η s is never of the form c 1 ξ s + c 2 B s with c 1 < 0. Checking this condition can be broken down into two steps: (i) examine if (η s , ξ s , B s ) pass through a case of linear dependence; (ii) if so, examine the sign of c 1 .
If one finds a time at which the stated vectors are linearly dependent with c 1 < 0 then the given trajectory does not lie on any invariant torus transverse to the given field ξ.
The conditions (i) and (ii) are stopping criteria for the integration of s(t) and η s(t) .To detect them, we follow the "general" formulation of [6].In the present context, this uses the magnetic flux-form β for (i) and a 1-form λ for (ii) such that λ(B) = 0 and λ(ξ) > 0 (except on zeroes of ξ).The reason that β suffices here (instead of the symplectic form on energy levels used in [6]) is that B belongs to the kernel of β, hence the triple product of (η s , ξ s , B s ) using the standard volume-form Ω = |B| −2 B ∧ β reduces to β(η s , ξ s ).Therefore, the general formulation translated to magnetic fields says that Theorem 3.1 Given initial conditions s 0 = s(0), choose initial vector η s 0 = ξ s 0 .If there is a time t such that at s = s(t): (i) β(η s , ξ s ) changes sign, and (ii) λ(η s ) < 0, then there is no invariant torus through s 0 transverse to ξ.
Thus, we can mark s 0 (and indeed its whole forward and backward orbit) as a point of the region of nonexistence of the desired class of invariant tori.In other words, we can eliminate s 0 (and its orbit) from being on invariant tori of the given class.
To apply this theorem, one needs to choose a 1-form λ with the required properties.There is some freedom here.In the case of ξ = ∇ψ for some choice of metric, following the steps from [6], we choose λ = dψ − kB , where k = B • ∇ψ/|B| 2 = B ψ /|B| 2 and both B and |B| 2 are defined using the chosen metric (B is the 1-form such that for all vectors u, B u = B • u).By construction, λ(B) = 0 and λ(ξ) > 0 everywhere except where ξ is parallel to B (by the Cauchy-Schwarz inequality).In our case, B φ > 0 implies that the only places where ξ is parallel to B are where ξ = 0, i.e., on the magnetic axis.
For our examples, it is simplest if the chosen metric is diagonal in adapted toroidal coordinates.It is best also if the metric makes dimensional sense.Thus we use which approximates Euclidean metric near the magnetic axis.
In vector calculus notation, the two quantities in Theorem 3.1 are expressed as β(η, ξ) = ξ • βη and λ(η) = η • (∇ψ − kB) for the aforementioned choice for λ, considering β as a matrix (as described in Appendix B), and • denoting the dot product with respect to the metric used.Also, note that η • ∇ψ = η ψ , independently of the metric.
Note that if we use V instead of B for our examples then V φ = 1 implies that an initially poloidal vector remains poloidal.Assuming ξ is chosen to be poloidal, this means firstly, the integration of tangent orbits does not require the φ-component to be represented, and secondly, β needs evaluating only on pairs of poloidal vectors, where it has the simple form dψ ∧ dϑ.We applied the method to both B and V , the first to demonstrate its general applicability, the second to speed up computations.

Using stellarator symmetry
The method of the previous subsection has a refinement for systems that admit a timereversal symmetry, as noted and used in [4].In the present context of magnetic field lines, the equivalent is the "stellarator symmetry" R : (r, θ, φ) −→ (r, −θ, −φ) [10], that translates to (ψ, ϑ, φ) → (ψ, −ϑ, −φ).The flow of a magnetic field B has stellarator symmetry if B is R-antisymmetric, i.e., BR = − dR B, or equivalently Bs = −B s , where B = dR B and s = R(s).Although this hypothesis narrows the magnetic fields that can be considered, it is a commonly assumed property in the fusion plasma literature and the design of stellarators.If the phases ζ mn are chosen zero then our examples have stellarator symmetry.
For magnetic fields with stellarator symmetry and initial conditions on a symmetry line (the half-lines of fixed points of R, i.e. θ = 0 or π, φ = 0 or π, r > 0) it is possible to simplify the test and speed up the computation by a factor of at least two, as the backward trajectory is a reflection by R of the forward one.Only now, we need to choose ξ to be R-symmetric, and to take as initial condition an R-antisymmetric vector η s 0 (not ξ s 0 as in the general formulation) on the symmetry semi-line, independent of B s 0 .
The two steps of the non-existence condition are then reduced to one, namely β(η s , ξ s ) = 0 for some t > 0. This is because if η s 0 evolves to η s(t) for some t > 0 with β(η s , ξ s ) = 0, then η s = c 1 ξ s + c 2 B s for some c 1 , c 2 , and c 1 = 0 since η was independent of B and independence is preserved by the evolution.But, by reflection, η s 0 also evolves backwards in time to ηs = −c 1 ξs + c 2 Bs .As previously explained, the change in sign of the component along ξ is incompatible with existence of an invariant torus through s 0 transverse to ξ. Theorem 3.2 Let B and ξ be R-antisymmetric and R-symmetric, respectively.Given initial conditions s 0 = s(0) on the symmetry lines, choose an initial R-antisymmetric vector η s 0 .If there is a time t such that β(η s , ξ s ) changes sign at a point s = s(t), then there is no invariant torus through s 0 = s(0) transverse to ξ.
To fix ideas, we choose ξ = ∂ ψ (equivalently, ψ∂ ψ to make it have a limit on the magnetic axis, as its magnitude doesn't matter) and η s 0 = ∂ ϑ .Then they both lie in a poloidal plane and so under the dynamics of V , η s remains in a poloidal plane and hence β can be simplified to dψ ∧ dϑ again.

Results
In this section, we apply the Converse KAM method laid out in Section 3 to integrable and non-integrable cases of magnetic fields of the type described in Section 2. In particular, we apply the method to find regions without invariant tori (that is, flux surfaces) transverse to the ψ-direction.Eliminated from being on such tori, they reveal magnetic islands and chaotic regions.
More specifically, Theorem 3.1 is applied to regular grids of initial conditions in symplectic coordinates (ỹ, z) (9) over the plane φ = 0 for the magnetic fields given by (7).The resolution of the grid is 160×160 initial conditions taken in each sample.Counting the initial conditions that are detected by the method allows to approximately bound from below the area (which in symplectic coordinates represents the toroidal flux) not occupied by tori transverse to the chosen direction.
Note that areas are the same when computed in symplectic coordinates (ỹ, z) or in (ψ, ϑ), because dỹ ∧ dz = dψ ∧ dϑ.Because of this, the area S of nonexistence in the plane φ = 0 can be estimated by counting the number of initial conditions detected by Theorem 3.1 on a regular grid over (ỹ, z, φ = 0).In other words, if S is the set of points detected by Theorem 3.1 on a N × N regular grid over [ỹ 0 − L, ỹ0 + L] × [z 0 − L, z0 + L], the area S is approximated by Also, Theorem 3.2 is applied to the same magnetic fields, however on a different set of initial conditions.As the method requires orbits starting from stellarator-symmetric lines, the initial conditions are taken uniformly in 2ψ/B 0 along the two semi-lines θ = 0, π (where z = 0) on the φ = 0 plane.In particular, it uses a regular partition of 200 points of the interval [−1, 1] in the ỹ-axis.We could also have taken initial conditions on the other two half-lines (θ = 0, π on φ = π), but for the choice of signs of ε mn that we use, we believe that φ = 0, θ = 0 is "dominant" in the sense that all the primary island chains have an elliptic point on it, and this tends to maximise the set of eliminated trajectories.
The figures that display the results of Theorem 3.2 are followed by Poincaré sections produced from the iteration of the selected initial points.If any point is detected for nonexistence then so is its whole trajectory; thus even though Theorem 3.2 is restricted to symmetric initial conditions, it has implications for a much larger set.However, estimating the areas occupied by the detected points from the results of this formulation is a more challenging problem, which we hope to address in the future.
To cater for the possibility that the termination condition is never reached on a trajectory, we choose a timeout t f .If the timeout is reached then the status of the chosen initial condition is undecided.This, of course, should include all initial conditions that are on invariant tori of the given class, but may include others for which more time would be required to detect the nonexistence.Depending on the implementation, the timeout values might not indicate how long were the trajectories.Thus, in the figures we also display the average of the last computed value of φ divided by 2π, i.e., the average number of toroidal laps.
For a trajectory, we denote by t * the time at which non-existence was detected or t f if it was not detected.As a measure of non-existence of tori of a given class, figures show in hues the relative time of detection, using the ratio q = t * /t f of time of detection to timeout.
In all the examples throughout this section and elsewhere, we take the following values and function for the vector potential ( 6) As previously mentioned, the results in all the forthcoming figures are presented over the poloidal plane φ = 0. Unless stated otherwise, the timeout used is t f = 200, which amounts to ∼ 32 laps around the z-axis.
The results from the general formulation of Theorem 3.1 applied to this case are shown in Figure 3.The choice of m = 2 allows an analytical expression for the separatrix delimiting the island.It is given by (10) with the limiting value of Ψ from (11).From this, we obtain the width ∆ψ for the island as a function of ϑ at given φ: Recall that w 2 corresponds to shear, so if f 1 = 0 we see the familiar behaviour ∆ψ ∼ 2 εf 0 ( n 2 − w 1 ) sin(ζ/2)/w 2 for small perturbation ε.The area S I of the island(s) can be computed by numerical integration of ∆ψ dϑ.By our choice of coordinates, this is equal to its toroidal flux.
Using (12), we calculate the area S of the nonexistence region detected by Theorem 3.1.In Figure 4, we see it as a function S = S(t f ) of timeout t f .Recall that timeout units correspond to t f /(2π) laps around the z-axis.As expected, the plot shows that the value of the estimated areas increase monotonically with t f up to a limiting value that agrees with the island area S I .
Note also that the time of first detection of the island can be predicted: it is the time for a tangent vector at the centre of the island to make one half of a poloidal revolution (this assumes that the rotation number in the island decreases as distance from the centre increases, else it would be detected earlier).In terms of Greene's residue R, which can be written as sin 2 (α/2) for eigenvalues e ±iα of the return map to a poloidal section, we see that for small R (i.e., approximating R = ( α 2 ) 2 ), the time t c of first detection should be asymptotically where T is the time for one toroidal revolution.
Figure 5 shows the numerically computed residue for the centre (and the x-point) of the island as a function of ε.Using the V -field, the period of the island centre is 4π, thus π 2 T ≈ 19.74, and comparing with Figure 4, we see that the formula (15) gives a reasonable prediction of the first time of detection of the island.Furthermore, for this example, we see from Figure 4 that more than half the area of the island has been detected within twice the time of first detection.
The area S of the nonexistence region detected by Theorem 3.1 as a function S = S(ε) of the perturbation parameter ε is presented in Figure 6 for different values of timeout t f .It shows that the error of the estimation is reduced for small values of ε as t f increases, but ultimately deviates for larger ε.The most likely explanation of this behaviour is that our regular grid does not have enough points in the magnetic island to give a reliable estimation for this region of the parameter.
Figure 7 shows the Converse KAM results now using Theorem 3.2.Since this theorem can only be applied on symmetric semi-lines, the results cannot be directly used to estimate the nonexistence area.The left picture shows the relative times of detection in terms of q = t * /t f for initial conditions on the partition of the semi-lines.The picture on the right shows the Poincaré plot obtained from iteration of the selected initial conditions for the given timeout t f .Compared to the one in Figure 3, it is worth noting that the time of first detection of the island is now half as much, as expected.

Non-integrable examples
Next we consider magnetic fields with more than one helical term, derived from (6).The Poincaré section in all the forthcoming examples displays features of typical near-integrable systems: tori of different classes and chaotic regions near the hyperbolic saddle of the resonances (magnetic islands).Both formulations, Theorems 3.1 and 3.2, yield closely aligned results.Using a radial direction field, they are able to identify and eliminate points (and in fact whole field lines) that do not lie on tori of the original class.They do not distinguish, however, between the ones lying on tori of another class or in chaotic regions.But if required, the use of a suitable foliation centered on the elliptic field lines of an island chain could differentiate between those two cases.

Example 2
The second example corresponds to the magnetic field derived from ( 6) for two modes now, namely the resonances 2/1 and 3/2 with same perturbation parameter value Following the same order as in previous example, the Converse KAM results using the formulation of Theorem 3.1 are shown in Figures 8-10.
Figure 8 shows the detection in symplectic coordinates, using the same color scheme as in Figure 3.As we can see on the right, the different q-hues suggest the location of the two magnetic islands corresponding to this example.Figure 9 shows the computed area S = S(t f ) of nonexistence from Theorem 3.1 for the present example for different values of the perturbation parameter ε.As in Figure 4, we see that the estimated areas increase monotonically with t f and seem to be approaching a limit.Figure 10 shows the estimated area S = S(ε) of nonexistence from Theorem 3.1 now as a function of the perturbation parameter ε for different values of timeout t f .The behaviour seems to be not as simple as in Figure 6 for the case of the single resonance.A possible explanation of this may be the interaction between the resonances as they grow with ε.The Converse KAM results using Theorem 3.2 are shown in Figure 11.The left plot, compared to the one in Figure 7, shows an asymmetrical distribution of the relative time of detection q, which is consistent with the resonances used in this example.

Example 3
The next example corresponds to the magnetic field obtained from (6) for the resonances 2/1 again and 5/4 now, with fixed value ε 21 = 0.001 and varying value of ε 54 = ε.That is, The Converse KAM results using the formulation of Theorem 3.1 are shown in Figures 12-14.Figure 12 shows the detection in symplectic coordinates for ε = 0.01, Figure 13 shows the computed area S = S(t f ) of nonexistence for the present example for different values of the perturbation parameter ε, and Figure 14 shows S = S(ε) for different values of timeout t f .The results behave as expected, except for the particularity that the islands seem to be detected at relatively different times.This is noticeable in the right plot of Figure 12, where the resonance 5/4 is seen mostly in the orange area, while the 2/1 has a green hue instead.Larger timeout t f is required for the method to detect magnetic islands of small amplitude (as quantified by the residue).
The corresponding Converse KAM results of Theorem 3.2 are shown in Figure 15.The results are similar to Figure 11 for Example 2, besides the different relative detection time q for each magnetic island that we see again.

Conclusions
The paper reports on numerical implementation of the Converse KAM method from [5] on some example magnetic fields.It has demonstrated that the method allows one to identify many of the points that do not belong to any flux surface of a given class.It has been shown to reach decisions in relatively short times on these examples.In an example with an integrable island, it detects a large fraction of the island in a fieldline flow time of order π 2 T / √ R (in the symmetric formulation, or twice this for the general formulation), where T is the time for one revolution around the z-axis and R is Greene's residue for the island.For fields with stellarator symmetry, it suffices to examine initial conditions just on the symmetry lines.
The method can be used to compute a good lower bound for the toroidal flux that is not on flux surfaces of the desired class.This is suitable for passing as an objective function to include in optimisation of the design of stellarator fields.One could also compute a lower bound on the volume not on flux surfaces of given class, by integrating the return time with respect to the toroidal flux over the detected points on a poloidal section.
A crucial component of this method is the selection of a suitable direction field transverse to the class of tori of interest.Our examples have a natural one, which made them a simple test case.For a more general magnetic field, one would need to determine the magnetic axis and a suitable direction field from it, but the freedom to choose the direction field means that the method could in principle be used for fields with bean-shaped cross-section, as in W7X.Also, one could be interested in survival of tori both around the magnetic axis and in a major island.For example, for perturbations of an integrable field with one helical mode, denoting by Ψ the conserved quantity for the integrable case, one could use the direction field ∇Ψ, with respect to some metric.For some discussion about how to choose the direction field in other contexts, see [11,6].
Although the examples treated here are simple, they already include ones that display the typical mix of islands and chaos.A next goal is to report on applications of the method to fields produced by the stellarator optimisation code SIMSOPT.These are designed to be close to integrable, so the aim is to detect and quantify the remaining deviations from integrability, which requires an efficient method, as we believe is ours.
It would be good to implement also the extension of the method called "killends", which uses bounds on the slope of invariant tori of given class to extend the region through which they cannot pass [5].This can eliminate more points without computing more trajectories.Indeed, it can result in a saving on the total length of trajectories computed, because in essence the trajectories from a grid of points in a transverse section are computed for one revolution, whereas in the general formulation used here, the trajectories from a grid were computed until non-existence detected or timeout, which typically takes many revolutions.
The formulations presented can be extended to guiding-centre motion, which is a 2parameter family of 3D-systems, parametrised by magnetic moment µ and energy E. The 3D space is the set of (x, v) ∈ R 3 × R satisfying 1  2 mv 2 + µ|B(x)| = E, where v represents parallel velocity.The flux-form β is replaced by eβ + md(vb ), where b is the unit vector along the magnetic field (e and m are the charge and mass of the particle).One would have to choose appropriate classes of tori for guiding-centre motion.The issues will be somewhat similar to those for the planar circular restricted three-body problem, treated in [6].
In the magnetic field context, a similar method ("phase rotation") was introduced by White [12] (see also Fig 6.8 in [13]) and applied to guiding-centre motion in a magnetic field, but in our opinion it needs some clarification.Firstly, it needs stating that the class of tori under consideration are the graphs of functions P = P ζ (θ) in the given coordinate system.Secondly, it is stated that the angle χ for the displacement vector between orbits on different invariant tori cannot rotate by more than π, but should specify that this means that relative to the vertical, χ has to remain in (−π, π).Our method can be considered to be the limiting case from displacement to tangent vectors, which may be more effective because the displacement vector from a trajectory in an island to a trajectory with a different rotation number might not rotate by more than π whereas an infinitesimal one will.Moreover, our method does not require a Poincaré section, and extends to other classes of tori.White's method has the advantage relative to [4] that it does not require shear, though it was realised some time ago that at least the 3D case of [4] did not require shear (leading eventually to [5]).
The method can be applied to other plasma physics problems too.For example, for an interface in a stepped pressure equilibrium to support a pressure jump, one needs an invariant torus of the "pressure-jump" Hamiltonian.Some conditions under which none exist were determined by [14], but it would be useful to extend them using the Converse KAM method.This would, for example, shed light on the work of [15].The standard direction field is the relevant one for this problem, so it suffices to use [4] rather than the current paper.Note that that paper applied the Converse KAM method to another plasma physics problem: the motion of a charged particle in the field of two electrostatic waves.Further plasma physics examples suggested by a reviewer include time-dependent fields such as those arising from resonances between particles and Alfvén modes [16] and the quasi-single-helicity states reported in [17].at a point.The relation between the flux-form and the field is where Ω is the volume-form (non-degenerate 3-form) with respect to which B is divergencefree (divB = 0 if and only if β is closed: dβ = 0).The formula says how to obtain β from B: β(ξ, η) = Ω(B, ξ, η) = B • (ξ × η) in vector calculus notation.
Conversely, given a 2-form β and a volume-form Ω in 3D, one can obtain vector field B at any point by noting that the kernel of β (ker β = {u : i u β = 0}) is non-zero because β is antisymmetric and dimension 3 is not even.Take any non-zero u ∈ ker β and extend to a basis (u, v, w); let B = β(v,w) Ω(u,v,w) u.Then the resulting B does not depend on the choices of (u, v, w) and satisfies i B Ω = β.
Note that representation of magnetic field lines by a closed 2-form β does not require the volume-form Ω.The fieldlines are the integral curves of ker β.All the volume-form does is to define the speed at which B goes along them.

A.2 Vector potential
It is commonplace that for a 3D vector field B, div B = 0 iff there exists a vector field A such that B = curl A. We explain here the related result for 2-forms.
First, we note that the above statement is not quite true for general 3D manifolds.One should strengthen the definition of a magnetic field from div B = 0 to S B • dS = 0 for every closed surface S (this is often called "absence of magnetic monopoles").The corresponding strengthening of the definition of a magnetic flux form β is that S β = 0 for all closed surfaces S.
It follows that there exists a 1-form α such that β = dα.Such an α is called a potential for β.Given a vector potential A for B, a potential for β is α = A , where A (ξ) = A•ξ for all tangents ξ.Conversely, any potential α defines a vector potential A. Just as A is non-unique up to addition of any gradient (even multi-valued), α is non-unique up to addition of any closed 1-form.The relation β = dα is equivalent to B = curl A, but is simpler because it makes no use of a Riemannian metric.

A.3 Integrable fields
We say a magnetic field B is integrable if there is a function Ψ with non-zero derivative almost everywhere such that B • ∇Ψ = 0.It follows that the surfaces of constant Ψ are invariant under the fieldline flow.
Integrability is nicely addressed at the level of continuous symmetries of potentials for the flux-form.If there is a vector field u such that the Lie derivative L u α = df for some function ζ then using L u = i u d + di u on differential forms one obtains i u β = dΨ, with Ψ = f − i u α.This says that i u i B Ω = dΨ, so in particular i B dΨ = 0 by antisymmetry of Ω.In vector calculus these relations are B × u = ∇Ψ and B • ∇Ψ = 0.If u is independent of B almost everywhere then dΨ = 0 almost everywhere.speed along ker ω E .
Analogously, given a volume-form Ω in 3D and a closed 2-form β, a speed for B is determined by i B Ω = β.

Figure 3 :
Figure 3: Converse KAM results using Theorem 3.1 for Example 1 with ε = 0.004 in symplectic coordinates (9).On the left, red = nonexistence, blue = undetermined.On the right, hues vary from fast detection (red) to no detection at all (blue) within timeout.The white curves are level sets of the invariant Ψ.

Figure 4 :
Figure 4: Nonexistence area S(t f ) detected by Theorem 3.1 of Converse KAM for different values of ε for Example 1.The dashed lines represent the corresponding island areas S I .

Figure 5 :
Figure 5: Greene's residue of the island centre (red) and negative residue of the island x-point (blue) for Example 1, as functions of ε.

Figure 6 :
Figure 6: Nonexistence area S(ε) detected by Theorem 3.1 of Converse KAM for different values of t f for Example 1, compared with the area S I (ε) of the island.

Figure 7 :
Figure 7: Converse KAM results using Theorem 3.2 for Example 1 with ε = 0.004, over the symmetrical semi-lines ϑ = 0, π (z = 0).Hues vary from fast detection (red) to no detection at all (blue) within timeout; level sets of the invariant Ψ are superimposed (right).Relative time of detection as a function of symmetrical initial position (left).

Figure 8 :
Figure 8: Converse KAM results using Theorem 3.1 for Example 2 with ε = 0.003 in symplectic coordinates (9).On the left, red = nonexistence, blue = undetermined.On the right, hues vary from fast detection (red) to no detection at all (blue) within timeout.

Figure 9 :
Figure 9: Nonexistence area S(t f ) detected by Theorem 3.1 of Converse KAM for different values of ε for Example 2.

Figure 10 :
Figure 10: Nonexistence area S(ε) detected by Theorem 3.1 of Converse KAM for different values of t f for Example 2.

Figure 11 :
Figure 11: Converse KAM results using Theorem 3.2 for Example 2 with ε = 0.003, over the symmetrical semi-lines ϑ = 0, π (z = 0).Hues vary from fast detection (red) to no detection at all (blue) within timeout (right).Relative time of detection as a function of symmetrical initial position (left).

Figure 12 :
Figure 12: Converse KAM results using Theorem 3.1 for Example 3 with ε = 0.01 in symplectic coordinates (9).On the left, red = nonexistence, blue = undetermined.On the right, hues vary from fast detection (red) to no detection at all (blue) within timeout.

Figure 13 :
Figure 13: Nonexistence area S(t f ) detected by Theorem 3.1 of Converse KAM for different values of ε for Example 3.

Figure 14 :
Figure 14: Nonexistence area S(ε) detected by Theorem 3.1 of Converse KAM for different values of t f for Example 3.

Figure 15 :
Figure 15: Converse KAM results using Theorem 3.2 for Example 3 with ε = 0.01, over the symmetrical semi-lines ϑ = 0, π (z = 0).Hues vary from fast detection (red) to no detection at all (blue) within timeout (right).Relative time of detection as a function of symmetrical initial position (left).