unphysical. An alternative is to smooth the coefficients artificially as,

$\displaystyle B_{1}^{\text{smooth}}(r)$
$\displaystyle=B_{1}(r)\text{sign}(n-m\iota(r))\frac{n-m\iota(r)}{\sqrt{(n-m\iota(r))^{2}+\epsilon}}$
(6.59a) $\displaystyle B_{2}^{\text{smooth}}(r)$
$\displaystyle=B_{2}(r)\frac{(n-m\iota(r))^{2}}{(n-m\iota(r))^{2}+\epsilon},$
(6.59b)

where $\epsilon\ll 1$ is a scalar chosen to account for the smoothing by
classical diffusion. When $\epsilon\rightarrow 0$, the Euler-Lagrange equation
remains unchanged. For small but finite $\epsilon$, the coefficients are only
modified in the vicinity of $r_{s}$. This is similar to a technique used in
the IPEC [181] code.

##### Analytic solution near singular surfaces

To study the solutions of the Euler-Lagrange equation with $m\neq 0$ and
$n\neq 0$ further, we consider a limit in which analytic solutions can be
obtained. We will take $p^{\prime}(\psi)=0$ and $\iota(r)=\iota_{1}r$ where
$\iota_{1}$ is a constant. In this case the force-balance equation (6.35)
gives us the following expression for the flux in terms of hypergeometric
functions,

$\displaystyle\psi(r)=\frac{r^{2}\psi_{0}\prescript{}{2}{F}_{1}\left(\frac{1}{2};\frac{3}{4};\frac{3}{2};-\frac{r^{4}\iota_{1}^{2}}{R_{0}^{2}}\right)}{\prescript{}{2}{F}_{1}\left(\frac{1}{2};\frac{3}{4};\frac{3}{2};-\frac{\iota_{1}^{2}}{R_{0}^{2}}\right)}.$
(6.60)

We define a variable $r_{s}=n/(m\iota_{1})$ such that a singular surface
occurs at $r=r_{s}$. The coefficients of the homogeneous problem can be
expressed as,

$\displaystyle B_{1}(r)$
$\displaystyle=\frac{3}{r-r_{s}}-\frac{5r_{s}}{r^{2}-rr_{s}}-\frac{6}{r+r^{5}\iota_{1}^{2}/R_{0}^{2}}-\frac{2R_{0}^{2}}{rR_{0}^{2}+r^{3}r_{s}^{2}\iota_{1}^{2}}$
(6.61a) $\displaystyle B_{2}(r)$
$\displaystyle=\frac{1+m^{2}}{r^{2}}+\frac{4}{r_{s}r-r^{2}}+\frac{m^{2}r_{s}^{2}\iota_{1}^{2}}{R_{0}^{2}}+\frac{2R_{0}^{2}(r+r_{s})}{r^{2}(r-r_{s})(R_{0}^{2}+r^{2}r_{s}^{2}\iota_{1}^{2})}.$
(6.61b)

In the limit of small shear, $\epsilon_{\iota}=\iota_{1}r_{s}^{2}/R_{0}\ll 1$,
we can approximate the coefficients as,

$\displaystyle B_{1}(r)$
$\displaystyle=\frac{3r_{s}-5r}{r(r-r_{s})}+\mathcal{O}\left(\epsilon_{\iota}^{2}\right)$
(6.62a) $\displaystyle B_{2}(r)$
$\displaystyle=\frac{m^{2}-1}{r^{2}}+\mathcal{O}\left(\epsilon_{\iota}^{2}\right).$
(6.62b)

In practice we choose a very small value for this expansion parameter
($\epsilon_{\iota}\sim 10^{-4}$) so that dropping the higher order terms is a
very good approximation. For the $m=2$, $n=1$ mode subject to a boundary
perturbation,

$\displaystyle\xi_{2,1}^{rc}(1)$ $\displaystyle=1$ $\displaystyle\delta
F^{2,1}_{rc}(r)$ $\displaystyle=0,$ (6.63)

we have the analytic solution,

$\displaystyle\xi^{rc}_{2,1}(r)=r\text{Re}\left[\frac{\prescript{}{2}{F}_{1}\left(3-\sqrt{7};3+\sqrt{7},5,\frac{r}{r_{s}}\right)}{\prescript{}{2}{F}_{1}\left(3-\sqrt{7};3+\sqrt{7},5,\frac{1}{r_{s}}\right)}\right].$
(6.64)

We first consider the case in which $r_{s}=2$ such that a singular surface
does not appear within the volume. We compare the numerical solution of the
Euler-Lagrange equation with a finite-difference calculation with VMEC. We
impose a boundary perturbation of the form,

$\displaystyle\delta R(\psi_{0},\theta_{b},\phi)$
$\displaystyle=\Delta\cos(3\theta_{b}-\phi)$ (6.65a) $\displaystyle\delta
Z(\psi_{0},\theta_{b},\phi)$ $\displaystyle=\Delta\sin(3\theta_{b}-\phi),$
(6.65b)

where $\phi$ is the geometric toroidal angle. The perturbed field is computed
with a two-point centered-difference stencil with amplitude $\Delta=10^{-4}$.
The results of the calculations are shown in Figure 6.7. We note that the
Euler-Lagrange solution agrees well with the analytic solution, with a volume-
averaged difference of $\Delta_{V}=1.86\times 10^{-3}$, but there is a small
discrepancy between the VMEC solution and the analytic solution near the edge,
with a volume-averaged difference of $\Delta_{V}=9.60\times 10^{-3}$. One
possible source of this error is the treatment of singularities by the VMEC
code. While recent results have indicated that VMEC equilibria can exhibit
$1/x$-like behavior near rational surfaces [144, 160], the numerical solution
is not truly singular on such surfaces, and very large numerical resolution is
necessary in order to see behavior resembling a singularity. Therefore, we do
not expect the displacement vector computed with finite-difference VMEC to
agree with the Euler-Lagrange solution. Although for this equilibrium, $\iota$
does not resonate with the harmonics of the displacement vector, it may
resonate with other modes present in the nonlinear equilibrium.

Figure 6.7: Benchmark of screw pinch $m=2$, $n=1$ mode with a boundary
perturbation (6.63). The solution of the Euler-Lagrange equation (6.36) with
coefficients (6.57) is compared with an analytic solution (6.64) and a finite-
difference calculation from VMEC. This equilibrium does not contain a resonant
surface within the volume.

Next we consider an equilibrium with a singular surface in the volume,
$r_{s}=0.5$. The Euler-Lagrange equation is solved with both the power-series
method, which captures the singular nature of the solution, and the
coefficient smoothing method (6.59) with several values of $\epsilon$. Again,
we compare with a finite-difference VMEC solution with a boundary perturbation
given by (6.65). With the power-series method, we find agreement between the
Euler-Lagrange and analytic solutions. As expected, the solutions with
smoothed coefficients do not reproduce the analytic expression. However,
neither of these approaches approximates the VMEC solution well. Although the
VMEC equilibrium is fairly well-resolved (701 flux surfaces, $10^{-12}$ force
tolerance, $m\leq 4$ poloidal modes, $|n|\leq 4$ toroidal modes), we do not
observe a response near $r=r_{s}$. We may need to consider a revised treatment
of the singularity to match the behavior from VMEC better.

Figure 6.8: Benchmark of screw pinch $m=2$, $n=1$ mode with a boundary
perturbation (6.63). The solution of the Euler-Lagrange equation (6.36) with
coefficients (6.57) is compared with an analytic solution (6.64) and a finite-
difference calculation from VMEC. This equilibrium contains a resonant surface
at $r=0.5$ ($\psi=0.20$).

###  6.4 Tokamak shape gradient

We will now demonstrate the linearized equilibrium technique to compute the
shape gradient of the vacuum magnetic well figure of merit discussed in
Chapter 5,

$\displaystyle f_{W}(S_{P})=\int_{V_{P}}d^{3}x\,w(\psi),$ (6.66)

with,

$\displaystyle
w(\psi)=\exp(-(\psi-\psi_{m,1})^{2}/\psi_{w}^{2})-\exp(-(\psi-\psi_{m,2})^{2}/\psi_{w}^{2}),$
(6.67)

where $\psi_{m,1}=0.9\psi_{0}$, $\psi_{m,2}=0.1\psi_{0}$, and
$\psi_{w}=0.05\psi_{0}$. The shape gradient of $f_{W}$ is obtained with an
adjoint approach by computing a perturbed equilibrium state corresponding to
the addition of a bulk force with no displacement of the boundary,

$\displaystyle\delta\textbf{x}\cdot\nabla\psi$ $\displaystyle=0$
$\displaystyle\delta\textbf{F}=-\nabla w(\psi).$ (6.68)

The resulting perturbed field, $\delta\textbf{B}[\bm{\xi}]$, is used to
compute the shape gradient,

$\displaystyle\mathcal{G}=\frac{\delta\textbf{B}[\bm{\xi}]\cdot\textbf{B}}{\mu_{0}}\bigg{\rvert}_{S_{P}}.$
(6.69)

We perform this calculation for an axisymmetric configuration with a plasma
boundary given by,

$\displaystyle R(\psi_{0},\theta_{b})$
$\displaystyle=R_{0}+a\cos(\theta_{b})+b\cos(2\theta_{b})$ (6.70a)
$\displaystyle Z(\psi_{0},\theta_{b})$
$\displaystyle=a\sin(\theta_{b})-b\sin(2\theta_{b}),$ (6.70b)

with $R_{0}=3$, $a=1$, and $b=0.1$. Owing to its toroidal symmetry, all of the
toroidal modes of the displacement vector decouple. Given the toroidal
symmetry of the bulk force perturbation, we only need to consider the $n=0$
modes. Therefore, the only singular point of the Euler-Lagrange equation is at
the origin. As before, the magnetic axis is not included on the computational
grid, and the coupled BVP is solved with the bvp4c routine. The radial
displacement vector is computed retaining modes $m\leq 30$.

The resulting shape gradient obtained from the Euler-Lagrange solution is
shown in Figure 6.9 along with that computed with a finite-difference VMEC
calculation,

$\displaystyle\delta p(\psi)=\Delta w(\psi).$ (6.71)

A two-point centered-difference derivative is computed with magnitude
$\Delta=10$. The surface-averaged fractional difference between the Euler-
Lagrange and VMEC solutions is computed to be $7.3\times 10^{-3}$.

(a)

(b)

(c)

Figure 6.9: The shape gradient of the vacuum magnetic well (6.66) is computed
for a tokamak equilibrium with triangularity (6.70) with the solution of the
Euler-Lagrange equation corresponding to the adjoint problem (6.68) and a
finite-difference approximation of the adjoint problem with VMEC (6.71).

###  6.5 Conclusions

We have demonstrated a variational method for computing perturbed equilibrium
states corresponding to the addition of a bulk force or boundary perturbation.
We considered the simplified geometry of a screw pinch to demonstrate the
behavior of each of the modes of the displacement vector. Numerical solutions
of the Euler-Lagrange equation are benchmarked with finite-difference
calculations of the nonlinear equilibrium code, VMEC, and with analytic
solutions in certain limits. Finally, we employed this approach to compute the
shape gradient of a figure of merit of interest for stellarator optimization
in toroidally symmetric geometry. We aim to apply this approach for computing
such shape gradients in stellarator geometry, though this task may be somewhat
more challenging. In fully 3D geometry, there may exist several singular
surfaces throughout a volume due to toroidal mode coupling, each of which
needs to be treated carefully,

While the Euler-Lagrange equation exhibits singular behavior at rational
surfaces, the equilibria computed with the VMEC code do not appear to exhibit
any singular response, as demonstrated in Section 6.3.4. If the goal is to
linearize about VMEC equilibria, we therefore may not want to solve the Euler-
Lagrange equation exactly, but to artificially smooth the coefficients
appearing in the ODE. As an alternative, artificial viscosity could be added
to the Euler-Lagrange system with the addition of a small term involving a
higher-order derivative. This technique, commonly used in the fluid dynamics
community [156, 67], turns a singular ODE into an ODE with a singular
perturbation. It remains to be demonstrated that the shape gradients obtained
from Euler-Lagrange solutions including such smoothing techniques can
reproduce the expected shape gradients computed with the VMEC code.

In addition to the demonstration for three-dimensional geometry, there are
several interesting extensions of the work discussed in this Chapter. As
discussed in Chapter 5, there are several figures of merit for which the
adjoint problem requires the addition of a perturbation to the prescribed
toroidal current profile. This would necessitate generalizing this formulation
to allow for perturbations to the magnetic field that vary the rotational
transform profile. While the work in this Chapter has been applied to compute
the shape gradient with respect to the plasma boundary, it may be possible to
couple perturbations of the boundary to coil perturbations in order to compute
the coil shape gradient. This may benefit from a method similar to that used
in the IPEC code, in which the virtual casing principle is applied to couple
boundary perturbations to changes in the external magnetic fields.

The further development of this linear equilibrium approach would enable the
shape gradient of many additional figures of merit to be computed with an
adjoint method. Even if an adjoint method is not applied, the linear
equilibrium approach could prove very fruitful for gradient-based, fixed-
boundary optimization. Replacing a finite-difference calculation by an
analytic derivative may reduce computational cost and noise associated with
the finite-difference step size, enabling more efficient sensitivity and
tolerance calculations for stellarator configurations.

## Chapter 7: Conclusions

In this Thesis, we have aimed to address fundamental challenges (Section
1.4.4) associated with stellarator optimization using the adjoint method and
shape sensitivity analysis:

  1. 1.

Coil complexity

  2. 2.

Non-convexity

  3. 3.

High-dimensionality

  4. 4.

Tight engineering tolerances.

The adjoint method allows us to efficiently compute derivatives in the context
of several problems of interest for stellarator optimization. These
derivatives enable navigation through high-dimensional, non-convex spaces with
gradient-based methods. We demonstrate gradient-based optimization with
adjoints in Chapter 3, for the design of coil shapes with minimal complexity.
Computing the shape gradient of coil metrics to perturbations of the winding
surface allows us to gain intuition about features of configurations that
enable simpler coils. We also demonstrate gradient-based optimization of the
local magnetic geometry for finite-collisionality neoclassical properties in
Chapter 5. While including such objective functions is typically prohibitively
expensive for non-convex, high-dimensional optimization, we demonstrate
convergence toward a local optimum with a minimal number of function
evaluations. With this adjoint method, we also gain intuition of the
sensitivity of the bootstrap current and particle fluxes to perturbations in
the field strength, informing engineering tolerances. Finally, in Chapter 5 we
demonstrate an adjoint method for computing the plasma surface and coil shape
gradient for functions that depend on MHD equilibrium solutions. Importantly,
the coil shape gradient can be used to evaluate engineering tolerances for
such figures of merit (Section 2.1.3). While it has not yet been demonstrated
in this Thesis, these shape gradients can also enable efficient adjoint-based
optimization, either in the space of the plasma boundary or coil shapes. As
discussed in Section 1.4, the direct optimization of coil shapes may result in
coils that can be more feasibly engineered than those resulting from the
traditional two-step optimization.

For several problems discussed in this Thesis, it is convenient to apply the
discrete adjoint method (Section 2.2.1). For the winding surface optimization
problem in Chapter 3, the forward problem is solved as a discrete linear
system, so the discrete adjoint operator can be obtained by simply taking the
matrix transpose. A similar discrete adjoint method was applied for
neoclassical optimization in Chapter 4, as the discretized form of the drift-
kinetic equation takes the form of a linear system in the SFINCS code.

Physical insight into the structure of the relevant equations can inform the
development of continuous adjoint methods (Section 2.2.2). For the
neoclassical application, the adjoint equation was obtained based on an inner
product similar to the free-energy norm from gyrokinetic theory. The self-
adjointness of the linear Fokker-Planck operator with respect to this inner
product enabled straightforward calculation of the adjoint operator. For the
MHD application, the adjoint equation is obtained by noting the self-
adjointness of the MHD force operator, generalized to allow for perturbations
of the rotational transform and currents in the vacuum region. Finally, in
Chapter 6, a variational method for solving the adjoint equations obtained in
Chapter 5 is presented. Here we are able to borrow a variational method from
MHD stability theory to efficiently compute the adjoint equilibrium problem.

###  7.1 Outlook

There are several natural extensions of the work presented in this Thesis.

####  7.1.1 Further development of adjoint methods

  * •

The advancement of the adjoint approach for functions of MHD equilibria
necessitates the further development of a linearized equilibrium code, as
outlined in Chapter 6. While we have demonstrated this technique for
axisymmetric equilibria, we plan to extend it to 3D equilibria. In this way,
adjoint methods for computing the shape gradient of the departure from quasi-
symmetry (Section 5.5.6), effective ripple (Section 5.5.5), and several
finite-collisionality neoclassical quantities (Section 5.5.7) could be
demonstrated.

  * •

In Chapter 3, we applied the adjoint method to compute derivatives with
respect to the winding surface parameters. Similarly, we can apply the adjoint
method to compute derivatives with respect to plasma surface parameters. This
would allow for the identification of plasma surfaces that do not require
overly-complex coils, facilitating the incorporation of coil considerations in
plasma configuration optimization [36]. Similar figures of merit (without
derivative information) have been used in the ROSE code [59].

####  7.1.2 Further application of derivatives

We have not yet taken full advantage of derivative information for stellarator
optimization problems.

  * •

The analysis of sensitivity and tolerances presented in this Thesis is based
on a local model, using a linear approximation of a function with first
derivative information. A more accurate global analysis can be computed from
Monte-Carlo sampling, which typically requires many function evaluations to
converge. Uncertainty quantification can be accelerated through the
application of a surrogate model of the design space [238] with the
incorporation of the uncertainty of the data. A surrogate model is an
approximation to an expensive simulation based on a small number of
evaluations of the function. The number of required evaluations to build the
surrogate is reduced with a gradient-enhanced Gaussian process regression
model [146]; thus the availability of adjoint-based gradients would enable
more accurate uncertainty quantification. In addition to sensitivity analysis,
once a surrogate is constructed, it can replace the expensive model during
optimization, allowing for more efficient local or global optimization.

  * •

In particular, one type of surrogate function of interest is a neural network,
which can be trained more efficiently using derivative information. Neural
networks with certain choices of activation functions are differentiable, and
can therefore be optimized with gradient-based optimization techniques.
Gradient-based shape optimization with neural networks has proven fruitful in
the field of aerodynamics [222].

  * •

Optimization under uncertainty methods optimize the expected value of an
objective function by performing a sample average over a distribution of
possible deviations. These techniques can improve the robustness of the
optimum by avoiding small local minima and obtaining solutions with reduced
risk. This technique has proven effective for the optimization of coil shapes
with increased tolerances [150, 151], using a Monte-Carlo approach. To avoid
the excessive cost of a Monte-Carlo method, a linear or quadratic
approximation can be made such that the expectation value and variance can be
computed with derivative information [3] obtained with an adjoint method.

We look forward to the adoption of adjoint methods and shape optimization
tools for many stellarator design problems.

## Appendix A: Toroidal coordinate systems

In this Appendix, we briefly review coordinate systems for describing scalar
and vector fields in toroidal systems. Comprehensive introductions to this
topic are provided in the textbook [54], the review article [97], and the
tutorial [121].

###  A.1 Toroidal coordinates

In this Thesis, we often want to describe surfaces of toroidal topology or the
volumes enclosed by such surfaces. We can describe the position on a toroidal
surface by two angles (Figure A.1). A poloidal angle, denoted by $\theta$,
increases by $2\pi$ upon one rotation the short way around the torus. A
toroidal angle, denoted by $\phi$, increases by $2\pi$ upon one rotation the
long way around the torus.

Figure A.1: The position on a toroidal surface, $S$, is described by the
toroidal and poloidal angles. Figure adapted from [121].

We will consider a volume, $V$, bounded by a toroidal surface, $S$. Suppose
that we use a set of continuously nested toroidal surfaces, $\Gamma(r)$, as a
radial coordinate $r$, such that the position within this volume can be
expressed as $\textbf{x}(r,\theta,\phi)$. A vector field, A can be expressed
in the basis of the gradients of the coordinates,

$\displaystyle\textbf{A}=A_{r}\nabla
r+A_{\theta}\nabla\theta+A_{\phi}\nabla\phi,$ (A.1)

the covariant form, or the derivatives of the position vectors with respect to
the coordinates,

$\displaystyle\textbf{A}=A^{r}\frac{\partial\textbf{x}}{\partial
r}+A^{\theta}\frac{\partial\textbf{x}}{\partial\theta}+A^{\phi}\frac{\partial\textbf{x}}{\partial\phi},$
(A.2)

the contravariant form. The two basis vectors can be related through the dual
relations,

$\displaystyle\frac{\partial\textbf{x}}{\partial x_{i}}=\frac{\nabla
x_{j}\times\nabla x_{k}}{\nabla x_{i}\cdot\nabla x_{j}\times\nabla x_{k}},$
(A.3)

where $(x_{i},x_{j},x_{k})=(r,\theta,\phi)$ or cyclic permutations. Such a
coordinate system is generally non-orthogonal, so $\partial\textbf{x}/\partial
x_{i}$ is not necessarily parallel to $\nabla x_{i}$. Several useful relations
in non-orthogonal coordinate systems are summarized in Table A.1. For a more
detailed discussion, refer to Chapter 2 in [54].

