1. Introduction

1.1 The Navier–Stokes Millennium Problem

The Navier–Stokes existence and smoothness problem is one of the seven Millennium Prize Problems identified by the Clay Mathematics Institute in the year 2000 [1]. The problem asks:

Given smooth, divergence-free initial velocity \( u_0 : \mathbb{R}^3 \to \mathbb{R}^3 \) with finite energy \( \int_{\mathbb{R}^3} |u_0|^2 \, dx < \infty \), does there exist a smooth solution \( u : \mathbb{R}^3 \times [0, \infty) \to \mathbb{R}^3 \) to the incompressible Navier–Stokes equations \[ \frac{\partial u}{\partial t} + (u \cdot \nabla) u = -\nabla p + \nu \Delta u, \qquad \nabla \cdot u = 0 \tag{1} \] for all \( t > 0 \)?

The question has resisted resolution for over eight decades, dating back to Leray's foundational work in 1934 [2], which established the existence of weak solutions (now called Leray–Hopf weak solutions) but left their uniqueness and smoothness unresolved. In two dimensions, Ladyzhenskaya [3] proved global existence and uniqueness of smooth solutions in 1969, exploiting the fact that vorticity is a scalar transported by the flow without stretching. In three dimensions, the situation is fundamentally different: the vortex stretching term couples velocity and vorticity in a way that permits potential finite-time singularity formation.

The history of the problem can be traced through several landmark results. Leray [2] showed that weak solutions exist globally but may develop singularities, and that the set of singular times has one-dimensional Hausdorff measure zero. Prodi [19] and Serrin [20] established conditional regularity criteria: solutions in \( L^p_t L^q_x \) with \( 2/p + 3/q \leq 1 \) are smooth, providing the first systematic understanding of which integrability conditions prevent singularity formation. Caffarelli, Kohn, and Nirenberg [4] proved that the set of space-time singularities has one-dimensional parabolic Hausdorff measure zero — the celebrated partial regularity theorem. Escauriaza, Seregin, and Šverák [18] extended this to prove that \( L^3 \)-solutions cannot develop singularities at the endpoint of the Prodi–Serrin scale, resolving a long-standing gap. Beale, Kato, and Majda [5] established that the solution remains smooth on \( [0, T] \) if and only if \[ \int_0^T \left\|\omega(\cdot, t)\right\|_{L^\infty} \, dt < \infty, \tag{2} \] where \( \omega = \nabla \times u \) is the vorticity. This criterion reduced the regularity question to a single integral condition on the maximum vorticity, but controlling that integral has proven extraordinarily difficult.

More recently, Tao [6] showed that an averaged version of the 3D Navier–Stokes equations — obtained by replacing the bilinear form \( B(u, u) \) with a carefully engineered operator that preserves the energy identity but breaks the cancellation structure — admits solutions that develop a finite-time singularity. Tao's construction is based on a “fluid computer” that programs the nonlinearity to build a self-similar cascade leading to singularity. Any argument about the true Navier–Stokes equations must therefore use specific structural properties of the actual nonlinearity, not merely energy-level estimates.

1.2 The Vortex Stretching Obstacle

The central difficulty is the vortex stretching term. Taking the curl of (1), we obtain the vorticity equation: \[ \frac{\partial \omega}{\partial t} + (u \cdot \nabla) \omega = (\omega \cdot \nabla) u + \nu \Delta \omega. \tag{3} \]

The term \( (\omega \cdot \nabla) u \) is the vortex stretching term. It has no analogue in two dimensions (where \( \omega \) is a scalar and the term vanishes identically). In three dimensions, it represents the tilting and stretching of vortex lines by the velocity gradient. Physically, as a vortex tube is stretched along its axis, conservation of angular momentum demands that the tube's cross-section shrink and the vorticity intensify. This is the mechanism by which tornadoes, waterspouts, and bathtub vortices concentrate angular momentum.

Mathematically, the stretching term creates a coupling between \( \omega \) and \( \nabla u \) that is supercritical with respect to the natural energy scaling. In the enstrophy evolution equation \[ \frac{d}{dt} \frac{1}{2} \int |\omega|^2 \, dx = -\nu \int |\nabla \omega|^2 \, dx + \int \omega_i \omega_j S_{ij} \, dx, \tag{4} \] where \( S_{ij} = \frac{1}{2}(\partial_i u_j + \partial_j u_i) \) is the strain rate tensor, the production term \( \int \omega_i \omega_j S_{ij} \, dx \) is cubic in the velocity gradient. The diffusion term \( -\nu \int |\nabla \omega|^2 \, dx \) is quadratic but involves higher-order derivatives. The competition between these two terms determines whether enstrophy remains bounded or blows up. Foias, Manley, Rosa, and Temam [17] provide a comprehensive treatment; Doering [16] gives a modern survey.

1.3 Survey of Existing Computational Approaches

Direct Numerical Simulation (DNS). DNS resolves all scales of turbulence by discretising the full equations on a grid fine enough to capture the Kolmogorov microscale \( \eta = (\nu^3/\varepsilon)^{1/4} \). The computational cost scales as \( Re^{9/4} \) in three dimensions [11], making DNS infeasible for high Reynolds numbers.

Large Eddy Simulation (LES). LES resolves only the large-scale motions and models the subgrid-scale stresses through closure models such as the Smagorinsky model [12]. LES is computationally cheaper than DNS but introduces modelling errors.

Reynolds-Averaged Navier–Stokes (RANS). RANS decomposes all quantities into mean and fluctuating parts and solves for the mean flow with turbulence models (\( k \)-\( \varepsilon \), \( k \)-\( \omega \), Reynolds stress transport) providing closure [11].

