Perturbing an axisymmetric magnetic equilibrium to obtain a quasi-axisymmetric stellarator

It is demonstrated that finite-pressure, approximately quasi-axisymmetric stellarator equilibria can be directly constructed (without numerical optimization) via perturbations of given axisymmetric equilibria. The size of such perturbations is measured in two ways, via the fractional external rotation and, alternatively, via the relative magnetic field strength, i.e. the average size of the perturbed magnetic field, divided by the unperturbed field strength. It is found that significant fractional external rotational transform can be generated by quasi-axisymmetric perturbations, with a similar value of the relative field strength, despite the fact that the former scales more weakly with the perturbation size. High mode number perturbations are identified as a candidate for generating such transform with local current distributions. Implications for the development of a general non-perturbative solver for optimal stellarator equilibria is discussed.


Introduction
Quasi-symmetry is a property of magnetic fields that ensures the confinement of collisionless particle orbits. Axisymmetric magnetic equilibria possess this property in a trivial sense, whereas the related class of stellarators, called quasi-axisymmetric (QAS), satisfy the symmetry in a way that is hidden to the naked eye [Nührenberg et al., 1994].
The close relationship between axisymmetry and QAS suggests that the second class may be continuously connected to the first, and in particular that QAS stellarators may be obtained by deformation of axisymmetric equilibria [Boozer, 2008]. It has also been suggested that modifying tokamak equilibria by non-axisymmetric shaping might help overcome the stability issues that plague them, and a previous study, using conventional numerical optimization, has demonstrated that suitable QAS may indeed be found as deformed tokamak equilibria [Ku and Boozer, 2009]. The idea of passively stabilizing a tokamak by non-axisymmetric perturbations is also supported by a number of experimental results, when the perturbation generates a sufficient "external" boost in rotational transform [W VIIA Team, 1980, Pandya et al., 2015.
Solving the MHD equilibrium problem for optimal stellarator equilibria, without the use of numerical optimization algorithms (i.e. "direct construction" of optimal solutions), is potentially beneficial due to the speedup offered . So far, the only ways to do this have involved approximations to the problem like small distance from the magnetic axis [Garren and Boozer, 1991a,b, Landreman and Sengupta, 2018, or small deviation from axi-symmetry [Plunk and Helander, 2018]. Solving these approximate problems can also lead to fundamental insights about the properties of equilibria, and the size of the solution space.
There are possible practical advantages of directly constructing QAS stellarator equilibria via perturbation of axisymmetric equilibria, as compared to conventional optimization. For instance, the general perturbation can be constructed as a sum of independent QAS modes with different toroidal mode numbers. After pre-computation of these modes, the corresponding space of QAS equilibria can be easily scanned, without further computational cost, whereas each step of a conventional optimizer involves solving the equilibrium problem anew. Also, there is no fundamental constraint on axisymmetric equilibrium measures, like aspect ratio, so these may be set arbitrarily to explore new areas of stellarator design space, which may have been inaccessible with conventional optimization.
In a previous paper [Plunk and Helander, 2018], it was proved that nearly axi-symmetric magnetic fields can be constructed to satisfy the condition of quasi-axisymmetry on a single magnetic surface. These solutions, however, apply only to vacuum conditions, where the plasma itself does not contribute significantly to the magnetic field. The present work considers the more general case of finite pressure equilibria. Formidable challenges are present in this general problem, starting with an increased complexity arising from the nonlinear coupling of multiple fields. The presence of singularities in the force balance equation makes the general problem of obtaining equilibria ill-posed, even without the requirement of satisfying a special symmetry. As we will show, the issue of force balance singularities may be overcome, at least at first order in the expansion, by suitable choice of the zeroth-order rotational transform profile. The complexity of the system, however, makes it more difficult to establish existence of solutions by the same methods employed by Plunk and Helander [2018]. We therefore turn to devising a method to numerically solve the system. This, as we find, gives evidence that the same problem as solved in the vacuum limit by Plunk and Helander [2018], namely the problem of finding a perturbation of specified toroidal mode number N that satisfies the condition of QAS on a single magnetic surface, is indeed well-posed, at least in some practical sense.
The contents of the paper are as follows. In section 2 the basic equations and notation are established, and the "inverse" MHD equilibrium problem formulation is described. In section 3, the expansion about axi-symmetry is performed, and the equations are given to find perturbations satisfying QAS on a specified magnetic surface. The issue of force balance singularities is discussed, and a strategy to overcome them is described. In section 4, a numerical method is described to solve the first order system, and a set of solutions are given, based on a zeroth order ITER-like equilibrium. The VMEC [Hirshman and Whitson, 1983] and BOOZ_XFORM [Sanchez et al., 2000] codes are used to demonstrate that the solutions can satisfy the appropriate level of QAS as predicted by the theory.

Preliminaries
The MHD equilibrium equations are We assume topologically toroidal magnetic surfaces, here labeled by the flux function ψ. To solve these equations, we use a similar approach as previous works Garren and Boozer [1991a,b], Hegna [2000], Boozer [2002], Weitzner [2014]. Boozer angles are denoted θ and ϕ. The contravariant form of B is written where ι is the rotational transform, and 2πψ is the toroidal flux. This form of B satisfies zero divergence, assuming flux-surface geometry. The covariant form is written This form is a consequence of j · ∇ψ = 0 (3), and Ampere's law (1); see e.g. Helander [2014]. The basic strategy to find an equilibrium is to assert B con = B cov together with force balance (3), relying on the fact that these forms of the magnetic field incorporate Eqns. 1 and 2 as well as the assumption of topologically toroidal magnetic surfaces. Either the magnetic coordinates ψ, θ, and ϕ can be considered as the unknown functions of spatial coordinates ("direct formulation"), or the coordinate mapping x(ψ, θ, ϕ) can be considered as the unknown function of magnetic coordinates ("inverse formulation"). Both formulations are used here.
It is convenient at zeroth order to solve the Grad-Shafranov equation (e.g. using the direct formulation). This means that we are able to specify the axisymmetric shape of the outer magnetic surface. We will also specify the current and pressure profiles at this stage, and consider them as fixed for the remainder of the calculation. We will use the indirect formulation for the problem at next order, i.e. the problem of QAS-preserving perturbations, as it casts the problem as a fixed boundary problem with QAS as the boundary condition.

Problem formulation
With the inverse formulation, the independent variables of the problem are the magnetic coordinates, and QAS is expressed as a simple constraint, ∂B/∂ϕ = 0. Instead of using magnetic flux as a coordinate, we will use a coordinate system based on a dimensionless radial coordinate ρ = ψ/ψ b , where ψ b denotes the value of ψ on the boundary surface. Note that most physical quantities are not analytic in ψ at the magnetic axis (ψ = 0), but can be expanded in ρ [Garren and Boozer, 1991a]. This idea is motivated by considering ρ and θ as polar coordinates and then assuming that a Taylor expansion can be made in the pseudo-cartesian coordinatesx = ρ cos(θ) andȳ = ρ sin(θ).
With the inverse formulation, the unknown of the theory is the coordinate mapping x(ρ, θ, ϕ), and the equilibrium equations are written in terms of various derivatives ∂x/∂ρ, and so forth. These equations can be translated into equations involving the metrics via the usual identities (reviewed in Appx. A). The equation B con = B cov becomes Force balance can be expressed as a scalar equation, since it only has a component in the ∇ρ direction. One uses j = µ −1 0 ∇ × B cov and take the scaler product of Eqn. 3 with ∇θ × ∇ϕ We introduce the following normalized quantities:Ḡ = G/(2ψ b ),Ī = I/(2ψ b ),K = ρK andp = p/(2ψ b ) 2 . Note that we defineJ = J/ρ in the limiting sense so that although J = (∇ρ·∇θ ×∇ϕ) −1 tends to zero with ρ,J does not. Finally, we will need the following expression for the regularized Jacobian:J Defining alsoB = B/(2ψ b ) we have the useful relationJ = (Ḡ + ι Ī )/B 2 (Eqn. 51) so that QAS can be expressed most conveniently here as In the present work, we look for solutions that satisfy this condition on a single magnetic surface; we will not consider here the question of whether this condition might, under special circumstances, be satisfied globally, i.e. uniformly in ρ.

The expansion about axisymmetry
We write the coordinate mapping x(ρ, θ, ϕ) as a series expansion in the small parameter , where x 0 corresponds to the zeroth-order axisymmetric equilibrium. We will consider the pressurep and currentsḠ andĪ as fixed to their zeroth order values (there is no loss of generality as any higher order variation in these functions can be absorbed into the zeroth order forms). This confines our attention to axisymmetry breaking perturbations. We must however allow the deformation to modify ι andK.
For a nearly axisymmetric equilibrium, it is sensible to take the components of Eqn. 20 along the cylindrical unit vectorsR,φ,ẑ (such thatR ×φ =ẑ).

Order 0
The zeroth order coordinate mapping is (see also appendix D) where the cylindrical unit vectors are functions of the geometric toroidal coordinate, related to the boozer angle by ϕ = φ + ν, and is expanded as so that φ 0 = ϕ − ν 0 and φ 1 = −ν 1 , etc., and, for simplicity, the unit vectors will be defined according to the zeroth order expression of the geometric toroidal angle, With these definitions, derivatives of the zeroth order coordinate mapping are evaluated as The zeroth order MHD constraint is where, Eqns. 17-19 can be substituted in and the equation projected along the unit vectorsR, φ andẑ to obtain three coupled equations. Note that we avoid explicitly writing the lengthy equations that result, and will do likewise with others that follow, especially when they do not give any useful insight. Force balance is

Inverting the Grad-Shafranov solution
It is convenient to use Grad-Shafranov (GS) theory to obtain the zeroth order equilibrium. This approach gives control of the axisymmetric plasma shape, and also benefits from existing understanding of the equation and its numerical solution. A solution of the GS equation is the poloidal flux function Ψ(R, z) is obtained from a given pressure function p, and the poloidal flux function G. From these quantities, the corresponding profiles I and ι 0 , the current potential K 0 and coordinate mapping components R 0 , Z 0 and ν 0 can be calculated. To perform the coordinate inversion, derivatives of the GS solution are computed, so a high degree of accuracy is needed. A method is described in appendix D.

O( 1 )
As in Plunk and Helander [2018], we do not modify the toroidal angle beyond zeroth order in the expansion (ν 1 = 0, etc., in Eqn. 14), but instead consider the corrections to the coordinate mapping to have a component in theφ direction, i.e.
from which it follows that As ϕ is an ignorable coordinate in the properly formulated first order equilibrium equations, we will assume with N = 0 an integer. The deformation is thus non-axisymmetric, and the axisymmetric (ϕaveraged) part of the local MHD constraint 58 is ι 1 (Ḡg 23 ) = 0, from which we conclude that The first order MHD constraint is then What is needed is the exp(iN ϕ) component of this equation, obtained by substituting Eqns. 26-29 into x 1 , Eqn. 22, and its derivatives, Eqns. 23-25. The further substitution of zeroth order expressions, Eqns. 17-19, and projection along the cylindrical unit vectors, then yields a set of three equations for the unknownsR 1 ,Ẑ 1 ,Φ 1 ,K 1 in terms of the known zeroth order solutions R 0 , Z 0 , ν 0 andK 0 . The system is completed with the force balance equation, The exp(iN ϕ) component of the first order Jacobian,Ĵ 1 , is obtained from Eqn. 8 by substituting Eqns 18-19 and Eqns 24-25, into the following expression: We note that QAS implies thatĴ 1 = 0, so force balance on any QAS surface reduces to iNK 1 + ι 0 ∂K 1 /∂θ = 0, which, assuming irrational ι 0 , implieŝ On magnetic surfaces where QAS is not satisfied, the possibility of resonances in Eqn. 32 must be considered. It is easy to see that the equation can be uniquely solved forK 1 , periodic in θ, if ι 0 is not equal to a rational number. Actually, some rational numbers are resonant, and some are not, in particular there are resonances at any magnetic surface where ι 0 satisfies for arbitrary integer m. One strategy to avoid resonances is to constrain ι 0 to lie between two neighboring singular values. In that case, force balance can be considered "soluble" throughout the plasma volume. Note that, assuming ι 0 ∼ 1, such "safe" ranges becomes increasingly narrow at large N , although resonances may be considered "high order" in this limit, and therefore less likely to pollute the solution.. Note that Eqn. 34 demonstrates that force balance is non-resonant on a QAS magnetic surface. Actually, this reflects a general non-perturbative property, which follows directly from exact force balance and the relationship between B andJ (Eqns. 7 and 51), namely that quasisymmetry drastically simplifies the source term in force balance (the Fourier series ofJ in θ and ϕ has only terms of a single helicity), so that an MHD equilibrium that is quasi-symmetric globally (on all magnetic surfaces) should be free from nontrivial resonances; see also Burby et al. [2019] and Rodríguez et al. [2020].

O( 2 )
Proceeding to the next order, the equations are quite similar as before, but now include terms that are quadratic in first order quantities. The second order coordinate mapping is and its derivatives are The appearance of the nonlinear terms occurs (in the MHD constraint and force balance equation) at toroidal mode numbers ±2N , and also in the axisymmetric component, which now must be solved to obtain ι 2 . We note that the appearance of toroidal mode numbers ±2N at second order implies a denser set of possible force balance resonances at higher orders in the expansion, i.e. ι 0 = 2N/m, which may justify further restriction on the chosen profile for ι 0 . Even if the problem will only be solved at first order, higher order resonances may occur in the exact force balance equation that must be satisfied by the full equilibrium.
In the vacuum case [Plunk and Helander, 2018], ι 2 was obtained as a solubility constraint of the axisymmetric component (ϕ-average) of the local MHD constraint, Eqn. 58. This result does not appear to generalize in a simple way, implying that the full system (MHD constraint plus force balance) must be solved to obtain ι 2 .

Numerical Solution
The task now is to solve the system composed of Eqn. 31 and Eqn. 32 for the unknownsR 1 , Z 1 ,Φ 1 ,K 1 , subject to QAS (K = 0) on a specified magnetic surface, typically the outermost magnetic surface (ρ = 1). Note that this boundary surface need not necessarily be taken to be located at the plasma edge. It has been suggested [Henneberg et al., 2019] that it may be optimal to satisfy QAS at some intermediate magnetic surface, which can be implemented here by redefinition of the coordinate ρ, or simply choosing a boundary ρ < 1. Henceforth we assume ρ = 1 for simplicity.
The "pseudo-cartesian" coordinatesx = ρ cos(θ) andȳ = ρ sin(θ) are used for numerical purposes, instead of the polar coordinates ρ and θ. These have the advantage that they do not possess the singularity of the polar coordinates (ρ, θ) as ρ → 0, and they do not require periodicity to be enforced in θ, or any analyticity at ρ = 0. The only advantage found in using ρθ coordinates is to explicitly observe development of non-analyticity in the solutions on resonant surfaces (satisfying ι 0 = N/m).
The finite element method is used, reformulating the problem as an equivalent "least squares" problem. The least squares finite element method offers better convergence and stability properties for systems of first order PDEs [Jiang, 1998]. To show how the problem is reformulated, we introduce the vector field u(x,ȳ) = [K 1 ,R 1 ,Ẑ 1 ,Φ 1 ] T , together with the inner product where v * denotes the complex conjugate of v, and the integral is performed over the computational domain, the unit disk Ω. The original first order system of equations (i.e. theR,φ, and z components of Eqn. 31, coupled with Eqn. 32) can be written as Lu = 0, where where u j denotes the jth component of u, ∂ 1 and ∂ 2 denote ∂/∂x and ∂/∂ȳ, respectively, and the tensors a ij and α ijk encode the coefficients of the system of equations. The adjoint of L is denoted as L † , and is given by With these definitions, our problem is transformed into solving the following eigenvalue problem subject to the QAS boundary condition u 1 = 0 on the boundary ∂Ω, for solutions with eigenvalue λ → 0. This system is generated by computer algebra, and not explicitly written down, due to its complexity.

Examples
To demonstrate that the above numerical method works, in practice, and give a flavor of possible solutions, we consider perturbations of two tokamak equilibria, based on ITER. The ITER-like equilibria have their outer surface shape defined by the Solev'ev equilibrium given in Pataki et al. [2013], but scaled up so that the magnetic axis has a radial position of 6. To independently evaluate the first order QAS numerical solutions, the outer surface shape can be generated and provided to the VMEC code [Hirshman and Whitson, 1983] as input for a fully nonlinear calculation, as was done in Plunk and Helander [2018], and the result then passed to the BOOZ_XFORM code [Sanchez et al., 2000] to check the level of QAS as predicted by the theory. To produce the surface shape, the perturbation amplitude is controlled via the arbitrary small parameter in x ≈ x 0 + x 1 . Three such surfaces are shown in Fig. 2.
The solutions are valid to first order in the expansion, and it can therefore be expected that the error in QAS should scale as 2 , the confirmation of which is shown for one case in Fig. 3, where the error is measured as follows: whereB mn is the Fourier coefficient of |B| calculated in the Boozer angles by the BOOZ_XFORM. It should be noted that not all of the cases reported here match so closely with the theoretical scaling. Some have only a limited range at larger values of where the quadratic scaling is observed, and exhibiting the weaker 1 scaling for smaller values of , associated with non-QAS perturbations. Although it is expected, for instance, that numerical error in the first order solution can introduce 1 scaling which must dominate at sufficiently small , it does not appear that the linear scaling observed here is related to numerical error in the three codes being used, of the type introduced by finite resolution. It is therefore suspected that a more fundamental issue is at fault, for instance (1) the presence of force balance singularities in the fully nonlinear calculation performed by VMEC (which are formally absent from our first order calculation of x 1 ), or (2) the possibility that the problem we are solving (QAS on a single surface of non-axisymmetric perturbations) is sometimes (or always) not well posed; this issue will be investigated in future work. Nevertheless, the low observed QAS error in solutions obtained so far indicate that the numerical method developed here should be practically useful. One issue encountered with using the inverse representation of a magnetic equilibrium is that the coordinate mapping is not generally invertible for the magnetic coordinates. Invertibility breaks down, when distinct points in the magnetic coordinate space, say (ρ 1 , θ 1 , ϕ 1 ) and (ρ 2 , θ 2 , ϕ 2 ), yield the same point in physical space, e.g. x(ρ 1 , θ 1 , ϕ 1 ) = x(ρ 2 , θ 2 , ϕ 2 ). The QAS solutions here, being based on known axisymmetric equilibria will not suffer from this problem if the perturbation amplitude is chosen to be sufficiently small. However, the problem can be reliably encountered at large values of , as demonstrated in Fig. 4. What is remarkable about this phenomenon, which is associated with the perturbation overwhelming the zeroth order mapping, is that the theoretical QAS scaling tends to hold even as the singular point is approached, as demonstrated by Fig. 3. Therefore, the VMEC solutions shown here are generally chosen to correspond to a value of close to the singular point, but not so large as to create sharp features in the outer magnetic surface that require more than 10 − 20 Fourier harmonics to properly resolve in VMEC.
Using the procedure described above, the QAS solutions, though formally only perturbative, can yield strongly shaped plasma equilibria with reasonable level of QAS, and finite "external" rotational transform, as measured by the difference between the total rotational transform and that of the original axisymmetric equilibrium. This is shown by table 1, where a total of nine cases are described, corresponding to three toroidal mode numbers applied to the equilibria of three different values of constant rotational transform. Each row of the table corresponds to a single value of (although a sequence of values were generally calculated to investigate scaling).  The fraction of external rotational transform generated by the perturbation is given in the third column: Next is the root-mean-squared value of the modulus of the perturbed magnetic field, divided by zeroth order field strength, with the average performed over θ and ϕ, denoted E: where we note that the above expression assumes the first order Jacobian to be zero (QAS). This measure gives some sense of how strong the perturbation is, and may be used to estimate the size of external current distributions needed to achieve the total field. The next column provides the QAS error, Q, defined in Eqn. 44. The chosen values of are somewhat arbitrary, so it is useful to calculate normalized values to compare the various solutions. For that reason we also give inferred values of the magnetic perturbation measures E 10 and E 15 , that would be obtained for external rotational transforms of 10% and 15%, respectively. Analogous quantities for QAS error are denoted Q 10 and Q 15 . These are calculated for each case by using the fact that δ ι − ι 0 scales theoretically as 2 (confirmed for all cases), as does Q, whereas E scales as 1 . We stress that these values, being obtained from first order solutions, and not benefitting from any further optimization, should only be taken as a indicator of what can be achieved by perturbing an axisymmetric equilibrium. However, what seems clear is that, despite the fact that the perturbation of the field E scales as 1 whereas external transform scales as 2 , significant values of the latter can still be achieved at modest values of the former, as for instance shown by the ι 0 = 0.43, N = 4 case where the rms field strength fraction is not much larger that the external rotational transform fraction. An interesting qualitative feature of the QAS perturbations is the tendency to "localize" to the inboard side (e.g. at lower values of the radial coordinate R), in the sense that the amplitude of the perturbed magnetic field is larger there than on the outboard. This is quantified in the final column of the table 1, labeled "Loc", where we calculate the ratio of the root-mean-squared value of the perturbed magnetic field δB on the outboard (defined such that θ = 0) to inboard (θ = π), where the average is done only over the toroidal angle ϕ. This feature is more pronounced at larger N and lower aspect ratio, as was also observed for the vacuum case (see the appendix of Plunk and Helander [2018]; the explanation here may be similar). We note that the high-N perturbations also only weakly penetrates radially into the plasma, as the rotational transform ι can be observed to fall rapidly, from the edge, to the unperturbed value ι 0 .

Conclusion
This work gives the first set of results showing the direct construction of QAS perturbations of non-trivial axisymmetric plasma equilibria. It has been demonstrated that, despite the perturbative nature of the calculation, relatively strongly shaped stellarator equilibria may be obtained with significant external rotational transform (10 − 15%), at a similar level of average perturbed magnetic field. This finding agrees with a previous study [Ku and Boozer, 2009] using conventional numerical optimization. However, the method of the present work allows for a more extensive exploration of the design space, as the general QAS perturbation corresponds to a sum of modes with suitably non-resonant toroidal mode numbers. This potentially opens new avenues for exploring the concept of a stellarator-tokamak hybrid device.
One interesting initial finding is that relatively high-N perturbations (here as high as N = 8) seem to efficiently produce finite external rotational transform (e.g. 10%), while diminishing strongly in amplitude both radially and polloidally, and showing very good satisfaction of QAS, much less than 1% error. Such perturbations may be generated by a more "modest" distribution of coils localized to the inboard side of the plasma. Additionally, with the perturbation localized to the high-field side of the device, it should only significantly affect the radial drift of barely trapped particles, rendering the overall neoclassical transport especially small.
One benefit of perturbative studies like the present is the ability to characterize the size of the solution space of optimal stellarators. Similar to what was found by Plunk and Helander [2018], we conclude that the freedom in QAS designs comes from (1) the zeroth order equilibrium, which is in this case includes plasma profiles in addition to the two-dimensional shaping (e.g. triangularity, elongation, aspect ratio, etc.), and (2) the solution space of the QAS perturbation. For the latter, there are also some differences: first, it appears that, for fixed toroidal mode number N , the solution is unique, whereas Plunk and Helander [2018] found that solutions came in pairs -the latter situation may stem from the symmetry of the ι 0 = 0 scenario; we note that this limit cannot be approached within the current framework, as the resonant values of ι become dense near ι = 0. Second, the choice of toroidal mode number N is constrained, at least formally, by the profile ι 0 (ρ), in the sense that resonances (ι 0 = N/m) must be avoided to guarantee smooth solutions at first order, with further resonances might be considered at higher order. Therefore, the realizable solution space for QAS perturbations of a given axisymmetric equilibrium may be limited to a small number of toroidal mode numbers. Such a small space might be rapidly explored to identify QAS equilibria that satisfy additional requirements.
The success demonstrated here in directly constructing QAS solutions with an inverse method, using Boozer coordinates, gives some hope that the fully nonlinear problem may be formulated and solved in a similar fashion, i.e. with a code similar to VMEC that would obtain quasisymmetric stellarator equilibria directly, and without approximations. To accomplish this, it is necessary to identify the appropriate amount of boundary information to yield a well-posed problem; the findings of this paper should provide a useful guide in this endeavor. J 2B2 = |e 3 + ι e 2 | 2 = g 33 + 2ι g 23 + ι 2 g 22 .
Eqns. 62-63 are linearly independent (i.e. G + ι I = 0 by Eqn. 51), so we have For obtaining the Grad Shafranov equation, it is convenient to use a mixed form of B, using the toroidal part of the covariant field, and the non-toroidal part of the contravariant field: where Ψ(ψ) is the poloidal magnetic flux, and dΨ/dψ = ι . Using this form, force balance and Ampere's law imply µ 0 ∇Ψdp/dΨ = (∇ × B mix ) × B mix , which immediately yields the Grad-Shafranov equation, where Although this is a complete specification of the axisymmetric field, we require the full set of coordinates to solve our the perturbative problem, so we must develop the more general representations, i.e. B cov and B con . The functions θ(R, Z) and ν(R, Z) (and I and ι as functions of Ψ) can be obtained from additional equations derived from components of the MHD constraint, B con = B cov : G(ψ)∇ϕ + I(ψ)∇θ + K(ψ, θ, ϕ)∇ψ = ∇ψ × ∇θ − ι ∇ψ × ∇ϕ, Theφ component gives where {A, B} = ∂ R A∂ Z B−∂ R B∂ Z A, and theR andẐ components can be combined to eliminate K, yielding Then ψ can be obtained from Finally, K can be obtained from the ∇ψ component of Eqn. 68 (i.e. the condition that B cov has no component pointing out of the magnetic surface), K(R, Z) = 1 |∇ψ| 2 (G∇ψ · ∇ν + I∇ψ · ∇θ). (72)