Jacobian | $\sqrt{g}=\left(\frac{\partial\textbf{x}}{\partial x_{i}}\times\frac{\partial\textbf{x}}{\partial x_{j}}\right)\cdot\frac{\partial\textbf{x}}{\partial x_{k}}=\left(\left(\nabla x_{i}\times\nabla x_{j}\right)\cdot\nabla x_{k}\right)^{-1}$
---|---
Differential volume | $d^{3}x=|\sqrt{g}|dx_{i}dx_{j}dx_{k}$
Differential length | $d\textbf{x}=\sum_{i=1}^{3}\frac{\partial\textbf{x}}{\partial x_{i}}dx_{i}$
Differential surface area (constant $x_{k}$) | $d^{2}x=|\sqrt{g}||\nabla x_{k}|dx_{i}dx_{j}$
Divergence of vector field | $\nabla\cdot\textbf{A}=\sum_{i=1}^{3}\frac{1}{\sqrt{g}}\frac{\partial}{\partial x_{i}}\left(\sqrt{g}A^{i}\right)$
Curl of vector field | $\nabla\times\textbf{A}=\sum_{k=1}^{3}\frac{1}{\sqrt{g}}\left(\frac{\partial A_{j}}{\partial x_{i}}-\frac{\partial A_{i}}{\partial x_{j}}\right)\frac{\partial\textbf{x}}{\partial x_{k}}$
Gradient of scalar | $\nabla q=\sum_{i=1}^{3}\frac{\partial q}{\partial x_{i}}\nabla x_{i}$
Table A.1: Summary of formulas used to describe the geometry of a non-
orthogonal coordinate system $(x_{1},x_{2},x_{3})$. In the above,
$\\{i,j,k\\}$ is a cyclic permutation of $\\{1,2,3\\}$. Table adapted from
[121].

###  A.2 Flux coordinates

If magnetic surfaces exist, indicating that the magnetic field is tangent to a
set of continuously nested toroidal surfaces, we can use the toroidal flux
through such surfaces as a coordinate, defined as,

$\displaystyle
2\pi\psi\equiv\int_{\mathcal{S}_{T}(\psi)}d^{2}x\,\textbf{B}\cdot\hat{\textbf{n}}.$
(A.4)

In the above expression, $\mathcal{S}_{T}(\psi)$ is an open surface such that
$\partial\mathcal{S}_{T}(\psi)$ is a loop on $\Gamma(\psi)$ that closes after
one poloidal rotation (Figure A.2). The unit normal is $\hat{\textbf{n}}$,
often chosen to point in the direction of increasing $\phi$.

Figure A.2: The plasma domain, $V_{P}$, is bounded by a toroidal surface,
$S_{P}$. We make the assumption that there exists a set of toroidal magnetic
surfaces, $\Gamma(\psi)$. The toroidal flux through each of these surfaces is
defined by (A.4) with $\mathcal{S}_{T}(\psi)$ an open surface bounded by a
poloidally closed curve on $\Gamma(\psi)$, $\partial\mathcal{S}_{T}(\psi)$.

Another choice for labeling magnetic surfaces is the poloidal flux function,
$\chi$,

$\displaystyle
2\pi\chi\equiv\int_{\mathcal{S}_{P}(\psi)}d^{2}x\,\textbf{B}\cdot\hat{\textbf{n}},$
(A.5)

where $\mathcal{S}_{P}(\psi)$ is an open surface such that
$\partial\mathcal{S}_{P}(\psi)$ is a loop on $\Gamma(\psi)$ that closes after
one toroidal rotation (Figure A.3).

Figure A.3: The poloidal flux through the magnetic surface, $\Gamma(\psi)$, is
defined by (A.5) with $\mathcal{S}_{P}(\psi)$ an open surface bounded by a
toroidally closed curve on $\Gamma(\psi)$, $\partial\mathcal{S}_{P}(\psi)$.

The rotational transform quantifies the number of poloidal turns of a field
line per toroidal turn,

$\displaystyle\iota\equiv\lim_{n\rightarrow\infty}\frac{\sum_{k=1}^{n}\left(\Delta\theta\right)_{k}}{2\pi
n}.$ (A.6)

Here $(\Delta\theta)_{k}$ is the change in poloidal angle in toroidal rotation
$k$ and $n$ counts the toroidal turns. If flux surfaces exist, then the
rotational transform can be computed from the derivative of the poloidal flux
with respect to the toroidal flux,

$\displaystyle\iota(\psi)=\chi^{\prime}(\psi),$ (A.7)

If a flux label, $\psi$, is used as one of the coordinates, known as a flux
coordinate system, then the contravariant form for the magnetic field
simplifies,

$\displaystyle\textbf{B}=B^{\theta}\frac{\partial\textbf{x}}{\partial\theta}+B^{\phi}\frac{\partial\textbf{x}}{\partial\phi},$
(A.8)

from the assumption that $\textbf{B}\cdot\nabla\psi=0$. Given
$\nabla\cdot\textbf{B}=0$ and using (A.3), we can express the magnetic field
as,

$\displaystyle\textbf{B}=\nabla\psi\times\nabla\left(\theta-\iota(\psi)\phi+\lambda(\psi,\theta,\phi)\right),$
(A.9)

where $\lambda(\psi,\theta,\phi)$ is $2\pi$-periodic in $\theta$ and $\phi$
(Section 11.1 in [121]).

In a flux-coordinate system, the flux-surface average,

$\displaystyle\langle
A\rangle_{\psi}=\frac{\int_{0}^{2\pi}d\theta\int_{0}^{2\pi}d\phi\,\sqrt{g}A}{V^{\prime}(\psi)},$
(A.10)

appears in many calculations, where

$\displaystyle
V^{\prime}(\psi)=\int_{0}^{2\pi}d\theta\int_{0}^{2\pi}d\phi\,\sqrt{g},$ (A.11)

is the differential volume associated with a change in flux. The flux-surface
average can be equivalently defined as the average over the infinitesimal
volume between flux surfaces,

$\displaystyle\left\langle A\right\rangle_{\psi}=\lim_{\Delta V\rightarrow
0}\frac{1}{\Delta V}\left(\int_{V_{P}(\psi)+\Delta
V}d^{3}x\,A-\int_{V_{P}(\psi)}d^{3}x\,A\right),$ (A.12)

where $V_{P}(\psi)$ is the volume enclosed by a surface labeled by $\psi$ and
$V_{P}(\psi)+\Delta V$ is the volume of a neighboring surface. The flux-
surface average is discussed in more detail in Section 4.9 of [54].

###  A.3 Magnetic coordinates

A flux coordinate system can be defined with many choices of poloidal and
toroidal angles. With some choices of these angles, the contravariant
expression for the magnetic field can simplify further. Given (A.9), the
definition of the poloidal and toroidal angles can be shifted to $\vartheta$
and $\varphi$ such that the magnetic field can be expressed as,

$\displaystyle\textbf{B}=\nabla\psi\times\nabla\left(\vartheta-\iota(\psi)\varphi\right).$
(A.13)

Such angles define a magnetic coordinate system. For example, one choice is
$\vartheta=\theta+\lambda(\psi,\theta,\phi)$ and $\varphi=\phi$. For any
choice of $\varphi$, there is a corresponding choice of $\vartheta$ that
defines a magnetic coordinate system. With this choice of angles, the magnetic
field lines are said to be straight in the $\vartheta-\varphi$ plane,

$\displaystyle\frac{d\vartheta(l)}{d\varphi(l)}=\frac{\textbf{B}\cdot\nabla\vartheta}{\textbf{B}\cdot\nabla\varphi}=\iota(\psi),$
(A.14)

with a slope given by the rotational transform. Here $l$ measures length along
a field line such that $df/dl=\hat{\textbf{b}}\cdot\nabla f$ for any quantity
$f$, where $\hat{\textbf{b}}=\textbf{B}/B$ is the unit vector in the direction
of the magnetic field.

From the covariant form for the magnetic field,

$\displaystyle\textbf{B}=B_{\vartheta}\nabla\vartheta+B_{\varphi}\nabla\varphi+B_{\psi}\nabla\psi,$
(A.15)

we can compute the net toroidal and poloidal currents enclosed by the surface
labeled by $\psi$,

$\displaystyle I_{T}(\psi)$
$\displaystyle\equiv\int_{\mathcal{S}_{T}(\psi)}d^{2}x\,\textbf{J}\cdot\hat{\textbf{n}}=\frac{1}{\mu_{0}}\oint_{\partial\mathcal{S}_{T}(\psi)}d\textbf{l}\cdot\textbf{B}=\frac{1}{\mu_{0}}\int_{0}^{2\pi}d\vartheta\,B_{\vartheta}$
(A.16a) $\displaystyle I_{P}(\psi)$
$\displaystyle\equiv\int_{\mathcal{S}_{P}(\psi)}d^{2}x\,\textbf{J}\cdot\hat{\textbf{n}}=\frac{1}{\mu_{0}}\oint_{\partial\mathcal{S}_{P}(\psi)}d\textbf{l}\cdot\textbf{B}=\frac{1}{\mu_{0}}\int_{0}^{2\pi}d\varphi\,B_{\varphi},$
(A.16b)

where $\mathcal{S}_{T}$ is defined in Figure A.2 and $\mathcal{S}_{P}$ is
defined in Figure A.3. Under the additional assumption that
$\textbf{J}\cdot\nabla\psi=0$, which follows from MHD force balance (1.3a)
with $p(\psi)$, we can write the covariant form as,

$\displaystyle\textbf{B}=I(\psi)\nabla\vartheta+G(\psi)\nabla\varphi+K(\psi,\vartheta,\varphi)\nabla\psi+\nabla
H(\psi,\vartheta,\varphi),$ (A.17)

where $I(\psi)=\mu_{0}I_{T}(\psi)/(2\pi)$ and
$G(\psi)=\mu_{0}I_{P}(\psi)/(2\pi)$. See Section 2.5 in [97], Section 9.2 in
[121], and Chapter 6.5 of [54] for details.

###  A.4 Boozer coordinates

As previously mentioned, there are many choices of magnetic coordinates
corresponding to different choices of toroidal angle, $\varphi$. Suppose we
begin with a system defined by $(\psi,\vartheta,\varphi)$ and want to
transform for a system defined by
$(\psi,\vartheta^{\prime},\varphi^{\prime})$. In order for the primed system
to remain a magnetic coordinate system, we must have
$\varphi^{\prime}=\varphi+\gamma(\psi,\vartheta,\varphi)$ and
$\vartheta^{\prime}=\vartheta+\iota(\psi)\gamma(\psi,\vartheta,\varphi)$,
where $\gamma(\psi,\vartheta,\varphi)$ is $2\pi$-periodic in $\vartheta$ and
$\varphi$. To construct the Boozer coordinate system [23], we will make a
particular choice for $\gamma$ to simplify the covariant form for the magnetic
field (A.17). The corresponding changes to the quantities appearing in the
covariant form (A.17) are

$\displaystyle H^{\prime}$
$\displaystyle=H-\left(\iota(\psi)I(\psi)+G(\psi)\right)\gamma(\psi,\vartheta,\varphi)$
(A.18a) $\displaystyle K^{\prime}$
$\displaystyle=K+\gamma(\psi,\vartheta,\varphi)\left(\iota(\psi)I^{\prime}(\psi)+G^{\prime}(\psi)\right).$
(A.18b)

Boozer coordinates are defined such that $H^{\prime}=0$, or
$\gamma(\psi,\vartheta,\varphi)=H(\psi,\vartheta,\varphi)/(\iota(\psi)I(\psi)+G(\psi))$.
With this choice of transformation, we will denote
$\vartheta_{B}=\vartheta+\iota\gamma$ and $\varphi_{B}=\varphi+\gamma$. The
covariant form becomes,

$\displaystyle\textbf{B}=I(\psi)\nabla\vartheta_{B}+G(\psi)\nabla\varphi_{B}+K(\psi,\vartheta_{B},\varphi_{B})\nabla\psi.$
(A.19)

By dotting the covariant with the contravariant form, we obtain an expression
for the Jacobian,

$\displaystyle\sqrt{g}=\frac{1}{\nabla\psi\times\nabla\vartheta_{B}\cdot\nabla\varphi_{B}}=\frac{G(\psi)+\iota(\psi)I(\psi)}{B^{2}}.$
(A.20)

We note that the Jacobian only varies on a surface through the magnetic field
strength; thus each of the contravariant and covariant components of the
magnetic field, except for $K(\psi,\vartheta_{B},\varphi_{B})$, possesses the
same property. (The radial covariant component,
$K(\psi,\vartheta_{B},\varphi_{B})$, is related to the field strength through
the MHD force balance equation (1.3a).) For this reason, the Boozer coordinate
system is extremely convenient for analyzing guiding center motion and
neoclassical transport, as we will in Chapter 4.

## Appendix B: Justification for current potential

In this Appendix, we justify the form for a continuous current density
supported on a toroidal surface, $S_{C}$,

$\displaystyle\textbf{J}_{C}(\theta,\phi)=\hat{\textbf{n}}\times\nabla\Phi,$
(B.1)

where $\hat{\textbf{n}}$ is the unit normal vector.

We consider an extension of $\textbf{J}_{C}$ in a neighborhood of $S_{C}$ of
width $\Delta b$,

$\displaystyle\widetilde{\textbf{J}}_{C}(b,\widetilde{\theta},\widetilde{\phi})=\textbf{J}_{C}(\theta,\phi),$
(B.2)

where we define extensions of $\theta$ and $\phi$ as,

$\displaystyle\widetilde{\theta}(\textbf{x})$
$\displaystyle=\theta(\textbf{x}-b(\textbf{x})\nabla b)$ (B.3a)
$\displaystyle\widetilde{\phi}(\textbf{x})$
$\displaystyle=\phi(\textbf{x}-b(\textbf{x})\nabla b),$ (B.3b)

or a normal projection onto $S_{C}$. We consider $b\in[-\frac{\Delta
b}{2},\frac{\Delta b}{2}]$ to be a “thickened” region of continuous current
density. We impose the constraint that
$\nabla\cdot\widetilde{\textbf{J}}_{C}=0$, expressed in the
$(b,\widetilde{\theta},\widetilde{\phi})$ coordinate system (Table A.1),

$\displaystyle\frac{1}{\sqrt{g}}\left(\frac{\partial\left(\sqrt{g}\widetilde{\textbf{J}}_{C}\cdot\nabla
b\right)}{\partial
b}+\frac{\partial\left(\sqrt{g}\widetilde{\textbf{J}}_{C}\cdot\nabla\widetilde{\theta}\right)}{\partial\widetilde{\theta}}+\frac{\partial\left(\sqrt{g}\widetilde{\textbf{J}}_{C}\cdot\nabla\widetilde{\phi}\right)}{\partial\widetilde{\phi}}\right)=0,$
(B.4)

where $\sqrt{g}=\partial\textbf{x}/\partial
b\cdot\left(\partial\textbf{x}/\partial\widetilde{\theta}\times\partial\textbf{x}/\partial\widetilde{\phi}\right)$
By the definition of our extension, the first term will vanish. In the limit
that $\Delta b\rightarrow 0$, the divergence-free condition is expressed as,

$\displaystyle\nabla_{\Gamma}\cdot\textbf{J}_{C}\equiv\frac{1}{\sqrt{g}}\left(\frac{\partial\left(\sqrt{g}J^{\theta}\right)}{\partial\theta}+\frac{\partial\left(\sqrt{g}J^{\phi}\right)}{\partial\phi}\right)=0,$
(B.5)

where we have expressed the current in the contravariant basis as
$\textbf{J}_{C}=J^{\theta}\partial\textbf{x}/\partial\theta+J^{\phi}\partial\textbf{x}/\partial\phi$
and $\nabla_{\Gamma}\cdot$ is the surface divergence (Appendix 3 in [229]).
For a continuous current density, Ampere’s law (1.3b) implies that
$\nabla\cdot\textbf{J}=0$. Thus the equivalent condition for a current
supported on a surface is $\nabla_{\Gamma}\cdot\textbf{J}_{C}=0$ [11]. The
surface divergence of a vector field tangent to a surface $\Gamma$
($\textbf{A}\cdot\hat{\textbf{n}}=0$ on $\Gamma$) defined in terms of a
general continuous extension, $\widetilde{\textbf{A}}$ in a neighborhood of
$\Gamma$ is,

$\displaystyle\nabla_{\Gamma}\cdot\textbf{A}\equiv\left(\nabla\cdot\widetilde{\textbf{A}}\right)\big{\rvert}_{\Gamma}-\hat{\textbf{n}}\cdot\left(\nabla\widetilde{\textbf{A}}\right)\big{\rvert}_{\Gamma}\cdot\hat{\textbf{n}}.$
(B.6)

In (B.2), we have defined our extension such that $\nabla
b\cdot\left(\nabla\widetilde{\textbf{J}}_{C}\right)=0$ such that the second
term in the above expression vanishes.

Given (B.5), we can write,

$\displaystyle J^{\theta}$
$\displaystyle=-\frac{1}{\sqrt{g}}\frac{\partial\Phi(\theta,\phi)}{\partial\phi}$
(B.7a) $\displaystyle J^{\phi}$
$\displaystyle=\frac{1}{\sqrt{g}}\frac{\partial\Phi(\theta,\phi)}{\partial\theta},$
(B.7b) where, $\displaystyle\Phi=\int d\theta\,\sqrt{g}J^{\phi}.$ (B.7c)

In other words,

$\displaystyle\textbf{J}_{C}=\hat{\textbf{n}}\times\nabla\Phi.$ (B.8)

## Appendix C: Adjoint derivative at fixed $J_{\text{max}}$

We enforce $J_{\text{max}}=$ constant in the REGCOIL solve in order to obtain
the regularization parameter $\lambda$ by requiring that the following
constraint be satisfied within a given tolerance,

$\displaystyle
G\left(\Omega,\overrightarrow{\bm{\Phi}}(\Omega,\lambda)\right)=J_{\text{max}}\left(\Omega,\overrightarrow{\bm{\Phi}}(\Omega,\lambda)\right)-J^{\text{target}}_{\text{max}}=0.$
(C.1)

Here $J^{\text{target}}_{\text{max}}$ is the target maximum current density
and $\overrightarrow{\bm{\Phi}}$ is chosen to satisfy the forward equation
(3.8),

$\displaystyle\overrightarrow{\textbf{F}}\left(\Omega,\overrightarrow{\bm{\Phi}},\lambda\right)=\overleftrightarrow{\textbf{A}}(\Omega,\lambda)\overrightarrow{\bm{\Phi}}-\overrightarrow{\textbf{b}}(\Omega,\lambda)=0.$
(C.2)

A log-sum-exponent function is used to approximate the maximum function,
similar to that used to approximate $d_{\text{coil-plasma}}$ (3.24),

$\displaystyle J_{\text{max}}\approx
J_{\text{max},\,\text{lse}}=\frac{1}{p}\log\left(\frac{\int_{S_{C}}d^{2}x\,\exp\left(pJ\right)}{A_{\text{coil}}}\right).$
(C.3)

We compute the total differential of $\overrightarrow{\textbf{F}}$,

$d\overrightarrow{\textbf{F}}(\Omega,\overrightarrow{\bm{\Phi}},\lambda)=\sum_{m,n}\left(\frac{\partial\overleftrightarrow{\textbf{A}}(\Omega,\lambda)}{\partial\Omega_{m,n}}\overrightarrow{\bm{\Phi}}-\frac{\partial\overrightarrow{\textbf{b}}(\Omega,\lambda)}{\partial\Omega_{m,n}}\right)d\Omega_{m,n}+\overleftrightarrow{\textbf{A}}d\overrightarrow{\bm{\Phi}}\\\
+\left(\overleftrightarrow{\textbf{A}}^{K}\overrightarrow{\bm{\Phi}}-\overrightarrow{\textbf{b}}^{K}\right)d\lambda=0.$
(C.4)

Here
$\overleftrightarrow{\textbf{A}}^{K}=\partial\overleftrightarrow{\textbf{A}}/\partial\lambda$
and
$\overrightarrow{\textbf{b}}^{K}=\partial\overrightarrow{\textbf{b}}/\partial\lambda$.
We left multiply by $\overleftrightarrow{\textbf{A}}^{-1}$ and solve for
$d\overrightarrow{\bm{\Phi}}$ such that
$d\overrightarrow{\textbf{F}}(\Omega,\overrightarrow{\bm{\Phi}},\lambda)=0$,

$\displaystyle
d\overrightarrow{\bm{\Phi}}=-\sum_{m,n}\overleftrightarrow{\textbf{A}}^{-1}\left(\frac{\partial\overleftrightarrow{\textbf{A}}(\Omega,\lambda)}{\partial\Omega_{m,n}}\overrightarrow{\bm{\Phi}}-\frac{\partial\overrightarrow{\textbf{b}}(\Omega,\lambda)}{\partial\Omega_{m,n}}\right)d\Omega_{m,n}-\overleftrightarrow{\textbf{A}}^{-1}\left(\overleftrightarrow{\textbf{A}}^{K}\overrightarrow{\bm{\Phi}}-\overrightarrow{\textbf{b}}^{K}\right)d\lambda.$
(C.5)

We also compute the total differential of $G$,

$\displaystyle dG(\Omega,\overrightarrow{\bm{\Phi}})=\sum_{m,n}\frac{\partial
G(\Omega,\overrightarrow{\bm{\Phi}})}{\partial\Omega_{m,n}}d\Omega_{m,n}+\frac{\partial
G(\Omega,\overrightarrow{\bm{\Phi}})}{\partial\overrightarrow{\bm{\Phi}}}\cdot
d\overrightarrow{\bm{\Phi}}=0.$ (C.6)

Using the form for $d\overrightarrow{\bm{\Phi}}$ (C.5), we compute $d\lambda$
in terms of $d\Omega_{m,n}$,

