Physics Model

NTX solves the monoenergetic neoclassical transport problem on a single flux surface using the Legendre-space formulation developed in Javier Escoto’s PhD thesis, Fast monoenergetic neoclassical transport coefficients in stellarators.

This page documents the physical model that NTX actually implements, the normalizations that appear in the code, and the source files in which each piece is assembled.

Scope Of The Model

NTX solves the monoenergetic drift-kinetic equation with:

• fixed particle speed

• Lorentz pitch-angle scattering

• Boozer-angle spatial coordinates on one flux surface

• an externally specified radial electric field through epsi_hat or er_hat

• two right-hand sides:

  • the radial-transport drive

  • the parallel-conductivity / bootstrap-current drive

The present code therefore computes the monoenergetic geometric coefficients

\[\hat D_{11},\quad \hat D_{31},\quad \hat D_{13},\quad \hat D_{33}, \quad \hat D_{33,\mathrm{Sp}}.\]

These are the quantities exposed in src/ntx/transport.py and surfaced in src/ntx/solver.py.

NTX does not attempt to solve the full multi-species ambipolar problem by itself. That layer is delegated to downstream tools such as NEOPAX, which consume the NTX monoenergetic tables.

Coordinate And Field Representation

On a single flux surface, the magnetic-field strength is represented in Boozer angles as

\[B(\theta,\zeta) = \sum_{m,n} B_{mn} \cos\!\left(m\theta + n N_\mathrm{fp}\zeta\right)\]

for stellarator-symmetric surfaces, or with additional sine coefficients when those are provided.

This is implemented in:

• src/ntx/geometry.py: evaluate_fourier_series(...)

• src/ntx/booz.py: Boozer-file loading

• src/ntx/vmec.py: VMEC wout loading through vmec_jax

The spatial geometry object stored on the angular grid is GeometryOnGrid, defined in src/ntx/geometry.py.

Geometric Quantities Used By The Solver

For Boozer-surface inputs, NTX constructs

\[\mathcal J = \frac{|B_\zeta + \iota B_\theta|}{B^2}, \qquad B^\theta = \frac{\iota}{\mathcal J}, \qquad B^\zeta = \frac{1}{\mathcal J}.\]

For VMEC-derived surfaces, NTX evaluates Fourier series for

• B

• \mathcal J

• B_\theta

• B_\zeta

• B^\theta

• B^\zeta

directly on the angular grid.

The magnetic-drift factor that enters the transport source is

\[\hat v_m(\theta,\zeta) = \frac{B_\theta \,\partial_\zeta B - B_\zeta \,\partial_\theta B} {\mathcal J B^3}.\]

This quantity is assembled in src/ntx/geometry.py inside _boozer_geometry_on_grid(...) and _vmec_geometry_on_grid(...).

Drift-Kinetic Equation Treated By NTX

NTX solves the local monoenergetic drift-kinetic equation for the non-adiabatic response on a fixed flux surface,

\[\left[ \xi \frac{1}{B}\left(B^\theta\partial_\theta+B^\zeta\partial_\zeta\right) + \frac{\hat E_\psi}{\mathcal J\langle B^2\rangle} \left(-B_\zeta\partial_\theta+B_\theta\partial_\zeta\right) - \frac{1-\xi^2}{2B^2} \left(B^\theta\partial_\theta B+B^\zeta\partial_\zeta B\right)\partial_\xi - C_L \right] f = s.\]

Here \xi = v_\parallel / v is pitch angle and C_L is the Lorentz pitch-angle-scattering operator,

\[C_L f = \frac{\hat\nu}{2} \partial_\xi\left[(1-\xi^2)\partial_\xi f\right].\]

The model is monoenergetic: speed is fixed, energy scattering is not retained, and the full species/profile problem is outside this local solve. NTX treats two source systems, one for radial transport and one for parallel flow. Those monoenergetic responses are later integrated or mapped by profile and NEOPAX workflows when a species-resolved bootstrap-current calculation is needed.

Monoenergetic Legendre System

The unknown distribution is expanded in Legendre modes of pitch angle,

\[f(\theta,\zeta,\xi) = \sum_{k=0}^{N_\xi} f^{(k)}(\theta,\zeta) P_k(\xi).\]

NTX then solves the block-tridiagonal hierarchy