Lattice Boltzmann Method (LBM). LBM evolves particle distribution functions on a lattice, recovering Navier–Stokes dynamics in the continuum limit through the Chapman–Enskog expansion [13]. Our framework incorporates the LBM concept of relaxation deficit (Level 18) as a diagnostic.

Galerkin Truncation. Galerkin projection onto a finite set of Fourier modes preserves the quadratic structure of the nonlinearity and the energy identity [8]. The 3D kernel of Section 2.7 is built this way. The shell system of Sections 2.3–2.6 is not: it is a hand-constructed ODE system on scalar wavenumbers.

Tao's Fluid Computer. Tao [6, 7] proposed viewing the Navier–Stokes nonlinearity as a “programmable” system that can be engineered to perform arbitrary computations. Our framework's \( H'' \) diagnostic, which measures prediction confidence, is directly motivated by this connection.

1.4 Motivation

Classical energy estimates seek an inequality of the form \[ \frac{d}{dt} \Omega \leq -c_1 \nu \Omega^{1+\delta} + c_2 \Omega^{3/2} \] where \( \Omega = \int |\omega|^2 \, dx \) is the enstrophy, and treat each quantity in isolation. The diagnostic \( H \) of this paper instead records many coupled trajectory statistics of a low-dimensional shell model at once. The construction borrows its composition rules from the Unified Adaptation Theorem [21].

1.5 How Our Scaffold Differs from Classical Estimates

The diagnostic differs from a classical energy estimate in three ways:

1. Coupled monitoring. Classical estimates bound individual quantities. \( H \) records 26 trajectory statistics of the shell model and their couplings.
2. Trajectory statistics. Classical estimates produce one-sided inequalities. \( H \) tracks the coupling among enstrophy \( L_1 \), stability margin \( L_4 \) and convergence score \( L_2 \) along a computed trajectory.
3. Dual-number derivatives. Each quantity is computed as a dual number — value plus derivative — via forward-mode automatic differentiation.

\( H \) is a diagnostic of the shell model's trajectories, not an estimate on the Navier–Stokes equations.

2. The Shell Model and the First 3D Kernel

2.1 The Full Navier–Stokes Vorticity Equation

We begin with the three-dimensional incompressible Navier–Stokes equations in vorticity form. Let \( u = (u_1, u_2, u_3) \) be the velocity field and \( \omega = \nabla \times u = (\omega_1, \omega_2, \omega_3) \) the vorticity. The vorticity equation is:

\[ \frac{\partial \omega}{\partial t} + (u \cdot \nabla) \omega = (\omega \cdot \nabla) u + \nu \Delta \omega. \tag{6} \]

Expanding in components, this is a system of three coupled PDEs involving advection, vortex stretching and tilting, and viscous diffusion. The stretching terms couple all three components of \( \omega \) to all three components of \( \nabla u \), creating a nine-component tensor interaction.

2.2 Fourier Expansion and the Two Systems

On the periodic domain \( [0, 2\pi]^3 \) each component of velocity and vorticity has a Fourier series:

\[ u(x, t) = \sum_{k \in \mathbb{Z}^3} \hat{u}(k, t) e^{ik \cdot x}, \qquad \omega(x, t) = \sum_{k \in \mathbb{Z}^3} \hat{\omega}(k, t) e^{ik \cdot x}. \]

A Galerkin truncation retains the modes with \( |k| \leq K \), keeps the quadratic nonlinearity, and conserves \( \sum_k |\hat{u}(k)|^2 \) in the inviscid limit. The 3D kernel of Section 2.7 is a truncation of this kind. The shell system of Sections 2.3–2.6 is a different object: a hand-constructed cubic ODE system on scalar wavenumbers \( k_i = i \), with a relaxation term, a cubic term and a signed cascade term. It is not a Galerkin truncation of Navier–Stokes and does not conserve energy in the inviscid limit.

2.3 The Coupling Structure

Velocity–Vorticity Coupling (\( \sigma \)).

\[ \frac{du_i}{dt} = \sigma(\omega_i - u_i) - \nu k_i^2 u_i. \tag{7} \]

Linear Vortex Stretching (\( \lambda \)).

\[ F_{\text{stretch}}^{(i)} = \lambda \sum_j c_{ij} \omega_j u_j. \tag{8} \]

Quadratic Stretching (\( \lambda_2 \)).

\[ F_{\text{quad}}^{(i)} = \lambda_2 |\omega|^2 \omega_i. \tag{9} \]

Forward Cascade (\( \beta \)).

\[ F_{\text{cascade}}^{(i)} = \sum_{j < i} \beta_{ji} \omega_j^2 \cdot \text{sgn}(\omega_i). \tag{10} \]

2.4 The 6-Mode Shell Model

The smallest shell model with nontrivial dynamics has 6 degrees of freedom: 3 velocity components \( (u_1, u_2, u_3) \) and 3 vorticity components \( (\omega_1, \omega_2, \omega_3) \) at wavenumbers \( k_1 = 1 \), \( k_2 = 2 \), \( k_3 = 3 \). The complete system:

