  Stellarators are a promising route to steady-state fusion power. However, to
achieve the required confinement, the magnetic geometry must be highly
optimized. This optimization requires navigating high-dimensional spaces, often
necessitating the use of gradient-based methods. The gradient of the
neoclassical fluxes is expensive to compute with classical methods, requiring
$O(N)$ flux computations, where $N$ is the number of parameters. To reduce the
cost of the gradient computation, we present an adjoint method for computing
the derivatives of moments of the neoclassical distribution function for
stellarator optimization. The linear adjoint method allows derivatives of
quantities which depend on solutions of a linear system, such as moments of the
distribution function, to be computed with respect to many parameters from the
solution of only two linear systems. This reduces the cost of computing the
gradient to the point that the finite-collisionality neoclassical fluxes can be
used within an optimization loop.
  With the neoclassical adjoint method, we compute solutions of the drift
kinetic equation and an adjoint drift kinetic equation to obtain derivatives of
neoclassical quantities with respect to geometric parameters. When the number
of parameters in the derivative is large ($\mathcal{O}(10^2)$), this adjoint
method provides up to a factor of 200 reduction in cost. We demonstrate
adjoint-based optimization of the field strength to obtain minimal bootstrap
current on a surface. With adjoint-based derivatives, we also compute the local
sensitivity to magnetic perturbations on a flux surface and identify regions
where tight tolerances on error fields are required for control of the
bootstrap current or radial transport. Furthermore, the solve for the ambipolar
electric field is accelerated using a Newton method with derivatives obtained
from the adjoint method.