$d\lambda=\left(\frac{\partial
G(\Omega,\overrightarrow{\bm{\Phi}})}{\partial\overrightarrow{\bm{\Phi}}}\cdot\left[\overleftrightarrow{\textbf{A}}^{-1}\left(\overleftrightarrow{\textbf{A}}^{K}\overrightarrow{\bm{\Phi}}-\overrightarrow{\textbf{b}}^{K}\right)\right]\right)^{-1}\\\
\times\sum_{m,n}\left(\frac{\partial
G(\Omega,\overrightarrow{\bm{\Phi}})}{\partial\Omega_{m,n}}-\frac{\partial
G(\Omega,\overrightarrow{\bm{\Phi}})}{\partial\overrightarrow{\bm{\Phi}}}\cdot\left[\overleftrightarrow{\textbf{A}}^{-1}\left(\frac{\partial\overleftrightarrow{\textbf{A}}(\Omega,\lambda)}{\partial\Omega_{m,n}}\overrightarrow{\bm{\Phi}}-\frac{\partial\overrightarrow{\textbf{b}}(\Omega,\lambda)}{\partial\Omega_{m,n}}\right)\right]\right)d\Omega_{m,n}.$
(C.7)

Using (C.5) and (C.7), the derivative of $\overrightarrow{\bm{\Phi}}$ with
respect to $\Omega_{m,n}$ subject to equations (C.1) and (C.2) is given by the
following expression,

$\frac{\partial\overrightarrow{\bm{\Phi}}(\Omega,\lambda(\Omega))}{\partial\Omega_{m,n}}=-\overleftrightarrow{\textbf{A}}^{-1}\left(\frac{\partial\overleftrightarrow{\textbf{A}}(\Omega,\lambda)}{\partial\Omega_{m,n}}\overrightarrow{\bm{\Phi}}-\frac{\partial\overrightarrow{\textbf{b}}(\Omega,\lambda)}{\partial\Omega_{m,n}}\right)-\frac{\overleftrightarrow{\textbf{A}}^{-1}\left(\overleftrightarrow{\textbf{A}}^{K}\overrightarrow{\bm{\Phi}}-\overrightarrow{\textbf{b}}^{K}\right)}{\frac{\partial
G(\Omega,\overrightarrow{\bm{\Phi}})}{\partial\overrightarrow{\bm{\Phi}}}\cdot\left[\overleftrightarrow{\textbf{A}}^{-1}\left(\overleftrightarrow{\textbf{A}}^{K}\overrightarrow{\bm{\Phi}}-\overrightarrow{\textbf{b}}^{K}\right)\right]}\\\
\times\left(\frac{\partial
G(\Omega,\overrightarrow{\bm{\Phi}})}{\partial\Omega_{m,n}}-\frac{\partial
G(\Omega,\overrightarrow{\bm{\Phi}})}{\partial\overrightarrow{\bm{\Phi}}}\cdot\left[\overleftrightarrow{\textbf{A}}^{-1}\left(\frac{\partial\overleftrightarrow{\textbf{A}}(\Omega,\lambda)}{\partial\Omega_{m,n}}\overrightarrow{\bm{\Phi}}-\frac{\partial\overrightarrow{\textbf{b}}(\Omega,\lambda)}{\partial\Omega_{m,n}}\right)\right]\right).$
(C.8)

Here $\overrightarrow{\bm{\Phi}}$ is understood to be a function of $\Omega$
and $\lambda$ through (C.2) and $\lambda$ is understood to be a function of
$\Omega$ through (C.1). We use the adjoint method to avoid solving a linear
system involving the operator $\overleftrightarrow{\textbf{A}}$ for each
$\Omega_{m,n}$,

$\frac{\partial\overrightarrow{\bm{\Phi}}(\Omega,\lambda(\Omega))}{\partial\Omega_{m,n}}=-\overleftrightarrow{\textbf{A}}^{-1}\left(\frac{\partial\overleftrightarrow{\textbf{A}}(\Omega,\lambda)}{\partial\Omega_{m,n}}\overrightarrow{\bm{\Phi}}-\frac{\partial\overrightarrow{\textbf{b}}(\Omega,\lambda)}{\partial\Omega_{m,n}}\right)-\frac{\overleftrightarrow{\textbf{A}}^{-1}\left(\overleftrightarrow{\textbf{A}}^{K}\overrightarrow{\bm{\Phi}}-\overrightarrow{\textbf{b}}^{K}\right)}{\frac{\partial
G(\Omega,\overrightarrow{\bm{\Phi}})}{\partial\overrightarrow{\bm{\Phi}}}\cdot\left[\overleftrightarrow{\textbf{A}}^{-1}\left(\overleftrightarrow{\textbf{A}}^{K}\overrightarrow{\bm{\Phi}}-\overrightarrow{\textbf{b}}^{K}\right)\right]}\\\
\times\left(\frac{\partial
G(\Omega,\overrightarrow{\bm{\Phi}})}{\partial\Omega_{m,n}}-\left[\left(\overleftrightarrow{\textbf{A}}^{T}\right)^{-1}\frac{\partial
G(\Omega,\overrightarrow{\bm{\Phi}})}{\partial\overrightarrow{\bm{\Phi}}}\right]\cdot\left(\frac{\partial\overleftrightarrow{\textbf{A}}(\Omega,\lambda)}{\partial\Omega_{m,n}}\overrightarrow{\bm{\Phi}}-\frac{\partial\overrightarrow{\textbf{b}}(\Omega,\lambda)}{\partial\Omega_{m,n}}\right)\right).$
(C.9)

We introduce a new adjoint vector $\overrightarrow{\widetilde{\textbf{q}}}$,
defined to be the solution of,

$\displaystyle\overleftrightarrow{\textbf{A}}^{T}\overrightarrow{\widetilde{\textbf{q}}}=\frac{\partial
G(\Omega,\overrightarrow{\bm{\Phi}})}{\partial\overrightarrow{\bm{\Phi}}}.$
(C.10)

Equation (C.9) is then used to compute the derivatives of $\chi^{2}_{B}$ with
respect to $\Omega_{m,n}$,

$\displaystyle\frac{\partial\chi^{2}_{B}\left(\Omega,\overrightarrow{\bm{\Phi}}(\Omega,\lambda(\Omega))\right)}{\partial\Omega_{m,n}}=\frac{\partial\chi^{2}_{B}(\Omega,\overrightarrow{\bm{\Phi}})}{\partial\Omega_{m,n}}+\frac{\partial\chi^{2}_{B}(\Omega,\overrightarrow{\bm{\Phi}})}{\partial\overrightarrow{\bm{\Phi}}}\cdot\frac{\partial\overrightarrow{\bm{\Phi}}(\Omega,\lambda(\Omega))}{\partial\Omega_{m,n}}.$
(C.11)

This result can be written in terms of both adjoint variables,
$\overrightarrow{\textbf{q}}$ and $\overrightarrow{\widetilde{\textbf{q}}}$,

$\frac{\partial\chi^{2}_{B}\left(\Omega,\overrightarrow{\bm{\Phi}}(\Omega,\lambda(\Omega))\right)}{\partial\Omega_{m,n}}=\frac{\partial\chi^{2}_{B}(\Omega,\overrightarrow{\bm{\Phi}})}{\partial\Omega_{m,n}}-\overrightarrow{\textbf{q}}\cdot\left(\frac{\partial\overleftrightarrow{\textbf{A}}(\Omega,\lambda)}{\partial\Omega_{m,n}}\overrightarrow{\bm{\Phi}}-\frac{\partial\overrightarrow{\textbf{b}}(\Omega,\lambda)}{\partial\Omega_{m,n}}\right)\\\
-\frac{\overrightarrow{\textbf{q}}\cdot\left(\overleftrightarrow{\textbf{A}}^{K}\overrightarrow{\bm{\Phi}}-\overrightarrow{\textbf{b}}^{K}\right)}{\overrightarrow{\widetilde{\textbf{q}}}\cdot\left(\overleftrightarrow{\textbf{A}}^{K}\overrightarrow{\bm{\Phi}}-\overrightarrow{\textbf{b}}^{K}\right)}\left(\frac{\partial
G(\Omega,\overrightarrow{\bm{\Phi}})}{\partial\Omega_{m,n}}-\overrightarrow{\widetilde{\textbf{q}}}\cdot\left(\frac{\partial\overleftrightarrow{\textbf{A}}(\Omega,\lambda)}{\partial\Omega_{m,n}}\overrightarrow{\bm{\Phi}}-\frac{\partial\overrightarrow{\textbf{b}}(\Omega,\lambda)}{\partial\Omega_{m,n}}\right)\right).$
(C.12)

The same method is used to compute derivatives of
$\mathinner{\\!\left\lVert\textbf{J}\right\rVert}_{2}$. So, to obtain the
derivatives at fixed $J_{\text{max}}$, we compute a solution to the two
adjoint equations, (3.22) and (C.10), in addition to the forward equation,
(3.8).

## Appendix D: Trajectory models

In the SFINCS coordinate system, the DKE can be written in the following way,

$\dot{\textbf{x}}\cdot\nabla f_{1s}+\dot{X}_{s}\frac{\partial f_{1s}}{\partial
X_{s}}+\dot{\xi}_{s}\frac{\partial
f_{1s}}{\partial\xi_{s}}-C_{s}(f_{1s})=-\left(\textbf{v}_{\text{m}s}\cdot\nabla\psi\right)\frac{\partial
f_{Ms}}{\partial\psi}.$ (D.1)

To obtain the trajectory coefficients ($\dot{\textbf{x}}$, $\dot{X}_{s}$, and
$\dot{\xi}_{s}$) several approximations are made. For example, any terms that
require radial coupling ($\psi$ derivatives of $f_{1s}$) cannot be retained,
as this would necessitate solving a five-dimensional system.

Under the full trajectory model, the trajectory coefficients are chosen such
that $\mu$ conservation is maintained as radial coupling is dropped,

$\displaystyle\dot{\textbf{x}}$
$\displaystyle=v_{||}\hat{\textbf{b}}+\frac{\Phi^{\prime}(\psi)}{B^{2}}\textbf{B}\times\nabla\psi$
(D.2a) $\displaystyle\dot{X}_{s}$
$\displaystyle=-\left(\textbf{v}_{\text{m}s}\cdot\nabla\psi\right)\frac{q_{s}}{2T_{s}X_{s}}\Phi^{\prime}(\psi)$
(D.2b) $\displaystyle\dot{\xi}_{s}$
$\displaystyle=-\frac{1-\xi_{s}^{2}}{2B\xi_{s}}v_{||}\hat{\textbf{b}}\cdot\nabla
B+\xi_{s}(1-\xi_{s}^{2})\frac{1}{2B^{3}}\Phi^{\prime}(\psi)\textbf{B}\times\nabla\psi\cdot\nabla
B.$ (D.2c)

Under the DKES trajectory model, the $\textbf{E}\times\textbf{B}$ velocity is
taken to be divergenceless,

$\displaystyle\textbf{v}_{E}^{\text{DKES}}=\frac{\textbf{B}\times\nabla\Phi}{\langle
B^{2}\rangle_{\psi}},$ (D.3)

where the flux surface average of a quantity is (4.8). Under the DKES
trajectory model, the trajectory coefficients are taken to be,

$\displaystyle\dot{\textbf{x}}$
$\displaystyle=v_{||}\hat{\textbf{b}}+\frac{1}{\langle
B^{2}\rangle_{\psi}}\Phi^{\prime}(\psi)\textbf{B}\times\nabla\psi$ (D.4a)
$\displaystyle\dot{X}_{s}$ $\displaystyle=0$ (D.4b)
$\displaystyle\dot{\xi}_{s}$
$\displaystyle=-\frac{1-\xi_{s}^{2}}{2B\xi_{s}}v_{||}\hat{\textbf{b}}\cdot\nabla
B.$ (D.4c)

These effective trajectories are adopted in the widely-used DKES code [113,
230].

## Appendix E: Adjoint collision operator

We want to find an adjoint collision operator, $C_{s}^{\dagger}$, that
satisfies the following relation,

$\displaystyle\left\langle\int
d^{3}v\,\frac{g_{1s}C_{s}(f_{1s})}{f_{Ms}}\right\rangle_{\psi}=\left\langle\int
d^{3}v\,\frac{f_{1s}C_{s}^{\dagger}(g_{1s})}{f_{Ms}}\right\rangle_{\psi}.$
(E.1)

The linearized Fokker-Planck collision operator can be written as,

$\displaystyle
C_{s}(f_{1s})=\sum_{s^{\prime}}C_{ss^{\prime}}^{L}(f_{1s},f_{1s^{\prime}})=\sum_{s^{\prime}}C_{ss^{\prime}}(f_{1s},f_{Ms^{\prime}})+C_{ss^{\prime}}(f_{Ms},f_{1s^{\prime}}),$
(E.2)

where $s^{\prime}$ sums over species. The first term on the right hand side of
(E.2) is referred to as the test-particle collision operator,
$C^{T}_{ss^{\prime}}(f_{1s})=C_{ss^{\prime}}(f_{1s},f_{Ms^{\prime}})$, and the
second the field-particle collision operator,
$C^{F}_{ss^{\prime}}(f_{1s^{\prime}})=C_{ss^{\prime}}(f_{Ms},f_{1s^{\prime}})$.
The test and field terms satisfy the following relations [198, 221],

$\displaystyle\int
d^{3}v\,\frac{g_{1s}C_{ss^{\prime}}(f_{1s},f_{Ms^{\prime}})}{f_{Ms}}$
$\displaystyle=\int
d^{3}v\,\frac{f_{1s}C_{ss^{\prime}}(g_{1s},f_{Ms^{\prime}})}{f_{Ms}}$ (E.3a)
$\displaystyle\int
d^{3}v\,\frac{g_{1s}C_{ss^{\prime}}(f_{Ms},f_{1s^{\prime}})}{f_{Ms}}$
$\displaystyle=\frac{T_{s^{\prime}}}{T_{s}}\int
d^{3}v\,\frac{f_{1s^{\prime}}C_{s^{\prime}s}(f_{Ms^{\prime}},g_{1s})}{f_{Ms^{\prime}}}.$
(E.3b)

For collisions between species of the same temperature, we see that
$C_{s}(f_{1s})$ is self-adjoint. The adjoint operator with respect to the
inner product (4.32) is thus,

$\displaystyle C_{s}^{\dagger}$
$\displaystyle=C_{s}^{T}+\sum_{s^{\prime}}\frac{f_{Ms}}{f_{Ms^{\prime}}}\frac{T_{s^{\prime}}}{T_{s}}C_{s^{\prime}s}^{F}.$
(E.4)

## Appendix F: Adjoint collisionless trajectories

We want to find an adjoint operator, $\mathbb{L}_{0s}^{\dagger}$, that
satisfies,

$\displaystyle\left\langle\int
d^{3}v\,\frac{g_{1s}\mathbb{L}_{0s}f_{1s}}{f_{Ms}}\right\rangle_{\psi}=\left\langle\int
d^{3}v\,\frac{f_{1s}\mathbb{L}_{0s}^{\dagger}g_{1s}}{f_{Ms}}\right\rangle_{\psi},$
(F.1)

for both trajectory models, where $\mathbb{L}_{0s}$ is defined in (4.10) with
(D.4) for the DKES trajectories model and (D.2) for the full trajectory model.
Throughout we use the velocity space element in SFINCS coordinates,
$d^{3}v=2\pi v_{ts}^{3}X_{s}^{2}d\xi_{s}dX_{s}$.

####  F.0.1 DKES trajectories

The operator under consideration is,

$\displaystyle\mathbb{L}_{0s}=v_{||}\hat{\textbf{b}}\cdot\nabla+\hat{\textbf{v}}_{E}^{\text{DKES}}\cdot\nabla-\frac{1-\xi_{s}^{2}}{2B\xi_{s}}v_{||}\hat{\textbf{b}}\cdot\nabla
B\frac{\partial}{\partial\xi_{s}}.$ (F.2)

Considering the contribution of the streaming term in (F.2) to the left hand
side of (F.1) we obtain,

$\displaystyle\left\langle\int
d^{3}v\,\frac{g_{1s}v_{||}\hat{\textbf{b}}\cdot\nabla
f_{1s}}{f_{Ms}}\right\rangle_{\psi}=-\left\langle\int
d^{3}v\,\frac{f_{1s}v_{||}\textbf{B}\cdot\nabla\left(g_{1s}/B\right)}{f_{Ms}}\right\rangle_{\psi}.$
(F.3)

Here the identity
$\langle\nabla\cdot\textbf{Q}\rangle_{\psi}=1/V^{\prime}(\psi)\partial/\partial\psi\left(V^{\prime}(\psi)\langle\textbf{Q}\cdot\nabla\psi\rangle_{\psi}\right)$
for any vector Q has been used. We next consider the contribution of the
$\textbf{E}\times\textbf{B}$ drift term in (F.2),

$\displaystyle\left\langle\int
d^{3}v\,\frac{g_{1s}\textbf{v}_{E}^{\text{DKES}}\cdot\nabla
f_{1s}}{f_{Ms}}\right\rangle_{\psi}$ $\displaystyle=-\left\langle\int
d^{3}v\,\frac{f_{1s}\textbf{v}_{E}^{\text{DKES}}\cdot\nabla
g_{1s}}{f_{Ms}}\right\rangle_{\psi}.$ (F.4)

Here we have used the identity,

$\displaystyle\left\langle\textbf{B}\times\nabla\psi\cdot\nabla
w\right\rangle_{\psi}=0,$ (F.5)

for any $w$. We consider the contribution of the mirror-force term in (F.2),

$\displaystyle\left\langle\int
d^{3}v\,\frac{g_{1s}\dot{\xi}_{s}}{f_{Ms}}\frac{\partial
f_{1s}}{\partial\xi_{s}}\right\rangle_{\psi}=-\left\langle\int
d^{3}v\,\frac{f_{1s}\dot{\xi}_{s}}{f_{Ms}}\frac{\partial
g_{1s}}{\partial\xi_{s}}\right\rangle_{\psi}-\left\langle\int
d^{3}v\,\frac{v_{||}}{B}\hat{\textbf{b}}\cdot\nabla
B\frac{g_{1s}f_{1s}}{f_{Ms}}\right\rangle_{\psi}.$ (F.6)

Combining (F.3-F.6), we obtain

$\displaystyle\left\langle\int
d^{3}v\,\frac{g_{1s}\mathbb{L}_{0s}f_{1s}}{f_{Ms}}\right\rangle_{\psi}=-\left\langle\int
d^{3}v\,\frac{f_{1s}\mathbb{L}_{0s}g_{1s}}{f_{Ms}}\right\rangle_{\psi}.$ (F.7)

Therefore, in the DKES trajectory model we obtain (4.51).

####  F.0.2 Full trajectories

The operator under consideration for the full model is,

$\mathbb{L}_{0s}=v_{||}\hat{\textbf{b}}\cdot\nabla+\textbf{v}_{E}\cdot\nabla+\frac{(1+\xi_{s}^{2})X_{s}}{2B}\textbf{v}_{E}\cdot\nabla
B\frac{\partial}{\partial X_{s}}\\\
-\frac{1-\xi_{s}^{2}}{2B\xi_{s}}v_{||}\hat{\textbf{b}}\cdot\nabla
B\frac{\partial}{\partial\xi_{s}}+\frac{\xi_{s}(1-\xi_{s}^{2})}{2B}\textbf{v}_{E}\cdot\nabla
B\frac{\partial}{\partial\xi_{s}}.$ (F.8)

The contribution to (F.1) from the streaming term in (F.8) is identical to
that in the case of the DKES trajectory model, (F.3). We next consider the
contribution from the $\textbf{E}\times\textbf{B}$ drift term in (F.8),

$\displaystyle\left\langle\int d^{3}v\,\frac{g_{1s}\textbf{v}_{E}\cdot\nabla
f_{1s}}{f_{Ms}}\right\rangle_{\psi}$ $\displaystyle=-\left\langle\int
d^{3}v\,\frac{f_{1s}B^{2}\textbf{v}_{E}\cdot\nabla\left(g_{1s}/B^{2}\right)}{f_{Ms}}\right\rangle_{\psi},$
(F.9)

again using (F.5). The contribution from the $\dot{X}_{s}$ term in (F.8) is,

$\left\langle\int d^{3}v\,\frac{g_{1s}\dot{X}_{s}}{f_{Ms}}\frac{\partial
f_{1s}}{\partial X_{s}}\right\rangle_{\psi}=-\left\langle\int
d^{3}v\,\frac{f_{1s}\dot{X}_{s}}{f_{Ms}}\frac{\partial g_{1s}}{\partial
X_{s}}\right\rangle_{\psi}\\\ -\left\langle\int
d^{3}v\,(3+2X_{s}^{2})(1+\xi_{s}^{2})\frac{g_{1s}f_{1s}}{2f_{Ms}B}\textbf{v}_{E}\cdot\nabla
B\right\rangle_{\psi}.$ (F.10)

The contribution from the mirror term in (F.8) is the same as in the case of
the DKES trajectories model (F.6). We consider the contribution from the final
term in (F.8),

$\left\langle\int
d^{3}v\,\frac{g_{1s}\xi_{s}(1-\xi_{s}^{2})\textbf{v}_{E}\cdot\nabla
B}{2Bf_{Ms}}\frac{\partial
f_{1s}}{\partial\xi_{s}}\right\rangle_{\psi}=-\left\langle\int
d^{3}v\,\frac{f_{1s}\xi_{s}(1-\xi_{s}^{2})\textbf{v}_{E}\cdot\nabla
B}{2Bf_{Ms}}\frac{\partial g_{1s}}{\partial\xi_{s}}\right\rangle_{\psi}\\\
-\left\langle\int d^{3}v\,(1-3\xi_{s}^{2})\textbf{v}_{E}\cdot\nabla
B\frac{f_{1s}g_{1s}}{2Bf_{M}}\right\rangle_{\psi}.$ (F.11)

