Grok Challenge #1, Paper 1 - Intrafere Research Group

← PAPER DIRECTORY

This is the first paper and rough draft answer generated for @Grok’s 2/28/2026 challenge to MOTO ASI by Intrafere Research Group. A fully autonomous super intelligence AI harness. This paper cost around $15 to produce and took approximately two hours. As of the time of this writing, MOTO is continuing to solve Grok’s challenge with successive papers, attempting recursive improvement and solution-set diversification. Company Website: https://intrafere.com/ Software GitHub that produced this paper: https://github.com/Intrafere/MOTO-Autonomous-ASI Grok Fusion Solution Challenge Link: https://x.com/grok/status/2027657401625690332 ================================================================================ AUTONOMOUS AI SOLUTION Disclaimer: This is an autonomous AI solution generated with the MOTO harness. This paper was not peer reviewed and was autonomously generated without user oversight or interaction beyond the original user prompt, therefore, this text may contain errors. These papers often contain ambitious content and/or extraordinary claims, all content should be viewed with extreme scrutiny. (EDITOR NOTE: This single paper does not attempt to solve the user’s prompt entirely, it is meant to be one piece toward the complex solution required for the users prompt – this paper is the first paper in the series – total solutions typically are achieved in later papers) User’s Research Prompt: Deliver a complete, engineering-ready blueprint for a compact stellarator fusion reactor achieving sustained Q>15 net gain by 2030—using only near-term materials, full MHD/plasma stability models, tritium breeding cycle, and <$5B build cost. Include all equations, sim code, and falsifiable tests. AI Model Authors: x-ai/grok-4, openai/gpt-5.2, z-ai/glm-5, moonshotai/kimi-k2.5 Generated: 2026-02-28 ================================================================================ TITLE: A Falsifiable 0D\\/+1D Systems Model for Compact Q>15 Stellarators: Transport-Operator Closures and Sensitivity Analysis Abstract We construct a deterministic, modular, and explicitly falsifiable 0D/+1D systems map for compact-stellarator performance studies aimed at high fusion gain, formulated as a residual problem $\mathcal{R}(\mathbf{u},\mathbf{p})=0$ with a separate feasibility-constraint layer $\mathcal{C}(\mathbf{u},\mathbf{p})\le 0$. The state $\mathbf{u}$ contains flux-surface-averaged profiles (at minimum $T_e(\rho),T_i(\rho)$, optionally $n_e(\rho)$) and derived scalars $(W,\tau_E,Q,\ldots)$; the output is required to include intermediate observables (profiles, fluxes, integrated power channels, margins) so that incorrect closures can be rejected on existing devices. The +1D backbone is a conservative two-temperature steady energy balance with a unit-consistent diffusive closure that admits an explicit quadrature representation and thus a direct “infer $\chi_{exp}(\rho)$” falsifier. A fast global closure is derived as a scalar root problem $P-\mathcal{L}(P)=0$ with checkable existence/uniqueness conditions, enabling power-scan falsification. The same transport structure is recast as a self-adjoint Sturm–Liouville operator with Robin edge loss, producing additional falsifiable observables: the principal decay rate $\lambda_1$, Rayleigh-quotient bounds, monotonicity tests, and modulation transfer functions. We specify swap-ready physics extensions (pinch, ambipolar $E_r$, energetic-particle deposition, particle/radiation/density-limit closures, bootstrap/island penalties, stability and SOL boundaries), a positivity-preserving finite-volume discretization with mandatory conservation diagnostics, exact adjoint sensitivities (including eigenvalue/resolvent derivatives), and a calibrate-or-reject statistical/OED wrapper. The contribution is a mathematically explicit framework and near-term rejection protocol, not a validated claim of $Q>15$. I. Introduction I.A. Problem setting: compact stellarator design loop with a falsifiable 0D/+1D model Compact stellarator concepts targeting high fusion gain (e.g. the informal goal $Q>15$) motivate end-to-end design loops that must evaluate many candidate configurations and operating points. In such loops, purely empirical scalar confinement scalings are fast but often under-determine *why* a prediction succeeds or fails, while high-fidelity equilibrium, transport, and kinetic calculations can be too expensive to support rapid iteration. This paper develops a middle-ground mathematical object: a deterministic coupled map from a parameter vector $\mathbf{p}$ (geometry proxies, operating set points, and closure coefficients) to a state $\mathbf{u}$ (profiles and derived scalars) and a vector of predicted observables $\widehat{\mathbf{y}}$. The central methodological constraint is *falsifiability on existing devices*. Because reactor-regime data do not exist for the target parameter range, the model must produce intermediate observables that can be compared to present stellarator diagnostics, so that incorrect closure assumptions can be rejected (or localized and replaced) before extrapolation. I.B. Design goals and falsifiability contract The document formalizes four design goals that guide the construction. 1. Fast end-to-end evaluation: the core closure is written in conservative 1D form and paired with reductions (including a scalar root formulation in a power-law operating neighborhood) that enable rapid steady evaluation. 2. Intermediate observables: the output $\widehat{\mathbf{y}}$ is required to include not only $Q$ but also profiles, fluxes, integrated power channels, constraint margins, and (when invoked) operator-level decay and modulation quantities. 3. Near-term rejection tests: each modeling choice is paired with explicit measurable comparisons (power scans, profile-based transport inference, modulation transfer functions, vacuum-field metrics, etc.) so that failure is detectable within realistic experimental time horizons. 4. Swap-ready upgrades: the coupled map is assembled from modules with declared interfaces, so that a component (equilibrium/geometry, neoclassical coefficients, turbulence surrogates, radiation tables, edge/SOL closure) can be replaced without altering the surrounding wrapper, the residual bookkeeping, or the statistical rejection protocol. These goals are encoded as a falsifiability contract (Section II): the steady coupled model is posed as a residual equation $\mathcal{R}(\mathbf{u},\mathbf{p})=0$ with a separate inequality constraint layer $\mathcal{C}(\mathbf{u},\mathbf{p})\le 0$, and every module is required to expose auditable diagnostics and (when used for optimization or OED) compatible derivatives. I.C. Central mathematical themes Three mathematical themes structure the analysis. 1. Transport closures in conservative +1D form. The core backbone is a flux-surface averaged two-temperature steady energy balance on $\rho\in[0,1]$, with explicit unit conventions and a conservative finite-volume discretization that enforces discrete power-balance identities and positivity requirements. 2. Transport-operator (Sturm–Liouville) closures. The same diffusive structure is re-expressed as a self-adjoint weighted operator $\mathcal{L}$ with Robin edge loss, producing additional observables: a principal decay rate $\lambda_1$, Rayleigh-quotient bounds, monotonicity properties, and resolvent-based modulation transfer functions $\Theta(\rho,\omega)$. These quantities provide falsifiers that are distinct from steady profile fits. 3. Exact sensitivity derivatives. The coupled discrete system is written in a form compatible with the implicit function theorem and adjoint differentiation, yielding exact gradients of objectives and constraint margins with respect to many parameters. Operator-level sensitivities (eigenvalue and transfer-function derivatives) are recorded in forms that depend only on module-supplied coefficient derivatives, consistent with the swap-ready requirement. I.D. Roadmap and organization Section II defines the coupled systems map and the calibrate-or-reject philosophy, including the explicit choice that scalar performance metrics are *derived* from a solved state rather than treated as primitive fits. Section III fixes notation, units, and the bookkeeping of stored energy, power channels, $\tau_E$, and $Q$. Section IV specifies a reduced (0D) equilibrium/geometry parameterization together with vacuum-field falsification hooks. Section V introduces the steady +1D transport backbone, including a direct quadrature representation for diffusion-only closures and a fast scalar-root closure under stated continuity/monotonicity hypotheses. Section VI develops the transport-operator formulation, bounds, monotonicity, and modulation predictions. Section VII lists extension modules needed for credible high-gain projections (pinch/advection, ambipolar $E_r$, energetic-particle deposition, particle balance and composition, radiation, density-limit surrogate, bootstrap/island penalties, stability constraints, and a minimal edge/SOL closure). Sections VIII–IX add a time-dependent 0D thermal stability wrapper (with a falsifiable small-signal eigenvalue $\gamma_T$ and an $\alpha$-heating emulation protocol) and a conservative positivity-preserving discretization with mandatory numerical diagnostics. Sections X–XI present the exact sensitivity/adjoint framework and a statistical calibration, model-discrimination, and optimal-experimental-design wrapper (including explicit measurement operators and a similarity-law falsifier). Section XII consolidates these elements into an integrated near-term falsification plan with module-level rejection logic, and Section XIII records limitations and regime-of-validity risks that must be externally verified before any extrapolative use. The overarching intent is not to assert verified $Q>15$ performance, but to provide a mathematically explicit, interface-stable, and rejection-oriented 0D/+1D framework in which each closure yields intermediate, testable predictions and can be replaced without collapsing the end-to-end structure. II. Model Architecture and Falsifiability Contract II.A. Definition of the coupled systems map We formalize the systems model as a parameter-to-observable map obtained by solving a coupled steady system. 1. Parameter vector. Let $\mathbf{p}\in\mathcal{P}\subset\mathbb{R}^m$ denote all inputs, partitioned as \[ \mathbf{p}=(\mathbf{p}_{geom},\mathbf{p}_{op},\boldsymbol{\theta}), \] where $\mathbf{p}_{geom}$ collects geometry/equilibrium proxies (e.g. a reduced boundary surface description and rotational-transform proxy parameters), $\mathbf{p}_{op}$ collects operational set-points (e.g. magnetic field, heating powers, boundary temperatures/densities, isotope mix assumptions), and $\boldsymbol{\theta}$ collects closure coefficients (transport, radiation, constraint margins) that are not directly measured. 2. State variables. Let $\rho\in[0,1]$ be a normalized flux label. The state $\mathbf{u}$ consists of (i) profile functions and (ii) scalars: \[ \mathbf{u}=\bigl(T_e(\rho),T_i(\rho),n_e(\rho),\ldots;\; W,\tau_E,Q,\ldots\bigr). \] In computation, profiles are represented on a grid $\{\rho_j\}_{j=0}^N$, so $\mathbf{u}\in\mathbb{R}^d$ for some $d$ after discretization. 3. Residual form. The coupled steady model is written as a residual equation \[ \mathcal{R}(\mathbf{u},\mathbf{p})=\mathbf{0},\qquad \mathcal{R}:\mathcal{U}\times\mathcal{P}\to\mathbb{R}^d, \] which compactly encodes the core transport balance, power balance, and any auxiliary closure relations used to evaluate sinks/sources. Constraints that represent feasibility limits are represented separately as inequalities (defined below). 4. Systems solution operator and observables. Whenever $\mathcal{R}(\cdot,\mathbf{p})=0$ has an admissible solution $\mathbf{u}(\mathbf{p})$, we define the solution operator \[ \mathcal{M}:\mathcal{P}\to\mathcal{U},\qquad \mathcal{M}(\mathbf{p})=\mathbf{u}(\mathbf{p}), \] and the observable prediction map \[ \mathcal{O}:\mathcal{U}\times\mathcal{P}\to\mathbb{R}^k,\qquad \widehat{\mathbf{y}}=\mathcal{O}(\mathbf{u}(\mathbf{p}),\mathbf{p}). \] The end-to-end prediction map is then \[ \mathcal{F}:\mathcal{P}\to\mathbb{R}^k,\qquad \mathcal{F}(\mathbf{p}) := \mathcal{O}(\mathcal{M}(\mathbf{p}),\mathbf{p}). \] Admissibility is enforced by a state constraint set $\mathcal{U}_{phys}\subset\mathcal{U}$ (e.g. positivity of temperatures and densities, monotonicity conditions when imposed, and boundedness of derived quantities). A parameter point $\mathbf{p}$ is deemed infeasible if no $\mathbf{u}\in\mathcal{U}_{phys}$ solves the residual. II.B. Modular decomposition (\u201cswap-ready\u201d interfaces) The architecture is defined as a composition of modules, each being a mathematical map with declared inputs/outputs. The key requirement is that a module can be replaced by another implementation with the same interface without altering the rest of the system. II.B.1. Geometry/equilibrium proxy module The geometry module is a map \[ \mathcal{G}:\mathbf{p}_{geom}\mapsto \mathbf{g}(\rho), \] where $\mathbf{g}(\rho)$ denotes flux-surface-averaged geometric factors required by transport and constraints. Examples include the volume derivative $V'(\rho)$, characteristic metric coefficients required to relate $d/d\rho$ to a physical radial coordinate, and neoclassical/turbulence control proxies such as an effective ripple-like function $\epsilon_{eff}(\rho)$ or related measures. The interface contract is purely mathematical: any upgraded equilibrium calculation is acceptable if it produces $\mathbf{g}(\rho)$ in consistent units and with consistent normalization. II.B.2. Core transport and power-balance module Given geometry factors and operational inputs, the transport module produces either (i) profile solutions directly or (ii) transport coefficients and a residual operator for the coupled profiles. Abstractly, \[ \mathcal{T}:(\mathbf{g},\mathbf{p}_{op},\boldsymbol{\theta}_T)\mapsto \mathcal{R}_T, \] where $\mathcal{R}_T$ is the portion of the residual encoding the steady flux-surface-averaged energy balances. The power-balance component evaluates volumetric sources and sinks and defines scalar performance metrics as derived quantities: \[ \mathcal{P}:(\mathbf{u},\mathbf{p})\mapsto (W,\tau_E,Q,\ldots). \] The model is organized so that (a) $Q$ is not a primitive but a derived function of a solved state, and (b) intermediate quantities entering $Q$ (e.g. profile-dependent losses) remain explicit observables. II.B.3. Constraint layer as inequalities and diagnostics Feasibility limits are collected as inequality constraints \[ \mathcal{C}(\mathbf{u},\mathbf{p})\le 0 \] where $\mathcal{C}:\mathcal{U}\times\mathcal{P}\to\mathbb{R}^\ell$. These represent limits such as density caps, stability margins, stochasticity indicators, engineering bounds, or edge/SOL heat flux limits. Each component $\mathcal{C}_i$ must be accompanied by a diagnostic target, i.e. a measurable proxy that permits empirical testing of whether $\mathcal{C}_i$ is correctly triggered. The role of the constraint layer is twofold: (i) It filters out parameter points even if the residual equations admit a solution (\u201csolution exists but is not credible\u201d). (ii) It produces explicit, separately falsifiable predicted quantities (e.g. a predicted limit threshold or a predicted margin to threshold). II.C. Falsifiability objects beyond scalar Q The falsifiability contract requires that the output $\widehat{\mathbf{y}}$ include intermediate observables that can be compared to existing diagnostics without burning-plasma operation. We therefore require $\widehat{\mathbf{y}}$ to contain, as available, 1. Profile-level observables: $T_e(\rho_j),T_i(\rho_j),n_e(\rho_j)$ on a prescribed grid, and boundary values (edge temperatures/densities) when treated as predicted rather than prescribed. 2. Flux-level observables: radial heat fluxes $q_e(\rho_j),q_i(\rho_j)$ and integrated conductive losses at the boundary, together with a global steady residual diagnostic \[ \mathfrak{r}_{PB}:=\bigl|P_{in}-P_{out}\bigr|, \] where $P_{in}$ and $P_{out}$ are computed from the same $\mathbf{u}$ and must agree to within stated numerical and measurement tolerances. 3. Transport-coefficient observables: effective diffusivities $\chi_e(\rho),\chi_i(\rho)$ and any neoclassical/turbulent decomposition the chosen transport closure provides. These are included so that experimental transport inference (or modulation analysis) can be compared to model-implied transport. 4. Operator-level observables (optional but admissible): if the transport closure is expressed via a linear(ized) operator $\mathcal{L}$ acting on perturbations of profiles, then the model may output spectral or transfer-function summaries (e.g. a principal decay rate, frequency response amplitude/phase). Such quantities are treated as standard observables: they must be computed from declared inputs and are subject to the same rejection tests as profile data. 5. Constraint-related observables: each inequality constraint $\mathcal{C}_i(\mathbf{u},\mathbf{p})\le 0$ must be paired with a predicted measurable proxy $\widehat{z}_i$ (e.g. a threshold value or margin), and the output vector may include both $\mathcal{C}_i$ and $\widehat{z}_i$ for transparency. This contract is designed so that failure can be localized: a disagreement in a profile or flux can refute a transport closure even when the computed scalar $Q$ happens to match one operating point. II.D. Explicit \u201ccalibrate-or-reject\u201d philosophy We formalize calibration as a constrained estimation problem in a limited parameter subspace, followed by validation that forbids additional tuning. 1. Limited calibration set. Choose a designated subset of closure coefficients $\boldsymbol{\theta}_{cal}\subset\boldsymbol{\theta}$ and hold all other components of $\mathbf{p}$ fixed at measured or externally specified values. The admissible calibration set must be declared a priori and should be low-dimensional. Mathematically, we define a calibration subspace \[ \Theta_{cal}\subset\mathbb{R}^{m_{cal}},\qquad \boldsymbol{\theta}_{cal}\in\Theta_{cal}, \] with explicit bounds representing physical plausibility. 2. Data model and misfit. Let a calibration dataset be $\mathcal{D}_{cal}=\{(\mathbf{y}^{(k)},\Sigma^{(k)},\mathbf{p}^{(k)})\}_{k=1}^{N_{cal}}$, where $\mathbf{y}^{(k)}$ is a vector of measured observables, $\Sigma^{(k)}$ an uncertainty covariance (diagonal is acceptable if correlations are unavailable), and $\mathbf{p}^{(k)}$ the measured inputs for shot/condition $k$. The weighted least-squares objective is \[ J(\boldsymbol{\theta}_{cal}) := \sum_{k=1}^{N_{cal}} \bigl(\mathbf{y}^{(k)}-\mathcal{F}(\mathbf{p}^{(k)};\boldsymbol{\theta}_{cal})\bigr)^T \bigl(\Sigma^{(k)}\bigr)^{-1}\bigl(\mathbf{y}^{(k)}-\mathcal{F}(\mathbf{p}^{(k)};\boldsymbol{\theta}_{cal})\bigr), \] where $\mathcal{F}(\mathbf{p}^{(k)};\boldsymbol{\theta}_{cal})$ indicates the dependence of $\mathcal{F}$ on $\boldsymbol{\theta}_{cal}$ with all other parameters fixed. A calibrated parameter vector is any minimizer \[ \boldsymbol{\theta}_{cal}^*\in\operatorname*{arg\,min}_{\boldsymbol{\theta}_{cal}\in\Theta_{cal}} J(\boldsymbol{\theta}_{cal}). \] If the minimum is not attained (e.g. due to noncompactness) or is non-identifiable (flat directions), the module is deemed non-calibratable with the declared dataset and must be rejected or augmented with additional observables. 3. Rejection criterion (validation). For an independent validation dataset $\mathcal{D}_{val}$ (disjoint from calibration conditions), define the validation statistic \[ \chi^2_{val}:=\sum_{k\in\mathcal{D}_{val}} \bigl(\mathbf{y}^{(k)}-\mathcal{F}(\mathbf{p}^{(k)};\boldsymbol{\theta}_{cal}^*)\bigr)^T \bigl(\Sigma^{(k)}\bigr)^{-1}\bigl(\mathbf{y}^{(k)}-\mathcal{F}(\mathbf{p}^{(k)};\boldsymbol{\theta}_{cal}^*)\bigr). \] The model (or a specified module) is rejected if $\chi^2_{val}$ exceeds a pre-declared threshold consistent with the assumed noise model and degrees of freedom, or if systematic residual structure is observed in any individual falsifiability object (e.g. consistent sign bias in predicted heat flux profiles across multiple conditions). 4. Localizing failure to modules. To support module-level rejection, we require that each module expose a set of internal \u201cinterface observables\u201d $\widehat{\mathbf{y}}_{int}$ that are (i) predicted by the module and (ii) either measured directly or inferable from measurements with stated uncertainty. A module is rejected if its own $\widehat{\mathbf{y}}_{int}$ fails validation regardless of whether the aggregate $Q$ prediction happens to agree at one operating point. 5. Non-negotiable numerical sanity checks. The computational realization of $\mathcal{F}$ must explicitly report: (i) whether a physical solution exists, (ii) whether $\mathbf{u}\in\mathcal{U}_{phys}$, and (iii) the global residual measures associated with its discrete equations. Any violation of positivity or global conservation identities (beyond declared tolerances) is treated as an immediate reject condition for that computed point, not as a tunable effect. III. Preliminaries, Notation, and Core Performance Metrics III.A. Flux label and geometry conventions We use a normalized flux label $\rho\in[0,1]$ with $\rho=0$ at the magnetic axis and $\rho=1$ at the last closed flux surface (LCFS). For reduced-model purposes we adopt a simple physical minor-radius coordinate \[ r := a\rho, \] where $a$ is an effective minor radius (m). Volume enclosed by the surface $\rho$ is $V(\rho)$ (m$^3$), and we define the key geometry factor \[ V'(\rho) := \frac{dV}{d\rho}\quad [\mathrm{m}^3]. \] For any scalar field $f$ defined on a flux surface, the flux-surface average is denoted by $\langle f\rangle(\rho)$. No specific coordinate system is assumed here; the only required geometric input to the +1D backbone is that all radial integrals can be expressed in the form \[ \int_{\text{volume}} f\,dV = \int_0^1 \langle f\rangle(\rho)\,V'(\rho)\,d\rho. \] We also use the surface area $A(\rho)$ (m$^2$) when a flux is naturally expressed per unit area; if an implementation does not compute $A$ directly, it must still provide a consistent mapping between a radial heat flux density (W/m$^2$) and the corresponding integrated power crossing the surface (W). III.B. Species, profiles, and conventions We consider at minimum electrons $e$ and a single effective ion species $i$; extension to multiple ion/impurity species is made by replacing $i$ by an index set $s\in\mathcal{S}$ and enforcing quasineutrality. 1. Profile unknowns. The core +1D state is expressed in terms of flux functions \[ T_e(\rho),\;T_i(\rho)\;[\mathrm{keV}]\qquad\text{and}\qquad n_e(\rho)\;[\mathrm{m}^{-3}], \] and, when needed, additional profiles (impurity density, helium ash fraction, etc.). 2. Density treatment (prescribed vs solved). Two cases are permitted: (i) Prescribed density: $n_e(\rho)$ is an input profile (measured or assumed), and the energy equations are solved for $T_e,T_i$. (ii) Closed particle balance: $n_e(\rho)$ (or a parameterized family) is part of the state and is determined by auxiliary constraints (e.g. global fueling balance, particle confinement time). In either case, the computational module must declare explicitly which components of $\mathbf{u}$ are solved and which are prescribed. 3. Quasineutrality and effective ion density. In a multi-species setting with ion species $s$ of charge $Z_s$, quasineutrality is \[ n_e(\rho)=\sum_{s\in\mathcal{S}_i} Z_s\,n_s(\rho). \] When a single effective ion species is used, we denote its density by $n_i$ and assume $n_i\approx n_e$ for a singly-charged dominant ion with small impurity content; any deviation (e.g. dilution) must be declared in the parameter set $\mathbf{p}$ and propagated consistently into fusion and radiation terms. III.C. Stored energy, power channels, and gain III.C.1. Thermal stored energy The thermal energy content (J) is defined as \[ W := \int_0^1 \left(\frac{3}{2} n_e(\rho)\,T_e^{(J)}(\rho) + \frac{3}{2} n_i(\rho)\,T_i^{(J)}(\rho)\right) V'(\rho)\,d\rho, \] where $T_s^{(J)}$ denotes temperature expressed as an energy per particle (J). (Equivalently, if one works with thermodynamic temperature in Kelvin, then $T_s^{(J)}=k_B T_s^{(K)}$.) In implementations that use $T$ in keV, we use the conversion \[ T^{(J)} = C_{keV}\,T\,[\mathrm{keV}],\qquad C_{keV}=10^3 e=1.602176634\times 10^{-16}\ \mathrm{J/keV}, \] so that \[ W = \frac{3}{2}\,C_{keV}\,\int_0^1 \bigl(n_e T_e + n_i T_i\bigr)\,V'(\rho)\,d\rho, \] when $n$ is in m$^{-3}$ and $T$ in keV. III.C.2. Radial conductive (transport) power Let $q_s(\rho)$ denote the outward radial heat flux density of species $s$ (W/m$^2$) associated with the chosen closure. The corresponding integrated conductive power flowing outward across surface $\rho$ is \[ P_{cond,s}(\rho) := A(\rho)\,q_s(\rho)\quad [\mathrm{W}], \] and the total conductive loss from the confined region is \[ P_{cond} := \sum_s P_{cond,s}(1). \] Any equivalent definition is acceptable provided it is consistent with the discrete conservation identity used by the transport solver (i.e. the divergence of the modeled flux equals the net volumetric source/sink in steady state). III.C.3. Radiation and other sinks We denote the total radiated power as \[ P_{rad} := \int_0^1 p_{rad}(\rho)\,V'(\rho)\,d\rho\quad [\mathrm{W}], \] where $p_{rad}(\rho)$ is the volumetric radiated power density (W/m$^3$) arising from bremsstrahlung and line radiation (or any surrogate thereof). Additional non-conductive sinks (charge exchange, mechanical work terms neglected in the reduced model, etc.) are collected into \[ P_{other} := \int_0^1 p_{other}(\rho)\,V'(\rho)\,d\rho\quad [\mathrm{W}], \] with the understanding that the model must either (i) provide an explicit closure for $p_{other}$ or (ii) declare it as a bounded, externally estimated residual term to be tracked as an explicit diagnostic. III.C.4. Heating, fusion power, and Q We distinguish between injected auxiliary power $P_{aux}$ and absorbed power $P_{aux,abs}$ (W). The absorbed auxiliary power density is $p_{aux,abs}(\rho)$ (W/m$^3$), with \[ P_{aux,abs} = \int_0^1 p_{aux,abs}(\rho)\,V'(\rho)\,d\rho. \] Fusion power is \[ P_{fus} := \int_0^1 p_{fus}(\rho)\,V'(\rho)\,d\rho, \] where $p_{fus}$ is computed from the assumed fuel composition and the reactivity $\langle\sigma v\rangle$. For a D–T mixture with ion densities $n_D,n_T$, one common form is \[ p_{fus}(\rho) = n_D(\rho)\,n_T(\rho)\,\langle\sigma v\rangle(T_i(\rho))\,E_{fus}, \] with $E_{fus}$ the energy per fusion reaction (J). The alpha-heating power is a fixed fraction of fusion power, \[ P_{\alpha} = f_{\alpha}\,P_{fus},\qquad f_{\alpha}:=\frac{E_{\alpha}}{E_{fus}}\approx \frac{3.5}{17.6}\approx 0.20, \] and we allow a deposited fraction $\eta_{\alpha,dep}\in[0,1]$ so that \[ P_{\alpha,dep} := \eta_{\alpha,dep}\,P_{\alpha}. \] The total heating power entering the core energy balance is \[ P_{heat} := P_{aux,abs}+P_{\alpha,dep}. \] The fusion gain is defined as \[ Q := \frac{P_{fus}}{P_{aux,abs}}, \] with the convention that $Q$ is reported only when $P_{aux,abs}>0$ and the solution is feasible. (For ignition-like conditions, one may equivalently use $P_{\alpha,dep}$ in the denominator or report a separate ignition margin; the chosen convention must be declared and used consistently.) III.C.5. Energy confinement time and steady power balance We define the energy confinement time by the usual global ratio \[ \tau_E := \frac{W}{P_{loss}},\qquad P_{loss}:=P_{cond}+P_{rad}+P_{other}. \] In steady state, the energy balance requires \[ P_{heat} = P_{loss}, \] and any computed solution must satisfy this identity to within the declared numerical tolerance (in addition to satisfying the discretized profile equations). For transient comparisons (e.g. in experimental datasets), the corresponding identity becomes \[ \frac{dW}{dt} = P_{heat}-P_{loss}, \] and the model must state explicitly whether it is being used in steady mode ($dW/dt=0$) or as a component inside a time-dependent wrapper. III.D. Unit-consistency requirements and dimensional analysis checks All modules in the coupled map must obey the following unit conventions (or provide explicit conversion factors at the interface): 1. Base units. – $\rho$ dimensionless, $r$ in m. – Densities $n$ in m$^{-3}$. – Temperatures $T$ in keV (with the conversion constant explicitly used when computing energies in J). – Geometry factor $V’$ in m$^3$, area $A$ in m$^2$. – Volumetric powers $p(\rho)$ in W/m$^3$, integrated powers $P$ in W. 2. Transport coefficients. If a diffusive closure is used later of the form $q_s=-n_s\chi_s\,dT_s/dr$ (or its $\rho$-form equivalent), then $\chi_s$ must have units m$^2$/s, and the model must ensure that the combination $n\chi\,dT/dr$ produces W/m$^2$ after the keV-to-J conversion. 3. Global consistency check (mandatory diagnostic). Given computed $\mathbf{u}$, each evaluation must report the dimensionless global residual \[ \mathfrak{r}_{PB} := \frac{|P_{heat}-P_{loss}|}{\max(P_{heat},P_{loss},P_{ref})}, \] where $P_{ref}>0$ is a fixed small reference power inserted to avoid division by zero in low-power tests. A computation that reports $Q$ must satisfy $\mathfrak{r}_{PB}\le\varepsilon_{PB}$ for a declared tolerance $\varepsilon_{PB}$ (numerical, not fitted). 4. Dimensional sanity checks. The following implications must hold in any physically admissible state: – If $n\ge 0$ and $T\ge 0$, then $W\ge 0$. – If $P_{aux,abs}\ge 0$, $P_{\alpha,dep}\ge 0$, and all sink terms are computed nonnegative by construction, then $P_{heat}\ge 0$, $P_{loss}\ge 0$, and $\tau_E$ is well-defined whenever $P_{loss}>0$. Any violation indicates either an interface unit inconsistency or a numerical failure and is treated as infeasibility of that computed point. IV. Magnetic Equilibrium and Geometry Parameterization (0D Backbone) This section specifies a reduced, explicitly parameterized representation of the magnetic configuration used by the systems map. The intent is not to compute a full 3D MHD equilibrium, but to (i) define a consistent set of geometry parameters $\mathbf{p}_{geom}$ and (ii) define the derived geometric factors $\mathbf{g}(\rho)$ required by the +1D transport backbone and feasibility constraints. IV.A. Boundary surface Fourier representation We represent the plasma boundary (e.g. the LCFS) as a smooth toroidal surface embedded in $\mathbb{R}^3$, parameterized by poloidal and toroidal angles $(\theta,\phi)\in[0,2\pi)^2$. In cylindrical coordinates, the surface is given by Fourier series \[ R(\theta,\phi)=R_0+\sum_{(m,n)\in\mathcal{K}}\Bigl(R_{mn}^c\cos(m\theta-n\phi)+R_{mn}^s\sin(m\theta-n\phi)\Bigr), \] \[ Z(\theta,\phi)=\sum_{(m,n)\in\mathcal{K}}\Bigl(Z_{mn}^c\cos(m\theta-n\phi)+Z_{mn}^s\sin(m\theta-n\phi)\Bigr), \] where $R_0>0$ is a major-radius scale and $\mathcal{K}\subset\mathbb{Z}_{\ge 0}\times\mathbb{Z}$ is a finite mode set. The geometry parameter vector may be taken as \[ \mathbf{p}_{geom}:=\bigl(R_0,\{R_{mn}^c,R_{mn}^s,Z_{mn}^c,Z_{mn}^s\}_{(m,n)\in\mathcal{K}},\;a,\;\text{(additional proxy parameters)}\bigr). \] Admissibility/regularity conditions. The surface map $(\theta,\phi)\mapsto (R\cos\phi,R\sin\phi,Z)$ must be an embedding (injective, $C^2$) and must satisfy a non-self-intersection condition. In practice, the reduced model enforces this by requiring a positive minimum Jacobian determinant for the parameterization and by rejecting any $\mathbf{p}_{geom}$ for which the computed minimum coil-to-plasma distance proxy (defined below) is nonpositive. Reduced mode set convention. For compact concepts, one typically restricts $\mathcal{K}$ to low $m$ and low-to-moderate $|n|$. In the present mathematical document, $\mathcal{K}$ is treated as part of the model specification: the choice is a modeling assumption that is falsifiable via vacuum-field comparisons (Section IV.E). IV.B. Rotational transform modeling The systems model requires a flux-surface label $\rho\in[0,1]$ and a rotational transform profile $\iota(\rho)$. In a reduced setting we treat $\iota(\rho)$ as a parametric flux function rather than as a quantity derived from full equilibrium solves. A minimal parameterization is a low-order even polynomial in $\rho$: \[ \iota(\rho)=\iota_0+\iota_2\rho^2+\iota_4\rho^4, \] with coefficients included in $\mathbf{p}_{geom}$. When a monotonicity convention is imposed (e.g. to exclude low-order resonances in the core), it is implemented as an inequality constraint \[ \min_{\rho\in[0,1]}\iota'(\rho)\ge -\delta_{\iota’} \] for a declared tolerance $\delta_{\iota’}\ge 0$. Two derived scalars are used frequently as confinement proxies: \[ \bar\iota := \int_0^1 \iota(\rho)\,w(\rho)\,d\rho,\qquad \bar s := \int_0^1 s(\rho)\,w(\rho)\,d\rho, \] where $w(\rho)\ge 0$ is a fixed weight with $\int_0^1 w=1$ and \[ s(\rho):=\frac{\rho}{\iota(\rho)}\,\frac{d\iota}{d\rho}(\rho) \] is a dimensionless shear-like quantity. The model must declare $w$ (e.g. uniform in $\rho$ or volume-weighted), since it affects comparisons of proxy trends. IV.C. Geometry-quality proxies The +1D transport and constraint layers require geometric information beyond $V'(\rho)$. In a swap-ready architecture, we treat these as outputs of the geometry module: \[ \mathcal{G}:\mathbf{p}_{geom}\mapsto \mathbf{g}(\rho)=\bigl(V'(\rho),A(\rho),g_{met}(\rho),\epsilon_{eff}(\rho),\dots\bigr), \] where $g_{met}(\rho)$ denotes any metric factor needed to convert $d/d\rho$ to a physical radial derivative. IV.C.1. Effective ripple proxy $\epsilon_{eff}(\rho)$ We assume the geometry module supplies a nonnegative, dimensionless function $\epsilon_{eff}(\rho)\ge 0$ intended to summarize neoclassical/energetic-particle sensitivity to 3D magnetic spectrum. In this document $\epsilon_{eff}$ is treated purely as an input-output object: no particular computation method is assumed. To ensure that $\epsilon_{eff}$ is used consistently across transport closures, we impose a basic admissibility requirement: \[ \epsilon_{eff}\in L^{\infty}([0,1]),\qquad 0\le \epsilon_{eff}(\rho)\le \epsilon_{max} \] for a declared bound $\epsilon_{max}$. Any design point violating the bound is infeasible. IV.C.2. Magnetic-well / related configuration metrics Many reduced stability constraints require a scalar or profile-like measure of favorable average curvature or well depth. We therefore allow an additional geometry proxy $\mathcal{W}(\rho)$ (dimensionless) that is interpreted as a “well strength,” with the convention that larger $\mathcal{W}$ indicates more favorable conditions. The only structural requirement in this paper is that $\mathcal{W}(\rho)$ be reported alongside any constraint that uses it, so that the constraint is empirically testable as a function of a measurable/inferrable configuration descriptor. IV.D. Coil/engineering constraints as falsifiable inequalities Engineering feasibility is represented by explicit inequality constraints of the form \[ \mathcal{C}_{eng}(\mathbf{u},\mathbf{p})\le 0. \] In this reduced setting we treat coil-related quantities as geometry-level outputs (or measurable inputs), not as consequences of the transport solve. IV.D.1. Coil–plasma distance proxy Let $\Gamma_{plasma}$ denote the plasma boundary surface from IV.A, and let $\Gamma_{coil}$ denote a prescribed coil-winding surface (or a set of filament curves) provided in the parameter set. Define the minimum distance \[ d_{cp}:=\inf_{\mathbf{x}\in\Gamma_{plasma}}\inf_{\mathbf{y}\in\Gamma_{coil}}\|\mathbf{x}-\mathbf{y}\|. \] We enforce a required clearance by the inequality \[ \mathcal{C}_{cp}:=d_{min}-d_{cp}\le 0, \] where $d_{min}>0$ is a declared design limit. The paired falsifiability object is the reported distance $d_{cp}$ itself. IV.D.2. Peak-field proxy and current/field limits Let $B_{coil}^{peak}$ denote a proxy for the peak magnetic field on coils (or winding surface). This may be computed from a vacuum-field Biot–Savart model or treated as an externally evaluated input. The feasibility constraint is \[ \mathcal{C}_{B}:=B_{coil}^{peak}-B_{max}\le 0 \] for a declared engineering bound $B_{max}$. The output must include $B_{coil}^{peak}$ (and $B_{max}$) to make the constraint auditable. IV.D.3. Coil complexity proxy To prevent unphysical optimization toward arbitrarily high Fourier content, we include a complexity penalty or constraint based on the boundary coefficients. One purely geometric (and implementation-independent) choice is a weighted $\ell^2$ norm on the boundary modes: \[ \mathcal{K}_{coil}:=\sum_{(m,n)\in\mathcal{K}} w_{mn}\Bigl((R_{mn}^c)^2+(R_{mn}^s)^2+(Z_{mn}^c)^2+(Z_{mn}^s)^2\Bigr), \] with nonnegative weights $w_{mn}$ increasing with mode number (e.g. $w_{mn}=1+m^2+n^2$). The corresponding constraint is \[ \mathcal{C}_{K}:=\mathcal{K}_{coil}-K_{max}\le 0. \] This object is falsifiable in the limited sense that it is a declared rule of admissibility rather than a physical law; its role is to prevent the reduced representation from leaving its regime of intended smoothness. IV.E. Vacuum-field falsification hooks Because the geometry module is often the least directly validated component, we include vacuum-field comparisons that can be conducted without high-performance plasma conditions. IV.E.1. Boundary-normal field metric Let $\Gamma_{target}$ be a fixed diagnostic surface (e.g. a proxy for the desired LCFS). Let $\mathbf{n}(\theta,\phi)$ be its outward unit normal and $\mathbf{B}_{vac}(\theta,\phi)$ the predicted vacuum magnetic field from the coil model. Define the boundary-normal component \[ B_n(\theta,\phi):=\mathbf{B}_{vac}(\theta,\phi)\cdot\mathbf{n}(\theta,\phi). \] A scalar discrepancy metric is the RMS norm \[ \|B_n\|_{RMS}:=\left(\frac{1}{|\Gamma_{target}|}\int_{\Gamma_{target}} B_n^2\,dA\right)^{1/2}. \] The model must report $\|B_n\|_{RMS}$ as an observable whenever a vacuum-field coil model is used, enabling direct comparison to probe or mapping measurements. IV.E.2. Small-perturbation response tests Suppose a controllable perturbation parameter $\delta\in\mathbb{R}$ modifies the vacuum field (e.g. trim-coil current). The model predicts a linearized response \[ \|B_n\|_{RMS}(\delta)\approx \|B_n\|_{RMS}(0)+\delta\,\left.\frac{d}{d\delta}\|B_n\|_{RMS}(\delta)\right|_{\delta=0}. \] The derivative may be computed analytically if the vacuum solver is linear in current, or numerically by finite difference with a declared step size. Agreement of both (i) the baseline $\|B_n\|_{RMS}(0)$ and (ii) the slope with measured perturbation response constitutes a direct falsification hook for the geometry/coil module: systematic mismatch indicates that the reduced representation is inadequate for the intended use, independent of transport modeling. IV.E.3. Monotone trend checks across small design sets When comparing a set of nearby configurations $\{\mathbf{p}_{geom}^{(j)}\}_{j=1}^J$ for which vacuum metrics are measured, the model can be tested by monotone trend agreement between predicted and measured scalars (e.g. $\|B_n\|_{RMS}$, $B_{coil}^{peak}$). As a concrete diagnostic, one may compute a rank-correlation coefficient $\varrho\in[-1,1]$ between predicted and measured values; persistent low $|\varrho|$ indicates that the reduced parameterization fails to capture the dominant dependencies in the vacuum response. This completes the specification of the 0D geometry backbone: a declared parameter vector $\mathbf{p}_{geom}$, a well-defined interface $\mathcal{G}$ producing $\mathbf{g}(\rho)$, and explicit inequality constraints and observables that make geometry assumptions empirically testable. V. Core +1D Transport Model: Flux-Surface Averaged Two-Temperature Closure This section defines the mathematical backbone that maps geometry factors $\mathbf{g}(\rho)$, operational inputs, and transport-closure parameters to steady profiles $T_e(\rho),T_i(\rho)$ (and, optionally, $n_e(\rho)$) together with derived fluxes and global power-balance diagnostics. V.A. Steady +1D energy equations on a $\rho$-grid We write the steady, flux-surface-averaged electron and ion thermal energy balances in conservative flux form. Let $P_s(\rho)$ denote the integrated outward conductive power crossing the surface $\rho$ for species $s\in\{e,i\}$, consistent with III.C.2: \[ P_s(\rho)=A(\rho)\,q_s(\rho),\qquad [P_s]=\mathrm{W}. \] Let $p_s(\rho)$ denote the net volumetric heating power density into species $s$ (W/m$^3$), i.e. source minus sink terms, including collisional exchange between species with equal and opposite contributions: \[ \int_0^1 p_e(\rho)\,V'(\rho)\,d\rho + \int_0^1 p_i(\rho)\,V'(\rho)\,d\rho = P_{aux,abs}+P_{\alpha,dep}-P_{rad}-P_{other}, \] with the conventions of III.C. The steady balance for each species is \[ \frac{dP_s}{d\rho}(\rho)=V'(\rho)\,p_s(\rho),\qquad \rho\in(0,1). \] We impose the regularity condition at the magnetic axis (no net power crossing the axis) \[ P_s(0)=0. \] At the edge $\rho=1$, the model must declare one of the following boundary condition types for each species: (i) Dirichlet (edge temperature prescribed): $T_s(1)=T_{s,edge}$. (ii) Neumann/flux (edge conductive power prescribed): $P_s(1)=P_{s,edge}$. (iii) Robin (edge conductance): $P_s(1)=h_s\,A(1)\,(T_s(1)-T_{s,wall})$ for some effective $h_s\ge 0$. Type (iii) is convenient when a minimal edge/SOL module supplies $h_s$ and $T_{s,wall}$, but the present core model only requires that the chosen boundary condition be explicit and auditable. Discrete representation. On a grid $0=\rho_0<\rho_1<\cdots<\rho_N=1$, a finite-volume discretization is naturally expressed in terms of face fluxes $P_{s,j+1/2}$ and cell-averaged sources $p_{s,j}$. The discrete conservative identity takes the form \[ P_{s,j+1/2}-P_{s,j-1/2}=\int_{\rho_{j-1/2}}^{\rho_{j+1/2}} V'(\rho)\,p_s(\rho)\,d\rho, \] which is the per-species version of the global residual $\mathfrak{r}_{PB}$ in III.D. V.B. Base diffusive closure and falsifiable integral constraints V.B.1. Diffusive constitutive law The base closure models the radial heat flux density as purely diffusive in the physical minor-radius coordinate $r=a\rho$: \[ q_s(\rho)=-n_s(\rho)\,\chi_s(\rho)\,\frac{dT_s^{(J)}}{dr}(\rho), \] where $\chi_s$ has units m$^2$/s, $n_s$ is m$^{-3}$, and $T_s^{(J)}$ is the temperature in Joules. Writing $T_s$ in keV as in III.B, we have $T_s^{(J)}=C_{keV}T_s$ with $C_{keV}=1.602176634\times 10^{-16}\,\mathrm{J/keV}$. Therefore \[ q_s(\rho)=-C_{keV}\,n_s(\rho)\,\chi_s(\rho)\,\frac{1}{a}\,\frac{dT_s}{d\rho}(\rho). \] Equivalently, in terms of the integrated conductive power, \[ P_s(\rho)=-\frac{C_{keV}}{a}\,A(\rho)\,n_s(\rho)\,\chi_s(\rho)\,\frac{dT_s}{d\rho}(\rho). \] This relation provides a direct route to inferring an effective diffusivity from measured profiles and reconstructed power flow. V.B.2. Direct integral mapping $\{\chi(\rho),V'(\rho)\}\to\{T(\rho),P(\rho)\}$ Assume for the moment that $A(\rho),V'(\rho),n_s(\rho)$ and $\chi_s(\rho)$ are known, and that $p_s(\rho)\in L^1([0,1])$ is specified. Then the conservative balance integrates to \[ P_s(\rho)=\int_0^{\rho} V'(\xi)\,p_s(\xi)\,d\xi, \] where the integration constant is fixed by $P_s(0)=0$. Combining with the diffusive closure yields an explicit quadrature for the temperature profile: \[ \frac{dT_s}{d\rho}(\rho)=-\frac{a}{C_{keV}}\,\frac{P_s(\rho)}{A(\rho)\,n_s(\rho)\,\chi_s(\rho)}, \] so for a Dirichlet boundary condition $T_s(1)=T_{s,edge}$, \[ T_s(\rho)=T_{s,edge}+\frac{a}{C_{keV}}\int_{\rho}^{1}\frac{P_s(\xi)}{A(\xi)\,n_s(\xi)\,\chi_s(\xi)}\,d\xi. \] Proposition V.1 (well-posedness of the diffusion quadrature). Assume $A,n_s,\chi_s\in L^{\infty}([0,1])$ with $A(\rho)\ge A_{min}>0$, $n_s(\rho)\ge n_{min}>0$, and $\chi_s(\rho)\ge \chi_{min}>0$ a.e. on $[0,1]$. If $p_s\in L^1([0,1])$, then the formula above defines an absolutely continuous $T_s$ on $[0,1]$ satisfying the steady conservative balance and the Dirichlet edge condition. Proof. Since $p_s\in L^1$ and $V’\in L^{\infty}$ by construction of the geometry proxy, $P_s$ is absolutely continuous with $P_s’ = V’ p_s$ in the distributional sense and $P_s(0)=0$. The integrand $P_s/(A n_s \chi_s)$ is in $L^1$ under the stated positivity bounds, hence $T_s$ defined by the quadrature is absolutely continuous and satisfies the differential relation almost everywhere. $\square$ Falsifiable integral constraint (inference of $\chi$). Given measured $T_s(\rho)$, density $n_s(\rho)$, and reconstructed net source $p_s(\rho)$ (from power deposition and sink diagnostics), one may compute the implied integrated power $P_s(\rho)=\int_0^{\rho}V’ p_s$ and thus infer an experimental effective diffusivity \[ \chi_{s,exp}(\rho):=-\frac{a}{C_{keV}}\,\frac{P_s(\rho)}{A(\rho)\,n_s(\rho)\,T_s'(\rho)} \] where $T_s’$ denotes a suitably regularized derivative (e.g. spline fit with reported smoothing). Agreement between $\chi_{s,exp}$ and the modeled $\chi_s$ over a radius band, within propagated uncertainty, is a direct falsifiability test of the diffusion-based closure. V.C. Neoclassical+turbulent additive transport with collisionality dependence We model the effective diffusivity as an additive decomposition \[ \chi_s(\rho)=\chi_{s,nc}(\rho)+\chi_{s,turb}(\rho), \] where each term is nonnegative by construction. The purpose of this paper is to define a unit-consistent and identifiable parameterization rather than to assert a unique first-principles formula; the coefficients below are therefore treated as closure parameters $\boldsymbol{\theta}_T$ subject to calibration-or-reject rules (II.D). V.C.1. Dimensionless collisionality Let $\nu_s(\rho)$ be an effective collision frequency (s$^{-1}$) used consistently within the transport closure. Define a dimensionless collisionality proxy \[ \nu_*^{(s)}(\rho):=\frac{\nu_s(\rho)\,R_0}{v_{ts}(\rho)},\qquad v_{ts}(\rho):=\sqrt{\frac{2T_s^{(J)}(\rho)}{m_s}}, \] where $R_0$ is the major-radius scale from $\mathbf{p}_{geom}$ and $m_s$ is the species mass. V.C.2. Neoclassical surrogate structure A unit-consistent neoclassical surrogate can be expressed as \[ \chi_{s,nc}(\rho)=C_{nc,s}\,\mathcal{H}_{nc}\bigl(\nu_*^{(s)}(\rho)\bigr)\,\mathcal{E}(\rho)\,\frac{\rho_s(\rho)^2\,v_{ts}(\rho)}{R_0}, \] where \[ \rho_s(\rho):=\frac{v_{ts}(\rho)}{\Omega_s},\qquad \Omega_s:=\frac{Z_s e B_0}{m_s}, \] $\mathcal{E}(\rho)\ge 0$ is a geometry quality proxy (e.g. a function of $\epsilon_{eff}(\rho)$), and $\mathcal{H}_{nc}\ge 0$ is a dimensionless regime-interpolation function. A minimal smooth choice that interpolates between an inverse-collisionality-like and plateau-like scaling is \[ \mathcal{H}_{nc}(\nu_*):=\left(\frac{1}{(\nu_*+\nu_{*,min})^{p}}+h_0\right)^{1/p}, \] with fixed $p\ge 1$, $\nu_{*,min}>0$ (regularization), and $h_0\ge 0$. The constants $C_{nc,s},h_0$ are elements of $\boldsymbol{\theta}_T$ and must be reported. V.C.3. Turbulence surrogate with critical-gradient activation Define the normalized temperature gradient \[ \kappa_s(\rho):=\frac{a}{L_{Ts}}(\rho):=-a\,\frac{1}{T_s(\rho)}\,\frac{dT_s}{dr}(\rho)=-\frac{1}{T_s(\rho)}\,\frac{dT_s}{d\rho}(\rho). \] A gyro-Bohm-normalized turbulent diffusivity surrogate is \[ \chi_{s,turb}(\rho)=C_{gB,s}\,\frac{\rho_s(\rho)^2\,c_s(\rho)}{a}\,g\bigl(\kappa_s(\rho);\kappa_{c,s},\beta_s\bigr), \] where $c_s:=\sqrt{T_s^{(J)}/m_s}$ and \[ g(\kappa;\kappa_c,\beta):=\max\left(0,\frac{\kappa-\kappa_c}{\kappa_c}\right)^{\beta} \] with $\kappa_c>0$ and $\beta\ge 1$ as closure parameters. V.C.4. Collisionality-dependent coupling between neoclassical and turbulent terms To encode the empirical observation that regimes can change with $\nu_*$ (without asserting a specific underlying mechanism), we allow a smooth modulation of the turbulent contribution by a dimensionless factor $S(\nu_*;\theta_S)\in[0,1]$: \[ \chi_s(\rho)=\chi_{s,nc}(\rho)+S\bigl(\nu_*^{(s)}(\rho);\theta_S\bigr)\,\chi_{s,turb}(\rho), \] where a convenient two-parameter family is \[ S(\nu_*;\nu_{*0},\eta)=\frac{1}{1+(\nu_*/\nu_{*0})^{\eta}},\qquad \nu_{*0}>0,\ \eta>0. \] This preserves nonnegativity and allows falsifiable trend tests in controlled density/temperature scans. V.D. Transport\u2013stability coupling variables Many reduced turbulence activations and feasibility constraints are written in terms of a pressure-gradient drive variable and a shear variable. Let the total thermal pressure (Pa) be \[ p(\rho):=n_e(\rho)\,T_e^{(J)}(\rho)+n_i(\rho)\,T_i^{(J)}(\rho). \] Using $r=a\rho$, define a stellarator-adapted dimensionless drive proxy \[ \alpha(\rho):=-\frac{2\mu_0 R_0}{B_0^2}\,\frac{1}{\iota(\rho)^2}\,\frac{dp}{dr}(\rho). \] Since typically $dp/dr<0$ outward, $\alpha(\rho)\ge 0$ in nominal confinement regimes. We use the shear-like quantity already defined in IV.B, \[ s(\rho):=\frac{\rho}{\iota(\rho)}\,\frac{d\iota}{d\rho}(\rho). \] Stability activation/constraint. The model may implement either: (i) a hard feasibility inequality $\alpha(\rho)\le \alpha_{crit}(\rho)$ for all $\rho$, where $\alpha_{crit}$ is a declared function of $s$ and geometry proxies; or (ii) a soft activation in the turbulent closure, e.g. \[ \chi_{s,turb}(\rho)\ \mapsto\ \chi_{s,turb}(\rho)\,\sigma\bigl(\alpha(\rho)-\alpha_{crit}(\rho)\bigr), \] with $\sigma(x)=\tfrac12(1+\tanh(x/\Delta_\alpha))$ and $\Delta_\alpha>0$ declared. Both implementations are falsifiable because they predict the intermediate observables $\alpha(\rho)$, $s(\rho)$, and the implied activation radius band where $\alpha$ crosses the threshold. V.E. Power-balance closure and fast evaluation This subsection isolates a commonly used steady global closure that can be evaluated without a fixed-point iteration, provided a local power-law surrogate is assumed for the global confinement time. V.E.1. Power-law transport-throughput surrogate (assumption) Assume that in a restricted operating neighborhood the (transport) energy-throughput time depends on total heating power by a power law \[ \tau_{tr}(P_{heat})=\tau_0\left(\frac{P_{heat}}{P_0}\right)^{-\gamma},\qquad \tau_0>0,\ P_0>0,\ \gamma\ge 0, \] where $\tau_0,\gamma$ (and possibly their dependence on geometry proxies) are declared closure parameters. Here $\tau_{tr}$ is used only to close the *conductive/transport* throughput via $P_{cond}\approx W/\tau_{tr}$. The global confinement time $\tau_E$ remains defined by III.C.5 as $\tau_E:=W/P_{loss}$ with $P_{loss}=P_{cond}+P_{rad}+P_{other}$, and is computed *after* the steady balance is solved. V.E.2. Scalar steady closure as a root problem Define the total modeled loss power functional \[ \mathcal{L}(P):=\frac{W(P)}{\tau_{tr}(P)}+P_{rad}(P)+P_{other}(P), \] where $W(P),P_{rad}(P),P_{other}(P)$ are the stored energy and sinks implied by the profile solution associated with heating level $P$ (for example, through a fixed-shape parameterization or through a fast direct solve of the diffusive backbone with sources proportional to $P$). A steady global balance is equivalent to the scalar equation \[ P-\mathcal{L}(P)=0. \] Proposition V.2 (sufficient conditions for existence and uniqueness). Assume $\mathcal{L}:[0,\infty)\to[0,\infty)$ is continuous. Fix $P_{max}>0$. If there exist $0<P_-<P_+\le P_{max}$ such that \[ \mathcal{L}(P_-)>P_- \quad\text{and}\quad \mathcal{L}(P_+)<P_+, \] then the scalar steady balance equation $P-\mathcal{L}(P)=0$ has at least one root $P^*\in(P_-,P_+)$. (If $\mathcal{L}(0)=0$, then $P=0$ is a trivial root; the sign-change condition above ensures the existence of a strictly positive root.) If, in addition, $\mathcal{L}$ is nondecreasing on $[P_-,P_+]$ and satisfies the contraction-type bound \[ \sup_{P\in[P_-,P_+]} \mathcal{L}'_{+}(P) < 1, \] where $\mathcal{L}'_{+}$ denotes the upper right Dini derivative, then the root in $[P_-,P_+]$ is unique. Proof. Let $f(P):=P-\mathcal{L}(P)$. The hypotheses imply $f(P_-)<0$ and $f(P_+)>0$, so by the intermediate value theorem there exists $P^*\in(P_-,P_+)$ with $f(P^*)=0$. For uniqueness, note that the Dini bound implies that for any $P_1<P_2$ in $[P_-,P_+]$, \[ \mathcal{L}(P_2)-\mathcal{L}(P_1) \le \alpha\,(P_2-P_1)\quad\text{for some}\quad \alpha<1, \] so \[ f(P_2)-f(P_1)=(P_2-P_1)-\bigl(\mathcal{L}(P_2)-\mathcal{L}(P_1)\bigr)\ge (1-\alpha)(P_2-P_1)>0. \] Hence $f$ is strictly increasing on $[P_-,P_+]$ and can have at most one zero there. $\square$ V.E.3. No-fixed-point computation path for $P_{aux,abs}$, $\tau_E$, and $Q$ Once the scalar root $P_{heat}^*$ is obtained (by bisection or safeguarded Newton using Proposition V.2 assumptions as checks), the reported performance metrics follow deterministically: \[ \tau_{tr}^*=\tau_{tr}(P_{heat}^*),\qquad W^*=W(P_{heat}^*),\qquad P_{loss}^*=\mathcal{L}(P_{heat}^*),\qquad \tau_E^*=\frac{W^*}{P_{loss}^*}. \] The core power-balance residual reported in II.C and III.D is \[ \mathfrak{r}_{PB}=\frac{|P_{heat}^*-P_{loss}^*|}{\max(P_{heat}^*,P_{loss}^*,P_{ref})}. \] If fusion power $P_{fus}$ is evaluated from the resulting profiles, then $Q$ is computed by III.C.4. The falsifiability content of this closure is that a power scan at fixed configuration produces a predicted one-parameter family $\{(P_{heat},W,P_{rad},\tau_E)\}$ constrained by a single scalar equation; systematic violations beyond uncertainty directly refute either the assumed power-law throughput surrogate $\tau_{tr}(P_{heat})$ or the modeled sink partition. VI. Transport-Operator Closures (Eigenvalue/Sturm–Liouville Formulation) This section introduces an operator-based closure for the +1D transport backbone. The closure produces intermediate observables beyond steady profiles, namely decay rates and modulation transfer functions, which are directly comparable to perturbative experiments (power steps and periodic modulation). The construction is purely mathematical: it is a re-expression of the same diffusive transport structure as an evolution operator whose spectrum can be bounded and differentiated exactly. VI.A. Operator form of the +1D energy transport equation For clarity we present the operator closure for a single temperature-like scalar $T(\rho,t)$ on $\rho\in[0,1]$; the extension to the two-temperature case is obtained by applying the same construction to each species (with additional coupling terms treated as part of the reaction/relaxation operator). Let $r=a\rho$. Consider a linear(ized) time-dependent flux-surface averaged energy balance of the form \[ M(\rho)\,\partial_t T(\rho,t) = \frac{1}{a^2}\,\frac{d}{d\rho}\left(K(\rho)\,\frac{dT}{d\rho}(\rho,t)\right) – H(\rho)\,T(\rho,t) + S(\rho,t), \] where: – $M(\rho)>0$ is a heat-capacitance weight (units chosen so that $M\partial_t T$ has units of power density); – $K(\rho)>0$ is an effective radial conduction coefficient induced by the diffusive closure (in the simplest case, $K\propto A(\rho)n(\rho)\chi(\rho)$, up to the keV-to-J conversion and geometric factors); – $H(\rho)\ge 0$ is a linearized sink coefficient (e.g. radiative damping linearized about an equilibrium, or effective loss to unmodeled channels); – $S$ is a forcing term representing modulated heating. Boundary conditions are chosen to match the core solver’s regularity and edge closure: – Axis regularity (no singular flux): \[ \frac{dT}{d\rho}(0,t)=0. \] – Edge Robin condition (effective edge conductance): \[ -\frac{1}{a}\,K(1)\,\frac{dT}{d\rho}(1,t)= h_{edge}\,\bigl(T(1,t)-T_{wall}(t)\bigr),\qquad h_{edge}\ge 0. \] When $h_{edge}=0$ this becomes a no-flux boundary; for large $h_{edge}$ it approaches a Dirichlet boundary. Define the weighted Hilbert space $\mathsf{H}:=L^2([0,1],M\,d\rho)$ with inner product \[ \langle f,g\rangle_M := \int_0^1 M(\rho)\,f(\rho)\,g(\rho)\,d\rho. \] On the domain \[ \mathsf{D}:=\Bigl\{\varphi\in H^1([0,1]) : \varphi'(0)=0\Bigr\}, \] introduce the bilinear form \[ a(\varphi,\psi) := \int_0^1 \frac{K(\rho)}{a^2}\,\varphi'(\rho)\,\psi'(\rho)\,d\rho + \int_0^1 H(\rho)\,\varphi(\rho)\,\psi(\rho)\,d\rho + h_{edge}\,\varphi(1)\,\psi(1), \] which is symmetric and nonnegative if $K>0$, $H\ge 0$, $h_{edge}\ge 0$. The associated operator $\mathcal{L}$ is defined (in the standard weak sense) by \[ \langle \mathcal{L}\varphi,\psi\rangle_M = a(\varphi,\psi)\quad\text{for all}\ \psi\in\mathsf{D}. \] Formally, $\mathcal{L}$ acts as the Sturm–Liouville-type differential expression \[ (\mathcal{L}\varphi)(\rho)= -\frac{1}{M(\rho)}\,\frac{1}{a^2}\,\frac{d}{d\rho}\left(K(\rho)\,\frac{d\varphi}{d\rho}(\rho)\right) + \frac{H(\rho)}{M(\rho)}\,\varphi(\rho) \] with the boundary conditions above encoded by the variational term $h_{edge}\,\varphi(1)\psi(1)$. VI.B. Principal eigenvalue as a confinement/decay predictor We consider the unforced homogeneous evolution \[ \partial_t T + \mathcal{L}T = 0. \] Under the assumptions \[ K\in L^{\infty},\ \ 0<K_{min}\le K(\rho)\le K_{max},\qquad M\in L^{\infty},\ \ 0<M_{min}\le M(\rho)\le M_{max}, \] with $H\in L^{\infty}$, $H\ge 0$, and $h_{edge}\ge 0$, the operator $\mathcal{L}$ is self-adjoint and nonnegative in $\mathsf{H}$ (in the sense of forms), hence has a real spectrum bounded below by $0$. If, additionally, either $h_{edge}>0$ or $H\not\equiv 0$, then the lowest eigenvalue is strictly positive. The eigenproblem associated with $\mathcal{L}$ is the generalized Sturm–Liouville system \[ -\frac{1}{a^2}\,\frac{d}{d\rho}\left(K(\rho)\,\varphi'(\rho)\right)+H(\rho)\,\varphi(\rho)=\lambda\,M(\rho)\,\varphi(\rho),\qquad \varphi'(0)=0, \] with edge condition \[ \frac{1}{a^2}K(1)\,\varphi'(1)+h_{edge}\,\varphi(1)=0. \] Let $0\le \lambda_1\le \lambda_2\le\cdots$ denote the eigenvalues (counted with multiplicity). For solutions expanded in the $\langle\cdot,\cdot\rangle_M$-orthonormal eigenfunctions $\{\varphi_k\}$, the evolution is \[ T(\rho,t)=\sum_{k\ge 1} c_k\,e^{-\lambda_k t}\,\varphi_k(\rho), \] so the slowest decay time is $\tau_{op}:=1/\lambda_1$ (when $\lambda_1>0$). This quantity is a predicted observable of the same transport coefficients $K,M,H,h_{edge}$ used in the steady closure. VI.B.1. Variational characterization and Rayleigh bounds The principal eigenvalue admits the Rayleigh–Ritz characterization \[ \lambda_1 = \inf_{\varphi\in\mathsf{D},\ \varphi\ne 0}\ \frac{a(\varphi,\varphi)}{\langle \varphi,\varphi\rangle_M} = \inf_{\varphi\in\mathsf{D},\ \varphi\ne 0}\ \frac{\int_0^1 \frac{K}{a^2}|\varphi’|^2\,d\rho + \int_0^1 H|\varphi|^2\,d\rho + h_{edge}|\varphi(1)|^2}{\int_0^1 M|\varphi|^2\,d\rho}. \] Two immediately computable (and hence falsifiable) inequalities follow from testing the Rayleigh quotient with specific $\varphi$. Upper bound (constant test function). Let $\varphi\equiv 1$. Then $\varphi’\equiv 0$, and \[ \lambda_1 \le \frac{\int_0^1 H(\rho)\,d\rho + h_{edge}}{\int_0^1 M(\rho)\,d\rho}. \] This bound is meaningful when $H$ and $h_{edge}$ represent measurable damping/edge-loss mechanisms. Lower bound (a coarse coercivity estimate). For any $\varphi\in\mathsf{D}$, the identity $\varphi(\rho)=\varphi(1)-\int_{\rho}^1 \varphi'(s)\,ds$ and Cauchy–Schwarz imply \[ \int_0^1 |\varphi(\rho)|^2\,d\rho \le 2|\varphi(1)|^2 + 2\int_0^1 |\varphi'(\rho)|^2\,d\rho. \] Using $\int_0^1 M|\varphi|^2\le M_{max}\int_0^1|\varphi|^2$ and $a(\varphi,\varphi)\ge \frac{K_{min}}{a^2}\int|\varphi’|^2 + h_{edge}|\varphi(1)|^2$, we obtain \[ \lambda_1\ge \frac{1}{2M_{max}}\,\min\left(\frac{K_{min}}{a^2},\ h_{edge}\right). \] This bound is intentionally conservative; its role is to provide a dimensionally transparent inequality that any physically plausible implementation must satisfy. VI.B.2. Monotonicity with respect to transport and edge-loss knobs The min–max principle yields monotone dependence of $\lambda_1$ on the coefficients. Lemma VI.1 (monotonicity). Let $(K_1,H_1,M_1,h_1)$ and $(K_2,H_2,M_2,h_2)$ be two coefficient sets with $M_1\equiv M_2$ and pointwise inequalities $K_2\ge K_1$, $H_2\ge H_1$, $h_2\ge h_1$. Then the principal eigenvalues satisfy $\lambda_1^{(2)}\ge \lambda_1^{(1)}$. Proof. For any admissible $\varphi$, the numerator $a(\varphi,\varphi)$ increases under $K,H,h$ increasing, while the denominator $\langle\varphi,\varphi\rangle_M$ is unchanged. Taking the infimum over $\varphi\ne 0$ preserves the inequality. $\square$ This provides sign-definite falsification tests when a module asserts that a configuration change increases/decreases $K$ or edge conductance. VI.C. Modulation/transfer-function predictions Consider periodic forcing at angular frequency $\omega$ around a steady baseline, \[ \partial_t T + \mathcal{L}T = f(\rho)\,e^{i\omega t}, \] with $f\in\mathsf{H}$. Seeking a time-harmonic response $T(\rho,t)=\Re\{\Theta(\rho,\omega)e^{i\omega t}\}$ yields the resolvent equation \[ (i\omega + \mathcal{L})\,\Theta(\cdot,\omega)=f. \] Expanding in the eigenbasis $\{\varphi_k\}$ gives the explicit transfer-function representation \[ \Theta(\rho,\omega)=\sum_{k\ge 1}\frac{\langle f,\varphi_k\rangle_M}{i\omega+\lambda_k}\,\varphi_k(\rho). \] Therefore, the amplitude and phase of the response at any radius are predicted from $\{\lambda_k,\varphi_k\}$ and the forcing projection $\langle f,\varphi_k\rangle_M$. In particular, at low frequency $\omega\ll\lambda_2$, the response is dominated by $k=1$ with approximately uniform phase and a single-pole roll-off, while at higher frequency multiple modes contribute, producing radius-dependent phase lags. These features are directly falsifiable with heat-pulse or ECH-modulation experiments whenever $f(\rho)$ is known (or can be reconstructed) and $T$ is measured as a function of radius. VI.C.1. Eigenvalue sensitivity formula (Hellmann–Feynman for generalized eigenproblems) Let $p$ be a scalar parameter that perturbs coefficients $K,H,M,h_{edge}$ smoothly (in a differentiable sense). Write the generalized eigenproblem in weak form: \[ a_p(\varphi,\psi)=\lambda\,m_p(\varphi,\psi),\qquad m_p(\varphi,\psi):=\int_0^1 M_p\,\varphi\psi\,d\rho. \] Assume $\lambda_1(p)$ is simple and choose eigenfunctions normalized by $m_p(\varphi_1,\varphi_1)=1$. Then \[ \frac{d\lambda_1}{dp}=\frac{\partial a_p}{\partial p}(\varphi_1,\varphi_1)-\lambda_1\,\frac{\partial m_p}{\partial p}(\varphi_1,\varphi_1). \] Writing out the derivatives yields the explicit sensitivity identity \[ \frac{d\lambda_1}{dp}=\int_0^1 \frac{1}{a^2}\,\frac{\partial K}{\partial p}\,|\varphi_1’|^2\,d\rho + \int_0^1 \frac{\partial H}{\partial p}\,|\varphi_1|^2\,d\rho + \frac{\partial h_{edge}}{\partial p}\,|\varphi_1(1)|^2 – \lambda_1\int_0^1 \frac{\partial M}{\partial p}\,|\varphi_1|^2\,d\rho. \] This formula is an exact derivative given the eigenpair $(\lambda_1,\varphi_1)$ and is compatible with the swap-ready philosophy: it depends only on module-level coefficient derivatives $\partial_p K$, $\partial_p H$, $\partial_p M$, $\partial_p h_{edge}$. VI.D. Falsification tests specific to the operator closure The operator closure produces falsifiable objects $\lambda_1$ (and optionally higher $\lambda_k$ and eigenfunctions) and the modulation response $\Theta(\rho,\omega)$. The following failures constitute direct rejection of the operator-level closure for the tested conditions. 1. Eigenvalue/decay-time mismatch. If a power-step experiment produces a measured dominant decay time $\tau_{decay}^{exp}$ (e.g. from a best-fit exponential to stored-energy or core-temperature relaxation), then the model prediction $\tau_{op}=1/\lambda_1$ must satisfy \[ \left|\frac{\tau_{op}-\tau_{decay}^{exp}}{\tau_{decay}^{exp}}\right|\le \varepsilon_{\tau} \] for a declared tolerance $\varepsilon_{\tau}$ after any permitted calibration. 2. Violation of Rayleigh-quotient bounds. Given measured or reconstructed coefficient proxies entering $H$ and $h_{edge}$ (or at minimum reasonable bounds on them), the inequality \[ \lambda_1 \le \frac{\int_0^1 H\,d\rho + h_{edge}}{\int_0^1 M\,d\rho} \] must not be violated by the inferred $\lambda_1$ beyond uncertainty. Because this bound depends only on integrals and not on eigenfunction resolution, it is robust against modest numerical noise. 3. Incorrect modulation transfer function. For known forcing $f(\rho)$ and measured complex response $\Theta_{exp}(\rho,\omega)$ (amplitude and phase), the model response $\Theta_{model}$ computed from $(i\omega+\mathcal{L})^{-1}f$ must match within declared uncertainty over a radius band. Systematic phase-lag errors (e.g. correct amplitude but wrong radial phase slope across multiple $\omega$) falsify the assumed $K(\rho)$ and/or $h_{edge}$ structure even when steady profiles can be fit. 4. Wrong trend under controlled parameter changes. If a scan changes only those inputs that, under the declared closure, monotonically increase $K$ or $h_{edge}$ while leaving $M$ fixed to leading order, then Lemma VI.1 implies a monotone increase of $\lambda_1$. Observing a statistically significant decrease of the measured decay rate in such a scan constitutes a direct falsification of the asserted coefficient dependence. These tests are intentionally operator-level: they can reject a transport surrogate without requiring burning-plasma conditions, because $\lambda_1$ and $\Theta(\rho,\omega)$ are accessible on existing stellarators via small perturbations and standard profile diagnostics. VII. Physics Extension Modules Required for Q>15 Credibility (All Swap-Ready) This section specifies extensions that can be coupled to the +1D backbone while preserving the interface structure of Section II. Each extension is defined as an explicit map (possibly implicit via a residual) from the current state and parameters to either (i) modified transport coefficients and boundary data, (ii) additional state variables, or (iii) feasibility inequalities $\mathcal{C}(\mathbf{u},\mathbf{p})\le 0$. The mathematical requirement is that any introduced closure expose intermediate observables sufficient to calibrate-or-reject (II.D). VII.A. Advection–diffusion (“heat pinch”) extension VII.A.1. Flux form and steady equation A minimal extension of the diffusive heat flux allows an advective (pinch) contribution proportional to $T$: \[ q(\rho) = -n(\rho)\,\chi(\rho)\,\frac{dT^{(J)}}{dr}(\rho) + n(\rho)\,V_T(\rho)\,T^{(J)}(\rho), \] where $\chi\ge 0$ has units m$^2$/s and $V_T$ has units m/s. In keV units, $T^{(J)}=C_{keV}T$ so the same form holds with a global factor $C_{keV}$ in $q$. Let $P(\rho)=A(\rho)q(\rho)$ be the integrated conductive power across $\rho$. In steady state, \[ \frac{dP}{d\rho}(\rho)=V'(\rho)\,p(\rho),\qquad P(0)=0. \] Assuming $A,n,\chi,V_T$ are given, the constitutive relation gives a first-order linear ODE for $T$: \[ \frac{dT}{d\rho} – \\underbrace{\frac{a\,V_T(\rho)}{\chi(\rho)}}_{:=\,\eta(\rho)}\,T(\rho)= -\frac{a}{C_{keV}}\,\frac{P(\rho)}{A(\rho)\,n(\rho)\,\chi(\rho)}. \] For a Dirichlet edge condition $T(1)=T_{edge}$, the solution is given by integrating factor: \[ T(\rho)=\exp\!\left(\int_{\rho}^{1}\eta(s)\,ds\right)T_{edge} +\frac{a}{C_{keV}}\int_{\rho}^{1}\exp\!\left(\int_{\rho}^{\xi}\eta(s)\,ds\right)\frac{P(\xi)}{A(\xi)n(\xi)\chi(\xi)}\,d\xi. \] This reduces to the purely diffusive quadrature of V.B.2 when $V_T\equiv 0$. VII.A.2. Well-posedness and admissibility Assume $A,n,\chi\in L^{\infty}([0,1])$ with $A\ge A_{min}>0$, $n\ge n_{min}>0$, $\chi\ge \chi_{min}>0$, and $\eta\in L^1([0,1])$. If $p\in L^1([0,1])$ and $V’\in L^{\infty}$, then $P(\rho)=\int_0^{\rho}V'(\xi)p(\xi)d\xi$ is absolutely continuous and the formula above defines an absolutely continuous $T$ satisfying the steady balance and the edge condition. The admissibility set $\mathcal{U}_{phys}$ is extended by requiring $V_T$ and $\chi$ to keep $T\ge 0$ (or $T\ge T_{min}>0$) on $[0,1]$. VII.A.3. Identifiability and modulation falsifier Because $V_T$ and $\chi$ enter through the ratio $\eta=aV_T/\chi$ and through $\chi$ in the source term, steady profiles alone may not uniquely determine both fields. A falsifiable requirement for using the pinch extension is therefore: – either $\eta(\rho)$ is parameterized with low dimension (e.g. $\eta(\rho)=\eta_0+\eta_2\rho^2$) and calibrated with declared datasets; – or periodic modulation data are used, since the forced resolvent $(i\omega+\mathcal{L}_{adv})^{-1}$ is sensitive to advective skewness, producing measurable phase/asymmetry changes distinct from pure diffusion. A module-level reject condition is: if adding $V_T$ cannot simultaneously improve (within uncertainties) both (i) steady power-balance closure and (ii) modulation phase profiles for a fixed calibrated parameter set, then the pinch term is rejected as either unnecessary or incorrectly parameterized. VII.B. Ambipolar radial electric field $E_r$ closure VII.B.1. Implicit closure as a root condition Stellarator neoclassical transport is often strongly dependent on $E_r$. We model $E_r(\rho)$ as an additional profile determined by an ambipolarity condition written abstractly as \[ \mathcal{A}(E_r(\rho);\mathbf{x}(\rho),\theta_{Er}) = 0,\qquad \rho\in[0,1], \] where $\mathbf{x}(\rho)$ denotes local state/geometry inputs (e.g. $n_s,T_s,\epsilon_{eff},\iota$) and $\theta_{Er}$ collects any surrogate coefficients. The map $\mathcal{A}$ is required to be reported as a diagnostic residual, \[ \mathfrak{r}_{amb}(\rho):=|\mathcal{A}(E_r(\rho);\mathbf{x}(\rho),\theta_{Er})|, \] so closure adequacy is testable even when $E_r$ itself is noisy. VII.B.2. Existence/uniqueness criterion (local) At each $\rho$, suppose $\mathcal{A}(\cdot;\mathbf{x}(\rho),\theta_{Er})$ is continuous and strictly monotone in $E_r$ on an interval $[E_{min},E_{max}]$ with opposite signs at the endpoints. Then there exists a unique $E_r(\rho)\in(E_{min},E_{max})$ solving the closure. Failure of monotonicity or sign-change is treated as a module-level infeasibility flag (the closure cannot determine $E_r$ without additional physics). VII.B.3. Swap-ready interface The neoclassical diffusivity component in V.C.2 is replaced by \[ \chi_{s,nc}(\rho) = \chi_{s,nc}\bigl(\nu_*^{(s)}(\rho),\epsilon_{eff}(\rho),E_r(\rho);\theta_{nc}\bigr), \] keeping the external interface as $\chi_{s,nc}(\rho)$ and reporting $E_r(\rho)$ and $\mathfrak{r}_{amb}(\rho)$ as falsifiability objects. Higher-fidelity replacements (e.g. tabulated neoclassical calculations) are admissible provided they preserve this interface. VII.C. Energetic particle / alpha confinement and deposition VII.C.1. Timescale-competition deposition fraction To avoid an unconstrained constant $\eta_{\alpha,dep}$, we define a deposited fraction from a competition between a slowing-down (energy deposition) time $\tau_{sd}$ and a loss time $\tau_{loss}$: \[ \eta_{\alpha,dep} := \mathcal{D}\!\left(\frac{\tau_{loss}}{\tau_{sd}}\right),\qquad 0\le \mathcal{D}(x)\le 1,\ \mathcal{D}'(x)\ge 0. \] A minimal admissible choice is $\mathcal{D}(x)=1-e^{-x}$, yielding $\eta_{\alpha,dep}\approx \tau_{loss}/\tau_{sd}$ for small $\tau_{loss}$ and $\eta_{\alpha,dep}\to 1$ for $\tau_{loss}\gg \tau_{sd}$. The model must report $\tau_{sd}$, $\tau_{loss}$, and $\eta_{\alpha,dep}$. VII.C.2. Reduced dependencies and falsifiability objects – $\tau_{sd}$ is modeled as a local functional of $n_e$ and $T_e$ (and composition) and is strictly decreasing in $n_e$ and increasing in $T_e$ under the surrogate assumptions. – $\tau_{loss}$ is modeled as a functional of geometry proxies (e.g. $\epsilon_{eff}$), magnetic field strength, and possibly clearance proxies (IV.D.1) if prompt orbit loss is included. A module is rejected if a configuration scan intended to change the loss proxy (e.g. varying $\epsilon_{eff}$ with other profiles held comparable) yields measured fast-ion loss indicators with an opposite monotone trend than the predicted $\partial\eta_{\alpha,dep}/\partial\epsilon_{eff}$ sign. VII.D. Particle balance, fueling, and composition dynamics VII.D.1. 0D inventory closure for species Let $N_s$ be the total particle inventory of species $s$ in the confined volume: \[ N_s := \int_0^1 n_s(\rho)\,V'(\rho)\,d\rho. \] A minimal steady particle balance is \[ 0 = S_s – \frac{N_s}{\tau_{p,s}} – L_s, \] where $S_s$ is a net source rate (particles/s), $\tau_{p,s}>0$ is a particle confinement time, and $L_s\ge 0$ collects any explicitly modeled sinks (e.g. pumping, charge-exchange out of the confined volume if modeled separately). Solving gives \[ N_s = \tau_{p,s}(S_s-L_s),\qquad \text{requiring } S_s\ge L_s. \] The implied fueling requirement is a directly falsifiable observable: $S_s$ needed to maintain a measured $N_s$ given an inferred $\tau_{p,s}$. VII.D.2. Helium ash fraction (steady) Let $S_{He}$ denote the helium production rate from fusion (particles/s), computed from $P_{fus}$ and the reaction energy. In steady state with negligible additional He sources, \[ N_{He} = \tau_{p,He}\,S_{He}. \] Define a global ash fraction by $f_{He}:=N_{He}/N_e$ (or another declared normalization). Then \[ f_{He} = \frac{\tau_{p,He}S_{He}}{N_e}. \] Because $S_{He}$ is predicted from the energy model, this yields a nonlinear coupling: increased $P_{fus}$ increases $f_{He}$, which can reduce fuel density and thus reduce $P_{fus}$. The module must report $f_{He}$ and a feasibility inequality such as \[ \mathcal{C}_{He}:=f_{He}-f_{He,max}\le 0, \] where $f_{He,max}$ is a declared limit. VII.D.3. Impurity dilution (minimal) Let $Z_{eff}$ be either measured input or a closure output. A minimal dilution factor for fusion-capable ions can be represented as a scalar $f_{dil}\in(0,1]$ multiplying the available fuel density in the fusion source: \[ n_D n_T \mapsto f_{dil}^2\,(n_D n_T). \] The module must report $f_{dil}$ as an observable and declare how it depends on ash/impurity inventories. VII.E. Radiation and $Z_{eff}$ closure (including impurity dilution feedback) We represent the radiated power density as \[ p_{rad}(\rho) = p_{brem}(\rho) + p_{line}(\rho),\qquad p_{rad}\ge 0. \] A unit-consistent bremsstrahlung surrogate has the generic form \[ p_{brem}(\rho) = C_{brem}\,Z_{eff}(\rho)\,n_e(\rho)^2\,\sqrt{T_e^{(J)}(\rho)}, \] with $C_{brem}>0$ a known dimensional constant (or an interface-provided constant if using a particular unit system). Line radiation is represented as \[ p_{line}(\rho) = n_e(\rho)^2\,\mathcal{L}\bigl(T_e(\rho);\theta_{line}\bigr), \] where $\mathcal{L}\ge 0$ is a cooling-curve surrogate with a low-dimensional parameter set $\theta_{line}$ (or a swap-ready lookup table). The falsifiability contract requires that the module expose: – the predicted profile $p_{rad}(\rho)$ and integral $P_{rad}$; – the decomposition fraction $P_{line}/P_{rad}$ (or analogous) when line radiation is active; – the calibrated constants (if any) and their allowed ranges. A reject condition is: if a controlled impurity-content change (with measured $Z_{eff}$ and comparable heating) cannot be matched in both $P_{rad}$ and the resulting profile response with a single fixed $\theta_{line}$, then the radiation closure is rejected. VII.F. Stellarator-specific density limit surrogate (radiative/Sudo-type) VII.F.1. Inequality form We impose a density limit as a feasibility constraint involving a scalar density measure (e.g. volume-average) \[ \bar n_e := \frac{1}{V(1)}\int_0^1 n_e(\rho)\,V'(\rho)\,d\rho \] and a predicted maximum $n_{max}$: \[ \mathcal{C}_{dens}:=\bar n_e – n_{max}(\mathbf{p},P_{heat};\theta_{dens})\le 0. \] The function $n_{max}$ is treated as a swap-ready surrogate. The module must report $n_{max}$ and the margin $n_{max}-\bar n_e$ as explicit observables. VII.F.2. Nonlinear coupling and multiple branches If $\bar n_e$ is itself solved (rather than prescribed) via particle balance, then $n_{max}$ and the energy balance couple nonlinearly through $P_{\alpha,dep}(n,T)$ and $P_{rad}(n,T,Z_{eff})$, potentially producing multiple steady solutions. Mathematically, if the reduced steady problem can be written as a scalar equation in $\bar n_e$, \[ F(\bar n_e)=0, \] then multiple branches correspond to multiple roots. A minimal falsifiability object in this case is the predicted number of admissible roots in an experimentally scanned parameter range (e.g. a density ramp), together with the predicted critical point where roots appear/disappear (fold). The model is rejected if observed behavior (single-valued vs hysteretic/multi-valued) disagrees robustly with the predicted branch structure under the same declared closure. VII.G. Bootstrap current $\to$ $\iota$ shift $\to$ resonance/islands $\to$ transport penalty loop VII.G.1. Algebraic bootstrap closure and $\Delta\iota$ We model the bootstrap current by an algebraic functional of profiles and geometry, \[ I_{bs} = \mathcal{B}(\mathbf{u},\mathbf{g};\theta_{bs}), \] and represent its effect on the rotational transform proxy as an additive shift \[ \iota(\rho)\ \mapsto\ \iota(\rho)+\Delta\iota(\rho),\qquad \Delta\iota(\rho)=\mathcal{S}_{\iota}(\rho)\,I_{bs}, \] where $\mathcal{S}_{\iota}(\rho)$ is an interface-provided sensitivity profile (or a low-dimensional parameterization). Both $I_{bs}$ and $\Delta\iota(\rho)$ are required observables. VII.G.2. Resonance indicator and feasibility Let $\mathcal{C}_{max}(\mathbf{u},\mathbf{p})\ge 0$ denote a scalar indicator of island overlap/stochasticity risk, computed from $\iota+\Delta\iota$ and any declared resonance data. Feasibility can be enforced either as a hard inequality \[ \mathcal{C}_{max}\le 1, \] or as a soft confinement penalty applied to the power-balance closure, \[ \tau_E \mapsto \tau_E\,\exp\!\left(-C_{stoch}\,\max(0,\mathcal{C}_{max}-1)^p\right),\qquad C_{stoch}>0,\ p\ge 1. \] This penalty is falsifiable: if a scan changes $\mathcal{C}_{max}$ substantially according to the closure while other dominant knobs are held comparable, the measured confinement trend must correlate with the predicted penalty factor. VII.H. Pressure-gradient (ballooning/interchange) stability constraint ($s$–$\alpha$ surrogate) We reuse the drive proxy $\alpha(\rho)$ and shear $s(\rho)$ defined in V.D. A minimal stability module specifies a critical curve $\alpha_{crit}(\rho)$ as an explicit functional \[ \alpha_{crit}(\rho)=\mathcal{K}_{\alpha}\bigl(s(\rho),\mathbf{g}(\rho);\theta_{stab}\bigr), \] with $\mathcal{K}_{\alpha}>0$. The feasibility constraint is \[ \mathcal{C}_{stab}(\rho):=\alpha(\rho)-\alpha_{crit}(\rho)\le 0\quad\text{for a.e. }\rho\in[0,1], \] or, in discretized form, at all grid points $\rho_j$. The module must report: – $\alpha(\rho_j)$, $s(\rho_j)$, $\alpha_{crit}(\rho_j)$; – the maximum margin $\max_j (\alpha(\rho_j)-\alpha_{crit}(\rho_j))$. A reject condition is: if instability-threshold scans show a systematic mismatch in the inferred threshold ordering (e.g. configurations predicted to have larger $\alpha_{crit}$ consistently disrupt/lose confinement earlier), beyond uncertainty, then $\mathcal{K}_{\alpha}$ (or its assumed dependencies) is rejected. VII.I. Minimal scrape-off-layer (SOL) / edge boundary module The core solver in V.A permits an effective edge Robin condition expressed via an edge conductance $h_{edge}$. A minimal SOL module supplies $h_{edge}$ and (optionally) $T_{wall}$ from a reduced parallel-conduction balance. VII.I.1. Edge closure as a power-throughput law Define the power leaving the confined region through the edge for a given channel as \[ P_{edge} := h_{edge}\,A(1)\,(T(1)-T_{wall}),\qquad h_{edge}\ge 0. \] This equation is not a new physical assertion; it is an interface: any edge/SOL model that can be reduced to a linearized relation between boundary temperature and exhaust power supplies $h_{edge}$. VII.I.2. A concrete reduced $h_{edge}$ from parallel conduction (one possible surrogate) If the SOL parallel heat transport is conduction-dominated with $q_{\parallel}\propto T^{7/2}/L_{\parallel}$ along a connection length $L_{\parallel}$, then linearizing around an operating point $T_{edge}$ yields an effective boundary conductance scaling \[ h_{edge} \approx \frac{\partial P_{edge}}{\partial T(1)}\bigg|_{T_{edge}}\frac{1}{A(1)} \ \propto\ \frac{T_{edge}^{5/2}}{L_{\parallel}}\,\mathcal{G}_{geo}, \] where $\mathcal{G}_{geo}$ is a declared geometry factor converting parallel exhaust to an effective cross-field boundary loss. The module must report $L_{\parallel}$, $\mathcal{G}_{geo}$, and $h_{edge}$. VII.I.3. Falsifiability objects The edge module is tested by comparing predicted and measured edge quantities under controlled scans: – $T_{edge}$ (from Thomson/ECE where available) versus the implied value from the core solution and $h_{edge}$; – wall/divertor heat flux proxies versus $P_{edge}$; – trend with geometry modifications entering $\mathcal{G}_{geo}$ and $L_{\parallel}$. Systematic mismatch of these intermediate observables constitutes rejection of the edge closure even if core profiles can be matched by retuning $\chi$, because $h_{edge}$ is precisely the interface that separates edge physics from core transport. VIII. Time-Dependent 0D Burning-Plasma Thermal Stability Module This section defines a minimal time-dependent module that wraps the steady closures of Sections III–VII into a lumped (0D) dynamical system for thermal stability. The goal is not to resolve profile evolution, but to (i) predict whether the steady operating point is dynamically attracting, (ii) detect parameter regions with multiple steady branches (hysteresis), and (iii) define a small-signal thermal eigenvalue $\gamma_T$ that is directly falsifiable by power-step experiments. VIII.A. Dynamic energy-balance ODE around steady solutions We use a scalar thermal state variable $W(t)\ge 0$ (J), interpreted as the total thermal stored energy as in III.C.1, but now time-dependent. The reduced energy balance is \[ \frac{dW}{dt} = P_{heat}(W;\mathbf{p}) – P_{loss}(W;\mathbf{p}), \tag{VIII.1} \] where $\mathbf{p}$ denotes fixed inputs/parameters (geometry, boundary conditions, densities/composition when treated as prescribed, and closure coefficients). 1. Heating functional. We write \[ P_{heat}(W;\mathbf{p}) := P_{aux,abs}(W;\mathbf{p}) + P_{\alpha,dep}(W;\mathbf{p}), \tag{VIII.2} \] where $P_{aux,abs}\ge 0$ may be prescribed as an external time-dependent control input or closed by a controller (VIII.D), and $P_{\alpha,dep}\ge 0$ is the deposited alpha power implied by the fusion model and the energetic-particle module (VII.C). 2. Loss functional. The reduced total loss power is \[ P_{loss}(W;\mathbf{p}) := P_{cond}(W;\mathbf{p}) + P_{rad}(W;\mathbf{p}) + P_{other}(W;\mathbf{p}), \tag{VIII.3} \] with each term defined consistently with the steady bookkeeping in III.C.2–III.C.5, but evaluated on the profile state implied by $W$ under the chosen profile-parameterization or fast backbone solve. 3. Consistency with steady solutions. Any steady state $\bar W$ of (VIII.1) satisfies \[ P_{heat}(\bar W;\mathbf{p})=P_{loss}(\bar W;\mathbf{p}), \tag{VIII.4} \] which is the global steady balance. Thus the dynamic module introduces no new steady equilibria beyond those already admitted by the steady closure; it only classifies their dynamical stability. Assumption VIII.A (regularity for linearization). For the local stability results below, assume $P_{heat}(\cdot;\mathbf{p})$ and $P_{loss}(\cdot;\mathbf{p})$ are continuously differentiable in a neighborhood of $\bar W$. This assumption is not automatic if the steady closure includes hard max/min switches (e.g. activation via $\max(0,\cdot)$); in that case the module must use either (i) smooth approximations (e.g. softplus/tanh) or (ii) one-sided derivatives and piecewise-linear stability conditions. VIII.B. Multiple equilibria (branching) as a falsifiable outcome Define the net power function \[ F(W;\mathbf{p}) := P_{heat}(W;\mathbf{p})-P_{loss}(W;\mathbf{p}). \tag{VIII.5} \] Equilibria are roots of $F$. Multiple branches occur when $F(\cdot;\mathbf{p})$ has more than one positive root. 1. Fold (saddle-node) condition. A necessary condition for the creation/annihilation of equilibria under a smooth one-parameter continuation $\mathbf{p}=\mathbf{p}(\mu)$ is the fold condition \[ F(\bar W;\mathbf{p})=0,\qquad \partial_W F(\bar W;\mathbf{p})=0. \tag{VIII.6} \] When $\partial_W^2 F(\bar W;\mathbf{p})\neq 0$ and $\partial_\mu F(\bar W;\mathbf{p})\neq 0$, standard implicit-function arguments imply that the equilibrium curve $W(\mu)$ turns at $\mu$ (a fold). The module must report, for any scanned control parameter (e.g. $P_{aux}$ or $\bar n_e$ if solved), whether it predicts: (i) a single equilibrium for all $\mu$ in the scan range, or (ii) the existence of a fold point characterized by (VIII.6). 2. Falsifiability object. The predicted number of admissible steady equilibria in a scan is itself an observable: if the module predicts two stable branches separated by an unstable branch (see VIII.C), then a quasi-static ramp in the control parameter should exhibit hysteresis in $W$ (or in any monotone proxy such as central temperature). Conversely, if the module predicts a unique equilibrium for all $\mu$ in a scan range, observation of repeatable hysteresis in otherwise comparable conditions is a direct reject condition for the reduced $F$ functional (typically indicating missing state variables or missing delays). VIII.C. Small-signal thermal eigenvalue $\gamma_T$: sign and magnitude Let $\bar W$ be a steady solution: $F(\bar W;\mathbf{p})=0$. Linearizing (VIII.1) about $\bar W$ gives \[ \frac{d}{dt}\,\delta W(t) = \gamma_T\,\delta W(t), \qquad \gamma_T := \partial_W F(\bar W;\mathbf{p}). \tag{VIII.7} \] Therefore: – $\gamma_T<0$ implies local asymptotic stability (exponential return with time constant $\approx 1/|\gamma_T|$). - $\gamma_T>0$ implies linear instability (thermal runaway within the reduced model). – $\gamma_T=0$ corresponds to neutral linear stability and is the fold/threshold condition (VIII.6). The module must report $\gamma_T$ as a primary intermediate observable. Magnitude calibration check. If a power-step experiment around an operating point yields a measured dominant relaxation time $\tau_T^{exp}$ in stored energy (or in a monotone proxy), and if the perturbation is sufficiently small that the linearization is valid, then the model prediction must satisfy \[ \left|\frac{1/|\gamma_T| – \tau_T^{exp}}{\tau_T^{exp}}\right|\le \varepsilon_T \tag{VIII.8} \] for a declared tolerance $\varepsilon_T$, after any allowed calibration. Systematic mismatch of the sign of $\gamma_T$ (stable predicted but unstable observed, or vice versa) is an immediate reject condition for the coupled $P_{heat}$/$P_{loss}$ dependence on $W$. Remark (relation to $\tau_E$). The reduced eigenvalue $\gamma_T$ is not the same object as the operator decay rate $\lambda_1$ in VI.B. $\lambda_1$ describes spatial diffusion/edge-loss decay of perturbations under a linear transport operator, while $\gamma_T$ is the net-feedback slope of global heating minus global losses (including temperature-dependent fusion and radiation). The two can be jointly tested: $\lambda_1$ by modulation/propagation experiments (Section VI) and $\gamma_T$ by global thermal relaxation under controlled heating steps. VIII.D. $\alpha$-heating emulator protocol (non-burning devices) Because true alpha heating is not available on present non-DT stellarators, we define an experimental emulation protocol in which auxiliary heating is closed-loop controlled to reproduce a chosen temperature dependence. 1. Emulator definition. Fix a measured scalar temperature proxy $\Theta(t)$ (e.g. a volume-average or a central $T_e$ proxy) and a reference value $\Theta_{ref}>0$. Define an emulator heating law \[ P_{aux,abs}(t)=P_0 + k\,\Bigl(\frac{\Theta(t)}{\Theta_{ref}}\Bigr)^{\nu}, \qquad P_0\ge 0,\ k\ge 0,\ \nu>0. \tag{VIII.9} \] This is a purely operational control law. The corresponding reduced dynamics is still (VIII.1), but now $P_{heat}$ inherits an explicit positive feedback component. The emulator parameters $(k,\nu)$ are declared inputs and therefore constitute a falsifiable experiment design. 2. Predicted stability shift. Suppose $\Theta$ is a differentiable function of $W$ locally, with $d\Theta/dW>0$. Then the emulator modifies the thermal eigenvalue by adding \[ \Delta\gamma_T = \partial_W P_{aux,abs}(\bar W) = k\,\nu\,\Theta_{ref}^{-\nu}\,\Theta(\bar W)^{\nu-1}\,\frac{d\Theta}{dW}(\bar W)\ \ge 0. \tag{VIII.10} \] Thus increasing $k$ (or $\nu$) monotonically increases $\gamma_T$ at fixed $\bar W$, potentially driving a stable operating point toward instability. This monotone effect is itself falsifiable: if the controller gain $k$ is increased in repeated experiments at comparable conditions, the measured relaxation rate must change with the sign predicted by (VIII.10), within uncertainty in $d\Theta/dW$. 3. Emulator-based falsifier. Choose $(k,\nu)$ so that the reduced model predicts $\gamma_T$ crossing zero at a reachable operating point. If experiments do not exhibit any qualitative change in stability behavior (e.g. onset of slow divergence or bistability/hysteresis) when the predicted $\gamma_T$ changes sign by a substantial margin, then either the reduced mapping $W\mapsto\Theta$ or the modeled $P_{loss}(W)$ dependence is inadequate. Conversely, observation of instability where the model predicts $\gamma_T\ll 0$ rejects the loss/heating slope closure. The emulator protocol makes the stability module testable without requiring DT plasmas: it converts the mathematically specified thermal feedback in $P_{\alpha,dep}(W)$ into a controllable auxiliary feedback that can be dialed and compared to the predicted $\gamma_T$ and branch structure. IX. Numerics for the +1D Backbone: Conservative, Positivity-Preserving Discretization This section specifies a conservative finite-volume discretization for the steady +1D transport backbone in V.A, together with mandatory discrete diagnostics. The purpose is to ensure that any implementation of the residual $\mathcal{R}_T$ (II.B.2) has (i) a discrete global conservation identity consistent with III.C–III.D, (ii) a systematic criterion for rejecting nonphysical discrete solutions (negativity and/or oscillations), and (iii) auditable numerical diffusion when advection (pinch) terms are present. IX.A. Finite-volume discretization of the two-temperature equations We discretize each species balance in conservative form on $\rho\in[0,1]$. Let \[ 0=\rho_0<\rho_1<\cdots<\rho_N=1 \] be cell centers, with faces $\rho_{j+1/2}$ satisfying $\rho_{j-1/2}<\rho_j<\rho_{j+1/2}$ and $\rho_{-1/2}=0$, $\rho_{N+1/2}=1$. Define physical radial spacing $\Delta r_{j+1/2}:=a(\rho_{j+1}-\rho_j)$ and cell volume weights \[ \Delta V_j := \int_{\rho_{j-1/2}}^{\rho_{j+1/2}} V'(\rho)\,d\rho>0. \] For each species $s\in\{e,i\}$, denote the unknown temperatures by $T_{s,j}\approx T_s(\rho_j)$. Face-conductive powers $P_{s,j+1/2}\approx P_s(\rho_{j+1/2})$ are defined as the primary conservative flux variables. Steady finite-volume balance. The steady equation $dP_s/d\rho = V’p_s$ is discretized as \[ P_{s,j+1/2}-P_{s,j-1/2} = \Delta V_j\,p_{s,j},\qquad j=0,1,\dots,N, \tag{IX.1} \] where $p_{s,j}$ is a cell-average approximation to $p_s$ on $[\rho_{j-1/2},\rho_{j+1/2}]$. Axis condition. The regularity condition $P_s(0)=0$ becomes \[ P_{s,-1/2}=0. \tag{IX.2} \] Edge boundary conditions. At the edge face $\rho_{N+1/2}=1$, we allow the same types as in V.A: 1. Dirichlet: $T_{s,N}=T_{s,edge}$ (implemented by replacing the equation at $j=N$ with the constraint). 2. Prescribed edge power: $P_{s,N+1/2}=P_{s,edge}$. 3. Robin: $P_{s,N+1/2}=h_s A(1)(T_{s,N}-T_{s,wall})$, with $h_s\ge 0$. Diffusive constitutive relation. For the diffusion-only closure (V.B.1), the integrated conductive power at a face is approximated by a two-point flux of the form \[ P_{s,j+1/2} = -\mathcal{K}_{s,j+1/2}\,(T_{s,j+1}-T_{s,j}), \tag{IX.3} \] with a nonnegative face transmissibility \[ \mathcal{K}_{s,j+1/2} := \frac{C_{keV}}{a^2}\,A_{j+1/2}\,\bigl(n_s\chi_s\bigr)_{j+1/2}\,\frac{a}{\rho_{j+1}-\rho_j} = \frac{C_{keV}}{a}\,A_{j+1/2}\,\bigl(n_s\chi_s\bigr)_{j+1/2}\,\frac{1}{\rho_{j+1}-\rho_j}, \tag{IX.4} \] where $A_{j+1/2}:=A(\rho_{j+1/2})$ and $(n_s\chi_s)_{j+1/2}$ is a declared averaging rule (e.g. harmonic average to preserve flux continuity when coefficients vary strongly). The factor $C_{keV}$ is required for unit consistency (III.C.1). With (IX.3), the steady discretization becomes a tridiagonal linear system for each species whenever $p_{s,j}$ does not depend on $T$; with temperature-dependent sources/sinks and $\chi_s(T,\nabla T)$, (IX.1)–(IX.4) define a nonlinear residual to be solved by Newton or a safeguarded fixed-point method. Two-temperature coupling. Electron-ion energy exchange terms should be discretized in a cell-local and energy-conserving way: if $p_{ei}(\rho)$ denotes the volumetric exchange from electrons to ions (so that $p_e$ includes $-p_{ei}$ and $p_i$ includes $+p_{ei}$), then the discretization must satisfy for each cell \[ \Delta V_j\,p_{e,ei,j} + \Delta V_j\,p_{i,ei,j} = 0. \tag{IX.5} \] Violation of (IX.5) constitutes an implementation error because it breaks exact internal energy conservation even if global power balance is otherwise enforced. IX.B. Discrete global conservation identity and the global power-balance residual Summing (IX.1) over all cells $j=0,\dots,N$ and using (IX.2) yields the exact discrete identity \[ P_{s,N+1/2} = \sum_{j=0}^{N} \Delta V_j\,p_{s,j}. \tag{IX.6} \] Summing over species gives \[ \sum_s P_{s,N+1/2} = \sum_s \sum_{j=0}^N \Delta V_j\,p_{s,j}. \tag{IX.7} \] The right-hand side is the discretized net input to the confined region (auxiliary absorption plus alpha deposition minus radiation and other sinks) partitioned by species. Thus (IX.7) is the discrete statement that the net volumetric power equals the net outward boundary throughput. Mandatory diagnostic residual. Every computed steady solution must report the dimensionless global power-balance residual \[ \mathfrak{r}_{PB}^{(disc)} := \frac{\left|\,P_{heat}^{(disc)}-P_{loss}^{(disc)}\,\right|}{\max\bigl(P_{heat}^{(disc)},P_{loss}^{(disc)},P_{ref}\bigr)}, \tag{IX.8} \] where \[ P_{loss}^{(disc)}:=\sum_s P_{s,N+1/2} + P_{rad}^{(disc)} + P_{other}^{(disc)} \tag{IX.9} \] with $P_{rad}^{(disc)}$ and $P_{other}^{(disc)}$ computed from the same discrete profiles. A reported $Q$ is admissible only if $\mathfrak{r}_{PB}^{(disc)}\le \varepsilon_{PB}$ (III.D). Module-level falsifier. Any case where a solver reports convergence of the profile residuals while (IX.6)–(IX.7) fail beyond tolerance indicates a broken conservative implementation (e.g. inconsistent discretizations of flux and volumetric terms). Such cases must be treated as numerical rejection, not as physics discrepancy. IX.C. Positivity / maximum-principle conditions and reject conditions The admissible state set $\mathcal{U}_{phys}$ (II.A) requires $T_{s,j}\ge 0$ (or $T_{s,j}\ge T_{min}>0$). The discretization must therefore (i) be constructed to preserve positivity whenever the continuous problem is positivity-preserving under the declared sources/sinks, and (ii) reject any computed state violating positivity. IX.C.1. Diffusion-only discrete maximum principle (linear case) Consider one species with diffusion-only flux (IX.3) and a linear sink term $-\sigma_j T_j$ with $\sigma_j\ge 0$, so that \[ \Delta V_j\,p_j = \Delta V_j\,(s_j-\sigma_j T_j) \tag{IX.10} \] with a nonnegative source $s_j\ge 0$. With Dirichlet or Robin boundary conditions that are nonnegative (e.g. $T_{edge}\ge 0$, $T_{wall}\ge 0$), the resulting linear system can be written \[ \mathbf{A}\mathbf{T}=\mathbf{b}. \tag{IX.11} \] A sufficient condition for a discrete maximum principle is that $\mathbf{A}$ be an M-matrix: diagonal entries positive, off-diagonal entries nonpositive, and row sums nonnegative, with strict diagonal dominance on at least one row (or irreducibility). For the tridiagonal diffusion operator formed from $\mathcal{K}_{j\pm 1/2}\ge 0$, these conditions hold if $\sigma_j\ge 0$ and Robin/Dirichlet terms enter with nonnegative coefficients. In this case $\mathbf{A}^{-1}\ge 0$ componentwise, hence $\mathbf{b}\ge 0\Rightarrow \mathbf{T}\ge 0$. Implementation requirement. When a module claims the diffusion-only closure with linear sinks in a particular operating mode, it must verify (and if needed, report) that the assembled matrix $\mathbf{A}$ is an M-matrix. If not, the discretization is not guaranteed positivity-preserving and the run must be flagged. IX.C.2. Nonlinear sources/sinks: positivity enforcement and rejection When $p_j$ depends nonlinearly on $T$ (fusion heating, radiation, temperature-dependent $\chi$, etc.), positivity is not automatic. The solver must therefore enforce one of the following admissible strategies: 1. Positivity-preserving iterate: use a Newton or fixed-point method with step control (damping/line search) that enforces $T_{s,j}^{(k+1)}\ge T_{min}$ at every iterate. 2. Variable transform: solve for $\theta_{s,j}=\log(T_{s,j})$ or $T_{s,j}=T_{min}+\exp(\theta_{s,j})$. Reject conditions (mandatory): after convergence, if any $T_{s,j}<0$ or if the solver required clipping of $T$ without accounting for it in the residual evaluation, then the computed point is infeasible and must not be used for $Q$, $\tau_E$, $\lambda_1$, or any falsifiability object. IX.D. Advection stabilization criterion and explicit reporting of numerical diffusion For the advection--diffusion extension (VII.A), the steady constitutive relation implies a flux that depends on both $T$ and $\nabla T$. A naive central differencing of the advective part can produce spurious oscillations and loss of positivity when advection dominates diffusion. IX.D.1. Cell-face Péclet number Define a face Péclet number for species/field $T$: \[ \mathrm{Pe}_{j+1/2} := \frac{|V_{T,j+1/2}|\,\Delta r_{j+1/2}}{\chi_{j+1/2}}, \tag{IX.12} \] where $\Delta r_{j+1/2}=a(\rho_{j+1}-\rho_j)$. Large $\mathrm{Pe}$ indicates advection-dominated transport. IX.D.2. Required stabilization rule An implementation must declare an advection discretization scheme and report $\max_{j}\mathrm{Pe}_{j+1/2}$. Two acceptable classes are: 1. Upwind (first-order) conservative flux: \[ P_{j+1/2} = -\mathcal{K}_{j+1/2}(T_{j+1}-T_j) + A_{j+1/2}C_{keV}\,n_{j+1/2}\,V_{T,j+1/2}\,T_{upwind}, \tag{IX.13} \] where $T_{upwind}=T_j$ if $V_{T,j+1/2}>0$ and $T_{upwind}=T_{j+1}$ if $V_{T,j+1/2}<0$. This is positivity-preserving under standard conditions but introduces numerical diffusion. 2. Exponential fitting (Scharfetter--Gummel-type) flux for the linear ODE structure in VII.A.1: interpret the face flux as the exact solution of a constant-coefficient advection--diffusion problem on the interval $[\rho_j,\rho_{j+1}]$. This yields a conservative flux that remains stable at large $\mathrm{Pe}$ and reduces numerical diffusion relative to (IX.13). Any such scheme must be documented by giving the resulting explicit face-flux formula. Numerical diffusion reporting (mandatory). If upwinding or any limiter is used, the module must report an equivalent added diffusivity $\chi_{num,j+1/2}\ge 0$ such that the total flux can be written as \[ P_{j+1/2}= -\frac{C_{keV}}{a}\,A_{j+1/2}(n\chi_{eff})_{j+1/2}\,\frac{T_{j+1}-T_j}{\rho_{j+1}-\rho_j} + \text{(remaining consistent terms)}, \tag{IX.14} \] with $\chi_{eff}=\chi+\chi_{num}$. The quantity $\chi_{num}$ must be included among reported diagnostics because it affects falsifiability comparisons of inferred $\chi_{exp}$ (V.B.2). IX.E. Verification checks tied to falsifiability The discretization is part of the falsifiability contract (II.C, II.D): predicted profiles/fluxes must be numerically trustworthy and self-consistent. Each run must therefore produce the following verification outputs. 1. Discrete conservation closure. Report (i) $\mathfrak{r}_{PB}^{(disc)}$ from (IX.8), and (ii) the per-species residuals \[ \mathfrak{r}_{s}^{(disc)}:=\frac{\left|\,P_{s,N+1/2}-\sum_{j}\Delta V_j p_{s,j}\,\right|}{\max\bigl(|P_{s,N+1/2}|,\sum_j \Delta V_j|p_{s,j}|,P_{ref}\bigr)}. \tag{IX.15} \] Large $\mathfrak{r}_s^{(disc)}$ invalidates any comparison of $q_s(\rho)$, $\chi_s(\rho)$, or operator-based objects $\lambda_1$ derived from the same discrete state. 2. Positivity/monotonicity flags. Report \[ T_{s,min}:=\min_j T_{s,j},\qquad n_{min}:=\min_j n_j \tag{IX.16} \] and a boolean flag indicating whether $T_{s,min}\ge T_{min}$ (and similarly for density if solved). If the model imposes monotonic temperature profiles (optional admissibility rule), report the discrete monotonicity violation magnitude $\max_j (T_{j+1}-T_j)_+$. 3. Grid convergence indicator. Because falsifiability comparisons may be sensitive to numerical resolution, an implementation must support at least two grid resolutions $N$ and $2N$ and report a self-consistency metric for a scalar of interest $\Psi$ (e.g. $W$, $P_{cond}$, $\lambda_1$): \[ \mathfrak{r}_{grid}(\Psi):=\frac{|\Psi_{2N}-\Psi_N|}{\max(|\Psi_{2N}|,|\Psi_N|,\Psi_{ref})}. \tag{IX.17} \] If $\mathfrak{r}_{grid}$ exceeds a declared tolerance for any reported falsifiability object, the run must be flagged as under-resolved. 4. Operator-consistency check (when VI is used). If the operator closure is computed from the same $K,M,H,h_{edge}$ that underlie the discrete transport model, the implementation must verify that the discrete operator used for eigenvalues/transfer functions is consistent with the finite-volume flux discretization (e.g. the same transmissibilities $\mathcal{K}_{j+1/2}$ and boundary terms). Otherwise, disagreement between steady and modulation predictions could be purely numerical. These numerical requirements ensure that disagreements with data can be attributed to modeling assumptions (transport/sink closures) rather than to hidden discretization artifacts. X. Exact Sensitivity Analysis and Design Optimization Layer This section derives exact (non-finite-difference) sensitivities of performance metrics and constraint functions with respect to parameters $\mathbf{p}$. The derivations are stated at the discretized level $\mathbf{u}\in\mathbb{R}^d$, because the coupled map is implemented numerically and must expose auditable derivatives consistent with the discrete conservation identities of Section IX. X.A. Variational / implicit-function formulation Let the steady coupled system be the residual equation from II.A, \[ \mathcal{R}(\mathbf{u},\mathbf{p})=\mathbf{0},\qquad \mathcal{R}:\mathbb{R}^d\times\mathbb{R}^m\to\mathbb{R}^d. \] Let $\Psi(\mathbf{u},\mathbf{p})\in\mathbb{R}$ be any scalar quantity of interest, such as $Q$, $W$, $\tau_E$, $P_{rad}$, $\lambda_1$, or a soft-penalty objective combining multiple terms. Assumption X.A (regularity and solvability). Fix a point $(\bar\mathbf{u},\bar\mathbf{p})$ with $\mathcal{R}(\bar\mathbf{u},\bar\mathbf{p})=0$. Assume: 1. $\mathcal{R}$ is continuously differentiable in a neighborhood of $(\bar\mathbf{u},\bar\mathbf{p})$. 2. The Jacobian $\mathcal{J}:=\partial_{\mathbf{u}}\mathcal{R}(\bar\mathbf{u},\bar\mathbf{p})\in\mathbb{R}^{d\times d}$ is nonsingular. Then, by the finite-dimensional implicit function theorem, there exists a neighborhood $\mathcal{N}$ of $\bar\mathbf{p}$ and a unique differentiable solution map $\mathbf{u}(\mathbf{p})$ on $\mathcal{N}$ satisfying $\mathcal{R}(\mathbf{u}(\mathbf{p}),\mathbf{p})=0$ and $\mathbf{u}(\bar\mathbf{p})=\bar\mathbf{u}$. Differentiating the identity $\mathcal{R}(\mathbf{u}(\mathbf{p}),\mathbf{p})\equiv 0$ yields \[ \mathcal{J}\,\frac{d\mathbf{u}}{d\mathbf{p}} = -\partial_{\mathbf{p}}\mathcal{R}(\bar\mathbf{u},\bar\mathbf{p}). \tag{X.1} \] Consequently, the total derivative of $\Psi(\mathbf{u}(\mathbf{p}),\mathbf{p})$ at $\bar\mathbf{p}$ is \[ \frac{d\Psi}{d\mathbf{p}} = \partial_{\mathbf{u}}\Psi\,\frac{d\mathbf{u}}{d\mathbf{p}} + \partial_{\mathbf{p}}\Psi. \tag{X.2} \] Remark (nonsmooth closures). If $\mathcal{R}$ contains kinks (e.g. $\max(0,\cdot)$ activations), then $\partial_{\mathbf{u}}\mathcal{R}$ may fail to exist at switching points. In that case one must either (i) replace switches by smooth approximations and treat the smoothing scale as a declared numerical parameter, or (ii) use piecewise-differentiation with one-sided derivatives and explicitly report when a derivative crosses a switching surface. Either approach is admissible only if the implementation reports the active-set state (which switches are active) as part of the diagnostic output. X.B. Adjoint method for exact gradients with many parameters Direct evaluation of (X.1) for all components of $\mathbf{p}$ requires solving $m$ linear systems with matrix $\mathcal{J}$, which is expensive when $m$ is large. The adjoint method reduces the cost to (typically) one linear solve per scalar objective $\Psi$. Let $\Psi$ be a scalar objective. Define the adjoint variable $\boldsymbol{\lambda}\in\mathbb{R}^d$ as the solution of the adjoint system \[ \mathcal{J}^T\,\boldsymbol{\lambda} = \bigl(\partial_{\mathbf{u}}\Psi(\bar\mathbf{u},\bar\mathbf{p})\bigr)^T. \tag{X.3} \] Then the gradient of the reduced objective $\psi(\mathbf{p}):=\Psi(\mathbf{u}(\mathbf{p}),\mathbf{p})$ at $\bar\mathbf{p}$ is \[ \frac{d\psi}{d\mathbf{p}} = \partial_{\mathbf{p}}\Psi(\bar\mathbf{u},\bar\mathbf{p}) - \boldsymbol{\lambda}^T\,\partial_{\mathbf{p}}\mathcal{R}(\bar\mathbf{u},\bar\mathbf{p}). \tag{X.4} \] Proof. From (X.2) and (X.1), \[ \frac{d\psi}{d\mathbf{p}} = \partial_{\mathbf{u}}\Psi\,\Bigl(-\mathcal{J}^{-1}\partial_{\mathbf{p}}\mathcal{R}\Bigr) + \partial_{\mathbf{p}}\Psi. \] Let $\boldsymbol{\lambda}$ satisfy (X.3), i.e. $\boldsymbol{\lambda}^T=\partial_{\mathbf{u}}\Psi\,\mathcal{J}^{-1}$. Substitution gives (X.4). $\square$ X.B.1. Boundary conditions and block structure In the discretizations of Section IX, $\mathbf{u}$ aggregates unknowns across radii and species, and boundary conditions are either embedded as rows of $\mathcal{R}$ (Dirichlet constraints) or as flux/balance equations including boundary terms (Neumann/Robin). The adjoint system (X.3) uses the transpose of the same assembled Jacobian, hence its boundary behavior is automatically consistent with the primal discretization. This is essential for falsifiability: an adjoint derived from a different operator than the primal can give formally correct sensitivities for a different numerical scheme. X.B.2. Separation of module derivatives (swap-ready gradient contract) The expression (X.4) depends on $\partial_{\mathbf{p}}\mathcal{R}$ and $\partial_{\mathbf{p}}\Psi$. To preserve the swap-ready architecture of II.B, we require that each module expose its own partial derivatives with respect to its declared inputs. Concretely, suppose the residual is assembled from module contributions \[ \mathcal{R}(\mathbf{u},\mathbf{p}) = \mathcal{R}^{(core)}(\mathbf{u};\mathbf{g},\mathbf{p}_{op},\boldsymbol{\theta}) + \sum_{r} \mathcal{R}^{(r)}(\mathbf{u};\cdots), \] where $\mathbf{g}=\mathcal{G}(\mathbf{p}_{geom})$ is the geometry output and $\mathcal{R}^{(r)}$ are extension residuals (e.g. additional balance equations). Then \[ \partial_{\mathbf{p}}\mathcal{R} = \partial_{\mathbf{g}}\mathcal{R}\,\partial_{\mathbf{p}}\mathbf{g} + \partial_{\mathbf{p}_{op}}\mathcal{R} + \partial_{\boldsymbol{\theta}}\mathcal{R}, \tag{X.5} \] where each term is supplied by the corresponding module. A replacement geometry module changes $\partial_{\mathbf{p}}\mathbf{g}$ but does not change the adjoint wrapper; a replacement transport module changes $\partial_{\mathbf{g}}\mathcal{R}$ and $\partial_{\boldsymbol{\theta}}\mathcal{R}$, but again leaves (X.4) unchanged. X.C. Sensitivity formulas tied to operator closures When the operator closure of Section VI is used, additional objectives and constraints can depend on operator-level quantities, such as the principal eigenvalue $\lambda_1$ or modulation response $\Theta(\rho,\omega)$. This subsection records exact derivative identities in a form compatible with the adjoint framework. X.C.1. Eigenvalue derivative (discrete and continuous) Let $\lambda_1$ be computed from a discretized generalized eigenproblem \[ \mathbf{A}(\mathbf{p})\,\boldsymbol{\varphi} = \lambda\,\mathbf{M}(\mathbf{p})\,\boldsymbol{\varphi}, \tag{X.6} \] where $\mathbf{A}$ and $\mathbf{M}$ are symmetric and $\mathbf{M}$ is positive definite (the matrix analog of the weighted form in VI.A). Assume $\lambda_1$ is simple and normalize $\boldsymbol{\varphi}_1^T\mathbf{M}\boldsymbol{\varphi}_1=1$. Then for any scalar parameter $p$, the Hellmann--Feynman identity yields \[ \frac{d\lambda_1}{dp} = \boldsymbol{\varphi}_1^T\Bigl(\frac{\partial \mathbf{A}}{\partial p} - \lambda_1\frac{\partial \mathbf{M}}{\partial p}\Bigr)\boldsymbol{\varphi}_1. \tag{X.7} \] The continuous coefficient-level formula in VI.C.1 is the direct analog of (X.7) at the variational level. In a swap-ready implementation, the operator module must provide $\partial_p \mathbf{A}$ and $\partial_p \mathbf{M}$ by differentiating the same coefficient assembly used to compute $\mathbf{A}$ and $\mathbf{M}$. X.C.2. Rayleigh-quotient sensitivity bounds Even when full derivatives $\partial_p\mathbf{A}$ are expensive, one can bound eigenvalue changes using Rayleigh quotients. For any admissible test vector $\mathbf{v}\ne 0$, define \[ \mathcal{R}_Q(\mathbf{v};\mathbf{p}) := \frac{\mathbf{v}^T\mathbf{A}(\mathbf{p})\mathbf{v}}{\mathbf{v}^T\mathbf{M}(\mathbf{p})\mathbf{v}}. \] Then $\lambda_1(\mathbf{p})\le \mathcal{R}_Q(\mathbf{v};\mathbf{p})$. If $\mathbf{v}$ is held fixed while $\mathbf{p}$ varies, differentiating $\mathcal{R}_Q$ gives an explicit bound on the directional derivative of $\lambda_1$ whenever monotonicity conditions (as in Lemma VI.1) apply. These inequalities are useful as falsifiable sign checks: if the code claims $\partial_p\mathbf{A}\succeq 0$ in a scan (in the sense of quadratic forms) and $\mathbf{M}$ fixed, then $\lambda_1$ must be nondecreasing with $p$. X.C.3. Modulation/transfer-function sensitivity (resolvent derivative) Let the complex modulation response be defined by the linear system \[ \bigl(i\omega\mathbf{M}+\mathbf{A}\bigr)\,\boldsymbol{\Theta}(\omega)=\mathbf{f}, \tag{X.8} \] where $\mathbf{f}$ represents a discretized forcing profile in the $\langle\cdot,\cdot\rangle_M$ inner product. For a parameter $p$, differentiating (X.8) gives \[ \bigl(i\omega\mathbf{M}+\mathbf{A}\bigr)\,\frac{d\boldsymbol{\Theta}}{dp} = -\Bigl(i\omega\frac{\partial\mathbf{M}}{\partial p}+\frac{\partial\mathbf{A}}{\partial p}\Bigr)\boldsymbol{\Theta} + \frac{\partial\mathbf{f}}{\partial p}. \tag{X.9} \] Thus sensitivities of $\boldsymbol{\Theta}$ reduce to solving one additional linear system with the same complex matrix as the forward modulation solve, plus module-supplied coefficient derivatives. When the objective depends on amplitude/phase at selected radii, (X.9) can be combined with an adjoint of the complex system to compute gradients with respect to many parameters efficiently. X.D. Analytic sensitivities for simplified 0D power-law confinement surrogates In the special case where the global closure is reduced to a scalar root as in V.E.2, one can obtain closed-form derivative expressions that serve as both computational accelerators and sign-definite falsifiability predictions. Let the steady scalar equation be \[ F(P;\mathbf{p}) := P-\mathcal{L}(P;\mathbf{p})=0, \tag{X.10} \] with $\mathcal{L}(P;\mathbf{p})$ continuous and differentiable in $P$ near a physical root $P^*$. Assume $\partial_P F(P^*;\mathbf{p})\ne 0$, i.e. \[ 1-\partial_P\mathcal{L}(P^*;\mathbf{p})\ne 0. \tag{X.11} \] Then by the implicit function theorem, \[ \frac{dP^*}{dp_i} = -\frac{\partial_{p_i}F(P^*;\mathbf{p})}{\partial_P F(P^*;\mathbf{p})} = \frac{\partial_{p_i}\mathcal{L}(P^*;\mathbf{p})}{1-\partial_P\mathcal{L}(P^*;\mathbf{p})}. \tag{X.12} \] If $\tau_E(P)=\tau_0(P/P_0)^{-\gamma}$ and $\mathcal{L}(P)=W(P)/\tau_E(P)+P_{rad}(P)+P_{other}(P)$, then $\partial_P\mathcal{L}$ contains the explicit contribution from $\tau_E$: \[ \partial_P\Bigl(\frac{W(P)}{\tau_E(P)}\Bigr)=\frac{W'(P)}{\tau_E(P)}+\frac{\gamma}{P}\,\frac{W(P)}{\tau_E(P)}. \tag{X.13} \] These expressions make explicit which assumptions control sensitivity amplification: as $\partial_P\mathcal{L}\to 1$, the denominator in (X.12) approaches zero and small parameter changes produce large changes in $P^*$, which is precisely the onset of fold-like behavior in scalar closures. Log-derivative ranking. For positive quantities $X(\mathbf{p})$, define the dimensionless sensitivity \[ S_{X,p_i}:=\frac{\partial \ln X}{\partial \ln p_i}=\frac{p_i}{X}\,\frac{\partial X}{\partial p_i}. \tag{X.14} \] When the objective is $Q=P_{fus}/P_{aux,abs}$ and $P_{aux,abs}$ is treated as an independent design/control variable, the ranking of geometric or closure parameters can be expressed through $S_{P_{fus},p_i}$, which in turn depends on sensitivities of the temperature profiles and thus of $P^*$ via (X.12) in reduced closures or via the full adjoint (X.4) in +1D closures. Cross-sensitivities. Second derivatives (Hessian-vector products) are obtainable by differentiating (X.3) and (X.4) with respect to $\mathbf{p}$. Because full Hessians are expensive in large $m$, an admissible optimization interface is the Hessian-vector product $\nabla^2\psi(\mathbf{p})\mathbf{v}$, which can be computed by a second-order adjoint or by differentiating the linearized system in direction $\mathbf{v}$. Any claimed parameter synergy (non-additivity) corresponds to off-diagonal Hessian structure and is testable by comparing predicted mixed partial signs $\partial^2\psi/\partial p_i\partial p_j$ against finite perturbation experiments. X.E. Constraint handling as augmented objectives Feasibility constraints are expressed as $\mathcal{C}(\mathbf{u},\mathbf{p})\le 0$ (II.B.3). For gradient-based design, two standard differentiable surrogates are admissible. X.E.1. Quadratic penalty Define a penalized objective \[ \Psi_{pen}(\mathbf{u},\mathbf{p}) := \Psi(\mathbf{u},\mathbf{p}) + \frac{\mu}{2}\,\|\bigl(\mathcal{C}(\mathbf{u},\mathbf{p})\bigr)_+\|_2^2, \tag{X.15} \] where $(\cdot)_+$ is the componentwise positive part and $\mu>0$ is a declared penalty weight. When $\mathcal{C}$ is differentiable, $\partial_{\mathbf{u}}\Psi_{pen}$ and $\partial_{\mathbf{p}}\Psi_{pen}$ are available and can be inserted into the adjoint system (X.3) to obtain gradients (X.4). If $\mathcal{C}$ is nondifferentiable at the boundary, one may replace $(x)_+$ by a smooth approximation and report the smoothing parameter. X.E.2. Log-barrier for strict interior optimization When maintaining strict feasibility is required, define the barrier objective on the interior $\{\mathcal{C}<0\}$: \[ \Psi_{bar}(\mathbf{u},\mathbf{p}) := \Psi(\mathbf{u},\mathbf{p}) - \mu\sum_{i=1}^{\ell} \ln\bigl(-\mathcal{C}_i(\mathbf{u},\mathbf{p})\bigr),\qquad \mu>0. \tag{X.16} \] The derivatives are \[ \partial_{\mathbf{u}}\Psi_{bar} = \partial_{\mathbf{u}}\Psi – \mu\sum_i \frac{1}{\mathcal{C}_i}\,\partial_{\mathbf{u}}\mathcal{C}_i, \qquad \partial_{\mathbf{p}}\Psi_{bar} = \partial_{\mathbf{p}}\Psi – \mu\sum_i \frac{1}{\mathcal{C}_i}\,\partial_{\mathbf{p}}\mathcal{C}_i. \tag{X.17} \] These expressions are well-defined only when $\mathcal{C}_i<0$. Therefore any implementation using a barrier must report the minimum constraint margin $\min_i (-\mathcal{C}_i)$ as a diagnostic and reject steps that attempt to cross feasibility boundaries. X.E.3. Adjoint gradients with constraints For either (X.15) or (X.16), the adjoint system remains (X.3) with $\Psi$ replaced by the chosen augmented objective. The gradient formula (X.4) then provides exact derivatives of the augmented objective with respect to all parameters using one adjoint solve. This is the primary mechanism by which the systems model supports optimization while preserving the calibrate-or-reject contract: the optimization layer cannot silently alter physics, but can only move within the declared parameter set $\mathbf{p}$ using derivatives computed from the declared residual equations and constraints. XI. Statistical Calibration, Model Discrimination, and Optimal Experimental Design (OED) This section formalizes the statistical wrapper that turns the deterministic map $\mathcal{F}:\mathbf{p}\mapsto\widehat{\mathbf{y}}$ into a testable (calibrate-or-reject) model against multi-diagnostic data. The emphasis is on (i) explicit measurement operators mapping state to diagnostic channels, (ii) constrained calibration in a low-dimensional parameter subset, (iii) rejection tests with uncertainty accounting, and (iv) experiment design criteria that rank which measurements best constrain the transport and constraint modules. XI.A. Calibrate-or-reject wrapper: limited fits followed by predictive validation Fix a partition of the full parameter vector $\mathbf{p}$ into measured/declared inputs $\mathbf{p}_{meas}$ and calibratable coefficients $\boldsymbol{\theta}_{cal}$. For each experimental condition $k$ (shot, time-window, configuration), let $\mathbf{p}_{meas}^{(k)}$ be the recorded inputs and let $\mathbf{y}^{(k)}\in\mathbb{R}^{k_y}$ be a vector of measured observables (profiles sampled at radii, integrated powers, modulation transfer-function values, etc.). The model predicts \[ \widehat{\mathbf{y}}^{(k)}(\boldsymbol{\theta}_{cal}) := \mathcal{F}\bigl(\mathbf{p}_{meas}^{(k)},\boldsymbol{\theta}_{cal}\bigr), \] where all non-calibrated coefficients are held fixed at declared values. We adopt a Gaussian error model with covariance $\Sigma^{(k)}\succ 0$ (estimated from diagnostics and reconstruction uncertainties), \[ \mathbf{y}^{(k)} = \widehat{\mathbf{y}}^{(k)}(\boldsymbol{\theta}_{cal}) + \boldsymbol{\varepsilon}^{(k)}, \qquad \boldsymbol{\varepsilon}^{(k)}\sim\mathcal{N}(0,\Sigma^{(k)}). \] The negative log-likelihood (up to additive constants) is the weighted least squares \[ J(\boldsymbol{\theta}_{cal})=\frac12\sum_{k=1}^{N}\bigl(\mathbf{y}^{(k)}-\widehat{\mathbf{y}}^{(k)}(\boldsymbol{\theta}_{cal})\bigr)^T\bigl(\Sigma^{(k)}\bigr)^{-1}\bigl(\mathbf{y}^{(k)}-\widehat{\mathbf{y}}^{(k)}(\boldsymbol{\theta}_{cal})\bigr). \tag{XI.1} \] A calibration produces an estimator $\boldsymbol{\theta}_{cal}^*$ (e.g. constrained minimizer of $J$ over a declared box $\Theta_{cal}$). Validation is performed on a disjoint dataset $\mathcal{D}_{val}$ with the calibrated coefficients frozen. XI.B. Chi-square rejection tests and uncertainty accounting XI.B.1. Chi-square statistic with multi-diagnostic aggregation For a validation dataset $\mathcal{D}_{val}\subset\{1,\dots,N\}$, define \[ \chi^2_{val}:=\sum_{k\in\mathcal{D}_{val}} \bigl(\mathbf{y}^{(k)}-\widehat{\mathbf{y}}^{(k)}(\boldsymbol{\theta}_{cal}^*)\bigr)^T\bigl(\Sigma^{(k)}\bigr)^{-1}\bigl(\mathbf{y}^{(k)}-\widehat{\mathbf{y}}^{(k)}(\boldsymbol{\theta}_{cal}^*)\bigr). \tag{XI.2} \] Under the stated Gaussian model and assuming the prediction map is correct and $\Sigma^{(k)}$ is accurate, $\chi^2_{val}$ concentrates near its degrees of freedom. In practice, the model is rejected if $\chi^2_{val}$ exceeds a pre-declared threshold or if residuals show structured bias in any declared falsifiability object (e.g. consistent radius-dependent sign in $T_e$ misfit or systematic phase-lag discrepancy in transfer functions). XI.B.2. Nuisance parameters and correlated errors Diagnostic reconstructions often contain correlated errors and nuisance degrees of freedom. We represent nuisance parameters by $\boldsymbol{\eta}\in\mathbb{R}^{m_\eta}$ and extend the prediction to \[ \widehat{\mathbf{y}}^{(k)}(\boldsymbol{\theta}_{cal},\boldsymbol{\eta})=\mathcal{H}^{(k)}\bigl(\mathbf{u}^{(k)}(\boldsymbol{\theta}_{cal}),\mathbf{p}_{meas}^{(k)},\boldsymbol{\eta}\bigr), \tag{XI.3} \] where $\mathcal{H}^{(k)}$ is a measurement operator (defined in XI.B.3) and $\mathbf{u}^{(k)}$ is the model state solving $\mathcal{R}(\mathbf{u},\mathbf{p})=0$. Two admissible treatments are: 1. Profiling (generalized least squares): minimize $J(\boldsymbol{\theta}_{cal},\boldsymbol{\eta})$ jointly over $(\boldsymbol{\theta}_{cal},\boldsymbol{\eta})$ with declared bounds/priors on $\boldsymbol{\eta}$. 2. Marginalization (Bayesian): introduce a prior $\pi(\boldsymbol{\eta})$ and use the marginal likelihood \[ \mathcal{L}_{marg}(\boldsymbol{\theta}_{cal}) \propto \int \exp\bigl(-J(\boldsymbol{\theta}_{cal},\boldsymbol{\eta})\bigr)\,\pi(\boldsymbol{\eta})\,d\boldsymbol{\eta}. \tag{XI.4} \] The falsifiability contract requires that any nuisance parameters be explicitly declared (and few in number); otherwise, the calibration can absorb genuine model error. XI.B.3. Measurement operators (synthetic diagnostics) A diagnostic does not measure $\mathbf{u}$ directly; it measures a transformed quantity. For each condition $k$, define an operator \[ \mathcal{H}^{(k)}:\mathcal{U}\times\mathcal{P}\times\mathbb{R}^{m_\eta}\to\mathbb{R}^{k_y}, \qquad \widehat{\mathbf{y}}^{(k)}=\mathcal{H}^{(k)}(\mathbf{u}^{(k)},\mathbf{p}_{meas}^{(k)},\boldsymbol{\eta}). \tag{XI.5} \] Examples of admissible $\mathcal{H}^{(k)}$ include: - Channel sampling: $T_e(\rho)\mapsto (T_e(\rho_{ch,1}),\dots,T_e(\rho_{ch,M}))$. - Line-integrated signals: $p_{rad}(\rho)\mapsto$ predicted bolometer chords via known geometry weights. - Diamagnetic energy: $\mathbf{u}\mapsto W(\mathbf{u})$ as in III.C.1 with explicitly recorded conversion factors. - Modulation response: $\mathcal{L}\mapsto\Theta(\rho,\omega)$ sampled in radius and frequency (VI.C), including a consistent mapping from modulated absorbed power to forcing $f$. Because $\mathcal{H}^{(k)}$ enters the statistical objective, it is part of the model and must itself be auditable; systematic mismatch attributable to an incorrect $\mathcal{H}^{(k)}$ is a rejection of the diagnostic model, not a license to retune physics closures. XI.C. Fisher information and identifiability ranking across measurements Let $\boldsymbol{\theta}\in\mathbb{R}^{m_\theta}$ denote a set of parameters to be estimated (typically $\boldsymbol{\theta}_{cal}$, possibly augmented by nuisance parameters). Define the stacked residual vector \[ \mathbf{r}(\boldsymbol{\theta}) := \begin{bmatrix} \Sigma^{(1)-1/2}(\mathbf{y}^{(1)}-\widehat{\mathbf{y}}^{(1)}(\boldsymbol{\theta}))\\ \vdots\\ \Sigma^{(N)-1/2}(\mathbf{y}^{(N)}-\widehat{\mathbf{y}}^{(N)}(\boldsymbol{\theta})) \end{bmatrix}. \tag{XI.6} \] At a reference $\bar{\boldsymbol{\theta}}$, the Gauss--Newton approximation to the Hessian of $2J$ is \[ \mathbf{F}(\bar{\boldsymbol{\theta}}) := \mathbf{J}_r(\bar{\boldsymbol{\theta}})^T\mathbf{J}_r(\bar{\boldsymbol{\theta}}), \qquad \mathbf{J}_r:=\partial_{\boldsymbol{\theta}}\mathbf{r}. \tag{XI.7} \] Under the Gaussian model, $\mathbf{F}$ is the Fisher information matrix. It provides: - Identifiability: small eigenvalues of $\mathbf{F}$ indicate poorly constrained parameter combinations. - Parameter ranking: for scalar $\theta_i$, the marginal variance lower bound is approximately $(\mathbf{F}^{-1})_{ii}$ (Cram\'er--Rao), provided the linearization is accurate. - Diagnostic value: adding a measurement channel updates $\mathbf{F}$ by an additive positive semidefinite term, enabling quantitative comparison of candidate diagnostics. Implementation note. The sensitivities $\partial_{\boldsymbol{\theta}}\widehat{\mathbf{y}}$ should be computed by the exact derivative mechanisms of Section X (adjoint or implicit differentiation) to avoid false identifiability conclusions due to finite-difference noise. XI.D. OED-guided minimal experiment sets An experiment design chooses a set of conditions (actuators, scan points, modulation frequencies, diagnostic channels) to maximize information about $\boldsymbol{\theta}$ subject to cost/feasibility constraints. Let a candidate experiment $e$ contribute a Fisher increment $\mathbf{F}_e(\boldsymbol{\theta})\succeq 0$. For a set $\mathcal{E}$ of experiments, \[ \mathbf{F}_{\mathcal{E}} = \sum_{e\in\mathcal{E}} \mathbf{F}_e. \tag{XI.8} \] Common scalar criteria include: - D-optimality: maximize $\log\det(\mathbf{F}_{\mathcal{E}})$ (minimizes volume of the linearized confidence ellipsoid). - A-optimality: minimize $\mathrm{tr}(\mathbf{F}_{\mathcal{E}}^{-1})$ (minimizes average variance). - E-optimality: maximize the smallest eigenvalue of $\mathbf{F}_{\mathcal{E}}$ (improves worst-case identifiability). A minimal, falsifiability-driven design objective is: choose $\mathcal{E}$ so that $\mathbf{F}_{\mathcal{E}}$ is well-conditioned on the declared calibratable subspace, while ensuring that the experiments also excite the operator-level observables ($\lambda_1$, modulation transfer functions) needed to discriminate diffusion-only from advection--diffusion closures. XI.E. Dimensionless-similarity reduction and similarity-law falsifiers Reduced transport surrogates often assert that responses depend primarily on a small set of dimensionless control parameters. This subsection defines a rigorous similarity check that can reject such assertions without requiring that a specific scaling exponent be correct. XI.E.1. Dimensionless grouping and invariance statement Let $\mathbf{z}\in\mathbb{R}^{m_z}$ denote a vector of dimensionless groups constructed from $\mathbf{p}_{meas}$ and selected outputs (e.g. normalized gyroradius proxy, collisionality proxy, normalized pressure proxy, geometry proxies). We do not assume a specific $\mathbf{z}$ here; we only require that the model declare (i) the mapping $\mathbf{p}\mapsto\mathbf{z}(\mathbf{p})$ and (ii) which observables are expected to be approximately invariant at fixed $\mathbf{z}$ (e.g. shape of normalized profiles, normalized transfer functions). A similarity hypothesis can be stated as: there exists a function $\Phi$ such that for the chosen observable vector $\mathbf{y}$, \[ \widehat{\mathbf{y}} \approx \Phi(\mathbf{z}) \tag{XI.9} \]\nup to declared residual terms capturing known symmetry-breaking effects (e.g. boundary conditions or specific constraint activations). XI.E.2. Similarity-law violation test Given two conditions $k$ and $\ell$ with $\|\mathbf{z}^{(k)}-\mathbf{z}^{(\ell)}\|$ small, define the normalized discrepancy \[ D_{k\ell} := \bigl(\mathbf{y}^{(k)}-\mathbf{y}^{(\ell)}\bigr)^T\bigl(\Sigma^{(k)}+\Sigma^{(\ell)}\bigr)^{-1}\bigl(\mathbf{y}^{(k)}-\mathbf{y}^{(\ell)}\bigr). \tag{XI.10} \] If the similarity hypothesis (XI.9) is correct for the selected observables and the chosen $\mathbf{z}$, then for a collection of near-similar pairs $(k,\ell)$ the statistics $D_{k\ell}$ should be consistent with noise. Persistent large $D_{k\ell}$ indicates that either (i) the chosen dimensionless group set $\mathbf{z}$ is incomplete, or (ii) the model is missing physics that breaks similarity in the tested regime. Either outcome is a falsifier for the declared similarity reduction. XI.E.3. Connection to model discrimination Similarity checks provide a diagnostic for whether discrepancies should be attributed to parameter mis-calibration versus structural model error. If calibration improves fit at one $\mathbf{z}$ but similarity violations persist across conditions with comparable $\mathbf{z}$, then the functional form of $\mathcal{F}$ (transport closure or constraint activation) is structurally inconsistent with the similarity hypothesis and must be revised rather than re-tuned. This completes the statistical wrapper: the coupled map is made empirically testable by explicit measurement operators, auditable uncertainty models, exact-derivative-based identifiability analysis, and experiment design criteria that prioritize measurements most capable of discriminating between competing transport-operator closures. XII. Integrated Falsification Plan (12–18 Month Horizon) This section consolidates the model-level and module-level falsification objects introduced previously into a concrete near-term test program. The goal is not to “validate” the model by fitting many coefficients, but to execute a limited calibration (II.D, XI.A) followed by predictive scans that can reject incorrect closures using multi-diagnostic comparisons. Throughout, let $\boldsymbol{\theta}_{cal}$ denote the declared low-dimensional calibration set, and let $\boldsymbol{\theta}_{cal}^*$ be the calibrated value obtained on a designated calibration dataset. All tests below are performed on disjoint validation data with $\boldsymbol{\theta}_{cal}$ frozen. Disagreements are quantified using (XI.2) and by structured residual tests on individual falsifiability objects. XII.A. Core transport falsifiers XII.A.1. Power scans at fixed configuration (global closure test) Protocol. Fix a configuration (geometry proxies $V',A,\epsilon_{eff},\iota$ and engineering constraints) and perform a scan in absorbed auxiliary power $P_{aux,abs}$ over a range where the steady assumption is defensible (small $|dW/dt|$ compared to $P$). For each scan point $k$, assemble the validation observable vector including at minimum \[ \mathbf{y}^{(k)} = \bigl(W^{exp},\ P_{rad}^{exp},\ P_{aux,abs}^{exp},\ T_e^{exp}(\rho_{ch}),\ T_i^{exp}(\rho_{ch})\bigr) \] with a declared covariance $\Sigma^{(k)}$ (XI.A). Predicted consistency constraints. 1. Discrete power balance: for each point, the solver must satisfy $\mathfrak{r}_{PB}^{(disc)}\le \varepsilon_{PB}$ (IX.B). Points violating this are numerical rejects and are not admitted as physical predictions. 2. If the scalar-root closure (V.E.2) is used in the given operating neighborhood, then the measured scan should lie on the predicted one-parameter family induced by $P-\mathcal{L}(P)=0$. Practically, this is tested by evaluating the chi-square statistic (XI.2) on the scan data with fixed $\boldsymbol{\theta}_{cal}^*$. Systematic departures (e.g. a curvature mismatch in $W$ vs $P_{aux,abs}$ not explainable by uncertainties) reject either the assumed $\tau_{tr}(P)$ structure or the modeled sink partition. Reject condition (global, example form). Reject the global closure on the scanned configuration if $\chi^2_{val}$ exceeds a pre-declared threshold, or if scan residuals show consistent sign-bias in at least one of $W$, $P_{rad}$, or $T_e(\rho)$ across a nontrivial fraction of scan points. XII.A.2. Profile-based transport inference (local closure test) Protocol. For steady windows with measured profiles $T_s^{exp}(\rho)$, density $n^{exp}(\rho)$, and reconstructed net sources $p_s^{exp}(\rho)$, compute the implied integrated power \[ P_s^{exp}(\rho)=\int_0^{\rho} V'(\xi)\,p_s^{exp}(\xi)\,d\xi \] (using the same geometry normalization as the model), and infer the experimental effective diffusivity (V.B.2) \[ \chi_{s,exp}(\rho):=-\frac{a}{C_{keV}}\,\frac{P_s^{exp}(\rho)}{A(\rho)\,n_s^{exp}(\rho)\,\partial_{\rho}T_s^{exp}(\rho)} \] with a declared smoothing/regularization procedure for $\partial_{\rho}T_s^{exp}$ and propagated uncertainty. Comparison. Compare $\chi_{s,exp}(\rho)$ against the modeled $\chi_{s,nc}(\rho)+\chi_{s,turb}(\rho)$ (V.C) using a banded misfit metric, e.g. \[ D_s := \int_{\rho_{min}}^{\rho_{max}} \frac{\bigl(\chi_{s,exp}(\rho)-\chi_{s,model}(\rho)\bigr)^2}{\sigma_{\chi,s}(\rho)^2}\,d\rho, \] where $\sigma_{\chi,s}(\rho)$ is a declared uncertainty envelope derived from diagnostic uncertainties and source reconstruction error. Reject condition. Reject the transport-coefficient closure on the tested regime if $D_s$ persistently exceeds a declared threshold across multiple conditions and the discrepancy has stable structure (e.g. the model systematically underpredicts mid-radius $\chi$ while matching edge values), indicating a structural closure error rather than calibration drift. XII.A.3. Collisionality scans (regime interpolation test) Protocol. Execute a scan varying $\nu_*$ (V.C.1) primarily through density and/or temperature changes while holding geometry fixed and attempting to control turbulence indicators (as much as diagnostics permit). For each scan point, compute modeled $\nu_*^{(s)}(\rho)$ and compare the observed trend of inferred $\chi_{s,exp}(\rho)$ with the modeled trend of $\chi_{s,nc}(\rho)$ through the interpolation function $\mathcal{H}_{nc}(\nu_*)$ and any $S(\nu_*)$ coupling. Reject condition. If over a scan where $\nu_*$ changes by an order-one factor and the model asserts a monotone dependence (through a declared $\partial\chi/\partial\nu_*$ sign in a radius band), but the inferred $\chi_{s,exp}$ shows the opposite statistically significant trend in that band, reject the regime-interpolation structure (not merely the fitted coefficient magnitude). XII.B. Modulation experiments XII.B.1. Transfer functions for operator closure validation Protocol. With periodic modulation of absorbed power at one or more frequencies $\omega$, reconstruct the complex temperature response $\Theta_{exp}(\rho,\omega)$ (amplitude and phase) from ECE/Thomson or other profile diagnostics using a declared measurement operator $\mathcal{H}$ (XI.B.3). Construct the corresponding modeled response $\Theta_{model}$ by solving the resolvent equation (VI.C) or its discrete analog (X.8). Comparison. Define a complex residual vector stacking selected radii and frequencies, \[ \mathbf{r}_{TF} = \Sigma_{TF}^{-1/2}\bigl(\Theta_{exp}-\Theta_{model}\bigr), \] and use $\|\mathbf{r}_{TF}\|_2^2$ as a dedicated validation statistic. Reject condition. Reject the operator closure if the model cannot match both amplitude and phase profiles over a frequency set using the frozen $\boldsymbol{\theta}_{cal}^*$, particularly if the phase-lag slope in radius is systematically wrong across $\omega$, indicating incorrect radial structure in $K(\rho)$ and/or the edge parameter $h_{edge}$. XII.B.2. Heat-pulse propagation as a nested-surface transport check Protocol. Perform a power step or short pulse and track the time-of-arrival and decay of the temperature perturbation at multiple radii. Compare with predictions from the same discrete transport operator used for steady solutions (IX.E.4). The primary falsifiability objects are (i) the dominant decay time (VI.D.1) and (ii) the inferred radial propagation pattern implied by the operator eigenfunction expansion (VI.C). Reject condition. If the propagation signature indicates non-diffusive behavior incompatible with any admissible $K,M,H,h_{edge}$ in the declared closure class (e.g. clear evidence of rapid nonlocal coupling or strong stochastic transport not representable as a self-adjoint diffusion operator), reject the nested-surface diffusion assumption for that regime, triggering module replacement rather than coefficient refit. XII.C. $E_r$ tests: ambipolar residual and sign-change checks Protocol. Where CXRS (or equivalent) permits inference of $E_r(\rho)$, compare measured $E_r^{exp}(\rho)$ against the ambipolar prediction $E_r^{model}(\rho)$ produced by the closure (VII.B). Independently, evaluate and report the closure residual $\mathfrak{r}_{amb}(\rho)=|\mathcal{A}(E_r;\mathbf{x},\theta_{Er})|$ on the modeled solution. Reject conditions. 1. If $E_r$ is predicted within a declared interval $[E_{min},E_{max}]$ by a unique-root condition (VII.B.2) but the measured $E_r$ persistently lies outside this interval by more than declared uncertainty, reject the $E_r$ closure in that regime. 2. If the closure predicts a sign change in $E_r$ (or any sharp structural feature) across a scan and the measured $E_r$ does not exhibit it within uncertainty, reject the functional form of $\mathcal{A}$ or its input dependencies. XII.D. Impurity/ash/radiation tests Protocol. Use bolometry and spectroscopy-derived proxies (including $Z_{eff}$ where available) to compare the predicted radiation profile $p_{rad}(\rho)$ and integrated $P_{rad}$ (VII.E) against measurements, under controlled perturbations such as impurity seeding or wall-condition changes. Reject condition. If, after the single allowed calibration of the declared line-radiation parameter set $\theta_{line}$ on the calibration dataset, the model fails to predict the direction and magnitude of $\Delta P_{rad}$ and the associated profile response $\Delta T_e(\rho)$ under impurity perturbations, reject the radiation closure (as a structural error) rather than increasing calibration freedom. XII.E. Density-limit tests Protocol. Perform controlled density ramps at roughly fixed geometry and heating scheme. Track collapse/detachment proxies and compare to the predicted density limit margin $n_{max}-\bar n_e$ (VII.F.1). If the coupled model predicts multiple admissible steady branches (VII.F.2, VIII.B), then quasi-static ramps should be checked for hysteresis in a chosen monotone proxy (e.g. $W$ or central temperature). Reject condition. Reject the density-limit module if it predicts a sharp boundary (or branch multiplicity) in a parameter range where experiments robustly show no such transition (or vice versa), beyond declared uncertainty in $\bar n_e$, $P_{rad}$, and absorbed power. XII.F. Bootstrap/island/stochasticity tests Protocol. In $\beta$-scans or profile scans where $I_{bs}$ is expected to vary, compare predicted $\Delta\iota(\rho)$ and the derived stochasticity indicator $\mathcal{C}_{max}$ (VII.G) to iota-profile reconstructions and island/stochasticity diagnostics. Additionally, compare measured confinement degradation $\tau_E$ trends against the predicted penalty factor when a soft penalty is used. Reject condition. Reject the bootstrap–stochasticity coupling if a scan predicted to cross $\mathcal{C}_{max}=1$ by a significant margin shows no corresponding confinement degradation (or shows degradation when $\mathcal{C}_{max}\ll 1$ is predicted), indicating that either $\Delta\iota$ is wrong or that $\mathcal{C}_{max}$ is not a valid indicator in the tested regime. XII.G. SOL/edge tests Protocol. Use edge diagnostics (edge temperature, wall/divertor heat loads, connection-length proxies) to test the edge boundary map $h_{edge},T_{wall}\mapsto P_{edge}$ (VII.I). The required intermediate observables are $h_{edge}$ (as computed by the edge module), $T(1)$ from the core solve, and the implied $P_{edge}=h_{edge}A(1)(T(1)-T_{wall})$. Reject condition. Reject the edge module if it cannot reproduce the measured scaling of wall/divertor heat load proxies with $P_{SOL}$ and geometry changes using frozen parameters, or if inferred $h_{edge}$ would require unphysical values (e.g. negative, or orders-of-magnitude swings inconsistent with the declared surrogate structure) to match the data. XII.H. Vacuum-field/coil-model tests Protocol. Measure vacuum-field normal component metrics (IV.E.1) and their response to controlled perturbations (IV.E.2). Compare $\|B_n\|_{RMS}$ and its derivatives with respect to perturbation knobs against model predictions. Reject condition. Reject the geometry/coil proxy layer if it fails both baseline and perturbative trend tests, because subsequent transport predictions depend on geometry factors that would then be untrustworthy even if transport closures were otherwise adequate. XII.I. Explicit rejection thresholds and partial-failure protocol 1. Threshold declaration. For each module $\mathcal{M}_r$, define a validation statistic $S_r$ (chi-square on module observables, or a dedicated norm such as $\|\mathbf{r}_{TF}\|_2^2$) and a pre-declared rejection threshold $\tau_r$. A run rejects module $r$ if $S_r>\tau_r$ on the validation dataset. Thresholds must be declared before validation and should be consistent with the assumed noise model (XI.B). 2. Partial failure handling. If only a subset of modules fail (e.g. edge module fails but core diffusion closure passes transport inference), the protocol is: (i) freeze the passing modules and their calibrated coefficients; (ii) replace or augment only the failing module(s) by a new closure with the same interface (II.B); (iii) re-calibrate only the coefficients local to the replaced module using the original calibration set augmented only by diagnostics relevant to that module; (iv) repeat validation on the same held-out datasets. This preserves falsifiability by preventing cascading retuning across the pipeline. 3. Whole-model rejection. If failures occur across multiple independent falsifiability objects that are not attributable to a single module (e.g. simultaneous systematic mismatch in vacuum metrics, core transport inference, and modulation phase), then the coupled map $\mathcal{F}$ is rejected in the tested regime and must be revised structurally (e.g. through replacement of the geometry backbone or the nested-surface transport assumption). XIII. Limitations, Gaps, and Required External Verification (Explicit) This section documents limitations inherent in the reduced 0D/+1D approach and identifies which coefficients and assumptions must be externally verified or tightly calibrated before extrapolation to a compact Q>15 stellarator regime can be treated as credible. XIII.A. Unverified/placeholder coefficients requiring external confirmation Several closures above introduce dimensional constants or exponents (e.g. in $\chi_{nc}$, $\chi_{turb}$, radiation surrogates, density-limit surrogates, and any edge conductance model). In the present document these are treated as elements of $\boldsymbol{\theta}$ rather than as established physical constants unless they are fundamental unit conversion factors (e.g. $C_{keV}$). Therefore: 1. Any coefficient not fixed by dimensional analysis (III.D) must be either: (i) calibrated within the declared $\boldsymbol{\theta}_{cal}$ budget and then validated out-of-sample, or (ii) replaced by a higher-fidelity external calculation (e.g. tabulated neoclassical transport, atomic physics databases) at the same interface. 2. Extrapolation claims (e.g. from present devices to reactor regime) are only admissible when the similarity reductions (XI.E) are demonstrated on the validation dataset; otherwise, apparent agreement at one operating point is not evidence of correct scaling. XIII.B. Regime-of-validity risks in the reduced PDE/ODE structure 1. Nested-flux-surface assumption. The +1D backbone assumes transport can be represented as a flux-surface averaged radial operator. In regimes with strong stochasticity, magnetic islands, or 3D edge effects, this assumption may fail qualitatively. The required response is not to “tune” $\chi$, but to reject the nested-surface backbone in that regime and activate a module that explicitly penalizes/flags stochasticity (VII.G) or replaces the transport model. 2. Closure smoothness and differentiability. Exact sensitivities (Section X) assume differentiability of $\mathcal{R}$ and objectives. Hard switches (e.g. $\max(0,\cdot)$) can violate these assumptions and lead to ill-defined gradients. Any implementation that uses nonsmooth activations must either smooth them (with a declared smoothing scale) or explicitly report active sets and use one-sided derivatives; otherwise, gradient-based optimization results are not trustworthy. 3. Diagnostic reconstruction dependence. Several falsifiers require reconstructed quantities (absorbed power deposition, radiation profiles, $E_r$, $\iota(\rho)$). Reconstruction systematics enter $\Sigma$ and nuisance parameters (XI.B.2). If a module can only be made to “pass” by introducing many nuisance degrees of freedom, the falsifiability contract is violated; this should be treated as an inability to test the model with the available diagnostics rather than as evidence that the model is correct. XIII.C. Missing physics not resolved by the current minimal set The reduced system omits, or treats only as surrogates, several effects that may be decisive in compact high-performance stellarators. Examples include: 1. Full finite-$\beta$ equilibrium changes and their feedback on transport proxies beyond the simple $\iota$ parameterization. 2. Detailed neutral/atomic physics (recycling, charge exchange, and nonlocal radiation transport), which can dominate edge power balance and influence core profiles through boundary conditions. 3. Kinetic electron effects and multi-scale turbulence phenomena not represented by a critical-gradient surrogate, especially if turbulence becomes nonlocal or if multiple instabilities coexist. These omissions are not “bugs” but modeling scope choices. They become falsifiable when the corresponding intermediate observables (edge losses, radiation partition, modulation response shapes) systematically contradict data, indicating that the missing physics is active in the tested regime. XIII.D. Swap-ready upgrade map (interface invariants) The swap-ready requirement means upgrades must preserve interface contracts. The following invariants must not change when modules are replaced: 1. Geometry interface invariants. Upgraded equilibrium/geometry solvers may change how $V'(\rho)$, $A(\rho)$, and proxy profiles are computed, but must preserve the normalization conventions of $\rho$, the units specified in III.D, and the ability to evaluate volume integrals $\int_0^1 \langle f\rangle V’\,d\rho$. 2. Transport interface invariants. Any upgraded transport model must still produce auditable profile-level and flux-level observables (II.C), satisfy the discrete conservation identities (IX.B) in steady state, and provide consistent coefficient/operator outputs if operator-level falsifiers are claimed (VI, IX.E.4). 3. Sensitivity interface invariants. A replacement module must supply the partial derivatives needed by (X.4) and (X.5) (or declare that it is non-differentiable and thus incompatible with gradient-based optimization). If derivatives are approximate (e.g. from surrogate models), their approximation error must be characterized because it directly affects Fisher information (XI.C) and OED decisions. 4. Statistical wrapper invariants. Changes to physics modules must not alter the statistical protocol: calibration parameters must remain explicitly bounded and limited in number (II.D), measurement operators must remain explicit (XI.B.3), and rejection thresholds must be declared before validation (XII.I). These limitations and interface constraints are part of the falsifiability contract: they ensure that increased model fidelity improves predictive content rather than increasing the degrees of freedom available to fit any dataset. Conclusion A. Summary of the falsifiable 0D/+1D systems model structure This document formulated a deterministic parameter-to-observable map $\mathcal{F}:\mathbf{p}\mapsto\widehat{\mathbf{y}}$ by posing the coupled stellarator “systems model” as a residual problem $\mathcal{R}(\mathbf{u},\mathbf{p})=0$ together with an explicit constraint layer $\mathcal{C}(\mathbf{u},\mathbf{p})\le 0$. The state $\mathbf{u}$ includes profile unknowns (at minimum $T_e(\rho),T_i(\rho)$, optionally $n_e(\rho)$) and derived scalar performance metrics $(W,\tau_E,Q,\ldots)$. The architecture was specified as a composition of swap-ready modules with declared interfaces: a geometry/equilibrium proxy module $\mathcal{G}$ producing the geometric factors needed by transport, a core +1D transport/power-balance module $\mathcal{T}$, and additional physics/feasibility modules contributing residual equations, boundary data, or inequalities. A central design choice was the falsifiability contract: the model is required to output intermediate observables beyond the scalar gain $Q$, including profiles, fluxes, integrated powers, explicit feasibility margins, and (when used) operator-level decay/transfer-function quantities. Calibration was formalized as a constrained fit in a low-dimensional coefficient set followed by out-of-sample validation with pre-declared rejection criteria, so that disagreement can be localized to a module rather than absorbed by unconstrained tuning. B. Summary of transport-operator closures and their falsifiable content For the core +1D backbone, the paper defined flux-surface averaged steady energy balances in conservative form and a base diffusive closure $q_s=-n_s\chi_s\,dT_s/dr$, including an explicit quadrature mapping from sources $p_s(\rho)$ to integrated power $P_s(\rho)$ and temperature profiles under stated positivity assumptions. This yielded a direct falsifier: from measured profiles and reconstructed sources, one can infer an experimental $\chi_{exp}(\rho)$ and compare it to the model-implied $\chi_{nc}+\chi_{turb}$ structure. An additive neoclassical+turbulent surrogate was specified with explicit unit-consistent parameterizations and a collisionality-dependent modulation of turbulent transport. In addition, a fast global closure was given in which steady balance reduces to a scalar root equation $P-\mathcal{L}(P)=0$ under continuity and monotonicity hypotheses, together with sufficient conditions for existence/uniqueness that can be checked in computation and falsified by power scans. The transport-operator closure re-expressed the (linearized) time-dependent transport equation as a self-adjoint Sturm–Liouville-type operator $\mathcal{L}$ on a weighted $L^2$ space, with Robin edge loss encoded by an effective conductance $h_{edge}$. The principal eigenvalue $\lambda_1$ was identified as a predicted decay rate ($\tau_{op}=1/\lambda_1$ when $\lambda_1>0$), endowed with Rayleigh-quotient bounds and monotonicity properties. These operator-level objects yield additional, near-term falsifiers: step-response decay times and modulation transfer functions $\Theta(\rho,\omega)$ must match measurements, and controlled parameter scans must obey sign-definite eigenvalue trends implied by the min–max principle. C. Summary of sensitivity, adjoint, and OED framework Exact sensitivity analysis was derived at the discretized level via the implicit-function theorem and adjoint equations. For any scalar objective $\Psi(\mathbf{u},\mathbf{p})$ (including $Q$, $\tau_E$, radiated power, constraint margins, or $\lambda_1$), the adjoint system $\mathcal{J}^T\boldsymbol{\lambda}=(\partial_{\mathbf{u}}\Psi)^T$ with $\mathcal{J}=\partial_{\mathbf{u}}\mathcal{R}$ yields gradients in the swap-ready form $d\psi/d\mathbf{p}=\partial_{\mathbf{p}}\Psi-\boldsymbol{\lambda}^T\partial_{\mathbf{p}}\mathcal{R}$. Operator-closure sensitivities were recorded through Hellmann–Feynman identities for generalized eigenproblems and resolvent differentiation for modulation response. On the statistical side, the deterministic map was wrapped in an explicit measurement-operator model $\mathcal{H}$, a Gaussian misfit functional, and chi-square rejection tests, with nuisance-parameter handling constrained to remain auditable and low-dimensional. Fisher information was used to quantify identifiability and to support optimal experimental design via standard scalar criteria. A dimensionless-similarity formalism was included to make extrapolation claims themselves falsifiable: if predicted invariances at fixed dimensionless groups do not hold within uncertainty across near-similar conditions, the asserted similarity reduction is rejected. D. Near-term experimental falsification path and decision logic The integrated falsification plan assembled concrete validation protocols using measurements feasible on present stellarators: power scans testing the scalar-root closure, profile-based inference of $\chi_{exp}(\rho)$, collisionality scans testing regime interpolation, modulation/heat-pulse experiments testing operator transfer functions and $\lambda_1$, and targeted checks of additional modules (ambipolar $E_r$ residuals, radiation closure under impurity perturbations, density-limit ramps and branch structure, bootstrap/island indicators and confinement penalties, SOL/edge closure via $h_{edge}$, and vacuum-field response metrics for the geometry/coil proxy layer). The decision logic was made explicit: if a module fails its own falsifiability objects on held-out data, it is to be replaced at the same interface rather than globally re-tuned; if multiple independent falsifiers fail in a way that cannot be localized, the coupled map is rejected in that regime. Finally, the limitations section emphasized that many coefficients and regime assumptions in the reduced closures remain placeholders unless externally verified or tightly calibrated, and that the nested-flux-surface backbone and differentiability requirements impose concrete regimes of validity. The paper’s contribution is therefore not a claim of verified Q>15 performance, but a mathematically explicit, interface-stable, and rejection-oriented framework in which each modeling assumption produces intermediate observables and quantitative tests that can invalidate it on a realistic experimental timescale. [HARD CODED END-OF-PAPER MARK — ALL CONTENT SHOULD BE ABOVE THIS LINE] ================================================================================ MODEL CREDITS This autonomous solution attempt was generated with the Intrafere LLC AI Harness, MOTO, and the following model(s): – x-ai/grok-4 (40 API calls) – openai/gpt-5.2 (33 API calls) – z-ai/glm-5 (27 API calls) – moonshotai/kimi-k2.5 (1 API calls) Total AI Model API Calls: 101 ================================================================================

More Info