Confinement physics: quasisymmetry, omnigenity, stability

The pages Equations and derivations and Algorithms derive the solver physics — the energy functional whose stationary point is an equilibrium, and the numerics that find it. This page derives the target functionals that a stellarator-design campaign puts in front of that equilibrium: the confinement and stability metrics of vmex.core.omnigenity, vmex.core.optimize, vmex.core.nyquist, vmex.core.stability and vmex.core.bootstrap. It is the first-principles companion to Objectives library, which catalogs the same quantities as ready-to-use optimizer terms.

The unifying object is the field strength \(|B|\) expressed in Boozer coordinates: quasisymmetry, omnigenity, and the neoclassical transport that they control are all statements about the angular structure of \(|B|\) at fixed flux label, and are cleanest in the coordinates in which field lines are straight and the parallel current is a flux function.

Boozer coordinates

Motivation

Guiding-center orbits — and therefore neoclassical transport, fast-ion confinement, and the bootstrap current — depend on the magnetic geometry almost entirely through the variation of \(|B|\) along and across field lines. Boozer coordinates \((s, \theta_B, \zeta_B)\) are the flux coordinates that make this variation maximally transparent: field lines are straight,

\[\frac{d\theta_B}{d\zeta_B} = \iota(s),\]

and the covariant representation of \(\mathbf B\) collapses to the two surface functions

\[\mathbf B = \nabla\psi\times\nabla\theta_B + \iota\,\nabla\zeta_B\times\nabla\psi = I(s)\,\nabla\theta_B + G(s)\,\nabla\zeta_B + \beta_*\,\nabla\psi,\]

where \(2\pi I(s)\) is the toroidal current inside \(s\), \(2\pi G(s)\) the poloidal current outside it (VMEC’s buco and bvco, surface_currents()), and \(\psi=\Phi/2\pi\) is the toroidal flux per radian. Because \(I\) and \(G\) are flux functions, the magnetic differential equation for the coordinate shift is linear and can be solved surface by surface.

Construction

Starting from the straight-field-line (VMEC) angles \((\theta,\zeta)\) with the renormalized poloidal angle \(u=\theta+\lambda\), the Boozer angles differ by a single periodic scalar \(\nu(s,\theta,\zeta)\):

\[\theta_B = \theta + \lambda + \iota\,\nu, \qquad \zeta_B = \zeta + \nu .\]

The shift \(\nu\) follows from the covariant field. Introduce the periodic part \(w\) of the Boozer generating potential, fixed by

\[\frac{\partial w}{\partial\theta} = B_\theta, \qquad \frac{\partial w}{\partial\zeta} = B_\zeta,\]

with \(B_\theta,B_\zeta\) the covariant components on the surface (magnetic_fields()). In boozer_bmnc_state() this is inverted spectrally: after an FFT, the non-axisymmetric (\(m\ne 0\)) harmonics of \(w\) come from \(B_\theta\) and the axisymmetric (\(m=0\)) harmonics from \(B_\zeta\), matching the mode split of the Fortran booz_xform. The coordinate shift is then

\[\nu = \frac{w - I(s)\,\lambda}{G(s) + \iota(s)\,I(s)} .\]

Finally the Boozer \(|B|\) spectrum is obtained by the angle-transform quadrature, weighting \(|B|\) (evaluated on the VMEC grid) by the Jacobian of the angle map:

\[\hat B^{B}_{mn}(s) = \Bigl\langle\,|B|\, \cos\!\bigl(m\,\theta_B - n\,\zeta_B\bigr)\, \frac{\partial(\theta_B,\zeta_B)}{\partial(\theta,\zeta)}\Bigr\rangle,\]

the angle bracket being the surface average over the original \((\theta,\zeta)\) grid. These are the bmnc_b coefficients (physical mode numbers xm_b, xn_b) consumed by every metric below.