Combining (F.3), (F.9), (F.10), (F.6), and (F.11), we obtain

$\left\langle\int
d^{3}v\,\frac{g_{1s}\mathbb{L}_{0s}f_{1s}}{f_{Ms}}\right\rangle_{\psi}=-\left\langle\int
d^{3}v\,\frac{f_{1s}\mathbb{L}_{0s}g_{1s}}{f_{Ms}}\right\rangle_{\psi}\\\
+\Phi^{\prime}(\psi)\frac{q_{s}}{T_{s}}\left\langle\int
d^{3}v\,\left(\textbf{v}_{\text{m}s}\cdot\nabla\psi\right)\frac{f_{1s}g_{1s}}{f_{Ms}}\right\rangle_{\psi}.$
(F.12)

Therefore, under the full trajectory model we obtain (4.52).

## Appendix G: Symmetry of the sensitivity function

In this Appendix we discuss several symmetry properties of the local
sensitivity function, $S_{\mathcal{R}}$, defined through (4.65). The arguments
that follow are similar to those in Appendix C of [138]. Throughout we will
assume that $B$ is stellarator symmetric and $N_{P}$ symmetric. We will show
that this implies $N_{P}$ symmetry of $S_{\mathcal{R}}$. In the limit that
$E_{r}\rightarrow 0$, then $S_{\mathcal{R}}$ also has stellarator symmetry.

####  G.0.1 Symmetry of $S_{\mathcal{R}}$ implied by Fourier derivatives

First we would like to show that $S_{\mathcal{R}}$ is stellarator symmetric if
and only if $\partial\mathcal{R}/\partial B_{m,n}^{s}=0$ for all $m$ and $n$,
where we express $B$ in a Fourier series,

$\displaystyle
B=\sum_{m,n}B_{m,n}^{c}\cos(m\vartheta_{B}-n\varphi_{B})+B_{m,n}^{s}\sin(m\vartheta_{B}-n\varphi_{B}).$
(G.1)

The perturbation, $\delta B$, is decomposed similarly. We begin with the “if”
portion of the argument. From (4.65) we have,

$\displaystyle\frac{\partial\mathcal{R}}{\partial B_{m,n}^{s}}$
$\displaystyle=V^{\prime}(\psi)^{-1}\int_{0}^{2\pi}d\vartheta_{B}\int_{0}^{2\pi}d\varphi_{B}\,\sqrt{g}S_{\mathcal{R}}\sin(m\vartheta_{B}-n\varphi_{B}).$
(G.2)

Suppose $\partial\mathcal{R}/\partial B_{m,n}^{s}=0$ for all $m$ and $n$. The
quantity $(\sqrt{g}S_{\mathcal{R}})$ can be represented as a Fourier series,

$\displaystyle\left(\sqrt{g}S_{\mathcal{R}}\right)=\sum_{m,n}A_{m,n}^{c}\cos(m\vartheta_{B}-n\varphi_{B})+A_{m,n}^{s}\sin(m\vartheta_{B}-n\varphi_{B}).$
(G.3)

From (G.2), we see that $A_{m,n}^{s}=0$ for all $m$ and $m$. Thus the quantity
$(\sqrt{g}S_{\mathcal{R}})$ must be even under the transformation
$(\vartheta_{B},\varphi_{B})\rightarrow(-\vartheta_{B},-\varphi_{B})$. We now
note that $\sqrt{g}$ must be even from (4.61) under the assumption that $B$ is
stellarator symmetric. Therefore $S_{\mathcal{R}}$ must be stellarator
symmetric, assuming that $\sqrt{g}$ does not vanish anywhere, which must be
the case for any well-defined coordinate transformation.

We continue with the “only if” portion of the argument. Suppose
$S_{\mathcal{R}}$ is stellarator symmetric. As $\sqrt{g}$ is also stellarator
symmetric, $(\sqrt{g}S_{\mathcal{R}})$ can be expressed in a Fourier series as
(G.3) with $A_{m,n}^{s}=0$ for all $m$ and $n$. Thus from (G.2)
$\partial\mathcal{R}/\partial B_{m,n}^{s}=0$ for all $m$ and $n$.

We next show that if $B$ is $N_{P}$ symmetric, then $S_{\mathcal{R}}$ is
$N_{P}$ symmetric if and only if $\partial\mathcal{R}/\partial B_{m,n}^{c}=0$
for all $n$ that are not integer multiples of $N_{P}$. We begin with the “if”
portion of the argument. From (4.65),

$\displaystyle\frac{\partial\mathcal{R}}{\partial
B_{m,n}^{c}}=V^{\prime}(\psi)^{-1}\int_{0}^{2\pi}d\vartheta_{B}\int_{0}^{2\pi}d\varphi_{B}\,\sqrt{g}S_{\mathcal{R}}\cos(m\vartheta_{B}-n\varphi_{B}).$
(G.4)

Suppose $\partial\mathcal{R}/\partial B_{m,n}^{c}=0$ for all $n$ which are not
integer multiples of $N_{P}$. Here $(\sqrt{g}S_{\mathcal{R}})$ can be
expressed in a Fourier series as (G.3) with $A_{m,n}^{s}=0$ for all $m$ and
$n$. Inserting the Fourier series into (G.4), we find that $A_{m,n}^{c}=0$ for
all $n$ that are not integer multiples of $N_{P}$. Thus
$(\sqrt{g}S_{\mathcal{R}})$ must be $N_{P}$ symmetric. As $\sqrt{g}$ must be
$N_{P}$ symmetric, this implies $S_{\mathcal{R}}$ possesses the same symmetry.

Next we consider the “only if” portion of the argument. Suppose that
$S_{\mathcal{R}}$ is $N_{P}$ symmetric. As $\sqrt{g}$ is also $N_{P}$
symmetric, then $(\sqrt{g}S_{\mathcal{R}})$ can be expressed in a Fourier
series as (G.3) where the sum includes $n$ that are integer multiples of
$N_{P}$. Inserting the Fourier series into (G.4), we find that
$\partial\mathcal{R}/\partial B_{m,n}^{c}=0$ for all $n$ that are not integer
multiples of $N_{P}$.

####  G.0.2 Symmetry of Fourier derivatives

To continue, we need to show that $\partial\mathcal{R}/\partial B_{m,n}^{s}=0$
for all $m$ and $n$ and $\partial\mathcal{R}/\partial B_{m,n}^{c}=0$ for all
$n$ which are not integer multiples of $N_{P}$. We begin with the $N_{P}$
symmetry argument. We consider the symmetry of $f_{1s}$ implied by (D.1).
Under the transformation $\varphi_{B}\rightarrow\varphi_{B}+2\pi/N_{P}$, we
find that each of the trajectory coefficients remain unchanged, as well as the
source term and collision operator. Therefore we can conclude that $f_{1s}$ is
$N_{P}$ symmetric. We can also note that each of the $\widetilde{\mathcal{R}}$
vectors are $N_{P}$ symmetric, as well as $\sqrt{g}$. We consider the
integrand that appears in the flux surface average in (4.34),

$\displaystyle D_{s}(\vartheta_{B},\varphi_{B})=\int
d^{3}v\,\frac{f_{1s}\widetilde{\mathcal{R}}_{s}^{f}\sqrt{g}}{f_{Ms}}.$ (G.5)

Here the superscript and subscript on $\widetilde{\mathcal{R}}$ denotes that
we consider the unknowns corresponding to the distribution function of species
$s$. We note that
$D_{s}(\vartheta_{B},\varphi_{B}+2\pi/N_{P})=D_{s}(\vartheta_{B},\varphi_{B})$.
The quantity $\mathcal{R}$ can be expressed in terms of $D_{s}$ as follows,

$\displaystyle\mathcal{R}=\sum_{s}V^{\prime}(\psi)^{-1}\int_{0}^{2\pi}d\varphi_{B}\int_{0}^{2\pi}d\vartheta_{B}\,D_{s}.$
(G.6)

Next we consider the functional derivative of $\mathcal{R}$ with respect to
$B$, defined as in (4.64). The derivative with respect to $B_{m,n}^{c}$ can be
thus defined as,

$\displaystyle\frac{\partial\mathcal{R}}{\partial
B_{m,n}^{c}}=V^{\prime}(\psi)^{-1}\int_{0}^{2\pi}d\varphi_{B}\int_{0}^{2\pi}d\vartheta_{B}\,\left(\sum_{s}\frac{\delta
D_{s}}{\delta B}-\mathcal{R}\frac{\delta\sqrt{g}}{\delta
B}\right)\cos(m\vartheta_{B}-n\varphi_{B}).$ (G.7)

As the functional derivative maintains the $N_{P}$ symmetry of $D_{s}$ and
$\sqrt{g}$, the quantity in parenthesis in (G.7) can be expressed in a Fourier
series containing only $n$ that are integer multiples of $N_{P}$. Thus we see
that the quantity $\partial\mathcal{R}/\partial B_{m,n}^{c}=0$ for all $n$
that are not integer multiples of $N_{P}$.

Next we consider a similar argument for stellarator symmetry. We begin by
considering the symmetry of $f_{1s}$ implied by (D.1) in the case $E_{r}=0$.
Under the transformation
$(\vartheta_{B},\varphi_{B},v_{||})\rightarrow(-\vartheta_{B},-\varphi_{B},-v_{||})$,
we see that both the collisionless trajectory operator and the collision
operator maintain the parity of $f_{1s}$, while the source term is odd.
Therefore, $f_{1s}$ must be odd under this transformation. In this case, we
can write $f_{1s}$ as,

$\displaystyle
f_{1s}=f^{-}_{a,s}(X_{s},\xi_{s})f^{+}_{b,s}(\vartheta_{B},\varphi_{B})+f^{+}_{a,s}(X_{s},\xi_{s})f^{-}_{b,s}(\vartheta_{B},\varphi_{B}),$
(G.8)

where $f^{-}_{a,s}(X_{s},-\xi_{s})=-f^{-}_{a,s}(X_{s},\xi_{s})$,
$f^{+}_{a,s}(X_{s},-\xi_{s})=f^{+}_{a,s}(X_{s},-\xi_{s})$, and analogous
expressions for $f_{b,s}^{+}$ and $f_{b,s}^{-}$.

We next note that each of the $\widetilde{\mathcal{R}}^{f}_{s}$ are odd under
the transformation
$(\vartheta_{B},\varphi_{B},v_{||})\rightarrow(-\vartheta_{B},-\varphi_{B},-v_{||})$.
As $\sqrt{g}$ is even, then we can express
$\widetilde{\mathcal{R}}_{s}^{f}\sqrt{g}$ in a similar way to (G.8),

$\displaystyle\widetilde{\mathcal{R}}_{s}^{f}\sqrt{g}=B_{a,s}^{-}(X_{s},\xi_{s})B_{b,s}^{+}(\vartheta_{B},\varphi_{B})+B_{a,s}^{+}(X_{s},\xi_{s})B_{b,s}^{-}(\vartheta_{B},\varphi_{B}).$
(G.9)

The integrand that appears in the flux surface average becomes,

$D_{s}=\int
d^{3}v\,f_{Ms}^{-1}\bigg{(}f_{a,s}^{-}(X_{s},\xi_{s})B_{a,s}^{-}(X_{s},\xi_{s})f_{b,s}^{+}(\vartheta_{B},\varphi_{B})B_{b,s}^{+}(\vartheta_{B},\varphi_{B})\\\
+f_{a,s}^{+}(X_{s},\xi_{s})B_{a,s}^{+}(X_{s},\xi_{s})f_{b,s}^{-}(\vartheta_{B},\varphi_{B})B_{b,s}^{-}(\vartheta_{B},\varphi_{B})\bigg{)}.$
(G.10)

We see that $D_{s}$ is even with respect to the transformation
$(\vartheta_{B},\varphi_{B})\rightarrow(-\vartheta_{B},-\varphi_{B})$. The
quantity $\mathcal{R}$ can be written as in (G.6) and the derivative with
respect to a stellarator asymmetric mode is

$\displaystyle\frac{\partial\mathcal{R}}{\partial
B_{m,n}^{s}}=V^{\prime}(\psi)^{-1}\int_{0}^{2\pi}d\varphi_{B}\int_{0}^{2\pi}d\vartheta_{B}\,\left(\sum_{s}\frac{\delta
D_{s}}{\delta B}-\mathcal{R}\frac{\delta\sqrt{g}}{\delta
B}\right)\sin(m\vartheta_{B}-n\varphi_{B}).$ (G.11)

The functional derivative with respect to $B$ does not change the parity of
$D_{s}$ or $\sqrt{g}$, thus we see that the quantity in parenthesis in the
above equation is even with respect to the transformation
$(\vartheta_{B},\varphi_{B})\rightarrow(-\vartheta_{B},-\varphi_{B})$.
Therefore, $\partial\mathcal{R}/\partial B_{m,n}^{s}=0$ for all $m$ and $n$. A
similar argument cannot be made if $E_{r}\neq 0$, as the inhomogeneous drive
term in (D.1) no longer has definite parity. However, according to the
arguments in [112] the transport coefficients do obey this symmetry property.

## Appendix H: Derivatives at ambipolarity

In this Appendix, we derive an expression for derivatives of moments of the
distribution function at fixed ambipolarity rather than fixed $E_{r}$ by
determining the relationship between geometry parameters, $\Omega$, and
$E_{r}$. We begin by assuming that the continuous adjoint approach outlined in
Section 4.3.1 is used. The approach taken here is analogous to that used in
Appendix C, in which an additional adjoint equation is used to compute
derivatives at a fixed constraint function for optimization of stellarator
coil shapes.

Consider the set of unknowns computed with SFINCS, $F$, which depends on
parameters $\Omega$ and $E_{r}$. The total differential of $F$ satisfies,

$\mathbb{L}dF(\Omega,E_{r})=\left(\frac{\partial\mathbb{S}(\Omega,E_{r})}{\partial
E_{r}}-\frac{\partial\mathbb{L}(\Omega,E_{r})}{\partial
E_{r}}F\right)dE_{r}\\\
+\sum_{i=1}^{N_{\Omega}}\left(\frac{\partial\mathbb{S}(\Omega,E_{r})}{\partial\Omega_{i}}-\frac{\partial\mathbb{L}(\Omega,E_{r})}{\partial\Omega_{i}}F\right)d\Omega_{i},$
(H.1)

which follows from (4.31). Consider $J_{r}(\Omega,F)$, which depends on
$E_{r}$ through $F$. The total differential of $J_{r}$ can be computed,

$\displaystyle
dJ_{r}(\Omega,F(\Omega,E_{r}))=\sum_{i=1}^{N_{\Omega}}\frac{\partial
J_{r}\left(\Omega,F\right)}{\partial\Omega_{i}}d\Omega_{i}+\left\langle\widetilde{J_{r}},dF(\Omega,E_{r})\right\rangle,$
(H.2)

which can be written using (H.1) and the solution to (4.74),

$dJ_{r}(\Omega,F(\Omega,E_{r}))=\left\langle\lambda^{J_{r}},\left(\frac{\partial\mathbb{L}(\Omega,E_{r})}{\partial
E_{r}}F-\frac{\partial\mathbb{S}(\Omega,E_{r})}{\partial
E_{r}}\right)\right\rangle dE_{r}\\\
+\sum_{i=1}^{N_{\Omega}}\left(\frac{\partial
J_{r}(\Omega,F)}{\partial\Omega_{i}}+\left\langle\lambda^{J_{r}},\left(\frac{\partial\mathbb{L}(\Omega,E_{r})}{\partial\Omega_{i}}F-\frac{\partial\mathbb{S}(\Omega,E_{r})}{\partial\Omega_{i}}\right)\right\rangle\right)d\Omega_{i}.$
(H.3)

By enforcing $dJ_{r}(\Omega,F(\Omega,E_{r}))=0$, we obtain the relationship
between $E_{r}$ and $\Omega$ at ambipolarity,

$\frac{\partial
E_{r}(\Omega)}{\partial\Omega_{i}}\bigg{\rvert}_{dJ_{r}=0}=-\left\langle\lambda^{J_{r}},\left(\frac{\partial\mathbb{L}(\Omega,E_{r})}{\partial
E_{r}}F-\frac{\partial\mathbb{S}(\Omega,E_{r})}{\partial
E_{r}}\right)\right\rangle^{-1}\\\ \left(\frac{\partial
J_{r}(\Omega,F)}{\partial\Omega_{i}}+\left\langle\lambda^{J_{r}},\left(\frac{\partial\mathbb{L}(\Omega,E_{r})}{\partial\Omega_{i}}F-\frac{\partial\mathbb{S}(\Omega,E_{r})}{\partial\Omega_{i}}\right)\right\rangle\right).$
(H.4)

Consider a moment of the distribution function,
$\mathcal{R}(\Omega,F(\Omega,E_{r}))$. The derivative with respect to
$\Omega_{i}$ at fixed ambipolarity can thus be computed,