\[L_k f^{(k-1)} + D_k f^{(k)} + U_k f^{(k+1)} = s^{(k)}.\]

This is the exact operator form implemented in src/ntx/operators.py and solved in src/ntx/solver.py.

The lower, diagonal, and upper blocks contain:

\[\begin{split}L_k &= \frac{k}{(2k-1)B} \left( B^\theta \partial_\theta + B^\zeta \partial_\zeta \right) \;+\; \frac{k(k-1)}{2(2k-1)B^2} \left(B^\zeta \partial_\zeta B + B^\theta \partial_\theta B\right), \\ D_k &= \frac{\hat E_\psi}{\mathcal J \langle B^2\rangle} \left(-B_\zeta \partial_\theta + B_\theta \partial_\zeta\right) \;+\; \frac{\hat \nu}{2}k(k+1), \\ U_k &= \frac{k+1}{(2k+3)B} \left( B^\theta \partial_\theta + B^\zeta \partial_\zeta \right) \;+\; \frac{-(k+1)(k+2)}{2(2k+3)B^2} \left(B^\zeta \partial_\zeta B + B^\theta \partial_\theta B\right).\end{split}\]

The exact coefficient arrays corresponding to these formulas are constructed in coefficients_for_k(...) in src/ntx/operators.py.

Source Systems

NTX solves two linear systems per monoenergetic case:

1. a transport-driven system used for D11 and D31

2. a parallel-current / conductivity system used for D13 and D33

In code these are the source arrays s1 and s3 from source_modes(...) in src/ntx/operators.py.

Their Legendre support is sparse:

• s1 occupies modes k = 0 and k = 2

• s3 occupies mode k = 1

That structure is central to the NTX solver design: only the low-order modes k = 0,1,2 are required in the backward substitution, even when the forward Schur elimination proceeds to much larger N_\xi.

Nullspace Constraint

The monoenergetic system has the usual constant-mode nullspace in the k = 0 equation. NTX removes it by replacing the first spatial row with

\[f^{(0)}(\theta_0,\zeta_0)=0.\]

That operation is implemented by apply_nullspace_condition(...) in src/ntx/operators.py.

Transport Coefficients

Once the low-order Legendre modes are known, NTX evaluates the monoenergetic coefficients through flux-surface averages. In the code:

\[\hat D_{11} = \frac{\left\langle -2 \hat v_{m,0} f_{1,0} -\frac{2}{5}\hat v_{m,2} f_{1,2} \right\rangle}{\Psi_\mathrm{coeff}^2},\]

\[\hat D_{31} = \frac{\left\langle \frac{2}{3} B f_{1,1}/B_0 \right\rangle} {\Psi_\mathrm{coeff}}, \qquad \hat D_{13} = \frac{\left\langle -2 \hat v_{m,0} f_{3,0} -\frac{2}{5}\hat v_{m,2} f_{3,2} \right\rangle} {\Psi_\mathrm{coeff} B_0},\]

\[\hat D_{33} = \frac{\left\langle \frac{2}{3} B f_{3,1}/B_0 \right\rangle}{B_0}.\]

NTX also computes the Spitzer reference contribution

\[\hat D_{33,\mathrm{Sp}} = \frac{2}{3\hat \nu} \left\langle \frac{B^2}{B_0^2} \right\rangle,\]

which is the expression quoted in the thesis appendix and implemented in coefficients_from_modes(...) in src/ntx/transport.py.

Onsager Symmetry

The monoenergetic formulation satisfies the expected Onsager relation

\[\hat D_{13} = -\hat D_{31}\]

under the assumptions of the model. NTX tracks the scalar residual

\[|\hat D_{13} + \hat D_{31}|\]

through onsager_error(...) in src/ntx/transport.py.

This quantity is written into the CLI NetCDF/NPZ/HDF5 outputs and exposed in the high-level TransportResult.

Electric-Field Normalization

NTX accepts either epsi_hat directly or er_hat together with a surface-dependent normalization.

For Boozer / DKES-style surfaces,

\[\hat E_\psi = \hat E_r / \psi_p.\]

For VMEC surfaces, NTX computes

\[r_n = \sqrt{\psi_n}, \qquad \hat r = a_\mathrm{minor} r_n, \qquad \frac{d\hat \psi}{d\hat r} = \frac{2 \hat \psi_a r_n}{a_\mathrm{minor}}, \qquad \hat E_\psi = \hat E_r \left(\frac{d\hat \psi}{d\hat r}\right)^{-1}.\]