\[ \begin{aligned} \frac{du_1}{dt} &= \sigma(\omega_1 - u_1) - \nu \cdot 1^2 \cdot u_1, \\ \frac{du_2}{dt} &= \sigma(\omega_2 - u_2) - \nu \cdot 2^2 \cdot u_2, \\ \frac{du_3}{dt} &= \sigma(\omega_3 - u_3) - \nu \cdot 3^2 \cdot u_3, \\ \frac{d\omega_1}{dt} &= -\nu \omega_1 + \lambda (c_{11} \omega_1 u_1 + c_{12} \omega_2 u_2 + c_{13} \omega_3 u_3) + \lambda_2 |\omega|^2 \omega_1, \\ \frac{d\omega_2}{dt} &= -4\nu \omega_2 + \lambda (c_{21} \omega_1 u_1 + c_{22} \omega_2 u_2 + c_{23} \omega_3 u_3) + \lambda_2 |\omega|^2 \omega_2 + \beta_{12} \omega_1^2 \cdot \text{sgn}(\omega_2), \\ \frac{d\omega_3}{dt} &= -9\nu \omega_3 + \lambda (c_{31} \omega_1 u_1 + c_{32} \omega_2 u_2 + c_{33} \omega_3 u_3) + \lambda_2 |\omega|^2 \omega_3 + \beta_{13} \omega_1^2 \cdot \text{sgn}(\omega_3) + \beta_{23} \omega_2^2 \cdot \text{sgn}(\omega_3). \end{aligned} \]

Here \( |\omega|^2 = \omega_1^2 + \omega_2^2 + \omega_3^2 \). The coupling coefficients satisfy \( c_{ij} = k_j / k_i \) for triad-connected modes. The cascade coefficients \( \beta_{ij} = \beta_0 k_j / k_i \).

2.5 General \( N \)-Mode Shell Model

For the general \( N \)-mode shell model with \( K = N/2 \) wavenumbers:

\[ \begin{aligned} \frac{du_i}{dt} &= \sigma(\omega_i - u_i) - \nu k_i^2 u_i, \\ \frac{d\omega_i}{dt} &= -\nu k_i^2 \omega_i + \lambda \sum_{j=1}^{K} c_{ij} \omega_j u_j + \lambda_2 |\omega|^2 \omega_i + \sum_{j < i} \beta_{ji} \omega_j^2 \cdot \text{sgn}(\omega_i), \end{aligned} \tag{11} \]

for \( i = 1, \ldots, K \), with wavenumbers \( k_i = i \) and \( N = 2K \) total ODEs.

2.6 Initial Conditions

The standard initial condition is sinusoidal with exponential decay:

\[ u_i(0) = A \cdot \frac{\sin(k_i)}{k_i}, \qquad \omega_i(0) = A \sin(k_i), \tag{12} \]

where \( A \) is the amplitude parameter.

2.7 The First 3D Kernel

Alongside the shell model, this record contains a first Fourier–Galerkin solver for the 3D equations on the periodic torus \( \mathbb{T}^3 \) with:

• All integer wavevectors \( \mathbf{k} = (k_1, k_2, k_3) \) with \( |\mathbf{k}| \leq N_{\max} \)
• Coupling coefficients from the trilinear form \( B(\mathbf{k}, \mathbf{p}, \mathbf{q}) = -(\hat{u}_\mathbf{p} \cdot \mathbf{q})\, P_\mathbf{k}(\hat{u}_\mathbf{q}) \)
• Leray projection \( P_\mathbf{k}(\mathbf{v}) = \mathbf{v} - (\mathbf{k} \cdot \mathbf{v} / |\mathbf{k}|^2)\mathbf{k} \) enforcing incompressibility
• Divergence-free initial conditions with energy decaying as \( 1/|\mathbf{k}| \)
• All triads \( (\mathbf{p}, \mathbf{q}) \) with \( \mathbf{p} + \mathbf{q} = \mathbf{k} \) precomputed and cached

Mode and triad counts of the 3D kernel at each truncation.

\( N_{\max} \)	Modes	Triads
2	32	426
3	122	6,642
4	256	30,360
5	514	122,472

This kernel stores real coefficients and omits the factor \( -i \) in the Fourier nonlinearity, so it injects energy; the defect was found and corrected on 2026-03-24. The mode enumeration, triad precomputation and Leray projection carry over unchanged to the corrected complex-coefficient solver. Every 3D excursion amplitude in this record was computed with the real-coefficient kernel and is reported only as a record of that computation.

3. The Diagnostic \( H \)

3.1 The 26 Diagnostic Levels

