Multi-region relaxed magnetohydrodynamics in plasmas with slowly changing boundaries --- resonant response of a plasma slab
Robert L. Dewar, Stuart R. Hudson, Amitava Bhattacharjee, Zensho Yoshida
MMulti-region relaxed magnetohydrodynamics in plasmas with slowly changingboundaries — resonant response of a plasma slab
R. L. Dewar ∗ Centre for Plasmas and Fluids, Research School of Physics & Engineering,The Australian National University, Canberra, ACT 2601, Australia
S. R. Hudson † and A. Bhattacharjee ‡ Princeton Plasma Physics Laboratory, PO Box 451, Princeton NJ 08543, USA
Z. Yoshida § Graduate School of Frontier Sciences, University of Tokyo, Kashiwa, Chiba 277-8561, Japan (Dated: October 15, 2018)The adiabatic limit of a recently proposed dynamical extension of Taylor relaxation, multi-regionrelaxed magnetohydrodynamics (MRxMHD) is summarized, with special attention to the appropriatedefinition of relative magnetic helicity. The formalism is illustrated using a simple two-region,sheared-magnetic-field model similar to the Hahm–Kulsrud–Taylor (HKT) rippled-boundary slabmodel. In MRxMHD a linear Grad–Shafranov equation applies, even at finite ripple amplitude. Theadiabatic switching on of boundary ripple excites a shielding current sheet opposing reconnectionat a resonant surface. The perturbed magnetic field as a function of ripple amplitude is calculatedby invoking conservation of magnetic helicity in the two regions separated by the current sheet. Atlow ripple amplitude “half islands” appear on each side of the current sheet, locking the rotationaltransform at the resonant value. Beyond a critical amplitude these islands disappear and therotational transform develops a discontinuity across the current sheet.
I. INTRODUCTION
The deficiencies of ideal magnetohydrodynamics(MHD) for describing typical fusion plasmas arise fromits assumption of infinite electrical conductivity, whichimplies “frozen-in” magnetic flux [1], and also its assump-tion of zero thermal conductivity, which implies frozen-inentropy (i.e. a thermodynamically adiabatic equation ofstate applying in each fluid element). The problem withfrozen-in entropy is obvious—thermal conductivity alongmagnetic field lines is in fact extremely high. The prob-lem with frozen-in flux is that it precludes changes inmagnetic-field-line topology through such reconnectionphenomena as the growth of magnetic islands at resonantmagnetic surfaces. Such islands may be excited by break-ing axisymmetry using external coils, as in stellarators ortokamaks with applied resonant magnetic perturbations(RMPs [2]), or through spontaneous tearing mode insta-bility [3].To allow for magnetic reconnection and parallel ther-mal equilibration, while retaining the non-dissipativecharacter of ideal MHD, we use a multiregion relaxation (MRx) model where complete Taylor relaxation [4] occursonly within subregions of the plasma. These relaxationregions are separated by a number (in principle, many) interfaces [5], or transport barriers, that act like thin lay-ers of ideal plasma where the ideal-MHD invariants are ∗ [email protected] † [email protected] ‡ [email protected] § [email protected] preserved, unlike the relaxation regions where the onlymagnetic invariants are helicity and total fluxes. Theseinterfaces frustrate the total relaxation postulated in theoriginal Taylor model, which cannot model the peakedcurrent and pressure profiles sought in fusion plasmaphysics. FIG. 1. Comparison of a smooth pressure ( p ) profile from aDIII-D reconstruction, using the STELLOPT code, and thestepped p -profile used in a SPEC calculation of the corre-sponding 3-D equilibrium. Also, plotted is the inverse rota-tional transform ≡ safety factor q . [Reproduced with permis-sion from Phys. Plasmas , 112502 (2012).] The MRx model is much more flexible, and has beenused as the theoretical basis for the three-dimensional(3-D) equilibrium code SPEC [6]. This has already hadsuccess [6, Sec. IV. E ] in modeling MHD equilibria us-ing experimental data from DIII-D tokamak shots, where a r X i v : . [ phy s i c s . p l a s m - ph ] M a r RMP coils were used to stochasticize the outer region ofthe plasma in order to ameliorate edge-localized modes(ELMs) [2]. Figure 1 reproduces Fig. 7 of [6], showingthe smooth pressure profile produced by a STELLOPT[7]/VMEC [8] reconstruction to fit the DIII-D data. Alsoshown is a closely approximating stepped-pressure profileused in a SPEC calculation with multiple relaxation re-gions and the same rippled boundary as used by VMEC.Evident in both are wide regions around the q = 2 and q = 3 surfaces where the pressure profile is flattened,presumably due to the presence of field-line chaos andislands generated by the RMP coils, as verified in theSPEC-produced Fig. 8 of [6]. However, while VMEC canproduce adequate macroscopic fits at a specific time inthe discharge, it is based on ideal MHD so it cannot re-solve field-line chaos and islands, and hence cannot modelthe development of equilibria exhibiting field-line chaos.The present paper presents the theoretical basis for be-lieving that SPEC can potentially do this.Taylor [9] postulates macroscopic relaxation to a force-free magnetic field B as due to turbulent fluctuationsof short ( microscale ) wavelength, scaling as the squareroot of the resistivity. Other fluctuation arguments forubiquity of relaxation may be advanced [10].However, the viewpoint we adopt in this paper is thatfield-line chaos (“stochasticity”) leads to relaxation byentangling [11] the microscopic flux tubes that each carrytheir own conserved magnetic helicity, combined with areconnection mechanism that leaves only the total mag-netic helicity conserved, as assumed by Taylor. A related,dynamical-systems-based, line of reasoning advanced byHudson et al. [6] in justification of our MRx model is thata partition of the plasma into relaxed regions invariantunder field-line flow is the only class of ideal-MHD solu-tion that avoids the mathematical pathologies in general3-D geometries identified by Grad [12]. - - - - - � / � � / � ψ / ψ � - FIG. 2. Boundaries and magnetic surfaces in a typi-cal Hahm–Kulsrud–Taylor (HKT) rippled-boundary case (seeSec. III B for details). The unperturbed boundaries are at x/a = ±
1. The shielding current sheet is along the y -axis,which separates the upper and lower relaxed-MHD regions.Note the half islands near the current sheet. As there are typically local tangential discontinuitiesacross our postulated interfaces (and a discontinuityin | B | if they support pressure differences), the inter-faces support globally extended current sheets . On the macroscale on which Taylor relaxation applies, these cur-rent sheets are of zero width—their net current may berepresented by a Dirac δ function. In developing the MRxapproach we have normally assumed that such long-livedtoroidal current sheets can exist in general 3-D equilibriaonly if the rotational transforms 1 /q on both their insideand outside faces are strongly irrational numbers. Thisassumption is based on a Hamilton–Jacobi constructionrelating the equilibrium surface currents on the two facescombined with Kolmogorov–Arnol’d–Moser (KAM) ar-guments [13–15]. The finite-wavelength stability of suchinterfaces in general 3-D geometry has yet to be inves-tigated. However, a criterion for ideal-MHD stability to localized variations may be established [16] by adaptingthe energy principle treatment of sharp-boundary equi-libria in [17] to show stability when there is no point ofzero magnetic shear (no tangential discontinuity) on asurface, suggesting jumps in rotational transform acrossinterfaces are favorable for stability.Our early development of the multi-region relaxationidea (see e.g. [6]) was based on Taylor’s minimum en-ergy variational principle [9], which produces magnetohy-dro static equilibria. The generalization to a fully fledgedfluid dynamics, Multi-region Relaxed Magnetohydrody-namics (MRxMHD) has only recently been enunciated[18]. This new formulation is a fully dynamical , time-dependent field theory whose self-consistency is ensuredby deriving it from an action principle rather than anenergy principle.The resulting model is simpler and more flexible thanideal MHD, and, we hope, is more physically applica-ble to fusion physics due to the aforementioned prob-lems with ideal MHD. The existence and stability con-siderations alluded to above are within the frameworkof the zero-Larmor-radius, dissipationless MRxMHD for-malism itself. At a minimum, internal consistency of themodel justifies using MRxMHD as a regularization anddiscretization of ideal MHD that is useful for numerical purposes. However, to argue that MRxMHD provides areasonable physical model in a specific physical contextwe need to go beyond its formal framework to considerseveral questions:Q1. Is there a microscopic mechanism for volume relax-ation to occur on a reasonably short timescale?Q2. Is there a mechanism for macroscopic current sheetsto form on a similar timescale?Q3. Are such current sheets robust enough to actas the transport-barrier interfaces postulated inMRxMHD?Q4. Is there experimental evidence that might help an-swer these questions?Current interest in modeling RMP penetration usingMHD equilibrium codes [19–23] makes it interesting toexplore whether dynamical MRxMHD can be used tomodel the development of shielding current sheets on ini-tially resonant rational surfaces as the perturbations areswitched on.In this paper, as well as developing the MRxMHDformalism for boundary-driven time-dependent pertur-bations, we illustrate it by considering the excitation ofa single resonant current sheet (to model perhaps the q = 2 or 3 surface in Fig. 1) by slowly ramping up theamplitude of a sinusoidal rippling of the boundary in avery simple geometry, namely the Hahm–Kulsrud–Taylor(HKT) rippled-boundary slab model [24] illustrated inFig. 2. We address the above questions as plausibly aswe are able in this exploratory stage of determining thedomain of applicability of the dynamical MRxMHD ap-proach.The much studied HKT model provides a macroscopicgeometry that is computationally simple and has a con-tinuous symmetry in the z direction, thus making the B -line dynamics macroscopically integrable (i.e. it hasgood flux surfaces everywhere). However, we assume, asin DIII-D [2] (addressing Q1 by appealing to Q4), thereare other RMPs simultaneously present that have littleeffect on the chosen RMP other than the crucial roleof providing the microscopic field-line chaos required tojustify invoking Taylor relaxation. We assume the onlyMRxMHD interface is a current sheet forming at the res-onant surface x = 0, thus dividing the plasma into twoequal-pressure relaxation regions, and address Q2 by cit-ing the work of Huang et al. [25] who show that currentsheets can form rapidly even in the presence of field-linechaos. Even if relaxation applies physically only near x = 0, little is lost by assuming relaxation throughouteach plasma subregion, as, in regions with constant equi-librium pressure, adiabatic relaxed MHD agrees with lin-earized ideal-marginal [26] MHD away from resonances[27, 28]. On the other hand, near the resonant currentsheet we find the linear response is weak and nonlinearterms become important, so relaxed MHD is certainlymore appropriate than linear ideal MHD [29, 30] and ar-guably more appropriate than fully nonlinear ideal MHD[31].While the concept of time is implicit in the present pa-per, we assume the switching on of boundary ripple to be adiabatic , i.e. to be sufficiently slow that the system canbe considered to evolve through a continuous sequence ofsteady states. It is our aim to add a quasi-dynamical di-mension to the static, equilibrium calculations of Loizu,Hudson et al. , [21, 22] so as to address the physical ac-cessibility of the equilibria with current sheets and dis-continuous rotational transform they calculated. Thiswork is also complementary to the ideal-MHD HKT sim-ulation study of Zhou et al. [31], who demonstrate theformation of a nonlinear ideal current sheet using a varia-tional integrator in Lagrangian labeling that enforces thefrozen-in-flux condition exactly.By definition, on resonant surfaces q is necessarily ra-tional, but the continuous symmetry of the HKT modelgets around the above-mentioned KAM existence ar- gument against rational interface q values because thecharacteristics of the Hamilton–Jacobi equation are inte-grable in this case. Furthermore we shall find that, abovea small threshold amplitude of the RMP, a discontinuityin q develops across the interface, making q on either sideno longer resonant. It has also been shown [32, 33] thatan interface located exactly at a rational surface sup-presses the associated tearing mode, essentially becausethe ideal-MHD invariants within the interface suppressreconnection.However, in a real plasma, current sheets are neitherzero width, nor perfectly conducting—at long times wemight expect [24] tearing instability to lead physicallyto island formation through internal reconnection withinthe current sheet. If this were so, the answer to Q3 wouldnot be as positive as we would like, as it would reduce thetimescale on which MRxMHD applies. But the situationmay be rescued by appeal to Q4—in the RMP experi-ments on DIII-D [2] there is a strong toroidal flow, andthe RMPs have nonzero toroidal mode number. Thuswe may invoke the flow-suppression of reconnection dis-covered by Parker and Dewar [34] in HKT geometry toargue that it is not physically unrealistic to assume theinterface current sheet to be robust on a long timescale.In Sec. II we summarize the general MRxMHD for-malism as presented in [18]. In addition we discuss thecorrect relative helicity to use in MRxMHD. We also ex-plain why we can ignore explicit consideration of flowin analyzing the RMP switch-on problem in the HKTmodel geometry, despite having invoked flow suppressionof reconnection above.The HKT model [24] is developed in Sec. III, where it isshown that MRxMHD leads to a linear Grad–Shafranovequation describing the magnetic field for ripple of ar-bitrary amplitude. In the HKT model there is reflec-tion symmetry about the resonant current sheet inter-face, which, in Cartesian coordinates x, y, z , we take tobe located on the plane x = 0. Due to the assumed re-flection symmetry, | B | is continuous across x = 0 so onlytangential discontinuities in B can arise.In Sec. IV we establish the general formalism for cal-culating solutions of the Beltrami equation in the Grad–Shafranov represention, including, in Sec. IV D, a Fourierdecomposition of Beltrami fields into plane waves. In Vwe give general expressions in Grad–Shafranov represen-tation for the magnetic energy, the vector potential andthe magnetic helicity.These formal developments are used to compute thespatially evanescent plasma response to boundary ripple,determining the Fourier coefficients from the boundaryconditions, flux and helicity constraints. It is this in-ternal disturbance that resonates at the x = 0 magneticsurface to excite the shielding-current-sheet states, whichare explored numerically in VI over ranges of initial mag-netic shear, ripple amplitude, and poloidal mode num-ber. In Sec. VI A studies are performed using amplitudessmall enough to need only the lowest spatial harmonic inthe Fourier expansion, as in [24]. Scalings with respectto initial magnetic shear and amplitude are determinedempirically. It is found that the excited current sheetconsists both of a k y (cid:54) = 0 ripple response and a net aver-age k y = 0 current. For amplitudes above the thresholdvalue at which the net average current (quadratic in am-plitude) begins to dominate the ripple currrent (linear inamplitude), a jump in rotational transform occurs acrossthe current sheet.In Sec. VI B further studies at higher amplitude areperformed by including more terms in the Fourier sums,so as to maintain a sinusoidal boundary ripple. This al-lows investigation of the dependence on poloidal modenumber, m , of the threshold amplitude for rotationaltransform discontinuity. It is found that the thresholdis highest at lowest m , presumably because the exponen-tial screening of the sheet current ripple is lowest in thiscase.Conclusions are given and directions for further workare indicated in Sec. VII, followed by Appendix A show-ing why loop integrals (cid:72) d l · A on interfaces must beincluded as MRxMHD constraint invariants, and Ap-pendix B illustrating why the simple (cid:82) A · B dV form ofthe magnetic helicity (i.e. with no vacuum helicity sub-tracted) is the correct helicity constraint invariant touse in MRxMHD under the constraint derived in Ap-pendix A.An electronic Supplement [Supp] is provided online. Itprovides further detail on deriving those equations be-low flagged by a citation to the Supplement. Also in theSupplement is Appendix C, giving a derivation of theGrad–Shafranov form of the Beltrami equation, alterna-tive to that given in Sec. III B, by deriving it directlyfrom the Woltjer–Taylor variational principle of extrem-izing magnetic energy subject to the constraint of con-stant magnetic helicity. II. THE DYNAMICAL MRXMHD MODEL
In [18] the equations for MRxMHD were derived asEuler–Lagrange equations from a Lagrangian L = (cid:88) i L i − (cid:90) Ω v B · B µ dV , (1)where the volume integration (cid:82) dV in the last term isover a vacuum region Ω v , with B denoting magnetic fieldand µ the permeability of free space. (However in theHKT model there is no vacuum region, so we do notneed this term in the present paper.) The sum (cid:80) i isover Lagrangians L i given by L i = (cid:90) Ω i L MHD dV + τ i ( S i − S i ) + µ i ( K i − K i ) . (2)Here Ω i denotes a plasma relaxation region and L MHD isthe standard MHD Lagrangian density [35, 36], ρv / − p/ ( γ − − B / (2 µ ), with ρ denoting mass density, p the plasma pressure, and γ the ratio of specific heats. As in ideal MHD, in MRxMHD mass is conserved mi-croscopically (i.e. in each fluid element dV ) by constrain-ing it holonomically [35, 36] to the strain field of La-grangian fluid element displacements.However, instead of using ideal-MHD holonomic La-grangian constraints on p and B to conserve entropy andflux microscopically we treat them as Eulerian fields. Theconstraint ∇· B = 0 is enforced by using the representa-tion B ≡ ∇× A , regarding the vector potential A as anindependently variable field, which is constrained onlyby conservation of total magnetic helicity 2 µ K i in each macroscopic subregion Ω i , and conservation of loop in-tegrals (cid:72) d l · A on the boundaries ∂ Ω i (see Appendix A).Likewise p is constrained only by conservation of totalsubregion entropies S i . These nonholonomic constraintsof constant K i and S i are then enforced through La-grange multipliers µ i ( t ) and τ i ( t ), respectively.The entropy S i , given in [18], is a functional of ρ and p but its specific form will not be needed in this paper.However the magnetic helicity constraint functional, K i ≡ (cid:90) Ω i A · B µ dV , (3)will play a critical role in our analysis of adiabatic re-sponse to ripple switch on. After the Euler–Lagrangeequation from variation of A is derived and solved ateach time, µ i is chosen so as to satisfy the helicity con-straint K i − K i = 0 (the subtracted constant K i beingthe initial value of K i ). As the problem we address inthis paper involves time-dependent geometric changes inthe boundaries we discuss below the constraints on thegauge of A required for K i to be truly invariant undersuch boundary changes, an issue not resolved in [18].We assume the plasma in each relaxation region is en-tirely enclosed by its boundary (the no-gap condition[37]), implying the tangentiality condition n i · B = 0 on ∂ Ω i , (4)where ∂ Ω i denotes the boundary of region Ω i and n i isthe unit normal at each point on ∂ Ω i .This boundary condition is intimately connected withthe question of invariance or otherwise of K i with respectto gauge changes A (cid:55)→ A + ∇ χ , which is equivalentto asking whether the volume integral (cid:82) Ω i ∇· ( B χ ) dV vanishes or not. Using Gauss’ theorem and Eq. (4)reduces this volume integral to a sum, over all topo-logically distinct cross sections S l , of surface integrals (cid:82) S l (cid:74) χ (cid:75) l n l · B dS , where (cid:74) χ (cid:75) l denotes the discontinuity(jump) in χ across S l .Thus K i is invariant under variations in χ if this gaugepotential is single valued , that is if (cid:74) χ (cid:75) l = 0 over the ν cross sections S l , where the genus ν is the number oftopologically distinct directions in Ω i [18]. (The genus inour annular tori Ω ± is 2, corresponding to the toroidaland poloidal directions.)However, single-valued gauge potentials do not exhaustthe topologically allowed possibilities: consider trans-formations of the form A (cid:55)→ A + ∇ χ l H i , where the ν non-single-valued functions χ l H i are harmonic func-tions (i.e. solutions of Laplace’s equation in Ω i ). Thejumps (periods) (cid:113) χ l H i (cid:121) l are constant over each S l , so (cid:82) S l (cid:113) χ l H i (cid:121) l n l · B dS = (cid:113) χ l H i (cid:121) l Φ l , where Φ l is the mag-netic flux through S l . Thus K i would not be gauge in-variant with respect to such transformations. However,such transformations during the evolution of the plasmaare ruled out by requiring constancy of loop integrals (cid:72) d l · A on the boundaries (the no-gaps condition, see Ap-pendix A), because allowing (cid:113) χ l H i (cid:121) l (cid:54) = 0 would changeone or more of these loop integrals on ∂ Ω i .Most of the loop integrals (cid:72) d l · A can be related to theinvariant fluxes Φ l within the plasma, but there remainsthe problem that the magnetic helicities Eq. (3), whileinvariant because of the no-gaps condition, are still notuniquely defined because of the initial gauge freedomarising from the unknown vacuum poloidal flux threadingthe toroidal vacuum-plasma interface. There are histori-cally two distinct approaches to fixing this problem, onebased on subtracting off products of line integrals (cid:72) d l · A on boundaries [9, 27, 38] and the other based on sub-tracting off corresponding vacuum helicities [39–41] toform relative helicities.Both methods involve magnetic fluxes (though repre-sented in different ways) and are both appropriate for fixed -boundary problems. However the present prob-lem involves varying boundaries and it is not clear thatthe vacuum helicity is invariant in such cases (see Ap-pendix B), casting doubt on the utility of the relative he-licity concept. This quandary is resolved in Appendix Ain favor of making invariance of boundary loop integrals (cid:72) d l · A a fundamental postulate of (no-gaps) MRxMHDbut working with the new relative helicity used in [18], K i − K i , which is relative to the initial helicity ratherthan to the vacuum helicity.Variation of A (holding µ i fixed) in Hamilton’s ActionPrinciple, δ S ≡ δ (cid:82) L dt = 0, gives a
Beltrami equation , ∇× B = µ i B (5)in each subregion Ω i , to be solved under the tangentialityboundary condition Eq. (4) on ∂ Ω i .The action principle also gives the fluid equationswithin each Ω i by varying p and the fluid positions ξ un-der the microscopic mass conservation constraint, whichin Eulerian form is ∂ρ∂t = − ∇· ( ρ v ) . (6)These variations give the compressible Euler fluid equa-tions for the mass velocity v , ρ (cid:18) ∂ v ∂t + v ·∇ v (cid:19) = − ∇ p , (7)and pressure p = τ i ρ [18].Because of the force-free nature of the magnetic fieldimplied by Eq. (5), there is no Lorentz force term in Eq. (7)— the plasma flow and magnetic field coupleonly at the interfaces . This peculiarity of dynamicalMRxMHD makes the HKT geometry particularly attrac-tive: in this geometry the deforming plasma boundariesare externally forced and the interface between the twomirror-image plasma regions is plane, thus allowing usto ignore flow in the subsequent analysis of our simpleillustrative case.In more general geometries the requirement that soundwaves not be excited also sets the slow timescale on whichan adiabatic analysis is appropriate. Suffice it to say herethat the plasma response to boundary ripple becomesincompressible in the very low frequency limit [33], so ρ and p are constant in space (and also time if the volumesof Ω i are kept constant during ripple switch on).Variation of fluid positions at the interface ∂ Ω i,j ≡ ∂ Ω i ∩ ∂ Ω j gives the force-balance condition across thecurrent sheet on this boundary (cid:115) p + B µ (cid:123) i,j = 0 , (8)the brackets (cid:74) · (cid:75) i,j denoting the jump in a quantity as theobservation point crosses the interface from the Ω i sideof to the Ω j side.However, due to the reflection symmetry about x = 0assumed in our simple HKT-like model, illustrated inFig. 2, the interface between the upper and lower re-laxation regions (which we denote by Ω + and Ω − , re-spectively) continues to be located on the x = 0 planethroughout the switching on of the RMP. Also, B remains an even function of x and hence continuous(though not necessarily differentiable) across the inter-face, and also (cid:74) p (cid:75) = 0. Thus Eq. (8) is trivially satisfiedin this paper. III. HKT-BELTRAMI SLAB MODELA. Unperturbed pseudo-toroidal equilibrium
In this subsection we limit attention to the initial, un-perturbed state of a slab plasma, before boundary rip-ple is switched on. Then all magnetic field lines canbe assumed to lie in parallel planar magnetic surfaces x = const.To relate slab geometry, as best we can, to that of atoroidal confinement device such as a tokamak, we as-sume the system to be topologically periodic in y and z , with periodic boundary condition lengths L pol = 2 πa and L tor = 2 πR , respectively. Here R is the nominal ma-jor radius of the device and a is a representative radialscale length, typically less than the mean minor radiusof an actual plasma. The y and z periodic variables arethen linearly related to the 2 π -periodic poloidal angle θ and toroidal angle ζ of a toroidal magnetic coordinatesystem, θ = ya , ζ = zR . (9)An unperturbed equilibrium field line passing throughthe point θ = θ , ζ = 0, on surface x = x is thendescribed by the line in θ , ζ space θ = θ + ι -( x ) ζ, ζ = zR , (10)where ι -( x ) [ ≡ /q ( x ), where q is the unperturbed“safety factor”] is the rotational transform on the fluxsurface. From Eq. (9), the line Eq. (10) is given in x , y space as y = aθ + ι -( x ) azR , (11)so that the general infinitesimal line element along a fieldline on an arbitrary surface x = const is d l ≡ dy e y + dz e z = (cid:104) ι -( x ) aR e y + e z (cid:105) dz , (12)by Eq. (11), with unit basis vectors e y ≡ ∇ y and e z ≡ ∇ z . Thus the equilibrium magnetic field is parallel to ι -( x ) a e y + R e z , so we may write B = B ( x ) ι -( x ) a e y + R e z (cid:112) ι -( x ) a + R = B ( x ) a e y + q ( x ) R e z (cid:112) a + q ( x ) R . (13)Thus B y /B z = ι - a/R and B z /B y = qR/a , giving thewell-known expressions ι -( x ) = RB y ( x ) aB z ( x ) , q ( x ) = aB z ( x ) RB y ( x ) . (14)Assuming an equilibrium with a sheared magneticfield, only an isolated magnetic surface(s) x = x res will resonate with a wavelike perturbation with poloidalmode number m such that k y ≡ m (2 π/L pol ) = m/a (15)and toroidal mode number n [ k z ≡ − n (2 π/L tor ) = − n/R ] [42] when the phase fronts coincide with fieldlines. That is, when k · B = 0, which, using Eq. (13),is the condition ι - ( x res ) m − n = 0, or ι - ( x res ) = nm , q ( x res ) = mn . (16)In the HKT model [24] x res = 0 and the boundaryripple is applied only in the poloidal direction, so n =0, k z = 0, and k y = m/a [43]. As | q ( x res ) | = ∞ wehenceforth use only ι -( x ) to characterize the pitch of theequilibrium field.We depart from [24] in taking both unperturbed andperturbed magnetic fields to obey Eq. (5), though theLagrange multiplier µ must change slightly with increas-ing ripple amplitude in order to satisfy the constant-magnetic-helicity constraint. (However µ will be the same in both Ω + and Ω − due to the assumed symme-try about x = 0.) Denoting the unperturbed value of µ by µ (not to be confused with the vacuum permeability µ ) we find the unperturbed solution of Eq. (5), B (0) ( x ) = B (sin µ x e y + cos µ x e z ) , (17)where B is a constant. Using Eq. (14) we find the rota-tional transform in the unperturbed state, ι -( x ) = Ra tan µ x . (18) B. Grad–Shafranov (GS) representation
We follow Hahm and Kulsrud [24] in using a flux-function representation for the full, perturbed magneticfield, defining ψ ( x, y ) such that B = F ( ψ ) e z + e z ×∇ ψ (19)where F ( ψ ) is B z expressed as a function of ψ .Note that Eq. (19) implies that B ·∇ ψ ≡
0, i.e. B is everywhere tangential to level surfaces ψ = const, sothat ψ has the property of being a label for magneticsurfaces . Figure 2 shows rippled magnetic surfaces givenby constructing representative contours of ψ ( x, y ) after m = 2 [see Eq. (15)] sinusoidal boundary ripple of am-plitude α = 0 .
21 [see Eq. (49)] has been switched on,starting from the “tokamak-relevant” case ( µ a = 0 . ψ is not unique, because it can be changedby a constant amount without changing the observable e z ×∇ ψ , the poloidal magnetic field . Such a baseline shiftalso changes the functional form of F to retain the invari-ance of the observable B z , the toroidal magnetic field .To remove this arbitrariness we set the baseline for ψ by fixing it on both boundaries x = ± x bdy ( y ) to be theconstant value ψ = ψ a : in the HKT model we assume ψ ( x, y ) to be even in x and to increase away from x = 0,as | y | increases, up to ψ a . Though ψ ( x, y ) is continuousacross the current sheet at x = 0, its derivative ∂ y ψ isin general discontinuous there, so it is sometimes con-venient to consider ψ ( x, y ) as defined on two Riemannsheets intersecting along the cut at x = 0.Note that Eq. (19) implies ∇× B = ∇ ψ e z − F (cid:48) ( ψ ) e z ×∇ ψ (20)so that, crossing Eq. (20) with Eq. (19) ,( ∇× B ) × B = − ( ∇ ψ + F F (cid:48) ) ∇ ψ . (21)For force-free fields, such as those described by Eq. (5),the left-hand side of Eq. (21) vanishes, leaving us withthe equation ∇ ψ + F F (cid:48) = 0 , (22)which is the Grad–Shafranov (GS) equation in slab ge-ometry in the special case p (cid:48) = 0. As will be shownbelow, this is a linear equation and reduces, in the limit µa →
0, to the Laplace equation assumed in [24]. TheGS representation of Beltrami solutions is also useful inaxisymmetric toroidal geometry [44].Substituting Eq. (19) and Eq. (20) in Eq. (5) we gettwo equations defining a Beltrami field in the GS repre-sentation, ∇ ψ = µF (23)and F (cid:48) ( ψ ) = − µ , (24)which are consistent with the GS equation, Eq. (22).Equations (22–24) are also derived variationally fromfirst principles in [Supp]:Appendix C by minimizing mag-netic energy at constant helicity.Integrating Eq. (24) gives F ( ψ ) = C − µψ (25)where C is a spatial constant, though it is not invariantunder application of ripple. Also, as the left-hand sideof Eq. (25), F = B z , is a physical observable, C mustcounterbalance the arbitrary baseline constant includedin ψ .In the following we find it useful to define the area-weighted average of an arbitrary function f over a surfaceof section across the upper relaxation region Ω + as f ≡ (cid:104) f (cid:105) ≡ A + (cid:90) πa dy (cid:90) x bdy ( y )0 dx f ( x, y ) , (26)where A + is the cross-sectional area over one topologicalperiodicity length, A + ≡ (cid:90) πa dy (cid:90) x bdy ( y )0 dx = 2 πa (cid:104) x bdy (cid:105) , (27)the total cross-sectional area across the whole plasma be-ing A ≡ A + + A − = 2 A + . (In general two-relaxation-region problems we would need to define separate aver-aging operators (cid:104) f (cid:105) ± in Ω + and Ω − , but the reflectionsymmetry assumed in the HKT model means these aver-ages are equal for even parity functions and the negativeof each other for odd parity functions.) Note we have in-troduced two equivalent averaging notations, (cid:104)· · ·(cid:105) beinga useful alternative to · · · for lengthy expressions.To decompose C into an invariant part and a geomet-rically dependent part we average Eq. (25) over a surfaceof section of Ω + , as in Eq. (26), to give C = F + µψ , (28)using which Eq. (25) becomes F = F − µ (cid:101) ψ , (29) where (cid:101) ψ ≡ ψ − ψ (30)is the deviation from the mean poloidal flux ψ . Deter-mination of the non-invariant quantity ψ (and hence C )will be discussed in Sec. IV B. ( a )- - � / � � � = � � ( b )- - � / � - - ι - = � - � ( c )- - � / � - - Φ � ( d )- - � / � Ψ � FIG. 3. Unperturbed (plane boundary) profiles in the RFP-relevant case µ a = 1 .
4. (a) Toroidal field B z = F = B cos µ x . (b) Rotational transform ι -. c) Toroidal fluxΦ. (d) Poloidal flux Ψ. Parameters and units are such that a = B a = 1, R = a cot µ a . ( a )- - � / � � � = � ( b )- - � / � - - ι - = � - � ( c )- - � / � - - - Φ ( d )- - � / � Ψ FIG. 4. Unperturbed (plane boundary) profiles in thetokamak-relevant case µ a = 0 .
2. (a) Toroidal field B z = F = B cos µ x . (b) Rotational transform ι -. c) Toroidal fluxΦ (0) ( ψ (0) ( x )). (d) Poloidal flux Ψ. Parameters and units areas in Fig. 3, the choices of µ in the two figures being discussedin Sec. III D. C. Fluxes
The poloidal flux Ψ( ψ ) between the current sheet,where ψ = ψ cut , and a magnetic surface ψ = ψ s is thesurface integral of the flux density e y · e z ×∇ ψ = ∂ x ψ over an area in any plane y = const bounded by the cur-rent sheet x = 0, the magnetic surface labeled by ψ s , andthe lines z = const and z = const + 2 πR ,Ψ ± ( ψ s ) ≡ (cid:90) πR dz (cid:90) x ± ( ψ s | y )0 ∂ψ∂x dx = 2 π ( ψ − ψ cut ) R , (31)where x = x ± ( ψ s | y ) denotes the upper (+) or lower ( − )branch of the solution to the equation ψ ( x, y ) = ψ s , forgiven y . (Where y is arbitrary for magnetic surfaces out-side half islands such as are seen in Fig. 2, but, for defin-ing the “private flux” within such an island, y must ob-viously be restricted to lie within the island.) The linearrelation between the poloidal flux Ψ and the function ψ justifies the terminology poloidal flux function for thelatter.The fact that Ψ( ψ ( x, y )) is an even function of x inthe HKT model is illustrated in Figs. 3 and 4 for theunperturbed case Eq. (17) (in which special case there isno current sheet, so Ψ is differentiable at x = 0).Assuming here the magnetic surface spans the fullpoloidal periodicity length 2 πa (i.e. it is not in a halfisland) we also define the toroidal flux Φ( ψ ) as a mag-netic surface quantity by integrating the toroidal mag-netic field B z = F ( ψ ) over one period in y between theresonant surface x = 0 and the given magnetic surface x = x ± ( ψ s | y ) in Ω ± ,Φ ± ( ψ s ) ≡ (cid:90) πa dy (cid:90) x ± ( ψ s | y )0 F ( ψ s ) dx . (32)[Note that Φ + ( ψ s ) = − Φ − ( ψ s ) so Φ( ψ ( x, y )) is an oddfunction of x , as illustrated in Figs. 3 and 4.] We gen-eralize the “safety factor” q , defined for the unperturbedfield in Eq. (14), as q ( ψ ) ≡ d Φ /d Ψ, which, in the GSrepresentation Eq. (19), can be written q ± ( ψ s ) = ± F ( ψ s )2 πR (cid:73) pol dl | ∇ ψ | (33)where dl ≡ ( dx + dy ) / is an element of length alonga contour ψ = ψ s running between y = 0 and y = 2 πa .This general definition applies equally to the perturbedand unperturbed system.However, as mentioned in Sec III, it is more convenientto work with the reciprocal of q , the rotational transform , ι - ± ( ψ ) = d Ψ ± d Φ ± . (34)As illustrated for the unperturbed case in Figs. 3 and 4, ι - is an odd function of x . In this special case it is con-tinuous at x = 0, but for large enough ripple amplitudewe shall find that it may be discontinuous there. Dotting both sides of Eq. (19) with e z and integratingover one wavelength of the cross section, we thus find ourfirst invariant Φ tor ≡ Φ + ( ψ a ) − Φ − ( ψ a ) = 2Φ + ( ψ a ), the total toroidal flux , to beΦ tor = A F . (35)The toroidal flux and magnetic helicity contained be-tween the everywhere perfectly conducting boundaries x = ± x bdy ( y ) are conserved throughout, from switch onto reconnection. Also, to avoid the external work re-quired to change the mean toroidal field F and pressure p we assume the rippling of the walls is done in such wayas to preserve area, A = A = 4 πa , (36)which from Eq. (27) is ensured by requiring (cid:104) x bdy ( y ) (cid:105) = a , (37)making the adiabatic plasma response incompressible[33].With area thus conserved, toroidal flux conservation isequivalent to invariance of F , which also applies sepa-rately in both upper and lower relaxation regions due tothe assumed reflection symmetry. Thus in both regionsthe toroidal flux conservation condition is equivalent tothe constraint F − F = 0 (38)during switch-on of the boundary ripple perturbation,where F denotes the unperturbed value of the meantoroidal field, which is calculated below. D. Unperturbed state in GS representation
From Eq. (17) the unperturbed toroidal magnetic fieldis F ( x ) ≡ B cos µ x . Thus F = B (cid:104) cos µ x (cid:105) , (39)where (cid:104) cos µ x (cid:105) = 1 a (cid:90) a dx cos µ x = sin µ aµ a ∼ − µ a O (cid:0) ( µ a ) (cid:1) . (40)Note the interesting fact that F has zeros at µ a = πn (41)for integer n (cid:54) = 0, corresponding to extreme reversed-fieldstates. However we shall not consider such large valuesof µ a in this paper.For use later in this paper it will be found useful todefine U , the plane-slab unit-vector solution of Eq. (5)with general µ , U ( x | µ ) ≡ sin µx e y + cos µx e z . (42)In terms of U , the unperturbed field B (0) ( x ), given aboveby Eq. (17), can be represented as B U ( x | µ ).In the GS representation, U can be represented by ψ U ( x | µ ) ≡ µ (1 − cos µx ) = 2 µ sin µx , (43a) F U ( x | µ ) ≡ cos µx = 1 − µψ U . (43b)Setting ψ (0) ( x ) = B ψ U ( x | µ ), F ( ψ (0) ) = B F U ( x | µ )in Eq. (19) verifies that B (0) ( x ) = B U ( x | µ ). Note thatwe have chosen the arbitrary constant in the sheared-fieldflux function to be such that ψ U = 0 on the y -axis. Notealso the useful identity ψ (cid:48) + µ ψ = 2 µψ U , (44)where ψ (cid:48) U ≡ ∂ψ U /∂x .For devices such as the reversed-field pinch (RFP), µ a is O (1) [9]. Such a case, close to the value π/ B is dominated bythe toroidal magnetic field B z = F ( ψ ), which can bemodeled in MRxMHD by taking µ a (cid:28)
1. Althoughlarge, B z is then approximately constant, the interest-ing physics being in the behavior of the poloidal field B y = ∂ x ψ . Thus, rather than specifying B directly, wefind it more convenient to specify the boundary poloidalfield , B a = ψ (cid:48) a ≡ B ψ (cid:48) U ( a | µ ). Using Eq. (43a) B is thengiven by B = B a sin µ a , (45)which diverges in the limit µ a → ψ a ≡ B ψ U ( a | µ ). Using Eqs. (43a) and 45, thisis given by [Supp] ψ a = B a µ tan µ a , (46)which approaches aB a as µ a → ι - (0) = ± x = ± a . Then Eq. (18) gives a/R =tan µ a ≈ µ a . For our standard tokamak-relevant refer-ence case we take µ a = 1 / a/R = tan µ a ≈ . R/a ≈ IV. RIPPLED STATES IN GSREPRESENTATIONA. Rippled boundary conditions
The perfectly-conducting boundary walls are then de-formed (rippled) by switching on, over a time short com-pared with the reconnection timescale for the result-ing long-lived current sheet, wavelike perturbations with(fundamental) poloidal wave number k y . Reflection sym-metry of both walls and plasma about the y -axis is as-sumed, so ψ remains an even function of x .From Eq. (42), U (0) = e z so the resonance condition k · B = 0 (see Sec. IV D) is satisfied at x = 0, i.e. along the y -axis, where a shielding current sheet of full width 2 πa initially forms to prevent island formation. This currentsheet cuts Ω into two disjoint subdomains, Ω + , between x = 0 and upper boundary x = x bdy ( y ), and Ω − , between x = 0 and the lower boundary x = − x bdy ( y ).We shall find the fully shielded state , immediately af-ter the boundary perturbation is switched on, by assum-ing Taylor relaxation occurs independently in Ω ± [a spe-cial case of the Multi-Region Relaxed MHD (MRxMHD)problem [6]], so the perturbed initial magnetic fields inΩ ± , before reconnection of the shielding current sets in,are found by solving Eq. (5) under the boundary andother conditions discussed below.On an equilibrium current sheet it can be shown [14,Appendix A] that the normal component of B must van-ish. In terms of the representation Eq. (19), ψ is thusconstant on both sides of the current sheet. Also ψ mustbe continuous across the current sheet.By definition, in the fully shielded state no poloidalflux has yet been reconnected through x = 0. Also, nopoloidal flux can escape through the perfectly conductingwalls. Thus the current sheet boundary condition ψ = ψ cut at x = 0 applies, with ψ cut fixed at its unperturbedvalue, which we have chosen to be ψ cut ≡ . (47)On the rippled walls the flux function remains ψ = ψ a toconserve poloidal flux. Likewise toroidal flux is trappedbetween the walls and current sheet, consistently withEq. (38).We consider two methods for defining the boundarywaveform function x bdy ( y ):Bdy-1. The indirect implicit boundary method [45] wherewe specify the boundary conditions on ψψ ( ± a, y ) − (cid:104) ψ ( ± a, y ) (cid:105) = 2 α ψ a cos k y y . (48)The factor ψ a , defined in Eq. (46), is introducedin Eq. (48) to make the ripple amplitude param-eter α dimensionless, the factor 2 being to make α the same as in Ref. 45 in the limit µ a → x bdy ( y ) is then defined by thecontour ψ = ψ a : x = x bdy ( y | α ), with0 ψ constructed so as to enforce Eq. (47), thearea/toroidal flux conservation equation Eq. (36),and the magnetic helicity constraint K i − K i = 0.The domain Ω = Ω − ∪ Ω + is now completely spec-ified, its complete boundary ∂ Ω being the unionof the two external boundaries x = ± x bdy ( y ) andthe internal boundary formed by the cut alongthe y -axis. Note that we do not linearize withrespect to α so x bdy ( y ) is not an exact sinusoidin this method.Bdy-2. The direct explicit boundary method , as used inFig. 2, where we prescribe x bdy to be exactly si-nusoidal, x bdy ( y ) = a (1 − α cos k y y ) . (49)The advantage of method Bdy-1 is that it allows a simpleclosed-form solution of the perturbed Beltrami equationsimilar to the type found by Hahm and Kulsrud [24].The disadvantage is that x bdy ( y ) becomes highly non-sinusoidal at quite moderate ripple amplitudes α and themethod breaks down as α increases further.Method Bdy-2, on the other hand, can treat large rip-ple amplitudes, as evidenced for instance in Fig. 2. Itsdisadvantage is that ψ must be expanded in an infiniteseries of higher harmonics if Eq. (49) is to be satisfiedexactly [see Eq. (64)]. However, in practice a good ap-proximation can be found with a reasonable number ofexpansion functions. For small α the two methods areequivalent. B. GS equation boundary conditions
In the present application of the GS formulation oursimple boundary condition on the current sheet, Eq. (47)allows us to identify C immediately as F (0), the toroidalmagnetic field on the current sheet at x = 0. However,for α (cid:54) = 0, F (0) is not known a priori but must be deter-mined along with ψ in the solution procedure.Unlike F (0), F is a known constant, from Eq. (38).Also (cid:101) ψ is independent of the arbitrary constant in ψ andis thus a better flux variable to work with. Like ψ itmust be constant on the boundaries and current sheets,obeying the boundary conditions (cid:101) ψ ( x bdy ( y ) , y ) = ψ a − ψ , ∀ y (50a) (cid:101) ψ (0 , y ) = ψ cut − ψ , (50b)Unlike ψ , neither of these boundary values is known a priori . Instead ψ ( α ) needs to be determined, alongwith µ , under the Taylor relaxation constraints and theconstraint implied by Eq. (30), (cid:68) (cid:101) ψ (cid:69) ≡ . (51) Substituting Eq. (29) in Eq. (23) gives a linear GSequation in the form of an inhomogeneous Helmholtzequation, ( ∇ + µ ) (cid:101) ψ = µF . (52)Averaging both sides of Eq. (52) and using Eq. (51) wefind (cid:68) ∇ (cid:101) ψ (cid:69) = µF with (cid:68) ∇ (cid:101) ψ (cid:69) = − A + (cid:90) ∂ Ω + n ·∇ (cid:101) ψ dl , (53)where the right-hand side is found by applying Gauss’theorem, with dl = ( dx + dy ) / an element of lengthalong two contours, the upper wall x = x bdy ( y ) and theupper side of the current sheet x = 0 between y = − πa and y = πa , n being the inward directed unit normal ateach point on ∂ Ω + .Noting from Eq. (43a) that ( ∇ + µ ) ψ U = µ , solvingEq. (52) can be reduced to the solution of a homogeneousequation using the ansatz (cid:101) ψ ( x, y ) = F ψ U ( x | µ ) + (cid:98) ψ ( x, y ) . (54)where (cid:98) ψ obeys the homogeneous Helmholtz equation ( ∇ + µ ) (cid:98) ψ = 0 (55)under the boundary and averaging conditions followingfrom Eqs. (50a - 51) (cid:98) ψ ( x bdy ( y ) , y ) = ψ a − ψ − F ψ U ( x bdy ( y ) | µ ) ∀ y , (56a) (cid:98) ψ (0 , y ) = ψ cut − ψ ∀ y ∈ cuts , (56b) (cid:68) (cid:98) ψ (cid:69) = − F (cid:104) ψ U (cid:105) . (56c)The parameter (cid:104) ψ U (cid:105) is a functional of the boundaryshape, which may or may not be known a priori depend-ing on whether we use method Bdy-1 or Bdy-2. Theunknown parameters to be solved for are µ , ψ , and coef-ficients of terms in the ansatz for (cid:98) ψ ( x, y ) to be discussedin Sec. IV D. C. Unperturbed state: include only k y = 0 In this section we calculate expressions for initial (un-perturbed) states before ripple is switched on. Theunperturbed state is defined by α = 0, with only x -dependent magnetic field B (0) ( x ) = B U ( x | µ ), where U is given by Eq. (42). The corresponding flux functionis ψ ( x ) ≡ ψ (0) ( x ) ≡ B ψ U ( x, µ ) , (57)where ψ U is defined in Eq. (43a).Applying the averaging operator defined in Eq. (26),with unperturbed boundary, we find [Supp] ψ = B µ (1 − (cid:104) cos µ x (cid:105) ) (58)1where (cid:104) cos µ x (cid:105) is defined in Eq. (40). Hence, inEq. (30), (cid:101) ψ ( x ) ≡ B ψ U ( x | µ ) − ψ = B µ (cid:18) sin µ aµ a − cos µ x (cid:19) = aB (cid:20)(cid:18) x a − (cid:19) µ a + O (cid:0) µ (cid:1)(cid:21) . (59)The decomposition Eq. (54), (cid:101) ψ ( x ) = F ψ U ( x, µ ) + (cid:98) ψ ( x ), implies (cid:98) ψ ( x ) ≡ (cid:101) ψ ( x ) − F ψ U ( x, µ )= F − B µ cos µ x = B µ (cid:18) sin µ aµ a − (cid:19) cos µ x , (60)which is in the kernel of ∇ + µ as required by Eq. (55). D. Rippled state: include k y (cid:54) = 0 terms Here we generalize the Hahm–Kulsrud [24] solutions byexpanding in a basis of plane-wave Beltrami solutions—a k y = 0 solution B U ( x | µ ), with B and µ to be determined,and “ripple” solutions that are periodic in the y directionand exponential in the x direction.The general solution of the Beltrami equation Eq. (5)is a superposition of divergence-free plane wave solutionswith wave vector k (cid:48) = ± k (cid:48) x e x ± k (cid:48) y e y such that k (cid:48) = µ .To satisfy 2 πa topological periodicity in the y directionwe introduce the poloidal mode number , m (cid:48) = 0 , , , . . . ,such that k (cid:48) y ≡ m (cid:48) /a . Thus Eq. (5) implies k (cid:48) x ≡ µ − k (cid:48) y = µ − m (cid:48) a . (61)The y -independent solutions considered in the previoussection correspond to m (cid:48) = 0, giving k (cid:48) x = ± µ . Ripplesolutions of Hahm–Kulsrud type require imaginary k (cid:48) x , sowe consider only the case | µ | < k (cid:48) y and set k (cid:48) x = ± iκ m ( µ ),where κ m (cid:48) ( µ ) = ( k (cid:48) y − µ ) / ≡ (cid:18) m (cid:48) a − µ (cid:19) / . (62)In the above, m (cid:48) is the fundamental poloidal modenumber m of the imposed ripple or a harmonic, m (cid:48) = lm ,where l = 0 , , , , . . . . We denote the fundamental ripplewavelength by λ m = 2 πam . (63)The requirement ∇· B = 0 is ensured by using the F, ψ representation Eq. (19), the most general 2 πa -periodic solution of Eq. (55), analytic on the halfplane x > (cid:98) ψ + ( x, y ) = c cos µx + ∞ (cid:88) l =1 c lm cos lmya cosh κ lm x + d sin µx + ∞ (cid:88) l =1 d lm cos lmya sinh κ lm x . (64)with the corresponding solution (cid:98) ψ − ( x, y ) on the halfplane x < (cid:98) ψ − ( x, y ) = (cid:98) ψ + ( | x | , y ).Our generalized Hahm–Kulsrud-type solutions are su-perpositions of the form Eq. (54), ψ = ψ + F ψ U + (cid:98) ψ on thecut x, y -plane, with the branch (cid:98) ψ = (cid:98) ψ ± being chosen ac-cording as x ≷ { c } , { d } , ψ , ψ cut and µ being chosen as described qualitatively in Secs. IV A.Method Bdy-1 uses only l = 0 and l = 1 while MethodBdy-2 in principle uses l = 0 , . . . , ∞ . V. ENERGY AND HELICITY IN GSREPRESENTATIONA. Relative magnetic energy density
Rather than use the total magnetic energy W = (cid:82) Ω B / µ dV we find it neater to define W , the averageenergy per unit volume, multiplied by µ . That is, W is defined as (cid:10) B (cid:11) /
2, which we henceforth refer to sim-ply as the energy density, there being no thermal energy,as we have taken p = 0, and kinetic energy being neg-ligible in adiabatic processes. Although we defined theaveraging operation in Eq. (26) to be over Ω + , the sameresults apply in Ω − due the assumed symmetry and thefact that the energy (and helicity—see below) densitiesare even functions of x .In GS representation Eq. (19), and using Eq. (29), theenergy density is given by W = 12 (cid:10) F + | ∇ ψ | (cid:11) = F W Σ , (65)where we have defined W Σ , the energy density in theinternally generated (i.e. non-vacuum) field, as, usingintegration by parts, [Supp] W Σ = 12 (cid:68) | ∇ (cid:101) ψ | + µ (cid:101) ψ (cid:69) = 12 ( ψ a − ψ ) µF + ψ a − ψ cut aλ m J + + µ (cid:68) (cid:101) ψ (cid:69) , (66)the ripple wavelength λ m being defined by Eq. (63). Inthe above manipulations we have used Eqs 30–53 and2the boundary conditions Eqs. 50a and 50b and have de-noted the line integrals on the top/bottom surface of thecurrent sheet cut over one ripple wavelength by J ± , J ± ≡ ± (cid:90) λ m / − λ m / (cid:101) ψ x (0 ± , y ) dy (67)where (cid:101) ψ x ( x, y ) ≡ ∂ x (cid:101) ψ ( x, y ). Taking into account thesymmetry about the y -axis we have J − = J + .To discover the physical meaning of J ± note that thestrength of a sheet current is j ∗ = (cid:114) n + ·∇ (cid:101) ψ (cid:122) / µ , where n + is the unit normal at each point on the upper surfaceof the current sheet and (cid:74) · (cid:75) denotes the jump in thisnormal direction. In our case the current sheet is at x = 0and n + = e x , so j ∗ ( y ) = 1 µ [ (cid:101) ψ x (0+ , y ) − (cid:101) ψ x (0 − , y )] . (68)Integrating Eq. (68) along the current sheet we see thatthe total current per ripple period in the current sheet is J/ µ where J ≡ J + + J − = 2 J + .Using Eq. (59) in Eq. (66) we find the unperturbedinternal magnetic energy, W Σ0 = B (cid:104) − (cid:104) cos µ x (cid:105) (cid:105) ∼ B a O ( a µ )] , (69)The latter form being found by using Eq. (40), eliminat-ing B using Eq. (45), and expanding in µ a .From Eq. (65) we see that, in the “tokamak-relevant”small- µ a limit (see Sec. III D), W / W Σ0 is dominated bythe large vacuum toroidal field energy term F / W Σ0 ,which, from Eq. (45), is seen to diverge like 1 /µ a in thislimit. However F is invariant under application of rip-ple because of our constant-volume constraint Eq. (36).Thus it is more instructive to work with the relative en-ergy density , ∆ W = W Σ − W Σ0 ≡ ∆ W Σ , (70)where the vacuum toroidal field energy has cancelled out. B. Relative magnetic helicity
It is readily verified, by calculating B = ∇× A andcomparing with Eq. (19), that A = − ψ e z + 1 µ e z ×∇ (cid:101) ψ (71)is a vector potential satisfying the requirement (see Ap-pendix A) of invariance of the loop integrals (cid:72) pol d l · A =2 πaψ cut and (cid:72) tor d l · A = 2 πRψ cut on the current sheet(topologically a torus). By eliminating the linear term in ψ from Eq. (71) usingEq. (25), an alternative form, A = B − C e z µ (72)is found that will be useful below for relating energy andhelicity.By analogy with W in Sec. V A we define the averagehelicity density K ≡ (cid:104) A · B (cid:105) / relative helicitydensity ∆ K ≡ (cid:104) A · B (cid:105) / − (cid:104) A · B (cid:105) /
2, so that, compar-ing with Eq. (3) the helicity constraint K i − K i = 0 isequivalent to ∆ K = 0 . (73)From Eq. (72), then Eq. (65) and e z · B = F , and elim-inating C with Eq. (28), we find the general expression[Supp] K ≡ (cid:104) A · B (cid:105) W Σ µ − ψ F , (74)expressing a linear relation between helicity and energy.The importance of the constant offset in this relation.arising in general from a surface term, (cid:82) ∂ Ω A × B · n dS , indiscussing minimum energy states was originally pointedout by Reiman [46, 47]. Below we use this result to obtainthe analytical version for K that is used in the numericalstudies presented in Sec. VI. As K = K after the un-knowns are solved for numerically, Eq. (74) is used againto calculate the numerical value of W Σ .The unperturbed helicity density K is simply a specialcase of Eq. (74). Subtracting it from K and using theinvariance of F and Eq. (66) gives∆ K = − ( ψ − ψ ) F + ψ a − ψ cut aµλ m J + + µ (cid:68) (cid:101) ψ (cid:69) − µ (cid:68) (cid:101) ψ (cid:69) . (75) VI. SHIELDED RMP SOLUTIONS
To explore the properties of this model quantitativelywe find Beltrami solutions satisfying rippled bound-ary conditions, appropriate to the methods discussed inSec. IV A, and the boundary conditions Eqs. (56b–56c)on the current sheet. To provide a complete set of equa-tions to solve numerically we also impose helicity conser-vation, Eq. (73). As Eq. (37) conserves cross-sectionalarea, toroidal flux conservation is equivalent to the con-servation of F , which is thus still given by Eq. (39). A. HKT-like rippled boundary condition
In this subsection we use the method Bdy-1 to specifythe rippled boundaries. Units such that 2 π/k y = a = 13 FIG. 5. Bdy-1, λ = a = m = 1: Plots of [ µ ( α, µ ) − µ ] / ( µ α ) vs . µ for selected values of α in the range [0 , . α -dependence hasbeen scaled out to high accuracy. have been used, as well as the choice λ ≡ π/k y = a .In this subsection we also assume L pol = a , so m = 1,but in the standard convention used elsewhere in thispaper, L pol = 2 πa , this would be equivalent to taking m = 2 π ≈ x = 0(denoted in [24] by subscript I, here denoted by subscript“sh”), is the special case of Eq. (64) (cid:98) ψ sh ( x, y ) ≡ αψ a sinh( κ m a ) (cid:18) | sinh κ m x | cos k y y + γ S κ m µ | sin µx | (cid:19) − ψ cos µx , (76)where, from Eq. (62), κ m ( µ ) = ( m /a − µ ) / , where µ = µ sh ( α, µ ), is to be determined. Comparing withEq. (64), we have set c and all l > l = 0 and l = 1 terms. Thecoefficient d = 2 αψ a / sinh( κ m a ) has been chosen so that ψ ( a, y ) automatically satisfies Eq. (48), the parameter α setting the ripple amplitude. Also the coefficient c = − ψ has been chosen so that ψ (0 , y ) automatically satisfiesEq. (56b), with the poloidal flux conservation conditionEq. (47). FIG. 6. Bdy-1, λ = a = m = 1: Plots of [ ψ ( α, µ ) − ψ (0 , µ )] /α (in units such that B a = 1) vs . µ for selectedvalues of α in the range [0 , . α in the legends is the same as that for thecorresponding curves, shown online by color.) By analogy with Eq. (13) of Ref. 45, in our expressionfor (cid:98) ψ sh above we have renormalized the amplitude d ofthe sin µx term by setting d = 2 αψ a γ S κ m /µ sinh( κ m a ), FIG. 7. Bdy-1, λ = a = m = 1: Plots of ∆ W Σ ( α, µ ) /α vs . µ for selected values of α in the range [0 , . λ = a = m = 1: Plots of γ S ( α, µ ) /α vs . µ for selected values of α in the range [0 , . the dimensionless parameter γ S adding a constant termto the sheet current on the x -axis—from Eq. (76), thejump in ∂ x (cid:101) ψ in the expression, Eq. (68), for the sheetcurrent j ∗ is given by (cid:114) ∂ x (cid:98) ψ sh (cid:122) ≡ αψ a κ m sinh κ m a (cos k y y + γ S ) , (77)so that the total-current parameter J , defined belowEq. (68), becomes J = 4 αψ a κ m λ m sinh κ m a γ S . (78)The boundary function x shbdy ( y | α ) is determined fromEq. (56a), the three parameters µ = µ sh ( α, µ ), ψ = ψ sh ( α, µ ), and γ S = γ shS ( α, µ ) being determined bysolving the 3 simultaneous equations Eq. (37), Eq. (51)and Eq. (73), with K given by Eq. (75). The average en-ergy density can then be found from Eq. (66) or Eq. (74).Numerical results showing the µ a -dependence of µ , ψ , W Σ , and γ S in the case k y = 2 π/a , also used inRef. 45, in units such that a = 1 , B a = 1 [see Eq. (45)],are given in Figs. 5–9. The fact that the scaled curvesfor different values of α are almost identical show thatthe small-amplitude scalings µ − ∝ α , ψ − ψ ∝ α , W Σ − W Σ0 ∝ α , and γ S ∝ α are a good approxima-tion for the range α < .
05 depicted (becoming exactin the limit α → µ a < µ , butvary more rapidly above this range as µ a approachesthe value π/ ≈ .
57 at which B (0) z = B cos µ x reversessign at x = ± a [cf. Fig. 3(a).].Quite apart from demonstrating a mathematical scal-ing law, the physics shown in Fig. 7 is worthy of remark4because of the sign of ∆ W Σ and its reversal at large val-ues of µ a . First note from Eq. (70) that ∆ W Σ and ∆ W are equal, so the negative values of ∆ W Σ at small-to-moderate µ a means the total magnetic field energy inthe plasma decreases as ripple is imposed. That is, theplasma does work on the boundary. It is tempting tointerpret this as implying that a slab plasma confined byMRxMHD interface current sheets at the vacuum-plasmaboundaries would be unstable toward spontaneous rip-pling.However, this conclusion would be unwarranted as thetotal energy includes not only the O ( α ) wave energy inthe m (cid:54) = 0 ripple, but also small, O ( α ), nonlinear correc-tions to the energy in the m = 0 background state. Thesign of such a background energy correction is dependenton the precise nature of the rippling process as it couldeasily be changed by allowing an O ( α ) change in thecross-sectional area A , rather than arbitrarily imposingits constancy through through the constraint Eq. (36).Furthermore, proper stability analysis of a free-boundaryMRxMHD plasma [33, 48] must include the change in thevacuum energy outside the plasma. Such issues will bediscussed further elsewhere. FIG. 9. Bdy-1, λ = a = m = 1: Plots of γ S ( α, µ ) /α vs . α for selected values of µ a within the restricted range [0 , γ S /α is approximately constant with respect to bothvariables in these ranges.FIG. 10. Bdy-1, λ = a = m = 1: Plots of the jump in thegradient of ψ , Eq. (77), vs . y/λ for µ = 1 . /a and selectedsmall values of α , showing the occurrence of current-densityreversal for the two smallest values. The linear α -dependence of the dimensionless param-eter γ S , shown more explicitly in Fig. 9, is particularlyinteresting, in the light of Eq. (77), as it means the m = 0response J scales as α , i.e. it is nonlinear . Thus, for small-enough values of α , the m = 0 response (the termin γ S ) is dominated by the linear response (the term incos k y y ). In this case the singular current density re-verses sign over a range of y . However, as shown inFig. 10, above the very small threshold value, α thr , atwhich γ S = 1, the m = 0 nonlinear response becomesincreasingly dominant and the sheet current becomes ofconstant sign. (Both figures use the same parameters, a = λ = 1, as the previous plots.) FIG. 11. Bdy-1, λ = a = m = 1: Level surfaces of ψ (mag-netic surfaces) in the case µ = 1 . /a , α = 0 . < α thr ,showing pairs of a small islands separated by the reversed-current section of the current sheet along the x -axis shown inFig. 10.FIG. 12. Bdy-1, λ = a = m = 1: Level surfaces of ψ in thecase µ = 1 . /a , α = 0 . > α thr , for which Fig. 10 showsthere is no current reversal and hence no magnetic islands. As the poloidal field at the current sheet is proportionalto the current-sheet strength, Eq. (68), such a transitionhas a profound effect on the topology of the magneticsurfaces close to the current sheet, as illustrated in Figs.5 and 6 of Ref. 45 and Figs 11 and 12 of the presentpaper. It is seen that small islands form near the currentsheet in the current-reversal case α < α thr [49], causingthe contour ψ = ψ cut = 0 to trifurcate into upper andlower magnetic surfaces x = x ± ψ ( y | ψ cut ) and the resonantsurface x = 0. In this case the toroidal flux function Φ,Eq. (32), jumps at x = 0 by the amount of “private”toroidal flux in the islands, as shown in Fig. 13. For α > α thr , Φ( x ) is continuous at x .In Fig. 14 the rotational transform, Eq. (34), is plotted, vs. flux surface label. In this plot the amplitude param-eter α = 0 . FIG. 13. Bdy-1, λ = a = m = 1, taking R = a : Toroidalflux vs. the magnetic surface label x , the x -value where a ψ -contour crosses the x -axis, in the case µ = 1 . /a , α =0 . < α thr . The dashed curve is for the unperturbed case α = 0. - - �� � � � - - ι - FIG. 14. Bdy-1, λ = a = m = 1, taking R = a : Rotationaltransform vs. x in the case µ = 1 . /a , α = 0 . > α thr (seetext). The dashed curve is for the unperturbed case α = 0. tational transform to change from the unperturbed res-onant value ι - = 0 at x = 0, jumping from a negativeto a positive value across the current sheet. For smallervalues of α the half-islands remove the discontinuity in ι -.This is because | ∇ ψ | = 0 at the current reversal points,so, from Eq. (33), q diverges as x →
0, locking ι - tozero at x = 0. However, the logarithmic nature of thesingularity means this approach to zero manifests itselfonly at extremely small x , so the slope of ι -( x ) at theorigin is so high that the plots for lower α look, to theeye, qualitatively the same as in Fig. 14 even at the veryfine resolution in x used in this figure. B. Sinusoidal rippled boundary condition
In Fig. 15 we compare boundaries generated by methodBdy-1, described in Sec. IV A, with the corresponding si-nusoidal boundaries defined by Eq. (49). For case (a),small amplitude ripple, method Bdy-1 produces a bound-ary indistinguishable from the target sinusoid, but for ( a )- - - � / � � / � ( b )- - - � / � � / � FIG. 15. Panels (a) and (b) show boundaries generated bymethod Bdy-1 (blue online) for ripple amplitudes α = 0 . .
2, respectively. For comparison, corresponding sinu-soidal boundaries (dashed, orange online), as used in methodBdy-2, are also shown. In both panels, m = 2, aµ = 0 . ( a )- - � / λ - - Δ � ��� α � % ( b )- - � / λ - - - Δ � ��� α � % FIG. 16. Percentage waveform errors using Bdy-2 to fit theprescribed sinusoidal boundary in case m = 2, aµ = 0 . α = 0 .
03; (b) at larger amplitude, α = 0 . larger amplitude, case (b), strong second harmonic erroris clear to the eye.In Fig. 16 we plot the difference between the pure sinu-soid defined by Eq. (49) and boundaries generated by theBdy-2 method, (a) for small-amplitude ripple, α = 0 . α = 0 .
21. The l -sum in Eq. (64)was truncated after l = 3 (hence the dominantly l = 4 er-ror) but even at α = 0 .
21 the percentage waveform errorof 1.5% is tolerable for graphical work. - ��� - ��� - ��� ��� ��� ��� � - ��� - ��������� 〚∂ � ψ〛 FIG. 17. Bdy-2, λ = 2 πa/
2: Plots of the jump in the gradientof ψ , Eq. (77), vs . y/λ for m = 2, µ a = 0 .
2, and the set ofamplitudes α given in the text, showing the occurrence ofcurrent-density reversal for all amplitudes in the set. Figures 17–19 plot the jump in the gradient of6 - ��� - ��� ��� ��� � �������������������� 〚∂ � ψ〛 FIG. 18. Bdy-2, λ = 2 πa/
3: Plots of the jump in the gradientof ψ , Eq. (77), vs . y/λ for m = 3, µ a = 0 .
2, and amplitudesgiven in the text, showing the occurrence of current-densityreversal for all but the highest amplitude. - ��� ��� � ������������������������ 〚∂ � ψ〛 FIG. 19. Bdy-2, λ = 2 πa/
4: Plots of the jump in the gradientof ψ , Eq. (77), vs . y/λ for m = 4, µ a = 0 . ψ , Eq. (77), in the tokamak-relevant case µ a =0 .
2, described after Eq. (46). The seven curves ineach plot are for the set of amplitude values α =0 . , . , . , . , . , . , .
21, whose correspondingvertical-axis intersections run from bottom to top (coloronline). The three figures are for three values of ripplewave number m . The figures reveal the dramatic effectof m on the phenomenon, illustrated in Figs. 10–14, ofthe locking of rotational transform at the resonant valueby half-island formation at small enough α , transition-ing beyond a threshold value of α to the removal of theresonance by the formation of a discontinuity in the ro-tational transform profile at the current sheet interface.The plots show there is qualitative transition in the ro-tational transform profile from strong resonance lockingfor m ≤ m ≥
4. (This is consistent withthe m ≈ α .)We interpret the stronger resonance locking at smaller m as due to the greater penetration of the ripple per-turbation from the boundary [ x = x bdy ( y )] to the in-terface [ x = 0] at longer wavelengths. This is a lin- ear, O ( α ), effect and can easily be seen from the factor1 / sinh κ m a ∼ exp( − κ m a ) in Eq. (77) and the scaling γ S ∝ α , evident from Fig. 9, showing γ S can be ignoredat linear order. On the other hand, it appears the O ( α )d.c. response from the γ S term is not so affected by theexponential decay of exp( − κ m a ) and begins to dominatethe linear response at larger m . VII. CONCLUSION
In this paper we have used numerical calculations togive an exploratory overview of a geometrically simpleapplication of dynamical MRxMHD in the adiabatic ap-proximation, as well as expanding on the general formu-lation in [18] regarding the fundamental question of theappropriate definition for magnetic helicity in MRxMHD.Our calculations have confirmed the physical accessi-bility of the static, equilibrium solutions for RMP reso-nant states found by Loizu, Hudson et al. , [21, 22], es-tablishing the existence of a threshold RMP amplitude atwhich rotational-transform jumps across the resonantlyexcited current sheets occur. Having established this the-oretical basis, the real test will be comparison with exper-iment, for instance revisiting the DIII-D reconstructionof RMP equilibria described in [6, Sec. IV. E ], but usingSPEC with appropriate entropy and magnetic helicityconstraints.While numerical calculations are a good way to explorethe implications of a theory without being tied to par-ticular parameter ranges, for a complete understandingthey need to be complemented by analytical work in ap-propriate asymptotic regimes. In particular, the simple α scalings found empirically in this paper give confidencethat an amplitude expansion could give an adequate un-derstanding of RMP screening dynamics. Such an ex-pansion procedure will be presented elsewhere.Another important area of research for which the HKTmodel is an ideal testbed for MRxMHD are investiga-tions of reconnection mechanisms giving rise to the tran-sition between perturbed states, fully shielded by the res-onantly excited current sheet, to states with fully devel-oped resonant islands. As this involves transfer of mass,entropy and magnetic flux between MRxMHD regions itis related to the field of helicity injection [39–41], whosestudy would require lifting the no-gaps restriction usedin this paper. APPENDICESAppendix A: Magnetic helicity conservation withmoving boundaries
In this Appendix we establish that the gauge constrainton A of conservation of surface loop integrals, mentionedin Sec. II, ensures time-independence of magnetic helic-ity in an ideal plasma region Ω t , with boundary ∂ Ω t de-7pendent on time t . As the guiding principle in formu-lating MRxMHD is to invoke only constraints that arealso appropriate in ideal MHD, this ideal constraint on A is inherited as one of the foundational postulates ofMRxMHD.We need only the no-gaps tangential boundary condi-tion Eq. (4) and the ideal Ohm’s Law E = − v × B , (A1)which, using Maxwell’s equations, gives the equations ofmotion for B and A ,∂ B ∂t = ∇× ( v × B ) , (A2)and ∂ A ∂t = v × B − ∇ ϕ , (A3)where A is a single-valued vector field and ϕ is a scalarpotential, not assumed to be single-valued at this point.Dotting both sides of Eq. (A3) with v and B we find twodifferential equations for ϕ , v ·∇ ϕ = − v · ∂ A ∂t , (A4)and B ·∇ ϕ = − B · ∂ A ∂t . (A5)Differentiating Eq. (3) and using the above assump-tions we find [Supp] the time derivative of the magnetichelicity functional K µ dKdt = (cid:90) ∂ Ω A · B n · v dS + (cid:90) Ω (cid:20) ∂ A ∂t · B + A · ∂ B ∂t (cid:21) dV = ν (cid:88) l =1 (cid:90) S l n · B (cid:74) ϕ (cid:75) dS (A6)It is thus seen that a sufficient condition for invarianceof K is that (cid:74) ϕ (cid:75) ≡
0, i.e. that ϕ be single valued every-where.We now test if single-valuedness is possible withoutcontradicting Galilean invariance and, if so, what restric-tions it places on the gauge of A . We assume the plasmais evolving under Eq. (A2) with a prescribed velocity field v from an initially integrable state with smoothly nestedmagnetic surfaces, which assumption will be hold for afinite time by the frozen-in flux argument [1].Consider first the Galilean invariance problem of a sta-tionary state in the LAB frame as viewed from a movingframe, so the origin of the LAB frame appears to be mov-ing with constant velocity v , hence v = v + v L , wheresubscripts L denote LAB frame fields viewed in the mov-ing frame, and ∂ t in LAB frame maps to D t ≡ ∂ t + v ·∇ in the frame of the observer. For example, Eq. (A2) be-comes D t B L = ∇× ( v L × B L ) , (A7)On the other hand, substituting v = v + v L in Eq. (A2)we find D t B = ∇× ( v L × B ) . (A8)Comparing Eq. (A7) and Eq. (A8) we see that B = B L ,verifying Galilean invariance of B in “pre-Maxwell” idealMHD. (A Lorentz-invariant generalization of helicity hasalso recently been developed [50].)Assuming D t A L = 0 and ϕ L single valued, the LABframe version of Eq. (A3) becomes ∇ ϕ L = v L × B L , (A9)It is easily seen from the LAB version of Eq. (A5) that ϕ L must be constant on each magnetic surface, so, fromthe LAB version of Eq. (A4), v L is, like B , a tangentialfield on each magnetic surface.Substituting v = v + v L in Eq. (A3) and usingEq. (A9) we find ∇ ( ϕ − ϕ L − v · A ) = − D t A . (A10)Taking line integrals of both sides on the magnetic sur-faces around topologically distinct loops C l , moving atthe LAB frame velocity v and cutting the correspond-ing surfaces of section S l , we have (cid:73) C l d l ·∇ ( ϕ − ϕ L − v · A ) = − (cid:73) C l d l · D t A . (A11)Assuming single-valuedness of ϕ (the other terms in theLHS integrand also being single-valued) the loop integralson the LHS vanish. As the contours C l are stationary inthe LAB frame we can commute D t outside the integralon the RHS to find ddt (cid:73) C l d l · A = 0 . (A12)Thus the loop integrals of A on magnetic surfaces aretime-invariant in all frames, which is consistent withGalilean invariance of A : ∇× A = ∇× A L is solved by A = A L + ∇ χ , where χ is an arbitrary but single-valued gauge potential. This confirms that single-valuedness of ϕ is consistent with Galilean invariance of K for systemsthat are stationary in the LAB frame.We now consider systems that are not stationary inany frame, i.e. v is an arbitrary function of timeand space, advecting the magnetic surfaces. Can weshow that single-valuedness of ϕ always implies time-invariance of (cid:72) C l d l · A around magnetic surfaces (partic-ularly the plasma boundary) even for loops not enclosingthe plasma? If so, this is the boundary condition consis-tent with conservation of K .8We begin with the advective [51] form of Eq. (A3) [cf.e.g. eq. (1b) of [36]]. d A dt = ( ∇ A ) · v − ∇ ϕ = − ( ∇ v ) · A − ∇ ( ϕ − v · A ) , (A13)where d/dt ≡ ∂ t + v ·∇ . We also need the advection of aline element d l = r ( r + d l , t ) − r ( r , t ) = d l ·∇ r ( r , t ),where d l is an infinitesimal displacement in the initialposition r of a fluid element some time before the present ddt d l = v ( r ( r + d l , t )) − v ( r ( r , t ))= d l ·∇ r ( r , t ) ∇ v ( r , t ))= d l ·∇ v , (A14)Then ddt (cid:73) C l d l · A = (cid:73) C l (cid:20)(cid:18) ddt d l (cid:19) · A + d l · d A dt (cid:21) = (cid:73) C l { d l · ( ∇ v ) · A − d l · [( ∇ v ) · A + ∇ ( ϕ − v · A )] } = 0 (A15)if and only if ϕ is single valued, which is also the conditionfor K to be time-invariant, so the full helicity constraintcondition in Ω i is conservation of K i and constancy of (cid:72) d l · A around all topologically distinct loops on each dis-joint component of the boundary ∂ Ω i .Note an additional loop integral constraint: If Ω i andΩ j are neighboring regions, the corresponding loop inte-grals (cid:72) ± d l · A on the two sides ± of the common boundaryΩ i,j are constrained to be equal because finiteness of B requires there be vanishing magnetic flux trapped withinthe common interface. Appendix B: Vacuum Helicity
In this Appendix we illustrate the fact that vacuumhelicity is not geometrically invariant by showing it is notinvariant even in slab geometry (if poloidal vacuum fieldis included). This shows that the Finn–Antonsen [41]form of the magnetic helicity (equivalent to the Jensen–Chu [39] relative helicity when there are no gaps in theperfectly conducting boundaries) is not appropriate inMRxMHD.Following [32] we define the harmonic (vacuum) com-ponent B H of a Beltrami field B in an annular toroidas the curl-free ( µ = 0) component carrying the toroidaland poloidal fluxes. Specifically, consider Ω + , for whichthe toroidal ( z -directed) flux is 2 πa F [ F being constantwhen µ = 0, from Eq. (24)] and the poloidal ( y -directed)flux is, from Eq. (31), 2 πψ a R . Then B H = F e z + e z ×∇ ψ H (B1) where the general form of ψ H ( x, y ) such that ∇ ψ = 0and ψ H (0 , y ) = 0 is [cf. Eq. (64) with µ = 0], ψ H ( x, y ) = d H0 | x | + ∞ (cid:88) l =1 d H lm cos lmya sinh (cid:12)(cid:12)(cid:12)(cid:12) lmxa (cid:12)(cid:12)(cid:12)(cid:12) . (B2)with the corresponding vector potential [cf. Eq. (71)] A H = − ψ H e z + e z ×∇ (cid:0) F x (cid:1) . (B3)The coefficients d H are to be chosen so that the Dirichletboundary condition ψ H ( x bdy ( y ) , y ) = ψ a ∀ y (B4)is satisfied, in order to conserve poloidal flux.Then we define the vacuum helicity [39–41] analo-gously to K i , Eq. (3), as [Supp] K H+ ≡ (cid:90) Ω + A H · B H µ dV = (cid:90) Ω + (cid:2) − F ψ H + (cid:0) ∇ F x (cid:1) ·∇ ψ H ) (cid:3) µ dV . (B5)Integration by parts then gives [Supp] K H+ = (cid:90) Ω + (cid:2) − F ψ H + F ∂ x ( xψ H ) (cid:3) µ dV = F µ (cid:34) − (cid:90) Ω + ψ H dV + 2 πa ψ a (cid:35) , (B6)using the area constraint Eq. (37).The term 2 πa ψ a is invariant under changes in x bdy ( y ).However, there appears no reason for the integral (cid:82) Ω + ψ H dV , evaluated from Eq. (B2) as [Supp] (cid:90) Ω + ψ H dV = (cid:90) πa − πa dy (cid:20) d H0 x ( y )+ ∞ (cid:88) l =1 d H lm cos (cid:18) lmya (cid:19) (cid:18) cosh (cid:12)(cid:12)(cid:12)(cid:12) lmx bdy ( y ) a (cid:12)(cid:12)(cid:12)(cid:12) − (cid:19)(cid:21) , (B7)to be invariant in general. As the vacuum helicityEq. (B6) includes this integral we conclude that K H+ isnot in general invariant and therefore not suitable fordefining a relative helicity that is conserved under defor-mations in boundary shape. Appendix C: Scalar Variational Principle
For a variational approach to deriving the Grad–Shafranov form of the Beltrami equation see the onlineSupplement [Supp].9
SUPPLEMENTARY MATERIAL
See online supplementary material [Supp] (provided asan ancillary file in this arXiv version) for further detailsrelevant to this paper: more detailed derivations and Ap-pendix C.
ACKNOWLEDGMENTS
One of the authors (RLD) gratefully acknowledges thesupport of Princeton Plasma Physics Laboratory andThe University of Tokyo during development of the con-cepts in this paper through collaboration visits over anumber of years, and some travel support from Aus-tralian Research Council grant DP11010288. The workof SRH and AB was supported under US DOE grant DE-AC02-09CH11466 and that of ZY was supported underJSPS grant KAKENHI 23224014. The numerical calcu-lations and plots were made using Mathematica 10, [52]. [1] W. A. Newcomb, Annals of Physics , 347 (1958).[2] T. E. Evans, R. A. Moyer, K. H. Burrell, M. E. Fenster-macher, I. Joseph, A. W. Leonard, T. H. Osborne, G. D.Porter, M. J. Schaffer, P. B. Snyder, P. R. Thomas, J. G.Watkins, and W. P. West, Nature Physics , 419 (2006).[3] R. B. White, D. A. Monticello, M. N. Rosenbluth, andB. V. Waddell, Physics of Fluids , 800 (1977).[4] J. B. Taylor, Phys. Rev. Lett. , 1139 (1974).[5] M. J. Hole, S. R. Hudson, and R. L. Dewar, J. PlasmaPhys. , 1167 (2006).[6] S. R. Hudson, R. L. Dewar, G. Dennis, M. J. Hole,M. McGann, G. von Nessi, and S. Lazerson, Phys. Plas-mas , 112502 (2012).[7] S. Lazerson, E. Lazarus, S. Hudson, N. Pablant,and D. Gates, in Proceedings: 39th EPS Con-ference on Plasma Physics and 16th Interna-tional Congress on Plasma Physics Stockholm,Sweden , europhysics conference abstracts, Vol.36F, edited by L. B. S. Ratynskaya and A. Fa-soli (European Physical Society, 2012) p. P4.077http://ocs.ciemat.es/epsicpp2012pap/pdf/P4.077.pdf.[8] S. P. Hirshman and O. Betancourt, J. Comput. Phys. ,99 (1991).[9] J. B. Taylor, Rev. Mod. Phys. , 741 (1986).[10] H. Qin, W. Liu, H. Li, and J. Squire, Phys. Rev. Lett. , 235001 (2012).[11] M. G. Rusbridge, Plasma Physics and Controlled Fusion , 1381 (1991).[12] H. Grad, Phys. Fluids , 137 (1967).[13] H. L. Berk, J. P. Freidberg, X. Llobet, P. J. Morrison,and J. A. Tataronis, Phys. Fluids , 3281 (1986).[14] M. McGann, S. R. Hudson, R. L. Dewar, and G. vonNessi, Phys. Letts. A , 3308 (2010).[15] M. McGann, Hamilton-Jacobi theory for connecting equi-librium magnetic fields across a toroidal surface support-ing a plasma pressure discontinuity , Ph.D. thesis, TheAustralian National University (2013).[16] D. Barmaz,
High- n stability of a pressure discontinuity ina three-dimensional plasma , Master’s thesis, Ecole Poly-technique F´ed´erale de Lausanne and The Australian Uni-versity (2011).[17] I. B. Bernstein, E. A. Frieman, M. D. Kruskal, and R. M.Kulsrud, Proc. Roy. Soc. London Ser. A , 17 (1958).[18] R. L. Dewar, Z. Yoshida, A. Bhattacharjee, and S. R.Hudson, J. Plasma Phys. , 515810604 (2015). [19] A. Reiman, N. M. Ferraro, A. Turnbull, J. K. Park,A. Cerfon, T. E. Evans, M. J. Lanctot, E. A. Lazarus,Y. Liu, G. McFadden, D. Monticello, and Y. Suzuki,Nucl. Fusion , 063026 (2015).[20] J. Loizu, S. Hudson, A. Bhattacharjee, and P. Helander,Phys. Plasmas , 022501 (2015).[21] J. Loizu, S. R. Hudson, A. Bhattacharjee, S. Lazerson,and P. Helander, Phys. Plasmas , 090704 (2015).[22] J. Loizu, S. R. Hudson, P. Helander, S. A. Lazerson, andA. Bhattacharjee, Phys. Plasmas , 055703 (2016).[23] S. A. Lazerson, J. Loizu, S. Hirshman, and S. R. Hudson,Physics of Plasmas , 012507 (2016).[24] T. S. Hahm and R. M. Kulsrud, Phys. Fluids , 2412(1985).[25] Y.-M. Huang, A. Bhattacharjee, and A. H. Boozer, As-trophys. J. , 106 (2014).[26] W. A. Newcomb, Annals of Physics , 232 (1960).[27] G. O. Spies, Phys. Plasmas , 3030 (2003).[28] R. Mills, M. J. Hole, and R. L. Dewar, J. Plasma Phys. , 637 (2009).[29] A. H. Boozer and N. Pomphrey, Phys. Plasmas ,110707 (2010).[30] R. B. White, Phys. Plasmas , 022105 (2013).[31] Y. Zhou, Y.-M. Huang, H. Qin, and A. Bhattacharjee,Phys. Rev. E , 023205 (2016).[32] Z. Yoshida and R. L. Dewar, J. Phys. A: Math. Gen. ,365502 (2012).[33] R. L. Dewar, L. H. Tuen, and M. J. Hole, Plasma Phys.Control. Fusion , 044009 (2017).[34] R. D. Parker and R. L. Dewar, Computer Physics Com-munications , 1 (1990).[35] W. A. Newcomb, Proceedings of the International Con-ference on Plasma Physics and Controlled Nuclear Fu-sion Research 1961 , Nucl. Fusion Suppl.