This conversion is implemented in:

• src/ntx/vmec.py for file-based VMEC inputs

• src/ntx/solver.py in MonoenergeticCase.resolved_epsi_hat(...)

What Downstream Tools Consume

Downstream bootstrap-current or ambipolar-transport tools typically consume D11, D13, and D33 over radial, collisionality, and electric-field scans.

NTX provides that through:

• solve_monoenergetic_scan(...)

• build_ntx_neopax_scan(...)

• scan_to_neopax_arrays(...)

in src/ntx/solver.py and src/ntx/neopax.py.

The public database bridge defaults to the raw D33 convention because that is the branch that preserves the integrated W7-X transfer gate. The D33_spitzer branch remains available as an explicit fixed-field closure stress option, but it is not a global runtime default.

Database Bridge And Bootstrap-Current Observable

The NTX-to-NEOPAX bridge is a normalization bridge between two well-defined physics objects:

1. NTX computes monoenergetic geometric coefficients on a flux surface.

2. The downstream energy convolution evaluates species transport matrices from those coefficients and the thermodynamic forces.

3. The momentum-restoring closure solves a finite Sonine moment system for the corrected parallel-flow moments.

4. The bootstrap-current density is assembled from the species parallel flows.

No fitted scalar is used in this chain.

Thermodynamic-Force Form

For species a, the reduced closure uses the linear-response form

\[\begin{split}\begin{pmatrix} \Gamma_a/n_a \\ Q_a/(n_a T_a) \\ U_{\parallel a}/n_a \end{pmatrix} = - \begin{pmatrix} L_{11,a} & L_{12,a} & L_{13,a} \\ L_{21,a} & L_{22,a} & L_{23,a} \\ L_{31,a} & L_{32,a} & L_{33,a} \end{pmatrix} \begin{pmatrix} A_{1,a} \\ A_{2,a} \\ A_3 \end{pmatrix}.\end{split}\]

Here A1 is the density/electrostatic drive, A2 is the temperature-gradient drive, and A3 is the parallel-electric-field drive. This separation matters for validation: the fixed-field bootstrap current is not determined by a single monoenergetic coefficient, but by how the density and temperature drives enter the row-3 response and the momentum-restoring Sonine solve.

For the finite-beta Redl stress audit, the analytic target is also stored as source terms on the same current observable,

\[\langle J\cdot B\rangle_\mathrm{Redl} = J_n + J_{T_e} + J_{T_i},\]

where J_n is the Redl density-gradient contribution and J_Te + J_Ti is the effective temperature-gradient target. The source-channel audit compares these terms to the frozen NTX+NEOPAX density/electric and effective-temperature responses without inserting any fitted multiplier into the runtime closure.

The no-momentum parallel-flow branch is therefore

\[U_{\parallel a}^{(0)} = -n_a\left(L_{31,a} A_{1,a} +L_{32,a} A_{2,a} +L_{33,a} A_3\right).\]

The momentum-restoring branch solves for Sonine coefficients c_{ak}. In the current finite basis the observable is fixed by the moment definition:

\[U_{\parallel a}=n_a c_{a0}.\]

This is why NTX keeps the observable map fixed in validation. Changing a fitted linear combination of c_{a0}, c_{a1}, ... can improve one benchmark but would change the physical moment definition. Any accepted closure improvement has to change the projected moment equations or collision blocks, not the current observable by a benchmark-tuned remap.

Monoenergetic Coefficient Scaling

The database bridge follows the dimensional role of each coefficient in the energy convolution:

\[D_{11}^{\mathrm{db}} = D_{11}^{\mathrm{NTX}}\left(\frac{dr}{ds}\right)^2, \qquad D_{13}^{\mathrm{db}} = D_{13}^{\mathrm{NTX}}\frac{dr}{ds}, \qquad D_{33}^{\mathrm{db}} = \nu D_{33}^{\mathrm{NTX}}.\]

The first factor appears twice in D11 because radial particle/heat transport contains two radial-gradient normalizations. The mixed D13 coefficient contains one radial and one parallel response, so it receives one dr/ds factor. The parallel conductivity coefficient is stored as nu * D33 because the consumer reconstructs the monoenergetic conductivity integrand by dividing by the collisionality used in the energy convolution. This is the exact mapping implemented in src/ntx/_neopax_bridge.py.

