An adjoint method for determining the sensitivity of island size to magnetic field variations
UUnder consideration for publication in J. Plasma Phys. An adjoint method for determining thesensitivity of island size to magnetic fieldvariations
Alessandro Geraldini † , M. Landreman , E. Paul Institute for Research in Electronics and Applied Physics, University of Maryland, CollegePark, MD 20742, USA Swiss Plasma Center, ´Ecole Polytechnique F´ed´erale de Lausanne, CH-1015 Lausanne,Switzerland Department of Astrophysical Sciences, Princeton University, Princeton, NJ 08544, USA(Received xx; revised xx; accepted xx)
An adjoint method to calculate the gradient of island width in stellarators is presentedand applied to a set of magnetic field configurations. The underlying method of calcula-tion of the island width is that of Cary & Hanson (1991) (with a minor modification),and requires that the residue of the island centre be small. Therefore, the gradient ofthe residue is calculated in addition. Both the island width and the gradient calculationsare verified using an analytical magnetic field configuration introduced in Reiman &Greenside (1986). The method is also applied to the calculation of the shape gradient ofthe width of a magnetic island in an NCSX vacuum configuration with respect to positionson a coil. A gradient-based optimization is applied to a magnetic field configurationstudied in Hanson & Cary (1984) to minimize stochasticity by adding perturbations toa pair of helical coils. Although only vacuum magnetic fields and an analytical magneticfield model are considered in this work, the adjoint calculation of the island widthgradient could also be applied to a magnetohydrodynamic (MHD) equilibrium if thederivative of the magnetic field with respect to the equilibrium parameters was known.Using the island width gradient calculation presented here, more general gradient-basedoptimization methods can be applied to design stellarators with small magnetic islands.Moreoever, the sensitivity of the island size may itself be optimized to ensure that coiltolerances with respect to island size are kept as high as possible.
1. Introduction
Stellarators (Spitzer Jr 1958) are promising candidates for a nuclear fusion devicewhose main advantage is to operate in an intrinsically steady state (Helander 2014). Inorder to avoid the need of a toroidal plasma current to produce a poloidal magnetic field,stellarators lack the continuous toroidal symmetry of the magnetic field vector which isa characteristic of the tokamak. Contrary to the tokamaks, however, stellarators havemagnetic fields that tend to be non-integrable and to develop magnetic islands whichbreak the otherwise nested toroidal magnetic surfaces (Rosenbluth et al. et al. † Email address for correspondence: [email protected] a r X i v : . [ phy s i c s . p l a s m - ph ] F e b A. Geraldini, M. Landreman, E. Paul degree of accuracy (Pedersen et al. et al. et al. et al. (1993) highlighted that the assumption of flux surfaces in a stellarator experiment isincorrect. Error fields have been studied in the CNT stellarator configuration (Hammond et al. et al. (2018) conclude that tools to gauge thesensitivity of island size are important as they allow trim coils to be used to reduce thevolume of magnetic islands. A recent paper by Zhu et al. (2019) addresses the issue of theidentification and removal of the error fields responsible for island size using a Hessianmatrix approach and making the simplifying assumptions that first derivatives are zeroand that variations in the magnetic coordinates can be ignored.There are several methods to calculate the width of a magnetic island. The most basicapproach relies on making a detailed scatter plot (Poincar´e plot) of the position at whicha magnetic field line intersects a given poloidal plane, and then measuring the island sizefrom the plot. This, however, is extremely time consuming and noisy, and is thereforeespecially inadequate if one wants to calculate the island width of a very large numberof configurations in a short time, let alone if one wants to obtain gradient information.Variations of this method employ automated algorithms to detect islands and calculatethe island width from integration of several magnetic field lines in an island (Pedersen et al. et al. (1990) and exploitsa Fourier decomposition of the magnetic field vector. In this work, we consider a measureof island width derived by Cary & Hanson (1991) that allows an efficient computation ofthe width of a magnetic island and also allows for the direct calculation of its gradient.This approach exploits the small-island approximation to calculate a measure of widththat depends only on the magnetic field line corresponding to the island centre and onthe equations for linearized displacements from the island centre.Adjoint methods have recently been applied in stellarators to obtain derivatives ofneoclassical fluxes (Paul et al. et al. et al. djoint method for sensitivity of island size to magnetic field variations
Reiman model ;(ii) the magnetic field in NCSX, specified by the position of a set of discrete points ona set of filamentary coils and by the current through each coil;(iii) a magnetic field produced by a pair of helical coils that was optimized in Hanson& Cary (1984) and Cary & Hanson (1986).A potential application of the gradient calculations is the fast calculation of coiltolerances with respect to island size. Another important application of the gradientcalculations developed here is optimization of stellarator surfaces. In order to minimizestochasticity and island size in stellarator vacuum magnetic fields, methods employedso far often minimize the magnitude of a quantity known as the residue (Greene 1968)of periodic field lines (Hanson & Cary 1984; Cary & Hanson 1986). This quantity iscalculated by linearising the equations for the magnetic field line about the island centreand calculating a matrix known as the full orbit tangent map. This is a linear maprelating the displacement of a magnetic field line from a nearby periodic field line, aftera full magnetic field line period, to the initial displacement. In general, this map isa two-dimensional matrix with unit determinant, and therefore has three degrees offreedom. The residue is related to the trace of this map, and provides a criterion fordetermining whether the closed field line is an O point (island centre) or an X point. Ifthe residue is zero there is no island chain (an unbroken rational flux surface), and if theresidue is small the size of the island chain is small compared to the length scale of themagnetic configuration. Therefore, the residue constitutes an extremely useful degree offreedom of the map. Minimizing the absolute value of the residue of a periodic field lineamounts to reducing the stochasticity in the magnetic field configuration and eventuallyalso reducing the volume occupied by the corresponding island chain in the magneticfield configuration. The gradient of the residue is used to find an optimal magnetic fieldconfiguration with small islands for a helical coil configuration previously optimized inCary & Hanson (1986).An aspect of the problem that is not considered in this work is the application tomagnetohydrodynamic (MHD) equilibrium configurations (Hegna 2012; Hudson et al.
2. Calculation of island size from periodic magnetic field trajectory
In this section, we review the calculation method of island widths introduced in Cary& Hanson (1991). In section 2.1 we assume the existence of toroidal flux surfaces and,upon considering an island-producing perturbation, we derive the equations for magneticfield lines in an island chain in the magnetic coordinates of the unperturbed system. Thelinearized motion near the island centre (O point) is analyzed in section 2.2 and the
A. Geraldini, M. Landreman, E. Paul equation for the displacement from the island centre is expressed in a frame rotatingwith the island centre poloidally around the magnetic axis. The equations describingmagnetic field line trajectories in cylindrical coordinates are obtained in section 2.3.Using the results of sections 2.2 and 2.3, in section 2.4 an expression for the island widthis obtained. 2.1.
Magnetic coordinates
It is often convenient to use magnetic coordinates (Helander 2014) when describingthe position along a magnetic field in systems with nested flux surfaces, such as the idealmagnetic configuration in fusion devices. Flux surfaces are closed toroidal surfaces wherethe enclosed toroidal magnetic flux 2 πψ through a surface of constant toroidal angle ϕ and the enclosed poloidal magnetic flux 2 πχ through a surface of constant poloidal angle θ are constant. The choice of magnetic coordinates is not unique: here we choose ϕ to bethe geometric toroidal angle which also constrains the poloidal angle θ (not generally ageometric angle). The magnetic field is given by B = ∇ ψ × ∇ θ + ∇ ϕ × ∇ χ ( ψ, θ, ϕ ) . (2.1)With nested flux surfaces, the poloidal flux χ is always a function of the toroidal flux ψ only, since both these quantities are constant on each flux surface. Then, ∇ χ and ∇ ψ are parallel to one another and the magnetic field line trajectory never crosses the fluxsurfaces since B · ∇ χ = B · ∇ ψ = 0.To study the small departure from the ideal configuration with nested flux surfaces, weuse magnetic coordinates ψ , θ and ϕ which correspond to the toroidal flux, poloidal angleand geometric toroidal angle of the unperturbed system. Equation (2.1) still describesthe magnetic field in such a system, but the function χ is comprised of a large piecewhich is equal to the poloidal flux of the unperturbed nested flux surfaces, χ ( ψ ), and asmall perturbation that breaks the flux surfaces, χ ( ψ, θ, ϕ ), χ ( ψ, θ, ϕ ) = χ ( ψ ) + χ ( ψ, θ, ϕ ). (2.2)As shown in appendix A, the magnetic field line trajectory is a Hamiltonian systemwhere the canonical co-ordinate q is θ , the canonical momentum p is ψ , the Hamiltonian H is χ and the time t is ϕ (Cary & Littlejohn 1983). Hamilton’s equations are thereforegiven by dψ/dϕ = − ∂χ/∂θ and dθ/dϕ = ∂χ/∂ψ . Since χ is small, to lowest order inthe perturbation the magnetic field line trajectory lies in the unperturbed nested fluxsurfaces, dψdϕ (cid:39)
0, (2.3) dθdϕ (cid:39) χ (cid:48) ( ψ ) ≡ ι ( ψ ). (2.4)Here, we have defined the rotational transform ι ( ψ ), which corresponds to the averagenumber of poloidal turns of a magnetic field line around the magnetic axis divided by theaverage number of toroidal turns. The magnetic field trajectory in magnetic coordinatesis therefore approximately given by ψ (cid:39) ψ (where ψ is a constant) and θ (cid:39) ι ( ψ ) ϕ + θ i ,where θ i is a constant.To calculate the effect of the perturbation to the function χ on the magnetic field linetrajectory, we re-express χ as a sum of its Fourier components, χ ( ψ, θ, ϕ ) = (cid:88) m,n χ m,n ( ψ ) exp ( imθ − inϕ ) , (2.5) djoint method for sensitivity of island size to magnetic field variations χ is periodic in θ and ϕ . Here, n is the toroidal mode numberand m is the poloidal mode number of the perturbation. As shown in appendix B, for anyunperturbed flux surface χ ( ψ ), the effect of χ is dominated by the pair of Fourier modesthat resonate with the rotational transform of the unperturbed flux surface, n/m = ι ( ψ )(Helander 2014; Cary & Littlejohn 1983). Thus, the effect of the perturbation is largestat rational flux surfaces, where the rotational transform is a rational number. In whatfollows, we assume that resonances at different rational flux surfaces do not overlap andinteract with each other and that higher order harmonics of the resonances have a smalleramplitude and can be neglected. Consider the rational flux surface where ι ( ψ ) = NM . (2.6)Here, N is the toroidal mode number and M is the poloidal mode number of the island-producing perturbation, which is readily re-expressed to χ ( ψ, θ, ϕ ) = (cid:15) ( ψ ) cos ( M ζ ( ψ ) + M θ − N ϕ ) . (2.7)An amplitude, (cid:15) ( ψ ) >
0, and a phase factor, ζ ( ψ ), of the resonant perturbation have beenintroduced to replace the pair of complex amplitudes χ M,N and χ − M, − N correspondingto the terms in χ that resonate with ι ( ψ ). The subscripts M and N on (cid:15) and ζ havebeen omitted for brevity. Introducing a new variable Θ , Θ = θ − N ϕM , (2.8) χ is re-expressed to χ ( ψ, Θ ) = (cid:15) ( ψ ) cos ( M ζ ( ψ ) + M Θ ) . (2.9)With this change of variable, the dependence on ϕ of the phase of χ has dropped.However, the function χ is no longer a Hamiltonian in the variables ( ψ, Θ ). As shown inappendix A.1, the Hamiltonian in the new variables is K ( ψ, Θ ) = χ ( ψ, Θ ) − ι ( ψ ) ψ , (2.10)and Hamilton’s equations are dΘ/dϕ = ∂K/∂ψ and dψ/dϕ = − ∂K/∂Θ . The Hamilto-nian K is independent of ϕ and is therefore conserved following the perturbed magneticfield lines.In what follows, K ( ψ, Θ ) is expanded close to the rational flux surface. The perturba-tion in ψ near the rational flux surface is ψ = ψ − ψ . (2.11)The function χ is expanded to χ ( ψ , Θ ) = χ ( ψ ) + χ (cid:48) ( ψ ) ψ + 12 χ (cid:48)(cid:48) ( ψ ) ψ + (cid:15) ( ψ ) cos ( M ζ ( ψ ) + M Θ )+ (cid:15) (cid:48) ( ψ ) ψ cos ( M ζ ( ψ ) + M Θ ) − M (cid:15) ( ψ ) ζ (cid:48) ( ψ ) ψ sin ( M ζ ( ψ ) + M Θ )+ O ( (cid:15) ˆ ψ , ψ ˆ ψ ), (2.12)where the normalized ψ perturbation is ˆ ψ = ψ ι (cid:48) ( ψ ) (equivalent to the ψ perturbationdivided by a measure of typical variations of ψ and χ across the toroidal configuration).Using equations (2.12), (2.10), and ι = χ (cid:48) ( ψ ), the Hamiltonian K becomes K ( ψ , Θ ) = χ ( ψ ) − ι ( ψ ) ψ + 12 ι (cid:48) ( ψ ) ψ + (cid:15) ( ψ ) cos ( M ζ ( ψ ) + M Θ ) A. Geraldini, M. Landreman, E. Paul + (cid:15) (cid:48) ( ψ ) ψ cos ( M ζ ( ψ ) + M Θ ) − M (cid:15) ( ψ ) ζ (cid:48) ( ψ ) ψ sin ( M ζ ( ψ ) + M Θ )+ O ( (cid:15) ˆ ψ , ψ ˆ ψ ). (2.13)Hamilton’s equations for the magnetic field line sufficiently close to the rational fluxsurface are thus dΘdϕ = ∂K∂ψ = ι (cid:48) ( ψ ) ψ + (cid:15) (cid:48) ( ψ ) cos ( M ζ ( ψ ) + M Θ ) − (cid:15) ( ψ ) M ζ (cid:48) ( ψ ) sin ( M ζ ( ψ ) + M Θ ) + O (ˆ (cid:15) ˆ ψ , ψ ˆ ψ ), (2.14)and dψ dϕ = − ∂K∂Θ = (cid:15) ( ψ ) M sin ( M ζ ( ψ ) + M Θ ) + O ( (cid:15) ˆ ψ ), (2.15)where ˆ (cid:15) = ι (cid:48) ( ψ ) (cid:15) ( ψ ).The fixed points of the magnetic field line flow in the ( ψ , Θ ) coordinates, whichrepresent closed magnetic field lines, occur at dΘ/dϕ = dψ /dϕ = 0. This correspondsto ( ψ , Θ ) = ( ¯ ψ , ¯ Θ ) such that sin (cid:0) M ζ ( ψ ) + M ¯ Θ (cid:1) = 0 and ¯ ψ = ± (cid:15) (cid:48) ( ψ ) /ι (cid:48) ( ψ ) = O ( (cid:15) ), where we have used cos (cid:0) M ζ ( ψ ) + M ¯ Θ (cid:1) = ∓
1. Introducing magnetic coordinatesrelative to a closed magnetic field line, δψ = ψ − ψ − ¯ ψ and δΘ = Θ − ¯ Θ , the equations(2.14) and (2.15) linearized near the closed magnetic field give ddϕ (cid:18) δψδΘ (cid:19) = (cid:18) O (ˆ (cid:15) ) ∓ (cid:15) ( ψ ) M + O ( (cid:15) ˆ ψ ) ι (cid:48) ( ψ ) + O ( ι (cid:48) ( ψ )ˆ (cid:15) ) O (ˆ (cid:15) ) (cid:19) (cid:18) δψδΘ (cid:19) . (2.16)As will be shown explicitly by expanding the Hamiltonian near the closed field line,the error terms in (2.16) arising from the diagonal elements of the matrix are negligiblebecause the only self-consistent ordering relating the characteristic sizes of δψ and δΘ is | ι (cid:48) | δψ ∼ M √ ˆ (cid:15)δΘ . Trajectories neighbouring the island centre (O point) are describedwhen the signs of the off-diagonal elements in the matrix in (2.16) are opposite, since inthis case the eigenvalues are purely imaginary and trajectories are periodic. The othersign choice corresponds to the trajectories passing close to the crossing point of the islandseparatrices (X point). Thus, we can replace ± everywhere with ± sgn ( ι (cid:48) ( ψ )), with thetop sign and bottom signs understood to correspond to O and X points respectively. Themotion sufficiently close to the centre is thus always described by ddϕ (cid:18) δψδΘ (cid:19) = (cid:18) O (ˆ (cid:15) ) − sgn ( ι (cid:48) ( ψ )) (cid:15) ( ψ ) M (1 + O (ˆ (cid:15) )) ι (cid:48) ( ψ )(1 + O (ˆ (cid:15) )) O (ˆ (cid:15) ) (cid:19) (cid:18) δψδΘ (cid:19) , (2.17)where we have set cos (cid:0) M ζ ( ψ ) + M ¯ Θ (cid:1) = − sgn ( ι (cid:48) ( ψ )) at an O point. Solving this linearsystem gives (cid:18) δψ ( ϕ ) δΘ ( ϕ ) (cid:19) = T (cid:18) δψ (0) δΘ (0) (cid:19) , (2.18)where the tangent map T is given by T (cid:39) (cid:18) cos( ωϕ ) − ( ω/ι (cid:48) ( ψ )) sin( ωϕ )( ι (cid:48) ( ψ ) /ω ) sin( ωϕ ) cos( ωϕ ) (cid:19) , (2.19)and the frequency at which neighbouring points rotate around the O point is ω (cid:39) M (cid:113) | ι (cid:48) ( ψ ) | (cid:15) ( ψ ) . (2.20)Note that, from equation (2.17), each element of the tangent map in (2.19) has an error djoint method for sensitivity of island size to magnetic field variations O (ˆ (cid:15) ) and the the frequency ω in (2.20) has an error of O (cid:16) ˆ (cid:15) (cid:112) | ι (cid:48) ( ψ ) | (cid:15) ( ψ ) (cid:17) (Cary &Hanson 1991).In order to relate the linearized motion along a field line neighbouring the islandcentre described by equations (2.18)-(2.20) to the island width, the motion along theisland separatrix is studied. The Hamiltonian on the separatrix is constant and equal toits value at the X point, K ( ¯ ψ, ¯ Θ ) = χ ( ψ ) − ι ( ψ ) ψ + sgn ( ι (cid:48) ( ψ )) (cid:15) ( ψ ) + O(ˆ (cid:15)(cid:15) ), (2.21)where we inserted ¯ ψ = − (cid:15) (cid:48) ( ψ ) / | ι (cid:48) ( ψ ) | (X point), sin (cid:0) M ζ ( ψ ) + M ¯ Θ (cid:1) = 0 (fixed point)and cos (cid:0) M ζ ( ψ ) + M ¯ Θ (cid:1) = sgn ( ι (cid:48) ( ψ )) (X point) in equation (2.13). From equation(2.13) and (cid:15) ( ψ ) >
0, the value of ψ is largest when cos ( M ζ ( ψ ) + M Θ ) = − sgn ( ι (cid:48) ( ψ ))on the separatrix, and so by evaluating K ( ψ, Θ ) at this point and equating it to equation(2.21) we obtain 12 | ι (cid:48) ( ψ ) | ψ = 2 (cid:15) ( ψ ) + O ( (cid:15) ˆ ψ , ˆ (cid:15)(cid:15) ). (2.22)Hence, the values of ψ at the separatrix at the point where the island is largest are ψ = ± (cid:112) (cid:15) ( ψ ) / | ι (cid:48) ( ψ ) | + O ( (cid:15) ) ∼ √ ˆ (cid:15)/ι (cid:48) ( ψ ), and the full island width, denoted Υ , is Υ = 4 (cid:115) (cid:15) ( ψ ) | ι (cid:48) ( ψ ) | . (2.23)Note that Υ is approximately equal to the full-island width in the magnetic coordinate ψ with an absolute error of O ( (cid:15) ), equivalent to a relative error of O (cid:0) ˆ (cid:15) / (cid:1) . Note also that | ¯ ψ | (cid:28) Υ , and so the island centre is equidistant from the two branches of the separatrixto lowest order in ˆ (cid:15) . From (2.20) and (2.23), the matrix T in (2.19) can be rewritten with Υ in place of ι (cid:48) to obtain T (cid:39) (cid:18) cos( ωϕ ) − sgn ( ι (cid:48) ( ψ )) ( M Υ/
4) sin( ωϕ )sgn ( ι (cid:48) ( ψ )) (4 /M Υ ) sin( ωϕ ) cos( ωϕ ) (cid:19) . (2.24)When linearising the equations for the magnetic field line near the island centre, theassociated Hamiltonian is equation (2.10) expanded near the island centre up to quadraticterms in δψ and δΘ . Retaining only the lowest order terms in (cid:15) gives K =const + 12 ι (cid:48) ( ψ ) δψ + 12 sgn ( ι (cid:48) ( ψ )) M (cid:15)δΘ + O (ˆ (cid:15)δψ ι (cid:48) ( ψ ) , M ˆ (cid:15)δψδΘ, M ˆ (cid:15)(cid:15)δΘ, M ˆ (cid:15)(cid:15)δΘ ), (2.25)where we have deduced that the only self-consistent ordering is ι (cid:48) ( ψ ) δψ ∼ M √ ˆ (cid:15)δΘ .Note that equation (2.25) can be made to be exactly quadratic: the linear term in theerror in (2.25) could be cancelled by calculating the first order correction in ˆ (cid:15) of the valueof Θ at the island centre and redefining δΘ relative to this more accurate coordinate.The neglected O ( M ˆ (cid:15)δψδΘ ) cross-term in (2.25) gives rise to the diagonal terms inthe matrix in equation (2.16), which have a negligible contribution once the ordering ι (cid:48) ( ψ ) δψ ∼ M √ ˆ (cid:15)δΘ is taken into account. The magnitude of the quadratic perturbationto the Hamiltonian, | δK | = δ u · K · δ u = (cid:0) δψ, δΘ (cid:1) (cid:18) | ι (cid:48) ( ψ ) | + O (ˆ (cid:15)ι (cid:48) ( ψ )) O ( M ˆ (cid:15) ) O ( M ˆ (cid:15) ) M (cid:15) + O ( M ˆ (cid:15)(cid:15) ) (cid:19) (cid:18) δψδΘ (cid:19) ,(2.26)is a scalar invariant, as it is conserved following the field lines neighbouring the O point. A. Geraldini, M. Landreman, E. Paul
To lowest order in ˆ (cid:15) , K is diagonal and thus the trajectories infinitesimally close to theisland centre in magnetic coordinates are ellipses that are approximately aligned with themagnetic coordinate directions and elongated in the Θ direction. The angle between thecharacteristic directions of the ellipses and the magnetic coordinate axes is small, O (ˆ (cid:15) ).Note that if we diagonalized the matrix K exactly, we would obtain additional error termsof the same order in the diagonal terms.2.2. Relating magnetic coordinates to lengths at the island centre
The relationship between the displacement from the island centre measured as a lengthin the poloidal cross-section at a given ϕ and the same displacement measured in magneticcoordinates must be linear if the displacement is infinitesimally small (as the relationshipbetween the two sets of coordinates must be locally described by a Taylor expansion).Thus, we define the local linear change of variables δ u = ( δψ, δΘ ) → δ ξ = ( δξ ⊥ , δξ (cid:107) ),where δξ ⊥ ( ϕ ) and δξ (cid:107) ( ϕ ) are displacements from the island centre in two orthogonaldirections (yet unspecified) in the poloidal plane with toroidal angle ϕ . The two sets ofcoordinates are related by δ u ( ϕ ) = Q ( ϕ ) · δ ξ ( ϕ ) for any ϕ , where Q ( ϕ ) = (cid:32) ∂ψ∂ξ ⊥ ( ϕ ) ∂ψ∂ξ (cid:107) ( ϕ ) ∂Θ∂ξ ⊥ ( ϕ ) ∂Θ∂ξ (cid:107) ( ϕ ) (cid:33) . (2.27)For each value of ϕ , the scalar invariant can be cast in the new coordinates, | δK | = δ ξ ( ϕ ) · Q (cid:124) ( ϕ ) · K · Q ( ϕ ) · δ ξ ( ϕ ), where Q (cid:124) denotes the transpose of Q and the locallaboratory invariant matrix Q (cid:124) ( ϕ ) · K · Q ( ϕ ) is symmetric since K is symmetric. Thus,the eigenvectors of the local laboratory invariant matrix are orthogonal and the pointsof intersection of a trajectory, infinitesimally close to the island centre, with any givenpoloidal plane ϕ lie on an ellipse.The dependence of the scalar invariant on ϕ expresses the fact that, when viewingthe motion continuously in different ϕ -planes, the poloidal displacement of a field lineinfinitesimally close to the island centre lies on a continuously varying ellipse . Thisis a consequence of the fact that in a stellarator the poloidal cross-section of fluxsurfaces taken at different toroidal angles is generally different, and thus the shape ofthe island continuously changes and rotates poloidally as the island centre is followedaround. Nonetheless, a set of equivalent flux surface sections always exists for values of ϕ that differ by an integer multiple of a field period, 2 π/n , where n (cid:62)
1. The fluxsurface sections are, in general, only equivalent to lowest order in the island-producingperturbation, (cid:15) , because an island chain may break the field periodicity ( n = 1 is aspecial case where the field periodicity is never broken). The closed field line intersectsthe set of approximately equivalent poloidal planes a finite number of times, L , beforereturning to the original position. Therefore, when snapshots of the position along afield line neighbouring the island centre are taken at ϕ = ϕ k + 2 πQL/n for any positiveinteger Q and initial toroidal angle ϕ k , the motion appears to be around the same ellipse .In this work, we choose one set of an infinite number of possible sets of symmetry planesby specifying the toroidal plane corresponding to ϕ = 0 and considering the set ofsymmetry planes given by ϕ k = 2 πk/n , where k is a positive integer. For stellaratorswith stellarator symmetry, the magnetic configuration in the plane ϕ = 0 is chosen to beup-down symmetric. Henceforth, any quantity that is a function of ϕ k will be denotedwith a subscript k , e.g. Q k = Q ( ϕ k ).As the matrix K is (approximately) diagonal, the diagonalization of the invariantmatrix Q (cid:124) k · K · Q k is achieved (approximately) by choosing the coordinates δξ (cid:107) ,k and djoint method for sensitivity of island size to magnetic field variations δξ ⊥ ,k such that ∂ψ/∂ξ (cid:107) ,k = ∂Θ/∂ξ ⊥ ,k = 0 for all k and thus δψ k = ∂ψ∂ξ ⊥ ,k δξ ⊥ ,k , (2.28) δΘ k = ∂Θ∂ξ (cid:107) ,k δξ (cid:107) ,k . (2.29)With this choice, the coordinates δξ ⊥ ,k and δξ (cid:107) ,k quantify the displacement (as a length)from the island centre in the poloidal plane ϕ = 2 πk/n , measured in the directionsassociated with δψ (across the flux surface) and δΘ (along the flux surface), respectively.In the new coordinates, the displacement from the fixed point satisfies the equation δ ξ k + q ( ϕ ) = S qR,k · δ ξ k , (2.30)where the tangent map in the rotating frame, S qR,k = Q − k + q · T k + q · Q k , is given by S qR,k (cid:39) ∂ξ ⊥ ,k + q ∂ψ ∂ψ∂ξ ⊥ ,k cos (cid:16) πωqn (cid:17) − ∂ξ ⊥ ,k + q ∂ψ ∂Θ∂ξ (cid:107) ,k MΥ sin (cid:16) πωqn (cid:17) ∂ξ (cid:107) ,k + q ∂Θ ∂ψ∂ξ ⊥ ,k MΥ sin (cid:16) πωqn (cid:17) ∂ξ (cid:107) ,k + q ∂Θ ∂Θ∂ξ (cid:107) ,k cos (cid:16) πωqn (cid:17) . (2.31)Note that ¯ w ⊥ ,k = Υ ∂ξ ⊥ ,k ∂ψ (2.32)is approximately equal to the island width measured at ϕ = 2 πk/n , with an absoluteerror of O ( Υ ∂ ξ ⊥ ,k /∂ψ ) giving rise to a relative error of O (ˆ (cid:15) / ). This relative error is tobe added to the equally large error in approximating the island width in the coordinate ψ as Υ in equation (2.23). Using equation (2.32), the off-diagonal elements of S qR,k explicitlyinclude the island width. Upon imposing δξ (cid:107) ,k = 0 as an initial condition, we obtain theequation δξ (cid:107) ,k + q δξ ⊥ ,k (cid:39) ∂ξ (cid:107) ,k + q ∂Θ M ¯ w ⊥ ,k sin (cid:18) πωqn (cid:19) , (2.33)where only the bottom-left element of (2.31) appears.Since the derivative in ∂ξ (cid:107) ,k + q /∂Θ is taken at fixed ϕ , the definition of Θ in (2.8)implies that ∂ξ (cid:107) ,k + q ∂θ = ∂ξ (cid:107) ,k + q ∂Θ . (2.34)Integrating in θ at fixed ϕ = ϕ k + q gives the circumference ¯ C around the unperturbedflux surface, ¯ C = (cid:90) π − π ∂ξ (cid:107) ,k + q ∂θ dθ. (2.35)Since each plane under consideration, given by ϕ = 2 π ( k + q ) /n , has the same unper-turbed flux surface indepently of the value of k + q , the circumference is independent ofthe index and thus ¯ C = ¯ C k + q . Consider the definition of Θ in equation (2.8). Followingthe island centre Θ k + q = Θ k = θ k − N ϕ k /M is conserved and therefore the poloidalangle changes according to the equation θ k + q − θ k = 2 πN q/ ( n M ). Since we are onlyconsidering the set of equivalent planes separated by intervals of 2 π/n in ϕ , the poloidalangles θ k + q correspond to the same flux surface poloidal cross-section. The field line firstreturns to the same poloidal location in an equivalent plane when N L/ ( n M ) = ¯ n , where0 A. Geraldini, M. Landreman, E. Paul L and ¯ n are the smallest possible integers satisfying this relation. Thus, the number ofdistinct islands crossed by a unique island centre magnetic field line in an equivalentpoloidal plane is L = ¯ nn N M. (2.36)The interval in poloidal angle θ , when defined between − π and π , between the fixedpoints in a given plane is 2 π/L , as there are L equally spaced fixed points. Thus, if L is sufficiently large, L (cid:29)
1, the integral in (2.35) can be replaced by a sum and thecircumference can be approximated by¯ C (cid:39) πL q + L − (cid:88) q = q ∂ξ (cid:107) ,k + q ∂θ (cid:39) q + L − (cid:88) q = q π ( δξ (cid:107) ,k + q /δξ ⊥ ,k ) M ¯ w ⊥ ,k L sin (2 πωq/n ) . (2.37)The error made in approximating the integral over the periodic function θ as a sum withequally spaced integration points is exponentially small in 1 /L , O (exp( − L )).The integer q can be chosen arbitrarily, although a convenient choice is made asfollows. If sin (2 πω ( k + q ) /n ) = 0, the denominator in the summand with index q is zero.However, no matter how precise, a numerical calculation of δξ (cid:107) ,k + q will not give exactlyzero, due to the fact that the elliptical motion described using the tangent map (2.24) hassmall errors both in the aspect ratio and in the axes directions (the higher order termsin ˆ (cid:15) that were neglected). To avoid the resulting divergence of the error, the first indexof the summand is chosen such that the the sum is centred around sin (2 πωq/n ) = 1,2 πω ( q + L/ /n = π/
2, giving q = int (cid:18) − L n ω (cid:19) , (2.38)where int denotes a function that rounds to the nearest integer. This amounts to followingthe linearized trajectory for a large number, q /L , of closed magnetic field line periods.With this choice, the linearized trajectory is followed until it has rotated by as close aspossible to π/ q in the middle of the summationinterval. Moreover, we assume that the frequency of rotation, ω in equation (2.20), issmall enough that the change in the quantity 2 πω ( k + q ) /n in the summation interval q (cid:54) q < q + L is also small, 2 πωL/n (cid:28)
1, giving sin (2 πωq/n ) (cid:39) q in the sum and thus C (cid:39) M π L q + L − (cid:88) q = q δξ (cid:107) ,k + q δξ ⊥ ,k ¯ w ⊥ ,k . (2.39)In Cary & Hanson (1991), this assumption is not made and the sine function in equation(2.37) is kept. Using equations (2.20) and the pessimistic ordering in which L/n islargest, L ∼ n M , this assumption requires the modestly stricter ordering 2 πM (cid:112) ι (cid:48) (cid:15) ∼ πM √ ˆ (cid:15) (cid:28) † † Cary & Hanson (1991) point out that — according to Greene (1979) — the rotation angleof points near the stable fixed point of a standard map (that resembles the full-orbit map M k considered here) has a threshold of about 60 ◦ above which chaos ensues. Therefore, if 2 πωL/n is not small, it is likely that the concept of island width has broken down. The footnote in page299 of Rosenbluth et al. (1966) argues that (cid:15) must decrease with M more rapidly than M − inorder to have non-overlapping islands, further justifying the stricter ordering. djoint method for sensitivity of island size to magnetic field variations Cylindrical coordinates
Since stellarators are toroidal devices, we choose the right-handed cylindrical coordi-nates (
R, ϕ, Z ) as a convenient fixed coordinate system. Here ϕ is the toroidal coordinate, R is the smallest distance of a point from the axis through the centre of the torus and Z is the smallest distance of a point from the mid-plane of the device.In what follows, we study the equations describing the magnetic field line poloidalposition, X = ( R, Z ), as a function of toroidal angle ϕ . Considering the magnetic fieldline as a streamline of the flow field B , the streamline trajectory satisfies RdϕB ϕ = dRB R = dZB Z , (2.40)and thus d X dϕ = V ( X , ϕ ) ≡ R B p ( X , ϕ ) B ϕ ( X , ϕ ) , (2.41)where we have denoted the component of the magnetic field vector in the poloidal planeas B p = ( B R , B Z ). To obtain the position of a magnetic field line at ϕ = ϕ k + q , weintegrate equation (2.41) in ϕ from an initial poloidal position X k = X ( ϕ k ), X k + q = F qk ( X k ) ≡ X k + (cid:90) ϕ k + q ϕ k V ( X ( ϕ ) , ϕ ) dϕ . (2.42)A closed magnetic field line ¯ X , such as the island centre, satisfies¯ X k + L = F Lk ( X k ) = ¯ X k . (2.43)For an initial condition that is infinitesimally close to the island centre, X ( ϕ k ) =¯ X k + δ X ( ϕ k ), the displacement from the closed field line as a function of toroidal anglesatisfies the linearized equation dδ X dϕ = (cid:2) ∇ X V (cid:0) ¯ X , ϕ (cid:1)(cid:3) (cid:124) · δ X , (2.44)where the Jacobian of the field-line-following equation is ∇ X V = ˆ e R B p ( ¯ X , ϕ ) B ϕ ( ¯ X , ϕ ) + R ∇ X B p ( ¯ X , ϕ ) B ϕ ( ¯ X , ϕ ) − R ( ∇ X B ϕ ( ¯ X , ϕ )) B p ( ¯ X , ϕ ) B ϕ ( ¯ X , ϕ ) . (2.45)We have denoted [ ∇ X V ] (cid:124) · δ ¯ X = δ ¯ R∂ V /∂R + δ ¯ Z∂ V /∂Z , and ˆ e R = ∇ X R . It is usefulto introduce a 2x2 matrix S k ( ϕ ) that solves a similar equation to (2.44), d S k dϕ = (cid:2) ∇ X V (cid:0) ¯ X , ϕ (cid:1)(cid:3) (cid:124) · S k ( ϕ ) , (2.46)with initial condition S k ( ϕ k ) = I . Then any solution to (2.44) can be found from δ X ( ϕ ) = S k ( ϕ ) δ X ( ϕ k ), as can be verified by substituting this expression into (2.44). Denoting S k ( ϕ k + q ) = S qk gives the equation δ X k + q = S qk δ X k , (2.47)so we can identify S qk with S qR,k in equations (2.30)-(2.31). The tangent map S qk is obtainedfrom the integral S qk ≡ I + (cid:90) ϕ k + q ϕ k (cid:2) ∇ X V (cid:0) ¯ X , ϕ (cid:1)(cid:3) (cid:124) · S k ( ϕ ) dϕ . (2.48)2 A. Geraldini, M. Landreman, E. Paul
In order to carry out this integration, the function ¯ X ( ϕ ) must be calculated separatelyfrom equation (2.41). As the tangent map is linear in δ X , it satisfies the property S qk ≡ S k + q − S k + q − . . . S k +1 S k . (2.49)The full-orbit tangent map is denoted M k = S Lk . (2.50)An important property of the full-orbit tangent map is that it has exactly unitdeterminant, det( M k ) = 1 for all k (Cary & Hanson 1986). This follows from theunderlying Hamiltonian nature of the magnetic field line trajectory (Meiss 1992). Thus,the characteristic equation for the eigenvalues of M k is λ − λ Tr( M k ) + 1 = 0, giving λ = 12 Tr( M k ) ± (cid:115)(cid:18)
12 Tr( M k ) (cid:19) −
1. (2.51)Hence, for | Tr( M k ) | < λ ± =exp ( ± iα ), with α = arccos (cid:18)
12 Tr( M k ) (cid:19) . (2.52)The angle α is the average angle of rotation around the island centre of a neighbouringtrajectory after the island centre comes back to its original poloidal position. For thisreason, it takes the same value irrespective of the fixed point k used to calculate it. Fromequation (2.19), the average angle of rotation after following around the island centreis 2 πωL/n . Hence, we obtain an expression for the average frequency of rotation oflinearized trajectories about the island centre, ω = n πL arccos (cid:18)
12 Tr( M k ) (cid:19) . (2.53)It is useful to recall that a closed magnetic field line is not necessarily an island centre.Since the relation det( M k ) = 1 always holds, the closed field line can be elliptic orhyperholic: it can be an O point or an X point respectively. When the closed field line isnot an O point, the concept of a rotation frequency breaks down. A quantity that canbe used to determine whether the closed field line is a centre or an X point is the residue(Greene 1968) of the full orbit tangent map, R = 12 −
14 Tr( M k ) . (2.54)For 0 < R <
1, a rotation frequency can be calculated from (2.53) and thus the closedfield line is a centre. If R < R > ω . In theisland width calculation of this section, it is not only assumed that 0 < R <
1, but alsothat
R (cid:28) πωL/n (cid:28) M k , satisfy the equation M T k σ M k = M k σ M T k = σ, (2.55)with σ = (cid:18) − (cid:19) . (2.56) djoint method for sensitivity of island size to magnetic field variations M T k ( σ M k ) M k = σ M k . Therefore, the matrix σ M k is aninvariant of the full-orbit tangent map. The symmetrized matrix W k = 12 ( σ M k + M (cid:124) k σ (cid:124) ) (2.57)satisfies the same property that the unsymmetrized counterpart σ M k satisfies, M T k · W k · M k = W k . (2.58)Consider the scalar invariant discussed after equation (2.27) calculated at ϕ = ϕ k , | δK | = δ ξ (cid:124) k · Q (cid:124) k · K · Q k · δ ξ k . This quantity is constant on a trajectory crossing the planes ϕ = 2 π ( k + LQ ) /n , where Q is an integer, as it directly follows from the quadraticperturbation to the Hamiltonian (in the magnetic coordinate analysis) about the islandcentre, | δK | = δ u (cid:124) · K · δ u . Similarly, the quantity δ X (cid:124) k · W k · δ X is constant on a trajectorycrossing the plane ϕ = 2 π ( k + LQ ) /n , since δ X (cid:124) k · M (cid:124) k · W k · M k · δ X k = δ X (cid:124) k · W k · δ X k .The vectors δ X k and δ ξ k are equivalent displacement vectors expressed in two differentcoordinate systems that are rotated with respect to one another. Since the scalar invariantremains unchanged after a rotation, the two invariants can only be related by an overallconstant, δ X (cid:124) k · W k · δ X k = γ k | δK | (Cary & Hanson 1991). Hence, the symmetric invariantmatrix W k has the same unit eigenvectors as the (approximately diagonal) matrix Q (cid:124) k · K · Q k = ι (cid:48) ( ψ ) (cid:16) ∂ψ∂ξ ⊥ (cid:17) + O (cid:0) (cid:15) ¯ C (cid:1) O (cid:16) M (cid:15) ¯ C (cid:17) O (cid:16) M (cid:15) ¯ C (cid:17) M (cid:15) (cid:16) ∂ξ (cid:107) ∂Θ (cid:17) + O (cid:16) M ˆ (cid:15)(cid:15) ¯ C (cid:17) , (2.59)where in the error terms we have assumed that | ∂ ξ /∂Θ | ∼ ¯ C and | ∂ ξ /∂ψ | ∼ ¯ Cι (cid:48) ( ψ ).Moreover, the eigenvalues of W k are equal to the eigenvalues of Q (cid:124) k · K · Q k multi-plied by a factor of γ k . From equation (2.59), the smallest eigenvalue of Q (cid:124) k · K · Q k is approximately equal to the bottom-right element, associated with the eigenvector( δξ ⊥ , δξ (cid:107) ) (cid:39) ( O (ˆ (cid:15) ) , O (ˆ (cid:15) )), and the largest eigenvalue is approximately equal to thetop-left element, associated with the eigenvector ( δξ ⊥ , δξ (cid:107) ) = (1 + O (ˆ (cid:15) ) , O (ˆ (cid:15) )). Therefore,the eigenvector of W k corresponding to the smallest eigenvalue is approximately tangentto the flux surface, and is denoted ˆ e k (cid:107) , and the eigenvector corresponding to the largesteigenvalue is approximately normal to the flux surface, and is denoted ˆ e k ⊥ . Since ˆ e k (cid:107) and ˆ e k ⊥ are just defined as eigenvectors of W k , their sign is yet unspecified. We imposethe constraint ˆ e (cid:107) ,k · M k · ˆ e ⊥ ,k > k to fix the relative sign of the eigenvectors.2.4. Island width
An expression for the island width at ϕ = ϕ k = 2 πk/n can be obtained by rearranging(2.39), ¯ w ⊥ ,k = 2 L ¯ CM π (cid:32) q + L − (cid:88) q = q δξ (cid:107) ,k + q δξ ⊥ ,k (cid:33) − . (2.60)Recall that δξ (cid:107) ,k = 0 and q is chosen for convenience according to equation (2.38), suchthat δξ ⊥ ,k + q (cid:39) q in the sum. An accurate calculation of the island width thusrests upon an accurate calculation of ¯ C and δξ (cid:107) ,k + q /δξ ⊥ ,k .To calculate δξ (cid:107) ,k + q /δξ ⊥ ,k , we follow the linearized magnetic field with initial condition δ X k · ˆ e (cid:107) k = 0, corresponding to δξ (cid:107) ,k (cid:39) δξ ⊥ ,k (cid:39) δ X k · ˆ e ⊥ ,k . Using equation (2.46),or repeated application of partial tangent maps as in (2.49), we obtain δ X k + q = S qk · δ X k .The final displacement from the fixed point in the direction parallel to the flux surface4 A. Geraldini, M. Landreman, E. Paul is δξ (cid:107) ,k + q (cid:39) ˆ e (cid:107) ,k + q · δ X k + q , giving δξ (cid:107) ,k + q δξ ⊥ ,k (cid:39) ˆ e (cid:107) ,k + q · S qk · ˆ e ⊥ ,k . (2.61)Equation (2.61) is defined to be positive for all k and q , which provides a second constraintto fix the sign of the vectors ˆ e ⊥ ,k and ˆ e (cid:107) ,k . With this additional constraint, the signs ofthe eigenvectors ˆ e ⊥ ,k , ˆ e (cid:107) ,k , ˆ e ⊥ ,k + q , ˆ e (cid:107) ,k + q for all k and q are specified with respect toone another. Note that the two constraints can only both be self-consistently applied ifthe small island width approximation is valid; if attempting to apply both constraintsleads to a contradiction, the islands are too wide. We then have q + L − (cid:88) q = q δξ (cid:107) ,k + q δξ ⊥ ,k (cid:39) Σ k ≡ q + L − (cid:88) q = q ˆ e (cid:107) ,k + q · S qk · ˆ e ⊥ ,k , (2.62)where we have defined the positive quantity Σ k .The circumference ¯ C can be approximated, for large enough number of fixed points,from the sum of the chords in the poloidal plane. However, to obtain the appropriatechords, the fixed points must first be reordered to proceed monotonically around theflux surface. To do so, we consider a reordering function ρ ( k ) such that p = ρ ( k ) is theindex labelling the reordered fixed points. We define ρ (0) = 0 such that p = 0 for k = 0,i.e. the set of reordered fixed points have the same point indexed as zero as the set offixed points ordered by following the field line. Moreover, we reorder the fixed pointsto monotonically circulate the magnetic axis in the same direction that they appear tocirculate when ordered following the magnetic field. The inverse of the reordering functionis denoted as ρ − and is defined such that ρ − ( ρ ( k )) = k , i.e. the field-line-following index k is returned for a given reordered index p . It can be shown that ρ ( k ) = ( n turns k ) mod L (2.63)and ρ − ( p ) = ( p mod n turns ) L + pn turns mod L, (2.64)where n turns is the total number of poloidal turns around the magnetic axis when followingthe field-line-ordered fixed points X k (i.e. the number of times that the magnetic axis appears to have been circulated poloidally). Here, k mod q indicates the remainder of k/q , equal to an integer between 0 and q −
1. The circumference ¯ C is approximated bythe sum of the chords obtained by joining the reordered fixed points,¯ C (cid:39) C ≡ L − (cid:88) k =0 | c k | , (2.65)where c k = X k − X k − , (2.66)and X k − is the fixed point preceding X k in the reordered set, such that k − = ρ − ( ρ ( k ) − . (2.67)The relative error in the chord approximation, i.e. ¯ C (cid:39) C , is O (1 /L ).Using equations (2.60), (2.62) and (2.65), the analytical small island width ¯ w ⊥ ,k can djoint method for sensitivity of island size to magnetic field variations w ⊥ ,k (cid:39) w ⊥ ,k , where w ⊥ ,k ≡ LCM πΣ k . (2.68)The value of M is obtained by counting how many of the L fixed points ¯ X k are also fixedpoints in the plane ϕ = 0, i.e. satisfy X (2 πL/n ) = X (0) = ¯ X k . The relative error in theisland width w ⊥ ,k compared to the true island width is O (ˆ (cid:15) / , L − ). The O (ˆ (cid:15) / ) piececomes from the error in approximating the width in the coordinate ψ using Υ in equation(2.23), and the additional error in relating changes in ψ to lengths in equation (2.32).The O ( L − ) piece comes from the chord approximation. Considering instead the relativeerror between w ⊥ ,k in equation (2.68) and ¯ w ⊥ ,k in equation (2.32), this is O (ˆ (cid:15), L − ).Here, the O (ˆ (cid:15) ) piece comes from the error in equation (2.61), which in turn comes fromthe O (ˆ (cid:15) ) error in the bottom-left matrix element of (2.31) due to the matrix K in (2.26)only being approximately diagonal.
3. Island size and residue variation using adjoint equations
Here we consider the variation of the different quantities — including the width —related to magnetic islands following the infinitesimally small variation of the magneticfield configuration from B ( R ) to B ( R ) + ∆ B ( R ). From equation (2.68), the variation ofthe island width is ∆w ⊥ ,k w ⊥ ,k (cid:39) ∆CC − ∆Σ k Σ k , (3.1)where ∆C is the variation of the circumference, approximated by the sum of chords inequation (2.65), and ∆Σ k is the variation of the sum in equation (2.62). Most of thissection is devoted to deriving expressions for ∆C and ∆Σ k in terms of ∆ B ( R ), using anadjoint method. In the final subsection, an expression for ∆ R is derived. Note that thevariation of the on-axis rotational transform is directly related to ∆ R if one considersthe magnetic axis as the periodic field line instead of an island chain O or X point.3.1. Circumference variation
A result of the change in magnetic configuration from B ( R ) to B ( R ) + ∆ B ( R ) is thatthe periodic field line position changes from ¯ X ( ϕ ) to ¯ X ( ϕ ) + ∆ ¯ X ( ϕ ). For the purpose ofthe island width calculation, the variation of the periodic field line position affects thecircumference. This changes from C in equation (2.65) to C + ∆C , where ∆C is givenby ∆C ( ¯ X ( ϕ ); ∆ ¯ X ( ϕ )) = L (cid:88) k =1 (cid:0) ∆ ¯ X k − ∆ ¯ X k − (cid:1) · ˆ c k , (3.2)where ˆ c k = c k / | c k | . Equation (3.2) can be re-expressed as ∆C ( ¯ X ( ϕ ); ∆ ¯ X ( ϕ )) = (cid:90) πL/n dϕ∆ ¯ X ( ϕ ) · (cid:34) L − (cid:88) k =0 ˆ c k (cid:0) δ ( ϕ − ϕ k ) − δ ( ϕ − ϕ k − ) (cid:1)(cid:35) . (3.3)Notice that the second sum in equation (3.3) can be re-cast as L − (cid:88) k (cid:48) =0 ˆ c k (cid:48) δ ( ϕ − ϕ k (cid:48)− ) = L − (cid:88) k =0 ˆ c k + δ ( ϕ − ϕ k ) , (3.4)6 A. Geraldini, M. Landreman, E. Paul where k + = ρ − ( ρ ( k ) + 1) , (3.5)since the order in which the sum is taken is not important. Re-expressing equation (3.3)using equation (3.4) gives the more convenient expression ∆C ( ¯ X ( ϕ ); ∆ ¯ X ( ϕ )) = (cid:90) πL/n dϕ∆ ¯ X ( ϕ ) · (cid:34) L − (cid:88) k =0 (cid:0) ˆ c k − ˆ c k + (cid:1) δ ( ϕ − ϕ k ) (cid:35) . (3.6)We impose the constraint that ¯ X ( ϕ ) + ∆ ¯ X ( ϕ ) remain a point along a closed magneticfield by introducing a Lagrangian L C = C + (cid:28) λ , d X dϕ − V ( X , ϕ ) (cid:29) , (3.7)where we have defined an inner product such that (cid:104) A , A (cid:105) = (cid:90) πL/n dϕ A · A . (3.8)Extremization of L with respect to λ leads to the constraint that X ( ϕ ) is the positionalong a magnetic field line, i.e. that X ( ϕ ) satisfies the differential equation (2.41). Wehave already found the periodic field line ¯ X ( ϕ ) satisfying (2.41) with periodic boundaryconditions. Now consider extrema of L with respect to changes in X ( ϕ ), ∆ L C = ∆C + (cid:28) λ , d∆ X dϕ − [ ∇ X V ] (cid:124) · ∆ X (cid:29) = 0 , (3.9)which can be re-expressed as ∆C − (cid:28) ∆ X , d λ dϕ + ∇ X V · λ (cid:29) + λ (2 πL/n ) · ∆ X (2 πL/n ) − λ (0) · ∆ X (0) = 0 . (3.10)By definition, the island centre maintains its periodicity after the perturbation, so weimpose ∆ X (2 πL/n ) = ∆ X (0). Additionally imposing the periodic boundary condition λ (0) = λ (2 πL/n ) one obtains the equation (cid:42) ∆ X , (cid:34) L − (cid:88) k =0 (cid:0) ˆ c k − ˆ c k + (cid:1) δ ( ϕ − ϕ k ) (cid:35) − d λ dϕ − ∇ X V · λ (cid:43) = 0 , (3.11)which leads to the differential equation d λ dϕ + ∇ X V · λ = L − (cid:88) k =0 (cid:0) ˆ c k − ˆ c k + (cid:1) δ ( ϕ − ϕ k ) . (3.12)Equation (3.12) is the adjoint equation useful for finding derivatives of the circumference.If ¯ X ( ϕ ) is a periodic solution of (2.41) and ¯ λ ( ϕ ) is a periodic solution of (3.12), then L C is stationary with respect to changes in ¯ X ( ϕ ) and ¯ λ ( ϕ ) arising from changes inthe magnetic field B ( R ). Therefore, L C is only affected by terms containing changes inthe magnetic field explicitly, ∆ L C = (cid:10) ¯ λ ( ϕ ) , − ∆ V ( ¯ X , ϕ ) (cid:11) . Furthermore, since ¯ X ( ϕ ) is amagnetic field line trajectory, the equation d ¯ X /dϕ = V ( ¯ X ( ϕ ) , ϕ ) is always satisfied and ∆C = ∆ L C (from equation (3.7)), giving ∆C = − (cid:90) πL/n dϕ ¯ λ · ∆ V ( ¯ X , ϕ ) , (3.13) djoint method for sensitivity of island size to magnetic field variations ∆ V ( ¯ X , ϕ ) = ¯ R∆ B p ( ¯ X , ϕ ) B ϕ ( ¯ X , ϕ ) − ¯ R B p ( ¯ X , ϕ ) ∆B ϕ ( ¯ X , ϕ ) B ϕ ( ¯ X , ϕ ) . (3.14)3.2. Variation of Σ k Upon varying the magnetic field from B ( R ) to B ( R ) + ∆ B ( R ), the quantity Σ k inequation (2.62) varies due to the variation of the eigenvectors ˆ e ⊥ k and ˆ e (cid:107) k + q and of thetangent map S qk , ∆Σ k = q + L − (cid:88) q = q (cid:0) ˆ e (cid:107) ,k + q · ∆ S qk · ˆ e ⊥ ,k + ∆e (cid:107) ,k + q ˆ e ⊥ k + q · S qk · ˆ e ⊥ ,k + ∆e ⊥ ,k ˆ e (cid:107) ,k + q · S qk · ˆ e (cid:107) ,k (cid:1) . (3.15)In equation (3.15), we have used the fact that the variation of a unit vector is perpen-dicular to the unit vector itself to re-express the variation of the normalized eigenvectorsof W k as ∆ ˆ e ⊥ ,k = ∆e ⊥ ,k ˆ e (cid:107) ,k and ∆ ˆ e (cid:107) ,k = ∆e (cid:107) ,k ˆ e ⊥ ,k . The slow rotation of nearby pointsaround the O point gives sin (2 πωk/n ) (cid:39) πωk/n ) (cid:39)
0. Thus, the matrixelements ˆ e ⊥ ,k + q · S qk · ˆ e ⊥ ,k and ˆ e (cid:107) ,k + q · S qk · ˆ e (cid:107) ,k , equal to the diagonal matrix elements in S qR,k (equation (2.31)), are both small in 2 πωL/n . This gives ∆Σ k (cid:39) q + L − (cid:88) q = q ˆ e (cid:107) ,k + q · ∆ S qk · ˆ e ⊥ ,k = q + L − (cid:88) q = q ˆ e (cid:107) ,k + q · ∆ s qk ; (3.16)here, we have introduced the variable s k ( ϕ ) = δ X ( ϕ ) / | δ X (2 πk/n ) | satisfying thedifferential equation (2.44), d s k dϕ = s k · ∇ X V , (3.17)with boundary condition s k (2 πk/n ) = s k = ˆ e ⊥ ,k , and defined the variation ∆ s k ( ϕ ). Were-express (3.16) to ∆Σ k (cid:39) Q − (cid:88) Q =0 (cid:90) π ( k + LQ + L ) /n π ( k + LQ ) /n ∆ s k ( ϕ ) · ˆ e (cid:107) ,k + q δ ( ϕ − ϕ k + q ) dϕ , (3.18)where the integer Q = (cid:98) q /L (cid:99) + 2, with the brackets (cid:98) and (cid:99) denoting the floor function,is chosen such that all values of q in the sum in (3.16) are counted. Introducing thenotation (cid:104) A · A (cid:105) k,Q = (cid:90) π ( k + LQ + L ) /n π ( k + LQ ) /n A · A dϕ, (3.19)we define the Lagrangian of Σ k , L Σ k = Σ k + Q − (cid:88) Q =0 (cid:32)(cid:28) λ k,Q , d X dϕ − V ( X , ϕ ) (cid:29) k,Q + (cid:28) µ k , d s k dϕ − s k · ∇ X V ( X , ϕ ) (cid:29) k,Q (cid:33) . (3.20)In (3.20) we have introduced two constraints using the adjoint variables λ k,Q ( ϕ ) and µ k ( ϕ ). With these constraints, the quantities L Σ k and Σ k are only equal to each other8 A. Geraldini, M. Landreman, E. Paul if X ( ϕ ) is a field line trajectory, satisfying equation (2.41), and s k satisfies the linearizedequation (3.17) for small displacements about the field line. † Extremization with respect to variations in X ( ϕ ) gives ∆ L Σ k = − Q − (cid:88) Q =0 (cid:28) ∆ X , d λ k,Q dϕ + ∇ X V · λ k,Q + s k · ∇ X ∇ X V · µ k (cid:29) + Q − (cid:88) Q =0 [ ∆ X (2 π ( k + QL + L ) /n ) · λ k,Q (2 π ( k + QL + L ) /n ) − ∆ X (2 π ( k + QL ) /n ) · λ k,Q (2 π ( k + QL ) /n )] = 0 , (3.21)where the Hessian of the field-line following equation (2.41) is ∇ X ∇ X V ( X , ϕ ) = 2ˆ e R ∇ X B p ( X , ϕ ) B ϕ ( X , ϕ ) − e R ( ∇ X B ϕ ( X , ϕ )) B p ( X , ϕ ) B ϕ ( X , ϕ ) + R ∇ X ∇ X B p ( X , ϕ ) B ϕ ( X , ϕ ) − R ( ∇ X B ϕ ( X , ϕ ))( ∇ X B p ( X , ϕ )) B ϕ ( X , ϕ ) − R ( ∇ X ∇ X B ϕ ( X p , ϕ )) B p ( X , ϕ ) B ϕ ( X , ϕ ) + 2 R ( ∇ X B ϕ ( X , ϕ ))( ∇ X B ϕ ( X , ϕ )) B p ( X , ϕ ) B ϕ ( X , ϕ ) . (3.22)We consider X ( ϕ ) to be a periodic field line ¯ X ( ϕ ) (an island centre) after the variationof the magnetic field configuration, such that ∆ X (2 π ( k + QL ) /n ) = ∆ X (2 π ( k + QL + L ) /n ). Therefore, all the boundary terms in equation (3.21) vanish if λ k,Q (2 π ( k + QL ) /n ) = λ k,Q (2 π ( k + QL + L ) /n ). The adjoint equations for λ k,Q ( ϕ ) are thus d λ k,Q dϕ + ∇ X V · λ k,Q + s k · ∇ X ∇ X V · µ k = 0 , (3.23)to be solved for periodic solutions ¯ λ k,Q ( ϕ ) in the interval 2 πQL/n (cid:54) ϕ − πk/n (cid:54) π ( QL + L ) /n for all k and Q .Extremizing L Σ k with respect to variations in s k ( ϕ ) gives an equation for µ k , ∆ L Σ k = Q − (cid:88) Q =0 (cid:42) ∆ s k , q + L − (cid:88) q = q ˆ e (cid:107) k + q δ ( ϕ − ϕ k + q ) − d µ k dϕ − ∇ X V · µ k (cid:43) Q + µ k (2 π ( k + LQ ) /n ) · ∆ s k (2 π ( k + LQ ) /n ) − µ k (2 πk/n ) · ∆ s k (2 πk/n ) = 0 , (3.24)where all boundary terms at ϕ = 2 π ( k + LQ ) /n for Q (cid:54) = 0 and Q (cid:54) = Q have vanishedbecause µ k is a continuous variable in the interval 2 πk/n (cid:54) ϕ (cid:54) π ( k + Q L ) /n . Whenconsidering variations in s k , the initial condition s k (2 πk/n ) = s k = ˆ e ⊥ k can be assumedto be unchanged and therefore ∆ s k (2 πk/n ) = 0. Hence, the boundary terms in (3.24)vanish by imposing µ k (2 π ( k + LQ ) /n ) = 0, resulting in the differential equation d µ k dϕ = q + L − (cid:88) q = q ˆ e (cid:107) k + q δ ( ϕ − ϕ k + q ) − ∇ X V · µ k . (3.25)Note that µ k (2 πk/n ) can be obtained from µ k (2 π ( k + LQ ) /n ) = 0 using equation(3.25). This can be solved by repeated applications of appropriate partial tangent maps † Note that the calculation could be carried out without replacing S k by s k if the vectoradjoint variable µ k were replaced by a matrix µ k . djoint method for sensitivity of island size to magnetic field variations ϕ k + q ; as shown in appendix C, the linear mappingfrom µ k to µ k +1 that follows from the homogeneous term on the right hand side of (3.25)is the adjoint of the partial tangent map S k , such that µ k (2 π ( k + q ) /L ) = (cid:16) S k + q (cid:17) (cid:124) · µ k (2 π ( k + q + 1) /L ). The adjoint variables λ k,Q ( ϕ ) can be obtained from equation (3.23)once the the adjoint variables µ k ( ϕ ) are obtained from (3.25).The magnetic field configuration is varied while considering ¯ X ( ϕ ) to be the magneticfield line trajectory at the island centre, satisfying equation (2.41) with a periodicboundary condition, and while constraining s k ( ϕ ) to satisfy equation (3.17) for linearizedtrajectories about the island centre. We thus conclude that Σ k = L Σ k and thus ∆Σ k = ∆ L Σ k . Therefore, the variation of Σ k is given by ∆Σ k = − Q − (cid:88) Q =0 (cid:90) π ( k + QL ) /n πk/n dϕ (cid:0) ¯ λ k,Q · ∆ V ( ¯ X , ϕ ) + s k · ∇ X ∆ V ( ¯ X , ϕ ) · µ k (cid:1) , (3.26)where ¯ λ k,Q is a periodic solution of (3.23) and µ k is the solution of (3.25) satisfying theboundary condition µ k (2 π ( k + LQ ) /n ) = 0. The function ∇ X ∆ V ( ¯ X , ϕ ) can be writtenexplicitly in terms of variations of the magnetic field and its gradients by differentiating(3.14), ∇ X ∆ V ( ¯ X , ϕ ) = ˆ e R ∆ B p ( ¯ X , ϕ ) B ϕ ( ¯ X , ϕ ) − ˆ e R B p ( ¯ X , ϕ ) ∆B ϕ ( ¯ X , ϕ ) B ϕ ( ¯ X , ϕ ) + R ∇ X ∆ B p ( ¯ X , ϕ ) B ϕ ( ¯ X , ϕ ) − R ( ∇ X B p ( ¯ X , ϕ )) ∆B ϕ ( ¯ X , ϕ ) B ϕ ( ¯ X , ϕ ) − R ∇ X B ϕ ( ¯ X , ϕ ) ∆ B p ( ¯ X , ϕ ) B ϕ ( ¯ X , ϕ ) − R ∇ X ( ∆B ϕ ( ¯ X , ϕ )) B p ( ¯ X , ϕ ) B ϕ ( ¯ X , ϕ ) + 2 R ( ∇ X B ϕ ) B p ∆B ϕ B ϕ . (3.27)3.3. Residue variation
The residue R varies due to the variation of the full orbit tangent map M k . We re-express the trace of the full orbit tangent map asTr( M ) = I : M , (3.28)where the operation denoted by the colon : corresponds to multiplying each matrixelement of one matrix with the corresponding matrix element of the other matrix andadding the results together. Note that k = 0 was chosen in (3.28), as the residue isindependent of k . From equation (2.54) and (3.28) we obtain ∆ R = − I : ∆ M . (3.29)Recall that M = S (2 πL/n ) and that the tangent map S satisfies the differentialequation (2.46) with boundary condition S (0) = I . We define the Lagrangian of R , L R = R + (cid:90) πL/n dϕ (cid:20) λ R · (cid:18) d X dϕ − V ( X , ϕ ) (cid:19) + µ R : (cid:18) d S dϕ − ( ∇ X V ( X , ϕ )) (cid:124) · S (cid:19)(cid:21) . (3.30)In (3.30) we have introduced two constraints using the adjoint variables λ R ( ϕ ) (a vector)and µ R ( ϕ ) (a matrix). With these constraints, the quantities L R and R are equal to eachother if X ( ϕ ) is a field line trajectory, satisfying equation (2.41), and S satisfies equation(2.46) for the tangent map.0 A. Geraldini, M. Landreman, E. Paul
Extremization with respect to variations in X ( ϕ ) gives ∆ L R = − (cid:28) ∆ X , d λ R dϕ + ∇ X V · λ R + ( ∇ X ∇ X V · µ R ) : S (cid:29) + ∆ X (2 πL/n ) · λ R (2 πL/n ) − ∆ X (0) · λ R (0) = 0 , (3.31)Again, all the boundary terms in equation (3.31) vanish if λ R (2 πL/n ) = λ R (0). Theadjoint equation for λ R ( ϕ ) is thus d λ R dϕ + ∇ X V · λ R + ( ∇ X ∇ X V · µ R ) : S = 0 , (3.32)to be solved in the interval 0 (cid:54) ϕ (cid:54) πL/n for periodic solutions, denoted ¯ λ R ( ϕ ).Extremizing L R with respect to variations in S ( ϕ ) gives an equation for µ R , ∆ L R = − I : ∆ M + (cid:28) µ R , d∆ S dϕ − ( ∇ X V ) (cid:124) · ∆ S (cid:29) = 0 . (3.33)Note that we have used the following definition for the inner product of two matrixquantities, (cid:104) A , A (cid:105) = (cid:90) πL/n dϕ A : A . (3.34)We re-express (3.33) to − I : ∆ M − (cid:28) ∆ S , dµ R dϕ + ∇ X V · µ R (cid:29) + µ R (2 πL/n ) : ∆S (2 πL/n ) − µ R (0) : ∆S (0) = 0 . (3.35)Noting that ∆ S (0) = 0 and ∆S (2 πL/n ) = ∆ M both hold true by definition, theboundary terms are zero provided µ R (2 πL/n ) = I is chosen. Thus, the adjoint equationfor µ R is dµ R dϕ = −∇ X V · µ R , (3.36)with µ R (2 πL/n ) = I as a boundary condition. From appendix C, an equivalent formof this boundary condition is µ R (0) = M (cid:124) .As before, the variation of the residue is given by ∆ R = − (cid:90) πL/n (cid:0) ¯ λ R · ∆ V ( ¯ X , ϕ ) + ( ∇ X ∆ V ( ¯ X , ϕ ) · µ R ) : S (cid:1) dϕ, (3.37)where ¯ λ R and µ R are the solutions of the adjoint equations (3.23) and (3.25) with theappropriate boundary conditions.To conclude this section, we briefly discuss a subsidiary result of the preceedinganalysis. The trace of the full orbit tangent map calculated at the magnetic axis isrelated the rotation angle of nearby trajectories about the magnetic axis. By definition,this is proportional to the on-axis rotational transform. Equating (2.52) (with k = 0) tothe product of the interval in toroidal angle, 2 π/n , and the on-axis rotational transform,¯ ι , and re-arranging for ¯ ι gives ¯ ι = n π arccos (cid:18)
12 Tr ( M ) (cid:19) . (3.38)Note that this equation is correct provided ¯ ι < n , a condition which is satisfied in most djoint method for sensitivity of island size to magnetic field variations ∆ ¯ ι = n π sin (2 π ¯ ι/n ) ∆ R . (3.39)
4. Numerical results
In this section, we present the numerical results obtained for the gradients of theisland width and of other properties of the periodic field line. In section 4.1, we brieflyexplain the numerical scheme used to obtain our results. We then present, in section 4.2numerical results for the island width, its gradient and the gradient of other island-relatedquantities in an analytical magnetic configuration studied in Reiman & Greenside (1986).In section 4.3, we present results for the shape gradient of the width of a magnetic islandin NCSX with respect to the positions on a type A modular coil. Finally, in section 4.4,we apply the gradient of the residue of a periodic field line to optimize a helical magneticconfiguration of the kind studied in Cary & Hanson (1986) and Hanson & Cary (1984).4.1.
Numerical scheme
A Runge-Kutta 4th order explicit scheme is used to integrate equations (2.41), (2.46),(3.12), (3.23) and (3.25) in toroidal angle ϕ . The number of grid points in ϕ per fieldperiod is denoted N ϕ , such that the number of grid points per toroidal turn is n N ϕ .A Newton method is used to search for periodic solutions of the magnetic field linesuch as the magnetic axis and a magnetic island centre. The search proceeds as follows.The position of the periodic field line is initially guessed as X (0) and equation (2.41)is integrated in ϕ from ϕ = 0 to ϕ = 2 πL/n , where L is an integer (= 1 for themagnetic axis). The tangent map is also integrated following the magnetic field line.The final position X (2 πL/n ) and tangent map S (2 πL/n ) are used to evaluate a nextguess for the magnetic axis as follows. The step X step (0) in position necessary to movecloser to the periodic solution at the next iteration, X i +1 (0) = X i (0) + X i step (0), iscalculated by imposing X i (0)+ X i step (0) = X i (2 πL/n )+ X i step (2 πL/n ) on the linearizedequations. This leads to X i (0) + X i step (0) = X i (2 πL/n ) + S i (2 πL/n ) · X i step (0) and,upon rearranging, to X i step (0) = (cid:0) I − S i (2 πL/n ) (cid:1) − · (cid:0) X i (2 πL/n ) − X i (0) (cid:1) . (4.1)The error is calculated from E i = (cid:12)(cid:12) X i (2 πL/n ) − X i (0) (cid:12)(cid:12) | X i (0) | (4.2)where the magnitude of the poloidal vector X = ( R, Z ) is defined as | X | = √ R + Z .A periodic field line solution is found if E i falls below a threshold value E thresh = 10 − ;we then consider ¯ X ( ϕ ) = X i ( ϕ ).Once a periodic field line is found, the values of the magnetic field and its first andsecond derivatives on the toroidal grid points and in the intermediate Runge Kutta stepsare stored to accelerate subsequent parts of the code such as the solutions of the adjointequations and the calculations of the gradients.2 A. Geraldini, M. Landreman, E. Paul
RZ R
Figure 1.
Poincar´e plots showing magnetic field lines near the island separatrix in the Reimanmagnetic field configuration with ι ax = 0 . ι (cid:48) ax = 0 .
38 and ε i = 0 for i (cid:54) = 6, for ε = 0 . ε = 0 .
01 (right).
Island width and gradient calculation: Reiman model
We study magnetic configurations of a form similar to the model field used in section 5of Reiman & Greenside (1986), given by equation (2.1) with ψ = 12 r , (4.3) χ = ι ax ψ + ι (cid:48) ax ψ − k max (cid:88) k =1 ε k (2 ψ ) k/ cos ( kθ − ϕ ) , (4.4)where k max is the largest value of k for which ε k (cid:54) = 0. Here, r and θ are chosen to be thepoloidal minor radius and the geometric poloidal angle, such that r = (cid:113) ( R − + Z , (4.5)and tan θ = ZR − . (4.6)The explicit expressions for the magnetic field and its derivatives are derived from equa-tions (2.1) and (4.3)-(4.6), and are given in Appendix D. The unperturbed configuration( ε k = 0 for all k ) is symmetric in the geometric poloidal and toroidal angles (it is notcurl-free and is therefore not a vacuum magnetic configuration). For the scope of thispaper, we focus on the set of parameters ι ax = 0 . ι (cid:48) ax = 0 .
38 and ε k = 0 for k (cid:54) = 6. Infigure 1, we show a Poincar´e plot of the island chain in the plane ϕ = 0 resulting fromthe parameters ε = 0 .
01 and ε = 0 . ε k = 0 for k (cid:54) = 6 in equation (4.4) results in equation (2.2) with χ = ι ax ψ + ι (cid:48) ax ψ and χ = ε (2 ψ ) cos ( kθ − ϕ + π ). The (unperturbed) minor radius ofthe island chain is obtained by calculating r = ¯ r at the resonance, ι res = 1 / ι ax + ι (cid:48) ax ¯ r ,giving ¯ r = (cid:112) ( ι res − ι ax ) /ι (cid:48) ax and ψ = ( ι res − ι ax ) / (2 ι (cid:48) ax ). Therefore, the circumference of djoint method for sensitivity of island size to magnetic field variations -10 -9 -8 -7 -6 -5 -4 -3 -2 -6 -5 -4 -3 -2 -1 w i d t h (a) w ⟂ w ⟂ ¯C/C¯w ⟂ -10 -9 -8 -7 -6 -5 -4 -3 -2 ε -10 -9 -8 -7 -6 -5 -4 -3 -2 -1 e rr o r (b) Figure 2.
Width of magnetic islands at the resonant flux surface with rotational transform ι res = 1 / ι ax = 0 . ι (cid:48) ax = 0 .
38 and ε i = 0 for i (cid:54) = 6, shown as a function of ε . (a) Uncorrected w ⊥ ( × ) and corrected w ⊥ ¯ C/C (+)computed values of width are compared with the analytical value ¯ w ⊥ in equation (4.9) (solidline). Here, w ⊥ is calculated from equation (2.68), C from equation (2.65) and ¯ C from equation(4.7). (b) Normalized error | − w ⊥ / ¯ w ⊥ | ( × ) and | − w ⊥ ¯ C/ ( ¯ w ⊥ C ) | (+). For ε > − , the errorin the corrected width decreases linearly with ε , as expected from the discussion at the end ofsection 2. For smaller values of ε this error changes sign and increases with ε − , most likelydue to roundoff error propagation. One in five markers are shown in both plots. The toroidalangle resolution was N ϕ = 80. the unperturbed resonant flux surface section in the Reiman model is¯ C = 2 π (cid:114) ι res − ι ax ι (cid:48) ax . (4.7)The width of the island chain expressed in the magnetic coordinate ψ , obtained byinserting (cid:15) ( ψ ) = ε (( ι res − ι ax ) /ι (cid:48) ax ) and ι (cid:48) ( ψ ) = 2 ι (cid:48) ax into equation (2.23), is Υ = 4 (cid:114) ε ι (cid:48) ax (cid:18) ι res − ι ax ι (cid:48) ax (cid:19) / . (4.8)From equation (2.32) and the equality ∂ξ ⊥ /∂ψ = dr/dψ = 1 / ¯ r , the expression for theisland width in terms of the parameters of the Reiman model is¯ w ⊥ = 4 (cid:114) ε ι (cid:48) ax ι res − ι ax ι (cid:48) ax . (4.9)Note that the subscript k in w ⊥ ,k is unnecessary here, since the poloidal symmetry ofthe Reiman model implies that the width of all islands in the chain is identical.4 A. Geraldini, M. Landreman, E. Paul
In figure 2 we plot the island width calculated using equation (2.68) and compare itwith that calculated using equation (4.9). The discrepancy between the two equations isalmost entirely due to the chord approximation for the circumference in equation (2.65). † This can be seen by comparing the island width calculated using the two methods with thequantity ¯ Cw ⊥ /C , a corrected island width where the sum of chords in (2.65) is replacedwith the more accurate measure of circumference in (4.7). The corrected island widthhas a near-perfect overlap with the width calculated from (4.9). The error between thetwo quantities decreases with ε — consistent with the discussion in the final paragraphof section 2 (recall that (4.9) comes from (2.32))— and is thus limited by the accuracy ofthe small-island analysis. Conversely, the uncorrected island width w ⊥ is a less accurateapproximation, as seen by the saturation of the error with decreasing ε . This is due tothe chord approximation of the circumference limiting the accuracy of the island widthevaluation, as the O ( L − ) error is dominant at small enough ε ∝ ˆ (cid:15) .The gradient of the island width with respect to the parameters κ ∈ { ι ax , ι (cid:48) ax , ε } is calculated using the method derived in this paper. When one such parameter isvaried infinitesimally, the infinitesimal magnetic field variation can be expressed as ∆ B = ∆κ∂ B /∂κ , where ∂ B /∂κ is the gradient of the magnetic field with respectto κ . The infinitesimal variation of the magnetic field gradient can be expressed as ∆ ∇ X B = ∆κ∂ ( ∇ X B ) /∂κ . Both ∂ B /∂κ and ∂ ( ∇ X B ) /∂κ are straightforwardly obtainedfrom the equations in appendix D. The variation of the circumference and of the tangentmap can also be expressed as ∆C = ∆κ∂C/∂κ and ∆Σ k = ∆κ∂Σ k /∂κ . Thus, thegradient of the circumference with respect to κ is ∂C∂κ = − (cid:90) πL/n dϕ ∂ V ∂κ · λ , (4.10)and the gradient of Σ k is ∂Σ k ∂κ = − Q − (cid:88) Q =0 (cid:90) π ( k + LQ + L ) /n π ( k + LQ ) /n (cid:18) ∂ V ∂κ · λ Q + s ⊥ k · (cid:18) ∂ ∇ X V ∂κ (cid:19) · µ (cid:19) , (4.11)where ∂ V ∂κ = RB ϕ ∂ B p ∂κ − R B p B ϕ ∂B ϕ ∂κ , (4.12) ∇ X ∂ V ∂κ = ˆ e R B ϕ ∂ B ∂κ − ˆ e R B p B ϕ ∂B ϕ ∂κ + R ∇ X B ϕ ∂ B p ∂κ − R ( ∇ X B p ) B ϕ ∂B ϕ ∂κ − R ( ∇ X B ϕ ) B ϕ ∂ B ∂κ − (cid:18) ∇ X ∂B ϕ ∂κ (cid:19) R B p B ϕ + 2 R ( ∇ X B ϕ ) B p B ϕ ∂B ϕ ∂κ . (4.13)Hence, and using equation (3.1), the gradient of the island width is given by ∂ ln w k ∂κ = ∂ ln C∂κ − ∂ ln Σ k ∂κ , (4.14)with ∂C/∂κ and ∂Σ k /∂κ given by equations (4.10) and (4.11) respectively. The gradientof the residue of a periodic field line follows from (3.37), ∂ R ∂κ = − (cid:90) πL/n (cid:18) λ R · ∂ V ∂κ + (cid:18) ∂ ∇ X V ∂κ · µ (cid:19) : S (cid:19) dϕ. (4.15) † Note that the other approximation (2.37) is exact in the Reiman model due to theunperturbed configuration being poloidally symmetric. djoint method for sensitivity of island size to magnetic field variations -8 -6 -4 -2 E R ( δ ) ι ax ι ′ax ε δ -10 -8 -6 -4 -2 E C ( δ ) -8 -6 -4 -2 E Σ ( δ ) -4 -3 -2 -1 finite difference step size -8 -6 -4 -2 E w ⟂ ( δ ) Figure 3.
Errors, relative to a centred difference approximation, in the gradients of residue R , circumference C , Σ and island width w ⊥ calculated with respect to the on-axis rotationaltransform ι ax , its first derivative ι (cid:48) ax and the amplitude of the resonant perturbation ε inthe Reiman model. The errors E ( δ ), defined in equation (4.18), are shown as a functionof a normalized finite difference step size R − ( ∂ R /∂κ ) δ . The configuration parameters are ι ax = 0 . ι (cid:48) ax = 0 . ε = 0 .
01 and ε i = 0 for i (cid:54) = 6. The resonant flux surface has rotationaltransform ι res = 1 /
6. The dashed line is E ( δ ) = δ . The toroidal angle resolution was N ϕ = 80. Considering the residue at the magnetic axis, the gradient of the on-axis rotationaltransform follows from (3.39), ∂ ¯ ι∂κ = n π sin (2 π ¯ ι/n ) ∂ R ∂κ . (4.16)where ¯ ι is obtained from (3.38) and is found to be equal to ι ax , as expected.The results of equations (4.10), (4.11), (4.14) and (4.15) can be compared witha numerical derivative of the quantities C , Σ k , w ⊥ ,k and R calculated using finitedifferences. To this end, we calculate centred difference (CD) derivatives of the form ∂ CD ln f∂κ ( κ, δ ) = ln [ f ( κ + δ )] − ln [ f ( κ − δ )]2 δ . (4.17)The error with respect to the adjoint calculation of the derivative in equation (4.17) is E f ( δ ) = (cid:12)(cid:12)(cid:12)(cid:12) ∂ CD ln f∂κ ( κ, δ ) − ∂ ln f∂κ ( κ ) (cid:12)(cid:12)(cid:12)(cid:12) (4.18)In equations (4.17)-(4.18), f ∈ {R , C, Σ k , w ⊥ ,k } . We expect the centred difference deriva-tive of the numerically calculated island width to approach the numerically calculated6 A. Geraldini, M. Landreman, E. Paul derivative as δ is decreased. Having ignored the pieces of the derivative of Σ k that aresmall in ˆ (cid:15) in equation (3.15), one might expect that E ( δ ) ∼ δ for δ larger than athreshold value and E ( δ ) ∼ ˆ (cid:15) for smaller δ . However, comparing the Reiman model tothe derivation of the island width, the quantity ζ ( ψ ) is independent of ψ in the Reimanmodel. Thus, the diagonal elements of the tangent map in (2.17), and correspondinglythe off-diagonal elements of the scalar invariant matrix in (2.26), are exactly zero. Hence,the eigenvectors ˆ e ⊥ and ˆ e (cid:107) are exactly aligned with ∇ X ψ and ∇ X θ respectively, for anyvalue of the parameters ι ax , ι (cid:48) ax and ε . Thus, the terms small in ˆ (cid:15) that were neglectedin the derivation of the gradient of the island width are exactly zero in the Reimanmodel, and E ( δ ) ∼ δ should hold for all values of δ . In figure 3 the quantity E f ( δ ) for f ∈ {R , C, Σ, w ⊥ } is shown for derivatives with respect to κ ∈ ( ι ax , ι (cid:48) ax , ε ) for ι ax = 0 . ι (cid:48) ax = 0 .
38 and ε = 0 .
01. The proportionality E f ( δ ) ∝ δ for all κ is strong evidencethat the gradient calculation is accurate for all quantities. In addition, the gradient ofthe on-axis rotational transform ¯ ι calculated from equation (4.16) is, as expected, unityfor κ = ι ax and zero for κ ∈ { ι (cid:48) ax , ε } .4.3. Shape gradient calculation: explicit coils
In this section, it will be useful to denote as R = ( X, Y, Z ) the position along a magneticfield line expressed in a set of right-handed Cartesian coordinates. The coordinate Z is thesame as in the two-dimensional vector X = ( R, Z ), while X = R cos ϕ and Y = R sin ϕ .The magnetic field produced by explicit coils is calculated using the Biot-Savart law.The number of inputs required for an explicit coil calculation is just the number offield periods and a set of (sufficiently resolved) positions along the coils. Since thecoils producing the magnetic field are continuous (even though they are numericallyapproximated as discrete), the magnetic field is a functional of the continuous periodicfunction r c ( l c ) = ( x c , y c , z c ) specifying the coil shape. Lowercase letters and a c subscriptdistinguish coil positions from positions along a magnetic field line. Here l c is the arclength along the coil measured from some reference point, l c ∈ [0 , L c ), where L c is thetotal coil length. The magnetic field is specified by the Biot-Savart law B coil = − N (cid:88) c =1 µ I c π (cid:90) dl c R − r c | R − r c | × d r c dl c . (4.19)The poloidal gradient of the magnetic field is ∇ X B coil = N (cid:88) c =1 µ I c π (cid:90) dl c | R − r c | (cid:20) − I p + 3 ( X − x c )( R − r c ) | R − r c | (cid:21) × d r c dl c . (4.20)where we introduced the poloidal identity matrix I p = ˆ e R ˆ e R +ˆ e Z ˆ e Z , with ˆ e Z = ∇ X Z , andthe position vector of the coil projected in the poloidal plane, x c = r c · (ˆ e R ˆ e R + ˆ e Z ˆ e Z ).The second derivative of the magnetic field is ∇ X ∇ X B coil = N (cid:88) c =1 µ I c π (cid:90) dl c | R − r c | (cid:20) (cid:18) I p − X − x c )( X − x c ) | R − r c | (cid:19) ( R − r c )+3 ( X − x c ) I p + 3 (ˆ e R ( X − x c ) ˆ e R + ˆ e Z ( X − x c ) ˆ e Z )] × d r c dl c . (4.21)The variation ∆ B of a coil-produced magnetic field is conveniently expressed in termsof the shape gradient G r c B of the magnetic field with respect to coils (Landreman &Paul 2018). The shape gradient G r c is a vector operator which, like ∇ X , is meaningful djoint method for sensitivity of island size to magnetic field variations R −0.6−0.4−0.20.00.20.40.6 Z Figure 4.
Poincar´e plot for the configuration produced by NCSX coils. The island for whichthe shape gradient of the width is calculated in figure 5 is shown enlarged in the inset; its widthusing (2.68) is w ⊥ = 0 . R = 0 . only when acting on a (scalar or tensor) quantity to its right. It is defined via ∆f = N (cid:88) c =1 (cid:90) dl c (cid:18) ∆ r c × d r c dl c (cid:19) · G r c f . (4.22)The shape gradient of the magnetic field can be extracted from the Biot-Savart law and(4.22) as shown in appendix E, giving G r c B coil = µ I c π | R − r c | (cid:20) I − R − r c ) ( R − r c ) | R − r c | (cid:21) , (4.23)where I = ˆ e X ˆ e X + ˆ e Y ˆ e Y + ˆ e Z ˆ e Z .The shape gradient of the island width with respect to the coils producing the magneticfield is G r c ln w k = G r c ln C − G r c ln Σ k . (4.24)Here G r c C is given by G r c C = − (cid:90) πL/n dϕ G r c V · λ , (4.25)with G r c V = R G r c B p B ϕ − R ( G r c B ϕ ) B p B ϕ , (4.26)8 A. Geraldini, M. Landreman, E. Paul ¯ G x c l n w ⟂ adjointcentred difference0 1 2 3 4 5 6 7−2000−100001000200030004000 ¯ G y c l n w ⟂ distance around coil l c −400−20002004006008001000 ¯ G z c l n w ⟂ Figure 5.
Shape gradient of the width of one of the magnetic islands in the NCSX configuration(figure 4) with respect to the positions along a type A modular coil (length L c = 7 . N ϕ = 30 with: the adjoint method (solid lines); the centred difference scheme with δκ c = 10 − (dashed lines). For each component, the mean residual between the two calculationsis about 2% of the mean absolute value. The adjoint calculation is over a hundred times faster. and G r c Σ k is given by G r c Σ k = − Q − (cid:88) Q =0 (cid:90) πL ( Q +1) /n πLQ/n dϕ (( G r c V ) · λ Q + s ⊥ k · ( ∇ X G r c V ) · µ ) , (4.27)with ∇ X G r c V = ˆ e R G r c B B ϕ − ˆ e R B B ϕ G r c B ϕ + R ∇ X G r c B B ϕ − R ( ∇ X B ) G r c B ϕ B ϕ − R ( ∇ X B ϕ ) G r c B B ϕ − R ( ∇ X G r c B ϕ ) B B ϕ + 2 R ( ∇ X B ϕ ) B G r c B ϕ B ϕ . (4.28)The spatial derivative (in the poloidal plane) of the shape gradient of the magnetic fieldis ∇ X G r c B = µ I c π | R − r c | (cid:20) X − x c ) (cid:18) I − R − r c ) ( R − r c ) | R − r c | (cid:19) +3ˆ e R ( R − r c ) ˆ e R + 3ˆ e R ˆ e R ( R − r c ) + 3ˆ e Z ( R − r c ) ˆ e Z + 3ˆ e Z ˆ e Z ( R − r c )] . (4.29)Using these equations, a numerical approximation to the shape gradient of the widthof an island in an NCSX vacuum configuration, shown in figure 4, was calculated for a djoint method for sensitivity of island size to magnetic field variations ι res = 1 /
4. The shape gradient calculation is compared in figure 5 with a finite differenceapproximation calculated by perturbing the discrete coil positions at certain locations.In figures 4 and 5, all plotted quantities are in SI units. For simplicity, in figure 5 we haveplotted an alternative definition of the shape gradient ¯ G c f = ( d r c /dl c ) × G r c f , whichsatisfies ∆f = N (cid:88) c =1 (cid:90) dl c ∆ r c · ¯ G r c f (4.30)and thus has a more intuitive geometric interpretation. In the adjoint calculation, theshape gradient is calculated by evaluating G r c f using the adjoint expressions outlined inthis section and a numerical approximation of d r c /dl c ,¯ G r c f = r c ( l c + δl c ) − r c ( l c − δl c )2 δl c × G r c f , (4.31)where δl c is the separation between equally spaced positions on the coil. Using a centreddifference approximation instead, the components of the shape gradient, ¯ G CD κ c f for κ c ∈{ x c , y c , z c } , are calculated from¯ G CD κ c f = f ( κ c + δκ c ) − f ( κ − δκ c )2 δκ c δl c . (4.32)For the coil considered in this study, the centred difference shape gradient taken with δκ c = 10 − has a good overlap with the shape gradient calculated using the adjointmethod, as shown in figure 5. The direct adjoint calculation is, however, over a hundredtimes faster. This is because the adjoint approach requires solving a few equations alongthe same periodic magnetic field line, so that the magnetic field and its gradients donot need to be re-evaluated when computing the shape gradient with respect to allcoil positions. Conversely, the finite difference calculation requires searching for a newperiodic field line after each finite difference step, which is expensive because it requiresre-computing the magnetic field at new locations for each coil perturbation.4.4. Optimization of helical coils via gradient of residue
We proceed to consider a model magnetic field studied by Hanson & Cary (1984) whichconsists of two pieces: a toroidal magnetic field generated by a long current-carrying wirepassing through the centre of the torus, and a magnetic field generated by a pair ofhelical coils with opposite current ( ± I heli ). The helical coil positions are specified by thepoloidal angle η ± as a function of the toroidal angle ϕ , η ± = n ϕl + k max (cid:88) k =0 (cid:20) A ± ,k cos (cid:18) kn ϕl (cid:19) + B ± ,k sin (cid:18) kn ϕl (cid:19)(cid:21) . (4.33)The position of the coils r ± = ( x ± , y ± , z ± ) is then obtained using the equations x ± = ( R + r cos η ± ) cos ϕ, (4.34) y ± = ( R + r cos η ± ) sin ϕ, (4.35) z ± = − r sin η ± . (4.36)0 A. Geraldini, M. Landreman, E. Paul
Here we have introduced the major and minor radius of the toroidal surface in which thehelical coils lie, R and r respectively. The magnetic field configuration is obtained byadding the toroidal magnetic field ˆ e ϕ R B t /R to the magnetic field obtained by applyingthe Biot-Savart law to the helical coils, B heli ( R ) = R B t R ˆ e ϕ + (cid:88) ± ± µ I heli π (cid:90) dϕ | R − r ± | d r ± dϕ × ( R − r ± ) . (4.37)Upon fixing R , r , I heli and B t , as done in Hanson & Cary (1984), the continuousparameters that can be used to perturb the magnetic field configuration are A k for0 (cid:54) k (cid:54) k max and B k for 1 (cid:54) k (cid:54) k max . Perturbing any one parameter of any onecoil (carrying a current ± I heli ), denoted κ ± , such that the variation in the coil positionis r ± ( ϕ ) → r ± ( ϕ ) + ∆κ ± ∂ r ± ( ϕ ) /∂κ ± causes the magnetic field to change such that B ( R ) → B ( R ) + ∆κ ± ∂ B ( R ) /∂κ ± . From the expressions (4.22) and (4.23), we get ∂ B ∂κ ± = ± µ I heli π (cid:90) dϕ | R − r ± | (cid:20) − I + 3( R − r ± ) ( R − r ± ) | R − r ± | (cid:21) · (cid:18) ∂ r ± ∂κ × d r ± dϕ (cid:19) . (4.38)The gradient of the magnetic field changes according to ∇ X B ( R ) → ∇ X B ( R ) + ∆κ ± ∂ ∇ X B ( R ) /∂κ ± ; from (4.22) and (4.29) we deduce ∂ ∇ X B ∂κ ± = ± µ I heli π (cid:90) dϕ | R − r ± | (cid:20) R − r ± ) I − R − r ± )( R − r ± ) ( R − r ± ) | R − r ± | + I ( R − r ± ) + ˆ e R ( R − r ± ) ˆ e R + ˆ e Z ( R − r ± ) ˆ e Z ] · (cid:18) ∂ r ± ∂κ × d r ± dϕ (cid:19) .(4.39)A gradient-based optimization scheme to demonstrate an application of the adjointgradient formulation was performed on the system with R = 1, r = 0 . µ I heli / π =0 . B t = 1. This was previously optimized in Cary & Hanson (1986) using aderivative-free algorithm. The fixed parameters are A + , = π , A − , = A − , = B + , = 0,and A ± ,k = B ± ,k = 0 for k (cid:62)
2. The vector p = ( A + , , B − , ) is composed of the onlyparameters which are optimized. Two fixed points are initially located and their residuesare used in the optimization. A priori, an infinite choice of appropriate objective functionsexist: the sensible properties are that they be functions of the residues of all the periodicfield lines and that they be zero when all the residues are zero. The island width can onlybe used in the objective function if it is known that the periodic field line is an O pointinstead of an X point. However, periodic field lines may switch from O to X points duringthe optimization as p is varied. Moreover, the measure of island width used in this workis only accurate if the island size is small. For these reasons, using the residues in theoptimization is preferable. We choose a linear combination P of the pairs of residues, P = 12 ( R + R ) . (4.40)One could in principle take the difference of the two residues or the mean square of theresidues as the objective function: the former converges to a non-optimal solution with R = R where the objective function is zero yet the residues have not been reducedenough; the latter has a much slower convergence. The objective function was minimizedusing gradient descent with a line search, with more details discussed in appendix F.The initial configuration with p = (0 . , .
0) has residues R = − .
143 and R = 79 . p = (0 . , . R = 0 . R = 0 . djoint method for sensitivity of island size to magnetic field variations R Z R Z Figure 6.
Poincar´e plots of unoptimized (top) and optimized (bottom) helical coil configuration.The squares indicate fixed points with residue R , while the crosses are those with residue R obtained by Cary & Hanson (1986), but give a slightly more optimized configuration:the small difference may be caused by different numerical resolutions in the Biot-Savartintegral around the helical coils. The optimized and unoptimized configurations are shownin figure 6, and the helical coils that produce them are shown in figure 7. The quadraticconvergence of the errors of the gradients of the residue for an intermediate step duringthe iteration is shown in figure 8.
5. Conclusion
In this paper we have derived the equations for the gradients of several quantitiesrelated to magnetic islands, including the island width and the residue of the periodic fieldline, using an adjoint method. The residue is a quantity that, when small and positive, isstrongly correlated with island size but that can be calculated for any periodic field line,even X points in stochastic regions. Thus, although it does not quantify the physicalwidth of an island, it is more versatile and can be used to minimize stochasticity inmagnetic configurations. The gradient of the island width was obtained by differentiatingthe measure of island width introduced by Cary & Hanson (1991). Although the islandwidth calculation is only accurate for small islands, optimized configurations usuallyhave small islands. Our adjoint approach thus provides an efficient and reliable methodto compute the sensitivity of the island size in optimized or near-optimized stellarators.We have performed and verified numerical gradient calculations on different magneticfield configurations. The analytical configuration of Reiman & Greenside (1986), shown2
A. Geraldini, M. Landreman, E. Paul
Figure 7.
Helical coils corresponding to the initial unoptimized configuration (dotted line) andto the optimized configuration (solid line). Coil 1 (red) carries the negative current ( − I heli ). -8 -6 -4 -2 E R ( δ ) B −,1 (30)A +,1 (30)δ B −,1 (80)A +,1 (80)10 -4 -3 -2 -1 finite difference step size -8 -6 -4 -2 E R ( δ ) Figure 8.
Error convergence of the centred difference approximation of the gradient of thetwo residues to the adjoint calculation when B − , = 0 . A + , = 0 . R − ( ∂ R /∂κ ) δ for κ ∈ { A + , , B − , } .The legend labels which parameter the gradient was taken with respect to for the differentsymbols, and in brackets the value of N ϕ . The higher resolution, N ϕ = 80, has an increasedaccuracy, although the optimization was carried out at N ϕ = 30 to make it faster. The decreaseof the converged error for larger N ϕ indicates discretization error in the adjoint-based approach. djoint method for sensitivity of island size to magnetic field variations Appendix A. Magnetic field line trajectory as a Hamiltonian system
A fundamental property of the magnetic field, namely that it is divergenceless, isimposed by defining B as the curl of a vector potential, B = ∇ × A . (A 1)Note that A is not uniquely defined: adding the gradient of a scalar function, ∇ g , calleda gauge transformation, leaves B unchanged. The form of B we have adopted in (2.1)corresponds to the vector potential A = ψ ∇ θ − χ ( ψ, θ, ϕ ) ∇ ϕ, (A 2)where we have chosen a particular gauge.Considering R ( ϕ ) as the position vector following a field line, equation (A 1) can bederived by extremizing the action S = (cid:90) ϕ ϕ L ( R ( ϕ ) , ϕ ) dϕ (A 3)with respect to R ( ϕ ), where the end-points along the path, R ( ϕ ) and R ( ϕ ), are fixed,and the Lagrangian L is L = A ( R ) · d R dϕ . (A 4)The choice (A 4) can be verified using the Euler-Lagrange equation resulting from the4 A. Geraldini, M. Landreman, E. Paul extremization of the action, ddϕ (cid:18) d L ∂ ˙ R (cid:19) = d L ∂ R , (A 5)which leads to the equation B × d R /dϕ = 0, implying that d R /dϕ is always parallel to B . Upon inserting equation (A 2) into (A 4), the action S in equation (A 3) becomes S = (cid:90) ϕ ϕ (cid:18) ψ dθdϕ − χ ( ψ, θ, ϕ ) (cid:19) dϕ. (A 6)The action expressed in the form (A 6) can be directly compared with the standard formfor the action of a Hamiltonian system, S = (cid:90) t t (cid:18) p dqdt − H ( q, p, t ) (cid:19) dt. (A 7)Hence, it follows that the trajectory of a magnetic field line as a function of toroidal angleconstitutes a Hamiltonian system where the canonical co-ordinate q is θ , the canonicalmomentum p is ψ , the Hamiltonian H is χ and the time t is ϕ .A.1. Hamiltonian after replacing θ with Θ We proceed to obtain the Hamiltonian after the change of variables ( θ, ψ ) → ( Θ, ψ ),where Θ = θ − ι ( ψ ) ϕ is the poloidal angle relative to the unperturbed magnetic fieldline at the flux surface ψ = ψ crossing ϕ = θ = 0. To this end, we re-express dθ/dϕ byapplying the chain rule to the equation θ = Θ + ι ( ψ ) ϕ , dθdϕ = (cid:18) ∂θ∂Θ (cid:19) ϕ dΘdϕ + (cid:18) ∂θ∂ϕ (cid:19) Θ = dΘdϕ + ι ( ψ ) . (A 8)Here, a subscript to the right of the parentheses indicates what variable is being keptconstant in the partial differentiation within the parentheses. By using (A 8), the actionin equation (A 6) becomes S = (cid:90) ϕ ϕ (cid:18) ψ dΘdϕ + ψι ( ψ ) − χ ( ψ, Θ ) (cid:19) dϕ. (A 9)Therefore, the Hamiltonian K in the new variables is given by equation (2.10). Appendix B. Dominance of resonant perturbation
The equations of the magnetic field line trajectory including the perturbation are dθdϕ = ι ( ψ ) + (cid:88) m,n χ (cid:48) m,n ( ψ ) exp ( imθ − inϕ ) , (B 1)and dψdϕ = − (cid:88) m,n χ m,n ( ψ ) im exp ( imθ − inϕ ) . (B 2)The solutions of equations (B 1) and (B 2) are assumed to be approximately given by thefunctions obtained from equations (2.3)-(2.4), ψ (cid:39) ψ + ψ ( ϕ ), (B 3) djoint method for sensitivity of island size to magnetic field variations θ (cid:39) θ i + ι ( ψ ) ϕ + θ ( ϕ ), (B 4)where ψ and θ are additional functions — assumed small — that depend on theperturbation, and θ i is the value of θ at ϕ = 0. Since | χ | (cid:28) χ , we assume that ψ (cid:39) ψ so that ψ and θ are small corrections, and calculate θ by neglecting it in the phase ofthe perturbation, dθ dϕ (cid:39) (cid:88) m,n χ (cid:48) m,n ( ψ ) exp ( imθ i + imι ϕ − inϕ ) . (B 5)This gives θ (cid:39) (cid:88) m,n χ (cid:48) m,n ( ψ ) imι − in exp ( imθ i + imι ϕ − inϕ ) , (B 6)which is divergent when n/m = ι . Hence, we deduce that the effect of the perturbationis dominated by the Fourier modes that coincide (resonate) with the rotational transformof the unperturbed flux surface, n/m = ι . Appendix C. Variation of periodic field line position
The periodic field line position changes due to the variation of the magnetic fieldconfiguration, affecting the circumference and in turn the island width as discussed insection 2. However, the gradient of the poloidal position vector ¯ X k is itself useful in orderto search faster and more reliably for the new position of the periodic field line after asmall change in parameters, e.g. during the optimization considered in section 4.4. Wethus proceed to derive the variation ∆ ¯ X k .We introduce the vector Lagrangian L ¯ X = ¯ X k + (cid:90) π ( k + L ) /n πk/n (cid:18) d X dϕ − V (cid:19) · µ ¯ X k dϕ . (C 1)In (3.30) we have introduced a constraint using the adjoint variable µ ¯ X k (a matrix).Extremization with respect to variations in X ( ϕ ) gives ∆ L ¯ X = ∆ ¯ X k + (cid:90) π ( k + L ) /n πk/n ∆ X · (cid:18) − dµ X k dϕ − ∇ X V · µ ¯ X k (cid:19) dϕ (C 2)+ ∆ X (2 π ( k + L ) /n ) · µ ¯ X k (2 π ( k + L ) /n ) − ∆ X (2 πk/n ) · µ ¯ X k (2 πk/n )) . (C 3)Considering ∆ X (2 π ( k + L ) /n ) = ∆ X (2 πk/n ) = ∆ ¯ X k , the adjoint equation is dµ ¯ X k dϕ = −∇ ¯ X V · µ ¯ X k , (C 4)and must be solved with the boundary condition µ ¯ X k (2 π ( k + L ) /n ) = µ ¯ X k (2 πk/n ) − I .The boundary condition can be simplified as follows. The transpose of (2.46) is d S (cid:124) k dϕ = S (cid:124) k · ∇ X V . (C 5)Premultiplying and postmultiplying by the adjoint of S k , S † k = ( S (cid:124) k ) − , and noting that6 A. Geraldini, M. Landreman, E. Paul S † k · d S (cid:124) k /dϕ = − (cid:16) d S † k /dϕ (cid:17) · S (cid:124) k (from the chain rule and S † k · S (cid:124) k = I ) results in d S † k dϕ = −∇ X V · S † k . (C 6)Comparing equations (C 4) and (C 6) and using that S k (2 π ( k + L ) /n ) = M k with S k (2 πk/n ) = I gives µ ¯ X k (2 π ( k + L ) /n ) = M † k · µ ¯ X (2 πk/n ) . (C 7)Hence, the appropriate boundary condition for equation (C 4) is µ ¯ X k (2 πk/n ) = (cid:104) I − M † k (cid:105) − . (C 8)The variation of ¯ X k is then given by ∆ ¯ X k = − (cid:90) π ( k + L ) /n πk/n ∆ V · µ ¯ X k dϕ. (C 9) Appendix D. Formulas for the Reiman magnetic field configuration
The analytical magnetic field configuration is adapted from a model magnetic fieldused in Reiman & Greenside (1986). The magnetic field components are, from equations(2.1) and (4.3)-(4.6), B R = ZR D + R − R R D , (D 1) B Z = − R − R R D + ZR D , (D 2) B ϕ = −
1, (D 3)where D = ι ax + ι (cid:48) ax r − k max (cid:88) k =1 kε k r k − cos ( kθ − ϕ ) , (D 4) D = k max (cid:88) k =1 kε k r k − sin ( kθ − ϕ ) . (D 5)From equations (D 1)-(D 3), the gradients of the magnetic field components with respectto the poloidal coordinates R and Z are ∂ R B R = − ZR D + ZR ∂ R D + R R D + R − R R ∂ R D , (D 6) ∂ Z B R = 1 R D + ZR ∂ Z D + R − R R ∂ Z D , (D 7) ∂ R B Z = − R R D − R − R R ∂ R D − ZR D + ZR ∂ R D , (D 8) ∂ Z B Z = − R − R R ∂ Z D + 1 R D + ZR ∂ Z D , (D 9) ∂ R B ϕ = ∂ Z B ϕ = 0, (D 10) djoint method for sensitivity of island size to magnetic field variations ∂ R D = 2 ι (cid:48) ax ( R − R ) − k max (cid:88) k =1 kε k r k − (( k −
2) ( R − R ) cos ( kθ − ϕ ) + kZ sin ( kθ − ϕ )) ,(D 11) ∂ R D = k max (cid:88) k =1 kε k r k − (( k −
2) ( R − R ) sin ( kθ − ϕ ) − kZ cos ( kθ − ϕ )) , (D 12) ∂ Z D = 2 ι (cid:48) ax Z − k max (cid:88) k =1 kε k r k − (( k − Z cos ( kθ − ϕ ) − k ( R − R ) sin ( kθ − ϕ )) ,(D 13) ∂ Z D = k max (cid:88) k =1 kε k r k − (( k − Z sin ( kθ − ϕ ) + k ( R − R ) cos ( kθ − ϕ )) . (D 14)From equations (D 6)-(D 10), the second derivatives of the magnetic field componentswith respect to R and Z are ∂ RR B R = 2 ZR D − R R D − ZR ∂ R D + 2 R R ∂ R D + ZR ∂ RR D + R − R R ∂ RR D ,(D 15) ∂ RZ B R = − R D − ZR ∂ Z D + 1 R ∂ R D + R R ∂ Z D + ZR ∂ RZ D + R − R R ∂ RZ D ,(D 16) ∂ ZZ B R = 2 R ∂ Z D + ZR ∂ ZZ D + R − R R ∂ ZZ D , (D 17) ∂ RR B Z = 2 R R D + 2 ZR D − R R ∂ R D − ZR ∂ R D − R − R R ∂ RR D + ZR ∂ RR D ,(D 18) ∂ RZ B Z = − R D − R R ∂ Z D − ZR ∂ Z D + 1 R ∂ R D − R − R R ∂ RZ D + ZR ∂ RZ D ,(D 19) ∂ ZZ B Z = 2 R ∂ Z D − R − R R ∂ ZZ D + ZR ∂ ZZ D , (D 20) ∂ RR B ϕ = ∂ RZ B ϕ = ∂ ZZ B ϕ = 0, (D 21)where ∂ RZ = ∂ ZR and ∂ RR D = 2 ι (cid:48) ax − k max (cid:88) k =1 kε k r k − (cid:2)(cid:0) ( k − (cid:0) r + ( k − R − R ) (cid:1) − k Z (cid:1) cos ( kθ − ϕ )+ kZ ( R − R )(2 k −
6) sin ( kθ − ϕ )] , (D 22) ∂ RZ D = − k max (cid:88) k =1 kε k r k − (cid:2) k (cid:0) r − ( k − R − R ) + ( k − Z (cid:1) sin ( kθ − ϕ )+ kZ ( R − R )(2 k − k + 8) cos ( kθ − ϕ ) (cid:3) , (D 23) ∂ ZZ D = 2 ι (cid:48) ax − k max (cid:88) k =1 kε k r k − (cid:2)(cid:0) ( k − (cid:0) r + ( k − Z (cid:1) − k ( R − R ) (cid:1) cos ( kθ − ϕ ) − kZ ( R − R )(2 k −
6) sin ( kθ − ϕ )] , (D 24)8 A. Geraldini, M. Landreman, E. Paul ∂ RR D = k max (cid:88) k =1 kε k r k − (cid:2)(cid:0) ( k − (cid:0) r + ( k − R − R ) (cid:1) − k Z (cid:1) sin ( kθ − ϕ ) − kZ ( R − R ) (2 k −
6) cos ( kθ − ϕ )] , (D 25) ∂ RZ D = k max (cid:88) k =1 kε k r k − (cid:2) k (cid:0) − r + ( k − R − R ) − ( k − Z (cid:1) cos ( kθ − ϕ )+ kZ ( R − R )(2 k − k + 8) sin ( kθ − ϕ ) (cid:3) , (D 26) ∂ ZZ D = k max (cid:88) k =1 kε k r k − (cid:2)(cid:0) ( k − (cid:0) r + ( k − Z (cid:1) − k ( R − R ) (cid:1) sin ( kθ − ϕ )+ kZ ( R − R )(2 k −
6) cos ( kθ − ϕ )] . (D 27)All the equations in this appendix are linear in the parameters κ ∈ ( ι ax , ι (cid:48) ax , ε ).Hence, the gradients of the magnetic field and its first poloidal derivatives can bestraightforwardly extracted. Appendix E. Shape gradient of coil-produced magnetic field
The magnetic field produced by a set of N coils, indexed c , is given by the Biot-Savartlaw (4.19). Perturbing the coils such that r c ( l c ) → r c ( l c ) + ∆ r c ( l c ) changes the magneticfield such that B ( R ) → B ( R ) + ∆ B ( R ), where the change in magnetic field is given by ∆ B coil = N (cid:88) c =1 µ I c π (cid:73) dl c (cid:18) ∆ r c · ∇ r c (cid:18) R − r c | R − r c | (cid:19) × d r c dl c + (cid:18) R − r c | R − r c | (cid:19) × d∆ r c dl c (cid:19) .(E 1)where ∇ r c (cid:18) R − r c | R − r c | (cid:19) = 1 | R − r c | (cid:18) − I + ( R − r c ) ( R − r c ) | R − r c | (cid:19) . (E 2)An important property of (E 2) is thatTr (cid:20) ∇ r c (cid:18) R − r c | R − r c | (cid:19)(cid:21) = 0 . (E 3)The second term in (E 1) can be re-expressed using integration by parts as (cid:73) dl c (cid:18) R − r c | R − r c | (cid:19) × d∆ r c dl c = − (cid:73) dl c d r c dl c · ∇ r c (cid:18) R − r c | R − r c | (cid:19) × ∆ r c , (E 4)where the term involving an integral of a total derivative has vanished due to the integralbeing periodic. Inserting this into the previous equation gives ∆ B coil = N (cid:88) c =1 µ I c π (cid:73) dl c (cid:18) ∆ r c · ∇ r c (cid:18) R − r c | R − r c | (cid:19) × d r c dl c − d r c dl c · ∇ r c (cid:18) R − r c | R − r c | (cid:19) × ∆ r c (cid:19) .(E 5)It is an exercise in vector identities to show that ∆ r c ·∇ r c (cid:18) R − r c | R − r c | (cid:19) × d r c dl c − d r c dl c · ∇ r c (cid:18) R − r c | R − r c | (cid:19) × ∆ r c = Tr (cid:20) ∇ r c (cid:18) R − r c | R − r c | (cid:19)(cid:21) ∆ r c × d r c dl c − ∇ r c (cid:18) R − r c | R − r c | (cid:19) · ∆ r c × d r c dl c . (E 6) djoint method for sensitivity of island size to magnetic field variations Appendix F. Optimization scheme
The iteration scheme used to search for solutions of P = 0 exploits the direction vector d n obtained from the gradient of P n with respect to the two variable parameters, d n = (cid:18) ∂P n ∂A − , , ∂P n ∂B + , (cid:19) . (F 1)Here, a subscipt n denotes the value of the quantity at the n th step in the iteration.Starting from an initial guess p = (0 . , . p n +1 = p n + p n, step , (F 2)where p n, step = − g n r P n d n | d n | , (F 3)and where g n > r is the number of times themove (F 2) is rejected at a particular iteration. The reference value of g n is denoted ¯ g n and is specified as follows. Initially we set ¯ g = 1. Then, if at any particular iteration thenumber of times r the move (F 2) gets rejected before being accepted exceeds a thresholdvalue r limit = 3, we set ¯ g n +1 = ¯ g n /
2, otherwise we keep the reference value unchanged,¯ g n +1 = ¯ g n . The move (F 2) is rejected if the step in parameter space leads to an increasein the objective function, P n +1 (cid:62) P n .There is an additional rejection criterion that is linked to a maximum value of thestep in the fixed point position, | ¯ X n +1 − ¯ X n | (the subscript labelling the fixed pointis omitted to avoid clutter, and the superscript is the iteration step number). This isnecessary in order to ensure that the iteration proceeds successfully. At the iterationstep n + 1 the scheme to search for a periodic field line is applied with the initial guessequal to ¯ X n − p n, step · ∇ p ¯ X n . The gradient of the fixed point position with respect tothe parameters, ∇ p ¯ X n = ∂ ¯ X n /∂ p , is calculated from the equations of Appendix C. Acommon way for the search for the new fixed point position X n +1 to fail is by returningthe magnetic axis (which is always a periodic solution of the magnetic field line equationsfor any toroidal interval which is a multiple of the the field periodicity). In order to avoidthis, the maximum step size is set to be a fraction u of the smallest distance between thefixed point and the magnetic axis (cid:12)(cid:12) ¯ X n − ¯ X n axis (cid:12)(cid:12) . The value u = 0 . P n +1 < P n and | ¯ X n +1 − ¯ X n | < u (cid:12)(cid:12) ¯ X n − ¯ X n axis (cid:12)(cid:12) ;reject move if P n +1 (cid:62) P n or | ¯ X n +1 − ¯ X n | (cid:62) u (cid:12)(cid:12) ¯ X n − ¯ X n axis (cid:12)(cid:12) . (F 4)Once a move (F 2) is accepted, it becomes permanent and the following move p n +1 , step is calculated. The optimization is stopped when | P n − − P n | < − . REFERENCESAntonsen, T., Paul, E. J. & Landreman, M.
Journal of PlasmaPhysics (2). A. Geraldini, M. Landreman, E. Paul
Cary, J. R. & Hanson, J. D.
The Physics of fluids (8),2464–2473. Cary, J. R. & Hanson, J. D.
Physics ofFluids B: Plasma Physics (4), 1006–1014. Cary, J. R. & Littlejohn, R. G.
Annals of Physics , 1–34.
Greene, J.
Journal of MathematicalPhysics , 1183. Greene, J. M.
Journal of MathematicalPhysics (5), 760–768. Hammond, K., Anichowski, A., Brenner, P., Pedersen, T. S., Raftopoulos, S.,Traverso, P. & Volpe, F.
Plasma Physics and Controlled Fusion (7), 074002. Hanson, J. D. & Cary, J. R.
The Physics offluids (4), 767–769. Hegna, C.
Physics of Plasmas (5), 056101. Helander, P.
Reportson Progress in Physics (8), 087001. Hudson, S., Monticello, D. & Reiman, A.
Physics of Plasmas (7), 3377–3381. Jaenicke, R., Ascasibar, E., Grigull, P., Lakicevic, I., Weller, A., Zippe, M., Hailer,H. & Schworer, K.
Nuclear Fusion (5), 687. Landreman, M. & Paul, E.
Nuclear Fusion (7), 076023. Lazerson, S. A., Bozhenkov, S., Israeli, B., Otte, M., Niemann, H., Bykov, V.,Endler, M., Andreeva, T., Ali, A., Drewelow, P. & others
Plasma Physics and Controlled Fusion (12), 124002. Lee, D., Harris, J. & Lee, G.
Nuclear fusion (10), 2177. Meiss, J.
Reviews of ModernPhysics (3), 795. Neilson, G., Gruber, C., Harris, J. H., Rej, D., Simmons, R. & Strykowsky, R.
IEEE transactions on plasma science (3), 320–327. Paul, E. J., Abel, I. G., Landreman, M. & Dorland, W.
Journal of Plasma Physics (5). Paul, E. J., Antonsen, T., Landreman, M. & Cooper, W. A.
Journal of Plasma Physics (1). Pedersen, T. S., Kremer, J., Lefrancois, R., Marksteiner, Q., Pomphrey, N.,Reiersen, W., Dahlgren, F. & Sarasola, X.
Fusion science and technology (3), 372–381. Pedersen, T. S., Otte, M., Lazerson, S., Helander, P., Bozhenkov, S., Biedermann,C., Klinger, T., Wolf, R. C. & Bosch, H.-S.
Nature communications (1),1–10. Reiman, A. & Greenside, H.
Computer Physics Communications (1), 157–167. Rosenbluth, M., Sagdeev, R., Taylor, J. & Zaslavski, G.
Nuclear Fusion (4), 297. Spitzer Jr, L.
The Physics of Fluids (4), 253–264. Strykowsky, R., Brown, T., Chrzanowski, J., Cole, M., Heitzenroeder, P., Neilson,G., Rej, D. & Viol, M. , pp. 1–4. IEEE.
Yamazaki, K., Yanagi, N., Ji, H., Kaneko, H., Ohyabu, N., Satow, T., Morimoto, S., djoint method for sensitivity of island size to magnetic field variations Yamamoto, J., Motojima, O. & LHD Design Group
Fusion Engineering and Design ,79 – 86. Zhu, C., Gates, D. A., Hudson, S. R., Liu, H., Xu, Y., Shimizu, A. & Okamura, S.
Nuclear Fusion59