Two implementations share these equations. The host driver vmex.core.boozer.run_booz_xform() calls booz_xform_jax on a wout_*.nc file and writes a standard boozmn_*.nc (used by vmec --booz and by quasi_isodynamic_residual_from_wout()). The traceable boozer_bmnc_state() evaluates the same transform in pure jax.numpy directly from the solver’s internal half-mesh field tables, so the Boozer spectrum — and any metric built on it — carries exact implicit gradients. The two agree to \(\sim10^{-6}\) on the dominant modes.

Quasisymmetry

Confinement rationale

A field is quasisymmetric when \(|B|\) depends on the two Boozer angles only through a single linear combination,

\[|B| = |B|\bigl(s,\; M\theta_B - N\zeta_B\bigr),\]

for a fixed integer helicity \((M,N)\). The magnitude, not the vector field, need only be symmetric: this is enough for the guiding-center Lagrangian to acquire an ignorable angle, hence a conserved canonical momentum and confined collisionless orbits — the same confinement a true axisymmetric tokamak enjoys, but in a compact 3D device. The helicity fixes the family (QuasisymmetryRatioResidual, helicity_n in units of nfp):

family	\((M,N)\)	\(|B|\) contours in Boozer angles
QA (quasi-axisymmetric)	\((1,0)\)	close poloidally (tokamak-like)
QH (quasi-helical)	\((1,\pm\mathrm{nfp})\)	close helically
QP (quasi-poloidal)	\((0,1)\)	close toroidally

The two-term residual

Rather than Fourier-filter the Boozer spectrum, vmex uses the Landreman–Paul two-term local residual, which is an exact pointwise diagnostic of the same condition and needs no mode truncation. On each requested surface, sampled on a uniform \((\theta,\phi)\) grid,

\[f_{\mathrm{QS}} = \frac{(\mathbf B\times\nabla B\cdot\nabla\psi)\,(N-\iota M) - (\mathbf B\cdot\nabla B)\,(M G + N I)}{B^{3}},\]