The holistic diagnostic \( H \) is built from 26 levels \( L_0, L_1, \ldots, L_{25} \), each measuring a different statistic of a shell-model trajectory. \( H \) is an in-sample classifier: its weights (Section 3.5) are learned from the outcomes it reports. Every level is computed as a dual number \( (v, v') \) where \( v \) is the value and \( v' = dv/d\theta \) is the derivative with respect to the state parameter \( \theta \).

We now define each level. For brevity we state the name, equation, physical meaning, and gate condition for each.

Level 0: State Vector \( \theta \). \( L_0(t) = (\theta(t), \dot\theta(t)) \), where \( \theta = (u_1, \ldots, u_K, \omega_1, \ldots, \omega_K) \in \mathbb{R}^{2K} \). Always active.

Level 1: Normalised Enstrophy.

\[ L_1(t) = \frac{\Omega(t)}{\Omega(0)}, \qquad \Omega(t) = \sum_{i=1}^{K} \omega_i(t)^2. \tag{13} \]

\( L_1 = 1 \) means enstrophy is at its initial level; \( L_1 > 100 \) is the excursion criterion of the protocol. Always active.

Level 2: Convergence Score.

\[ L_2(t) = S(t) = 1 - \frac{D_{\text{late}}(t)}{D_{\text{early}}(t)}. \tag{14} \]

\( S \to 1 \): converging to steady state. \( S < 0 \): diverging. Gate: \( t > t_{\min} \).

Level 3: PID Control on Convergence. Proportional \( S(t) \), Integral \( \frac{1}{t}\int_0^t S(\tau)\,d\tau \), Derivative \( S'(t) \). Gate: \( t > 2t_{\min} \).

Level 4: Stability Margin.

\[ L_4(t) = M(\theta, \varepsilon) = \frac{\left\|\theta_\varepsilon(t) - \theta(t)\right\|}{\varepsilon}. \tag{15} \]

Amplification factor of a small perturbation. Always active.

Level 5: Fragility Rate. \( L_5(t) = dM/dt \). Detects onset of instability.

Level 6: Perturbation Scaling. \( L_6(t) = dM/d\varepsilon \). Measures superlinear perturbation response. Gate: \( M > 2 \).

Level 7: Optimal Probe. \( L_7 = \varepsilon^* = \arg\min_\varepsilon \text{Var}[M(\theta, \varepsilon)] \). Learned via meta-gradient. Gate: \( t > 5{,}000 \).

Level 8: Parameter Sensitivity. \( L_8 = (\partial T_{\text{blow}}/\partial A,\; \partial T_{\text{blow}}/\partial \nu) \). Gate: \( t = T \).

Level 9: Coupling Sensitivity. \( L_9 = (\partial T_{\text{blow}}/\partial \lambda,\; \partial^2 T_{\text{blow}}/\partial \lambda \partial A) \). Gate: \( t = T \).

Level 10: Trajectory Acceleration. \( L_{10}(t) = \partial T_{\text{blow}} / \partial \alpha_\nu \). Response to time-varying viscosity. Gate: sweep active.

Level 11: Enstrophy Acceleration.

\[ L_{11}(t) = \frac{d^2 \Omega}{dt^2} = 2 \sum_i (\dot\omega_i^2 + \omega_i \ddot\omega_i). \tag{16} \]

Gate: \( d\Omega/dt > 0 \).

Level 12: Deceleration Ratio. \( L_{12}(t) = (d^2\Omega/dt^2)/(d\Omega/dt)^2 \). Distinguishes self-similar runaway growth from exponential growth.

Level 13: Enstrophy Jerk. \( L_{13}(t) = d^3\Omega/dt^3 \). Triple-positive signature \( (L_{11}, L_{13}, d\Omega/dt \text{ all } > 0) \) indicates imminent runaway growth.

Level 14: Saturation Time Predictor. \( L_{14}(t) = -\frac{d\Omega/dt}{d^2\Omega/dt^2} \). Predicted time to enstrophy saturation. Gate: \( d^2\Omega/dt^2 < 0 \).

Level 15: Enstrophy Doubling Time.

\[ \tau_d(n) = t_n - t_{n-1}, \qquad \text{where } \Omega(t_n) = 2^n \cdot \Omega(0). \tag{17} \]

Ratio \( \tau_d(n)/\tau_d(n-1) > 1 \): safe. \( < 0.85 \): runaway growth. Gate: 2+ doublings.

Level 16: Reynolds Fluctuation Ratio. \( L_{16}(t) = |\Omega(t) - \langle\Omega\rangle_W| / \langle\Omega\rangle_W \). Inspired by RANS. Gate: \( t > W \).

Level 17: Turbulent Timescale Ratio. \( L_{17}(t) = (k(t)/\varepsilon(t))/\tau_d(t) \). Compares \( k \)-\( \varepsilon \) timescale with doubling time.

Level 18: Relaxation Deficit. \( L_{18}(t) = |F_{\text{stretch}}(t)| / |F_{\text{diff}}(t)| \). Inspired by LBM. \( L_{18} = 1 \) is exact balance. Always active.

Level 19: Symmetry Breaking.

\[ L_{19}(t) = \frac{|\omega_x - \omega_y| + |\omega_y - \omega_z| + |\omega_z - \omega_x|}{|\omega_x| + |\omega_y| + |\omega_z|}. \tag{18} \]

\( L_{19} = 0 \): perfect isotropy. \( L_{19} \to 2 \): extreme anisotropy. Always active.

Level 20: Directional Alignment. \( L_{20}(t) = |\dot\omega \cdot \omega| / (|\dot\omega| \cdot |\omega|) \). Cosine of angle between vorticity and its rate of change. \( L_{20} \approx 1 \): dimension collapse (vortex tube stretching).

Level 21: Stretching Efficiency.

\[ L_{21}(t) = \frac{\omega^T S \omega}{|\omega|^2 \cdot |S|}. \tag{19} \]

Alignment between vorticity and strain. Earliest discriminator between bounded and runaway trajectories in the shell model (diverges at step ~6,000). Always active.

Level 22: Back-Reaction. \( L_{22}(t) = \sum_i |u_i(t) - u_i^{(\text{passive})}(t)| / \sum_i |u_i(t)| \). How much vorticity dynamics modified velocity. Gate: \( t > 1{,}000 \).

Level 23: Energy Transfer Flux. \( L_{23}(t) = \Pi(k_c, t) \). Energy flux across cutoff wavenumber \( k_c = K/2 \). Gate: \( K \geq 4 \).

Level 24: Phase Coherence. \( L_{24}(t) = |\sum_i \omega_i(t) e^{ik_i\phi(t)}| / \sum_i |\omega_i(t)| \). Degree of phase-locking across modes. Gate: \( K \geq 3 \).

Level 25: Spectral Slope.

\[ L_{25}(t) = -\frac{d \ln E(k)}{d \ln k}\bigg|_{\text{fit}}. \tag{20} \]

In Kolmogorov inertial range, slope is \( 5/3 \approx 1.67 \). Shallower slope indicates energy accumulation at small scales. Gate: \( K \geq 4 \).

Summary of the 26 Diagnostic Levels

Level	Name	Equation	Dual Part Meaning	Gate
\( L_0 \)	State vector	\( \theta = (u, \omega) \)	Velocity in phase space	Always
\( L_1 \)	Enstrophy	\( \Omega / \Omega(0) \)	Growth rate	Always
\( L_2 \)	Convergence	\( 1 - D_{\text{late}}/D_{\text{early}} \)	Convergence trend	\( t > t_{\min} \)
\( L_3 \)	PID control	\( (S, \bar{S}, S') \)	Convergence dynamics	\( t > 2t_{\min} \)
\( L_4 \)	Stability margin	\( \|\delta\theta\|/\varepsilon \)	Fragility rate	Always
\( L_5 \)	Fragility rate	\( dM/dt \)	Fragility acceleration	Always
\( L_6 \)	Perturbation scaling	\( dM/d\varepsilon \)	Nonlinearity trend	\( M > 2 \)
\( L_7 \)	Optimal probe	\( \varepsilon^* \)	Probe adaptation rate	\( t > 5000 \)
\( L_8 \)	Parameter sensitivity	\( \partial T/\partial (A, \nu) \)	Sensitivity trend	\( t = T \)
\( L_9 \)	Coupling sensitivity	\( \partial T/\partial \lambda \), cross	Interaction dynamics	\( t = T \)
\( L_{10} \)	Trajectory accel.	\( \partial T/\partial \alpha_\nu \)	Sweep effectiveness	Sweep active
\( L_{11} \)	Enstrophy accel.	\( d^2\Omega/dt^2 \)	Jerk	\( \dot\Omega > 0 \)
\( L_{12} \)	Deceleration ratio	\( \ddot\Omega / \dot\Omega^2 \)	Runaway type	\( \dot\Omega > 0 \)
\( L_{13} \)	Jerk	\( d^3\Omega/dt^3 \)	Snap	\( \dot\Omega > 0 \)
\( L_{14} \)	Saturation time	\( -\dot\Omega/\ddot\Omega \)	Predicted settling	\( \ddot\Omega < 0 \)
\( L_{15} \)	Doubling time	\( \tau_d \) and ratio	Trend	2+ doublings
\( L_{16} \)	Reynolds ratio	\( |\Omega - \bar\Omega|/\bar\Omega \)	Fluctuation trend	\( t > W \)
\( L_{17} \)	Timescale ratio	\( (k/\varepsilon)/\tau_d \)	Equilibration	\( \tau_d \) defined
\( L_{18} \)	Relaxation deficit	\( |F_{\text{str}}|/|F_{\text{diff}}| \)	Balance trend	Always
\( L_{19} \)	Symmetry breaking	\( \sum|\omega_i - \omega_j|/\sum|\omega_i| \)	Anisotropy growth	Always
\( L_{20} \)	Alignment	\( |\dot\omega \cdot \omega|/(|\dot\omega||\omega|) \)	Collapse trend	\( |\omega| \) large
\( L_{21} \)	Stretching eff.	\( \omega^T S \omega/(|\omega|^2|S|) \)	Alignment trend	Always
\( L_{22} \)	Back-reaction	\( |u - u^{\text{pass}}|/|u| \)	Growth trend	\( t > 1000 \)
\( L_{23} \)	Energy flux	\( \Pi(k_c) \)	Flux trend	\( K \geq 4 \)
\( L_{24} \)	Phase coherence	\( |\sum \omega_i e^{ik_i\phi}|/\sum|\omega_i| \)	Coherence trend	\( K \geq 3 \)
\( L_{25} \)	Spectral slope	\( -d\ln E/d\ln k \)	Slope evolution	\( K \geq 4 \)

3.2 The Feedback Loop

The enstrophy production equation is the starting point:

\[ \frac{d\Omega}{dt} = \underbrace{-2\nu \sum_i k_i^2 \omega_i^2}_{D(t)} + \underbrace{2\lambda \sum_i \omega_i \sum_j c_{ij} \omega_j u_j}_{P_1(t)} + \underbrace{2\lambda_2 |\omega|^2 \sum_i \omega_i^2}_{P_2(t)} + \underbrace{2\sum_i \omega_i F_{\text{cascade}}^{(i)}}_{C(t)}. \tag{21} \]

In the shell model \( P_2(t) = 2\lambda_2 \Omega^2 \) grows as the square of enstrophy, while \( D(t) \) grows only linearly; the cubic term is what makes runaway growth possible in this ODE system.

The coupling among \( L_1 \) (enstrophy), \( L_4 \) (stability margin) and \( L_2 \) (convergence score) is summarised by a loop gain \( G = (\partial L_4/\partial L_1) \cdot (\partial L_2/\partial L_4) \cdot (\partial L_1/\partial L_2) \), measured from the dual parts at enstrophy doubling. The following table records it for the 8-mode shell model.

Loop Gain as a Function of Amplitude

Loop gain \( G(A) \) for the 8-mode shell model at various amplitudes; \( A_T = 0.290 \) is the finite-window excursion amplitude of Section 3.3.

\( A \)	\( A/A_T \)	\( \partial L_4/\partial L_1 \)	\( \partial L_2/\partial L_4 \)	\( \partial L_1/\partial L_2 \)	\( G \)
0.15	0.52	0.82	−0.31	−0.45	0.11
0.20	0.69	1.24	−0.48	−0.58	0.35
0.25	0.86	1.85	−0.62	−0.73	0.84
0.28	0.97	2.15	−0.72	−0.81	1.25
0.30	1.03	2.40	−0.78	−0.85	1.59
0.35	1.21	3.10	−0.91	−0.94	2.65

3.3 The Finite-Window Excursion Amplitude

\( A_T \) is a property of the protocol, not of the equations: it depends on the window length (\( A_T \propto T^{-0.89} \)), on the excursion multiple (here 100), and on the solver. It is not a regularity threshold, and the tables in Sections 5 and 7 report it only as a record of what each system did under the stated protocol.

3.4 Time Gating

Each level \( L_i \) has a gate function \( g_i(t) \in \{0, 1\} \). The holistic diagnostic is:

\[ H(t) = \frac{\sum_{i=0}^{25} g_i(t) \cdot w_i \cdot s_i(t)}{\sum_{i=0}^{25} g_i(t) \cdot w_i}. \tag{24} \]

The denominator normalises by the active weight sum, ensuring \( H \in [0, 1] \) regardless of how many gates are open.

3.5 Self-Adapting Weights (\( H' \))

The weights are learned online:

\[ w_i^{(n+1)} = w_i^{(n)} + \eta \cdot \text{acc}^{(n)} \cdot s_i^{(n)}. \tag{25} \]

This Hebbian-like rule upweights levels that are active during correct predictions. The “correct” label is the outcome of the same trajectory, so the weights are fitted in-sample and the classifier is evaluated on the data that trained it.

Learned weights (most influential): \( L_{15} \) (doubling time): 0.35; \( L_{21} \) (stretching efficiency): 0.22; \( L_{18} \) (relaxation deficit): 0.15; \( L_{14} \) (saturation time): 0.12; \( L_4 \) (stability margin): 0.08.

3.6 Confidence Tracking (\( H'' \))

\( H'' = dC/dt \). Its value is reported at \( T = 50{,}000 \) and \( T = 100{,}000 \) steps in the shell-model runs.

4. The Enstrophy Doubling-Time Lemma

The lemma is elementary and is the one statement of this record that stands as a theorem. It says nothing about whether a given system's doubling times shrink.

Classification Accuracy vs Threshold \( c \)

Classification accuracy of the shrinking-doubling-time rule vs threshold \( c \) (6-mode shell model, in-sample).

Threshold \( c \)	True Positive Rate	True Negative Rate	Accuracy
0.70	100%	72%	80.0%
0.75	100%	78%	83.3%
0.80	100%	81%	85.0%
0.85	100%	82%	86.1%
0.90	95%	86%	83.3%
0.95	87%	90%	80.0%

5. The Viscosity Sweep

5.1 The Solid–Fluid–Gas Spectrum

Viscosity sweep (8-mode shell model, \( A = 0.3 \)). \( A_T \) is the finite-window excursion amplitude of Section 3.3.

\( \nu \)	Regime	\( A_T \)	\( \max\Omega/\Omega(0) \)	Outcome
\( 10^{-4} \)	Gas	0.052	>100	Excursion
\( 5 \times 10^{-4} \)	Gas	0.121	>100	Excursion
\( 10^{-3} \)	Low fluid	0.184	>100	Excursion
\( 5 \times 10^{-3} \)	Fluid	0.290	>100	Excursion (marginal)
\( 10^{-2} \)	Fluid	0.412	45.2	Bounded
\( 5 \times 10^{-2} \)	Viscous	0.783	8.4	Bounded
\( 10^{-1} \)	Solid	1.240	2.1	Bounded

In the shell model \( A_T \) increases monotonically with \( \nu \) over this sweep, with an empirical fit \( A_T \propto \nu^{0.4} \).

5.2 Resonance Discovery

When viscosity is swept periodically (\( \nu(t) = \nu_0 + \Delta\nu \sin(2\pi t/T_\nu) \)), a resonance emerges at \( T_\nu \approx 10{,}000 \) steps, with enstrophy response amplified 3–5×. The resonance period:

\[ T_{\text{res}} = \frac{2\pi}{\omega_{\text{nat}}}, \qquad \omega_{\text{nat}} = \sqrt{2\lambda_2 \Omega_{\text{eq}} \cdot 2\nu k_1^2 - (\nu k_1^2)^2}. \tag{27} \]

5.3 Hysteresis and Ratchet

Sweeping viscosity upward then downward reveals hysteresis with ratio \( R_H = \Omega_{\text{down}}(\nu_{\text{mid}})/\Omega_{\text{up}}(\nu_{\text{mid}}) \approx 78 \). The ratchet effect: after one complete viscosity cycle, net enstrophy increases by a factor of ~1.23. The ratchet criterion \( J \) of the diagnostic reads this asymmetry.

6. Forward and Backward Automatic Differentiation

6.1 Forward AD: Sensitivity of Future to Present

Forward AD sensitivities of \( H(T) \) to model parameters (6-mode shell model, \( A = 0.9 \)).

Parameter	\( \partial H(T)/\partial p \)	Relative Sensitivity
\( A \) (amplitude)	+12.4	1.00 (reference)
\( \nu \) (viscosity)	−8.7	0.70
\( \lambda \) (linear stretch)	+5.2	0.42
\( \lambda_2 \) (quad stretch)	+15.8	1.27
\( \sigma \) (vel-vort coupling)	+1.3	0.10
\( \beta_0 \) (cascade)	+2.1	0.17

The quadratic stretching coefficient \( \lambda_2 \) has the highest sensitivity — it is the dominant driver of enstrophy growth in the shell model.

6.2 Backward AD (Adjoint): Influence of Past on Present

Backward-mode AD solves the adjoint equation:

\[ \frac{d\lambda}{dt} = -J(\theta(t))^T \lambda, \qquad \lambda(T) = \frac{\partial H(T)}{\partial \theta(T)}. \tag{28} \]

6.3 Vorticity vs Velocity Attribution

Component-wise backward attribution at \( t = 0 \) (8-mode shell model).

Component	\( |\partial H(T)/\partial \theta_i(0)| \)	Relative	Type
\( \omega_1 \) (\( k=1 \))	8.42	13.2	Vorticity
\( \omega_2 \) (\( k=2 \))	6.17	9.7	Vorticity
\( \omega_3 \) (\( k=3 \))	5.83	9.1	Vorticity
\( \omega_4 \) (\( k=4 \))	5.11	8.0	Vorticity
\( u_1 \) (\( k=1 \))	0.64	1.0	Velocity
\( u_2 \) (\( k=2 \))	0.58	0.9	Velocity
\( u_3 \) (\( k=3 \))	0.47	0.7	Velocity
\( u_4 \) (\( k=4 \))	0.31	0.5	Velocity

Vorticity components carry 8–13× more attribution than velocity components, in this shell model.

7. Tables

The tables below record what the shell models did under the stated protocol (\( \nu = 0.005 \), \( \lambda_2 = 5.0 \), \( T = 50{,}000 \) steps, excursion multiple 100). \( A_T \) is the finite-window excursion amplitude; \( A_{T,P} \) is the amplitude at which the shrinking-doubling-time rule first fires; Score is the in-sample classification count of that rule.

7.1 Mode Table

Shell-model runs by mode count. \( N \) = total modes, \( K \) = wavenumbers.

Model	\( N \)	\( K \)	\( k_{\max} \)	\( \nu k_{\max}^2 \)	\( A_T \)	\( A_{T,P} \)	\( A_{T,P}/A_T \)	Score
6-mode	6	3	3	0.045	1.136	1.020	89.8%	17/20
8-mode	8	4	4	0.080	0.290	0.277	95.5%	16/16
10-mode	10	5	5	0.125	0.302	0.290	96.0%	13/14
12-mode	12	6	6	0.180	0.328	0.307	93.6%	14/14
16-mode	16	8	8	0.320	0.347	0.328	94.6%	14/14
20-mode	20	10	10	0.500	0.347	0.328	94.6%	14/14
24-mode	24	12	12	0.720	0.347	0.328	94.6%	14/14

7.2 Notes by Model

6-Mode Shell Model. \( K = 3 \), \( k_{\max} = 3 \). \( A_T = 1.136 \), much higher than the other models. Classification: 17/20 = 85%, with 3 false negatives near \( A_T \).

8-Mode Shell Model. \( K = 4 \), \( k_{\max} = 4 \). \( A_T = 0.290 \); the fourth wavenumber adds cascade terms. Classification: 16/16.

10-Mode and 12-Mode Shell Models. \( A_T = 0.302 \) and \( 0.328 \). Classification 13/14 and 14/14.

16, 20, and 24-Mode Shell Models. Diffusion at \( k_{\max} = 12 \) (24-mode) is \( 0.72 \), sixteen times the 6-mode value. These values are properties of the shell model under this window and threshold; the amplification factors reported in the same runs vary between mode counts, and no scaling exponent is fitted.

8. Discussion

8.1 Relation to the Unified Adaptation Theorem

The diagnostic \( H \), its self-adapting weights \( H' \) and its confidence tracker \( H'' \) are composed according to the update rules of the Unified Adaptation Theorem [21]. That composition governs how the diagnostic's own weights settle; it says nothing about the dynamics of the shell model or of the fluid equations.

8.2 Possible Uses of the Diagnostic Levels

1. Runaway-growth flagging. The shrinking-doubling-time rule (\( L_{15} \)) is cheap to evaluate along any computed trajectory; its classification counts above are in-sample.
2. Adaptive refinement. Stretching efficiency (\( L_{21} \)) could be used to focus resolution in a solver; this has not been tested.
3. Subgrid indicators. The relaxation deficit (\( L_{18} \)) is a physically motivated ratio that could be read by a closure model; this has not been tested.

8.3 Implementation in Simplex

The shell model and the diagnostic are implemented in Simplex [15], which provides native dual-number automatic differentiation: a dual number type with exact derivative propagation, forward-mode AD in a single pass, and compilation to LLVM.

9. Limitations

1. Shell model, not Navier–Stokes. The 6–24-mode systems are hand-constructed cubic ODE systems on scalar wavenumbers. They are not Galerkin truncations of the equations and do not conserve energy in the inviscid limit. No conclusion about the 3D Navier–Stokes equations follows from their behaviour.
2. The 3D kernel injects energy. The real-coefficient kernel of Section 2.7 omits the factor \( -i \) in the Fourier nonlinearity. Every 3D excursion amplitude in this record was computed with it. The corrected complex-coefficient solver is described in the companion record on energy conservation.
3. \( A_T \) is a protocol quantity. It depends on the window length (\( A_T \propto T^{-0.89} \)), on the excursion multiple, and on the solver. It is not a regularity threshold and does not converge to one.
4. \( H \) is fitted in-sample. Its weights are learned from the outcomes it reports; the classification scores are not out-of-sample accuracies.
5. Reproducibility of the shell-model tables. The equations and parameters printed here do not by themselves reproduce the tabulated excursion amplitudes; the published Simplex sources are the definition of the protocol.
6. Specific parameters and initial conditions. Results are at \( \nu = 0.005 \), \( \lambda_2 = 5.0 \), sinusoidal initial data, \( T = 50{,}000 \) steps (some runs 100,000), in 64-bit arithmetic.

10. Conclusion

This record documents a cubic shell model, a multi-level trajectory diagnostic \( H \), and a first real-coefficient Fourier–Galerkin kernel for 3D Navier–Stokes, together with the tables of what each did under the stated protocol. What stands is the doubling-time lemma, the mode-enumeration, triad-precomputation and Leray-projection structure of the kernel, which is the basis of the corrected energy-conserving solver, and the tables as a record; no conclusion about the 3D Navier–Stokes equations follows from this paper.

11. Reproducibility

All shell-model experiments are written in the Simplex programming language [15] and compiled to LLVM IR via the sxc compiler. The only dependencies are clang (for linking) and a C runtime providing sin, cos, sqrt, exp, ln.

To reproduce any experiment:

git clone https://github.com/senuamedia/lab
cd lab
./build.sh
./build/sxc <experiment>.sx -o /tmp/out.ll
clang /tmp/out.ll runtime/standalone_runtime.c \
  -o /tmp/out -lm -lssl -lcrypto
/tmp/out

The complete experimental suite, source code, and lab documentation are available at lab.senuamedia.com and GitHub.

12. References

1. C. L. Fefferman, “Existence and smoothness of the Navier–Stokes equation,” Clay Mathematics Institute Millennium Prize Problems, 2000.
2. J. Leray, “Sur le mouvement d'un liquide visqueux emplissant l'espace,” Acta Math., vol. 63, pp. 193–248, 1934.
3. O. A. Ladyzhenskaya, The Mathematical Theory of Viscous Incompressible Flow, Gordon and Breach, 2nd ed., 1969.
4. L. Caffarelli, R. Kohn, and L. Nirenberg, “Partial regularity of suitable weak solutions of the Navier–Stokes equations,” Comm. Pure Appl. Math., vol. 35, pp. 771–831, 1982.
5. J. T. Beale, T. Kato, and A. Majda, “Remarks on the breakdown of smooth solutions for the 3-D Euler equations,” Comm. Math. Phys., vol. 94, pp. 61–66, 1984.
6. T. Tao, “Finite time blowup for an averaged three-dimensional Navier–Stokes equation,” J. Amer. Math. Soc., vol. 29, pp. 601–674, 2016.
7. T. Tao, “Searching for singularities in the Navier–Stokes equations,” Nature Reviews Physics, vol. 1, pp. 418–419, 2019.
8. R. Temam, Navier–Stokes Equations: Theory and Numerical Analysis, North-Holland, 1984.
9. J. C. Robinson, Infinite-Dimensional Dynamical Systems, Cambridge University Press, 2001.
10. A. N. Kolmogorov, “The local structure of turbulence in incompressible viscous fluid for very large Reynolds numbers,” Dokl. Akad. Nauk SSSR, vol. 30, pp. 299–303, 1941.
11. S. B. Pope, Turbulent Flows, Cambridge University Press, 2000.
12. P. Sagaut, Large Eddy Simulation for Incompressible Flows, Springer, 3rd ed., 2006.
13. S. Succi, The Lattice Boltzmann Equation for Fluid Dynamics and Beyond, Oxford University Press, 2001.
14. W. T. Ashurst, A. R. Kerstein, R. M. Kerr, and C. H. Gibson, “Alignment of vorticity and scalar gradient with strain rate in simulated Navier–Stokes turbulence,” Phys. Fluids, vol. 30, pp. 2343–2353, 1987.
15. R. Higgins, “Simplex: A programming language with native dual-number automatic differentiation,” github.com/senuamedia/lab, 2026.
16. C. R. Doering, “The 3D Navier–Stokes problem,” Ann. Rev. Fluid Mech., vol. 41, pp. 109–128, 2009.
17. C. Foias, O. Manley, R. Rosa, and R. Temam, Navier–Stokes Equations and Turbulence, Cambridge University Press, 2001.
18. L. Escauriaza, G. A. Seregin, and V. Šverák, “\( L_{3,\infty} \)-solutions of Navier–Stokes equations and backward uniqueness,” Russian Math. Surveys, vol. 58, no. 2, pp. 211–250, 2003.
19. G. Prodi, “Un teorema di unicità per le equazioni di Navier–Stokes,” Ann. Mat. Pura Appl., vol. 48, pp. 173–182, 1959.
20. J. Serrin, “On the interior regularity of weak solutions of the Navier–Stokes equations,” Arch. Rational Mech. Anal., vol. 9, pp. 187–195, 1962.
21. R. Higgins, “Unified Adaptation Theorem: Convergence of Composed Adaptive Systems via Interaction Matrices and Higher-Order Convergence Diagnostics,” 2026.
22. W. Lohmiller and J.-J. E. Slotine, “On contraction analysis for non-linear systems,” Automatica, vol. 34, no. 6, pp. 683–696, 1998.
23. E. Hopf, “Über die Anfangswertaufgabe für die hydrodynamischen Grundgleichungen,” Math. Nachr., vol. 4, pp. 213–231, 1951.
24. R. Temam, Navier–Stokes Equations: Theory and Numerical Analysis, North-Holland, 1977.