$\displaystyle\frac{\partial\mathcal{R}(\Omega,F(\Omega,E_{r}(\Omega))}{\partial\Omega_{i}}=\frac{\partial\mathcal{R}(\Omega,F)}{\partial\Omega_{i}}+\left\langle\widetilde{\mathcal{R}},\frac{\partial
F(\Omega,E_{r}(\Omega))}{\partial\Omega_{i}}\right\rangle,$ (H.5)

where $E_{r}$ is viewed as a function of $\Omega$ through (H.4). The first
term corresponds to the explicit dependence on $\Omega_{i}$, while the second
contains dependence through $F$. Here $\partial
F(\Omega,E_{r}(\Omega))/\partial\Omega_{i}$ satisfies,

$\mathbb{L}\frac{\partial
F(\Omega,E_{r}(\Omega))}{\partial\Omega_{i}}=\left(\frac{\partial\mathbb{S}(\Omega,E_{r})}{\partial\Omega_{i}}-\frac{\partial\mathbb{L}(\Omega,E_{r})}{\partial\Omega_{i}}F\right)\\\
-\left(\frac{\partial\mathbb{S}(\Omega,E_{r})}{\partial
E_{r}}-\frac{\partial\mathbb{L}(\Omega,E_{r})}{\partial
E_{r}}F\right)\left\langle\lambda^{J_{r}},\left(\frac{\partial\mathbb{L}(\Omega,E_{r})}{\partial
E_{r}}F-\frac{\partial\mathbb{S}(\Omega,E_{r})}{\partial
E_{r}}\right)\right\rangle^{-1}\\\ \times\left(\frac{\partial
J_{r}(\Omega,F)}{\partial\Omega_{i}}+\left\langle\lambda^{J_{r}},\left(\frac{\partial\mathbb{L}(\Omega,E_{r})}{\partial\Omega_{i}}F-\frac{\partial\mathbb{S}(\Omega,E_{r})}{\partial\Omega_{i}}\right)\right\rangle\right),$
(H.6)

from (H.1) using (H.4). Using (H.6) and (4.44), we find

$\frac{\partial\mathcal{R}(\Omega,F(\Omega,E_{r}(\Omega))}{\partial\Omega_{i}}=\frac{\partial\mathcal{R}(\Omega,F)}{\partial\Omega_{i}}+\left\langle\lambda^{\mathcal{R}},\left(\frac{\partial\mathbb{L}(\Omega,E_{r})}{\partial\Omega_{i}}F-\frac{\partial\mathbb{S}(\Omega,E_{r})}{\partial\Omega_{i}}\right)\right\rangle\\\
-\left\langle\lambda^{\mathcal{R}},\left(\frac{\partial\mathbb{L}(\Omega,E_{r})}{\partial
E_{r}}F-\frac{\partial\mathbb{S}(\Omega,E_{r})}{\partial
E_{r}}\right)\right\rangle\\\ \times\frac{\left(\frac{\partial
J_{r}(\Omega,F)}{\partial\Omega_{i}}+\left\langle\lambda^{J_{r}},\left(\frac{\partial\mathbb{L}(\Omega,E_{r})}{\partial\Omega_{i}}F-\frac{\partial\mathbb{S}(\Omega,E_{r})}{\partial\Omega_{i}}\right)\right\rangle\right)}{\left\langle\lambda^{J_{r}},\left(\frac{\partial\mathbb{L}(\Omega,E_{r})}{\partial
E_{r}}F-\frac{\partial\mathbb{S}(\Omega,E_{r})}{\partial
E_{r}}\right)\right\rangle}.$ (H.7)

An analogous expression can be obtained using the discrete approach,

$\frac{\partial\mathcal{R}\left(\Omega,\overrightarrow{\textbf{F}}\left(\Omega,E_{r}(\Omega)\right)\right)}{\partial\Omega_{i}}=\frac{\partial\mathcal{R}\left(\Omega,\overrightarrow{\textbf{F}}\right)}{\partial\Omega_{i}}+\left\langle\overrightarrow{\bm{\lambda}}^{\mathcal{R}},\left(\frac{\partial\overrightarrow{\textbf{S}}(\Omega,E_{r})}{\partial\Omega_{i}}-\frac{\partial\overleftrightarrow{\textbf{L}}(\Omega,E_{r})}{\partial\Omega_{i}}\overrightarrow{\textbf{F}}\right)\right\rangle\\\
-\left\langle\overrightarrow{\bm{\lambda}}^{\mathcal{R}},\left(\frac{\partial\overrightarrow{\textbf{S}}(\Omega,E_{r})}{\partial
E_{r}}-\frac{\partial\overleftrightarrow{\textbf{L}}(\Omega,E_{r})}{\partial
E_{r}}\overrightarrow{\textbf{F}}\right)\right\rangle\\\
\times\frac{\left(\frac{\partial
J_{r}\left(\Omega,\overrightarrow{\textbf{F}}\right)}{\partial\Omega_{i}}+\left\langle\overrightarrow{\bm{\lambda}}^{J_{r}},\left(\frac{\partial\overrightarrow{\textbf{S}}(\Omega,E_{r})}{\partial\Omega_{i}}-\frac{\partial\overleftrightarrow{\textbf{L}}(\Omega,E_{r})}{\partial\Omega_{i}}\overrightarrow{\textbf{F}}\right)\right\rangle\right)}{\left\langle\overrightarrow{\bm{\lambda}}^{J_{r}},\left(\frac{\partial\overrightarrow{\textbf{S}}(\Omega,E_{r})}{\partial
E_{r}}-\frac{\partial\overleftrightarrow{\textbf{L}}(\Omega,E_{r})}{\partial
E_{r}}\overrightarrow{\textbf{F}}\right)\right\rangle},$ (H.8)

where (4.75) has been used.

## Appendix I: Derivation of generalized MHD self-adjointness relation

The quantity $U_{P}=U_{P_{1}}+U_{P_{2}}$ consists of two terms, accounting for
changes to the vector potential due to MHD perturbations,

$\displaystyle
U_{P_{1}}=\int_{V_{P}}d^{3}x\,\left(\delta\textbf{J}_{1}\cdot\bm{\xi}_{2}\times\textbf{B}-\delta\textbf{J}_{2}\cdot\bm{\xi}_{1}\times\textbf{B}\right),$
(I.1)

and changes to the rotational transform,

$\displaystyle
U_{P_{2}}=\int_{V_{P}}d^{3}x\,\left(\delta\chi_{1}(\psi)\delta\textbf{J}_{2}\cdot\nabla\varphi-\delta\chi_{2}(\psi)\delta\textbf{J}_{1}\cdot\nabla\varphi\right).$
(I.2)

The quantity $U_{P_{1}}$ can be expressed by using (5.26) and applying the
divergence theorem to the pressure gradient terms,

$U_{P_{1}}=\int_{V_{P}}d^{3}x\,\bm{\xi}_{2}\cdot\left(\textbf{J}\times\delta\textbf{B}_{1}+\nabla
p\left(\nabla\cdot\bm{\xi}_{1}\right)-\textbf{F}_{1}\right)\\\
-\int_{V_{P}}d^{3}x\,\bm{\xi}_{1}\cdot\left(\textbf{J}\times\delta\textbf{B}_{2}+\nabla
p\left(\nabla\cdot\bm{\xi}_{2}\right)-\textbf{F}_{2}\right).$ (I.3)

We will define
$\delta\widetilde{\textbf{B}}_{1,2}=\nabla\times\left(\bm{\xi}_{1,2}\times\textbf{B}\right)$
such that
$\delta\textbf{B}_{1,2}=\delta\widetilde{\textbf{B}}_{1,2}-\nabla\delta\chi_{1,2}(\psi)\times\nabla\varphi$.
The terms in (I.3) due to $\delta\widetilde{\textbf{B}}_{1,2}$ can be
evaluated using
$\textbf{J}=J_{||}\hat{\textbf{b}}+\hat{\textbf{b}}\times\nabla p/B$ and
(5.10),

$\int_{V_{P}}d^{3}x\,\left(\bm{\xi}_{2}\cdot\textbf{J}\times\delta\widetilde{\textbf{B}}_{1}-\bm{\xi}_{1}\cdot\textbf{J}\times\delta\widetilde{\textbf{B}}_{2}\right)=\int_{V_{P}}d^{3}x\,\frac{J_{||}}{B}\nabla\cdot\left(\left(\bm{\xi}_{1}\times\textbf{B}\right)\times\left(\bm{\xi}_{2}\times\textbf{B}\right)\right)\\\
+\int_{V_{P}}d^{3}x\,\frac{1}{B}\left(\left(\bm{\xi}_{2}\cdot\nabla
p\right)\hat{\textbf{b}}\cdot\delta\widetilde{\textbf{B}}_{1}-\left(\bm{\xi}_{1}\cdot\nabla
p\right)\hat{\textbf{b}}\cdot\delta\widetilde{\textbf{B}}_{2}\right).$ (I.4)

The first term in (I.4) can be simplified using $\nabla\cdot\textbf{J}=0$ and
noting that the perturbation can be written as
$\bm{\xi}_{1,2}=\xi^{\psi}_{1,2}\nabla\psi+\xi^{\perp}_{1,2}\hat{\textbf{b}}\times\nabla\psi$.
Applying the identity
$\textbf{B}\cdot\delta\widetilde{\textbf{B}}_{1,2}=-B^{2}\nabla\cdot\bm{\xi}_{1,2}-\bm{\xi}_{1,2}\cdot\nabla
B^{2}-\mu_{0}\bm{\xi}_{1,2}\cdot\nabla p$ to the second term, the following
expression can be obtained,

$\int_{V_{P}}d^{3}x\,\left(\bm{\xi}_{2}\cdot\textbf{J}\times\delta\widetilde{\textbf{B}}_{1}-\bm{\xi}_{1}\cdot\textbf{J}\times\delta\widetilde{\textbf{B}}_{2}\right)=\\\
\int_{V_{P}}d^{3}x\,\left(\left(\nabla\cdot\bm{\xi}_{2}\right)\bm{\xi}_{1}\cdot\nabla
p-\left(\nabla\cdot\bm{\xi}_{1}\right)\bm{\xi}_{2}\cdot\nabla p\right).$ (I.5)

Hence we obtain the following expression for $U_{P_{1}}$,

$U_{P_{1}}=\int_{V_{P}}d^{3}x\,\left(-\bm{\xi}_{2}\cdot\textbf{F}_{1}+\bm{\xi}_{1}\cdot\textbf{F}_{2}\right)\\\
-\int_{V_{P}}d^{3}x\,\left(\delta\chi_{1}^{\prime}(\psi)\bm{\xi}_{2}\cdot\nabla\psi-\delta\chi_{2}^{\prime}(\psi)\bm{\xi}_{1}\cdot\nabla\psi\right)\textbf{J}\cdot\nabla\varphi.$
(I.6)

We now consider $U_{P_{2}}$ defined in (I.2). Applying (5.24) for the change
in toroidal current, integrating by parts in $\psi$, and combining the
expressions for $U_{P_{1}}$ (I.3) and $U_{P_{2}}$ (I.2), we obtain,

$U_{P}=\int_{V_{P}}d^{3}x\,\left(-\bm{\xi}_{2}\cdot\textbf{F}_{1}+\bm{\xi}_{1}\cdot\textbf{F}_{2}\right)+2\pi\int_{V_{P}}d\psi\,\left(\delta\chi_{1}(\psi)\delta
I_{T,2}^{\prime}(\psi)-\delta\chi_{2}(\psi)\delta
I_{T,1}^{\prime}(\psi)\right)\\\
-\int_{S_{P}}d^{2}x\,\left(\delta\chi_{1}(\psi)\bm{\xi}_{2}-\delta\chi_{2}(\psi)\bm{\xi}_{1}\right)\cdot\hat{\textbf{n}}\textbf{J}\cdot\nabla\varphi.$
(I.7)

Next we combine $U_{P}$ (I.7) with $U_{B}$ (5.31) and $U_{C}$ (5.32) to obtain
the free-boundary adjoint relation (5.33).

To obtain the fixed-boundary adjoint relation, the integral over the plasma
volume (5.29) can be related to a surface integral by applying the divergence
theorem to arrive at 5.35. Using (5.19) and applying several vector
identities,

$U_{P}=-\frac{1}{\mu_{0}}\int_{S_{P}}d^{2}x\,\hat{\textbf{n}}\cdot\left(\bm{\xi}_{1}\delta\textbf{B}_{2}-\bm{\xi}_{2}\delta\textbf{B}_{1}\right)\cdot\textbf{B}\\\
-\frac{1}{\mu_{0}}\int_{S_{P}}d^{2}x\,\left(\delta\chi_{2}(\psi)\delta\textbf{B}_{1}-\delta\chi_{1}(\psi)\delta\textbf{B}_{2}\right)\cdot\nabla\varphi\times\hat{\textbf{n}}.$
(I.8)

Using (I.7) and expressing the second term in (I.8) as a perturbed current
using (5.24), the fixed boundary adjoint relation (5.36) is obtained.

## Appendix J: Alternate derivation of fixed-boundary adjoint relation

The MHD force operator,

$\displaystyle\textbf{F}[\bm{\xi}_{1,2}]=\textbf{J}\times\left(\nabla\times\left(\bm{\xi}_{1,2}\times\textbf{B}\right)\right)+\frac{\nabla\times\left(\nabla\times\left(\bm{\xi}_{1,2}\times\textbf{B}\right)\right)\times\textbf{B}}{\mu_{0}}+\nabla\left(\bm{\xi}_{1,2}\cdot\nabla
p\right),$ (J.1)

possesses the following self-adjointness property [20, 83],

$\displaystyle\int_{V_{P}}d^{3}x\,\left(\bm{\xi}_{2}\cdot\textbf{F}[\bm{\xi}_{1}]-\bm{\xi}_{1}\cdot\textbf{F}[\bm{\xi}_{2}]\right)=\frac{1}{\mu_{0}}\int_{S_{P}}d^{2}x\,\hat{\textbf{n}}\cdot\left(\bm{\xi}_{1}\textbf{B}\cdot\delta\widetilde{\textbf{B}}_{2}-\bm{\xi}_{2}\textbf{B}\cdot\delta\widetilde{\textbf{B}}_{1}\right),$
(J.2)

where
$\delta\widetilde{\textbf{B}}_{1,2}=\nabla\times\left(\bm{\xi}_{1,2}\times\textbf{B}\right)$
is the perturbed field corresponding to the MHD perturbations. As we consider
linearized equilibrium states that preserve $p(\psi)$, the perturbed pressure
satisfies $\delta p(\psi)=-\bm{\xi}\cdot\nabla p$. The force operator we adopt
(J.1) is the $\gamma\rightarrow 0$ limit of the more general form of the force
operator (5.9), which sometimes includes the term $\nabla\left(\gamma
p\nabla\cdot\bm{\xi}\right)$.

For perturbations described by 5.19, 5.20, 5.23, 5.24, 5.25 and 5.26, the
force operator satisfies,

$\displaystyle\textbf{F}[\bm{\xi}_{1,2}]=\textbf{J}\times\left(\nabla\delta\chi_{1,2}(\psi)\times\nabla\varphi\right)+\frac{\nabla\times\left(\nabla\delta\chi_{1,2}(\psi)\times\nabla\varphi\right)\times\textbf{B}}{\mu_{0}}-\delta\textbf{F}_{1,2}.$
(J.3)

Using J.3 and several vector identities, the left hand side of J.2 can be
written as

$\int_{V_{P}}d^{3}x\,\left(\bm{\xi}_{2}\cdot\textbf{F}[\bm{\xi}_{1}]-\bm{\xi}_{1}\cdot\textbf{F}[\bm{\xi}_{2}]\right)=\int_{V_{P}}d^{3}x\,\left(\delta\chi_{1}^{\prime}(\psi)\bm{\xi}_{2}-\delta\chi_{2}^{\prime}(\psi)\bm{\xi}_{1}\right)\cdot\nabla\psi\textbf{J}\cdot\nabla\varphi\\\
-\frac{1}{\mu_{0}}\int_{V_{P}}d^{3}x\,\nabla\psi\times\nabla\varphi\cdot\left(\delta\chi_{1}^{\prime}(\psi)\delta\widetilde{\textbf{B}}_{2}-\delta\chi_{2}^{\prime}(\psi)\delta\widetilde{\textbf{B}}_{1}\right)\\\
-\frac{1}{\mu_{0}}\int_{S_{P}}d^{2}x\,\left(\bm{\xi}_{2}\delta\chi_{1}^{\prime}(\psi)-\bm{\xi}_{1}\delta\chi_{2}^{\prime}(\psi)\right)\cdot\hat{\textbf{n}}\left(\nabla\psi\times\nabla\varphi\cdot\textbf{B}\right)\\\
-\int_{V_{P}}d^{3}x\,\left(\bm{\xi}_{2}\cdot\delta\textbf{F}_{1}-\bm{\xi}_{1}\cdot\delta\textbf{F}_{2}\right).$
(J.4)

In arriving at (J.4), we use $\textbf{J}\cdot\nabla\psi=0$, which follow from
MHD force balance (5.10). Using 5.24 to re-express the first two terms on the
right-hand side,

$\int_{V_{P}}d^{3}x\,\left(\bm{\xi}_{2}\cdot\textbf{F}[\bm{\xi}_{1}]-\bm{\xi}_{1}\cdot\textbf{F}[\bm{\xi}_{2}]\right)=2\pi\int_{V_{P}}d\psi\,\left(\delta
I_{T,2}(\psi)\delta\chi_{1}^{\prime}(\psi)-\delta
I_{T,1}(\psi)\delta\chi_{2}^{\prime}(\psi)\right)\\\
-\frac{1}{\mu_{0}}\int_{S_{P}}d^{2}x\,\left(\bm{\xi}_{2}\delta\chi_{1}^{\prime}(\psi)-\bm{\xi}_{1}\delta\chi_{2}^{\prime}(\psi)\right)\cdot\hat{\textbf{n}}\left(\nabla\psi\times\nabla\varphi\cdot\textbf{B}\right)\\\
-\int_{V_{P}}d^{3}x\,\left(\bm{\xi}_{2}\cdot\delta\textbf{F}_{1}-\bm{\xi}_{1}\cdot\delta\textbf{F}_{2}\right).$
(J.5)

Using 5.19 and J.2 we obtain 5.36.

## Appendix K: Interpretation of the displacement vector

For MHD perturbations such that
$\delta\textbf{B}=\nabla\times\left(\bm{\xi}\times\textbf{B}\right)$ the
displacement can be interpreted as a vector describing the motion of a field
lines. Thus a normal perturbation to the surface of the plasma as in (5.4) can
be expressed in terms of the displacement vector,

$\displaystyle\delta
f(S_{P};\bm{\xi})=\int_{S_{P}}d^{2}x\,\mathcal{G}\bm{\xi}\cdot\hat{\textbf{n}}.$
(K.1)

For perturbations that allow for changes in the rotational transform it
remains to be shown that a similar relation can be found.

As we require that $\psi$ remain a flux surface label in the perturbed
equilibrium, the Lagrangian perturbation to $\psi$ at fixed position is

$\displaystyle\delta\psi=-\delta\textbf{x}\cdot\nabla\psi.$ (K.2)

The perturbed magnetic field,
$\textbf{B}^{\prime}=\textbf{B}+\delta\textbf{B}$ must remain tangent to
$\psi^{\prime}=\psi+\delta\psi$ surfaces; thus to first order in the
perturbation,

$\displaystyle
0=\textbf{B}^{\prime}\cdot\nabla\psi^{\prime}=\textbf{B}\cdot\nabla\delta\psi+\delta\textbf{B}\cdot\nabla\psi.$
(K.3)

Applying the form for the perturbed field allowing for changes in the
rotational transform,
$\delta\textbf{B}=\nabla\times\left(\bm{\xi}\times\textbf{B}-\delta\chi(\psi)\nabla\varphi\right)$,
and using several vector identities, the following condition is obtained

$\displaystyle\textbf{B}\cdot\nabla\left(\delta\textbf{x}\cdot\nabla\psi\right)=\textbf{B}\cdot\nabla\left(\bm{\xi}\cdot\nabla\psi\right).$
(K.4)

This implies that
$\delta\textbf{x}\cdot\nabla\psi=\bm{\xi}\cdot\nabla\psi+F(\psi)$, where
$F(\psi)$ is some flux function which can be determined by requiring that the
perturbation to the toroidal flux as a function of $\psi$ vanishes,
$\delta\Psi_{T}(\psi)=0$.

The perturbed toroidal flux through a surface labeled by $\psi$ contains two
terms, corresponding to the flux of the unperturbed field through the
perturbed surface and the perturbed field through the unperturbed surface,

$\displaystyle\delta\Psi_{T}(\psi)=\int_{\partial
S_{T}(\psi)}d\vartheta\,\sqrt{g}\delta\textbf{x}\cdot\nabla\psi\textbf{B}\cdot\nabla\varphi+\int_{S_{T}(\psi)}d\psi
d\vartheta\,\sqrt{g}\delta\textbf{B}\cdot\nabla\varphi.$ (K.5)

Using the form for $\delta\textbf{B}$, applying the divergence theorem, and
noting that $\textbf{B}\cdot\nabla\varphi=\sqrt{g}^{-1}$, the following
condition is obtained,

$\displaystyle\delta\Psi_{T}(\psi)=\int_{0}^{2\pi}d\vartheta\,\left(\delta\textbf{x}\cdot\nabla\psi-\bm{\xi}\cdot\nabla\psi\right).$
(K.6)

By requiring that $\delta\Psi_{T}(\psi)=0$, we find that $F(\psi)=0$. Thus we
can express shape gradients in the form of (K.1) even when the rotational
transform is allowed to vary.

## Appendix L: Details of axis ripple calculation

In this Appendix, we compute the shape derivative of the finite-pressure
magnetic well figure of merit from (5.101) and show that if we impose an
adjoint perturbation of the form (5.102), the shape gradient is given by
(5.106).

We use the expression for the perturbation to the field strength (5.62) and
$\delta\psi=-\bm{\xi}_{1}\cdot\nabla\psi$ with (5.101) to obtain,

$\delta
f_{R}(S_{P};\bm{\xi}_{1})=\int_{S_{P}}d^{2}x\,\bm{\xi}_{1}\cdot\hat{\textbf{n}}\widetilde{f_{R}}-\int_{V_{P}}d^{3}x\,\frac{\partial\widetilde{f_{R}}}{\partial\psi}\bm{\xi}_{1}\cdot\nabla\psi\\\
-\int_{V_{P}}d^{3}x\,\frac{\partial\widetilde{f_{R}}}{\partial
B}\frac{1}{B}\left(B^{2}\nabla\cdot\bm{\xi}_{1}+\bm{\xi}_{1}\cdot\nabla\left(B^{2}+\mu_{0}p\right)+\delta\chi_{1}^{\prime}(\psi)\textbf{B}\cdot\left(\nabla\psi\times\nabla\varphi\right)\right).$
(L.1)

The third term can be integrated by parts to obtain,

$\delta
f_{R}(S_{P};\bm{\xi}_{1})=\int_{S_{P}}d^{2}x\,\bm{\xi}_{1}\cdot\hat{\textbf{n}}\left(\widetilde{f_{R}}-\frac{\partial\widetilde{f_{R}}}{\partial
B}B\right)+\int_{V_{P}}d^{3}x\,\left(\frac{\partial^{2}\widetilde{f_{R}}}{\partial
B\partial\psi}B-\frac{\partial\widetilde{f_{R}}}{\partial\psi}\right)\bm{\xi}_{1}\cdot\nabla\psi\\\
+\int_{V_{P}}d^{3}x\,\left(-\frac{\partial\widetilde{f_{R}}}{\partial
B}B\bm{\xi}_{1}\cdot\bm{\kappa}+B\frac{\partial^{2}\widetilde{f_{R}}}{\partial
B^{2}}\bm{\xi}_{1}\cdot\nabla
B+\delta\chi_{1}^{\prime}(\psi)\frac{\partial\widetilde{f_{R}}}{\partial
B}\hat{\textbf{b}}\cdot\left(\nabla\varphi\times\nabla\psi\right)\right),$
(L.2)

where the expression for the curvature in an equilibrium field (5.105) has
been applied.

We compute one term that appears in the fixed-boundary adjoint relation (5.36)
using the prescribed adjoint bulk force perturbation (5.102a),

$\int_{V_{P}}d^{3}x\,\bm{\xi}_{1}\cdot\textbf{F}_{2}=\int_{V_{P}}d^{3}x\,\left(-\frac{\partial^{2}p_{||}}{\partial
B\partial\psi}B+\frac{\partial
p_{||}}{\partial\psi}\right)\bm{\xi}_{1}\cdot\nabla\psi\\\
+\int_{V_{P}}d^{3}x\,\left(\frac{\partial p_{||}}{\partial
B}B\bm{\xi}_{1}\cdot\bm{\kappa}-B\frac{\partial^{2}p_{||}}{\partial
B^{2}}\bm{\xi}_{1}\cdot\nabla B\right),$ (L.3)

where we have applied the parallel force balance condition (5.103). Therefore,
if we impose $p_{||}=\widetilde{f_{R}}$, we obtain the following expression
for the shape derivative of $f_{R}$,

$\delta
f_{R}(S_{P};\bm{\xi}_{1})=\int_{S_{P}}d^{2}x\,\bm{\xi}_{1}\cdot\hat{\textbf{n}}\left(\widetilde{f_{R}}-\frac{\partial\widetilde{f_{R}}}{\partial
B}B\right)-\int_{V_{P}}d^{3}x\,\bm{\xi}_{1}\cdot\textbf{F}_{2}\\\
+\int_{V_{P}}d^{3}x\,\delta\chi_{1}^{\prime}(\psi)\frac{\partial\widetilde{f_{R}}}{\partial
B}\hat{\textbf{b}}\cdot\left(\nabla\varphi\times\nabla\psi\right).$ (L.4)

Upon application of the fixed-boundary adjoint relation we obtain (5.106) with
(5.102).

## Appendix M: Details of effective ripple in the $1/\nu$ regime calculation

Neoclassical transport in the $1/\nu$ collisionality regime is discussed in
many references including [65], [42], and [116]. In this Appendix we sketch
the computation of $\epsilon_{\text{eff}}^{3/2}$ originally introduced in
[168] and compute linear perturbations of $f_{\epsilon}$ (5.112), showing them
to take the form of (5.113).

In the $1/\nu$ regime, the distribution function is ordered in the parameter
$\nu_{*}=\nu/(v_{t}/L)\ll 1$, where $\nu$ is the collision frequency, the
thermal speed is $v_{t}=\sqrt{2T/m}$ for mass $m$ and temperature $T$, and $L$
is a macroscopic scale length,

$\displaystyle f_{1}=f_{1}^{-1}+f_{1}^{0}+\mathcal{O}(\nu_{*}).$ (M.1)

In velocity space we use a pitch angle coordinate
$\lambda=v_{\perp}^{2}/(v^{2}B)$, energy coordinate $\epsilon=v^{2}/2$, and
$\sigma=\text{sign}(v_{||})$, where $v_{\perp}=\sqrt{v^{2}-v_{||}^{2}}$ is the
perpendicular velocity and $v_{||}=\textbf{v}\cdot\hat{\textbf{b}}$ is the
parallel velocity. We use the field line label, $\alpha$, and length along a
field line, $l$, to describe location on a constant $\psi$ surface. In the
$1/\nu$ regime the $\textbf{E}\times\textbf{B}$ precession frequency is
assumed to be small relative to the collision frequency, so the drift kinetic
equation (4.2) becomes,

$\displaystyle v_{||}\frac{\partial f_{1}}{\partial
l}=C(f_{1})-\textbf{v}_{\text{m}}\cdot\nabla\psi\frac{\partial
f_{0}}{\partial\psi},$ (M.2)

where the Maxwellian with density $n$ is,

$\displaystyle f_{0}=n\pi^{-3/2}v_{t}^{-3}e^{-v^{2}/v_{t}^{2}},$ (M.3)

and the radial magnetic drift is,

$\displaystyle\textbf{v}_{\text{m}}\cdot\nabla\psi=(v^{2}+v_{||}^{2})\frac{m}{2qB^{3}}\nabla\psi\times\textbf{B}\cdot\nabla
B,$ (M.4)

for charge $q$. The drift kinetic equation to $\mathcal{O}(\nu_{*}^{-1})$ is,

$\displaystyle v_{||}\frac{\partial f_{1}^{-1}}{\partial l}=0.$ (M.5)

In the trapped portion of phase space, this implies that
$f_{1}^{-1}=f_{1}^{-1}(\psi,\alpha,\epsilon,\lambda)$, and in the passing
portion of phase space, this implies that
$f_{1}^{-1}=f_{1}^{-1}(\psi,\epsilon,\lambda,\sigma)$. The drift kinetic
equation to $\mathcal{O}(\nu_{*}^{0})$ is,

$\displaystyle v_{||}\frac{\partial f_{1}^{0}}{\partial
l}=C(f_{1}^{-1})-\textbf{v}_{\text{m}}\cdot\nabla\psi\frac{\partial
f_{0}}{\partial\psi}.$ (M.6)

In the passing region, this implies that $f_{1}^{-1}$ is a Maxwellian, so it
can be taken to vanish. We employ a pitch-angle scattering operator,

$\displaystyle
C=\frac{2\nu(\epsilon)v_{||}}{B\epsilon}\frac{\partial}{\partial\lambda}\left(\lambda
v_{||}\frac{\partial}{\partial\lambda}\right).$ (M.7)

The parallel streaming term in (M.6) is annihilated by the bounce averaging
operation,

$\displaystyle 0=\langle
C(f_{1}^{-1})\rangle_{b}-\langle\textbf{v}_{\text{m}}\cdot\nabla\psi\rangle_{b}\frac{\partial
f_{0}}{\partial\psi},$ (M.8)

where the bounce average of a quantity $A$ is $\langle
A\rangle_{b}=\tau^{-1}\oint dl\,A/v_{||}$ and the bounce time is $\tau=\oint
dl\,v_{||}^{-1}.$ The bounce-averaged equation (M.8) can be expressed in terms
of the parallel adiabatic invariant $J=\oint dl\,v_{||}$ using the relation,

$\displaystyle\langle\textbf{v}_{\text{m}}\cdot\nabla\psi\rangle_{b}=\frac{m}{q\tau}\frac{\partial
J}{\partial\alpha}.$ (M.9)

Integrating (M.8) with respect to $\lambda$ we obtain,

$\displaystyle\frac{\partial
f_{1}^{-1}}{\partial\lambda}=\frac{m\epsilon}{2q\lambda\nu(\epsilon)}\frac{\partial
f_{0}}{\partial\psi}\left(\oint
dl\,\frac{v_{||}}{B}\right)^{-1}\int_{1/B_{\max}}^{\lambda}d\lambda^{\prime}\,\frac{\partial
J}{\partial\alpha}.$ (M.10)

Here $B_{\max}$ is the maximum value of the field strength on the surface
labeled by $\psi$. We have used the boundary condition $\left(\oint
dl\,v_{||}/B\right)\partial
f_{1}^{-1}/\partial\lambda\rvert_{\lambda=1/B_{\max}}=0$, as there is no flux
in pitch-angle from the passing region. The integration with respect to
$\lambda$ is performed to obtain,

$\displaystyle\frac{\partial
f_{1}^{-1}}{\partial\lambda}=-\frac{m}{6q\lambda\nu(\epsilon)}\frac{\partial
f_{0}}{\partial\psi}\left(\oint
dl\,\frac{v_{||}}{B}\right)^{-1}\frac{\partial}{\partial\alpha}\left(\oint
dl\,\frac{v_{||}^{3}}{B}\right).$ (M.11)

The particle flux from $f_{1}^{-1}$ is obtained by multiplying (M.6) by
$f_{1}^{-1}(\partial f_{0}/\partial\psi)^{-1}$, integrating over velocity
space, and flux surface averaging,

$\displaystyle\left\langle\bm{\Gamma}\cdot\nabla\psi\right\rangle_{\psi}$
$\displaystyle\equiv\left\langle\int
d^{3}v\,f_{1}^{-1}\textbf{v}_{\text{m}}\cdot\nabla\psi\right\rangle_{\psi}=\left\langle\int
d^{3}v\,f_{1}^{-1}C(f_{1}^{-1})\left(\frac{\partial
f_{0}}{\partial\psi}\right)^{-1}\right\rangle_{\psi}.$ (M.12)

The velocity space integration is performed using the velocity-space Jacobian
$d^{3}v=2\pi\sum_{\sigma}B\epsilon/|v_{||}|d\lambda d\epsilon$. Upon
integration by parts in $\lambda$ and applying (M.11), the following
expression is obtained,

$\left\langle\bm{\Gamma}\cdot\nabla\psi\right\rangle_{\psi}=\\\
-\frac{4\sqrt{2}\pi}{V^{\prime}(\psi)}\left(\frac{m}{3q}\right)^{2}\int_{0}^{\infty}d\epsilon\,\left(\frac{\partial
f_{0}}{\partial\psi}\right)\frac{\epsilon^{5/2}}{\nu(\epsilon)}\int_{1/B_{\max}}^{1/B_{\min}}\frac{d\lambda}{\lambda}\,\int_{0}^{2\pi}d\alpha\,\sum_{i}\frac{(\frac{\partial}{\partial\alpha}\hat{K}_{i}(\alpha,\lambda))^{2}}{\hat{I}_{i}(\alpha,\lambda)},$
(M.13)

where the bounce integrals are defined by (5.111). The sum in (M.13) is taken
over trapping regions for particles with pitch angle $\lambda$ on a field line
labeled by $\alpha$ for left bounce points $\varphi_{-,i}\in[0,2\pi)$.

The parameter $\epsilon_{\text{eff}}^{3/2}$ quantifies the geometric
dependence of the $1/\nu$ particle flux. It is defined in terms of the radial
particle flux in the following way [168],

$\displaystyle\langle\bm{\Gamma}\cdot\nabla\psi\rangle_{\psi}=-32\langle|\nabla\psi|\rangle_{\psi}^{2}\left(\frac{m}{3q}\right)^{2}\frac{1}{B_{0}^{2}R^{2}}\epsilon_{\text{eff}}^{3/2}\int_{0}^{\infty}d\epsilon\,\left(\frac{\partial
f_{0}}{\partial\psi}\right)\frac{\epsilon^{5/2}}{\nu(\epsilon)}.$ (M.14)

We take our normalizing length and field values to be such that
$B_{0}R=\epsilon_{\text{ref}}^{-1}\langle|\nabla\psi|\rangle_{\psi}$, where
$\epsilon_{\text{ref}}$ is a reference aspect ratio. Comparing (M.13) with
(M.14) we obtain the expression for $\epsilon_{\text{eff}}^{3/2}$ (5.110). The
corresponding expression (29) in [168] is obtained by noting that
$\hat{H}^{\text{Nemov}}=-(\partial\hat{K}/\partial\alpha)\lambda^{1/2}B_{0}^{3/2}$
and $\hat{I}=2\hat{I}^{\text{Nemov}}$, where $\hat{H}^{\text{Nemov}}$ and
$\hat{I}^{\text{Nemov}}$ are given in (30)-(31) of [168].

The shape derivative of $f_{\epsilon}$ (5.112) is computed to be,

$\displaystyle\delta
f_{\epsilon}(S_{P};\bm{\xi}_{1})=\int_{V_{P}}d\psi\,w(\psi)\delta(V^{\prime}(\psi)\epsilon_{\text{eff}}^{3/2}(\psi)).$
(M.15)

The perturbation to the bounce integrals is computed using the following
identity for the perturbation of a line integral
$Q_{L}=\int_{l_{0}}^{l_{L}}dl\,Q$ due to displacement of the integration curve
by vector field $\delta\textbf{x}$ [9, 138],

$\displaystyle\delta
Q_{L}=\int_{l_{0}}^{l_{L}}dl\,\left(\delta\textbf{x}\cdot\left(-\bm{\kappa}Q+\left(\textbf{I}-\hat{\textbf{t}}\hat{\textbf{t}}\right)\cdot\nabla
Q\right)+\delta Q\right)+Q(l_{L})\delta l_{L}-Q(l_{0})\delta l_{0},$ (M.16)

where $\delta Q$ is the perturbation to the integrand at fixed position,
$\hat{\textbf{t}}=\textbf{x}^{\prime}(l)$ is the unit tangent vector,
$\bm{\kappa}=\textbf{x}^{\prime\prime}(l)$ is the curvature, and $\delta
l_{L}$ and $\delta l_{0}$ are perturbations to the bounds of the integral.

We compute the perturbation to the bounce integrals to be,

$\displaystyle\delta\hat{I}_{i}$ $\displaystyle=\oint
dl\,\left(-\frac{v_{||}}{vB}\bm{\kappa}\cdot\delta\textbf{x}-\left(\frac{\lambda
v}{2Bv_{||}}+\frac{v_{||}}{B^{2}v}\right)\left(\delta\textbf{x}\cdot\nabla
B+\delta B\right)\right)$ (M.17a) $\displaystyle\delta\hat{K}_{i}$
$\displaystyle=\oint
dl\,\left(-\frac{v_{||}^{3}}{v^{3}B}\bm{\kappa}\cdot\delta\textbf{x}-\left(\frac{3\lambda
v_{||}}{2Bv}+\frac{v_{||}^{3}}{B^{2}v^{3}}\right)\left(\delta\textbf{x}\cdot\nabla
B+\delta B\right)\right),$ (M.17b)

where $\delta B$ is the perturbation to the field strength (5.62) and
$\delta\textbf{x}$ is given by (5.22). We note that
$\delta\textbf{x}\cdot\hat{\textbf{b}}=0$ such that the perpendicular
projection, $(\textbf{I}-\hat{\textbf{t}}\hat{\textbf{t}})$, is not needed.
There is no contribution due to the perturbation of the bounce points, as the
integrand vanishes at these points. The expressions (5.113)-(5.115b) can now
be obtained by writing (M.15) in terms of the perturbations of the bounce
integrals, using $\bm{\xi}_{1}\cdot\nabla B+\delta
B=-B\left(\textbf{I}-\hat{\textbf{b}}\hat{\textbf{b}}\right)\mathrel{\mathop{\ordinarycolon}}\nabla\bm{\xi}_{1}-\delta\chi_{1}^{\prime}(\psi)\hat{\textbf{b}}\cdot(\nabla\psi\times\nabla\varphi)$
and
$\bm{\kappa}\cdot\bm{\xi}_{1}=-\hat{\textbf{b}}\hat{\textbf{b}}\mathrel{\mathop{\ordinarycolon}}\nabla\bm{\xi}_{1}$.

## Appendix N: Details of departure from quasi-symmetry calculation

In this Appendix we compute the shape derivative of $f_{QS}$ (5.121) to obtain
(5.126)-(5.127c) by expressing each term in (5.125) in the desired form. The
second term in (5.125) is expressed using
$\delta\psi=-\bm{\xi}_{1}\cdot\nabla\psi$,

$\displaystyle\frac{1}{2}\int_{V_{P}}d^{3}x\,w^{\prime}(\psi)\delta\psi\mathcal{M}^{2}=-\frac{1}{2}\int_{V_{P}}d^{3}x\,\mathcal{M}^{2}\bm{\xi}_{1}\cdot\nabla
w(\psi).$ (N.1)

The third term in (5.125) is computed upon application of (5.20), the
divergence theorem, and noting that
$\mathcal{M}=\textbf{B}\cdot\bm{\mathcal{A}}$,

$\int_{V_{P}}d^{3}x\,w(\psi)\mathcal{M}\delta\textbf{B}\cdot\bm{\mathcal{A}}=-\int_{S_{P}}d^{2}x\,\bm{\xi}_{1}\cdot\textbf{{n}}w(\psi)\mathcal{M}^{2}-\int_{V_{P}}d^{3}x\,w(\psi)\delta\chi^{\prime}_{1}(\psi)\mathcal{M}\nabla\psi\times\nabla\varphi\cdot\bm{\mathcal{A}}\\\
+\int_{V_{P}}d^{3}x\,\bm{\xi}_{1}\cdot\left(w(\psi)\mathcal{M}\left(\textbf{B}\times\left(\nabla\times\bm{\mathcal{A}}\right)\right)-\bm{\mathcal{A}}w(\psi)\textbf{B}\cdot\nabla\mathcal{M}+\mathcal{M}\nabla\left(w(\psi)\mathcal{M}\right)\right).$
(N.2)

The quantity $\bm{\mathcal{A}}$ can be projected into the perpendicular
direction as $\bm{\xi}_{1}\cdot\hat{\textbf{b}}=0$, noting that,

$\displaystyle\hat{\textbf{b}}\times\left(\bm{\mathcal{A}}\times\hat{\textbf{b}}\right)=-(\hat{\textbf{b}}\times\nabla\psi)\nabla_{||}B-F(\psi)\nabla_{\perp}B.$
(N.3)

Similarly, any terms in (N.2) involving $\bm{\xi}_{1}\cdot\nabla$ can be
expressed as $\bm{\xi}_{1}\cdot\nabla_{\perp}$. The corresponding terms in
(5.127a) are obtained using the expression for the curvature in an equilibrium
field. The fourth term in (5.125) is expressed in the following way upon
application of (5.62), the divergence theorem, and noting that
$\bm{\mathcal{S}}\cdot\nabla\psi=\nabla\cdot\bm{\mathcal{S}}=0$,

$\int_{V_{P}}d^{3}x\,w(\psi)\mathcal{M}\bm{\mathcal{S}}\cdot\nabla\delta
B=\int_{S_{P}}d^{2}x\,\bm{\xi}_{1}\cdot\mathbf{n}Bw(\psi)\bm{\mathcal{S}}\cdot\nabla\mathcal{M}-\int_{V_{P}}d^{3}x\,\bm{\xi}_{1}\cdot\left[B\nabla\left(w(\psi)\bm{\mathcal{S}}\cdot\nabla\mathcal{M}\right)\right]\\\
+\int_{V_{P}}d^{3}x\,w(\psi)(\bm{\mathcal{S}}\cdot\nabla\mathcal{M})\left(\delta\chi^{\prime}_{1}(\psi)\hat{\textbf{b}}\cdot(\nabla\psi\times\nabla\varphi)+B\bm{\xi}_{1}\cdot\bm{\kappa}\right).$
(N.4)

We express terms involving $\bm{\xi}_{1}\cdot\nabla$ as
$\bm{\xi}_{1}\cdot\nabla_{\perp}$ to obtain the corresponding terms in
(5.127a). The fifth term in (5.125) is expressed in the following way upon
application of $\delta\psi=-\bm{\xi}_{1}\cdot\nabla\psi$, the divergence
theorem, and several vector identities,

$\int_{V_{P}}d^{3}x\,w(\psi)\mathcal{M}\textbf{B}\times\nabla\delta\psi\cdot\nabla
B=-\int_{S_{P}}d^{2}x\,\bm{\xi}_{1}\cdot\hat{\textbf{n}}w(\psi)\mathcal{M}\nabla
B\times\textbf{B}\cdot\nabla\psi\\\
-\int_{V_{P}}d^{3}x\,\bm{\xi}_{1}\cdot\nabla\psi\nabla
B\cdot\nabla\times\left(w(\psi)\mathcal{M}\textbf{B}\right).$ (N.5)

The sixth term in (5.125) upon application of (5.124) is,

$-\int_{V_{P}}d^{3}x\,\frac{\delta
G(\psi)w(\psi)\mathcal{M}\textbf{B}\cdot\nabla B}{\iota(\psi)-(N/M)}=\\\
\frac{1}{4\pi^{2}}\int_{S_{P}}d^{2}x\,\frac{w(\psi)V^{\prime}(\psi)\langle\mathcal{M}\textbf{B}\cdot\nabla
B\rangle_{\psi}}{(\iota(\psi)-(N/M))}\left(\textbf{B}\cdot\nabla\psi\times\nabla\vartheta\right)\bm{\xi}_{1}\cdot\hat{\textbf{n}}\\\
-\frac{1}{4\pi^{2}}\int_{V_{P}}d^{3}x\,\bm{\xi}_{1}\cdot\nabla\left(\frac{w(\psi)V^{\prime}(\psi)\langle\mathcal{M}\textbf{B}\cdot\nabla
B\rangle_{\psi}}{(\iota(\psi)-(N/M))}\right)\textbf{B}\cdot\nabla\psi\times\nabla\vartheta\\\
+\frac{1}{4\pi^{2}}\int_{V_{P}}d^{3}x\,\frac{w(\psi)V^{\prime}(\psi)\langle\mathcal{M}\textbf{B}\cdot\nabla
B\rangle_{\psi}}{\iota(\psi)-(N/M)}\left(\bm{\xi}_{1}\cdot\left(\nabla\psi\nabla\cdot\left(\textbf{B}\times\nabla\vartheta\right)-\textbf{B}\times\nabla\times\left(\nabla\psi\times\nabla\vartheta\right)\right)\right)\\\
-\frac{1}{4\pi^{2}}\int_{V_{P}}d^{3}x\,\delta\chi^{\prime}_{1}(\psi)\frac{w(\psi)V^{\prime}(\psi)\langle\mathcal{M}\textbf{B}\cdot\nabla
B\rangle_{\psi}}{\sqrt{g}^{2}(\iota(\psi)-(N/M))}\frac{\partial\textbf{x}}{\partial\varphi}\cdot\frac{\partial\textbf{x}}{\partial\vartheta}.$
(N.6)

In obtaining the corresponding terms in (5.127a), terms involving
$\bm{\xi}_{1}\cdot\nabla$ are expressed as $\bm{\xi}_{1}\cdot\nabla_{\perp}$.
The seventh term in (5.125) is expressed using
$\delta\psi=-\bm{\xi}_{1}\cdot\nabla\psi$. Combining all terms, we obtain
(5.126)-(5.127c).

## Appendix O: Details of neoclassical figures of merit calculation

In this Section we compute the shape derivative of $f_{NC}$ (5.130) to obtain
(5.136)-(5.137c) by expressing each term in (5.135) in the desired form.
Throughout Boozer coordinates will be assumed.

The second term in (5.135) is expressed using
$\delta\psi=-\bm{\xi}_{1}\cdot\nabla\psi$. The third term in (5.135) can be
computed using (5.124), noting that
$V^{\prime}(\psi)/(4\pi^{2}\sqrt{g})=B^{2}/\langle B^{2}\rangle_{\psi}$ in
Boozer coordinates and applying the divergence theorem,

$\int_{V_{P}}d^{3}x\,w(\psi)\frac{\partial\mathcal{R}(\psi)}{\partial
G(\psi)}\delta
G(\psi)=-\int_{V_{P}}d^{3}x\,w(\psi)\frac{B^{2}\sqrt{g}}{\langle
B^{2}\rangle_{\psi}}\frac{\partial\mathcal{R}(\psi)}{\partial
G(\psi)}\bm{\xi}_{1}\cdot\nabla\psi(\nabla\times\textbf{B})\cdot\nabla\vartheta\\\
+\int_{V_{P}}d^{3}x\,\left(\bm{\xi}_{1}\cdot\nabla\left(\frac{\partial\mathcal{R}(\psi)}{\partial
G(\psi)}\frac{w(\psi)}{\langle
B^{2}\rangle_{\psi}}\right)B^{2}G(\psi)+\frac{w(\psi)}{\langle
B^{2}\rangle_{\psi}}\frac{\partial\mathcal{R}(\psi)}{\partial
G(\psi)}\bm{\xi}_{1}\cdot\textbf{B}\times\nabla\times\left(\frac{\partial\textbf{x}}{\partial\varphi}B^{2}\right)\right)\\\
+\int_{V_{P}}d^{3}x\,\frac{w(\psi)\delta\chi_{1}^{\prime}(\psi)B^{2}}{\sqrt{g}\langle
B^{2}\rangle_{\psi}}\frac{\partial\mathcal{R}(\psi)}{\partial
G(\psi)}\frac{\partial\textbf{x}}{\partial\varphi}\cdot\frac{\partial\textbf{x}}{\partial\vartheta}-\int_{S_{P}}d^{2}x\,w(\psi)\frac{B^{2}}{\langle
B^{2}\rangle_{\psi}}\frac{\partial\mathcal{R}(\psi)}{\partial
G(\psi)}G(\psi)\bm{\xi}_{1}\cdot\hat{\textbf{n}}.$ (O.1)

The fifth term in (5.135) can be computed using (5.62), the divergence
theorem, and the expression for the curvature in an equilibrium field (5.105),

$\int_{V_{P}}d^{3}x\,w(\psi)\langle S_{\mathcal{R}}\delta
B\rangle_{\psi}=\int_{V_{P}}d^{3}x\,\left(\bm{\xi}_{1}\cdot\nabla\left(w(\psi)S_{\mathcal{R}}\right)B-BS_{\mathcal{R}}w(\psi)\bm{\xi}_{1}\cdot\bm{\kappa}\right)\\\
-\int_{V_{P}}d^{3}x\,\delta\chi_{1}^{\prime}(\psi)S_{\mathcal{R}}w(\psi)\hat{\textbf{b}}\cdot\nabla\psi\times\nabla\varphi-\int_{S_{P}}d^{2}x\,w(\psi)S_{\mathcal{R}}B\bm{\xi}_{1}\cdot\hat{\textbf{n}}.$
(O.2)

The resulting terms can be combined to write the shape derivative in the form
of (5.136), noting that any terms involving $\bm{\xi}_{1}\cdot\nabla$ can be
expressed as $\bm{\xi}_{1}\cdot\nabla_{\perp}$.

## Appendix P: Linearized equilibrium energy functional and coefficient
matrices

###  P.1 Further simplification of energy functional

We will now further simplify the energy functional (6.11) using a magnetic
coordinate system. Each of the contravariant components of the perturbed
magnetic field are evaluated to be,

$\displaystyle Q^{\psi}$
$\displaystyle\equiv\delta\textbf{B}[\bm{\xi}]\cdot\nabla\psi=\frac{1}{\sqrt{g}}\left(\frac{\partial\xi^{\psi}}{\partial\varphi}+\iota\frac{\partial\xi^{\psi}}{\partial\vartheta}\right)$
(P.1a) $\displaystyle Q^{\vartheta}$
$\displaystyle\equiv\delta\textbf{B}[\bm{\xi}]\cdot\nabla\vartheta=\frac{1}{\sqrt{g}}\left(\frac{\partial\xi^{\alpha}}{\partial\varphi}-\frac{\partial\xi^{\psi}\iota}{\partial\psi}\right)$
(P.1b) $\displaystyle Q^{\varphi}$
$\displaystyle\equiv\delta\textbf{B}[\bm{\xi}]\cdot\nabla\varphi=-\frac{1}{\sqrt{g}}\left(\frac{\partial\xi^{\alpha}}{\partial\vartheta}+\frac{\partial\xi^{\psi}}{\partial\psi}\right).$
(P.1c)

We also express the current density in the contravariant basis as,

$\displaystyle\textbf{J}=J^{\vartheta}\frac{\partial\textbf{x}}{\partial\vartheta}+J^{\varphi}\frac{\partial\textbf{x}}{\partial\varphi}.$
(P.2)

The first term in the energy functional is expressed as,

$\displaystyle W_{1}$
$\displaystyle\equiv-\frac{1}{\mu_{0}}\int_{V_{P}}d^{3}x\,\delta\textbf{B}[\bm{\xi}]\cdot\delta\textbf{B}[\bm{\xi}]$
(P.3)
$\displaystyle=-\frac{1}{\mu_{0}}\int_{V_{P}}d^{3}x\,\Bigg{[}\left(Q^{\psi}\right)^{2}g_{\psi\psi}+\left(Q^{\vartheta}\right)^{2}g_{\vartheta\vartheta}+\left(Q^{\varphi}\right)^{2}g_{\varphi\varphi}+2Q^{\psi}Q^{\vartheta}g_{\psi\vartheta}\Bigg{]},$

where $g_{x_{i}x_{j}}=\partial\textbf{x}/\partial
x_{i}\cdot\partial\textbf{x}/\partial x_{j}$ are the metric coefficients. Here
we have assumed that $\varphi=\phi$, the geometric toroidal angle, such that
$g_{\vartheta\varphi}=g_{\psi\varphi}=0$.

The second term in the energy functional is expressed as,

$\displaystyle W_{2}$
$\displaystyle\equiv\int_{V_{P}}d^{3}x\,\bm{\xi}\cdot\textbf{J}\times\delta\textbf{B}[\bm{\xi}]$
(P.4)
$\displaystyle=\int_{V_{P}}d^{3}x\,\sqrt{g}\left(\xi^{\psi}\left(J^{\vartheta}Q^{\varphi}-J^{\varphi}Q^{\vartheta}\right)+Q^{\psi}\left(\xi^{\vartheta}J^{\varphi}-\xi^{\varphi}J^{\vartheta}\right)\right).$

Here we can note that the radial component of MHD force balance yields
$p^{\prime}(\psi)=J^{\vartheta}-\iota(\psi)J^{\varphi}$ to write,

$\displaystyle W_{2}$
$\displaystyle=\int_{V_{P}}d^{3}x\,\sqrt{g}\left(\xi^{\psi}\left(J^{\vartheta}Q^{\varphi}-J^{\varphi}Q^{\vartheta}\right)+Q^{\psi}\left(\xi^{\alpha}J^{\varphi}-p^{\prime}(\psi)\xi^{\varphi}\right)\right).$
(P.5)

The third term in the energy functional can be expressed as,

$\displaystyle W_{3}$
$\displaystyle\equiv\int_{V_{P}}d^{3}x\,\bm{\xi}\cdot\nabla\left(\bm{\xi}\cdot\nabla
p\right)$ (P.6)
$\displaystyle=\int_{V_{P}}d^{3}x\,\left[\xi^{\psi}\frac{\partial(\xi^{\psi}p^{\prime}(\psi))}{\partial\psi}+p^{\prime}(\psi)\left(\xi^{\alpha}\frac{\partial\xi^{\psi}}{\partial\vartheta}+\sqrt{g}Q^{\psi}\xi^{\varphi}\right)\right].$

Combining $W_{2}$ and $W_{3}$, we see that the energy functional indeed only
depends on $\xi^{\alpha}$ and $\xi^{\psi}$,

$W_{2}+W_{3}=\int_{V_{P}}d^{3}x\,\Bigg{(}\sqrt{g}\xi^{\psi}\left(J^{\vartheta}Q^{\varphi}-J^{\varphi}Q^{\vartheta}\right)+\xi^{\alpha}\textbf{J}\cdot\nabla\xi^{\psi}+\xi^{\psi}\frac{\partial(\xi^{\psi}p^{\prime}(\psi))}{\partial\psi}\Bigg{)}.$
(P.7)

We now can apply the divegernce theorem, noting that
$\nabla\cdot\textbf{J}=\textbf{J}\cdot\nabla\psi=0$, to obtain,

$\displaystyle
W_{2}+W_{3}=\int_{V_{P}}d^{3}x\,\Bigg{(}\xi^{\psi}\left(J^{\varphi}\iota^{\prime}(\psi)\xi^{\psi}-2\textbf{J}\cdot\nabla\xi^{\alpha}+\xi^{\psi}p^{\prime\prime}(\psi)\right)\Bigg{)}.$
(P.8)

We now see that the first three terms of the energy functional only depend on
$\xi^{\alpha}$ through its $\vartheta$ and $\varphi$ derivatives. Furthermore,
given the restriction of $\delta F_{\alpha}$ discussed in Appendix Q, the
$m=0$, $n=0$ mode of $\xi^{\alpha}$ will not enter the variational principle.

###  P.2 Explicit forms of coefficient matrices

We can now express the linear operators that couple the Fourier components of
$\xi^{\alpha}$, $\xi^{\psi}$, and $\partial\xi^{\psi}/\partial\psi$ given the
simplifications of the energy functional in the previous Section:

$\displaystyle\textbf{A}_{\psi^{\prime}\psi^{\prime}}$
$\displaystyle=-\frac{V^{\prime}(\psi)}{\mu_{0}}\left\langle\frac{1}{\left(\sqrt{g}\right)^{2}}\left(g_{\varphi\varphi}+\iota(\psi)^{2}g_{\vartheta\vartheta}\right)\bm{\mathcal{F}}^{\psi}\bm{\mathcal{F}}^{\psi}\right\rangle_{\psi}$
(P.9a) $\displaystyle\textbf{A}_{\psi\psi}$
$\displaystyle=\frac{V^{\prime}(\psi)}{\mu_{0}}\Bigg{\langle}\frac{1}{\left(\sqrt{g}\right)^{2}}\Bigg{[}-g_{\psi\psi}\Bigg{(}\frac{\partial\bm{\mathcal{F}}^{\psi}}{\partial\varphi}\frac{\partial\bm{\mathcal{F}}^{\psi}}{\partial\varphi}+\iota(\psi)^{2}\frac{\partial\bm{\mathcal{F}}^{\psi}}{\partial\vartheta}\frac{\partial\bm{\mathcal{F}}^{\psi}}{\partial\vartheta}\Bigg{)}$
(P.9b) $\displaystyle\hskip 14.22636pt-
g_{\psi\psi}\iota(\psi)\left(\frac{\partial\bm{\mathcal{F}}^{\psi}}{\partial\vartheta}\frac{\partial\bm{\mathcal{F}}^{\psi}}{\partial\varphi}+\frac{\partial\bm{\mathcal{F}}^{\psi}}{\partial\varphi}\frac{\partial\bm{\mathcal{F}}^{\psi}}{\partial\vartheta}\right)$
$\displaystyle\hskip
14.22636pt+\left(\mu_{0}\left(\sqrt{g}\right)^{2}\left(J^{\varphi}\iota^{\prime}(\psi)+p^{\prime\prime}(\psi)\right)-g_{\vartheta\vartheta}\left(\iota^{\prime}(\psi)\right)^{2}\right)\bm{\mathcal{F}}^{\psi}\bm{\mathcal{F}}^{\psi}$
$\displaystyle\hskip
14.22636pt+g_{\psi\vartheta}\iota^{\prime}(\psi)\left(\left(\frac{\partial\bm{\mathcal{F}}^{\psi}}{\partial\varphi}+\iota(\psi)\frac{\partial\bm{\mathcal{F}}^{\psi}}{\partial\vartheta}\right)\bm{\mathcal{F}}^{\psi}+\bm{\mathcal{F}}^{\psi}\left(\frac{\partial\bm{\mathcal{F}}^{\psi}}{\partial\varphi}+\iota(\psi)\frac{\partial\bm{\mathcal{F}}^{\psi}}{\partial\vartheta}\right)\right)\Bigg{]}\Bigg{\rangle}_{\psi}$
$\displaystyle\textbf{A}_{\psi\psi^{\prime}}$
$\displaystyle=\frac{V^{\prime}(\psi)}{\mu_{0}}\left\langle\frac{2\iota(\psi)}{\left(\sqrt{g}\right)^{2}}\Bigg{[}-\bm{\mathcal{F}}^{\psi}g_{\vartheta\vartheta}\iota^{\prime}(\psi)+g_{\psi\vartheta}\left(\frac{\partial\bm{\mathcal{F}}^{\psi}}{\partial\varphi}+\iota(\psi)\frac{\partial\bm{\mathcal{F}}^{\psi}}{\partial\vartheta}\right)\Bigg{]}\bm{\mathcal{F}}^{\psi}\right\rangle_{\psi}$
(P.9c) $\displaystyle\textbf{A}_{\alpha\alpha}$
$\displaystyle=-\frac{V^{\prime}(\psi)}{\mu_{0}}\left\langle\frac{1}{\left(\sqrt{g}\right)^{2}}\Bigg{[}g_{\vartheta\vartheta}\frac{\partial\bm{\mathcal{F}}^{\alpha}}{\partial\varphi}\frac{\partial\bm{\mathcal{F}}^{\alpha}}{\partial\varphi}+g_{\varphi\varphi}\frac{\partial\bm{\mathcal{F}}^{\alpha}}{\partial\vartheta}\frac{\partial\bm{\mathcal{F}}^{\alpha}}{\partial\vartheta}\Bigg{]}\right\rangle_{\psi}$
(P.9d) $\displaystyle\textbf{A}_{\alpha\psi^{\prime}}$
$\displaystyle=\frac{2V^{\prime}(\psi)}{\mu_{0}}\left\langle\,\frac{1}{\left(\sqrt{g}\right)^{2}}\Bigg{[}g_{\vartheta\vartheta}\iota\frac{\partial\bm{\mathcal{F}}^{\alpha}}{\partial\varphi}-g_{\varphi\varphi}\frac{\partial\bm{\mathcal{F}}^{\alpha}}{\partial\vartheta}\Bigg{]}\bm{\mathcal{F}}^{\psi}\right\rangle_{\psi}$
(P.9e) $\displaystyle\textbf{A}_{\alpha\psi}$
$\displaystyle=-\frac{2V^{\prime}(\psi)}{\mu_{0}}\Bigg{\langle}\frac{1}{\left(\sqrt{g}\right)^{2}}\Bigg{[}\left(-g_{\vartheta\vartheta}\iota^{\prime}(\psi)\frac{\partial\bm{\mathcal{F}}^{\alpha}}{\partial\varphi}+\mu_{0}\left(\sqrt{g}\right)^{2}\textbf{J}\cdot\nabla\bm{\mathcal{F}}^{\alpha}\right)\bm{\mathcal{F}}^{\psi}$
(P.9f) $\displaystyle\hskip
14.22636pt+g_{\psi\vartheta}\frac{\partial\bm{\mathcal{F}}^{\alpha}}{\partial\varphi}\left(\frac{\partial\bm{\mathcal{F}}^{\psi}}{\partial\varphi}+\iota(\psi)\frac{\partial\bm{\mathcal{F}}^{\psi}}{\partial\vartheta}\right)\Bigg{]}\Bigg{\rangle}_{\psi}$
$\displaystyle\textbf{I}_{\psi}$
$\displaystyle=2V^{\prime}(\psi)\left\langle\bm{\mathcal{F}}^{\psi}\delta
F_{\psi}\right\rangle_{\psi}$ (P.9g) $\displaystyle\textbf{I}_{\alpha}$
$\displaystyle=2V^{\prime}(\psi)\left\langle\bm{\mathcal{F}}^{\alpha}\delta
F_{\alpha}\right\rangle_{\psi},$ (P.9h)

where $\left\langle...\right\rangle_{\psi}$ is the flux-surface average
(A.10).

###  P.3 Invertibility of $\textbf{A}_{\alpha\alpha}$

Obtaining the Euler-Lagrange solution for $\xi^{\alpha}$ requires inverting
$\textbf{A}_{\alpha\alpha}$. We now show that this matrix is, in fact,
negative definite and thus invertible. For any non-zero vector
$\bm{\Xi}_{\alpha}$, we can write the inner product with
$\textbf{A}_{\alpha\alpha}$ as,

$\displaystyle\bm{\Xi}_{\alpha}\cdot\left(\textbf{A}_{\alpha\alpha}\bm{\Xi}_{\alpha}\right)=-\frac{1}{\mu_{0}}\int_{0}^{2\pi}d\vartheta\int_{0}^{2\pi}d\varphi\,\left[\frac{g_{\vartheta\vartheta}}{\sqrt{g}}\left(\frac{\partial\xi^{\alpha}}{\partial\varphi}\right)^{2}+\frac{g_{\varphi\varphi}}{\sqrt{g}}\left(\frac{\partial\xi^{\alpha}}{\partial\vartheta}\right)^{2}\right].$
(P.10)

We note that for a well-defined coordinate system, $g_{\vartheta\vartheta}>0$,
$g_{\varphi\varphi}>0$, and $\sqrt{g}>0$. While either
$\partial\xi^{\alpha}/\partial\varphi$ or
$\partial\xi^{\alpha}/\partial\vartheta$ may vanish, they will not vanish
simultaneously throughout the integrand as we have excluded the $n=0$, $m=0$
mode. Therefore, the integrand will only vanish at isolated points. Thus the
above integral is negative definite, and $\textbf{A}_{\alpha\alpha}$ is
invertible throughout the volume.

## Appendix Q: Constraint on bulk force perturbation

As shown in Appendix P, the first three terms in the energy functional (6.11)
only depend on $\xi^{\alpha}$ through its derivatives with respect to
$\vartheta$ and $\varphi$. In this Appendix, we show that it is always
possible to choose the in-surface component of the bulk force perturbation,
$\delta F_{\alpha}$, such that the final term in the energy functional,

$\displaystyle W_{4}\equiv\int_{V_{P}}d^{3}x\,\xi^{\alpha}\delta F_{\alpha},$
(Q.1)

does not depend on $\xi^{\alpha
c}_{0,0}=\frac{1}{(2\pi)^{2}}\int_{0}^{2\pi}d\vartheta\int_{0}^{2\pi}d\varphi\,\xi^{\alpha}$.
As $\xi^{\alpha c}_{0,0}$ does not enter our variational principle, we can
take it to vanish. The condition that $\xi^{\alpha c}_{0,0}$ does not enter
$W_{4}$ is equivalent to requiring that,

$\displaystyle\langle\delta F_{\alpha}\rangle_{\psi}=0,$ (Q.2)

on every surface, where $\langle\dots\rangle_{\psi}$ is the flux-surface
average (A.10). This follows from the surface-averaged in-surface component of
the linearized force-balance equation (6.2),

$\displaystyle\left\langle\frac{\partial\textbf{x}}{\partial\vartheta}\cdot\textbf{F}[\bm{\xi}]\right\rangle_{\psi}=0.$
(Q.3)

This property of the MHD force operator holds for any equilibrium field that
satisfies MHD force balance (6.1). To see this we note that the flux-surface
average can be defined in terms of an average over the infinitesimal volume
between flux surfaces $\Delta V$ (A.12). We can now apply the self-adjointness
relation (6.9) to simplify (Q.3),

$\left\langle\frac{\partial\textbf{x}}{\partial\vartheta}\cdot\textbf{F}[\bm{\xi}]\right\rangle_{\psi}=\left\langle\bm{\xi}\cdot\textbf{F}\left[\frac{\partial\textbf{x}}{\partial\vartheta}\right]\right\rangle_{\psi}\\\
+\lim_{\Delta V\rightarrow 0}\frac{1}{\mu_{0}\Delta
V}\left(\int_{\partial(V_{P}+\Delta
V)}d^{2}x\,\hat{\textbf{n}}\cdot\bm{\xi}\textbf{B}\cdot\delta\textbf{B}\left[\frac{\partial\textbf{x}}{\partial\vartheta}\right]-\int_{\partial(V_{P})}d^{2}x\,\hat{\textbf{n}}\cdot\bm{\xi}\textbf{B}\cdot\delta\textbf{B}\left[\frac{\partial\textbf{x}}{\partial\vartheta}\right]\right),$
(Q.4)

where we have noted that
$\hat{\textbf{n}}\cdot\frac{\partial\textbf{x}}{\partial\vartheta}=0$, as
$\hat{\textbf{n}}\propto\nabla\psi$. The quantity
$\delta\textbf{B}\left[\partial\textbf{x}/\partial\vartheta\right]=\nabla\times\left(\partial\textbf{x}/\partial\vartheta\times\textbf{B}\right)$
is shown to vanish by expressing B in contravariant form and using the dual
relations (A.3) between the contravariant and covariant basis vectors. The
remaining flux-surface averaged term can also be shown to vanish,

$\left\langle\bm{\xi}\cdot\textbf{F}\left[\frac{\partial\textbf{x}}{\partial\vartheta}\right]\right\rangle_{\psi}=\\\
\left\langle\bm{\xi}\cdot\left(\textbf{J}\times\delta\textbf{B}\left[\frac{\partial\textbf{x}}{\partial\vartheta}\right]+\frac{\left(\nabla\times\delta\textbf{B}\left[\frac{\partial\textbf{x}}{\partial\vartheta}\right]\right)\times\textbf{B}}{\mu_{0}}+\nabla\left(\frac{\partial\textbf{x}}{\partial\vartheta}\cdot\nabla
p\right)\right)\right\rangle_{\psi},$ (Q.5)

as $\frac{\partial\textbf{x}}{\partial\varphi}\cdot\nabla\psi=0$ and
$\delta\textbf{B}\left[\partial\textbf{x}/\partial\vartheta\right]=0$.

Therefore, we see that in order to satisfy linear force balance, $\delta
F_{\alpha}$ must be chosen to satisfy the condition (Q.2). However, this
property can always be imparted on a bulk force arising from the adjoint
formulation. Consider the fixed-boundary adjoint relation (5.36) without
perturbations to the rotational transform,

$\displaystyle\int_{V_{P}}d^{3}x\,\left(\bm{\xi}_{1}\cdot\textbf{F}_{2}-\bm{\xi}_{2}\cdot\textbf{F}_{1}\right)-\frac{1}{\mu_{0}}\int_{S_{P}}d^{2}x\,\hat{\textbf{n}}\cdot\left(\bm{\xi}_{2}\delta\textbf{B}[\bm{\xi}_{1}]\cdot\textbf{B}-\bm{\xi}_{1}\delta\textbf{B}[\bm{\xi}_{2}]\cdot\textbf{B}\right)=0.$
(Q.6)

As $\delta\textbf{B}[\bm{\xi}]$ does not depend on $\xi^{\alpha c}_{0,0}$, we
can choose to define the displacement vector such that $\xi^{\alpha
c}_{0,0}=0$. This is analogous to our convention that
$\bm{\xi}\cdot\textbf{B}=0$, as $\delta\textbf{B}[\bm{\xi}]$ does not depend
on the parallel component of $\bm{\xi}$. Given this convention for the
displacement vector, we can note that $\langle\delta
F_{\alpha,2}\rangle_{\psi}$ and $\langle\delta F_{\alpha,1}\rangle_{\psi}$ do
not enter the above adjoint relation. Therefore, we are free to choose our
bulk force such that the desired constraint (Q.2) is satisfied.

## Appendix R: Near-axis expansion of screw pinch equilibria

The MHD force-balance equation for a screw pinch is,

$\displaystyle\frac{d}{dr}\left(\mu_{0}p(r)+\frac{1}{2r^{2}}\left(\psi^{\prime}(r)\right)^{2}\right)+\frac{\iota(r)\psi^{\prime}(r)}{R_{0}^{2}r}\frac{d}{dr}\left(r\iota(r)\psi^{\prime}(r)\right)=0.$
(R.1)

We note that (R.1) remains unchanged under the transformation
$r\rightarrow-r$, so $\psi(r)$ must be even in $r$. Thus near the origin we
can express the flux function as,

$\displaystyle\psi(r)=\frac{\psi_{2}}{2}r^{2}+\frac{\psi_{4}}{4!}r^{4}+\frac{\psi_{6}}{6!}r^{6}+\mathcal{O}(r^{8}),$
(R.2)

under the assumption that $\psi(0)=0$. We similarly express the rotational
transform and pressure profiles in a power series near the axis,

$\displaystyle\iota(\psi(r))$
$\displaystyle=\iota_{0}+\iota_{1}\psi(r)+\frac{\iota_{2}}{2}\psi(r)^{2}+\frac{\iota_{3}}{3!}\psi(r)^{3}+\mathcal{O}(\psi(r)^{4})$
(R.3a) $\displaystyle p(\psi(r))$
$\displaystyle=p_{0}+p_{1}\psi(r)+\frac{p_{2}}{2}\psi(r)^{2}+\frac{p_{3}}{3!}\psi(r)^{3}+\mathcal{O}(\psi(r)^{4}).$
(R.3b)

The force-balance equation to $\mathcal{O}(r)$ becomes,

$\displaystyle\mu_{0}p_{1}\psi_{2}+\frac{2\iota_{0}^{2}\psi_{2}^{2}}{R_{0}^{2}}+\frac{\psi_{2}\psi_{4}}{3}=0,$
(R.4)

and to $\mathcal{O}(r^{3})$ it is,

$\displaystyle\frac{\mu_{0}p_{2}\psi_{2}^{2}}{2}+\frac{3\iota_{0}\iota_{1}\psi_{2}^{3}}{R_{0}^{2}}+\frac{\mu_{0}p_{1}\psi_{4}}{6}+\frac{\iota_{0}^{2}\psi_{2}\psi_{4}}{R_{0}^{2}}+\frac{\psi_{4}^{2}}{18}+\frac{\psi_{2}\psi_{6}}{30}=0.$
(R.5)

In order to determine the power series expansion of $\psi$, we match the
solution near the axis with a numerical solution for $\psi(r)$ at some chosen
boundary location near the axis, $r_{b}$. To perform an expansion to
$\mathcal{O}(r^{2})$, $\psi_{2}$ is chosen such that

$\displaystyle\psi_{2}=\frac{2\psi(r_{b})}{r_{b}^{2}}.$ (R.6)

To perform an expansion to $\mathcal{O}(r^{4})$, (R.4) is used to express
$\psi_{4}$ in terms of $\psi_{2}$, and $\psi_{2}$ is chosen such that
$\psi_{2}r_{b}^{2}/2+\psi_{4}r_{b}^{4}/4!=\psi(r_{b})$,

$\displaystyle\psi_{2}$
$\displaystyle=\frac{-\mu_{0}p_{1}r_{b}^{4}-8\psi(r_{b})}{2r_{b}^{2}\left(\frac{r_{b}^{2}\iota_{0}^{2}}{R_{0}^{2}}-2\right)}$
(R.7a) $\displaystyle\psi_{4}$
$\displaystyle=-3\left(\mu_{0}p_{1}+\frac{2\iota_{0}^{2}\psi_{2}}{R_{0}^{2}}\right).$
(R.7b)

To perform an expansion to $\mathcal{O}(r^{6})$, (R.4) and (R.5) are used to
express $\psi_{4}$ and $\psi_{6}$ in terms of $\psi_{2}$, and $\psi_{2}$ is
chosen such that
$\psi_{2}r_{b}^{2}/2+\psi_{4}r_{b}^{4}/4!+\psi_{6}r_{b}^{6}/6!=\psi(r_{b})$.
The resulting equation for $\psi_{2}$ is quadratic, but only one solution is
allowed in practice to ensure that
$(\psi_{6}r_{b}^{6}/6!)/(\psi_{2}r_{b}^{2}/2+\psi_{4}r_{b}^{4}/4!)\sim
r_{b}^{2}$ in the limit that $r_{b}\ll 1$,

$\displaystyle\psi_{2}$
$\displaystyle=-\frac{R_{0}^{2}}{12\iota_{0}\iota_{1}r_{b}^{6}}\Bigg{(}-24r_{b}^{2}+\frac{12r_{b}^{4}\iota_{0}^{2}}{R_{0}^{2}}+r_{b}^{6}\left(\mu_{0}p_{2}-\frac{8\iota_{0}^{4}}{R_{0}^{4}}\right)$
(R.8a) $\displaystyle\hskip
14.22636pt+r_{b}^{2}\Bigg{[}\left(-24+\mu_{0}p_{2}r_{b}^{4}+\frac{12r_{b}^{2}\iota_{0}^{2}}{R_{0}^{2}}-\frac{8r_{b}^{4}\iota_{0}^{4}}{R_{0}^{4}}\right)^{2}$
$\displaystyle\hskip
14.22636pt+\frac{48r_{b}^{2}\iota_{0}\iota_{1}}{R_{0}^{2}}\left(p_{1}r_{b}^{4}\left(-3+\frac{2r_{b}^{2}\iota_{0}^{2}}{R_{0}^{2}}\right)-24\psi(r_{b})\right)\Bigg{]}^{1/2}\Bigg{)}$
$\displaystyle\psi_{4}$
$\displaystyle=-3\left(\mu_{0}p_{1}+2\frac{\iota_{0}^{2}}{R_{0}^{2}}\psi_{2}\right)$
(R.8b) $\displaystyle\psi_{6}$
$\displaystyle=15\left(\frac{4\mu_{0}p_{1}\iota_{0}^{2}}{R_{0}^{2}}-\mu_{0}p_{2}\psi_{2}+\frac{8\iota_{0}^{4}\psi_{2}}{R_{0}^{4}}-\frac{6\iota_{0}\iota_{2}\psi_{2}^{2}}{R_{0}^{2}}\right).$
(R.8c)

We compare the resulting solution for $\psi$ to a numerical solution of (R.1)
using MATLAB’s bvp4c routine. The solution is computed for $r\in[0,1]$ with a
boundary condition of $\psi(0)=0$ and $\psi(1)=\psi_{0}$. The same profiles
are used as described in Section 6.3.1. The axis expansion solution is matched
with the numerical solution at $r_{b}=10^{-2}$. In Figure R.1 we present a
comparison between the numerical solution and axis expansion of $\psi(r)$. As
expected, the error in the axis expansion to $\mathcal{O}(r^{p})$ scales as
$\sim|r-r_{b}|^{p+2}$ as one moves away from $r=r_{b}$.

(a)

(b)

Figure R.1: (a) The axis expansion solutions to $\mathcal{O}(r^{2})$,
$\mathcal{O}(r^{4})$, and $\mathcal{O}(r^{6})$ are compared with the numerical
solution of $\psi(r)$ near the axis. (b) The absolute error in the expansion
is shown, $\rvert\sum_{n}\psi_{n}r^{n}/n!-\psi(r)\rvert$ where $\psi(r)$ is
the numerical solution. As expected, the error in the axis expansion to
$\mathcal{O}(r^{p})$ scales as $|r-r_{b}|^{p+2}$ near $r=r_{b}$.

## References

  * [1] Princeton plasma physics laboratory - timeline.  URL https://www.pppl.gov/about/history/timeline.  date accessed: 01/03/2019.
  * Abel et al. [2013] I. Abel, G. Plunk, E. Wang, M. Barnes, S. Cowley, W. Dorland, and A. Schekochihin.  Multiscale gyrokinetics for rotating tokamak plasmas: fluctuations, transport and energy flows.  _Reports on Progress in Physics_ , 76(11):116201, 2013.
  * Alexanderian et al. [2017] A. Alexanderian, N. Petra, G. Stadler, and O. Ghattas.  Mean-variance risk-averse optimal control of systems governed by PDEs with random parameter fields using quadratic approximations.  _SIAM/ASA Journal on Uncertainty Quantification_ , 5(1):1166–1192, 2017.
  * Allaire [2015] G. Allaire.  A review of adjoint methods for sensitivity analysis, uncertainty quantification and optimization in numerical codes.  _Ingénieurs de l’Automobile_ , 836:33, 2015.
  * Almagri et al. [1998] A. F. Almagri, D. T. Anderson, and S. F. B. Anderson.  Design and construction of HSX: A helically symmetric stellarator.  In _Helical System Research_. 1998.
  * Anderson [2019] D. Anderson.  Personal communication, 9 2019.
  * Anderson et al. [1990] D. V. Anderson, W. Cooper, R. Gruber, S. Merazzi, and U. Schwenn.  Methods for the efficient calculation of the (MHD) magnetohydrodynamic stability properties of magnetically confined fusion plasmas.  _The International Journal of Supercomputing Applications_ , 4(3):34, 1990.
  * Anderson et al. [1995] F. S. B. Anderson, A. F. Almagri, D. T. Anderson, P. G. Matthews, J. N. Talmadge, and J. L. Shohet.  The Helically Symmetric eXperiment,(HSX) goals, design and status.  _Fusion Technology_ , 27(3T):273–277, 1995.
  * Antonsen and Lee [1982] T. Antonsen and Y. Lee.  Electrostatic modification of variational principles for anisotropic plasmas.  _Physics of Fluids_ , 25(1):132, 1982.
  * Antonsen et al. [2019] T. Antonsen, E. J. Paul, and M. Landreman.  Adjoint approach to calculating shape gradients for three-dimensional magnetic confinement equilibria.  _Journal of Plasma Physics_ , 85(2), 2019.
  * Arnoldus [2006] H. F. Arnoldus.  Conservation of charge at an interface.  _Optics Communications_ , 265(1):52–59, 2006\.
  * Bader et al. [2019] A. Bader, M. Drevlak, D. Anderson, B. Faber, C. Hegna, K. Likin, J. Schmitt, and J. Talmadge.  Stellarator equilibria with reactor relevant energetic particle losses.  _Journal of Plasma Physics_ , 85(5), 2019.
  * Barnes et al. [2010] M. Barnes, I. Abel, W. Dorland, T. Görler, G. Hammett, and F. Jenko.  Direct multiscale coupling of a transport code to gyrokinetic turbulence codes.  _Physics of Plasmas_ , 17(5):056109, 2010.
  * Bauer et al. [2012] F. Bauer, O. Betancourt, and P. Garabedian.  _A Computational Method in Plasma Physics_.  Springer Science & Business Media, 2012.
  * Beidler et al. [1990] C. Beidler, G. Grieger, F. Herrnegger, E. Harmeyer, W. Lotz, H. Maassberg, P. Merkel, J. Nührenberg, F. Rau, J. Sapper, F. Sardei, R. Scardovelli, A. Schlüter, and H. Wobig.  Physics and engineering design for Wendelstein VII-X.  _Fusion Technology_ , 17(1):148, 1990.
  * Beidler et al. [2011] C. Beidler, K. Allmaier, M. Y. Isaev, S. Kasilov, W. Kernbichler, G. Leitold, H. Maassberg, D. Mikkelsen, S. Murakami, M. Schmidt, et al.  Benchmarking of the mono-energetic transport coefficients—results from the International Collaboration on Neoclassical Transport in Stellarators (ICNTS).  _Nuclear Fusion_ , 51(7):076001, 2011.
  * Beidler and D’haeseleer [1995] C. D. Beidler and W. D. D’haeseleer.  A general solution of the ripple-averaged kinetic equation (GSRAKE).  _Plasma Physics and Controlled Fusion_ , 37(4):463, 1995.
  * Belli and Candy [2015] E. A. Belli and J. Candy.  Neoclassical transport in toroidal plasmas with nonaxisymmetric flux surfaces.  _Plasma Physics and Controlled Fusion_ , 57(5):054012, 2015.
  * Berkl et al. [1968] E. Berkl et al.  Plasma physics and controlled nuclear fusion research 1968.  In _Proceedings of the 3rd International Conference Novosibirsk_ , volume 1, 1968.
  * Bernstein et al. [1958] I. Bernstein, E. Frieman, M. Kruskal, and R. Kulsrud.  An energy principle for hydromagnetic stability problems.  _Proceedings of the Royal Society A_ , 244(1236):17, 1958.
  * Boozer [1981] A. Boozer.  Plasma equilibrium with rational magnetic surfaces.  _The Physics of Fluids_ , 24(11):1999, 1981.
  * Boozer [1995] A. Boozer.  Quasi-helical symmetry in stellarators.  _Plasma Physics and Controlled Fusion_ , 37(11A):A103, 1995.
  * Boozer [1980] A. H. Boozer.  Guiding center drift equations.  _The Physics of Fluids_ , 23(5):904, 1980.
  * Boozer [1983] A. H. Boozer.  Transport and isomorphic equilibria.  _The Physics of Fluids_ , 26(2):496, 1983.
  * Boozer [2000] A. H. Boozer.  Stellarator coil optimization by targeting the plasma configuration.  _Physics of Plasmas_ , 7(8):3378, 2000.
  * Boozer [2015] A. H. Boozer.  Non-axisymmetric magnetic fields and toroidal plasma confinement.  _Nuclear Fusion_ , 55(2):025001, 2015.
  * Boozer [2019] A. H. Boozer.  Stellarators as a fast path to fusion energy.  _arXiv preprint arXiv:1912.06289_ , 2019.
  * Boozer and Nührenberg [2006] A. H. Boozer and C. Nührenberg.  Perturbed plasma equilibria.  _Physics of Plasmas_ , 13(10):102501, 2006.
  * Boyd and Vandenberghe [2004] S. Boyd and L. Vandenberghe.  _Convex Optimization_.  Cambridge University Press, 2004.
  * Brent [2013] R. P. Brent.  _Algorithms for Minimization Without Derivatives_.  Courier Corporation, 2013.
  * Brooks and Reiersen [2003] A. Brooks and W. Reiersen.  Coil tolerance impact on plasma surface quality for NCSX.  In _20th IEEE/NPSS Symposium on Fusion Engineering, 2003_ , page 553\. IEEE, 2003.
  * Brown et al. [2015] T. Brown, J. Breslau, D. Gates, N. Pomphrey, and A. Zolfaghari.  Engineering optimization of stellarator coils lead to improvements in device maintenance.  In _IEEE 26th Symposium on Fusion Engineering (SOFE)_ , Austin, Texas, 2015.
  * Calvo et al. [2017] I. Calvo, F. I. Parra, J. L. Velasco, and J. A. Alonso.  The effect of tangential drifts on neoclassical transport in stellarators close to omnigeneity.  _Plasma Physics and Controlled Fusion_ , 59(5):055014, 2017.
  * Calvo et al. [2018] I. Calvo, J. L. Velasco, F. I. Parra, J. A. Alonso, and J. M. García-Regaña.  Electrostatic potential variations on stellarator magnetic surfaces in low collisionality regimes.  _Journal of Plasma Physics_ , 84(4), 2018.
  * Canik et al. [2007] J. Canik, D. Anderson, F. Anderson, C. Clark, K. Likin, J. Talmadge, and K. Zhai.  Reduced particle and heat transport with quasisymmetry in the Helically Symmetric Experiment.  _Physics of Plasmas_ , 14(5):056107, 2007.
  * Carlton-Jones et al. [2019] A. Carlton-Jones, E. Paul, and W. Dorland.  Computing the shape gradient of coil complexity with respect to the plasma boundary with an adjoint method.  _Bulletin of the American Physical Society_ , 64, 2019.
  * Carreras et al. [1996] B. Carreras, V. Lynch, and A. Ware.  Configuration studies for a small-aspect-ratio tokamak stellarator hybrid.  Technical report, Oak Ridge National Lab., 1996.
  * Cary and Hanson [1991] J. R. Cary and J. D. Hanson.  Simple method for calculating island widths.  _Physics of Fluids B: Plasma Physics_ , 3(4):1006, 1991.
  * Cary and Shasharina [1997] J. R. Cary and S. G. Shasharina.  Omnigenity and quasihelicity in helical plasma confinement systems.  _Physics of Plasmas_ , 4(9):3323, 1997.
  * Choi and Kim [2006] K. K. Choi and N.-H. Kim.  _Structural Sensitivity Analysis and Optimization 1: Linear Systems_.  Springer Science & Business Media, 2006.
  * Coddington and Levinson [1955] E. A. Coddington and N. Levinson.  _Theory of Ordinary Differential Equations_.  Tata McGraw-Hill Education, 1955.
  * Connor and Hastie [1974] J. Connor and R. Hastie.  Neoclassical diffusion in an $l=3$ stellarator.  _Physics of Fluids_ , 17(114):114, 1974.
  * Cooper et al. [1992] W. Cooper, S. Hirshman, S. Merazzi, and R. Gruber.  3D magnetohydrodynamic equilibria with anisotropic pressure.  _Computer Physics Communications_ , 72(1):1, 1992\.
  * Cooper et al. [2005] W. Cooper, S. Hirshman, T. Yamaguchi, Y. Narushima, S. Okamura, S. Sakakibara, C. Suzuki, K. Watanabe, H. Yamada, and K. Yamazaki.  Three-dimensional anisotropic pressure equilibria that model balanced tangential neutral beam injection effects.  _Plasma Physics and Controlled Fusion_ , 47(3):561, 2005.
  * Cooper et al. [2006] W. Cooper, J. Graves, S. Hirshman, T. Yamaguchi, Y. Narushima, S. Okamura, S. Sakakibara, C. Suzuki, K. Watanabe, H. Yamada, et al.  Anisotropic pressure bi-Maxwellian distribution function model for three-dimensional equilibria.  _Nuclear Fusion_ , 46(7):683, 2006.
  * Coor et al. [1958] T. Coor, S. Cunningham, R. Ellis, M. Heald, and A. Kranz.  Experiments on the ohmic heating and confinement of plasma in a stellarator.  _The Physics of Fluids_ , 1(5):411, 1958.
  * Dekeyser [2014] W. Dekeyser.  _Optimal Plasma Edge Configurations for Next-Step Fusion Reactors_.  PhD thesis, Katholieke Universiteit Leuven, 2014.
  * Dekeyser et al. [2012] W. Dekeyser, D. Reiter, and M. Baelmans.  Divertor design through shape optimization.  _Contributions to Plasma Physics_ , 52(5):544, 2012.
  * Dekeyser et al. [2014a] W. Dekeyser, D. Reiter, and M. Baelmans.  Automated divertor target design by adjoint shape sensitivity analysis and a one-shot method.  _Journal of Computational Physics_ , 278:117, 2014a.
  * Dekeyser et al. [2014b] W. Dekeyser, D. Reiter, and M. Baelmans.  Optimal shape design for divertors.  _International Journal of Computational Science and Engineering 2_ , 9(5-6):397, 2014b.
  * Dekeyser et al. [2014c] W. Dekeyser, D. Reiter, and M. Baelmans.  A one shot method for divertor target shape optimization.  _Proceedings in Applied Mathematics and Mechanics_ , 14(1):1017, 2014c.
  * Delfour and Zolésio [2011] M. C. Delfour and J.-P. Zolésio.  _Shapes and Geometries_.  Society for Industrial and Applied Mathematics, 2011.
  * Dewar and Hudson [1998] R. Dewar and S. Hudson.  Stellarator symmetry.  _Physica D: Nonlinear Phenomena_ , 112(1):275–280, 1998.
  * D’haeseleer et al. [1991] W. D. D’haeseleer, W. N. Hitchon, J. D. Callen, and J. L. Shohet.  _Flux Coordinates and Magnetic Field Structure: A Guide to a Fundamental Tool of Plasma Theory_.  Springer, 1991.
  * Dinklage et al. [2018] A. Dinklage, C. Beidler, P. Helander, G. Fuchert, H. Maaßberg, K. Rahbarnia, T. S. Pedersen, Y. Turkin, R. Wolf, A. Alonso, et al.  Magnetic configuration effects on the Wendelstein 7-X stellarator.  _Nature Physics_ , 14(8):855–860, 2018.