\subsection{Long-Term Secular Stability, Planetary KAM Torus Preservation, and Historical Celestial Encounters}
\label{sec:long_term_stability}

The ultimate benchmark of any geometric symplectic integrator in celestial mechanics lies in its capacity to reconcile two antithetical requirements: preserving the global, secular phase-space topology of planetary systems over astronomical timescales (decades to millions of orbits) while simultaneously handling violent, localized close encounters without energy blowup or artificial dissipation.
In this section, we present the empirical validation of our normal-cone symplectic architecture across two complementary benchmark series:
\begin{enumerate}[label=(\roman*)]
    \item \textbf{Historical Encounter Benchmarks (\S\ref{sub:historical_encounters_validation}): Canonical Few-Body Models of Documented Historical Encounters}, evaluating the solver's adaptive behavior during deep gravitational flybys, sub-Roche tidal disruptions, and atlas chart transitions;
    \item \textbf{Secular Stability Benchmarks (\S\ref{sub:secular_stability_laplace}): Long-Term Secular Invariant Preservation and Frequency Analysis}, proving the absence of linear energy drift, the boundedness of the resonant Jovian Laplace angle, and the preservation of KAM invariant tori via Laskar's Numerical Analysis of Fundamental Frequencies (\textit{NAFF}) \cite{laskar1993frequency}.
\end{enumerate}

% -----------------------------------------------------------------------------
% SUBSECTION 7.2.1 : HISTORICAL ENCOUNTERS & COVARIANT ATLAS GAUGE CONTINUITY
% -----------------------------------------------------------------------------
\subsubsection{Canonical Few-Body Models of Documented Historical Encounters and Covariant Atlas Gauge Continuity}
\label{sub:historical_encounters_validation}

To test our architecture under realistic astronomical mass and distance ratios without conflating symplectic integration errors with empirical force-model uncertainties, we formulate in dual precision (\texttt{Float64} and 35-digit \texttt{Float128} under strictly identical integration parameters) three canonical, unsoftened Newtonian point-mass few-body benchmarks ($N \in \{4, 5\}$, $\varepsilon = 0, d_0 = 0$) calibrated to the characteristic encounter scales ($r_{\min}, v_\infty, a, e$) of celebrated events from historical ephemerides (intentionally omitting planetary zonal harmonics $J_2$, non-gravitational cometary outgassing, and hydrodynamical continuum dissipation so that exact conservative Hamiltonian and Lie--Bargmann invariants remain well-defined).