with \(M=\) helicity_m, \(N=\) helicity_n\(\times\)nfp, and \(G,I\) the Boozer covariant averages bvco/buco. The key algebraic fact is that \(f_{\mathrm{QS}}\) vanishes identically iff \(|B|\) is quasisymmetric with helicity \((M,N)\); there is no residual symmetry-breaking harmonic left to penalize. The flux-surface sum \(\sum f_{\mathrm{QS}}^2\), weighted by the surface measure \(\sqrt{\mathrm{nfp}\,\Delta\theta\,\Delta\phi\,|\sqrt g|/V'}\), reproduces simsopt’s QuasisymmetryRatioResidual A/B bit-for-bit. Kept in Gauss–Newton (per-point) form, it feeds the least-squares driver as an exact residual vector rather than a pre-summed scalar. The metric is evaluated from the parity-proven wout tables (vmex.core.nyquist) and also exposes a traceable residuals_state lane, so the same term optimizes under both jac=None and jac="implicit" (see Objectives library).

Omnigenity and quasi-isodynamicity

Omnigenity generalizes quasisymmetry: it asks only that the bounce-averaged radial drift of trapped particles vanish, without requiring a symmetry of \(|B|\). Equivalently (Cary & Shasharina 1997), the second adiabatic invariant

\[\mathcal J_\parallel(s,\alpha;B^*) = \oint v_\parallel\,d\ell, \qquad v_\parallel \propto \sqrt{1 - B/B^*},\]

is a flux function — independent of the field-line label \(\alpha\) at every trapping level \(B^*\). Every quasisymmetric field is omnigenous; the converse is not true, which leaves omnigenity the larger (and for poloidally-closed contours, bootstrap-suppressing) design space. A quasi- isodynamic (QI) field is an omnigenous field whose \(|B|\) contours close poloidally (\(M=0\)), so that the trapped-particle precession is purely poloidal and the bootstrap current is small by construction — the target of the nfp1–nfp4 decks in examples/data/.

The constructed-QI target

Because \(\mathcal J_\parallel\) uniformity is awkward to differentiate, QIResidual implements the equivalent level-set conditions of the constructed-QI-target method of Goodman et al. (2023). On each surface, \(|B|\) is sampled along Boozer field lines \(\theta_B=\alpha+\iota\phi_B\) over one field period, and three families of residual — each an exact zero of an exactly QI field — are stacked:

• Bounce-distance uniformity (well_weight). For every trapping level \(B^*\), the distance \(\delta(\alpha,B^*)\) between the two monotone branches of the magnetic well, minus its field-line average. This is the Cary–Shasharina condition (constant well width at fixed \(B^*\) across field lines) that Goodman’s “shuffle” step enforces; the branch envelopes are built from smooth running-maximum occupancy integrals so the term is differentiable.

• Extremum alignment (extremum_weight). The per-field-line \(B_{\min}\) and \(B_{\max}\) minus their field-line averages — poloidal closure of the extremal \(|B|\) contours (Goodman’s “align the maxima” step; also the flat-\(B_{\max}\) condition of Dudt et al. 2024).

• Single-well monotonicity (squash_weight). The pointwise distance between \(|B|\) and its monotone branch envelopes — Goodman’s “squash” distance, which penalizes side wells (more than one magnetic well per period).

Every operation (sigmoid occupancies, running maxima, level-space quadrature) is smooth or piecewise-smooth, so the residual is jit/grad/jvp-transparent and QI optimization runs with the exact implicit adjoint, exactly like the QS residual.

Metric fidelity: report the wout-based residual

The traceable QIResidual is built for gradient flow, not for labeling a configuration. Its level-space envelopes and finite sampling make it an excellent descent direction but an optimistic absolute scale: driven hard, it can report values well below what an independent Boozer analysis confirms. When quoting “how QI is this equilibrium?”, use the wout/Boozer-based quasi_isodynamic_residual_from_wout() (host booz_xform_jax, finite-difference-only) — the same construction evaluated on the fully resolved Boozer spectrum of the written wout. The recommended workflow is therefore: optimize with the traceable residual under jac="implicit", then report the wout-lane residual. The two agree in ranking; they diverge in absolute value precisely where the traceable form is being pushed below its trustworthy range.

Ideal MHD stability: Mercier and the magnetic well

Mercier criterion

The Mercier criterion is the ideal-MHD stability condition for localized interchange modes resonant on a rational surface — the toroidal generalization of the Suydam criterion. Following Glasser–Greene–Johnson, the surface is Mercier-stable where

\[D_{\mathrm{Merc}} = D_{\mathrm{shear}} + D_{\mathrm{well}} + D_{\mathrm{curr}} + D_{\mathrm{geod}} > 0 ,\]

a decomposition into four flux-surface integrals with distinct physical origin. mercier_and_jxb() ports VMEC2000 mercier.f term for term (with \({}' \equiv d/d\psi\), \(p\) the pressure, \(\iota\) the transform, \(V'\) the differential volume, and \(I_\varphi\) the enclosed toroidal current):

term	sign	origin
\(D_{\mathrm{shear}} = \tfrac14(\iota')^2\)	\(\ge 0\)	magnetic shear (always stabilizing)
\(D_{\mathrm{well}}\)	either	pressure gradient \(\times\) magnetic well \(V''\)
\(D_{\mathrm{curr}}\)	either	parallel-current / kink drive \(\propto \iota'\)
\(D_{\mathrm{geod}}\)	\(\le 0\)	geodesic curvature (always destabilizing)

The well term carries the sign of \(p'\,(V'' - p'\langle\cdots\rangle)\): with the pressure decreasing outward (\(p'<0\)), a magnetic well (\(V''<0\), volume that decreases toward the edge) is stabilizing. The geodesic term is a manifestly non-positive Schwarz-inequality remainder. Because the individual pieces involve radial derivatives of surface averages, the two surfaces nearest the axis and the edge carry the usual numerical noise; a practical objective penalizes min(DMerc[2:-1], 0). vmex exposes the profile as d_merc(), evaluated through the parity-proven wout engine (host NumPy, hence jac=None); it is validated against VMEC2000 golden wout files.

Magnetic well

The dominant stabilizing ingredient of \(D_{\mathrm{well}}\) — the sign of \(V''(s)\) — is also useful on its own as a cheap, fully traceable proxy. magnetic_well() returns the finite-difference vacuum-well measure

\[W = \frac{V'(0) - V'(1)}{V'(0)},\]

with \(V'=dV/ds\) extrapolated from the half-mesh differential volume \(vp\) (VMEC bcovar.f). Positive \(W\) means \(V'\) decreases outward — a magnetic well, favorable for interchange stability — matching simsopt’s vacuum_well. Being a pure (state, runtime) function it carries exact implicit gradients and is the traceable stand-in whenever the full DMerc profile is too expensive or not differentiable. Near-axis analytic context for both measures is in Landreman–Jorge (2020) and Kim–Jorge–Dorland (2021); see References.

Ideal ballooning

Interchange stability bounds only the \(n\to\infty\), radially-localized limit. The complementary ballooning limit — high toroidal mode number, extended along the field line — is provided by vmex.core.stability as a fully differentiable eigenvalue objective (a JAX port of the COBRA solve in the Gaur et al. formulation). It solves the self-adjoint field-line ODE eigenproblem

\[\frac{d}{d\eta}\!\Bigl(g\,\frac{dX}{d\eta}\Bigr) + c\,X = \lambda\,f\,X, \qquad X(\pm\eta_b)=0,\]

along the straight-field-line angle \(\eta\), where \(g\) is the line-bending term, \(c\) the pressure/curvature drive, and \(f>0\) the inertia; \(\lambda = (\gamma a_N/v_A)^2 > 0\) flags instability. ballooning_growth_rate() reduces the batched eigenvalues to a smooth softmax scalar built to be driven negative as a stable-by-construction constraint. The full coefficient definitions are in the stability module docstring and the usage recipe in Objectives library.

Bootstrap current (Redl)

The self-consistent parallel current a stellarator generates from its own pressure gradient — the bootstrap current — sets the achievable \(\beta\) and, in a QI device, must be kept small. vmex.core.bootstrap provides a differentiable evaluation and a self-consistency loop (reproducing Landreman–Buller–Drevlak, arXiv:2205.02914).

Two independent estimates of \(\langle\mathbf J\cdot\mathbf B\rangle\) are compared. The equilibrium value follows from the exact MHD identity

\[\langle\mathbf J\cdot\mathbf B\rangle(s) = \frac{\langle B^2\rangle\,I'(s) + \mu_0\,I(s)\,p'(s)}{2\pi\,\psi_a}, \qquad I(s) = \mathrm{signgs}\,\frac{2\pi}{\mu_0}\,\mathrm{buco}(s),\]

(vmec_j_dot_B()). The kinetic value is the Redl et al. (2021) analytic closure (j_dot_B_redl()), a fit in the effective trapped fraction \(f_t\) and the Sauter collisionalities, with the quasisymmetry isomorphism \(\iota\to\iota-\mathrm{nfp}\,\mathrm{helicity\_n}\) applied as in simsopt. The trapped fraction itself uses the singularity-removing substitution \(y=\sqrt{1-\lambda B_{\max}}\) in

\[f_t = 1 - \tfrac34\langle B^2\rangle \int_0^{1/B_{\max}}\! \frac{\lambda\,d\lambda}{\langle\sqrt{1-\lambda B}\rangle},\]

evaluated with fixed-order Gauss–Legendre quadrature so the whole chain stays differentiable. Their normalized mismatch is the residual RedlBootstrapMismatch (the exact formula and the finite-beta profile conventions are in Equations and derivations); driving it to zero, optionally with current_dofs freed, yields a current profile consistent with the plasma the equilibrium describes.