D33_spitzer is not a fitted alternate normalization. It is the analytic Spitzer-conductivity contribution

\[\hat D_{33,\mathrm{Sp}} = \frac{2}{3\hat \nu} \left\langle\frac{B^2}{B_0^2}\right\rangle,\]

which is the constant-field conductivity limit of the monoenergetic system. The fixed-field QA/QH stress audit can select this branch explicitly because the low-order momentum-restoring closure is a parallel-flow conductivity model. The public bridge still defaults to raw D33, because the integrated W7-X database transfer is validated on that convention.

Current Assembly And SFINCS Observable

The species current density from the closure is

\[j_{\parallel a}=Z_a e\,U_{\parallel a}.\]

The archived fixed-field reference stores flux-surface-averaged parallel-flow observables. In the SFINCS output convention,

\[\mathrm{FSABjHat} = \sum_a Z_a\,\mathrm{FSABFlow}_a, \qquad \mathrm{FSABjHatOverB0} = \frac{\mathrm{FSABjHat}}{\mathrm{B0OverBBar}}.\]

The comparison therefore uses the archived B0OverBBar and the imported Boozer handedness convention to convert the closure current into the same flux-surface-averaged observable before comparing with the archived profile. Both pieces are read from the reference output or from the imported geometry; neither is fitted to the QA/QH current.

The other important implementation point is semantic rather than numerical: get_Neoclassical_Fluxes_With_Momentum_Correction(...) returns the total corrected U_parallel. It is not an increment to add to the no-momentum U_parallel. Treating it as an increment double-counts the no-momentum branch and produces an unphysical current. The shipped examples and validation scripts now use the corrected return directly.

Why The Fixed-Field Gate Is A Stress Gate

The precise-QS fixed-field comparison now passes the total-current stress gate using only physics-derived operations:

• exact NTX coefficient normalizations,

• explicit D33_spitzer conductivity branch for the fixed-field stress model,

• SFINCS observable conversion from archived B0OverBBar,

• total corrected U_parallel semantics,

• no fitted constants and no species-dependent bridge factors.

It is still not promoted as species-current parity. Continuum drift-kinetic solvers with fuller linearized Fokker-Planck physics can resolve species-current components that a low-order momentum-restoring closure built from monoenergetic coefficients is not guaranteed to reproduce coefficient by coefficient. The promoted statement is therefore deliberately narrower: the total bootstrap-current profile is inside the documented fixed-field stress gate while the independent Redl comparison and the integrated W7-X raw-branch transfer remain separate validation gates.

Closure-Model Gates

The monoenergetic coefficient pipeline and the downstream momentum-restoring closure have to be validated separately.

For the closure side, the physically relevant gates are:

1. Coefficient-side invariants stay closed first. The monoenergetic solver must continue to satisfy the Onsager relation and the established database normalization bridge before any closure change is considered.

2. The observable map stays fixed. For the present Sonine/Laguerre basis, the corrected parallel flow is

   \[U_{\parallel a} = n_a c_{a0},\]

   so higher-order closure work must change the solved moment system, not the final observable by an empirical remap.

3. Finite-order closure must recover the current three-moment system exactly at P=2. Any arbitrary-order implementation has to reproduce the legacy P=2 Sonine basis, normalization factors, and source projections before it is used to interpret benchmark differences.

4. The closure must preserve symmetry properties implied by the projected theory. In practice this means retaining the Onsager/ambipolar structure emphasized by Sugama–Horton and Sugama–Nishimura, rather than inserting benchmark-fit constants into selected matrix entries.

5. The projected collision model must preserve the collisional invariants. Particle number, parallel momentum, and energy must remain exact invariants of the finite-order collisional system, and the weighted collisional form should retain the self-adjoint structure that underlies Onsager symmetry and non-negative entropy production.

6. Benchmark transfer is mandatory. A closure change is only acceptable if it improves the precise-QS fixed-field QA/QH current benchmark without regressing the already-validated integrated W7-X workflow.

These gates are now the design constraints for any higher-order closure work in the imported NTX+NEOPAX path. The runnable summary of the artifact-backed gates is in physics-gates.md and python scripts/check_physics_gates.py.