1 Introduction
While the stellarator is a promising magnetic configuration for the realization of steady-state fusion, the geometry of a stellarator must be carefully designed. This is because collisionless charged particle trajectories are not automatically confined in a stellarator, as they are in axisymmetric configurations. Consequently, the quality of confinement depends sensitively on the shape of the confining magnetic field. To optimize a configuration, figures of merit quantifying confinement, along with other physics criteria such as magnetohydrodynamic (MHD) stability, must be considered in numerical optimization of the MHD equilibrium. These figures of merit describing a configuration depend on the shape of the outer plasma boundary or the shape of the electromagnetic coils. It is thus desirable to obtain derivatives with respect to these shapes for optimization of equilibria or identification of sensitivity information. These so-called shape derivatives can be computed by directly perturbing the shape, recomputing the equilibrium and computing the resulting change to a figure of merit that depends on the equilibrium solution. However, this direct finite-difference approach requires recomputing the equilibrium for each possible perturbation of the shape. For stellarators whose geometry is described by a set of 
$N_{\unicode[STIX]{x1D6FA}}\sim 10^{2}$ parameters, this requires 
$N_{\unicode[STIX]{x1D6FA}}$ solutions to the MHD equilibrium equations. Despite this computational complexity, gradient-based optimization of stellarators has proceeded with the direct approach (e.g. Reiman et al. Reference Reiman, Fu, Hirshman, Ku, Monticello, Mynick, Redi, Spong, Zarnstorff and Blackwell1999; Ku et al. Reference Ku, Garabedian, Lyon, Turnbull, Grossman, Mau and Zarnstorff2008; Proll et al. Reference Proll, Mynick, Xanthopoulos, Lazerson and Faber2015).
The shape gradient quantifies the change in a figure of merit associated with any perturbation to a shape. Thus, if the shape gradient can be obtained, the shape derivative with respect to any perturbation is known (more precise definitions of the shape derivative and gradient are given in § 2). In this work, we provide demonstration of an adjoint approach for computing the shape gradient for functions of MHD equilibria that could be considered within a stellarator configuration optimization. The adjoint approach does not require direct perturbation of a shape, but rather only the solution of one additional force balance equation which depends on the figure of merit of interest. Thus, when 
$N_{\unicode[STIX]{x1D6FA}}$ is very large, as is generally the case for stellarator geometry, the adjoint approach is very advantageous.
In an accompanying work (Antonsen, Paul & Landreman Reference Antonsen, Paul and Landreman2019), two adjoint relations are derived: one involving perturbations to the plasma boundary, referred to as the fixed-boundary adjoint relation, and the other involving perturbations to currents in the vacuum region, known as the free-boundary adjoint relation. These can be considered generalizations of the self-adjointness of the force operator that arises in linearized MHD (Bernstein et al. Reference Bernstein, Frieman, Kruskal and Kulsrud1958). A summary of these results is presented in § 3. These adjoint relations are applied to obtain expressions for the shape gradient of several figures of merit in terms of solutions to an adjoint force balance equation.
Historically, stellarator optimization has been conducted in two stages: in the first, the plasma boundary is varied to optimize an MHD equilibrium for desired physical properties (Nührenberg & Zille Reference Nührenberg and Zille1988). As a second step, the coils are then optimized to provide the desired plasma boundary. The fixed-boundary adjoint relation provides a means to obtain the shape gradient with respect to the plasma boundary and can be used for the traditional optimization route. It is also advantageous to consider coupling the coil design with the physics optimization (Strickler et al. Reference Strickler, Hirshman, Spong, Cole, Lyon, Nelson, Williamson and Ware2004; Drevlak et al. Reference Drevlak, Beidler, Geiger, Helander and Turkin2018; Hudson et al. Reference Hudson, Zhu, Pfefferlé and Gunderson2018), with the aim of obtaining configurations which do not require overly complex coils. The free-boundary adjoint relation allows the computation of the shape gradient of equilibrium figures of merit with respect to coil geometry, allowing for direct optimization of coils. This approach can also be used to efficiently compute coil tolerances (Landreman & Paul Reference Landreman and Paul2018).
Although the adjoint relations are based on the equations of linearized MHD, we perform numerical calculations in this work with nonlinear MHD solutions with the addition of a small perturbation. In the accompanying paper, numerical calculations of the shape gradient with the adjoint approach were obtained for simple figures of merit that did not require modification of the Variational Moments Equilibrium Code (VMEC) (Hirshman & Whitson Reference Hirshman and Whitson1983). In this work, we demonstrate that the adjoint approach can be used to compute the shape gradient for other figures of merit that are relevant for the optimization of stellarator equilibria. We obtain expressions for the shape gradients of the vacuum magnetic well (§ 4), magnetic ripple (§ 5), effective ripple in the 
$1/\unicode[STIX]{x1D708}$ neoclassical regime (Nemov et al. Reference Nemov, Kasilov, Kernbichler and Heyn1999) where 
$\unicode[STIX]{x1D708}$ is the collision frequency (§ 6), departure from quasisymmetry (§ 7) and moments of the neoclassical distribution function (§ 8) in terms of the solution to an adjoint force balance equation. We present calculations of the shape gradient with the adjoint approach for the vacuum magnetic well, which does not require modification to VMEC. The calculation for the magnetic ripple is computed with a minor modification of the Anisotropic Neumann Inverse Moments Equilibrium Code (ANIMEC) (Cooper et al. Reference Cooper, Hirshman, Merazzi and Gruber1992). The adjoint force balance equations needed to compute the shape gradient for the other figures of merit require the addition of a bulk force that will necessitate further modification of an equilibrium or linearized MHD code. Numerical calculations for these figures of merit will, therefore, not be presented in this work.
2 Shape calculus fundamentals
We now introduce several definitions and relations from the field of shape calculus which will prove useful for calculations in this work. Consider a function, 
$F(S_{P})$, which depends implicitly on the plasma boundary, 
$S_{P}$, through the solution to the MHD equilibrium equations with boundary condition 
$\boldsymbol{B}\boldsymbol{\cdot }\boldsymbol{n}|_{S_{P}}=0$ where 
$\boldsymbol{n}$ is the outward unit normal on 
$S_{P}$ and 
$\boldsymbol{B}$ is the magnetic field. We define a functional integrated over the plasma volume, 
$V_{P}$,
$$\begin{eqnarray}\displaystyle f(S_{P})=\int _{V_{P}}\text{d}^{3}x\,F(S_{P}), & & \displaystyle\end{eqnarray}$$ where 
$S_{P}$ is the boundary of 
$V_{P}$. Consider a vector field describing displacements of the surface, 
$\unicode[STIX]{x1D6FF}\boldsymbol{r}$, and a displaced surface 
$S_{P,\unicode[STIX]{x1D716}}=\{\boldsymbol{r}_{0}+\unicode[STIX]{x1D716}\unicode[STIX]{x1D6FF}\boldsymbol{r}:\boldsymbol{r}_{0}\in S_{P}\}$. The shape derivative of 
$F$ is defined as
$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6FF}F(S_{P};\unicode[STIX]{x1D6FF}\boldsymbol{r})=\lim _{\unicode[STIX]{x1D716}\rightarrow 0}\frac{F(S_{P,\unicode[STIX]{x1D716}})-F(S_{P})}{\unicode[STIX]{x1D716}}. & & \displaystyle\end{eqnarray}$$ The shape derivative of 
$f$ is defined by the same expression with 
$F\rightarrow f$. Under certain assumptions of smoothness of 
$\unicode[STIX]{x1D6FF}F$ with respect to 
$\unicode[STIX]{x1D6FF}\boldsymbol{r}$, the shape derivative of the volume-integrated quantity, 
$f$, can be written in the following way (Delfour & Zolésio Reference Delfour and Zolésio2011a),
$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6FF}f(S_{P};\unicode[STIX]{x1D6FF}\boldsymbol{r})=\int _{V_{P}}\text{d}^{3}x\,\unicode[STIX]{x1D6FF}F(S_{P};\unicode[STIX]{x1D6FF}\boldsymbol{r})+\int _{S_{P}}\text{d}^{2}x\,\unicode[STIX]{x1D6FF}\boldsymbol{r}\boldsymbol{\cdot }\boldsymbol{n}F. & & \displaystyle\end{eqnarray}$$ The first term accounts for the Eulerian perturbation to 
$F$ while the second accounts for the motion of the boundary. This is referred to as the transport theorem for domain functionals and will be used throughout to compute the shape derivatives of figures of merit of interest.
According to the Hadamard–Zolésio structure theorem (Delfour & Zolésio Reference Delfour and Zolésio2011b), the shape derivative of a functional of 
$S_{P}$ (not restricted to the form of (2.1)) can be written in the following form,
$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6FF}f(S_{P};\unicode[STIX]{x1D6FF}\boldsymbol{r})=\int _{S_{P}}\text{d}^{2}x\,\unicode[STIX]{x1D6FF}\boldsymbol{r}\boldsymbol{\cdot }\boldsymbol{n}{\mathcal{G}}, & & \displaystyle\end{eqnarray}$$ assuming 
$\unicode[STIX]{x1D6FF}f$ exists for all 
$\unicode[STIX]{x1D6FF}\boldsymbol{r}$ and is sufficiently smooth. In the above expression, 
${\mathcal{G}}$ is the shape gradient. This is an instance of the Riesz representation theorem, which states that any linear functional can be expressed as an inner product with an element of the appropriate space (Rudin Reference Rudin2006). As the shape derivative of 
$f$ is linear in 
$\unicode[STIX]{x1D6FF}\boldsymbol{r}$, it can be written in the form of (2.4). Intuitively, the shape derivative does not depend on tangential perturbations to the surface. The shape gradient can be computed from derivatives with respect to the set of parameters, 
$\unicode[STIX]{x1D6FA}$, used to discretize 
$S_{P}$,
$$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x2202}f}{\unicode[STIX]{x2202}\unicode[STIX]{x1D6FA}_{i}}=\int _{S_{P}}\text{d}^{2}x\,\frac{\unicode[STIX]{x2202}\boldsymbol{r}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D6FA}_{i}}\boldsymbol{\cdot }\boldsymbol{n}{\mathcal{G}}. & & \displaystyle\end{eqnarray}$$ For example, 
$\unicode[STIX]{x1D6FA}=\{R_{mn}^{c},Z_{mn}^{s}\}$ could be assumed, where these are the Fourier coefficients in a cosine and sine representation of the cylindrical coordinates 
$(R,Z)$ of 
$S_{P}$. Upon discretization of the right-hand side on a surface, the above takes the form of a linear system that can be solved for 
${\mathcal{G}}$ (Landreman & Paul Reference Landreman and Paul2018). However, this approach requires performing at least one additional equilibrium calculation for each parameter with a finite-difference approach.
The shape gradient can also be computed with respect to perturbations of currents in the vacuum region. We now consider 
$f$ to depend on the shape of a set of filamentary coils, 
$C=\{C_{k}\}$, through a free-boundary solution to the MHD equilibrium equations. We consider a vector field of displacements to the coils, 
$\unicode[STIX]{x1D6FF}\boldsymbol{r}_{C}$. The shape derivative of 
$f$ can also be written in shape gradient form,
$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6FF}f(C;\unicode[STIX]{x1D6FF}\boldsymbol{r}_{C})=\mathop{\sum }_{k}\int _{C_{k}}\text{d}l\,\boldsymbol{S}_{k}\boldsymbol{\cdot }\unicode[STIX]{x1D6FF}\boldsymbol{r}_{C_{k}}, & & \displaystyle\end{eqnarray}$$ where 
$\boldsymbol{S}_{k}$ is the shape gradient for coil 
$k$, 
$C_{k}$ is the line integral along coil 
$k$ and the sum is taken over coils. Again, 
$\boldsymbol{S}_{k}$ can be computed from derivatives with respect to a set of a parameters describing coil shapes, analogous to (2.5).
To avoid the cost of direct computation of the shape gradient, we apply an adjoint approach. The shape gradient is thus obtained without perturbing the plasma surface or coil shapes directly, but instead by solving an additional adjoint equation that depends on the figure of merit of interest. We perform the calculation with the direct approach to demonstrate that the same derivative information is computed with either approach.
3 Adjoint relations for MHD equilibria
Here we summarize the model for perturbed MHD equilibria and the adjoint relations from the accompanying paper. Throughout we assume the existence of magnetic coordinates such that the magnetic field can be written in the contravariant form as,
$$\begin{eqnarray}\displaystyle \boldsymbol{B}=\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}\times \unicode[STIX]{x1D735}\unicode[STIX]{x1D703}-\unicode[STIX]{x1D704}(\unicode[STIX]{x1D713})\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}\times \unicode[STIX]{x1D735}\unicode[STIX]{x1D701}, & & \displaystyle\end{eqnarray}$$ where 
$2\unicode[STIX]{x03C0}\unicode[STIX]{x1D713}$ is the toroidal flux, 
$\unicode[STIX]{x1D703}$ is a poloidal angle, 
$\unicode[STIX]{x1D701}$ is a toroidal angle and 
$\unicode[STIX]{x1D704}(\unicode[STIX]{x1D713})$ is the rotational transform. The equilibrium magnetic field, 
$\boldsymbol{B}$, is assumed to be in force balance,
$$\begin{eqnarray}\displaystyle \frac{\boldsymbol{J}\times \boldsymbol{B}}{c}=\unicode[STIX]{x1D735}p, & & \displaystyle\end{eqnarray}$$ where 
$\boldsymbol{J}$ is the current density, 
$p(\unicode[STIX]{x1D713})$ is the plasma pressure and 
$c$ is the speed of light. The current density satisfies Ampere’s law,
$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D735}\times \boldsymbol{B}=\frac{4\unicode[STIX]{x03C0}}{c}\boldsymbol{J}. & & \displaystyle\end{eqnarray}$$ We will consider a fixed-boundary calculation such that the equilibrium equations (3.2)–(3.3) are solved with a specified value of toroidal flux 
$2\unicode[STIX]{x03C0}\unicode[STIX]{x1D713}_{0}$ on a given surface 
$S_{P}$, and free-boundary calculations such that they are solved with specified currents in the vacuum region. Two free functions of flux must also be specified, which we take to be 
$p(\unicode[STIX]{x1D713})$ and the rotational transform 
$\unicode[STIX]{x1D704}(\unicode[STIX]{x1D713})$ or the toroidal current contained within a flux surface,
$$\begin{eqnarray}\displaystyle I_{T}(\unicode[STIX]{x1D713})=\frac{c}{8\unicode[STIX]{x03C0}^{2}}\int _{0}^{2\unicode[STIX]{x03C0}}\text{d}\unicode[STIX]{x1D703}\int _{0}^{2\unicode[STIX]{x03C0}}\text{d}\unicode[STIX]{x1D701}\,\sqrt{g}\boldsymbol{B}\boldsymbol{\cdot }(\unicode[STIX]{x1D735}\unicode[STIX]{x1D701}\times \unicode[STIX]{x1D735}\unicode[STIX]{x1D713}), & & \displaystyle\end{eqnarray}$$ where the Jacobian is 
$\sqrt{g}=(\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}\times \unicode[STIX]{x1D735}\unicode[STIX]{x1D703}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D701})^{-1}$. Fixing the toroidal current is more common in the context of stellarator optimization: for 
$\unicode[STIX]{x1D6FD}=0$, 
$I_{T}$ can be taken to vanish while at finite 
$\unicode[STIX]{x1D6FD}$ it can be computed to be self-consistent with the neoclassical bootstrap current (Spong et al. Reference Spong, Hirshman, Berry, Lyon, Fowler, Strickler, Cole, Nelson, Williamson and Ware2001; Shimizu et al. Reference Shimizu, Liu, Isobe, Okamura, Nishimura, Suzuki, Xu, Zhang, Liu and Huang2018). We note that specification of an equilibrium state via (3.1)–(3.4) is not always possible, as magnetic surfaces may not exist and the necessary periodicity constraints on rational surfaces may not be satisfied in a general three-dimensional system. At this point we neglect these issues and proceed assuming that (3.1)–(3.4) are sufficient to define an equilibrium state.
We now consider a linearization about this equilibrium state resulting from a perturbation to the plasma boundary, 
$S_{P}$, the coil shapes, the scalar profiles (
$p(\unicode[STIX]{x1D713})$, 
$\unicode[STIX]{x1D704}(\unicode[STIX]{x1D713})$ or 
$I_{T}(\unicode[STIX]{x1D713})$), or the addition of a bulk force. From (3.1), the perturbed magnetic field can be expressed in terms of the perturbations to the magnetic coordinates, (
$\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D713}$, 
$\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D703}$, 
$\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D701}$),
$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6FF}\boldsymbol{B}=\unicode[STIX]{x1D735}\times [\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D713}(\unicode[STIX]{x1D735}\unicode[STIX]{x1D703}-\unicode[STIX]{x1D704}(\unicode[STIX]{x1D713})\unicode[STIX]{x1D735}\unicode[STIX]{x1D701})+\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}(\unicode[STIX]{x1D704}(\unicode[STIX]{x1D713})\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D701}-\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D703})-\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D6F7}(\unicode[STIX]{x1D713})\unicode[STIX]{x1D735}\unicode[STIX]{x1D701}], & & \displaystyle\end{eqnarray}$$ where 
$\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D6F7}(\unicode[STIX]{x1D713})$ is the perturbation to the poloidal flux profile such that 
$\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D6F7}^{\prime }(\unicode[STIX]{x1D713})=\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D704}(\unicode[STIX]{x1D713})$ is the perturbation to the rotational transform profile. We can express (3.5) in terms of the displacement vector,
$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6FF}\boldsymbol{B}=\unicode[STIX]{x1D735}\times (\unicode[STIX]{x1D743}\times \boldsymbol{B}-\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D6F7}(\unicode[STIX]{x1D713})\unicode[STIX]{x1D735}\unicode[STIX]{x1D701}), & & \displaystyle\end{eqnarray}$$with
$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D743}=\frac{\boldsymbol{B}}{B^{2}}\times [\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D713}(\unicode[STIX]{x1D735}\unicode[STIX]{x1D703}-\unicode[STIX]{x1D704}(\unicode[STIX]{x1D713})\unicode[STIX]{x1D735}\unicode[STIX]{x1D701})+\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}(\unicode[STIX]{x1D704}(\unicode[STIX]{x1D713})\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D701}-\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D703})]. & & \displaystyle\end{eqnarray}$$For perturbations which fix the rotational transform profile, the familiar expression for the perturbed magnetic field from ideal MHD stability theory is recovered.
We define a vector field which defines the displacement of a field line, 
$\unicode[STIX]{x1D6FF}\boldsymbol{r}$, such that the perturbation to the field line label 
$\unicode[STIX]{x1D6FC}=\unicode[STIX]{x1D703}-\unicode[STIX]{x1D704}(\unicode[STIX]{x1D713})\unicode[STIX]{x1D701}$ and toroidal flux satisfy,
$$\begin{eqnarray}\displaystyle \displaystyle \unicode[STIX]{x1D6FF}\unicode[STIX]{x1D713}+\unicode[STIX]{x1D6FF}\boldsymbol{r}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D713} & = & \displaystyle 0\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle \displaystyle \unicode[STIX]{x1D6FF}\unicode[STIX]{x1D6FC}+\unicode[STIX]{x1D6FF}\boldsymbol{r}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D6FC} & = & \displaystyle 0,\end{eqnarray}$$ and 
$\unicode[STIX]{x1D6FF}\boldsymbol{r}\boldsymbol{\cdot }\boldsymbol{B}=0$. Noting that 
$\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D6FC}=\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D703}-\unicode[STIX]{x1D704}(\unicode[STIX]{x1D713})\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D701}-(\unicode[STIX]{x1D704}^{\prime }(\unicode[STIX]{x1D713})\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D713}+\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D6F7}^{\prime }(\unicode[STIX]{x1D713}))\unicode[STIX]{x1D701}$, we find that
$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6FF}\boldsymbol{r}=\unicode[STIX]{x1D743}+\frac{\boldsymbol{b}\times \unicode[STIX]{x1D735}\unicode[STIX]{x1D6F7}(\unicode[STIX]{x1D713})}{B}\unicode[STIX]{x1D701}. & & \displaystyle\end{eqnarray}$$The linearized force balance equation is,
$$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x1D6FF}\boldsymbol{J}_{1,2}\times \boldsymbol{B}+\boldsymbol{J}\times \unicode[STIX]{x1D6FF}\boldsymbol{B}_{1,2}}{c}-\unicode[STIX]{x1D735}\unicode[STIX]{x1D6FF}p_{1,2}+\unicode[STIX]{x1D6FF}\boldsymbol{F}_{1,2}=0, & & \displaystyle\end{eqnarray}$$ where 
$\unicode[STIX]{x1D6FF}\boldsymbol{J}_{1,2}=(c/4\unicode[STIX]{x03C0})\unicode[STIX]{x1D735}\times \unicode[STIX]{x1D6FF}\boldsymbol{B}_{1,2}$ is the perturbed current density, 
$\unicode[STIX]{x1D6FF}p_{1,2}$ is the perturbed pressure and 
$\unicode[STIX]{x1D6FF}\boldsymbol{F}_{1,2}$ is an additional bulk force. For perturbations that fix the pressure profile, 
$p(\unicode[STIX]{x1D713})$, the change to the pressure at fixed position is
$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6FF}p_{1,2}=-\unicode[STIX]{x1D743}_{1,2}\boldsymbol{\cdot }\unicode[STIX]{x1D735}p, & & \displaystyle\end{eqnarray}$$ which follows from (3.7). Quantities with subscript 1 are called the direct perturbation, and those with subscript 2 are called the adjoint perturbation. Direct perturbations correspond to a specified perturbation to the boundary, 
$\unicode[STIX]{x1D6FF}\boldsymbol{r}\boldsymbol{\cdot }\boldsymbol{n}|_{S_{P}}=\unicode[STIX]{x1D743}_{1}\boldsymbol{\cdot }\boldsymbol{n}|_{S_{P}}$, or to the coil shapes, 
$\unicode[STIX]{x1D6FF}\boldsymbol{r}_{C}$, with no additional bulk force or perturbation to the profiles. Adjoint perturbations satisfy a modified force balance equation which depends on the figure of merit of interest. We will discuss several examples of adjoint perturbations in the following sections.
Upon application of the self-adjointness relation for the MHD force operator (Bernstein et al. Reference Bernstein, Frieman, Kruskal and Kulsrud1958) augmented by the introduction of the perturbed poloidal flux, the following fixed-boundary adjoint relation is obtained,
$$\begin{eqnarray}\displaystyle & & \displaystyle \int _{V_{P}}\text{d}^{3}x\,(-\unicode[STIX]{x1D743}_{1}\boldsymbol{\cdot }\unicode[STIX]{x1D6FF}\boldsymbol{F}_{2}+\unicode[STIX]{x1D743}_{2}\boldsymbol{\cdot }\unicode[STIX]{x1D6FF}\boldsymbol{F}_{1})-\frac{2\unicode[STIX]{x03C0}}{c}\int _{V_{P}}\text{d}\unicode[STIX]{x1D713}\,(\unicode[STIX]{x1D6FF}I_{T,2}(\unicode[STIX]{x1D713})\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D704}_{1}(\unicode[STIX]{x1D713})-\unicode[STIX]{x1D6FF}I_{T,1}(\unicode[STIX]{x1D713})\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D704}_{2}(\unicode[STIX]{x1D713}))\nonumber\\ \displaystyle & & \displaystyle \qquad -\,\frac{1}{4\unicode[STIX]{x03C0}}\int _{S_{P}}\text{d}^{2}x\,\boldsymbol{n}\boldsymbol{\cdot }(\unicode[STIX]{x1D743}_{2}\unicode[STIX]{x1D6FF}\boldsymbol{B}_{1}\boldsymbol{\cdot }\boldsymbol{B}-\unicode[STIX]{x1D743}_{1}\unicode[STIX]{x1D6FF}\boldsymbol{B}_{2}\boldsymbol{\cdot }\boldsymbol{B})=0.\end{eqnarray}$$ For perturbed MHD equilibria, the displacement vector 
$\unicode[STIX]{x1D743}$ describes the motion of the boundary 
$(\unicode[STIX]{x1D743}\boldsymbol{\cdot }\boldsymbol{n}|_{S_{P}}=\unicode[STIX]{x1D6FF}\boldsymbol{r}\boldsymbol{\cdot }\boldsymbol{n}|_{S_{P}})$. Thus the shape derivative with respect to the plasma boundary can be expressed with the replacement 
$\unicode[STIX]{x1D6FF}\boldsymbol{r}\rightarrow \unicode[STIX]{x1D743}$ (appendix C of Antonsen et al. (Reference Antonsen, Paul and Landreman2019)). Therefore we see that the boundary term in (3.13) is already in the form of a shape gradient (2.4). The task thus remains to express the shape derivative of a given figure of merit in terms of the first term in (3.13) and convert it to shape gradient form using the fixed-boundary relation.
A similar relation is obtained for perturbation of currents in the vacuum region rather than displacements of the plasma surface,
$$\begin{eqnarray}\displaystyle & & \displaystyle \int _{V_{P}}\text{d}^{3}x\,(-\unicode[STIX]{x1D743}_{1}\boldsymbol{\cdot }\unicode[STIX]{x1D6FF}\boldsymbol{F}_{2}+\unicode[STIX]{x1D743}_{2}\boldsymbol{\cdot }\unicode[STIX]{x1D6FF}\boldsymbol{F}_{1})+\frac{2\unicode[STIX]{x03C0}}{c}\int _{V_{P}}\text{d}\unicode[STIX]{x1D713}\left(\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D6F7}_{1}(\unicode[STIX]{x1D713})\frac{\text{d}\unicode[STIX]{x1D6FF}I_{T,2}(\unicode[STIX]{x1D713})}{\text{d}\unicode[STIX]{x1D713}}-\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D6F7}_{2}(\unicode[STIX]{x1D713})\frac{\text{d}\unicode[STIX]{x1D6FF}I_{T,1}(\unicode[STIX]{x1D713})}{\text{d}\unicode[STIX]{x1D713}}\right)\nonumber\\ \displaystyle & & \displaystyle \qquad +\,\frac{1}{c}\mathop{\sum }_{k}\left(I_{C_{k}}\int _{C_{k}}\text{d}l\,(\unicode[STIX]{x1D6FF}\boldsymbol{r}_{C_{1,k}}(\boldsymbol{x})\boldsymbol{\cdot }\boldsymbol{t}\times \unicode[STIX]{x1D6FF}\boldsymbol{B}_{2}-\unicode[STIX]{x1D6FF}\boldsymbol{r}_{C_{2,k}}(\boldsymbol{x})\boldsymbol{\cdot }\boldsymbol{t}\times \unicode[STIX]{x1D6FF}\boldsymbol{B}_{1})\right)=0,\end{eqnarray}$$ where we have made the assumption that the currents are confined to filamentary coils, and the coil shapes are perturbed without perturbations to their currents, 
$I_{C_{k}}$. Expressions which do not make these assumptions are provided in Antonsen et al. (Reference Antonsen, Paul and Landreman2019). The unit tangent vector along the coil is 
$\boldsymbol{t}$. We can note that the third term in the above expression is in the form of a coil shape gradient (2.6). Thus, this adjoint relation can be applied by expressing the shape derivative of a figure of merit in the form of the first two terms, and the shape gradient is computed from the solution to an adjoint equation.
These relations, equations (3.13) and (3.14), will now be applied to compute the shape gradients for several figures of merit with the adjoint approach.
4 Vacuum magnetic well
The averaged radial (i.e. normal to a flux surface) curvature is an important metric for MHD stability (Freidberg Reference Freidberg2014),
$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D705}_{\unicode[STIX]{x1D713}}\equiv \left\langle \unicode[STIX]{x1D73F}\boldsymbol{\cdot }\left(\frac{\unicode[STIX]{x2202}\boldsymbol{r}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}}\right)_{\unicode[STIX]{x1D6FC},l}\right\rangle _{\unicode[STIX]{x1D713}}=\left\langle \frac{1}{2B^{2}}\left(\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}}(8\unicode[STIX]{x03C0}p+B^{2})\right)_{\unicode[STIX]{x1D6FC},l}\right\rangle _{\unicode[STIX]{x1D713}}, & & \displaystyle\end{eqnarray}$$ where the curvature is 
$\unicode[STIX]{x1D73F}=\boldsymbol{b}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\boldsymbol{b}$, 
$\boldsymbol{b}=\boldsymbol{B}/B$ is a unit vector in the direction of the magnetic field, 
$\unicode[STIX]{x1D6FC}=\unicode[STIX]{x1D703}-\unicode[STIX]{x1D704}(\unicode[STIX]{x1D713})\unicode[STIX]{x1D701}$ is a field line label such that 
$\boldsymbol{B}=\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}\times \unicode[STIX]{x1D735}\unicode[STIX]{x1D6FC}$ and 
$l$ measures length along a field line. Subscripts in the above expression indicate quantities held fixed while computing the derivative. The flux surface average of a quantity 
$A$ is
$$\begin{eqnarray}\displaystyle \langle A\rangle _{\unicode[STIX]{x1D713}}=\frac{\displaystyle \int _{-\infty }^{\infty }\frac{\text{d}l}{B}\,A}{\displaystyle \int _{-\infty }^{\infty }\frac{\text{d}l}{B}}=\frac{\displaystyle \int _{0}^{2\unicode[STIX]{x03C0}}\text{d}\unicode[STIX]{x1D703}\int _{0}^{2\unicode[STIX]{x03C0}}\text{d}\unicode[STIX]{x1D701}\,\sqrt{g}A}{V^{\prime }(\unicode[STIX]{x1D713})}. & & \displaystyle\end{eqnarray}$$ Here, 
$V(\unicode[STIX]{x1D713})$ is the volume enclosed by the surface labelled by 
$\unicode[STIX]{x1D713}$. The average radial curvature appears in the ideal MHD potential energy functional for interchange modes, and it provides a stabilizing effect when 
$p^{\prime }(\unicode[STIX]{x1D713})\unicode[STIX]{x1D705}_{\unicode[STIX]{x1D713}}<0$. As typically 
$p^{\prime }(\unicode[STIX]{x1D713})<0$, 
$\unicode[STIX]{x1D705}_{\unicode[STIX]{x1D713}}>0$ is desirable for MHD stability. In a vacuum field, the expression for the averaged radial curvature reduces to
$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D705}_{\unicode[STIX]{x1D713}}=-\frac{V^{\prime \prime }(\unicode[STIX]{x1D713})}{V^{\prime }(\unicode[STIX]{x1D713})}. & & \displaystyle\end{eqnarray}$$ Thus, as volume increases with flux, 
$V^{\prime \prime }(\unicode[STIX]{x1D713})<0$ is advantageous (Helander Reference Helander2014). The quantity 
$p^{\prime }(\unicode[STIX]{x1D713})V^{\prime \prime }(\unicode[STIX]{x1D713})$ also appears in the Mercier criterion for ideal MHD interchange stability (Mercier & Luc Reference Mercier and Luc1974). Known as the vacuum magnetic well, 
$V^{\prime \prime }(\unicode[STIX]{x1D713})$ has been employed in the optimization of several stellarator configurations (e.g. Hirshman et al. Reference Hirshman, Spong, Whitson, Nelson, Batchelor, Lyon, Sanchez, Brooks, Y.-Fu and Goldston1999; Henneberg et al. Reference Henneberg, Drevlak, Nührenberg, Beidler, Turkin, Loizu and Helander2019).
We consider the following figure of merit
$$\begin{eqnarray}\displaystyle f_{W}=\int _{V_{P}}\text{d}\unicode[STIX]{x1D713}\,w(\unicode[STIX]{x1D713})V^{\prime }(\unicode[STIX]{x1D713}), & & \displaystyle\end{eqnarray}$$ where 
$w(\unicode[STIX]{x1D713})$ is a radial weight function which will be chosen so that (4.4) approximates 
$V^{\prime \prime }(\unicode[STIX]{x1D713})$. This can equivalently be written as
$$\begin{eqnarray}\displaystyle f_{W}=\int _{V_{P}}\text{d}^{3}x\,w(\unicode[STIX]{x1D713}). & & \displaystyle\end{eqnarray}$$4.1 Fixed-boundary shape gradient
We consider perturbations about an equilibrium with fixed toroidal current. For the direct perturbation, we have,
$$\begin{eqnarray}\displaystyle \displaystyle \unicode[STIX]{x1D6FF}\boldsymbol{F}_{1} & = & \displaystyle 0\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle \displaystyle \unicode[STIX]{x1D743}_{1}\boldsymbol{\cdot }\boldsymbol{n}|_{S_{P}} & = & \displaystyle \unicode[STIX]{x1D6FF}\boldsymbol{r}\boldsymbol{\cdot }\boldsymbol{n}|_{S_{P}}\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle \displaystyle \unicode[STIX]{x1D6FF}I_{T,1}(\unicode[STIX]{x1D713}) & = & \displaystyle 0,\end{eqnarray}$$ for a specified boundary perturbation 
$\unicode[STIX]{x1D6FF}\boldsymbol{r}\boldsymbol{\cdot }\boldsymbol{n}$. The shape derivative of 
$f_{W}$ is computed upon application of the transport theorem (2.3), noting that 
$\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D713}=-\unicode[STIX]{x1D743}_{1}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}$,
$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6FF}f_{W}(S_{P};\unicode[STIX]{x1D743}_{1})=-\int _{V_{P}}\text{d}^{3}x\,\unicode[STIX]{x1D743}_{1}\boldsymbol{\cdot }\unicode[STIX]{x1D735}w(\unicode[STIX]{x1D713})+\int _{S_{P}}\text{d}^{2}x\,\unicode[STIX]{x1D743}_{1}\boldsymbol{\cdot }\boldsymbol{n}w(\unicode[STIX]{x1D713}), & & \displaystyle\end{eqnarray}$$ where we have assumed 
$w(\unicode[STIX]{x1D713})$ to be differentiable. We recast the first term in (4.9) as a surface integral by applying the fixed-boundary adjoint relation (3.13) and prescribing the adjoint perturbation to satisfy the following,
$$\begin{eqnarray}\displaystyle \displaystyle \unicode[STIX]{x1D6FF}\boldsymbol{F}_{2} & = & \displaystyle -\unicode[STIX]{x1D735}(\unicode[STIX]{x1D6E5}_{P}w(\unicode[STIX]{x1D713}))\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle \displaystyle \unicode[STIX]{x1D743}_{2}\boldsymbol{\cdot }\boldsymbol{n}|_{S_{P}} & = & \displaystyle 0\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle \displaystyle \unicode[STIX]{x1D6FF}I_{T,2}(\unicode[STIX]{x1D713}) & = & \displaystyle 0.\end{eqnarray}$$ Strictly speaking, the adjoint perturbation is the linear response to the bulk force (4.10). Rather than solve a linearized force balance equation, we note that the adjoint bulk force takes the form of the gradient of a scalar. This is implemented by perturbing the pressure profile by 
$\unicode[STIX]{x1D6E5}_{P}w(\unicode[STIX]{x1D713})$, where 
$\unicode[STIX]{x1D6E5}_{P}$ is a constant chosen judiciously. Thus a small perturbation is applied to the pressure profile, the nonlinear equilibrium is computed, and the change in the fields are recorded. Accordingly 
$\unicode[STIX]{x1D6E5}_{P}$ must be small enough that nonlinear effects are not important, but large enough that round-off error does not dominate.
Upon application of (3.13) we obtain the following expression for the shape gradient which depends on the adjoint solution, 
$\unicode[STIX]{x1D6FF}\boldsymbol{B}_{2}$,
$$\begin{eqnarray}\displaystyle {\mathcal{G}}_{W}=\left(w(\unicode[STIX]{x1D713})+\frac{\unicode[STIX]{x1D6FF}\boldsymbol{B}_{2}\boldsymbol{\cdot }\boldsymbol{B}}{4\unicode[STIX]{x03C0}\unicode[STIX]{x1D6E5}_{P}}\right)_{S_{P}}. & & \displaystyle\end{eqnarray}$$ In figure 1 we present the computation of 
${\mathcal{G}}_{W}$ for the NCSX LI383 equilibrium (Zarnstorff et al. Reference Zarnstorff, Berry, Brooks, Fredrickson, Fu, Hirshman, Hudson, Ku, Lazarus and Mikkelsen2001) using the adjoint and direct approaches. We use a weight function
$$\begin{eqnarray}\displaystyle w(\unicode[STIX]{x1D713})=\exp (-(\unicode[STIX]{x1D713}-\unicode[STIX]{x1D713}_{m,1})^{2}/\unicode[STIX]{x1D713}_{w}^{2})-\exp (-(\unicode[STIX]{x1D713}-\unicode[STIX]{x1D713}_{m,2})^{2}/\unicode[STIX]{x1D713}_{w}^{2}) & & \displaystyle\end{eqnarray}$$ (see figure 1c) such that 
$f_{W}$ remains smooth while it approximates 
$V^{\prime }(\unicode[STIX]{x1D713}_{m,1})-V^{\prime }(\unicode[STIX]{x1D713}_{m,2})$ where 
$\unicode[STIX]{x1D713}_{m,1}=0.8\unicode[STIX]{x1D713}_{0}$, 
$\unicode[STIX]{x1D713}_{m,2}=0.1\unicode[STIX]{x1D713}_{0}$ and 
$\unicode[STIX]{x1D713}_{w}=0.05\unicode[STIX]{x1D713}_{0}$. We note that 
$f_{W}$ can be interpreted as measuring the change in volume due to the interchange of two flux tubes centred at 
$\unicode[STIX]{x1D713}_{m,1}$ and 
$\unicode[STIX]{x1D713}_{m,2}$. If 
$f_{W}>0$, this indicates that moving a flux tube radially outward will cause it to expand and lower the potential energy.
All equilibrium calculations are performed with the VMEC code. For the direct approach, derivatives with respect to the Fourier discretization of the boundary (
$R_{mn}^{c}$ and 
$Z_{mn}^{s}$) are computed for 
$m\leqslant 20$ and 
$|n|\leqslant 10$ using an 8-point centred difference stencil with a polynomial-fitting technique. The direct approach requires 6889 calls to VMEC while the adjoint approach requires two calls. It is clear from figure 1 that the adjoint approach yields the same gradient information as the finite-difference approach, at much lower computational cost. The small difference between figures 1(a) and 1(b) can be quantified as follows,
$$\begin{eqnarray}\displaystyle S_{\text{residual}}=\frac{|S_{\text{adjoint}}-S_{\text{direct}}|}{\displaystyle \sqrt{\left.\int _{S_{P}}\text{d}^{2}x\,S_{\text{adjoint}}^{2}\right/\int _{S_{P}}\text{d}^{2}x}}. & & \displaystyle\end{eqnarray}$$ The surface-averaged value of 
$S_{\text{residual}}$ is 
$3.8\times 10^{-2}$. We note that the number of required equilibrium calculations for the direct shape gradient calculation depends on the Fourier resolution and finite-difference stencil chosen. In this work we present the number of function evaluations required in order for the adjoint and direct shape gradient calculations to agree within a few per cent. As the Fourier resolution is increased, the results of the adjoint and direct methods converge to each other.
The residual difference is non-zero due to several sources of error, including discretization error in VMEC. As a result of the assumption of nested magnetic surfaces, MHD force balance (3.2) is not satisfied exactly, but a finite force residual is introduced. Error is also introduced by computing 
$\unicode[STIX]{x1D6FF}\boldsymbol{B}_{2}$ with the addition of a small perturbation to a nonlinear equilibrium calculation rather than from a linearized MHD solution.
4.2 Coil shape gradient
The shape derivative of 
$f_{W}$ can also be computed with respect to a perturbation of the coil shapes. We consider perturbations about an equilibrium with fixed toroidal current,
$$\begin{eqnarray}\displaystyle \displaystyle \unicode[STIX]{x1D6FF}\boldsymbol{F}_{1} & = & \displaystyle 0\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle \displaystyle \unicode[STIX]{x1D6FF}I_{T,1}(\unicode[STIX]{x1D713}) & = & \displaystyle 0,\end{eqnarray}$$ with specified perturbation to the coil shapes, 
$\unicode[STIX]{x1D6FF}\boldsymbol{r}_{C_{1}}\times \boldsymbol{t}$. We prescribe the following adjoint perturbation
$$\begin{eqnarray}\displaystyle \displaystyle \unicode[STIX]{x1D6FF}\boldsymbol{F}_{2} & = & \displaystyle -\unicode[STIX]{x1D735}(\unicode[STIX]{x1D6E5}_{P}w(\unicode[STIX]{x1D713}))\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle \displaystyle \unicode[STIX]{x1D6FF}I_{T,2}(\unicode[STIX]{x1D713}) & = & \displaystyle 0,\end{eqnarray}$$ with 
$\unicode[STIX]{x1D6FF}\boldsymbol{r}_{C_{2}}\times \boldsymbol{t}=0$. The same weight function (4.14) is applied, which decreases sufficiently fast that we can approximate 
$w(\unicode[STIX]{x1D713}_{0})=0$. Upon application of the free-boundary adjoint relation (3.14), we obtain the following coil shape gradient,
$$\begin{eqnarray}\displaystyle \boldsymbol{S}_{k}=\left.\frac{I_{C_{k}}\boldsymbol{t}\times \unicode[STIX]{x1D6FF}\boldsymbol{B}_{2}}{c\unicode[STIX]{x1D6E5}_{P}}\right|_{C_{k}}. & & \displaystyle\end{eqnarray}$$ The calculation of 
$\boldsymbol{S}_{k}$ for each of the 3 unique coil shapes from the NCSX C09R00 coil setFootnote 1 (Williamson et al. Reference Williamson, Brooks, Brown, Chrzanowski, Cole, Fan, Freudenberg, Fogarty, Hargrove and Heitzenroeder2005) is shown in figure 2. The field is computed in the vacuum region for the evaluation of 
$\unicode[STIX]{x1D6FF}\boldsymbol{B}_{2}$ using the DIAGNO code (Gardner Reference Gardner1990; Lazerson Reference Lazerson2012) with a 2-point centred difference stencil. The shape gradient is also computed with a direct approach. The Cartesian components of each coil are Fourier discretized (
$\boldsymbol{X}_{m}^{s}$, 
$\boldsymbol{X}_{m}^{c}$), and derivatives are computed with respect to 
$m\leqslant 40$ with a 4-point centred difference stencil. The fractional difference between the results obtained with the two approaches is
$$\begin{eqnarray}\displaystyle S_{\text{residual},k}^{l}=\frac{|S_{\text{adjoint},k}^{l}-S_{\text{direct},k}^{l}|}{\displaystyle \sqrt{\left.\int _{C_{k}}\text{d}l\,(S_{\text{adjoint},k}^{l})^{2}\right/\int _{C_{k}}\text{d}l}}. & & \displaystyle\end{eqnarray}$$ The line-averaged value of 
$S_{\text{residual},k}^{l}$ is 
$4.1\times 10^{-2}$. The direct approach required 2917 VMEC calls while the adjoint only required three.
Figure 2. The coil shape gradient for 
$f_{W}$ is calculated for each of the 3 unique NCSX coil shapes. The arrows indicate the direction of 
$\boldsymbol{S}_{k}$ (4.20), and their lengths indicate the magnitude scaled according to the legend.
5 Ripple on magnetic axis
We now consider a figure of merit which quantifies the ripple near the magnetic axis (Carreras, Lynch & Ware Reference Carreras, Lynch and Ware1996; Drevlak et al. Reference Drevlak, Geiger, Helander and Turkin2014, Reference Drevlak, Beidler, Geiger, Helander and Turkin2018). As all physical quantities must be independent of the poloidal angle on the magnetic axis, this quantifies the departure from quasi-helical or quasi-axisymmetry near the magnetic axis.
We define the magnetic ripple to be,
$$\begin{eqnarray}\displaystyle f_{R}=\int _{V_{P}}\text{d}^{3}x\,\widetilde{f}_{R} & & \displaystyle\end{eqnarray}$$with
$$\begin{eqnarray}\displaystyle & \displaystyle \widetilde{f_{R}}={\textstyle \frac{1}{2}}w(\unicode[STIX]{x1D713})(B-\overline{B})^{2} & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \overline{B}=\frac{\displaystyle \int _{V_{P}}\text{d}^{3}x\,w(\unicode[STIX]{x1D713})B}{\displaystyle \int _{V_{P}}\text{d}^{3}x\,w(\unicode[STIX]{x1D713})} & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle w(\unicode[STIX]{x1D713})=\exp (-\unicode[STIX]{x1D713}^{2}/\unicode[STIX]{x1D713}_{w}^{2}) & & \displaystyle\end{eqnarray}$$ with 
$\unicode[STIX]{x1D713}_{w}=0.1\unicode[STIX]{x1D713}_{0}$.
5.1 Fixed-boundary shape gradient
We compute perturbations about an equilibrium with fixed rotational transform
$$\begin{eqnarray}\displaystyle \displaystyle \unicode[STIX]{x1D6FF}\boldsymbol{F}_{1} & = & \displaystyle 0\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle \displaystyle \unicode[STIX]{x1D743}_{1}\boldsymbol{\cdot }\boldsymbol{n}|_{S_{P}} & = & \displaystyle \unicode[STIX]{x1D6FF}\boldsymbol{r}\boldsymbol{\cdot }\boldsymbol{n}|_{S_{P}}\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle \displaystyle \unicode[STIX]{x1D6FF}\unicode[STIX]{x1D704}_{1}(\unicode[STIX]{x1D713}) & = & \displaystyle 0.\end{eqnarray}$$Noting that the local perturbation to the field strength is
$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6FF}B=-\frac{1}{B}(B^{2}\unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D743}_{1}+\unicode[STIX]{x1D743}_{1}\boldsymbol{\cdot }\unicode[STIX]{x1D735}(B^{2}+4\unicode[STIX]{x03C0}p)+\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D704}_{1}(\unicode[STIX]{x1D713})\boldsymbol{B}\boldsymbol{\cdot }(\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}\times \unicode[STIX]{x1D735}\unicode[STIX]{x1D701})), & & \displaystyle\end{eqnarray}$$from (3.5), the shape derivative is computed with the transport theorem (2.3),
$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6FF}f_{R}(S_{P};\unicode[STIX]{x1D743}_{1})=\int _{S_{P}}\text{d}^{2}x\,\unicode[STIX]{x1D743}_{1}\boldsymbol{\cdot }\boldsymbol{n}\widetilde{f_{R}}+\int _{V_{P}}\text{d}^{3}x\,\left(\frac{\unicode[STIX]{x2202}\widetilde{f_{R}}}{\unicode[STIX]{x2202}B}\unicode[STIX]{x1D6FF}B+\frac{\unicode[STIX]{x2202}\widetilde{f_{R}}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}}\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D713}\right), & & \displaystyle\end{eqnarray}$$ where the partial derivative with respect to 
$B$ is performed at constant 
$\unicode[STIX]{x1D713}$. We prescribe the following adjoint perturbation,
$$\begin{eqnarray}\displaystyle \displaystyle \unicode[STIX]{x1D6FF}\boldsymbol{F}_{2} & = & \displaystyle -\unicode[STIX]{x1D6E5}_{P}\unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D64B}\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle \displaystyle \unicode[STIX]{x1D743}_{2}\boldsymbol{\cdot }\boldsymbol{n}|_{S_{P}} & = & \displaystyle 0\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle \displaystyle \unicode[STIX]{x1D6FF}\unicode[STIX]{x1D704}_{2}(\unicode[STIX]{x1D713}) & = & \displaystyle 0,\end{eqnarray}$$ where 
$\unicode[STIX]{x1D6E5}_{P}$ is again a constant scale factor. The bulk force perturbation required for the adjoint problem is written as the divergence of an anisotropic pressure tensor, 
$\unicode[STIX]{x1D64B}=p_{\bot }\unicode[STIX]{x1D644}+(p_{\Vert }-p_{\bot })\boldsymbol{b}\boldsymbol{b}$ where 
$\unicode[STIX]{x1D644}$ is the identity tensor. The parallel and perpendicular pressures are related by the parallel force balance condition,
$$\begin{eqnarray}\displaystyle \left.\frac{\unicode[STIX]{x2202}p_{\Vert }}{\unicode[STIX]{x2202}B}\right|_{\unicode[STIX]{x1D713}}=\frac{p_{\Vert }-p_{\bot }}{B}, & & \displaystyle\end{eqnarray}$$ which follows from the requirement that 
$\boldsymbol{b}\boldsymbol{\cdot }\unicode[STIX]{x1D6FF}\boldsymbol{F}_{2}=0$ (3.11). We take the parallel pressure to be
$$\begin{eqnarray}\displaystyle p_{\Vert }=\widetilde{f_{R}}. & & \displaystyle\end{eqnarray}$$Upon application of the fixed-boundary adjoint relation and the expression for the curvature in an equilibrium field, we obtain the following shape gradient,
$$\begin{eqnarray}\displaystyle {\mathcal{G}}_{R}=\left(p_{\bot }+\frac{\unicode[STIX]{x1D6FF}\boldsymbol{B}_{2}\boldsymbol{\cdot }\boldsymbol{B}}{4\unicode[STIX]{x03C0}\unicode[STIX]{x1D6E5}_{P}}\right)_{S_{P}}. & & \displaystyle\end{eqnarray}$$If instead the toroidal current is held fixed in the direct perturbation as in (4.6)–(4.8), then the required adjoint current perturbation is given by
$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6FF}I_{T,2}(\unicode[STIX]{x1D713})=\frac{c\unicode[STIX]{x1D6E5}_{P}}{2\unicode[STIX]{x03C0}}V^{\prime }(\unicode[STIX]{x1D713})\left\langle \frac{\unicode[STIX]{x2202}\widetilde{f}_{R}}{\unicode[STIX]{x2202}B}\boldsymbol{b}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D701}\times \unicode[STIX]{x1D735}\unicode[STIX]{x1D713}\right\rangle _{\unicode[STIX]{x1D713}}, & & \displaystyle\end{eqnarray}$$with the shape gradient unchanged. See appendix A for details of the calculation.
To compute the adjoint perturbation (5.9)–(5.15), we consider the addition of an anisotropic pressure tensor to the nonlinear force balance equation,
$$\begin{eqnarray}\displaystyle \frac{\boldsymbol{J}^{\prime }\times \boldsymbol{B}^{\prime }}{c}=\unicode[STIX]{x1D735}p^{\prime }+\unicode[STIX]{x1D6E5}_{P}\unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D64B}(\unicode[STIX]{x1D713}^{\prime },B^{\prime }), & & \displaystyle\end{eqnarray}$$ where 
$\unicode[STIX]{x1D64B}(\unicode[STIX]{x1D713}^{\prime },B^{\prime })=p_{\bot }(\unicode[STIX]{x1D713}^{\prime },B^{\prime })\unicode[STIX]{x1D644}+(p_{\Vert }(\unicode[STIX]{x1D713}^{\prime },B^{\prime })-p_{\bot }(\unicode[STIX]{x1D713}^{\prime },B^{\prime }))\boldsymbol{b}^{\prime }\boldsymbol{b}^{\prime }$. Here, primes indicate the perturbed quantities (i.e. 
$B^{\prime }=B+\unicode[STIX]{x1D6FF}B$) where unprimed quantities satisfy (3.2). As in § 4, the perturbation has a scale set by 
$\unicode[STIX]{x1D6E5}_{P}$ which is chosen to be small enough that the response is linear. Enforcing parallel force balance from (5.16) results in the following condition,
$$\begin{eqnarray}\displaystyle \left.\frac{\unicode[STIX]{x2202}p_{\Vert }}{\unicode[STIX]{x2202}B^{\prime }}\right|_{\unicode[STIX]{x1D713}^{\prime }}=\frac{p_{\Vert }-p_{\bot }}{B^{\prime }}. & & \displaystyle\end{eqnarray}$$ If we furthermore assume that 
$\unicode[STIX]{x1D6E5}_{P}\unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D64B}$ is small compared with the other terms in (5.16), we can consider it to be a perturbation to the base equilibrium (3.2). In this way, we can apply the perturbed force balance equation (3.11) with 
$\unicode[STIX]{x1D6FF}\boldsymbol{F}_{2}=-\unicode[STIX]{x1D6E5}_{P}\unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D64B}(\boldsymbol{B})$, where 
$\unicode[STIX]{x1D64B}$ is now evaluated with the equilibrium field which satisfies (3.2). Thus the desired pressure tensor (5.13) can be implemented by evaluating 
$p_{\Vert }$ at the perturbed field such that (5.17) is satisfied.
The pressure tensor defined by (5.12)–(5.13) has been implemented in the ANIMEC code (Cooper et al. Reference Cooper, Hirshman, Merazzi and Gruber1992), which modifies the VMEC variational principle to allow three-dimensional equilibrium solutions with anisotropic pressures to be computed. The ANIMEC code has been used to model equilibria with energetic particle species using pressure tensors based on bi-Maxwellian (Cooper et al. Reference Cooper, Graves, Hirshman, Yamaguchi, Narushima, Okamura, Sakakibara, Suzuki, Watanabe and Yamada2006) and slowing-down (Cooper et al. Reference Cooper, Hirshman, Yamaguchi, Narushima, Okamura, Sakakibara, Suzuki, Watanabe, Yamada and Yamazaki2005) distribution functions. The variational principle assumes that 
$p_{\Vert }$ only varies on a surface through 
$B$ and can, therefore, be used to include the required adjoint bulk force.
In figure 3, we present the computation of 
${\mathcal{G}}_{R}$ for the NCSX LI383 equilibrium using the adjoint and direct approaches. For the direct approach, derivatives with respect to the Fourier discretization of the boundary are computed for 
$m\leqslant 11$ and 
$|n|\leqslant 7$ using an 8-point centred difference stencil. The direct approach required 2761 calls to VMEC while the adjoint approach required two calls. The surface-averaged value of 
$S_{\text{residual}}$ (4.15) is 
$3.3\times 10^{-2}$.
6 Effective ripple in the 
$1/\unicode[STIX]{x1D708}$ regime
The effective ripple in the 
$1/\unicode[STIX]{x1D708}$ regime (Nemov et al. Reference Nemov, Kasilov, Kernbichler and Heyn1999) is a figure of merit which has proven valuable for neoclassical optimization (e.g. Zarnstorff et al. Reference Zarnstorff, Berry, Brooks, Fredrickson, Fu, Hirshman, Hudson, Ku, Lazarus and Mikkelsen2001; Ku et al. Reference Ku, Garabedian, Lyon, Turnbull, Grossman, Mau and Zarnstorff2008; Henneberg et al. Reference Henneberg, Drevlak, Nührenberg, Beidler, Turkin, Loizu and Helander2019). This quantity characterizes the geometric dependence of the neoclassical particle flux under the assumption of low collisionality such that 
$\unicode[STIX]{x1D716}_{\text{eff}}$ is analogous to the helical ripple amplitude, 
$\unicode[STIX]{x1D716}_{h}$, that appears in the expression of the 
$1/\unicode[STIX]{x1D708}$ particle flux for a classical stellarator (Galeev & Sagdeev Reference Galeev and Sagdeev1979). The following expression is obtained for the effective ripple,
$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D716}_{\text{eff}}^{3/2}(\unicode[STIX]{x1D713})=\frac{\unicode[STIX]{x03C0}}{4\sqrt{2}V^{\prime }(\unicode[STIX]{x1D713})\unicode[STIX]{x1D716}_{\text{ref}}^{2}}\int _{1/B_{\max }}^{1/B_{\min }}\frac{\text{d}\unicode[STIX]{x1D706}}{\unicode[STIX]{x1D706}}\,\int _{0}^{2\unicode[STIX]{x03C0}}\text{d}\unicode[STIX]{x1D6FC}\,\mathop{\sum }_{i}\frac{\displaystyle \left(\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D6FC}}\hat{K}_{i}(\unicode[STIX]{x1D6FC},\unicode[STIX]{x1D706})\right)^{2}}{\hat{I} _{i}(\unicode[STIX]{x1D6FC},\unicode[STIX]{x1D706})}. & & \displaystyle\end{eqnarray}$$ Here, 
$\unicode[STIX]{x1D706}=v_{\bot }^{2}/(v^{2}B)$ is the pitch angle, 
$\unicode[STIX]{x1D6FC}=\unicode[STIX]{x1D703}-\unicode[STIX]{x1D704}(\unicode[STIX]{x1D713})\unicode[STIX]{x1D701}$ is a field line label, 
$B_{\min }$ and 
$B_{\max }$ are the minimum and maximum values of the field strength on a surface labelled by 
$\unicode[STIX]{x1D713}$ and 
$\unicode[STIX]{x1D716}_{\text{ref}}$ is a reference aspect ratio. We have defined the bounce integrals
$$\begin{eqnarray}\displaystyle \displaystyle \hat{I} _{i}(\unicode[STIX]{x1D6FC},\unicode[STIX]{x1D706}) & = & \displaystyle \oint \text{d}l\,\frac{v_{\Vert }}{Bv}\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle \displaystyle \hat{K}_{i}(\unicode[STIX]{x1D6FC},\unicode[STIX]{x1D706}) & = & \displaystyle \oint \text{d}l\,\frac{\mathop{v}_{\Vert }^{3}}{Bv^{3}},\end{eqnarray}$$ where the notation 
$\oint \text{d}l=\sum _{\unicode[STIX]{x1D70E}}\unicode[STIX]{x1D70E}\int _{\unicode[STIX]{x1D701}_{-}}^{\unicode[STIX]{x1D701}_{+}}\text{d}\unicode[STIX]{x1D701}/\boldsymbol{b}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D701}$ indicates integration at constant 
$\unicode[STIX]{x1D706}$ and 
$\unicode[STIX]{x1D6FC}$ between successive bounce points where 
$v_{\Vert }(\unicode[STIX]{x1D701}_{+})=v_{\Vert }(\unicode[STIX]{x1D701}_{-})=0$ and 
$\unicode[STIX]{x1D70E}=\text{sign}(v_{\Vert })$. The sum in (6.1) is taken over wells at constant 
$\unicode[STIX]{x1D706}$ and 
$\unicode[STIX]{x1D6FC}$ for 
$\unicode[STIX]{x1D701}_{-,i}\in [0,2\unicode[STIX]{x03C0})$.
We consider an integrated figure of merit
$$\begin{eqnarray}\displaystyle f_{\unicode[STIX]{x1D716}}=\int _{V_{P}}\text{d}^{3}x\,w(\unicode[STIX]{x1D713})\unicode[STIX]{x1D716}_{\text{eff}}^{3/2}(\unicode[STIX]{x1D713}), & & \displaystyle\end{eqnarray}$$ where 
$w(\unicode[STIX]{x1D713})$ is a radial weight function. We perturb about an equilibrium with fixed toroidal current (4.6)–(4.8). The shape derivative of 
$f_{\unicode[STIX]{x1D716}}$ is computed to be
$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6FF}f_{\unicode[STIX]{x1D716}}(S_{P};\unicode[STIX]{x1D743}_{1})=\int _{V_{P}}\text{d}^{3}x\,(\unicode[STIX]{x1D64B}_{\unicode[STIX]{x1D716}}:\unicode[STIX]{x1D735}\unicode[STIX]{x1D743}_{1}+\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D704}_{1}(\unicode[STIX]{x1D713}){\mathcal{I}}_{\unicode[STIX]{x1D716}}), & & \displaystyle\end{eqnarray}$$ where the double dot (:) indicates contraction between dyadic tensors 
$\unicode[STIX]{x1D63C}$ and 
$\unicode[STIX]{x1D63D}$ as
$\unicode[STIX]{x1D63C}:\unicode[STIX]{x1D63D}=\sum _{i,j}A_{ij}B_{ji}$, with
$$\begin{eqnarray}\displaystyle {\mathcal{I}}_{\unicode[STIX]{x1D716}} & = & \displaystyle \frac{\unicode[STIX]{x03C0}w(\unicode[STIX]{x1D713})}{2\sqrt{2}\unicode[STIX]{x1D716}_{\text{ref}}^{2}}\int _{1/B_{\max }}^{1/B}\frac{\text{d}\unicode[STIX]{x1D706}}{\unicode[STIX]{x1D706}}\nonumber\\ \displaystyle & & \displaystyle \times \left[\frac{\displaystyle \left(\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D6FC}}\hat{K}(\unicode[STIX]{x1D6FC},\unicode[STIX]{x1D706},\unicode[STIX]{x1D701})\right)^{2}}{\hat{I} ^{2}(\unicode[STIX]{x1D6FC},\unicode[STIX]{x1D706},\unicode[STIX]{x1D701})}\left(-\unicode[STIX]{x1D701}\boldsymbol{B}\times \unicode[STIX]{x1D735}\unicode[STIX]{x1D713}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\left(\frac{|v_{\Vert }|}{vB^{2}}\right)+\boldsymbol{B}\times \unicode[STIX]{x1D735}\unicode[STIX]{x1D713}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D701}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}B}\left(\frac{|v_{\Vert }|}{vB}\right)\right)\right.\nonumber\\ \displaystyle & & \displaystyle +\,2\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D6FC}}\left(\frac{\displaystyle \frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D6FC}}\hat{K}(\unicode[STIX]{x1D6FC},\unicode[STIX]{x1D706},\unicode[STIX]{x1D701})}{\hat{I} (\unicode[STIX]{x1D6FC},\unicode[STIX]{x1D706},\unicode[STIX]{x1D701})}\right)\nonumber\\ \displaystyle & & \displaystyle \times \left.\left(-\unicode[STIX]{x1D701}\boldsymbol{B}\times \unicode[STIX]{x1D735}\unicode[STIX]{x1D713}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\left(\frac{|v_{\Vert }|^{3}}{v^{3}B^{2}}\right)+\boldsymbol{B}\times \unicode[STIX]{x1D735}\unicode[STIX]{x1D713}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D701}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}B}\left(\frac{|v_{\Vert }|^{3}}{v^{3}B}\right)\right)\right]\end{eqnarray}$$ and 
$\unicode[STIX]{x1D64B}_{\unicode[STIX]{x1D716}}=p_{\Vert }\boldsymbol{b}\boldsymbol{b}+p_{\bot }(\unicode[STIX]{x1D644}-\boldsymbol{b}\boldsymbol{b})$ with
$$\begin{eqnarray}\displaystyle p_{\Vert } & = & \displaystyle -\frac{\unicode[STIX]{x03C0}w(\unicode[STIX]{x1D713})}{2\sqrt{2}\unicode[STIX]{x1D716}_{\text{ref}}^{2}}\int _{1/B_{\max }}^{1/B}\frac{\text{d}\unicode[STIX]{x1D706}}{\unicode[STIX]{x1D706}}\,\left(\frac{\displaystyle \left(\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D6FC}}\hat{K}(\unicode[STIX]{x1D6FC},\unicode[STIX]{x1D706},\unicode[STIX]{x1D701})\right)^{2}}{\hat{I} ^{2}(\unicode[STIX]{x1D6FC},\unicode[STIX]{x1D706},\unicode[STIX]{x1D701})}\frac{|v_{\Vert }|}{v}+2\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D6FC}}\left(\frac{\displaystyle \frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D6FC}}\hat{K}(\unicode[STIX]{x1D6FC},\unicode[STIX]{x1D706},\unicode[STIX]{x1D701})}{\hat{I} (\unicode[STIX]{x1D6FC},\unicode[STIX]{x1D706},\unicode[STIX]{x1D701})}\right)\frac{|v_{\Vert }|^{3}}{v^{3}}\right)\nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle p_{\bot } & = & \displaystyle -\frac{\unicode[STIX]{x03C0}w(\unicode[STIX]{x1D713})}{2\sqrt{2}\unicode[STIX]{x1D716}_{\text{ref}}^{2}}\int _{1/B_{\max }}^{1/B}\frac{\text{d}\unicode[STIX]{x1D706}}{\unicode[STIX]{x1D706}}\,\left(\frac{\displaystyle \left(\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D6FC}}\hat{K}(\unicode[STIX]{x1D6FC},\unicode[STIX]{x1D706},\unicode[STIX]{x1D701})\right)^{2}}{\hat{I} ^{2}(\unicode[STIX]{x1D6FC},\unicode[STIX]{x1D706},\unicode[STIX]{x1D701})}\left(\frac{\unicode[STIX]{x1D706}vB}{2|v_{\Vert }|}+\frac{|v_{\Vert }|}{v}\right)\right.\nonumber\\ \displaystyle & & \displaystyle +\left.2\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D6FC}}\left(\frac{\displaystyle \frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D6FC}}\hat{K}(\unicode[STIX]{x1D6FC},\unicode[STIX]{x1D706},\unicode[STIX]{x1D701})}{\hat{I} (\unicode[STIX]{x1D6FC},\unicode[STIX]{x1D706},\unicode[STIX]{x1D701})}\right)\left(\frac{3\unicode[STIX]{x1D706}|v_{\Vert }|B}{2v}+\frac{|v_{\Vert }|^{3}}{v^{3}}\right)\right).\end{eqnarray}$$ Derivatives are computed assuming 
$\unicode[STIX]{x1D716}_{\text{ref}}$ is held constant. The bounce integrals are defined with respect to 
$\unicode[STIX]{x1D701}$ such that 
$\hat{I} (\unicode[STIX]{x1D6FC},\unicode[STIX]{x1D706},\unicode[STIX]{x1D701})=\hat{I} _{i}$ if 
$\unicode[STIX]{x1D701}\in [\unicode[STIX]{x1D701}_{-,i},\unicode[STIX]{x1D701}_{+,i}]$ and 
$\hat{I} (\unicode[STIX]{x1D6FC},\unicode[STIX]{x1D706},\unicode[STIX]{x1D701})=0$ if 
$\unicode[STIX]{x1D706}B(\unicode[STIX]{x1D6FC},\unicode[STIX]{x1D701})>1$. The same convention is used for 
$\hat{K}(\unicode[STIX]{x1D6FC},\unicode[STIX]{x1D706},\unicode[STIX]{x1D701})$. We prescribe the following adjoint perturbation
$$\begin{eqnarray}\displaystyle \displaystyle \unicode[STIX]{x1D6FF}\boldsymbol{F}_{2} & = & \displaystyle -\unicode[STIX]{x1D6E5}_{P}\unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D64B}_{\unicode[STIX]{x1D716}}\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle \displaystyle \unicode[STIX]{x1D743}_{2}\boldsymbol{\cdot }\boldsymbol{n}|_{S_{P}} & = & \displaystyle 0\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle \displaystyle \unicode[STIX]{x1D6FF}I_{T,2}(\unicode[STIX]{x1D713}) & = & \displaystyle \frac{c}{2\unicode[STIX]{x03C0}}V^{\prime }(\unicode[STIX]{x1D713})\unicode[STIX]{x1D6E5}_{P}\langle {\mathcal{I}}_{\unicode[STIX]{x1D716}}\rangle _{\unicode[STIX]{x1D713}}.\end{eqnarray}$$The adjoint bulk force must be consistent with parallel force balance from (3.11), which is equivalent to the condition
$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D735}_{\Vert }p_{\Vert }=\frac{\unicode[STIX]{x1D735}_{\Vert }B}{B}(p_{\Vert }-p_{\bot }). & & \displaystyle\end{eqnarray}$$ This can be shown to be satisfied by (6.7)–(6.8), noting that the 
$\unicode[STIX]{x1D706}$ integrand vanishes at 
$1/B$ such that there is no contribution from the parallel gradient acting on the bounds of the integral. There is also no contribution to the parallel gradient from the bounce integrals, as 
$|v_{\Vert }|$ vanishes at points of non-zero gradient of 
$\hat{I} (\unicode[STIX]{x1D6FC},\unicode[STIX]{x1D706},\unicode[STIX]{x1D701})$ and 
$\hat{K}(\unicode[STIX]{x1D6FC},\unicode[STIX]{x1D706},\unicode[STIX]{x1D701})$.
Upon application of the fixed-boundary adjoint relation (3.13) and integration by parts, we obtain the following expression for the shape gradient
$$\begin{eqnarray}\displaystyle {\mathcal{G}}_{\unicode[STIX]{x1D716}}=\left(p_{\bot }+\frac{\unicode[STIX]{x1D6FF}\boldsymbol{B}\boldsymbol{\cdot }\boldsymbol{B}}{4\unicode[STIX]{x03C0}\unicode[STIX]{x1D6E5}_{P}}\right)_{S_{P}}. & & \displaystyle\end{eqnarray}$$ See appendix B for details of the calculation. The approach demonstrated in this section could be extended to compute the shape gradients of other figures of merit involving bounce integrals, such as the 
$\unicode[STIX]{x1D6E4}_{c}$ metric for energetic particle confinement (Nemov et al. Reference Nemov, Kasilov, Kernbichler and Leitold2005) or the variation of the parallel adiabatic invariant on a flux surface (Drevlak et al. Reference Drevlak, Geiger, Helander and Turkin2014).
7 Departure from quasisymmetry
Quasisymmetry is desirable as it ensures collisionless confinement of guiding centres. This property follows when the field strength depends on a linear combination of the Boozer angles, 
$B(\unicode[STIX]{x1D713},\unicode[STIX]{x1D703}_{B},\unicode[STIX]{x1D701}_{B})=B(\unicode[STIX]{x1D713},M\unicode[STIX]{x1D703}_{B}-N\unicode[STIX]{x1D701}_{B})$ for fixed integers 
$M$ and 
$N$ (Nührenberg & Zille Reference Nührenberg and Zille1988; Boozer Reference Boozer1995). Several stellarator configurations have been optimized to be close to quasisymmetry (e.g. Reiman et al. Reference Reiman, Fu, Hirshman, Ku, Monticello, Mynick, Redi, Spong, Zarnstorff and Blackwell1999; Drevlak et al. Reference Drevlak, Brochard, Helander, Kisslinger, Mikhailov, Nührenberg, Nührenberg and Turkin2013; Liu et al. Reference Liu, Shimizu, Isobe, Okamura, Nishimura, Suzuki, Xu, Zhang, Liu and Huang2018; Henneberg et al. Reference Henneberg, Drevlak, Nührenberg, Beidler, Turkin, Loizu and Helander2019) by minimizing the amplitude of symmetry-breaking Fourier harmonics of the field strength. We will consider a figure of merit that does not require a Boozer coordinate transformation; instead, we use a general set of magnetic coordinates 
$(\unicode[STIX]{x1D713},\unicode[STIX]{x1D703},\unicode[STIX]{x1D701})$ to define our figure of merit.
In Boozer coordinates (Boozer Reference Boozer1981; Helander Reference Helander2014) (
$\unicode[STIX]{x1D713},\unicode[STIX]{x1D703}_{B},\unicode[STIX]{x1D701}_{B}$) the covariant form for the magnetic field is
$$\begin{eqnarray}\displaystyle \boldsymbol{B}=I(\unicode[STIX]{x1D713})\unicode[STIX]{x1D735}\unicode[STIX]{x1D703}_{B}+G(\unicode[STIX]{x1D713})\unicode[STIX]{x1D735}\unicode[STIX]{x1D701}_{B}+K(\unicode[STIX]{x1D713},\unicode[STIX]{x1D703}_{B},\unicode[STIX]{x1D701}_{B})\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}. & & \displaystyle\end{eqnarray}$$ Here, 
$G(\unicode[STIX]{x1D713})=(2/c)I_{P}(\unicode[STIX]{x1D713})$, where 
$I_{P}(\unicode[STIX]{x1D713})$ is the poloidal current outside the 
$\unicode[STIX]{x1D713}$ surface. The poloidal current can be computed using Ampere’s law and expressed as an integral over a surface labelled by 
$\unicode[STIX]{x1D713}$, 
$S_{P}(\unicode[STIX]{x1D713})$,
$$\begin{eqnarray}\displaystyle I_{P}(\unicode[STIX]{x1D713}) & = & \displaystyle \frac{c}{4\unicode[STIX]{x03C0}}\int _{0}^{2\unicode[STIX]{x03C0}}\text{d}\unicode[STIX]{x1D701}\,\boldsymbol{B}\boldsymbol{\cdot }\frac{\unicode[STIX]{x2202}\boldsymbol{r}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D701}}\nonumber\\ \displaystyle & = & \displaystyle -\frac{c}{8\unicode[STIX]{x03C0}^{2}}\int _{S_{P}(\unicode[STIX]{x1D713})}\text{d}^{2}x\,\boldsymbol{B}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D703}\times \boldsymbol{n}.\end{eqnarray}$$ The quantity 
$I(\unicode[STIX]{x1D713})=(2/c)I_{T}(\unicode[STIX]{x1D713})$, where 
$I_{T}(\unicode[STIX]{x1D713})$ is the toroidal current inside the 
$\unicode[STIX]{x1D713}$ surface (3.4). We quantify the departure from quasisymmetry in the following way,
$$\begin{eqnarray}\displaystyle f_{QS}=\frac{1}{2}\int _{V_{P}}\text{d}^{3}x\,w(\unicode[STIX]{x1D713})(\boldsymbol{B}\times \unicode[STIX]{x1D735}\unicode[STIX]{x1D713}\boldsymbol{\cdot }\unicode[STIX]{x1D735}B-F(\unicode[STIX]{x1D713})\boldsymbol{B}\boldsymbol{\cdot }\unicode[STIX]{x1D735}B)^{2}. & & \displaystyle\end{eqnarray}$$ Here, 
$w(\unicode[STIX]{x1D713})$ is a radial weight function and
$$\begin{eqnarray}\displaystyle F(\unicode[STIX]{x1D713})=\frac{(M/N)G(\unicode[STIX]{x1D713})+I(\unicode[STIX]{x1D713})}{(M/N)\unicode[STIX]{x1D704}(\unicode[STIX]{x1D713})-1}. & & \displaystyle\end{eqnarray}$$ If 
$f_{QS}=0$, then the field is quasisymmetric with mode numbers 
$M$ and 
$N$ (Helander Reference Helander2014), which can be shown using the covariant (3.1) and contravariant (7.1) representations of the magnetic field assuming 
$B=B(\unicode[STIX]{x1D713},M\unicode[STIX]{x1D703}_{B}-N\unicode[STIX]{x1D701}_{B})$ for fixed 
$M$ and 
$N$. Note that 
$f_{QS}$ quantifies the symmetry in Boozer coordinates but can be evaluated in any flux coordinate system.
We consider perturbation about an equilibrium with fixed toroidal current (4.6)–(4.8). The perturbations to the Boozer poloidal covariant component is computed using the transport theorem (2.3),
$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6FF}G(\unicode[STIX]{x1D713})=-\frac{1}{4\unicode[STIX]{x03C0}^{2}}\int _{S_{P}(\unicode[STIX]{x1D713})}\text{d}^{2}x\,(\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\boldsymbol{B}\times \unicode[STIX]{x1D735}\unicode[STIX]{x1D703})\unicode[STIX]{x1D743}_{1}\boldsymbol{\cdot }\boldsymbol{n}+\unicode[STIX]{x1D6FF}\boldsymbol{B}\times \unicode[STIX]{x1D735}\unicode[STIX]{x1D703}\boldsymbol{\cdot }\boldsymbol{n}). & & \displaystyle\end{eqnarray}$$ In arriving at (7.5) we have used the fact that spatial derivatives commute with shape derivatives. The first term accounts for the unperturbed current density through the perturbed boundary, and the second accounts for the perturbed current density through the unperturbed boundary. The contribution from the perturbation to the poloidal angle can be shown to vanish. Upon application of (3.5) we obtain, noting that 
$\int _{S_{P}(\unicode[STIX]{x1D713})}\text{d}^{2}x\,A=V^{\prime }(\unicode[STIX]{x1D713})\langle A|\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}|\rangle _{\unicode[STIX]{x1D713}}$ for any quantity 
$A$,
$$\begin{eqnarray}\displaystyle & & \displaystyle \unicode[STIX]{x1D6FF}G(\unicode[STIX]{x1D713})=-\frac{V^{\prime }(\unicode[STIX]{x1D713})}{4\unicode[STIX]{x03C0}^{2}}\nonumber\\ \displaystyle & & \displaystyle \quad \times \left\langle \unicode[STIX]{x1D743}_{1}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\boldsymbol{B}\times \unicode[STIX]{x1D735}\unicode[STIX]{x1D703})-\frac{1}{\sqrt{g}}\frac{\unicode[STIX]{x2202}\boldsymbol{r}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D701}}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\times (\unicode[STIX]{x1D743}_{1}\times \boldsymbol{B})-\frac{\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D704}_{1}(\unicode[STIX]{x1D713})}{\sqrt{g}^{2}}\frac{\unicode[STIX]{x2202}\boldsymbol{r}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D701}}\boldsymbol{\cdot }\frac{\unicode[STIX]{x2202}\boldsymbol{r}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D703}}\right\rangle _{\unicode[STIX]{x1D713}}.\quad\end{eqnarray}$$ Applying the transport theorem (2.3), the shape derivative of 
$f_{QS}$ takes the form,
$$\begin{eqnarray}\displaystyle & & \displaystyle \unicode[STIX]{x1D6FF}f_{QS}(S_{P};\unicode[STIX]{x1D743}_{1})=\frac{1}{2}\int _{S_{P}}\text{d}^{2}x\,\unicode[STIX]{x1D743}_{1}\boldsymbol{\cdot }\boldsymbol{n}{\mathcal{M}}^{2}w(\unicode[STIX]{x1D713})+\frac{1}{2}\int _{V_{P}}\text{d}^{3}x\,w^{\prime }(\unicode[STIX]{x1D713})\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D713}{\mathcal{M}}^{2}\nonumber\\ \displaystyle & & \displaystyle \quad +\,\int _{V_{P}}\text{d}^{3}x\,w(\unicode[STIX]{x1D713}){\mathcal{M}}\left(\unicode[STIX]{x1D6FF}\boldsymbol{B}\boldsymbol{\cdot }\boldsymbol{{\mathcal{A}}}+\boldsymbol{{\mathcal{S}}}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D6FF}B+\boldsymbol{B}\times \unicode[STIX]{x1D735}\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D713}\boldsymbol{\cdot }\unicode[STIX]{x1D735}B-\frac{\unicode[STIX]{x1D6FF}G(\unicode[STIX]{x1D713})\boldsymbol{B}\boldsymbol{\cdot }\unicode[STIX]{x1D735}B}{\unicode[STIX]{x1D704}(\unicode[STIX]{x1D713})-(N/M)}\right)\nonumber\\ \displaystyle & & \displaystyle \quad +\,\int _{V_{P}}\text{d}^{3}x\,w(\unicode[STIX]{x1D713}){\mathcal{M}}\left(\frac{F(\unicode[STIX]{x1D713})}{\unicode[STIX]{x1D704}(\unicode[STIX]{x1D713})-(N/M)}\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D704}_{1}(\unicode[STIX]{x1D713})\boldsymbol{B}\boldsymbol{\cdot }\unicode[STIX]{x1D735}B-\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D713}F^{\prime }(\unicode[STIX]{x1D713})\boldsymbol{B}\boldsymbol{\cdot }\unicode[STIX]{x1D735}B\right),\end{eqnarray}$$ where 
${\mathcal{M}}=\boldsymbol{B}\times \unicode[STIX]{x1D735}\unicode[STIX]{x1D713}\boldsymbol{\cdot }\unicode[STIX]{x1D735}B-F(\unicode[STIX]{x1D713})\boldsymbol{B}\boldsymbol{\cdot }\unicode[STIX]{x1D735}B$, 
$\boldsymbol{{\mathcal{A}}}=\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}\times \unicode[STIX]{x1D735}B-F(\unicode[STIX]{x1D713})\unicode[STIX]{x1D735}B$ and 
$\boldsymbol{{\mathcal{S}}}=\boldsymbol{B}\times \unicode[STIX]{x1D735}\unicode[STIX]{x1D713}-F(\unicode[STIX]{x1D713})\boldsymbol{B}$. After several steps outlined in appendix C, the shape derivative can be written in the following way,
$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6FF}f_{QS}(S_{P};\unicode[STIX]{x1D743}_{1})=\int _{V_{P}}\text{d}^{3}x\,(\unicode[STIX]{x1D743}_{1}\boldsymbol{\cdot }\boldsymbol{{\mathcal{F}}}_{QS}+\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D704}_{1}(\unicode[STIX]{x1D713}){\mathcal{I}}_{QS})+\int _{S_{P}}\text{d}^{2}x\,\unicode[STIX]{x1D743}_{1}\boldsymbol{\cdot }\boldsymbol{n}{\mathcal{B}}_{QS}, & & \displaystyle\end{eqnarray}$$with
$$\begin{eqnarray}\displaystyle \boldsymbol{{\mathcal{F}}}_{QS} & = & \displaystyle \frac{1}{2}\unicode[STIX]{x1D735}_{\bot }(w(\unicode[STIX]{x1D713}){\mathcal{M}}^{2})+((\boldsymbol{b}\times \unicode[STIX]{x1D735}\unicode[STIX]{x1D713})\unicode[STIX]{x1D735}_{\Vert }B+F(\unicode[STIX]{x1D713})\unicode[STIX]{x1D735}_{\bot }B)w(\unicode[STIX]{x1D713})\boldsymbol{B}\boldsymbol{\cdot }\unicode[STIX]{x1D735}{\mathcal{M}}\nonumber\\ \displaystyle & & \displaystyle +\,\boldsymbol{B}\times (\unicode[STIX]{x1D735}\times (\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}\times \unicode[STIX]{x1D735}B))w(\unicode[STIX]{x1D713}){\mathcal{M}}-B\unicode[STIX]{x1D735}_{\bot }(w(\unicode[STIX]{x1D713})\boldsymbol{{\mathcal{S}}}\boldsymbol{\cdot }\unicode[STIX]{x1D735}{\mathcal{M}})+\unicode[STIX]{x1D73F}Bw(\unicode[STIX]{x1D713})\boldsymbol{{\mathcal{S}}}\boldsymbol{\cdot }\unicode[STIX]{x1D735}{\mathcal{M}}\nonumber\\ \displaystyle & & \displaystyle -\,\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}\unicode[STIX]{x1D735}B\boldsymbol{\cdot }\unicode[STIX]{x1D735}\times (w(\unicode[STIX]{x1D713}){\mathcal{M}}\boldsymbol{B})\nonumber\\ \displaystyle & & \displaystyle +\,\frac{1}{4\unicode[STIX]{x03C0}^{2}}\left(-\unicode[STIX]{x1D735}_{\bot }\left(\frac{w(\unicode[STIX]{x1D713})V^{\prime }(\unicode[STIX]{x1D713})\langle {\mathcal{M}}\boldsymbol{B}\boldsymbol{\cdot }\unicode[STIX]{x1D735}B\rangle _{\unicode[STIX]{x1D713}}}{(\unicode[STIX]{x1D704}(\unicode[STIX]{x1D713})-(N/M))}\right)(\boldsymbol{B}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}\times \unicode[STIX]{x1D735}\unicode[STIX]{x1D703})\right.\nonumber\\ \displaystyle & & \displaystyle +\left.\frac{w(\unicode[STIX]{x1D713})V^{\prime }(\unicode[STIX]{x1D713})\langle {\mathcal{M}}\boldsymbol{B}\boldsymbol{\cdot }\unicode[STIX]{x1D735}B\rangle _{\unicode[STIX]{x1D713}}}{\unicode[STIX]{x1D704}(\unicode[STIX]{x1D713})-(N/M)}(\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\boldsymbol{B}\times \unicode[STIX]{x1D735}\unicode[STIX]{x1D703})-\boldsymbol{B}\times \unicode[STIX]{x1D735}\times (\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}\times \unicode[STIX]{x1D735}\unicode[STIX]{x1D703}))\right)\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle {\mathcal{B}}_{QS} & = & \displaystyle -\frac{1}{2}w(\unicode[STIX]{x1D713}){\mathcal{M}}^{2}+Bw(\unicode[STIX]{x1D713})\boldsymbol{{\mathcal{S}}}\boldsymbol{\cdot }\unicode[STIX]{x1D735}{\mathcal{M}}-w(\unicode[STIX]{x1D713}){\mathcal{M}}\unicode[STIX]{x1D735}B\times \boldsymbol{B}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}\nonumber\\ \displaystyle & & \displaystyle +\,\frac{w(\unicode[STIX]{x1D713})V^{\prime }(\unicode[STIX]{x1D713})\langle {\mathcal{M}}\boldsymbol{B}\boldsymbol{\cdot }\unicode[STIX]{x1D735}B\rangle _{\unicode[STIX]{x1D713}}}{4\unicode[STIX]{x03C0}^{2}(\unicode[STIX]{x1D704}(\unicode[STIX]{x1D713})-(N/M))}(\boldsymbol{B}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}\times \unicode[STIX]{x1D735}\unicode[STIX]{x1D703})\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle {\mathcal{I}}_{QS} & = & \displaystyle -w(\unicode[STIX]{x1D713}){\mathcal{M}}\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}\times \unicode[STIX]{x1D735}\unicode[STIX]{x1D701}\boldsymbol{\cdot }\boldsymbol{{\mathcal{A}}}+w(\unicode[STIX]{x1D713})(\boldsymbol{{\mathcal{S}}}\boldsymbol{\cdot }\unicode[STIX]{x1D735}{\mathcal{M}})\boldsymbol{b}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}\times \unicode[STIX]{x1D735}\unicode[STIX]{x1D701}\nonumber\\ \displaystyle & & \displaystyle +\,\frac{w(\unicode[STIX]{x1D713}){\mathcal{M}}\boldsymbol{B}\boldsymbol{\cdot }\unicode[STIX]{x1D735}B}{\unicode[STIX]{x1D704}(\unicode[STIX]{x1D713})-(N/M)}\left(F(\unicode[STIX]{x1D713})-\left\langle \frac{V^{\prime }(\unicode[STIX]{x1D713})}{4\unicode[STIX]{x03C0}^{2}\sqrt{g}^{2}}\frac{\unicode[STIX]{x2202}\boldsymbol{r}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D701}}\boldsymbol{\cdot }\frac{\unicode[STIX]{x2202}\boldsymbol{r}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D703}}\right\rangle _{\unicode[STIX]{x1D713}}\right).\end{eqnarray}$$ In (7.9), 
$\unicode[STIX]{x1D735}_{\Vert }=\boldsymbol{b}\boldsymbol{\cdot }\unicode[STIX]{x1D735}$ and 
$\unicode[STIX]{x1D735}_{\bot }=\unicode[STIX]{x1D735}-\boldsymbol{b}\unicode[STIX]{x1D735}_{\Vert }$ are the parallel and perpendicular gradients. We note that 
$\boldsymbol{{\mathcal{F}}}_{QS}$ satisfies the parallel force balance condition (
$\boldsymbol{b}\boldsymbol{\cdot }\boldsymbol{{\mathcal{F}}}_{QS}=0$) implied by (3.11).
We can now prescribe an adjoint perturbation which satisfies,
$$\begin{eqnarray}\displaystyle \displaystyle \unicode[STIX]{x1D6FF}\boldsymbol{F}_{2} & = & \displaystyle \unicode[STIX]{x1D6E5}_{QS}\boldsymbol{{\mathcal{F}}}_{QS}\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle \displaystyle \unicode[STIX]{x1D743}\boldsymbol{\cdot }\boldsymbol{n}|_{S_{P}} & = & \displaystyle 0\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle \displaystyle \unicode[STIX]{x1D6FF}I_{T,2}(\unicode[STIX]{x1D713}) & = & \displaystyle \frac{c\unicode[STIX]{x1D6E5}_{QS}}{2\unicode[STIX]{x03C0}}V^{\prime }(\unicode[STIX]{x1D713})\langle {\mathcal{I}}_{QS}\rangle _{\unicode[STIX]{x1D713}}.\end{eqnarray}$$Upon application of the fixed-boundary adjoint relation we obtain the following shape gradient,
$$\begin{eqnarray}\displaystyle {\mathcal{G}}_{QS}=\left(\frac{\unicode[STIX]{x1D6FF}\boldsymbol{B}_{2}\boldsymbol{\cdot }\boldsymbol{B}}{4\unicode[STIX]{x03C0}\unicode[STIX]{x1D6E5}_{QS}}+{\mathcal{B}}_{QS}\right)_{S_{P}}. & & \displaystyle\end{eqnarray}$$8 Neoclassical figures of merit
In § 6, we considered a figure of merit that quantifies the geometric dependence of the neoclassical particle flux in the 
$1/\unicode[STIX]{x1D708}$ regime. In applying this model, several assumptions are imposed, such as a small radial electric field, 
$E_{r}$, low collisionality and a simplified pitch-angle scattering collision operator. In this section, we consider a more general neoclassical figure of merit arising from a moment of the local drift kinetic equation, allowing for optimization at finite collisionality and 
$E_{r}$. It is assumed here that the collision time is comparable to the bounce time but shorter than the time needed to complete a magnetic drift orbit. Recently an adjoint method has been demonstrated for obtaining derivatives of neoclassical figures of merit (Paul et al. Reference Paul, Abel, Landreman and Dorland2019) with respect to local geometric quantities on a flux surface. The adjoint method described in this section will extend these results, such that shape derivatives with respect to the plasma boundary can be computed.
Consider the following figure of merit,
$$\begin{eqnarray}\displaystyle f_{NC}=\int _{V_{P}}\text{d}^{3}x\,w(\unicode[STIX]{x1D713}){\mathcal{R}}(\unicode[STIX]{x1D713}). & & \displaystyle\end{eqnarray}$$ Here 
${\mathcal{R}}(\unicode[STIX]{x1D713})$ is a flux surface-averaged moment of the neoclassical distribution function, 
$f_{1}$, which satisfies the local drift kinetic equation (DKE),
$$\begin{eqnarray}\displaystyle (v_{\Vert }\boldsymbol{b}+\boldsymbol{v}_{E})\boldsymbol{\cdot }\unicode[STIX]{x1D735}f_{1}-C(f_{1})=-\boldsymbol{v}_{m}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}\frac{\unicode[STIX]{x2202}f_{0}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}}, & & \displaystyle\end{eqnarray}$$ where 
$\boldsymbol{v}_{E}=\boldsymbol{E}\times \boldsymbol{B}/B^{2}$ is the 
$\boldsymbol{E}\times \boldsymbol{B}$ drift velocity, 
$\boldsymbol{v}_{m}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}$ is the radial magnetic drift velocity (B 4), 
$f_{0}$ is a Maxwellian (B 3) and 
$C$ is the linearized Fokker-Planck operator. For example, 
${\mathcal{R}}$ can be taken to be the bootstrap current,
$$\begin{eqnarray}\displaystyle J_{b}=\mathop{\sum }_{s}\frac{\displaystyle \left\langle B\int \text{d}^{3}v\,f_{1s}v_{\Vert }\right\rangle _{\unicode[STIX]{x1D713}}}{n_{s}\langle B^{2}\rangle _{\unicode[STIX]{x1D713}}^{1/2}}, & & \displaystyle\end{eqnarray}$$ where the sum is taken over species. We note that the geometric dependence that enters the DKE when written in Boozer coordinates only arises through the quantities 
$\{B,G(\unicode[STIX]{x1D713}),I(\unicode[STIX]{x1D713}),\unicode[STIX]{x1D704}(\unicode[STIX]{x1D713})\}$. Thus for simplicity, Boozer coordinates will be assumed throughout this section.
The perturbation to 
${\mathcal{R}}(\unicode[STIX]{x1D713})$ at fixed toroidal current (4.6)–(4.8) can be written as,
$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6FF}{\mathcal{R}}(\unicode[STIX]{x1D713})=\langle S_{{\mathcal{R}}}\unicode[STIX]{x1D6FF}B\rangle _{\unicode[STIX]{x1D713}}+\frac{\unicode[STIX]{x2202}{\mathcal{R}}(\unicode[STIX]{x1D713})}{\unicode[STIX]{x2202}G(\unicode[STIX]{x1D713})}\unicode[STIX]{x1D6FF}G(\unicode[STIX]{x1D713})+\frac{\unicode[STIX]{x2202}{\mathcal{R}}(\unicode[STIX]{x1D713})}{\unicode[STIX]{x2202}\unicode[STIX]{x1D704}(\unicode[STIX]{x1D713})}\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D704}_{1}(\unicode[STIX]{x1D713}). & & \displaystyle\end{eqnarray}$$ Here, 
$S_{{\mathcal{R}}}$ is a local sensitivity function which quantifies the change to 
${\mathcal{R}}$ associated with a perturbation of the field strength 
$\unicode[STIX]{x1D6FF}B$ defined in the following way. Consider the perturbation to 
${\mathcal{R}}$ resulting from a change in the field strength at fixed 
$G(\unicode[STIX]{x1D713})$, 
$I(\unicode[STIX]{x1D713})$, and 
$\unicode[STIX]{x1D704}(\unicode[STIX]{x1D713})$. The functional derivative of 
${\mathcal{R}}(\unicode[STIX]{x1D713})$ with respect to 
$B(\boldsymbol{r})$ can be expressed as,
$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6FF}{\mathcal{R}}(\unicode[STIX]{x1D6FF}B;B(\boldsymbol{r}))=\langle S_{{\mathcal{R}}}\unicode[STIX]{x1D6FF}B(\boldsymbol{r})\rangle _{\unicode[STIX]{x1D713}}. & & \displaystyle\end{eqnarray}$$ This is another instance of the Riesz representation theorem: 
$\unicode[STIX]{x1D6FF}{\mathcal{R}}$ is a linear functional of 
$\unicode[STIX]{x1D6FF}B$, with the inner product taken to be the flux surface average. Thus 
$S_{{\mathcal{R}}}$ can be thought of as analogous to the shape gradient (2.4).
The quantities 
$\{S_{{\mathcal{R}}},\unicode[STIX]{x2202}{\mathcal{R}}(\unicode[STIX]{x1D713})/\unicode[STIX]{x2202}G(\unicode[STIX]{x1D713}),\unicode[STIX]{x2202}{\mathcal{R}}(\unicode[STIX]{x1D713})/\unicode[STIX]{x2202}\unicode[STIX]{x1D704}(\unicode[STIX]{x1D713})\}$ can be computed with a related adjoint method (Paul et al. Reference Paul, Abel, Landreman and Dorland2019) with the SFINCS code (Landreman et al. Reference Landreman, Smith, Mollén and Helander2014). Here we consider SFINCS to be run on a set of surfaces such that (8.1) can be computed numerically. The derivatives computed by SFINCS will appear in the additional bulk force required for the adjoint perturbed equilibrium. The shape derivative of 
$f_{NC}$ can be computed on application of the transport theorem (2.3),
$$\begin{eqnarray}\displaystyle & & \displaystyle \unicode[STIX]{x1D6FF}f_{NC}(S_{P};\unicode[STIX]{x1D743}_{1})=\int _{S_{P}}\text{d}^{2}x\,\unicode[STIX]{x1D743}_{1}\boldsymbol{\cdot }\boldsymbol{n}w(\unicode[STIX]{x1D713}){\mathcal{R}}(\unicode[STIX]{x1D713})+\int _{V_{P}}\text{d}^{3}x\,\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D713}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}}(w(\unicode[STIX]{x1D713}){\mathcal{R}}(\unicode[STIX]{x1D713}))\nonumber\\ \displaystyle & & \displaystyle \quad +\,\int _{V_{P}}\text{d}^{3}x\,w(\unicode[STIX]{x1D713})\left(\frac{\unicode[STIX]{x2202}{\mathcal{R}}(\unicode[STIX]{x1D713})}{\unicode[STIX]{x2202}G(\unicode[STIX]{x1D713})}\unicode[STIX]{x1D6FF}G(\unicode[STIX]{x1D713})+\frac{\unicode[STIX]{x2202}{\mathcal{R}}(\unicode[STIX]{x1D713})}{\unicode[STIX]{x2202}\unicode[STIX]{x1D704}(\unicode[STIX]{x1D713})}\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D704}_{1}(\unicode[STIX]{x1D713})+\langle S_{R}\unicode[STIX]{x1D6FF}B\rangle _{\unicode[STIX]{x1D713}}\right).\end{eqnarray}$$After several steps outlined in appendix D, the shape derivative is written in the following form,
$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6FF}f_{NC}(S_{P};\unicode[STIX]{x1D743}_{1})=\int _{V_{P}}\text{d}^{3}x\,(\unicode[STIX]{x1D743}_{1}\boldsymbol{\cdot }\boldsymbol{{\mathcal{F}}}_{NC}+\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D704}_{1}(\unicode[STIX]{x1D713}){\mathcal{I}}_{NC})+\int _{S_{P}}\text{d}^{3}x\,\unicode[STIX]{x1D743}_{1}\boldsymbol{\cdot }\boldsymbol{n}{\mathcal{B}}_{NC} & & \displaystyle\end{eqnarray}$$with
$$\begin{eqnarray}\displaystyle \boldsymbol{{\mathcal{F}}}_{NC} & = & \displaystyle -\unicode[STIX]{x1D735}({\mathcal{R}}(\unicode[STIX]{x1D713})w(\unicode[STIX]{x1D713}))-\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}(\unicode[STIX]{x1D735}\times \boldsymbol{B})\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D703}\frac{\unicode[STIX]{x2202}{\mathcal{R}}(\unicode[STIX]{x1D713})}{\unicode[STIX]{x2202}G(\unicode[STIX]{x1D713})}w(\unicode[STIX]{x1D713})\frac{B^{2}\sqrt{g}}{\langle B^{2}\rangle _{\unicode[STIX]{x1D713}}}\nonumber\\ \displaystyle & & \displaystyle +\,\frac{w(\unicode[STIX]{x1D713})}{\langle B^{2}\rangle _{\unicode[STIX]{x1D713}}}\frac{\unicode[STIX]{x2202}{\mathcal{R}}(\unicode[STIX]{x1D713})}{\unicode[STIX]{x2202}G(\unicode[STIX]{x1D713})}\boldsymbol{B}\times \unicode[STIX]{x1D735}\times \left(\frac{\unicode[STIX]{x2202}\boldsymbol{r}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D701}}B^{2}\right)+G(\unicode[STIX]{x1D713})B^{2}\unicode[STIX]{x1D735}\left(\frac{w(\unicode[STIX]{x1D713})}{\langle B^{2}\rangle _{\unicode[STIX]{x1D713}}}\frac{\unicode[STIX]{x2202}{\mathcal{R}}(\unicode[STIX]{x1D713})}{\unicode[STIX]{x2202}G(\unicode[STIX]{x1D713})}\right)\nonumber\\ \displaystyle & & \displaystyle -\,\unicode[STIX]{x1D73F}w(\unicode[STIX]{x1D713})S_{{\mathcal{R}}}B+B\unicode[STIX]{x1D735}_{\bot }(w(\unicode[STIX]{x1D713})S_{{\mathcal{R}}})\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle \displaystyle {\mathcal{B}}_{NC} & = & \displaystyle w(\unicode[STIX]{x1D713}){\mathcal{R}}(\unicode[STIX]{x1D713})-\frac{w(\unicode[STIX]{x1D713})B^{2}}{\langle B^{2}\rangle _{\unicode[STIX]{x1D713}}}\frac{\unicode[STIX]{x2202}{\mathcal{R}}(\unicode[STIX]{x1D713})}{\unicode[STIX]{x2202}G(\unicode[STIX]{x1D713})}G(\unicode[STIX]{x1D713})-w(\unicode[STIX]{x1D713})S_{{\mathcal{R}}}B\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle \displaystyle {\mathcal{I}}_{NC} & = & \displaystyle \frac{\unicode[STIX]{x2202}{\mathcal{R}}(\unicode[STIX]{x1D713})}{\unicode[STIX]{x2202}G(\unicode[STIX]{x1D713})}\frac{w(\unicode[STIX]{x1D713})B^{2}}{\langle B^{2}\rangle _{\unicode[STIX]{x1D713}}\sqrt{g}}\frac{\unicode[STIX]{x2202}\boldsymbol{r}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D701}}\boldsymbol{\cdot }\frac{\unicode[STIX]{x2202}\boldsymbol{r}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D703}}+w(\unicode[STIX]{x1D713})\frac{\unicode[STIX]{x2202}{\mathcal{R}}(\unicode[STIX]{x1D713})}{\unicode[STIX]{x2202}\unicode[STIX]{x1D704}(\unicode[STIX]{x1D713})}-w(\unicode[STIX]{x1D713})S_{{\mathcal{R}}}\boldsymbol{b}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}\times \unicode[STIX]{x1D735}\unicode[STIX]{x1D701}.\quad\end{eqnarray}$$ The adjoint bulk force 
$\boldsymbol{{\mathcal{F}}}_{NC}$ is chosen to satisfy parallel force balance required by (3.11).
We consider the following adjoint perturbation,
$$\begin{eqnarray}\displaystyle \displaystyle \unicode[STIX]{x1D6FF}\boldsymbol{F}_{2} & = & \displaystyle \unicode[STIX]{x1D6E5}_{NC}\boldsymbol{{\mathcal{F}}}_{NC}\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle \displaystyle \unicode[STIX]{x1D743}_{2}\boldsymbol{\cdot }\boldsymbol{n}|_{S_{P}} & = & \displaystyle 0\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle \displaystyle \unicode[STIX]{x1D6FF}I_{T,2}(\unicode[STIX]{x1D713}) & = & \displaystyle \frac{c\unicode[STIX]{x1D6E5}_{NC}}{2\unicode[STIX]{x03C0}}V^{\prime }(\unicode[STIX]{x1D713})\langle {\mathcal{I}}_{NC}\rangle _{\unicode[STIX]{x1D713}}.\end{eqnarray}$$Upon application of the fixed-boundary adjoint relation we obtain the shape gradient,
$$\begin{eqnarray}\displaystyle {\mathcal{G}}_{NC}=\left({\mathcal{B}}_{NC}+\frac{\unicode[STIX]{x1D6FF}\boldsymbol{B}_{2}\boldsymbol{\cdot }\boldsymbol{B}}{4\unicode[STIX]{x03C0}\unicode[STIX]{x1D6E5}_{NC}}\right)_{S_{P}}. & & \displaystyle\end{eqnarray}$$9 Conclusions
We have demonstrated that the self-adjointness relations (§ 3) can be implemented to efficiently compute the shape gradient of figures of merit relevant for stellarator configuration optimization. The shape gradient is obtained by solving an adjoint perturbed force balance equation that depends on the figure of merit of interest. For the vacuum well parameter (§ 4), the additional bulk force required for the adjoint problem is simply the gradient of a function of flux, and so can be implemented by adding a perturbation to the pressure profile. For the magnetic ripple on axis (§ 5), the required bulk force takes the form of the divergence of a pressure tensor that only varies on a surface through the field strength. As this type of pressure tensor is currently treated by the ANIMEC code, this adjoint bulk force is implemented with a minor modification to the code. Computing the shape gradient of 
$\unicode[STIX]{x1D716}_{\text{eff}}^{3/2}$ with the adjoint approach also requires the addition of the divergence of a pressure tensor. However, this pressure tensor varies on a surface through the field line label due to the bounce integrals that appear (6.8)–(6.7). Thus the variational principle used by the ANIMEC code cannot be easily extended for this application. Similarly, the shape gradients for the quasisymmetry (§ 7) and neoclassical (§ 8) figures of merit require an adjoint bulk force that is not in the form of the divergence of a pressure tensor. This provides an impetus for the development of a flexible perturbed MHD equilibrium code that could enable these calculations. While several three-dimensional ideal MHD stability codes exist (Anderson et al. Reference Anderson, Cooper, Gruber, Merazzi and Schwenn1990; Schwab Reference Schwab1993; Strumberger & Günter Reference Strumberger and Günter2016), only the CAS3D code has been modified in order to perform perturbed equilibrium calculations (Nührenberg & Boozer Reference Nührenberg and Boozer2003; Boozer & Nührenberg Reference Boozer and Nührenberg2006). We hope to take advantage of linear MHD calculations for future work through the proper modification of an existing code or development of a new code.
To date, the numerical verification of this adjoint approach for MHD equilibria has been performed with the VMEC and ANIMEC codes, which solve the nonlinear force balance equations, equations (3.2) and (5.16). The adjoint perturbed force balance equation (3.11) is approximated by adding a perturbative bulk force or toroidal current profile, whose characteristic magnitudes are scaled by the 
$\unicode[STIX]{x1D6E5}$ constants (e.g. 
$\unicode[STIX]{x1D6E5}_{P}$ in (4.10)). As demonstrated in Antonsen et al. (Reference Antonsen, Paul and Landreman2019), these parameters must be small enough that nonlinear effects do not become important yet large enough that round-off error does not dominate. This furthermore motivates the development of a perturbed equilibrium code that could eliminate this source of noise.
As demonstrated, this adjoint approach for functions of MHD equilibria is quite flexible and can be applied to many quantities of interest. Because of the demonstrated efficiency in comparison with the direct approach to computing shape gradients, we anticipate many further applications of this method.
Acknowledgements
The authors would like to acknowledge productive discussions with R. Nochetto on shape calculus. This work was supported by the US Department of Energy through grants DE-FG02-93ER-54197 and DE-FC02-08ER-54964. The computations presented in this paper have used resources at the National Energy Research Scientific Computing Center (NERSC).
Appendix A. Details of axis ripple calculation
In this appendix we compute the shape derivative of the finite pressure magnetic well figure of merit from (5.8) and show that if we impose an adjoint perturbation of the form (5.9)–(5.11), the shape gradient is given by (5.14) with (5.9)–(5.11).
We use the expression for the perturbation to the field strength (5.7) and 
$\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D713}=-\unicode[STIX]{x1D743}_{1}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}$ with (5.8) to obtain
$$\begin{eqnarray}\displaystyle & & \displaystyle \unicode[STIX]{x1D6FF}f_{R}(S_{P};\unicode[STIX]{x1D743}_{1})=\int _{S_{P}}\text{d}^{2}x\,\unicode[STIX]{x1D743}_{1}\boldsymbol{\cdot }\boldsymbol{n}\widetilde{f_{R}}-\int _{V_{P}}\text{d}^{3}x\,\frac{\unicode[STIX]{x2202}\widetilde{f_{R}}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}}\unicode[STIX]{x1D743}_{1}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}\nonumber\\ \displaystyle & & \displaystyle \quad -\,\int _{V_{P}}\text{d}^{3}x\,\frac{\unicode[STIX]{x2202}\widetilde{f_{R}}}{\unicode[STIX]{x2202}B}\frac{1}{B}(B^{2}\unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D743}_{1}+\unicode[STIX]{x1D743}_{1}\boldsymbol{\cdot }\unicode[STIX]{x1D735}(B^{2}+4\unicode[STIX]{x03C0}p)+\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D704}_{1}(\unicode[STIX]{x1D713})\boldsymbol{B}\boldsymbol{\cdot }(\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}\times \unicode[STIX]{x1D735}\unicode[STIX]{x1D701})).\quad\end{eqnarray}$$The third term can be integrated by parts to obtain
$$\begin{eqnarray}\displaystyle & & \displaystyle \unicode[STIX]{x1D6FF}f_{R}(S_{P};\unicode[STIX]{x1D743}_{1})=\int _{S_{P}}\text{d}^{2}x\,\unicode[STIX]{x1D743}_{1}\boldsymbol{\cdot }\boldsymbol{n}\left(\widetilde{f_{R}}-\frac{\unicode[STIX]{x2202}\widetilde{f_{R}}}{\unicode[STIX]{x2202}B}B\right)+\int _{V_{P}}\text{d}^{3}x\,\left(\frac{\unicode[STIX]{x2202}^{2}\widetilde{f_{R}}}{\unicode[STIX]{x2202}B\unicode[STIX]{x2202}\unicode[STIX]{x1D713}}B-\frac{\unicode[STIX]{x2202}\widetilde{f_{R}}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}}\right)\unicode[STIX]{x1D743}_{1}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}\nonumber\\ \displaystyle & & \displaystyle \quad +\,\int _{V_{P}}\text{d}^{3}x\,\left(-\frac{\unicode[STIX]{x2202}\widetilde{f_{R}}}{\unicode[STIX]{x2202}B}B\unicode[STIX]{x1D743}_{1}\boldsymbol{\cdot }\unicode[STIX]{x1D73F}+B\frac{\unicode[STIX]{x2202}^{2}\widetilde{f_{R}}}{\unicode[STIX]{x2202}B^{2}}\unicode[STIX]{x1D743}_{1}\boldsymbol{\cdot }\unicode[STIX]{x1D735}B+\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D704}_{1}(\unicode[STIX]{x1D713})\frac{\unicode[STIX]{x2202}\widetilde{f_{R}}}{\unicode[STIX]{x2202}B}\boldsymbol{b}\boldsymbol{\cdot }(\unicode[STIX]{x1D735}\unicode[STIX]{x1D701}\times \unicode[STIX]{x1D735}\unicode[STIX]{x1D713})\right),\end{eqnarray}$$where the expression for the curvature in an equilibrium field has been applied.
We compute one term that appears in the fixed-boundary adjoint relation (3.13) using the prescribed adjoint bulk force perturbation (5.9)
$$\begin{eqnarray}\displaystyle \frac{1}{\unicode[STIX]{x1D6E5}_{P}}\int _{V_{P}}\text{d}^{3}x\,\unicode[STIX]{x1D743}_{1}\boldsymbol{\cdot }\unicode[STIX]{x1D6FF}\boldsymbol{F}_{2} & = & \displaystyle \int _{V_{P}}\text{d}^{3}x\,\left(\frac{\unicode[STIX]{x2202}^{2}p_{\Vert }}{\unicode[STIX]{x2202}B\unicode[STIX]{x2202}\unicode[STIX]{x1D713}}B-\frac{\unicode[STIX]{x2202}p_{\Vert }}{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}}\right)\unicode[STIX]{x1D743}_{1}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}\nonumber\\ \displaystyle & & \displaystyle +\,\int _{V_{P}}\text{d}^{3}x\,\left(-\frac{\unicode[STIX]{x2202}p_{\Vert }}{\unicode[STIX]{x2202}B}B\unicode[STIX]{x1D743}_{1}\boldsymbol{\cdot }\unicode[STIX]{x1D73F}+B\frac{\unicode[STIX]{x2202}^{2}p_{\Vert }}{\unicode[STIX]{x2202}B^{2}}\unicode[STIX]{x1D743}_{1}\boldsymbol{\cdot }\unicode[STIX]{x1D735}B\right),\end{eqnarray}$$ where we have applied the parallel force balance condition (5.12). Therefore, if we impose 
$p_{\Vert }=\widetilde{f_{R}}$, we obtain the following expression for the shape derivative of 
$f_{R}$,
$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6FF}f_{R}(S_{P};\unicode[STIX]{x1D743}_{1}) & = & \displaystyle \int _{S_{P}}\text{d}^{2}x\,\unicode[STIX]{x1D743}_{1}\boldsymbol{\cdot }\boldsymbol{n}\left(\widetilde{f_{R}}-\frac{\unicode[STIX]{x2202}\widetilde{f_{R}}}{\unicode[STIX]{x2202}B}B\right)+\frac{1}{\unicode[STIX]{x1D6E5}_{P}}\int _{V_{P}}\text{d}^{3}x\,\unicode[STIX]{x1D743}_{1}\boldsymbol{\cdot }\unicode[STIX]{x1D6FF}\boldsymbol{F}_{2}\nonumber\\ \displaystyle & & \displaystyle +\,\int _{V_{P}}\text{d}^{3}x\,\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D704}_{1}(\unicode[STIX]{x1D713})\frac{\unicode[STIX]{x2202}\widetilde{f_{R}}}{\unicode[STIX]{x2202}B}\boldsymbol{b}\boldsymbol{\cdot }(\unicode[STIX]{x1D735}\unicode[STIX]{x1D701}\times \unicode[STIX]{x1D735}\unicode[STIX]{x1D713}).\end{eqnarray}$$Upon application of the fixed-boundary adjoint relation we obtain (5.14) with (5.9)–(5.11).
Appendix B. Details of effective ripple in the $1/\unicode[STIX]{x1D708}$ regime calculation
Neoclassical transport in the 
$1/\unicode[STIX]{x1D708}$ collisionality regime is discussed in many references including Frieman (Reference Frieman1970), Connor & Hastie (Reference Connor and Hastie1974) and Ho & Kulsrud (Reference Ho and Kulsrud1987). In this appendix we sketch the computation of 
$\unicode[STIX]{x1D716}_{\text{eff}}^{3/2}$ originally introduced in Nemov et al. (Reference Nemov, Kasilov, Kernbichler and Heyn1999) and compute linear perturbations of 
$f_{\unicode[STIX]{x1D716}}$ (6.4), showing them to take the form of (6.5).
In the 
$1/\unicode[STIX]{x1D708}$ regime, the distribution function is ordered in the parameter 
$\unicode[STIX]{x1D708}_{\ast }=\unicode[STIX]{x1D708}/(v_{t}/L)\ll 1$, where 
$\unicode[STIX]{x1D708}$ 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,
$$\begin{eqnarray}\displaystyle f_{1}=f_{1}^{-1}+f_{1}^{0}+O(\unicode[STIX]{x1D708}_{\ast }). & & \displaystyle\end{eqnarray}$$ In velocity space we use a pitch-angle coordinate 
$\unicode[STIX]{x1D706}=v_{\bot }^{2}/(v^{2}B)$, energy coordinate 
$\unicode[STIX]{x1D716}=v^{2}/2$ and 
$\unicode[STIX]{x1D70E}=\text{sign}(v_{\Vert })$, where 
$v_{\bot }=\sqrt{v^{2}-v_{\Vert }^{2}}$ is the perpendicular velocity and 
$v_{\Vert }=\boldsymbol{v}\boldsymbol{\cdot }\boldsymbol{b}$ is the parallel velocity. We use the field line label, 
$\unicode[STIX]{x1D6FC}$, and length along a field line, 
$l$, to describe location on a constant 
$\unicode[STIX]{x1D713}$ surface. In the 
$1/\unicode[STIX]{x1D708}$ regime the 
$\boldsymbol{E}\times \boldsymbol{B}$ precession frequency is assumed to be small relative to the collision frequency, and the drift kinetic equation is,
$$\begin{eqnarray}\displaystyle v_{\Vert }\frac{\unicode[STIX]{x2202}f_{1}}{\unicode[STIX]{x2202}l}=C(f_{1})-\boldsymbol{v}_{m}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}\frac{\unicode[STIX]{x2202}f_{0}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}}, & & \displaystyle\end{eqnarray}$$ where the Maxwellian with density 
$n$ is
$$\begin{eqnarray}\displaystyle f_{0}=n\unicode[STIX]{x03C0}^{-3/2}v_{t}^{-3}e^{-v^{2}/v_{t}^{2}} & & \displaystyle\end{eqnarray}$$and the radial magnetic drift is
$$\begin{eqnarray}\displaystyle \boldsymbol{v}_{\text{m}}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}=(v^{2}+v_{\Vert }^{2})\frac{mc}{2qB^{3}}\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}\times \boldsymbol{B}\boldsymbol{\cdot }\unicode[STIX]{x1D735}B, & & \displaystyle\end{eqnarray}$$ for charge 
$q$. The drift kinetic equation to 
$O(\unicode[STIX]{x1D708}_{\ast }^{-1})$ is
$$\begin{eqnarray}\displaystyle v_{\Vert }\frac{\unicode[STIX]{x2202}f_{1}^{-1}}{\unicode[STIX]{x2202}l}=0. & & \displaystyle\end{eqnarray}$$ In the trapped portion of phase space, this implies that 
$f_{1}^{-1}=f_{1}^{-1}(\unicode[STIX]{x1D713},\unicode[STIX]{x1D6FC},\unicode[STIX]{x1D716},\unicode[STIX]{x1D706})$, and in the passing portion of phase space, this implies that 
$f_{1}^{-1}=f_{1}^{-1}(\unicode[STIX]{x1D713},\unicode[STIX]{x1D716},\unicode[STIX]{x1D706},\unicode[STIX]{x1D70E})$. The drift kinetic equation to 
$O(\unicode[STIX]{x1D708}_{\ast }^{0})$ is
$$\begin{eqnarray}\displaystyle v_{\Vert }\frac{\unicode[STIX]{x2202}f_{1}^{0}}{\unicode[STIX]{x2202}l}=C(f_{1}^{-1})-\boldsymbol{v}_{\text{m}}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}\frac{\unicode[STIX]{x2202}f_{0}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}}. & & \displaystyle\end{eqnarray}$$ 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,
$$\begin{eqnarray}\displaystyle C=\frac{2\unicode[STIX]{x1D708}(\unicode[STIX]{x1D716})v_{\Vert }}{B\unicode[STIX]{x1D716}}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D706}}\left(\unicode[STIX]{x1D706}v_{\Vert }\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D706}}\right). & & \displaystyle\end{eqnarray}$$The parallel streaming term in (B 6) is annihilated by the bounce averaging operation,
$$\begin{eqnarray}\displaystyle 0=\langle C(f_{1}^{-1})\rangle _{b}-\langle \boldsymbol{v}_{\text{m}}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}\rangle _{b}\frac{\unicode[STIX]{x2202}f_{0}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}}, & & \displaystyle\end{eqnarray}$$ where the bounce average of a quantity 
$A$ is 
$\langle A\rangle _{b}=\unicode[STIX]{x1D70F}^{-1}\oint \text{d}l\,A/v_{\Vert }$ and the bounce time is 
$\unicode[STIX]{x1D70F}=\oint \text{d}l\,v_{\Vert }^{-1}$. The bounce-averaged equation (B 8) can be expressed in terms of the parallel adiabatic invariant 
$J=\oint \text{d}l\,v_{\Vert }$ using the relation
$$\begin{eqnarray}\displaystyle \langle \boldsymbol{v}_{m}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}\rangle _{b}=\frac{mc}{q\unicode[STIX]{x1D70F}}\frac{\unicode[STIX]{x2202}J}{\unicode[STIX]{x2202}\unicode[STIX]{x1D6FC}}. & & \displaystyle\end{eqnarray}$$ Integrating (B 8) with respect to 
$\unicode[STIX]{x1D706}$ we obtain,
$$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x2202}f_{1}^{-1}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D706}}=\frac{mc\unicode[STIX]{x1D716}}{2q\unicode[STIX]{x1D706}\unicode[STIX]{x1D708}(\unicode[STIX]{x1D716})}\frac{\unicode[STIX]{x2202}f_{0}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}}\left(\oint \text{d}l\,\frac{v_{\Vert }}{B}\right)^{-1}\int _{1/B_{\max }}^{\unicode[STIX]{x1D706}}\text{d}\unicode[STIX]{x1D706}^{\prime }\,\frac{\unicode[STIX]{x2202}J}{\unicode[STIX]{x2202}\unicode[STIX]{x1D6FC}}. & & \displaystyle\end{eqnarray}$$ Here, 
$B_{\max }$ is the maximum value of the field strength on the surface labelled by 
$\unicode[STIX]{x1D713}$. We have used the boundary condition 
$(\oint \text{d}l\,v_{\Vert }/B)\unicode[STIX]{x2202}f_{1}^{-1}/\unicode[STIX]{x2202}\unicode[STIX]{x1D706}|_{\unicode[STIX]{x1D706}=1/B_{\max }}=0$, as there is no flux in pitch angle from the passing region. The integration with respect to 
$\unicode[STIX]{x1D706}$ is performed to obtain,
$$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x2202}f_{1}^{-1}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D706}}=-\frac{mc}{6q\unicode[STIX]{x1D706}\unicode[STIX]{x1D708}(\unicode[STIX]{x1D716})}\frac{\unicode[STIX]{x2202}f_{0}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}}\left(\oint \text{d}l\,\frac{v_{\Vert }}{B}\right)^{-1}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D6FC}}\left(\oint \text{d}l\,\frac{\mathop{v}_{\Vert }^{3}}{B}\right). & & \displaystyle\end{eqnarray}$$ The particle flux from 
$f_{1}^{-1}$ is obtained by multiplying (B 6) by 
$f_{1}^{-1}(\unicode[STIX]{x2202}f_{0}/\unicode[STIX]{x2202}\unicode[STIX]{x1D713})^{-1}$, integrating over velocity space and flux surface averaging,
$$\begin{eqnarray}\displaystyle \langle \unicode[STIX]{x1D71E}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}\rangle _{\unicode[STIX]{x1D713}}\equiv \left\langle \int \text{d}^{3}v\,f_{1}^{-1}\boldsymbol{v}_{\text{m}}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}\right\rangle _{\unicode[STIX]{x1D713}}=\left\langle \int \text{d}^{3}v\,f_{1}^{-1}C(f_{1}^{-1})\left(\frac{\unicode[STIX]{x2202}f_{0}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}}\right)^{-1}\right\rangle _{\unicode[STIX]{x1D713}}. & & \displaystyle\end{eqnarray}$$ The velocity space integration is performed using the Jacobian 
$\text{d}^{3}v=2\unicode[STIX]{x03C0}\sum _{\unicode[STIX]{x1D70E}}B\unicode[STIX]{x1D716}/|v_{\Vert }|\text{d}\unicode[STIX]{x1D706}\,\text{d}\unicode[STIX]{x1D716}$. Upon integration by parts in 
$\unicode[STIX]{x1D706}$ and applying (B 11), the following expression is obtained,
$$\begin{eqnarray}\displaystyle \langle \unicode[STIX]{x1D71E}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}\rangle _{\unicode[STIX]{x1D713}} & = & \displaystyle -\frac{4\sqrt{2}\unicode[STIX]{x03C0}}{V^{\prime }(\unicode[STIX]{x1D713})}\left(\frac{mc}{3q}\right)^{2}\int _{0}^{\infty }\text{d}\unicode[STIX]{x1D716}\,\left(\frac{\unicode[STIX]{x2202}f_{0}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}}\right)\frac{\unicode[STIX]{x1D716}^{5/2}}{\unicode[STIX]{x1D708}(\unicode[STIX]{x1D716})}\nonumber\\ \displaystyle & & \displaystyle \times \int _{1/B_{\max }}^{1/B_{\min }}\frac{\text{d}\unicode[STIX]{x1D706}}{\unicode[STIX]{x1D706}}\,\int _{0}^{2\unicode[STIX]{x03C0}}\text{d}\unicode[STIX]{x1D6FC}\,\mathop{\sum }_{i}\frac{\displaystyle \left(\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D6FC}}\hat{K}_{i}(\unicode[STIX]{x1D6FC},\unicode[STIX]{x1D706})\right)^{2}}{\hat{I} _{i}(\unicode[STIX]{x1D6FC},\unicode[STIX]{x1D706})},\end{eqnarray}$$ where the bounce integrals are defined by (6.2) and (6.3). The sum in (B 13) is taken over trapping regions for particles with pitch angle 
$\unicode[STIX]{x1D706}$ on a field line labelled by 
$\unicode[STIX]{x1D6FC}$. The sum is carried out for left bounce points 
$\unicode[STIX]{x1D701}_{-,i}\in [0,2\unicode[STIX]{x03C0})$.
The parameter 
$\unicode[STIX]{x1D716}_{\text{eff}}^{3/2}$ quantifies the geometric dependence of the 
$1/\unicode[STIX]{x1D708}$ particle flux. It is defined in terms of the radial particle flux in the following way (Nemov et al. Reference Nemov, Kasilov, Kernbichler and Heyn1999),
$$\begin{eqnarray}\displaystyle \langle \unicode[STIX]{x1D71E}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}\rangle _{\unicode[STIX]{x1D713}}=-32\langle |\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}|\rangle _{\unicode[STIX]{x1D713}}^{2}\left(\frac{mc}{3q}\right)^{2}\frac{1}{B_{0}^{2}R^{2}}\unicode[STIX]{x1D716}_{\text{eff}}^{3/2}\int _{0}^{\infty }\text{d}\unicode[STIX]{x1D716}\,\left(\frac{\unicode[STIX]{x2202}f_{0}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}}\right)\frac{\unicode[STIX]{x1D716}^{5/2}}{\unicode[STIX]{x1D708}(\unicode[STIX]{x1D716})}. & & \displaystyle\end{eqnarray}$$ We take our normalizing length and field values to be such that 
$B_{0}R=\unicode[STIX]{x1D716}_{\text{ref}}^{-1}\langle |\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}|\rangle _{\unicode[STIX]{x1D713}}$, where 
$\unicode[STIX]{x1D716}_{\text{ref}}$ is a reference aspect ratio. Comparing (B 13) with (B 14) we obtain the expression for 
$\unicode[STIX]{x1D716}_{\text{eff}}^{3/2}$ (6.1). The corresponding expression (29) in Nemov et al. (Reference Nemov, Kasilov, Kernbichler and Heyn1999) is obtained by noting that 
${\hat{H}}^{\text{Nemov}}=-(\unicode[STIX]{x2202}\hat{K}/\unicode[STIX]{x2202}\unicode[STIX]{x1D6FC})\unicode[STIX]{x1D706}^{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 Nemov et al. (Reference Nemov, Kasilov, Kernbichler and Heyn1999).
The shape derivative of 
$f_{\unicode[STIX]{x1D716}}$ (6.4) is computed to be
$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6FF}f_{\unicode[STIX]{x1D716}}(S_{P};\unicode[STIX]{x1D743}_{1})=\int _{V_{P}}\text{d}\unicode[STIX]{x1D713}\,w(\unicode[STIX]{x1D713})\unicode[STIX]{x1D6FF}(V^{\prime }(\unicode[STIX]{x1D713})\unicode[STIX]{x1D716}_{\text{eff}}^{3/2}(\unicode[STIX]{x1D713})). & & \displaystyle\end{eqnarray}$$ 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}}\text{d}l\,Q$ due to displacement of the integration curve by vector field 
$\unicode[STIX]{x1D6FF}\boldsymbol{r}$ (Antonsen & Lee Reference Antonsen and Lee1982; Landreman & Paul Reference Landreman and Paul2018),
$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6FF}Q_{L}=\int _{l_{0}}^{l_{L}}\text{d}l\,(\unicode[STIX]{x1D6FF}\boldsymbol{r}\boldsymbol{\cdot }(-\unicode[STIX]{x1D73F}Q+(\unicode[STIX]{x1D644}-\boldsymbol{t}\boldsymbol{t})\boldsymbol{\cdot }\unicode[STIX]{x1D735}Q)+\unicode[STIX]{x1D6FF}Q)+Q(l_{L})\unicode[STIX]{x1D6FF}l_{L}-Q(l_{0})\unicode[STIX]{x1D6FF}l_{0}, & & \displaystyle\end{eqnarray}$$ where 
$\unicode[STIX]{x1D6FF}Q$ is the perturbation to the integrand at fixed position, 
$\boldsymbol{t}=\unicode[STIX]{x2202}\boldsymbol{r}/\unicode[STIX]{x2202}l$ is the unit tangent vector, 
$\unicode[STIX]{x1D73F}=\unicode[STIX]{x2202}^{2}\boldsymbol{r}/\unicode[STIX]{x2202}l^{2}$ is the curvature and 
$\unicode[STIX]{x1D6FF}l_{L}$ and 
$\unicode[STIX]{x1D6FF}l_{0}$ are perturbations to the bounds of the integral.
We compute the perturbation to the bounce integrals to be
$$\begin{eqnarray}\displaystyle \displaystyle \unicode[STIX]{x1D6FF}\hat{I} _{i} & = & \displaystyle \oint \text{d}l\,\left(-\frac{v_{\Vert }}{vB}\unicode[STIX]{x1D73F}\boldsymbol{\cdot }\unicode[STIX]{x1D6FF}\boldsymbol{r}-\left(\frac{\unicode[STIX]{x1D706}v}{2Bv_{\Vert }}+\frac{v_{\Vert }}{B^{2}v}\right)(\unicode[STIX]{x1D6FF}\boldsymbol{r}\boldsymbol{\cdot }\unicode[STIX]{x1D735}B+\unicode[STIX]{x1D6FF}B)\right)\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle \displaystyle \unicode[STIX]{x1D6FF}\hat{K}_{i} & = & \displaystyle \oint \text{d}l\,\left(-\frac{\mathop{v}_{\Vert }^{3}}{v^{3}B}\unicode[STIX]{x1D73F}\boldsymbol{\cdot }\unicode[STIX]{x1D6FF}\boldsymbol{r}-\left(\frac{3\unicode[STIX]{x1D706}v_{\Vert }}{2Bv}+\frac{\mathop{v}_{\Vert }^{3}}{B^{2}v^{3}}\right)(\unicode[STIX]{x1D6FF}\boldsymbol{r}\boldsymbol{\cdot }\unicode[STIX]{x1D735}B+\unicode[STIX]{x1D6FF}B)\right),\end{eqnarray}$$ where 
$\unicode[STIX]{x1D6FF}B$ is the perturbation to the field strength (5.7) and 
$\unicode[STIX]{x1D6FF}\boldsymbol{r}$ is given by (3.10). We note that 
$\unicode[STIX]{x1D6FF}\boldsymbol{r}\boldsymbol{\cdot }\boldsymbol{b}=0$ such that the perpendicular projection, 
$(\unicode[STIX]{x1D644}-\boldsymbol{t}\boldsymbol{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 (6.5)–(6.8) can now be obtained by writing (B 15) in terms of the perturbations of the bounce integrals, using 
$\unicode[STIX]{x1D743}_{1}\boldsymbol{\cdot }\unicode[STIX]{x1D735}B+\unicode[STIX]{x1D6FF}B=-B(\unicode[STIX]{x1D644}-\boldsymbol{b}\boldsymbol{b}):\unicode[STIX]{x1D735}\unicode[STIX]{x1D743}_{1}-\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D704}_{1}\boldsymbol{b}\boldsymbol{\cdot }(\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}\times \unicode[STIX]{x1D735}\unicode[STIX]{x1D701})$ and 
$\unicode[STIX]{x1D73F}\boldsymbol{\cdot }\unicode[STIX]{x1D743}_{1}=-\boldsymbol{b}\boldsymbol{b}:\unicode[STIX]{x1D735}\unicode[STIX]{x1D743}_{1}$.
Appendix C. Details of departure from quasisymmetry calculation
In this appendix we compute the shape derivative of 
$f_{QS}$ (7.3) to obtain (7.8)–(7.11) by expressing each term in (7.7) in the desired form. The second term in (7.7) is expressed using 
$\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D713}=-\unicode[STIX]{x1D743}_{1}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}$,
$$\begin{eqnarray}\displaystyle \frac{1}{2}\int _{V_{P}}\text{d}^{3}x\,w^{\prime }(\unicode[STIX]{x1D713})\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D713}{\mathcal{M}}^{2}=-\frac{1}{2}\int _{V_{P}}\text{d}^{3}x\,{\mathcal{M}}^{2}\unicode[STIX]{x1D743}_{1}\boldsymbol{\cdot }\unicode[STIX]{x1D735}w(\unicode[STIX]{x1D713}). & & \displaystyle\end{eqnarray}$$ The third term in (7.7) is computed upon application of (3.5), the divergence theorem and noting that 
${\mathcal{M}}=\boldsymbol{B}\boldsymbol{\cdot }\boldsymbol{{\mathcal{A}}}$,
$$\begin{eqnarray}\displaystyle & & \displaystyle \int _{V_{P}}\text{d}^{3}x\,w(\unicode[STIX]{x1D713}){\mathcal{M}}\unicode[STIX]{x1D6FF}\boldsymbol{B}\boldsymbol{\cdot }\boldsymbol{{\mathcal{A}}}\nonumber\\ \displaystyle & & \displaystyle \quad =-\int _{S_{P}}\text{d}^{2}x\,\unicode[STIX]{x1D743}_{1}\boldsymbol{\cdot }\boldsymbol{n}w(\unicode[STIX]{x1D713}){\mathcal{M}}^{2}-\int _{V_{P}}\text{d}^{3}x\,w(\unicode[STIX]{x1D713})\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D704}_{1}(\unicode[STIX]{x1D713}){\mathcal{M}}\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}\times \unicode[STIX]{x1D735}\unicode[STIX]{x1D701}\boldsymbol{\cdot }\boldsymbol{{\mathcal{A}}}\nonumber\\ \displaystyle & & \displaystyle \qquad +\int _{V_{P}}\text{d}^{3}x\,\unicode[STIX]{x1D743}_{1}\boldsymbol{\cdot }(w(\unicode[STIX]{x1D713}){\mathcal{M}}(\boldsymbol{B}\times (\unicode[STIX]{x1D735}\times \boldsymbol{{\mathcal{A}}}))-\boldsymbol{{\mathcal{A}}}w(\unicode[STIX]{x1D713})\boldsymbol{B}\boldsymbol{\cdot }\unicode[STIX]{x1D735}{\mathcal{M}}+{\mathcal{M}}\unicode[STIX]{x1D735}(w(\unicode[STIX]{x1D713}){\mathcal{M}})).\nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$ The quantity 
$\boldsymbol{{\mathcal{A}}}$ can be projected into the perpendicular direction as 
$\unicode[STIX]{x1D743}_{1}\boldsymbol{\cdot }\boldsymbol{b}=0$, noting that
$$\begin{eqnarray}\displaystyle \boldsymbol{b}\times (\boldsymbol{{\mathcal{A}}}\times \boldsymbol{b})=-(\boldsymbol{b}\times \unicode[STIX]{x1D735}\unicode[STIX]{x1D713})\unicode[STIX]{x1D735}_{\Vert }B-F(\unicode[STIX]{x1D713})\unicode[STIX]{x1D735}_{\bot }B. & & \displaystyle\end{eqnarray}$$ Similarly, any terms in (C 2) involving 
$\unicode[STIX]{x1D743}_{1}\boldsymbol{\cdot }\unicode[STIX]{x1D735}$ can be expressed as 
$\unicode[STIX]{x1D743}_{1}\boldsymbol{\cdot }\unicode[STIX]{x1D735}_{\bot }$. The corresponding terms in (7.9) are obtained using the expression for the curvature in an equilibrium field. The fourth term in (7.7) is expressed in the following way upon application of (5.7), the divergence theorem, and noting that 
$\boldsymbol{{\mathcal{S}}}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}=\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{{\mathcal{S}}}=0$,
$$\begin{eqnarray}\displaystyle & & \displaystyle \int _{V_{P}}\text{d}^{3}x\,w(\unicode[STIX]{x1D713}){\mathcal{M}}\boldsymbol{{\mathcal{S}}}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D6FF}B\nonumber\\ \displaystyle & & \displaystyle \quad =\int _{S_{P}}\text{d}^{2}x\,\unicode[STIX]{x1D743}_{1}\boldsymbol{\cdot }\boldsymbol{n}Bw(\unicode[STIX]{x1D713})\boldsymbol{{\mathcal{S}}}\boldsymbol{\cdot }\unicode[STIX]{x1D735}{\mathcal{M}}-\int _{V_{P}}\text{d}^{3}x\,\unicode[STIX]{x1D743}_{1}\boldsymbol{\cdot }[B\unicode[STIX]{x1D735}(w(\unicode[STIX]{x1D713})\boldsymbol{{\mathcal{S}}}\boldsymbol{\cdot }\unicode[STIX]{x1D735}{\mathcal{M}})]\nonumber\\ \displaystyle & & \displaystyle \qquad +\int _{V_{P}}\text{d}^{3}x\,w(\unicode[STIX]{x1D713})(\boldsymbol{{\mathcal{S}}}\boldsymbol{\cdot }\unicode[STIX]{x1D735}{\mathcal{M}})(\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D704}_{1}(\unicode[STIX]{x1D713})\boldsymbol{b}\boldsymbol{\cdot }(\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}\times \unicode[STIX]{x1D735}\unicode[STIX]{x1D701})+B\unicode[STIX]{x1D743}_{1}\boldsymbol{\cdot }\unicode[STIX]{x1D73F}).\end{eqnarray}$$ We express terms involving 
$\unicode[STIX]{x1D743}_{1}\boldsymbol{\cdot }\unicode[STIX]{x1D735}$ as 
$\unicode[STIX]{x1D743}_{1}\boldsymbol{\cdot }\unicode[STIX]{x1D735}_{\bot }$ to obtain the corresponding terms in (7.9). The fifth term in (7.7) is expressed in the following way upon application of 
$\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D713}=-\unicode[STIX]{x1D743}_{1}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}$, the divergence theorem and several vector identities,
$$\begin{eqnarray}\displaystyle \int _{V_{P}}\text{d}^{3}x\,w(\unicode[STIX]{x1D713}){\mathcal{M}}\boldsymbol{B}\times \unicode[STIX]{x1D735}\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D713}\boldsymbol{\cdot }\unicode[STIX]{x1D735}B & = & \displaystyle -\int _{S_{P}}\text{d}^{2}x\,\unicode[STIX]{x1D743}_{1}\boldsymbol{\cdot }\boldsymbol{n}w(\unicode[STIX]{x1D713}){\mathcal{M}}\unicode[STIX]{x1D735}B\times \boldsymbol{B}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}\nonumber\\ \displaystyle & & \displaystyle -\,\int _{V_{P}}\text{d}^{3}x\,\unicode[STIX]{x1D743}_{1}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}\unicode[STIX]{x1D735}B\boldsymbol{\cdot }\unicode[STIX]{x1D735}\times (w(\unicode[STIX]{x1D713}){\mathcal{M}}\boldsymbol{B}).\qquad\end{eqnarray}$$The sixth term in (7.7) upon application of (7.6) is,
$$\begin{eqnarray}\displaystyle & & \displaystyle -\int _{V_{P}}\text{d}^{3}x\,\frac{\unicode[STIX]{x1D6FF}G(\unicode[STIX]{x1D713})w(\unicode[STIX]{x1D713}){\mathcal{M}}\boldsymbol{B}\boldsymbol{\cdot }\unicode[STIX]{x1D735}B}{\unicode[STIX]{x1D704}(\unicode[STIX]{x1D713})-(N/M)}\nonumber\\ \displaystyle & & \displaystyle \quad =\frac{1}{4\unicode[STIX]{x03C0}^{2}}\int _{S_{P}}\text{d}^{2}x\,\frac{w(\unicode[STIX]{x1D713})V^{\prime }(\unicode[STIX]{x1D713})\langle {\mathcal{M}}\boldsymbol{B}\boldsymbol{\cdot }\unicode[STIX]{x1D735}B\rangle _{\unicode[STIX]{x1D713}}}{(\unicode[STIX]{x1D704}(\unicode[STIX]{x1D713})-(N/M))}(\boldsymbol{B}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}\times \unicode[STIX]{x1D735}\unicode[STIX]{x1D703})\unicode[STIX]{x1D743}_{1}\boldsymbol{\cdot }\boldsymbol{n}\nonumber\\ \displaystyle & & \displaystyle \qquad -\,\frac{1}{4\unicode[STIX]{x03C0}^{2}}\int _{V_{P}}\text{d}^{3}x\,\unicode[STIX]{x1D743}_{1}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\left(\frac{w(\unicode[STIX]{x1D713})V^{\prime }(\unicode[STIX]{x1D713})\langle {\mathcal{M}}\boldsymbol{B}\boldsymbol{\cdot }\unicode[STIX]{x1D735}B\rangle _{\unicode[STIX]{x1D713}}}{(\unicode[STIX]{x1D704}(\unicode[STIX]{x1D713})-(N/M))}\right)\boldsymbol{B}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}\times \unicode[STIX]{x1D735}\unicode[STIX]{x1D703}\nonumber\\ \displaystyle & & \displaystyle \qquad +\,\frac{1}{4\unicode[STIX]{x03C0}^{2}}\int _{V_{P}}\text{d}^{3}x\,\frac{w(\unicode[STIX]{x1D713})V^{\prime }(\unicode[STIX]{x1D713})\langle {\mathcal{M}}\boldsymbol{B}\boldsymbol{\cdot }\unicode[STIX]{x1D735}B\rangle _{\unicode[STIX]{x1D713}}}{\unicode[STIX]{x1D704}(\unicode[STIX]{x1D713})-(N/M)}\nonumber\\ \displaystyle & & \displaystyle \quad \times \,(\unicode[STIX]{x1D743}_{1}\boldsymbol{\cdot }(\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\boldsymbol{B}\times \unicode[STIX]{x1D735}\unicode[STIX]{x1D703})-\boldsymbol{B}\times \unicode[STIX]{x1D735}\times (\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}\times \unicode[STIX]{x1D735}\unicode[STIX]{x1D703})))\nonumber\\ \displaystyle & & \displaystyle \qquad -\,\frac{1}{4\unicode[STIX]{x03C0}^{2}}\int _{V_{P}}\text{d}^{3}x\,\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D704}_{1}(\unicode[STIX]{x1D713})\frac{w(\unicode[STIX]{x1D713})V^{\prime }(\unicode[STIX]{x1D713})\langle {\mathcal{M}}\boldsymbol{B}\boldsymbol{\cdot }\unicode[STIX]{x1D735}B\rangle _{\unicode[STIX]{x1D713}}}{\sqrt{g}^{2}(\unicode[STIX]{x1D704}(\unicode[STIX]{x1D713})-(N/M))}\frac{\unicode[STIX]{x2202}\boldsymbol{r}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D701}}\boldsymbol{\cdot }\frac{\unicode[STIX]{x2202}\boldsymbol{r}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D703}}.\end{eqnarray}$$ In obtaining the corresponding terms in (7.9), terms involving 
$\unicode[STIX]{x1D743}_{1}\boldsymbol{\cdot }\unicode[STIX]{x1D735}$ are expressed as 
$\unicode[STIX]{x1D743}_{1}\boldsymbol{\cdot }\unicode[STIX]{x1D735}_{\bot }$. The seventh term in (7.7) is expressed using 
$\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D713}=-\unicode[STIX]{x1D743}_{1}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}$. Combining all terms, we obtain (7.8)–(7.11).
Appendix D. Details of neoclassical figures of merit calculation
In this section we compute the shape derivative of 
$f_{NC}$ (8.1) to obtain (8.7)–(8.10) by expressing each term in (8.6) in the desired form. Throughout, Boozer coordinates will be assumed.
The second term in (8.6) is expressed using 
$\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D713}=-\unicode[STIX]{x1D743}_{1}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}$. The third term in (8.6) can be computed using (7.6), noting that 
$V^{\prime }(\unicode[STIX]{x1D713})/(4\unicode[STIX]{x03C0}^{2}\sqrt{g})=B^{2}/\langle B^{2}\rangle$ in Boozer coordinates and applying the divergence theorem,
$$\begin{eqnarray}\displaystyle & & \displaystyle \int _{V_{P}}\text{d}^{3}x\,w(\unicode[STIX]{x1D713})\frac{\unicode[STIX]{x2202}{\mathcal{R}}(\unicode[STIX]{x1D713})}{\unicode[STIX]{x2202}G(\unicode[STIX]{x1D713})}\unicode[STIX]{x1D6FF}G(\unicode[STIX]{x1D713})=-\int _{V_{P}}\text{d}^{3}x\,w(\unicode[STIX]{x1D713})\frac{B^{2}\sqrt{g}}{\langle B^{2}\rangle _{\unicode[STIX]{x1D713}}}\frac{\unicode[STIX]{x2202}{\mathcal{R}}(\unicode[STIX]{x1D713})}{\unicode[STIX]{x2202}G(\unicode[STIX]{x1D713})}\unicode[STIX]{x1D743}_{1}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}(\unicode[STIX]{x1D735}\times \boldsymbol{B})\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D703}\nonumber\\ \displaystyle & & \displaystyle \quad +\,\int _{V_{P}}\text{d}^{3}x\,\left(\unicode[STIX]{x1D743}_{1}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\left(\frac{\unicode[STIX]{x2202}{\mathcal{R}}(\unicode[STIX]{x1D713})}{\unicode[STIX]{x2202}G(\unicode[STIX]{x1D713})}\frac{w(\unicode[STIX]{x1D713})}{\langle B^{2}\rangle _{\unicode[STIX]{x1D713}}}\right)B^{2}G(\unicode[STIX]{x1D713})\right.\nonumber\\ \displaystyle & & \displaystyle \quad +\left.\frac{w(\unicode[STIX]{x1D713})}{\langle B^{2}\rangle _{\unicode[STIX]{x1D713}}}\frac{\unicode[STIX]{x2202}{\mathcal{R}}(\unicode[STIX]{x1D713})}{\unicode[STIX]{x2202}G(\unicode[STIX]{x1D713})}\unicode[STIX]{x1D743}_{1}\boldsymbol{\cdot }\boldsymbol{B}\times \unicode[STIX]{x1D735}\times \left(\frac{\unicode[STIX]{x2202}\boldsymbol{r}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D701}}B^{2}\right)\right)\nonumber\\ \displaystyle & & \displaystyle \quad +\,\int _{V_{P}}\text{d}^{3}x\,\frac{w(\unicode[STIX]{x1D713})\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D704}_{1}(\unicode[STIX]{x1D713})B^{2}}{\sqrt{g}\langle B^{2}\rangle _{\unicode[STIX]{x1D713}}}\frac{\unicode[STIX]{x2202}{\mathcal{R}}(\unicode[STIX]{x1D713})}{\unicode[STIX]{x2202}G(\unicode[STIX]{x1D713})}\frac{\unicode[STIX]{x2202}\boldsymbol{r}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D701}}\boldsymbol{\cdot }\frac{\unicode[STIX]{x2202}\boldsymbol{r}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D703}}\nonumber\\ \displaystyle & & \displaystyle \quad -\int _{S_{P}}\text{d}^{2}x\,w(\unicode[STIX]{x1D713})\frac{B^{2}}{\langle B^{2}\rangle _{\unicode[STIX]{x1D713}}}\frac{\unicode[STIX]{x2202}{\mathcal{R}}(\unicode[STIX]{x1D713})}{\unicode[STIX]{x2202}G(\unicode[STIX]{x1D713})}G(\unicode[STIX]{x1D713})\unicode[STIX]{x1D743}_{1}\boldsymbol{\cdot }\boldsymbol{n}.\end{eqnarray}$$The fifth term in (8.6) can be computed using (5.7), the divergence theorem and the expression for the curvature in an equilibrium field,
$$\begin{eqnarray}\displaystyle & & \displaystyle \int _{V_{P}}\text{d}^{3}x\,w(\unicode[STIX]{x1D713})\langle S_{{\mathcal{R}}}\unicode[STIX]{x1D6FF}B\rangle _{\unicode[STIX]{x1D713}}=\int _{V_{P}}\text{d}^{3}x\,(\unicode[STIX]{x1D743}_{1}\boldsymbol{\cdot }\unicode[STIX]{x1D735}(w(\unicode[STIX]{x1D713})S_{{\mathcal{R}}})B-BS_{{\mathcal{R}}}w(\unicode[STIX]{x1D713})\unicode[STIX]{x1D743}_{1}\boldsymbol{\cdot }\unicode[STIX]{x1D73F})\nonumber\\ \displaystyle & & \displaystyle \quad -\,\int _{V_{P}}\text{d}^{3}x\,\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D704}_{1}(\unicode[STIX]{x1D713})S_{{\mathcal{R}}}w(\unicode[STIX]{x1D713})\boldsymbol{b}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}\times \unicode[STIX]{x1D735}\unicode[STIX]{x1D701}-\int _{S_{P}}\text{d}^{2}x\,w(\unicode[STIX]{x1D713})S_{{\mathcal{R}}}B\unicode[STIX]{x1D743}_{1}\boldsymbol{\cdot }\boldsymbol{n}.\end{eqnarray}$$ The resulting terms can be combined to write the shape derivative in the form of (8.7), noting that any terms involving 
$\unicode[STIX]{x1D743}_{1}\boldsymbol{\cdot }\unicode[STIX]{x1D735}$ can be expressed as 
$\unicode[STIX]{x1D743}_{1}\boldsymbol{\cdot }\unicode[STIX]{x1D735}_{\bot }$.