\medskip\noindent\textbf{Benchmark H1 : Comet D/1770 L1 (Lexell) Deep Jovian Passage and Covariant Atlas Gauge Continuity ($T_{A \to B}^{(\mathrm{PBM4})}$, Figure~\ref{fig:historical_encounters}).}
Comet Lexell underwent a historic close encounter with Jupiter in 1779, plunging deep inside the Jovian sphere of influence to a perijove of $r_{\min} = 1.4483 \times 10^{-3}\text{ AU} = 3.031\,R_J$ (well inside Io's orbital radius $a_{\mathrm{Io}} = 5.90\,R_J$). We integrate the 4-body point-mass model (Sun, Jupiter, Earth, Comet Lexell) across $t \in [0, 12.0]$ ($\approx 1.91\text{ yr}$) under a strictly constant macro-step $\Delta t = 0.08$ ($N_{\mathrm{steps}} = 150$ uniform macro-steps, maximum dyadic level $k_{\max} = 11$, zero dyadic saturation, with native solver substep observation enabled for $k_{\mathrm{sym}} \ge 4$) starting outside Jupiter's Hill sphere at $r_{0,J} = 0.7000\text{ AU} = 1{,}465\,R_J > R_{\mathrm{Hill},J} = 0.3551\text{ AU} = 743\,R_J$ in the heliocentric Jacobi chart $\mathcal{T}_{\mathrm{helio}} = [((0,2),1),3]$. As tracked in Figure~\ref{fig:historical_encounters}a, the Voronoi Hill admissibility criterion ($\gamma = 1.0$, Sections~\ref{sub:agglomerative_tree_building} and~\ref{sub:delaunay_jacobi_admissibility}) triggers two dynamic topological chart flips:
\begin{equation}
    \mathcal{T}_{\mathrm{helio}} = [((0,2),1),3] \;\underset{t_2 = 8.48}{\overset{t_1 = 2.56}{\longleftrightarrow}}\; \mathcal{T}_{\mathrm{Jovi}} = [(0,2),(1,3)],
\end{equation}
entering the Jovicentric chart at $t_{\mathrm{flip},1} = 2.56$ ($\eta_{\mathrm{Hill}} = 0.5809 < 1.0$), autonomously refining the symmetric dyadic rung along a strictly unimodal staircase $k_{\mathrm{sym}}(t) : 0 \to 1 \to 2 \to 3 \to 5 \to 6 \to 8 \to 11 \to 8 \to 5 \to 3 \to 2 \to 1$ across the deep perijove at $t_{\mathrm{peri}} = 5.4584$ ($r_{\min} = 3.031\,R_J$, resolved with $2^{11} = 2{,}048$ dyadic substeps of size $h_{\mathrm{sub}} = 3.91 \times 10^{-5}$), and flipping back to the heliocentric chart on egress at $t_{\mathrm{flip},2} = 8.48$ ($\eta_{\mathrm{Hill}} = 0.9969$, Jovicentric dwell interval $\Delta t_{\mathrm{Jovi}} = 5.92$).

At each chart transition and dyadic step change, the canonical non-commutative cotangent-lifted gauge connection established in Section~\ref{sub:pbm4_atlas_gauge_transitions} and Section~\ref{sub:canonical_gauge_connection} (Theorems~\ref{thm:pbm4_gauge_operator_exact} and~\ref{thm:pbm4_gauge_operator_sec6}):
\begin{equation}
\begin{aligned}
    T_{(A,h_A) \to (B,h_B)}^{(\mathrm{PBM4})} &= \bigl(\mathcal{Z}_4^{\mathrm{BM4},(B)}(h_B)\bigr)^{-1} \circ \mathbf{M}_{A\to B} \circ \mathcal{Z}_4^{\mathrm{BM4},(A)}(h_A), \\
    \mathcal{Z}_4^{\mathrm{BM4},(\alpha)} &:= \mathcal{Z}_4^{\mathrm{BM4}}(h_\alpha;\mathcal{T}_\alpha),
\end{aligned}
\end{equation}
maps the shadow state seamlessly across the chart boundary, where $\mathbf{M}_{A\to B} \in \operatorname{SL}_2(\mathbb{R}) \subset \operatorname{Sp}(6N,\mathbb{R})$ denotes the symplectic Jacobi coordinate transition morphism:
\begin{equation}
\begin{aligned}
    \tilde{\mathbf{z}}_B &= T_{(A,h_A) \to (B,h_B)}^{(\mathrm{PBM4})}(\tilde{\mathbf{z}}_A), \\
    \mathbf{z}_{\mathrm{phys}} &= \mathcal{Z}_4^{\mathrm{BM4},(B)}(h_B)\bigl(\tilde{\mathbf{z}}_B\bigr) = \mathbf{M}_{A\to B}\bigl(\mathcal{Z}_4^{\mathrm{BM4},(A)}(h_A)(\tilde{\mathbf{z}}_A)\bigr).
\end{aligned}
\end{equation}

As shown in Figure~\ref{fig:historical_encounters}b, the interfacial jump in the comet's specific orbital energy at $t_{\mathrm{flip},1} = 2.56$ is bounded to $|\Delta\mathcal{E}_3|/|\mathcal{E}_{3,0}| = 0.00$ in \texttt{Float64} and $4.48 \times 10^{-35}$ in \texttt{Float128} (with total 4-body Hamiltonian jump $|\Delta H_{\mathrm{flip}}|/|H_0| = 2.36 \times 10^{-35}$ in \texttt{Float128}). Over the entire encounter $t \in [0, 12.0]$, relative Hamiltonian error remains bounded to $\max_t |\Delta H(t)|/|H_0| = 5.48 \times 10^{-13}$ in \texttt{Float64} and $7.13 \times 10^{-17}$ in \texttt{Float128}: in accordance with backward error analysis (Theorem~\ref{thm:shadow_alignment} and Proposition~\ref{prop:dyadic_shadow_continuity}), \texttt{Float128} eliminates floating-point roundoff across the $k_{\max} = 11$ perijove substeps and isolates the bounded $\mathcal{O}(h_{\mathrm{sub}}^6 + \varepsilon\,h_{\mathrm{sub}}^4)$ shadow-Hamiltonian truncation oscillation at $\Delta t = 0.08$, whereas the exact quadratic Lie--Bargmann Noether invariants (Propositions~\ref{prop:strang_bch_parity_bargmann} and~\ref{prop:lie_bargmann_conservation_spectrum}) exhibit zero truncation error and reach the raw arithmetic floor $\|\Delta\mathbf{P}_{\mathrm{tot}}\|_2 \le 1.20 \times 10^{-16}$ / $5.57 \times 10^{-38}$ and $\|\Delta\mathbf{L}_{\mathrm{tot}}\|_2 \le 6.05 \times 10^{-16}$ / $7.74 \times 10^{-35}$ in \texttt{Float64} / \texttt{Float128}.

\begin{figure}[!htbp]
    \centering
    \includegraphics[width=\textwidth]{chapters/ch07_benchmarks_and_lean4/figures/fig_ch07_h1_lexell_atlas_gauge.png}
    \caption{\textbf{Benchmark H1 --- Comet D/1770 L1 (Lexell, 1770--1779): deep Jovian encounter ($3.031\,R_J$), dynamic Hill chart flips, and covariant gauge continuity ($N=4$).} \textbf{(a)}~Jovicentric distance $r_J(t)/R_J$ across $\mathcal{T}_{\mathrm{helio}} \leftrightarrow \mathcal{T}_{\mathrm{Jovi}}$ chart flips ($t = 2.56, 8.48$) and unimodal dyadic rung $k_{\mathrm{sym}}(t) \in \{0,\dots,11\}$. \textbf{(b)}~Continuous Hamiltonian error $|\Delta H(t)|/|H_0|$ under $T_{A\to B}^{(\mathrm{PBM4})}$ and Lie--Bargmann invariants $\|\Delta\mathbf{L}_{\mathrm{tot}}\|_2, \|\Delta\mathbf{P}_{\mathrm{tot}}\|_2$ (\texttt{Float64} and \texttt{Float128}).}
    \label{fig:h1_lexell_atlas_gauge}
    \label{fig:historical_encounters}
\end{figure}

\medskip\noindent\textbf{Benchmark H2 : Comet D/1993 F2 (Shoemaker--Levy 9) Sub-Roche Tidal Disruption and Synchronous Supernode Promotion (Figure~\ref{fig:h2_shoemaker_levy_roche}).}
During its July 7, 1992 perijove pass, Comet Shoemaker--Levy~9 penetrated deep inside Jupiter's fluid Roche limit ($r_{\mathrm{Roche}} = 2.45\,R_J = 1.1708 \times 10^{-3}\text{ AU}$). We model this regime via a 5-body point-mass system (Sun, Jupiter, and $K = 3$ self-gravitating progenitor fragments $\{2,3,4\}$ of mass $m_f = 10^{-11} M_\odot$ separated by $\delta r_0 = 2.0 \times 10^{-5}\text{ AU}$) starting at $r_{0,J} = 0.0120\text{ AU} = 25.11\,R_J$ across $N_{\mathrm{steps}} = 160$ macro-steps ($\Delta t = 1.8 \times 10^{-4}$, $t_{\max} = 0.0288 > t_{\mathrm{peri}} = 0.02178$, maximum dyadic level $k_{\max} = 4$, zero saturation) in both \texttt{Float64} and \texttt{Float128}. Figure~\ref{fig:h2_shoemaker_levy_roche}a tracks both the Jovicentric distance $r_J(t)/R_J$ and the dimensionless cubic tidal force ratio:
\begin{equation}
    \Lambda_S(t) = \frac{\|\nabla W_{\mathrm{apex}}^{\mathrm{screen}}\|}{\|\nabla H_{\mathrm{leaf}}\|} \approx \frac{M_J}{3 M_S}\Bigl(\frac{\bar{\rho}_S(t)}{r_J(t)}\Bigr)^3 = \frac{1}{3}\,\eta_S(t)^3,
\end{equation}
where $\eta_S = \bar{\rho}_S / R_{\mathrm{Hill}}(S)$ denotes the linear Hill ratio of Section~\ref{sub:synchronous_supernode_anti_zeno}.

Initially outside the Roche zone ($\Lambda_S(0) = 0.0267 < \Lambda_{\mathrm{crit}} = 0.25$), the progenitor cluster evolves as standard binary leaves until $\Lambda_S(t)$ crosses the Roche disruption threshold $\Lambda_{\mathrm{crit}} = 0.25$, dynamically promoting $\{2,3,4\}$ into a synchronous supernode across $84/160$ steps (canonical sub-Roche tidal Channel~A, Section~\ref{sub:synchronous_supernode_anti_zeno}). At the sub-Roche perijove $t_{\mathrm{peri}} = 0.02178$, the cluster reaches $r_{\min} = 6.3839 \times 10^{-4}\text{ AU} = 1.336\,R_J < 2.45\,R_J$ with $\Lambda_S^{\max} = 16{,}690$ ($\eta_S^{\max} = 36.86$, unimodal symmetric dyadic rung pyramid $k_{\mathrm{sym}} : 0 \to 1 \to 2 \to 3 \to 4 \to 3 \to 2 \to 1 \to 0$). As shown in Figure~\ref{fig:h2_shoemaker_levy_roche}b--c, differential Jovian tidal shear stretches the fragment chain into the ``string-of-pearls'' configuration with span ratio $S(t)/S_0 = \|\mathbf{q}_4(t) - \mathbf{q}_2(t)\|_2 / S_0$ peaking at $4.61\times$ just after perijove and reaching $3.82\times$ at $t_{\max} = 0.0288$.

In Figure~\ref{fig:h2_shoemaker_levy_roche}c, the \texttt{Float128} energy curve starts at $\sim 10^{-21}$ during the outer approach (where \texttt{Float64} sits at the $\sim 10^{-15}$ IEEE~754 roundoff floor) and rises smoothly to coincide with \texttt{Float64} at perijove ($\max_t |\Delta H(t)|/|H_0| = 3.99 \times 10^{-15}$ in \texttt{Float64} and $4.40 \times 10^{-15}$ in \texttt{Float128}), reflecting the reversible $\mathcal{O}(h_{\mathrm{sub}}^6)$ shadow-Hamiltonian truncation oscillation under extreme sub-Roche tidal shear ($\Lambda_S^{\max} = 16{,}690$ at $h_{\mathrm{sub}} = \Delta t / 2^4 = 1.125 \times 10^{-5}$), whereas the exact Lie--Bargmann Noether invariants remain unaffected by tidal truncation at $\|\Delta\mathbf{L}_{\mathrm{tot}}\|_2 \le 5.64 \times 10^{-18}$ / $1.10 \times 10^{-35}$ and $\|\Delta\mathbf{P}_{\mathrm{tot}}\|_2 \le 2.89 \times 10^{-18}$ / $1.19 \times 10^{-37}$ in \texttt{Float64} / \texttt{Float128}.

\begin{figure}[!htbp]
    \centering
    \includegraphics[width=\textwidth]{chapters/ch07_benchmarks_and_lean4/figures/fig_ch07_h2_shoemaker_levy_roche.png}
    \caption{\textbf{Benchmark H2 --- Comet D/1993 F2 (Shoemaker--Levy 9, July 1992): sub-Roche perijove ($1.336\,R_J < 2.45\,R_J$), dynamic tidal promotion ($\Lambda_S \ge 0.25$), and string-of-pearls elongation ($N=5$).} \textbf{(a)}~Jovicentric distance $r_J(t)/R_J$ and cubic tidal force ratio $\Lambda_S(t) = \frac{1}{3}\eta_S(t)^3$ ($\Lambda_S^{\max} = 16{,}690$). \textbf{(b)}~Fragment chain span ratio $S(t)/S_0$ ($4.61\times$ peak) and dyadic rungs $k_{\mathrm{sym}}(t) \in \{0,\dots,4\}$. \textbf{(c)}~Shadow-Hamiltonian oscillation $|\Delta H(t)|/|H_0|$ and Lie--Bargmann invariants (\texttt{Float64} and \texttt{Float128}).}
    \label{fig:h2_shoemaker_levy_roche}
\end{figure}
\medskip\noindent\textbf{Benchmark H3 : Asteroid (99942) Apophis 2029 Close Earth Flyby, Aten $\to$ Apollo Deflection, and Full-Flyby Palindromic Reversibility (Figure~\ref{fig:h3_apophis_flyby_palindrome}).}
On April 13, 2029, asteroid (99942) Apophis will pass inside Earth's geostationary satellite belt ($r_{\mathrm{GEO}} = 6.61\,R_\oplus$). We integrate the 4-body Sun--Earth--Moon--Apophis point-mass model across $N_{\mathrm{steps}} = 80$ macro-steps ($\Delta t = 2.5 \times 10^{-4}$, $t_{\max} = 0.0200 > t_{\mathrm{peri}} = 0.00950$) in both \texttt{Float64} and \texttt{Float128} (Figure~\ref{fig:h3_apophis_flyby_palindrome}). Approaching on its pre-flyby Aten-class orbit ($a_{\mathrm{pre}} = 0.9224\text{ AU} < 1.0\text{ AU}$, $e_{\mathrm{pre}} = 0.219$), Apophis grazes perigee at $r_{\min} = 2.5480 \times 10^{-4}\text{ AU} = 5.976\,R_\oplus \approx 38{,}115\text{ km}$ at $t_{\mathrm{peri}} = 0.00950$ (Figure~\ref{fig:h3_apophis_flyby_palindrome}a) and is deflected outward by Earth's trailing-side gravity assist into a post-flyby Apollo-class orbit ($a_{\mathrm{post}} = 1.1163\text{ AU} > 1.0\text{ AU}$, $e_{\mathrm{post}} = 0.217$), producing a physical semi-major axis jump of $\Delta a = a_{\mathrm{post}} - a_{\mathrm{pre}} = +0.1939\text{ AU} > 0$ (Figure~\ref{fig:h3_apophis_flyby_palindrome}b).

Across the bilateral pericenter guard $|t_{\mathrm{peri}}(t)| \le 2\,\Delta t_{\mathrm{macro}}$ (Definitions~\ref{def:conormal_timescales} and~\ref{def:conormal_timescale_triad}), the symmetric dyadic selector (Propositions~\ref{prop:skeel_gear_boundary_sync} and~\ref{prop:macroblock_palindromic_symplecticity}):
\[
    k_{\mathrm{sym}}(\mathbf{z}_n, \Delta t) = \max\Bigl(k(\mathbf{z}_n),\, k\bigl(\Phi_{\Delta t}^{(0)}(\mathbf{z}_n)\bigr)\Bigr), \qquad
    k_{\mathrm{sym}} : 0 \to 1 \to 2 \to 3 \to 2 \to 1 \to 0,
\]
climbs along the unimodal pyramid peaking at $k_{\max} = 3$ at perigee (step $m = 38$, Figure~\ref{fig:h3_apophis_flyby_palindrome}c), bounding the forward-flyby Hamiltonian excursion to $\max_t |\Delta H(t)|/|H_0| = 3.33 \times 10^{-14}$ in \texttt{Float64} and $2.95 \times 10^{-22}$ ($\mathcal{O}(h_{\mathrm{sub}}^6)$ shadow oscillation) in \texttt{Float128} (with Noether conservation $\|\Delta\mathbf{L}_{\mathrm{tot}}\|_2 \le 5.12 \times 10^{-20}$ / $6.42 \times 10^{-39}$).

Executing a full-flyby forward--backward reversibility test $\Phi_{-\Delta t}^{N_h} \circ \Phi_{+\Delta t}^{N_h}$ across perigee ($N_h \in \{20, 40, 60, 80\}$ steps, Figure~\ref{fig:h3_apophis_flyby_palindrome}c--d) confirms exact point-by-point palindromic dyadic sequence matching $k_{\mathrm{fwd}}(m) \equiv k_{\mathrm{bwd}}(N-1-m)$ in both \texttt{Float64} and \texttt{Float128} and annihilates the conservative shadow-Hamiltonian oscillation upon returning to $t = 0$, yielding 80-step return position error $\|\Delta\mathbf{q}_{\mathrm{FB}}\|_2 = 1.72 \times 10^{-16}\text{ AU}$ in \texttt{Float64} and $3.28 \times 10^{-33}\text{ AU}$ in \texttt{Float128}, and return energy error $|\Delta H_{\mathrm{FB}}|/|H_0| = 0.0$ across all four horizons in \texttt{Float64} and $8.46 \times 10^{-33}$ in \texttt{Float128}.
\pythonrefright{py_bench_gen04}{{benchmarks/run, figures/plot}_generator_04_*.py}{Float64 / Float128}

\begin{figure}[!htbp]
    \centering
    \includegraphics[width=\textwidth]{chapters/ch07_benchmarks_and_lean4/figures/fig_ch07_h3_apophis_flyby_palindrome.png}
    \caption{\textbf{Benchmark H3 --- Asteroid (99942) Apophis (April 13, 2029): sub-geostationary perigee ($5.976\,R_\oplus$), Aten $\to$ Apollo deflection, and full-flyby dyadic palindrome ($N=4$).} \textbf{(a)}~Geocentric hyperbolic trajectory inside the geostationary belt ($r_{\mathrm{GEO}} = 6.61\,R_\oplus$). \textbf{(b)}~Heliocentric osculating $a(t)$ and $e(t)$ ($\Delta a = +0.1939\text{ AU}$). \textbf{(c)}~Bilateral perigee horizon $|t_{\mathrm{peri}}|/\Delta t_{\mathrm{macro}} \le 2$ and palindromic dyadic rungs $k_{\mathrm{fwd}}(m) \equiv k_{\mathrm{bwd}}(N-1-m)$ ($k_{\max} = 3$). \textbf{(d)}~Multi-horizon round-trip reversibility error $\Phi_{-\Delta t}^{N_h} \circ \Phi_{+\Delta t}^{N_h}$ ($N_h \in \{20, 40, 60, 80\}$).}
    \label{fig:h3_apophis_flyby_palindrome}
\end{figure}

\begin{table}[!htbp]
    \centering
    \caption{\textbf{Certified dual-precision (\texttt{Float64} and 35-digit \texttt{Float128}) quantitative invariants and observables across documented historical astronomical encounters (\S\ref{sub:historical_encounters_validation}; $\gamma = 1.0, \varepsilon = 0, d_0 = 0$).}}
    \label{tab:campaign_II_kpi_table}
    \vspace{0.3em}
    \renewcommand{\arraystretch}{1.2}
    \resizebox{\textwidth}{!}{%
    \begin{tabular}{lcccccc}
        \toprule
        \makecell{\textbf{Historical}\\\textbf{Benchmark}} & \makecell{\textbf{Closest Approach}\\$r_{\min}$} & \makecell{\textbf{Max Rung}\\$k_{\max}$} & \makecell{\textbf{Topological / Tidal}\\\textbf{Regime}} & \makecell{\textbf{Max $|\Delta H|/|H_0|$}\\\textbf{(\texttt{Float64} / \texttt{Float128})}} & \makecell{\textbf{Noether $\|\Delta\mathbf{L}_{\mathrm{tot}}\|_2$}\\\textbf{(\texttt{Float64} / \texttt{Float128})}} & \makecell{\textbf{Certified Physical \&}\\\textbf{Algebraic Observables}} \\
        \midrule
        \makecell[l]{\textbf{H1}: Comet Lexell\\(1770--1779, $N=4$)} & \makecell{$3.031\,R_J$\\($1.4483 \times 10^{-3}\text{ AU}$)} & $11$ & \makecell{2 Hill chart flips\\$\mathcal{T}_{\mathrm{helio}} \leftrightarrow \mathcal{T}_{\mathrm{Jovi}}$\\($t = 2.56,\,8.48$)} & \makecell{$5.48 \times 10^{-13}$ /\\$7.13 \times 10^{-17}$} & \makecell{$6.05 \times 10^{-16}$ /\\$7.74 \times 10^{-35}$} & \makecell{Flip jump under $T_{A\to B}^{(\mathrm{PBM4})}$:\\$|\Delta\mathcal{E}_3|/|\mathcal{E}_{3,0}| = 0.00$ (\texttt{F64}),\\$|\Delta H_{\mathrm{flip}}|/|H_0| = 2.36 \times 10^{-35}$ (\texttt{F128})} \\
        \addlinespace[0.25em]
        \makecell[l]{\textbf{H2}: Shoemaker--Levy 9\\(July 1992, $N=5$)} & \makecell{$1.336\,R_J < 2.45\,R_J$\\($6.3839 \times 10^{-4}\text{ AU}$)} & $4$ & \makecell{Sub-Roche promotion\\$\Lambda_S(t) \ge 0.25$\\($84/160$ steps)} & \makecell{$3.99 \times 10^{-15}$ /\\$4.40 \times 10^{-15}$} & \makecell{$5.64 \times 10^{-18}$ /\\$1.10 \times 10^{-35}$} & \makecell{$\Lambda_S(0) = 0.0267 \to \Lambda_S^{\max} = 16{,}690$,\\string-of-pearls span $S/S_0$: $4.61\times$ peak\\($3.82\times$ at egress)} \\
        \addlinespace[0.25em]
        \makecell[l]{\textbf{H3}: (99942) Apophis\\(April 2029, $N=4$)} & \makecell{$5.976\,R_\oplus < 6.61\,R_\oplus$\\($2.5480 \times 10^{-4}\text{ AU}$)} & $3$ & \makecell{Bilateral perigee guard\\$|t_{\mathrm{peri}}| \le 2\,\Delta t_{\mathrm{macro}}$\\($80$ macro-steps)} & \makecell{$3.33 \times 10^{-14}$ /\\$2.95 \times 10^{-22}$\\(FB: $0.0$ / $8.46 \times 10^{-33}$)} & \makecell{$5.12 \times 10^{-20}$ /\\$6.42 \times 10^{-39}$} & \makecell{Aten ($0.9224\text{ AU}$) $\to$ Apollo ($1.1163\text{ AU}$),\\$\Delta a = +0.1939\text{ AU}$, exact dyadic palindrome,\\$\|\Delta\mathbf{q}_{\mathrm{FB}}\|_2 = 1.72 \times 10^{-16}$ / $3.28 \times 10^{-33}\text{ AU}$} \\
        \bottomrule
    \end{tabular}%
    }
\end{table}

% -----------------------------------------------------------------------------
% SUBSECTION 7.2.2 : LONG-TERM SECULAR STABILITY & LAPLACE RESONANCE
% -----------------------------------------------------------------------------
\subsubsection{Long-Term Secular Stability, 9-Generator Lie--Bargmann Metrology, Laplace 4:2:1 Resonance, and Laskar NAFF Frequency Confinement}
\label{sub:secular_stability_laplace}

Complementing the localized, high-eccentricity close encounters of Section~\ref{sub:historical_encounters_validation}, the second benchmark series evaluates long-term secular stability under canonical Voronoi chart admissibility ($\gamma = 1.0$, $n_{\mathrm{flips}} = 0$) and active adaptive dyadic substepping ($0$ saturation events) across all five symplectic splitting and processed integrators ($\mathcal{S}_2$ with $s_{\mathrm{eval}} = 1$, $\mathcal{Y}_4$ with $s_{\mathrm{eval}} = 3$, bare $\mathrm{BM}_4^{(7)}$ with $s_{\mathrm{eval}} = 6$, homologically processed $\mathrm{PBM}_4(\mathcal{Z}_4^{\mathrm{BM4}})$ with $s_{\mathrm{eval}} = 6$, and bare $\mathrm{BM}_6^{(11)}$ with $s_{\mathrm{eval}} = 10$ force evaluations per substep).

Whereas the four bare splitting schemes ($\mathcal{S}_2, \mathcal{Y}_4, \mathrm{BM}_4^{(7)}, \mathrm{BM}_6^{(11)}$) evolve on their native bare symplectic manifolds, the processed integrator $\mathrm{PBM}_4(\mathcal{Z}_4^{\mathrm{BM4}})$ evolves natively in the shadow manifold with non-destructive physical coordinate readout at observation epochs and automatic covariant gauge transport across dyadic step-size junctions $h_A \to h_B$:
\[
    \mathbf{z}_{\mathrm{phys}}(t_k) = \mathcal{Z}_4^{\mathrm{BM4}}(h_k)(\tilde{\mathbf{z}}_k), \qquad
    T_{(A,h_A)\to(B,h_B)}^{(\mathrm{PBM4})} = \bigl(\mathcal{Z}_4^{\mathrm{BM4},(B)}(h_B)\bigr)^{-1} \circ \mathbf{M}_{A\to B} \circ \mathcal{Z}_4^{\mathrm{BM4},(A)}(h_A).
\]
Applying $T_{(A,h_A)\to(B,h_B)}^{(\mathrm{PBM4})}$ exclusively to $\mathrm{PBM}_4(\mathcal{Z}_4^{\mathrm{BM4}})$ is essential: injecting $\mathcal{Z}_4^{\mathrm{BM4}}$ gauge kicks into a bare integrator such as $\mathrm{BM}_6^{(11)}$ across dyadic junctions would perturb its bare manifold by uncompensated $\mathcal{O}(\varepsilon h^4)$ impulses. Both conservative point-mass benchmarks ($N \in \{4, 5\}$, $\varepsilon = 0, d_0 = 0$) are executed over their full secular horizon in IEEE-754 \texttt{Float64} and paired with equal-parameter 35-digit \texttt{Float128} integrations to separate bounded $\mathcal{O}(h^p)$ shadow-Hamiltonian truncation oscillations from $\mathcal{O}(\epsilon_{\mathrm{mach}}\sqrt{s_{\mathrm{eval}}\,N_{\mathrm{steps}}})$ floating-point roundoff diffusion on the 9 Galilean Lie--Bargmann Noether invariants (Figures~\ref{fig:s1_jovian_laplace}--\ref{fig:s2_solar_system_brouwer_naff} and Table~\ref{tab:campaign_III_kpi_table}).

\medskip\noindent\textbf{Protocol S1 : Resonant Jovian $4\mathbin{:}2\mathbin{:}1$ Laplace Satellite System, Methodological Libration Ablation, and Dual-Precision Lie--Bargmann Metrology ($N=4$, Figure~\ref{fig:s1_jovian_laplace}).}
We model the four-body Jovian Galilean system as a conservative Newtonian point-mass hierarchy comprising Jupiter ($m_0 = 1.0$) and the three inner moons Io ($m_1 = 4.70 \times 10^{-5}, a_1 = 1.0$), Europa ($m_2 = 2.53 \times 10^{-5}, a_2 = 1.5874$), and Ganymede ($m_3 = 7.80 \times 10^{-5}, a_3 = 2.5198$). To avoid exciting artificial free-eccentricity transients, the system is initialized at the forced resonant Laplace apsidal equilibrium ($e_{1,\mathrm{forced}} = 0.0041, \varpi_1 = 0$; $e_{2,\mathrm{forced}} = 0.0094, \varpi_2 = \pi$; $e_{3,\mathrm{forced}} = 0.0006, \varpi_3 = 0$) in a single Jovicentric Jacobi chart ($n_{\mathrm{flips}} = 0$) and integrated over $200$ Io orbits ($T_{\max} = 1{,}256.75\text{ rad}$, macro-step $\Delta t = 0.25 \approx T_{\mathrm{Io}}/25.13$, $N_{\mathrm{steps}} = 5{,}027$, observation stride $5$, where the symmetric dyadic scheduler selects $k_{\mathrm{sym}} \equiv 3$, i.e., $h_{\mathrm{sub}} = \Delta t / 2^3 = 0.03125$ with $0$ saturation events), tracking three-body resonant locking via the osculating mean-longitude critical angle:
\begin{equation}
    \label{eq:laplace_angle}
    \theta_{\mathrm{Lap}}(t) = \bigl(\lambda_1(t) - 3\lambda_2(t) + 2\lambda_3(t)\bigr) \bmod 360^\circ \approx 180^\circ, \qquad \lambda_i(t) = \varpi_i(t) + M_i(t).
\end{equation}
The four panels of Figure~\ref{fig:s1_jovian_laplace} establish three distinct physical and numerical properties:
\begin{itemize}
    \item \textit{Shadow-Manifold Telescoping, \texttt{Float64} vs.\ \texttt{Float128} Time Series, and Resonant Phase Hierarchy (Panels~\ref{fig:s1_jovian_laplace}a--b)}: Under adaptive dyadic substepping ($k_{\mathrm{sym}} = 3$, $h_{\mathrm{sub}} = 0.03125$, $40{,}216$ substeps over $200$ Io orbits), the $2^3 = 8\times$ step reduction prevents the second-order Strang splitting $\mathcal{S}_2$ ($\max|\Delta H|/|H_0| = 3.98 \times 10^{-7}$, $\mathrm{RMS} = 2.47 \times 10^{-7}$) from breaking the resonance into full $360^\circ$ circulation; however, because its $\mathcal{O}(h_{\mathrm{sub}}^2)$ shadow Hamiltonian perturbation ($h_{\mathrm{sub}}^2 \approx 9.77 \times 10^{-4}$) shifts the modified Keplerian mean motions $n_i = \partial\widetilde{H}/\partial\Lambda_i$ slightly off exact $4\mathbin{:}2\mathbin{:}1$ commensurability, $\mathcal{S}_2$ still exhibits a $+7^\circ$ secular phase drift from $180^\circ$ to $187^\circ$ ($\theta_{\mathrm{Lap}} = 183.41^\circ \pm 3.53^\circ$, Figure~\ref{fig:s1_jovian_laplace}b, top).

    By contrast, for the fourth-order Yoshida scheme $\mathcal{Y}_4$ ($\max|\Delta H|/|H_0| = 1.98 \times 10^{-10}$), the $8\times$ dyadic reduction compresses its $\mathcal{O}(h_{\mathrm{sub}}^4)$ shadow frequency error by $8^4 = 4{,}096\times$, completely eliminating the secular phase drift observed at unsubstepped $k=0$ and locking $\mathcal{Y}_4$ onto the exact resonant libration curve $\theta_{\mathrm{Lap}}(t) = 179.85^\circ \pm 0.80^\circ$ alongside $\mathrm{BM}_4^{(7)}$ ($6.96 \times 10^{-13}$), $\mathrm{PBM}_4(\mathcal{Z}_4^{\mathrm{BM4}})$, and $\mathrm{BM}_6^{(11)}$. Plotting the logarithmic phase residual $|\theta_{\mathrm{Lap}}(t) - \theta_{\mathrm{Lap}}^{\mathrm{BM6}}(t)|$ in the bottom sub-panel of Figure~\ref{fig:s1_jovian_laplace}b resolves the underlying six-order-of-magnitude accuracy hierarchy across the integrators: the maximum resonant angle error drops by $1{,}445\times$ from $\mathcal{S}_2$ ($6.13^\circ$) to $\mathcal{Y}_4$ ($4.24 \times 10^{-3\,\circ}$) and by \textbf{$730{,}000\times$} down to $8.78 \times 10^{-6\,\circ}$ ($\mathrm{BM}_4^{(7)}$) and $8.40 \times 10^{-6\,\circ}$ ($\mathrm{PBM}_4$), directly demonstrating KAM symplectic phase conjugacy ($\mathrm{PBM}_4 \equiv \mathrm{BM}_4^{(7)}$) on the resonant torus.

    In energy (Figure~\ref{fig:s1_jovian_laplace}a, plotting both \texttt{Float64} solid and 35-digit \texttt{Float128} dashed time series over the full $200$-orbit horizon), the \texttt{Float64} and \texttt{Float128} trajectories of $\mathcal{S}_2$, $\mathcal{Y}_4$, and $\mathrm{BM}_4^{(7)}$ superimpose ($3.98 \times 10^{-7}$, $1.98 \times 10^{-10}$, and $6.96 \times 10^{-13}$ vs.\ $6.67 \times 10^{-13}$), whereas in \texttt{Float64} $\mathrm{PBM}_4$ and $\mathrm{BM}_6^{(11)}$ encounter the IEEE-754 Brouwer roundoff floor ($\max = 1.20 \times 10^{-13}$, $\mathrm{RMS} = 5.77 \times 10^{-14}$, $5.80\times$ max / $7.47\times$ RMS gain for $\mathrm{PBM}_4$, and $\max = 8.14 \times 10^{-14}$, $\mathrm{RMS} = 2.88 \times 10^{-14}$ for $\mathrm{BM}_6^{(11)}$); in 35-digit \texttt{Float128}, roundoff is eliminated and the dashed curves stay on their true periodic shadow-Hamiltonian bands across all $200$ Io orbits ($\mathrm{PBM}_4$: $\max = 3.80 \times 10^{-14}$, $\mathrm{RMS} = 1.25 \times 10^{-14}$, i.e., a \textbf{$17.56\times$ max / $32.62\times$ RMS} shadow gain over $\mathrm{BM}_4^{(7)}$; $\mathrm{BM}_6^{(11)}$: $\max = 1.90 \times 10^{-18}$, $\mathrm{RMS} = 1.02 \times 10^{-18}$).
    \item \textit{Methodological Ablation of Readout Coordinates and Initial Conditions (Panel~\ref{fig:s1_jovian_laplace}c)}: The two-tier layout of Figure~\ref{fig:s1_jovian_laplace}c under $\mathrm{PBM}_4(\mathcal{Z}_4^{\mathrm{BM4}})$ separates two classical metrological artifacts when measuring Equation~\eqref{eq:laplace_angle}: (i) in the top zoomed sub-panel ($[174.8^\circ, 185.2^\circ]$), substituting true polar angles $\vartheta_i = \operatorname{atan2}(y_i, x_i)$ for osculating mean longitudes $\lambda_i = \varpi_i + M_i$ superimposes the short-period equation-of-the-center modulation $\vartheta_i - \lambda_i = 2e_i\sin M_i + \mathcal{O}(e_i^2)$, inflating the apparent half-amplitude $5.1\times$ from $\pm 0.80^\circ$ to $\pm 4.05^\circ$; and (ii) in the bottom sub-panel, initializing the moons on unforced circular orbits ($e_i(0) = 0$) instead of the forced Laplace equilibrium ($10^3 e_i(t) \approx (4.1, 9.4, 0.6)$ plotted on the right axis) displaces the trajectory away from the resonant fixed point and excites a $39.7\times$ larger free-eccentricity libration of $\pm 31.78^\circ$.
    \item \textit{Paired Dual-Precision (\texttt{Float64} vs.\ \texttt{Float128}) Energy Truncation and 9-Generator Lie--Bargmann Metrology (Panel~\ref{fig:s1_jovian_laplace}d)}: Comparing all five integrators over the full $200$-orbit equal-parameter dual-precision benchmark ($N_{\mathrm{steps}} = 5{,}027$ at $\Delta t = 0.25, k_{\mathrm{sym}} = 3$) in \texttt{Float64} (solid columns) and 35-digit \texttt{Float128} (hatched columns) directly illustrates the structural distinction between shadow-\allowbreak Hamiltonian truncation and exact Noether conservation. For $\mathcal{S}_2$, $\mathcal{Y}_4$, and $\mathrm{BM}_4^{(7)}$, the purple energy columns $|\Delta H|/|H_0|$ have identical height in \texttt{Float64} and \texttt{Float128} ($\mathcal{S}_2$: $3.98 \times 10^{-7}$; $\mathcal{Y}_4$: $1.98 \times 10^{-10}$; $\mathrm{BM}_4^{(7)}$: $6.96 \times 10^{-13}$ vs.\ $6.67 \times 10^{-13}$), while $\mathrm{PBM}_4$ drops from $1.20 \times 10^{-13}$ to $3.80 \times 10^{-14}$ and $\mathrm{BM}_6^{(11)}$ drops from its \texttt{Float64} roundoff ceiling ($8.14 \times 10^{-14}$) down to its true $\mathcal{O}(h_{\mathrm{sub}}^6)$ shadow level ($1.90 \times 10^{-18}$) in \texttt{Float128}.

    Conversely, because all leaf and apex sub-flows depend exclusively on Euclidean pairwise distances $r_{ij} = \|\mathbf{q}_i - \mathbf{q}_j\|_2$ (Propositions~\ref{prop:strang_bch_parity_bargmann} and~\ref{prop:lie_bargmann_conservation_spectrum}), the 9 Galilean Lie--Bargmann Noether generators $(\mathbf{P}_{\mathrm{tot}}, \mathbf{L}_{\mathrm{tot}}, \mathbf{G}_{\mathrm{Gal}})$ have identically zero truncation error at any step size $h > 0$: their residual drift drops by $\sim 18\text{--}20$ orders of magnitude from \texttt{Float64} (solid columns: $\|\Delta\mathbf{P}\|_2/M_{\mathrm{tot}} \le 7.16 \times 10^{-18}$, $\|\Delta\mathbf{L}\|_2/\|\mathbf{L}_0\|_2 \le 5.28 \times 10^{-14}$, $\|\Delta\mathbf{G}\|_2/M_{\mathrm{tot}} \le 5.09 \times 10^{-15}$) down to the 35-digit \texttt{Float128} arithmetic floor (hatched columns: $\|\Delta\mathbf{P}\|_2/M_{\mathrm{tot}} \in [1.55 \times 10^{-38}, 7.83 \times 10^{-38}]$, $\|\Delta\mathbf{L}\|_2/\|\mathbf{L}_0\|_2 \in [3.74 \times 10^{-32}, 4.15 \times 10^{-31}]$, and $\|\Delta\mathbf{G}\|_2/M_{\mathrm{tot}} \le 3.55 \times 10^{-34}$).
\end{itemize}
\pythonrefright{py_bench_gen05}{benchmarks/{run_S1_laplace_resonance, run_S2_brouwer_naff}.py}{Float64 / Float128}

\begin{figure}[!htbp]
    \centering
    \includegraphics[width=\textwidth]{chapters/ch07_benchmarks_and_lean4/figures/fig_ch07_s1_jovian_laplace_resonance.png}
    \caption{\textbf{Protocol S1: Resonant Jovian Galilean $4\mathbin{:}2\mathbin{:}1$ Laplace satellite system, dual-precision time series, resonant phase hierarchy, and methodological ablation ($N=4$, $200$ Io orbits, $\Delta t = 0.25, k_{\mathrm{sym}} \equiv 3$).} \textbf{(a)}~Relative Hamiltonian error $|\Delta H(t)|/|H_0|$ in \texttt{Float64} (solid) and \texttt{Float128} (dashed). \textbf{(b)}~Mean-longitude Laplace angle $\theta_{\mathrm{Lap}}(t) = \lambda_1 - 3\lambda_2 + 2\lambda_3$ (top) and phase residual $|\theta_{\mathrm{Lap}}(t) - \theta_{\mathrm{Lap}}^{\mathrm{BM6}}(t)|$ (bottom). \textbf{(c)}~Readout coordinate ablation (top) and circular-vs.-forced initialization ablation (bottom). \textbf{(d)}~\texttt{Float64} vs.\ \texttt{Float128} energy and 9-generator Lie--Bargmann residuals.}
    \label{fig:s1_jovian_laplace}
    \label{fig:secular_brouwer_laplace}
\end{figure}
\medskip\noindent\textbf{Protocol S2 : Inner Solar System Secular Stability, Covariant Shadow Telescoping Across 264 Dyadic Junctions ($\propto t^0$) vs.\ Brouwer Noether Diffusion ($\propto t^{1/2}$), and Authentic Laskar NAFF Frequency Metrology ($N=5$, Figure~\ref{fig:s2_solar_system_brouwer_naff}).}
We integrate a five-body Newtonian point-mass model of the inner Solar System comprising the Sun ($m_0 = 1.0$) and the four terrestrial planets Mercury ($m_1 = 1.6601 \times 10^{-7}, a_1 = 0.3871\text{ AU}, e_1 = 0.2056$), Venus ($m_2 = 2.4478 \times 10^{-6}, a_2 = 0.7233\text{ AU}, e_2 = 0.0068$), Earth--Moon ($m_3 = 3.0404 \times 10^{-6}, a_3 = 1.0000\text{ AU}, e_3 = 0.0167$), and Mars ($m_4 = 3.2272 \times 10^{-7}, a_4 = 1.5237\text{ AU}, e_4 = 0.0934$) in mutually inclined 3D orbits over $N_{\mathrm{steps}} = 4{,}000$ macro-steps of size $\Delta t = 0.05\text{ rad} \approx 2.90\text{ days}$ ($T_{\max} = 200.0\text{ rad} = 31.83\text{ yr} \approx 132.2$ Mercury orbits), recording the physical state at every macro-step ($N_{\mathrm{samples}} = 4{,}001$, i.e., $30.27\text{ samples/Mercury orbit}$, Figure~\ref{fig:s2_solar_system_brouwer_naff}):
\begin{itemize}
    \item \textit{Heliocentric Atlas Invariance, Covariant Gauge Telescoping Across 264 Dyadic Junctions, and Bounded Shadow-Hamiltonian Plateau ($\propto t^0$, Panels~\ref{fig:s2_solar_system_brouwer_naff}a, d)}: Throughout the $31.83\text{-yr}$ integration, the Hill-weighted Delaunay--Jacobi admissibility guard $\gamma_{\mathrm{Hill}}$ at canonical Voronoi threshold $\gamma = 1.0$ prevents spurious heliocentric-to-synodic leaf flips during planetary inferior conjunctions, maintaining $n_{\mathrm{flips}} = 0$ and regular secular eccentricity hodographs $(\Delta h_i, \Delta k_i) = \bigl(\Delta(e_i\cos\varpi_i), \Delta(e_i\sin\varpi_i)\bigr)$ (Figure~\ref{fig:s2_solar_system_brouwer_naff}d). Simultaneously, because Mercury's eccentric orbit ($e_1 = 0.2056$) modulates its local dynamical timescale between perihelion and aphelion across each of its $\sim 132.2$ revolutions, the symmetric dyadic scheduler dynamically switches between $k_{\mathrm{sym}} = 2$ ($h_{\mathrm{sub}} = 0.0125$) and $k_{\mathrm{sym}} = 3$ ($h_{\mathrm{sub}} = 0.00625$), executing \textbf{$264$ dyadic step-size junctions} ($0$ saturation events).

    At every junction $h_A \to h_B$, the macro-step propagator applies the covariant gauge connection $T_{(A,h_A)\to(A,h_B)}^{(\mathrm{PBM4})} = \bigl(\mathcal{Z}_4^{\mathrm{BM4},(A)}(h_B)\bigr)^{-1} \circ \mathcal{Z}_4^{\mathrm{BM4},(A)}(h_A)$, preserving exact shadow-manifold telescoping across step-size transitions: the running RMS relative energy error $\operatorname{RMS}_{[0,t]}(\Delta H/H_0)$ in Figure~\ref{fig:s2_solar_system_brouwer_naff}a (shown in both \texttt{Float64} solid and \texttt{Float128} dashed curves) saturates onto a horizontal periodic plateau ($\propto t^0$), where $\mathrm{PBM}_4(\mathcal{Z}_4^{\mathrm{BM4}})$ reduces the cumulative RMS energy error by \textbf{$22.42\times$ in \texttt{Float64} and $27.33\times$ in \texttt{Float128}} over bare $\mathrm{BM}_4^{(7)}$ ($8.28 \times 10^{-14}$ vs.\ $1.86 \times 10^{-12}$ in \texttt{Float64}, and $6.74 \times 10^{-14}$ vs.\ $1.84 \times 10^{-12}$ in \texttt{Float128}; \textbf{$19.18\times$ \texttt{Float64} / $22.70\times$ \texttt{Float128}} in maximum error, $1.55 \times 10^{-13}$ vs.\ $2.97 \times 10^{-12}$ in \texttt{Float64} and $1.31 \times 10^{-13}$ vs.\ $2.96 \times 10^{-12}$ in \texttt{Float128}), while $\mathrm{BM}_6^{(11)}$ reaches $\mathrm{RMS} = 1.19 \times 10^{-14}$ ($\max = 3.27 \times 10^{-14}$) in \texttt{Float64} and $\mathrm{RMS} = 7.15 \times 10^{-18}$ ($\max = 1.18 \times 10^{-17}$) in \texttt{Float128}. Expressed in physical time units via $\beta_{E,\mathrm{yr}^{-1}} = 2\pi\,\beta_{E,\mathrm{rad}^{-1}}$, the linear least-squares secular energy slope is bounded to:
    \begin{equation}
        \beta_E(\mathrm{PBM}_4) = \frac{\mathrm{d}}{\mathrm{d}t}\Bigl(\frac{\Delta H}{H_0}\Bigr) = +1.05 \times 10^{-15}\text{ yr}^{-1}, \qquad \beta_E(\mathrm{BM}_6^{(11)}) = -1.67 \times 10^{-16}\text{ yr}^{-1},
    \end{equation}
    confirming that covariant gauge telescoping across dyadic step-size transitions eliminates secular energy drift.
    \item \textit{Unbiased Brouwer Stochastic Diffusion ($\propto t^{1/2}$) and Sub-Stage Arithmetic Scaling on Exact Lie--Bargmann Invariants (Panel~\ref{fig:s2_solar_system_brouwer_naff}b)}: In contrast to the Hamiltonian, the 9 Galilean Lie--Bargmann Noether invariants have zero truncation error and isolate pure IEEE-754 floating-point roundoff accumulation. Figure~\ref{fig:s2_solar_system_brouwer_naff}b compares the unclamped running-RMS drift of all three invariant families ($\|\Delta\mathbf{L}\|_2/\|\mathbf{L}_0\|_2$ in red, $\|\Delta\mathbf{G}\|_2/M_{\mathrm{tot}}$ in green, and $\|\Delta\mathbf{P}\|_2/M_{\mathrm{tot}}$ in blue) between $\mathrm{PBM}_4$ (solid curves, $s_{\mathrm{eval}} = 6$) and $\mathrm{BM}_6^{(11)}$ (dashed curves, $s_{\mathrm{eval}} = 10$), revealing two structural properties: (i)~both schemes occupy the exact same three invariant bands ($\sim 10^{-14}$ for $\Delta\mathbf{L}$, $\sim 10^{-17}$ for $\Delta\mathbf{G}$, and $\sim 10^{-19}$ for $\Delta\mathbf{P}$, dropping to $\le 3.84 \times 10^{-32}$ and $\le 6.36 \times 10^{-39}$ over the full $4{,}000$-step \texttt{Float128} integration), where the four- to five-order-of-magnitude suppression of $\|\Delta\mathbf{P}\|_2/M_{\mathrm{tot}}$ and $\|\Delta\mathbf{G}\|_2/M_{\mathrm{tot}}$ below $\epsilon_{\mathrm{mach}} = 2.22 \times 10^{-16}$ stems from the decoupling of the root barycenter node in the Jacobi tree and the weighting of planetary velocity roundoff by the planet-to-Sun mass ratios $m_i/M_{\mathrm{tot}} \sim 10^{-6}$; and

    (ii)~all curves follow Brouwer's unbiased random-walk scaling \cite{hairer2006geometric}:
    \[
        \sigma_{\mathrm{Noether}}(t) \propto \epsilon_{\mathrm{mach}}\sqrt{s_{\mathrm{eval}}\,n_{\mathrm{steps}}} \propto t^{1/2},
    \]
    with $\mathrm{BM}_6^{(11)}$ reaching $\max\|\Delta\mathbf{L}\|_2/\|\mathbf{L}_0\|_2 = 1.45 \times 10^{-14}$, $\max\|\Delta\mathbf{G}\|_2/M_{\mathrm{tot}} = 1.55 \times 10^{-17}$, and $\max\|\Delta\mathbf{P}\|_2/M_{\mathrm{tot}} = 1.71 \times 10^{-19}$ ($\mathrm{PBM}_4$: $1.33 \times 10^{-14}$, $2.27 \times 10^{-17}$, $1.68 \times 10^{-19}$; $\mathcal{S}_2$: $9.05 \times 10^{-15}$, $1.20 \times 10^{-17}$, $8.26 \times 10^{-20}$).
    \item \textit{Authentic Laskar NAFF Frequency Metrology and KAM Symplectic Conjugacy Invariance (Panel~\ref{fig:s2_solar_system_brouwer_naff}c)}: To quantify phase-space frequency preservation without synthetic proxies, we apply Laskar's Numerical Analysis of Fundamental Frequencies (\textit{NAFF}) with order-$p=3$ Hanning filtering (Section~\ref{sub:laskar_dual_convergence}) directly to Mercury's complex heliocentric orbital signal ($N_{\mathrm{samples}} = 4{,}001$, $30.27\text{ samples/orbit}$, well above the Nyquist threshold):
    \[
        z_{\mathrm{orb}}(t_k) = \bigl(x_1(t_k) - x_0(t_k)\bigr) + i\bigl(y_1(t_k) - y_0(t_k)\bigr),
    \]
    against a $2\times$-finer dyadic $\mathrm{BM}_6^{(11)}$ reference trajectory ($\Delta t_{\mathrm{ref}} = 0.025\text{ rad}$, $8{,}000$ macro-steps, $\nu_{\mathrm{ref}} = 4.15210957695\text{ rad/time unit}$). Figure~\ref{fig:s2_solar_system_brouwer_naff}c reports three complementary observables: the full-horizon fundamental frequency error $|\nu_{\mathrm{NAFF}} - \nu_{\mathrm{ref}}|/\nu_{\mathrm{ref}}$, the complex waveform $L^2$ error $\|z_{\mathrm{orb}} - z_{\mathrm{orb}}^{\mathrm{ref}}\|_2 / \|z_{\mathrm{orb}}^{\mathrm{ref}}\|_2$, and the two-window numerical frequency drift between $W_1 = [0, T_{\max}/2]$ and $W_2 = [T_{\max}/2, T_{\max}]$ after subtracting the true physical secular planetary precession shift $\Delta\nu_{1\to 2}^{\mathrm{phys}}/\nu_1 = 1.78 \times 10^{-8}$:
    \begin{equation}
        \frac{\delta\nu_{1\to 2}^{\mathrm{num}}}{\nu_{\mathrm{ref}}} = \frac{\bigl|(\nu_{W_2} - \nu_{W_1}) - (\nu_{W_2}^{\mathrm{ref}} - \nu_{W_1}^{\mathrm{ref}})\bigr|}{\nu_{\mathrm{ref}}}.
    \end{equation}
    Across the splitting hierarchy under adaptive dyadic substepping ($k_{\mathrm{sym}} \in \{2, 3\}$), the fundamental frequency error drops by \textbf{$2{,}120\times$} from $\mathcal{S}_2$ ($7.85 \times 10^{-10}$, $\delta\nu_{1\to 2}^{\mathrm{num}}/\nu_{\mathrm{ref}} = 1.17 \times 10^{-8}$, $L^2 = 2.35 \times 10^{-6}$) and from $\mathcal{Y}_4$ ($2.00 \times 10^{-12}$, improved by $166\times$ over unsubstepped $k=0$, $\delta\nu_{1\to 2}^{\mathrm{num}}/\nu_{\mathrm{ref}} = 4.17 \times 10^{-12}$, $L^2 = 2.36 \times 10^{-9}$) right down to the $\sim 3.7\text{--}4.4 \times 10^{-13}$ spectral resolution floor of the $4{,}001$-sample order-3 Hanning window on $\mathrm{BM}_4^{(7)}$ ($4.40 \times 10^{-13}$, $\delta\nu_{1\to 2}^{\mathrm{num}}/\nu_{\mathrm{ref}} = 1.96 \times 10^{-13}$, $L^2 = 7.21 \times 10^{-11}$), $\mathrm{PBM}_4(\mathcal{Z}_4^{\mathrm{BM4}})$ ($3.70 \times 10^{-13}$, $\delta\nu_{1\to 2}^{\mathrm{num}}/\nu_{\mathrm{ref}} = 9.90 \times 10^{-14}$, $L^2 = 3.78 \times 10^{-11}$), and $\mathrm{BM}_6^{(11)}$ ($4.23 \times 10^{-13}$, $\delta\nu_{1\to 2}^{\mathrm{num}}/\nu_{\mathrm{ref}} = 6.67 \times 10^{-14}$, $L^2 = 6.83 \times 10^{-11}$), with:
    \begin{equation}
        \frac{|\nu_{\mathrm{NAFF}}(\mathrm{PBM}_4) - \nu_{\mathrm{NAFF}}(\mathrm{BM}_4^{(7)})|}{\nu_{\mathrm{ref}}} = 7.02 \times 10^{-14}.
    \end{equation}
    This near-identity directly reflects the \textbf{KAM symplectic conjugacy property} \cite{laskar1993frequency,hairer2006geometric}: because the processed propagator
    \[
        \hat{\Phi}_h = \mathcal{Z}_4^{\mathrm{BM4}}(h) \circ \mathcal{S}_h \circ \bigl(\mathcal{Z}_4^{\mathrm{BM4}}(h)\bigr)^{-1}
    \]
    is canonically conjugate to the bare kernel $\mathcal{S}_h$ via the near-identity symplectomorphism $\mathcal{Z}_4^{\mathrm{BM4}}(h)$ and gauge connection $T_{(A,h_A)\to(A,h_B)}^{(\mathrm{PBM4})}$, both maps share identical asymptotic rotation numbers on invariant KAM tori while $\mathcal{Z}_4^{\mathrm{BM4}}(h)$ deforms the shadow torus into $22.42\times$ (\texttt{Float64}) and $27.33\times$ (\texttt{Float128}) closer RMS alignment with the constant-energy hypersurface.
\end{itemize}
\pythonrefright{py_bench_gen05}{{benchmarks/run, figures/plot}_generator_05_*.py}{Float64 / Float128}

\begin{figure}[!htbp]
    \centering
    \includegraphics[width=\textwidth]{chapters/ch07_benchmarks_and_lean4/figures/fig_ch07_s2_solar_system_brouwer_naff.png}
    \caption{\textbf{Protocol S2: Inner Solar System secular stability, covariant shadow telescoping across $264$ dyadic junctions ($\propto t^0$) vs.\ Brouwer Noether diffusion ($\propto t^{1/2}$), and authentic Laskar NAFF frequency confinement ($N=5$, $4{,}000$ macro-steps, $31.83\text{ yr} \approx 132.2$ Mercury orbits).} \textbf{(a)}~Running RMS energy error $\operatorname{RMS}_{[0,t]}(\Delta H/H_0)$ in \texttt{Float64} (solid) and \texttt{Float128} (dashed). \textbf{(b)}~Unclamped running-RMS drift of Lie--Bargmann invariants ($\|\Delta\mathbf{L}\|_2/\|\mathbf{L}_0\|_2$, $\|\Delta\mathbf{G}\|_2/M_{\mathrm{tot}}$, $\|\Delta\mathbf{P}\|_2/M_{\mathrm{tot}}$). \textbf{(c)}~Order-3 Hanning Laskar NAFF frequency error, two-window drift, and waveform $L^2$ error on Mercury's $z_{\mathrm{orb}}(t)$. \textbf{(d)}~Secular planetary eccentricity hodographs $(\Delta h_i, \Delta k_i)$ ($0$ synodic flips).}
    \label{fig:s2_solar_system_brouwer_naff}
\end{figure}

\begin{table}[!htbp]
    \centering
    \caption{\textbf{Certified quantitative long-term secular stability, 9-generator Lie--Bargmann metrology, and Laskar NAFF frequency metrics under canonical Voronoi admissibility ($\gamma = 1.0$) and active dyadic substepping (\S\ref{sub:secular_stability_laplace}).} Protocol~S1 integrates the 4-body Jovian $4\mathbin{:}2\mathbin{:}1$ Laplace resonance ($200$ Io orbits, $k_{\mathrm{sym}} \equiv 3$), and Protocol~S2 integrates the 5-body Inner Solar System ($31.83\text{ yr}$, $264$ dyadic junctions, $0$ synodic flips).}
    \label{tab:campaign_III_kpi_table}
    \vspace{0.3em}
    \renewcommand{\arraystretch}{1.22}
    \resizebox{\textwidth}{!}{%
    \begin{tabular}{lcccccc}
        \toprule
        \textbf{Integrator} & \makecell{\textbf{S1 Laplace Error}\\$\max|\Delta H|/|H_0|$ \textbf{F64 / F128}\\(\textbf{RMS F64 / F128}, $k=3$)} & \makecell{\textbf{S1 Laplace Angle}\\$\theta_{\mathrm{Lap}} = \lambda_1 - 3\lambda_2 + 2\lambda_3$\\(Mean Longitudes)} & \makecell{\textbf{S2 Solar System}\\$\operatorname{RMS}(\Delta H/H_0)$ \textbf{F64 / F128}\\($\max|\Delta H|/|H_0|$ \textbf{F64 / F128})} & \makecell{\textbf{S2 Secular Slope}\\$\beta_E\text{ [yr}^{-1}\text{]}$\\($2\pi\,\beta_{E,\mathrm{rad}^{-1}}$)} & \makecell{\textbf{S2 Laskar NAFF Error}\\$|\nu_{\mathrm{NAFF}} - \nu_{\mathrm{ref}}|/\nu_{\mathrm{ref}}$\\($\delta\nu_{1\to 2}^{\mathrm{num}}/\nu_{\mathrm{ref}}$)} & \makecell{\textbf{Lie--Bargmann (\texttt{Float64})}\\$\max\|\Delta\mathbf{L}\|_2/\|\mathbf{L}_0\|_2\text{ (S1 / S2)}$\\($\max\|\Delta\mathbf{P}\|_2/M_{\mathrm{tot}}\text{ S1 / S2}$)} \\
        \midrule
        $\mathcal{S}_2$ (Strang, $s=1$) & \makecell{$3.98 \times 10^{-7}$ / $3.98 \times 10^{-7}$\\($2.47 \times 10^{-7}$ / $2.47 \times 10^{-7}$)} & \makecell{$183.41^\circ \pm 3.53^\circ$\\(Librated at $k=3$)} & \makecell{$1.05 \times 10^{-6}$ / $1.05 \times 10^{-6}$\\($1.79 \times 10^{-6}$ / $1.79 \times 10^{-6}$)} & $+1.61 \times 10^{-10}$ & \makecell{$7.85 \times 10^{-10}$\\($1.17 \times 10^{-8}$)} & \makecell{$2.51 \times 10^{-14}$ / $9.05 \times 10^{-15}$\\($2.32 \times 10^{-18}$ / $8.26 \times 10^{-20}$)} \\
        \addlinespace[0.2em]
        $\mathcal{Y}_4$ (Yoshida, $s=3$) & \makecell{$1.98 \times 10^{-10}$ / $1.98 \times 10^{-10}$\\($1.20 \times 10^{-10}$ / $1.20 \times 10^{-10}$)} & \makecell{$179.85^\circ \pm 0.80^\circ$\\(Locked)} & \makecell{$5.67 \times 10^{-10}$ / $5.67 \times 10^{-10}$\\($9.39 \times 10^{-10}$ / $9.39 \times 10^{-10}$)} & $-6.43 \times 10^{-14}$ & \makecell{$2.00 \times 10^{-12}$\\($4.17 \times 10^{-12}$)} & \makecell{$1.89 \times 10^{-14}$ / $1.30 \times 10^{-14}$\\($4.79 \times 10^{-18}$ / $1.05 \times 10^{-19}$)} \\
        \addlinespace[0.2em]
        $\mathrm{BM}_4^{(7)}$ (Bare, $s=6$) & \makecell{$6.96 \times 10^{-13}$ / $6.67 \times 10^{-13}$\\($4.31 \times 10^{-13}$ / $4.09 \times 10^{-13}$)} & \makecell{$179.85^\circ \pm 0.80^\circ$\\(Locked)} & \makecell{$1.86 \times 10^{-12}$ / $1.84 \times 10^{-12}$\\($2.97 \times 10^{-12}$ / $2.96 \times 10^{-12}$)} & $-1.31 \times 10^{-15}$ & \makecell{$4.40 \times 10^{-13}$\\($1.96 \times 10^{-13}$)} & \makecell{$1.49 \times 10^{-14}$ / $2.19 \times 10^{-14}$\\($7.05 \times 10^{-18}$ / $1.53 \times 10^{-19}$)} \\
        \addlinespace[0.2em]
        $\mathbf{PBM_4(\mathcal{Z}_4^{\mathrm{BM4}})}$ ($s=6$) & \makecell{$\mathbf{1.20 \times 10^{-13}}$ / $\mathbf{3.80 \times 10^{-14}}$\\($\mathbf{5.77 \times 10^{-14}}$ / $\mathbf{1.25 \times 10^{-14}}$)} & \makecell{$\mathbf{179.85^\circ \pm 0.80^\circ}$\\(Locked)} & \makecell{$\mathbf{8.28 \times 10^{-14}}$ / $\mathbf{6.74 \times 10^{-14}}$\\($\mathbf{1.55 \times 10^{-13}}$ / $\mathbf{1.31 \times 10^{-13}}$)} & $\mathbf{+1.05 \times 10^{-15}}$ & \makecell{$\mathbf{3.70 \times 10^{-13}}$\\($\mathbf{9.90 \times 10^{-14}}$, KAM conj.)} & \makecell{$\mathbf{5.28 \times 10^{-14}}$ / $\mathbf{1.33 \times 10^{-14}}$\\($\mathbf{2.49 \times 10^{-18}}$ / $\mathbf{1.68 \times 10^{-19}}$)} \\
        \addlinespace[0.2em]
        $\mathrm{BM}_6^{(11)}$ (Bare, $s=10$) & \makecell{$8.14 \times 10^{-14}$ / $1.90 \times 10^{-18}$\\($2.88 \times 10^{-14}$ / $1.02 \times 10^{-18}$)} & \makecell{$179.85^\circ \pm 0.80^\circ$\\(Locked)} & \makecell{$1.19 \times 10^{-14}$ / $7.15 \times 10^{-18}$\\($3.27 \times 10^{-14}$ / $1.18 \times 10^{-17}$)} & $-1.67 \times 10^{-16}$ & \makecell{$4.23 \times 10^{-13}$\\($6.67 \times 10^{-14}$)} & \makecell{$2.23 \times 10^{-14}$ / $1.45 \times 10^{-14}$\\($7.16 \times 10^{-18}$ / $1.71 \times 10^{-19}$)} \\
        \bottomrule
    \end{tabular}%
    }
\end{table}
\clearpage