Solving the 3D MHD equilibrium equations in toroidal geometry by Newton's method

We describe a novel form of Newton's method for computing 3D MHD equilibria. The method has been implemented as an extension to the hybrid spectral/finite-difference Princeton Iterative Equilibrium Solver (PIES) which normally uses Picard iteration on the full nonlinear MHD equilibrium equations. Computing the Newton functional derivative numerically is not feasible in a code of this type but we are able to do the calculation analytically in magnetic coordinates by considering the response of the plasma's Pfirsch-Schlüter currents to small changes in the magnetic field. Results demonstrate a significant advantage over Picard iteration in many cases, including simple finite-β stellarator equilibria. The method shows promise in cases that are difficult for Picard iteration, although it is sensitive to resolution and imperfections in the magnetic coordinates, and further work is required to adapt it to the presence of magnetic islands and stochastic regions.