% --- Chapter 5 Preamble & Section 5.1 Content ---
Having established the continuous geometric foundations of the Fulton--MacPherson compactification $\mathbf{Q}[N]$, the $\mathcal{C}^{1,1}$-regularized screened apex potential $W_{\mathrm{apex}}^{\mathrm{screen}}$, and the synchronous supernode regularization in Section~\ref{ch:ch04_fulton_macpherson_splitting}, we now pass from continuous Hamiltonian geometry to discrete symplectic time-stepping. The overarching objective of this chapter is to construct a hierarchy of unsoftened ($\varepsilon = 0, d_0 = 0$), structure-preserving geometric integrators on the extended cotangent bundle $\mathcal{P}_{\mathrm{ext}} = T^*\mathbf{Q}[N] \times \mathbb{R}_t \times \mathbb{R}_E$ that conserve all ten Lie--Bargmann first integrals to machine precision while traversing exact binary and cluster collisions in a finite number of physical steps.

Our exposition proceeds in six logical stages. Section~\ref{sec:strang_propagator} constructs the foundational 7-operator hierarchical palindromic Strang propagator $\mathcal{S}_2(h)$ and establishes its unconditional extended symplecticity and discrete Galilean invariance. Section~\ref{sec:suzuki_barrier_mcgehee} confronts the Sheng--Suzuki impossibility barrier ($w_i < 0$ for order $p \ge 3$) and proves that the real-analytic invertibility of universal Sundman--Stumpff flows coupled with McGehee's unitary collision bounce involution ($\mathcal{B}_C^2 = \operatorname{Id}_{\mathbf{M}_C}$) renders retrograde sub-steps completely singularity-free, enabling Yoshida's fourth-order scheme $\mathcal{S}_4(h)$. To suppress the large backward excursion of minimal triplets, Section~\ref{sec:bm4_seven_stages} deploys the spectrally optimized 7-stage Blanes--Moan scheme $\mathrm{BM}_4$ and proves general order-$2p$ exponential shadow Hamiltonian confinement under zero softening. Section~\ref{sec:bm6_eleven_stages} extends this composition architecture to the 11-stage sixth-order scheme $\mathrm{BM}_6^{(11)}$ and establishes the bi-graded Lie algebra filtration $\mathcal{L}_{n,m}(A_{\mathrm{JK}}, W_{\mathrm{apex}}^{\mathrm{screen}})$ governing hierarchical error commutators. Section~\ref{sec:symplectic_correctors} then investigates symplectic boundary processing: we prove the \emph{Universal Obstruction Theorem} prohibiting global order-6 configuration correctors on $T^*\mathbf{Q}[N]$, construct the domain-restricted tidal corrector $K_4$ ($\mathrm{PBM}_4$) on $\mathcal{D}_{\mathrm{reg}}$ with its irreducible quadratic floor $\mathcal{O}(\varepsilon^2 h^4)$, and derive the covariant cotangent-lifted gauge transition operator $T_{A \to B}^{(\mathrm{PBM4})}$ across Jacobi charts. Finally, Section~\ref{sec:dyadic_multirate_scheduling} synthesizes these propagators into a time-reversible dyadic multirate scheduler governed by Skeel--Gear macro-boundary synchronization and the \emph{Leaf/Apex Decoupling Theorem}.

\subsection{The 7-Operator Hierarchical Palindromic Strang Propagator ($\mathcal{S}_2$)}
\label{sec:strang_propagator}

The foundational building block of our geometric numerical integrator on the compactified Fulton--MacPherson configuration space $\mathbf{Q}[N]$ (Section~\ref{ch:ch04_fulton_macpherson_splitting}) is the second-order symmetric Strang splitting propagator~\cite{strang1968construction} (independently developed by Marchuk~\cite{marchuk1968applications} and rooted in Trotter's product formula~\cite{trotter1959product}). In classical celestial mechanics, operator splitting for unsoftened $N$-body dynamics has historically foundered on close encounters: when pairwise separations $\rho_C \to 0$, naive kinetic-potential splittings suffer unbounded truncation error spikes. While advanced splitting schemes for perturbed systems $H = A + \varepsilon B$ (most notably Laskar and Robutel~\cite{laskar2001large} and Blanes, Casas, Farr{\'e}s, Laskar, Makazaga and Murua~\cite{blanes2013optimized}) achieve extraordinary efficiency in regular, quasi-circular secular regimes with a uniform distance barrier $\inf_{t} r_{ij}(t) \ge r_{\min} > 0$, unsoftened close encounters ($\rho_C \to 0$) violate $\varepsilon \ll 1$ and drive catastrophic Lie commutator explosions. On the compactified cotangent bundle $\mathcal{P}_{\mathrm{ext}} = T^*\mathbf{Q}[N] \times \mathbb{R}_t \times \mathbb{R}_E$, however, the continuous Hamiltonian established in Section~\ref{ch:ch04_fulton_macpherson_splitting} decomposes into autonomous, regular, and desingularized vector fields. 

In this section, we construct the discrete 7-operator hierarchical palindromic Strang scheme $\mathcal{S}_2(h)$, prove its unconditional extended symplecticity on $(\mathcal{P}_{\mathrm{ext}}, \Omega_{\mathrm{ext}})$, demonstrate the exact cancellation of all odd-order dissipative error terms in its Baker--Campbell--Hausdorff expansion, and verify the exact discrete conservation (with identically vanishing truncation error) of all ten generators of the centrally extended Lie--Bargmann symmetry algebra.

% ------------------------------------------------------------------------------
\subsubsection{Extended Phase Space and Elementary Lie Propagators}
\label{sub:strang_phase_space_propagators}

Recall from Section~\ref{sub:lie_bargmann_and_palindromy} (Proposition~\ref{prop:extended_palindromic_symplecticity}) that the dynamics is formulated on the $(6N+2)$-dimensional extended Poincaré--Sundman phase space:
\begin{equation}
    \mathcal{P}_{\mathrm{ext}} = T^*\big(\mathbf{Q}[N] \times \mathbb{R}^3_{\mathrm{cdm}}\big) \times \mathbb{R}_t \times \mathbb{R}_E \cong T^*\mathbf{Q}[N] \times T^*\mathbb{R}^3_{\mathrm{cdm}} \times \mathbb{R}_t \times \mathbb{R}_E,
    \label{eq:strang_extended_phase_space}
\end{equation}
which identifies smoothly with the standard cotangent bundle $(T^*(\mathbb{R}^{3N} \setminus \Delta)) \times \mathbb{R}_t \times \mathbb{R}_E$ on the dense open collisionless stratum $\mathbf{Q}[N]^\circ$, while rigorously incorporating the boundary strata $\partial \mathbf{Q}[N]$ associated with binary and multiple collisions.
Equipped with the canonical configuration-momentum conjugate pair $(q_t, p_t) = (t, -E)$, where $t \in \mathbb{R}$ represents the physical coordinate time and $E \in \mathbb{R}$ is the conjugate energy variable, the extended Liouville $1$-form is:
\begin{equation}
    \lambda_{\mathrm{ext}} = \sum_{i=1}^N \mathbf{p}_i \cdot \mathrm{d}\mathbf{q}_i + p_t\,\mathrm{d}t = \sum_{i=1}^N \mathbf{p}_i \cdot \mathrm{d}\mathbf{q}_i - E\,\mathrm{d}t,
    \label{eq:strang_extended_liouville}
\end{equation}
inducing the canonical closed non-degenerate extended Darboux $2$-form:
\begin{equation}
    \Omega_{\mathrm{ext}} \coloneqq -\mathrm{d}\lambda_{\mathrm{ext}} = \Omega_0 - \mathrm{d}p_t \wedge \mathrm{d}t = \Omega_0 - (-\mathrm{d}E) \wedge \mathrm{d}t = \Omega_0 + \mathrm{d}E \wedge \mathrm{d}t, \quad \mathrm{d}\Omega_{\mathrm{ext}} = 0,
    \label{eq:strang_extended_darboux_form}
\end{equation}
where $\Omega_0 = \sum_{i=1}^N \sum_{\alpha=1}^3 \mathrm{d}q_i^\alpha \wedge \mathrm{d}p_i^\alpha$ is the standard spatial Darboux form on $T^*\mathbb{R}^{3N}$. The canonical extended Poisson brackets satisfy:
\begin{equation}
    \{q_i^\alpha, p_j^\beta\} = \delta_{ij} \delta^{\alpha\beta}, \qquad \{t, -E\} = 1 \iff \{t, E\} = -1,
    \label{eq:strang_canonical_brackets}
\end{equation}
all other fundamental brackets vanishing identically. The autonomous extended Hamiltonian is given by:
\begin{equation}
    \mathcal{H}_{\mathrm{ext}}(\mathbf{q}, \mathbf{p}, t, E) = H(\mathbf{q}, \mathbf{p}) - E,
    \label{eq:strang_extended_hamiltonian}
\end{equation}
and physical motions correspond to trajectories confined to the regular zero-energy level hypersurface $\Sigma_0 = \{ (\mathbf{q}, \mathbf{p}, t, E) \in \mathcal{P}_{\mathrm{ext}} \mid \mathcal{H}_{\mathrm{ext}} = 0 \}$.

In any adapted Jacobi chart $U_{\mathcal{T}} \subset \mathbf{Q}[N]$ associated with a maximal flag tree $\mathcal{T} \in \operatorname{Trees}(N)$ (Section~\ref{sub:maximal_flag_jacobi_trees}), the physical Hamiltonian is rigorously structured as a \emph{Fulton--MacPherson hierarchical splitting}:
\begin{equation}
    H(\mathbf{r}, \boldsymbol{\pi}) = A_{\mathrm{hier}}(\mathbf{r}, \boldsymbol{\pi}) + W_{\mathrm{apex}}^{\mathrm{screen}}(\mathbf{r}) = A_{\mathrm{JK}}(\mathbf{r}, \boldsymbol{\pi}) + W_{\mathrm{apex}}^{\mathrm{screen}}(\mathbf{r}) + V_{\mathrm{kepler}}^{\mathrm{screened}}(\mathbf{r}),
    \label{eq:strang_five_operator_splitting}
\end{equation}
where the unperturbed hierarchical Keplerian operator $A_{\mathrm{JK}}$ aggregates all decoupled kinematic and Keplerian degrees of freedom:
\begin{equation}
    A_{\mathrm{JK}}(\mathbf{r}, \boldsymbol{\pi}) \coloneqq T_0(\boldsymbol{\pi}_0) + T_{\mathrm{comp}}(\boldsymbol{\pi}) + H_{\mathrm{leaf}}(\mathbf{r}, \boldsymbol{\pi}) + H_{\mathrm{super}}(\mathbf{r}, \boldsymbol{\pi}).
    \label{eq:conormal_ajk_operator_def}
\end{equation}
Crucially, in canonical Jacobi coordinates, the relative pairs and barycentric channels are mututally orthogonal under the Darboux form $\Omega_0 = \sum \mathrm{d}\mathbf{r}_k \wedge \mathrm{d}\boldsymbol{\pi}_k$, ensuring that all constituent channels of $A_{\mathrm{JK}}$ Poisson-commute pairwise:
\begin{equation}
    \{H_k^{\mathrm{Kepler}}, H_l^{\mathrm{Kepler}}\} \equiv 0 \quad (\forall k \ne l) \implies \{A_{\mathrm{JK}}, A_{\mathrm{JK}}\} \equiv 0.
    \label{eq:ajk_exact_abelian_commutativity}
\end{equation}
Consequently, the unperturbed flow $\exp(\delta t X_{A_{\mathrm{JK}}})$ factorizes into the exact direct product of its autonomous channels without any operator splitting error.

\begin{rembox}[Geometric Grounding: Calaque--Safronov Cotangent Specialization and Hierarchical Splitting]
\label{rem:calaque_safronov_grounding}
In the foundational framework of Fulton--Verdier~\cite{fulton1984intersection} and its shifted cotangent lift established by Calaque \& Safronov~\cite{calaque2024deformation}, the deformation to the normal cone $\mathcal{D}_\Delta Q \to \mathbb{A}^1$ induces a relative symplectic family $T^*(\mathcal{D}_\Delta Q / \mathbb{A}^1) \to \mathbb{A}^1$ whose special fiber at $t = 0$ specializes to the cotangent bundle of the normal cone $T^*(N_\Delta Q) \cong T^*\Delta \times T^*(N_\Delta Q / \Delta)$. 
Under this geometric lens, our operator splitting~\eqref{eq:strang_five_operator_splitting} realizes the canonical base-fiber polarization:
\begin{equation}
    H = H_{\mathrm{base}} + H_{\mathrm{normal}}^{\mathrm{apex}} + H_{\mathrm{super}}^{\mathrm{sync}},
\end{equation}
where $H_{\mathrm{base}} = T_0 + T_{\mathrm{comp}}$ generates the regular horizontal drift on the base stratum ($t = 1$), $H_{\mathrm{normal}}^{\mathrm{apex}} = H_{\mathrm{leaf}} + \mathcal{K}^{\mathrm{apex}}$ governs the decoupled two-body dynamics and vertical cotangent interaction kicks along the fibers ($t = 0$), and $H_{\mathrm{super}}^{\mathrm{sync}}$ regulates multiple collision loci on the deeper strata of the Fulton--MacPherson tree.

\noindent\textbf{Hierarchical Structure and Error Grading:} An important consequence of this geometry is that on the regular hierarchical domain $\mathcal{D}_{\mathrm{reg}}$, the unperturbed operator $A_{\mathrm{JK}}$ is completely integrable by exact analytic quadratures (universal Kepler--Sundman flow), while the perturbation $W_{\mathrm{apex}}^{\mathrm{screen}}$ is scaled by the dimensionless tidal parameter $\varepsilon \sim (r_{\mathrm{leaf}} / R_{\mathrm{ext}})^3 \le 0.125 \ll 1$. This structure breaks fundamentally from classical Cartesian splittings $T(\mathbf{p}) + V(\mathbf{q})$: because internal two-body dynamics is analytically absorbed into $A_{\mathrm{JK}}$, all pure velocity error commutators $L_T^4(V) \in \mathcal{P}_4$ vanish identically, and the Baker--Campbell--Hausdorff error series is bi-graded in powers $(h^{2n}, \varepsilon^m)$ with no residual term in $\mathcal{O}(\varepsilon^0 h^k)$.
\end{rembox}

\noindent The five elementary Hamiltonian sub-flows governing the Fulton--MacPherson decomposition~\eqref{eq:strang_five_operator_splitting}--\eqref{eq:conormal_ajk_operator_def} are defined explicitly as follows:
\begin{enumerate}
    \item \textbf{Center-of-Mass Kinetic Drift ($T_0$):} Governed by $T_0(\boldsymbol{\pi}_0) = \frac{\|\boldsymbol{\pi}_0\|^2}{2 M_{\mathrm{tot}}}$ with total mass $M_{\mathrm{tot}} = \sum_{i=1}^N m_i$. Its exact Hamiltonian flow for a step $\delta t$ is a linear shear:
    \begin{equation}
        \Phi_{\delta t}^{\mathrm{drift}, (0)} : \left( \mathbf{R}_{\mathrm{cdm}}, \boldsymbol{\pi}_0 \right) \longmapsto \left( \mathbf{R}_{\mathrm{cdm}} + \delta t \, \frac{\boldsymbol{\pi}_0}{M_{\mathrm{tot}}}, \, \boldsymbol{\pi}_0 \right),
        \label{eq:strang_flow_cdm}
    \end{equation}
    leaving all internal Jacobi relative coordinates and $(t, E)$ invariant.
    \item \textbf{Composite-Node Kinematic Drift ($T_{\mathrm{comp}}$):} Governed by the sum of kinetic energies of internal cluster barycenters $T_{\mathrm{comp}}(\boldsymbol{\pi}) = \sum_{v \in \mathcal{V}_{\mathrm{comp}}(\mathcal{T})} \frac{\|\boldsymbol{\pi}_v\|^2}{2 \mu_v}$, where $\mu_v = \frac{M_{L(v)} M_{R(v)}}{M_{L(v)} + M_{R(v)}}$ is the reduced mass. Within the Fulton--MacPherson hierarchical operator $A_{\mathrm{JK}}$, intermediate cluster barycenters advance along exact straight-line Galilean shears $\Phi_{\delta t}^{\mathrm{drift}, (\mathrm{comp})}$, while their mutual gravitational attractions are smoothly absorbed into the screened apex potential $W_{\mathrm{apex}}^{\mathrm{screen}}$:
    \begin{equation}
        \Phi_{\delta t}^{\mathrm{drift}, (\mathrm{comp})} : \mathbf{r}_v \longmapsto \mathbf{r}_v + \delta t \, \frac{\boldsymbol{\pi}_v}{\mu_v}, \quad \boldsymbol{\pi}_v \longmapsto \boldsymbol{\pi}_v \quad \forall v \in \mathcal{V}_{\mathrm{comp}}(\mathcal{T}).
        \label{eq:strang_flow_comp}
    \end{equation}
    Since $\mu_v > 0$ and $T_{\mathrm{comp}}$ is purely quadratic, this flow is smooth and non-singular across $\mathbf{r}_v = \mathbf{0}$. We denote the aggregate kinematic drift by $\Phi_{\delta t}^{\mathrm{drift}} := \Phi_{\delta t}^{\mathrm{drift}, (0)} \circ \Phi_{\delta t}^{\mathrm{drift}, (\mathrm{comp})}$.
    \item \textbf{Leaf Keplerian Pairs Flow ($H_{\mathrm{leaf}}$):} Governed by the decoupled sum of isolated two-body Keplerian pairs $H_{\mathrm{leaf}} = \sum_{C \in \mathcal{V}_{\mathrm{leaf}}(\mathcal{T})} \left( \frac{\|\boldsymbol{\pi}_C\|^2}{2 \mu_C} - \frac{G m_i m_j}{\|r_C\|} \right)$. By Proposition~\ref{prop:disjoint_cotangent_commutation}, disjoint sub-clusters commute identically under the Poisson bracket ($[\mathcal{X}_{H_C}, \mathcal{X}_{H_{C'}}] = 0$), so the collective flow factorizes as:
    \begin{equation}
        \Phi_{\vec{\Delta \tau}}^{\mathrm{leaf}} = \prod_{C \in \mathcal{V}_{\mathrm{leaf}}(\mathcal{T})} \Phi_{\Delta \tau_C}^{(C)},
        \label{eq:strang_flow_leaf}
    \end{equation}
    where each $\Phi_{\Delta \tau_C}^{(C)}$ is integrated analytically in universal Stumpff series under Sundman's eccentric anomaly clock $\mathrm{d}t_C = \rho_C \mathrm{d}\tau_C$ and extends smoothly across collision $\mathbf{M}_C$ via the McGehee bounce involution $\mathcal{B}_C^2 = \mathrm{id}$.
    \item \textbf{Synchronous Supernodes Flow ($H_{\mathrm{super}}$):} Governed by compact multi-body sub-clusters ($\rho_S \le \delta_{\mathrm{cluster}}^{\mathrm{nat}}$, dynamically determined by Proposition~\ref{prop:intrinsic_compactness_threshold}) under Sundman's unfolded local clock $\mathrm{d}t_S = \chi_S^2 \mathrm{d}\tau_S$ (Section~\ref{sec:synchronous_supernodes}). The collective flow is:
    \begin{equation}
        \Phi_{\vec{\Delta \tau}_S}^{\mathrm{super}} = \prod_{S \in \mathcal{S}} \Phi_{\Delta \tau_S}^{(S)},
        \label{eq:strang_flow_super}
    \end{equation}
    guaranteeing finite Anti-Zeno crossing ($N_{\mathrm{sub}} < \infty$, Theorem~\ref{thm:anti_zeno_regularization}) and preserving extended Darboux area.
    \item \textbf{Screened Apex and Keplerian Potential Kick ($\mathcal{K}^{\mathrm{apex}}$):} Governed by the configuration potential $V_{\mathrm{kick}}(\mathbf{q}) = W_{\mathrm{apex}}^{\mathrm{screen}}(\mathbf{q}) + V_{\mathrm{kepler}}^{\mathrm{screened}}(\mathbf{q})$, whose apex component is universally $\mathcal{C}^{1,1}$-bounded (Theorem~\ref{thm:universal_c11_apex}) and whose screened Keplerian attraction is smooth away from collisions. Because $V_{\mathrm{kick}}$ depends purely on configuration coordinates $\mathbf{q}$, its Hamiltonian vector field $\mathcal{X}_{V_{\mathrm{kick}}}$ is vertical along the cotangent fibers $T^*_{\mathbf{q}} \mathbf{Q}[N]$ (spanning the vertical subbundle $V(T^*\mathbf{Q}[N]) = \ker \mathrm{d}\pi$). In particular, any two configuration kick components $V_a, V_b$ commute identically: $[\mathcal{X}_{V_a}, \mathcal{X}_{V_b}] = 0$. Its exact Hamiltonian flow for a time step $\delta t$ is a vertical cotangent shear:
    \begin{equation}
        \mathcal{K}_{\delta t}^{\mathrm{apex}} : \mathbf{q}_i \longmapsto \mathbf{q}_i, \quad \mathbf{p}_i \longmapsto \mathbf{p}_i - \delta t \, \nabla_{\mathbf{q}_i} \big( W_{\mathrm{apex}}^{\mathrm{screen}}(\mathbf{q}) + V_{\mathrm{kepler}}^{\mathrm{screened}}(\mathbf{q}) \big) \quad (i = 1, \dots, N).
        \label{eq:strang_flow_apex}
    \end{equation}
\end{enumerate}

% ------------------------------------------------------------------------------
\subsubsection{The 7-Operator Hierarchical Palindromic Composition}
\label{sub:strang_7_stage_composition}

Combining these elementary Lie flows into a symmetric palindromic sequence yields the complete second-order normal cone splitting propagator (in the cotangent lift framework of Calaque--Safronov~\cite{calaque2024deformation}).

\begin{defbox}[The 7-Operator Hierarchical Palindromic Strang Propagator]{def:strang_7_stage_propagator}
Let $h > 0$ be the macroscopic time step. The \emph{7-operator hierarchical palindromic Strang propagator} $\mathcal{S}_2(h) : \mathcal{P}_{\mathrm{ext}} \to \mathcal{P}_{\mathrm{ext}}$ is the composite operator:
\begin{equation}
    \mathcal{S}_2(h) = \Phi_{h/2}^{\mathrm{drift}} \circ \Phi_{h/2}^{\mathrm{super}} \circ \Phi_{h/2}^{\mathrm{leaf}} \circ \mathcal{K}_h^{\mathrm{apex}} \circ \Phi_{h/2}^{\mathrm{leaf}} \circ \Phi_{h/2}^{\mathrm{super}} \circ \Phi_{h/2}^{\mathrm{drift}},
    \label{eq:strang_7_stage_definition}
\end{equation}
augmented by the physical time advancement $t \mapsto t + h$, while the conjugate energy variable $E$ is held strictly invariant ($E \mapsto E$).
\end{defbox}

The ordering of subflows in Equation~\eqref{eq:strang_7_stage_definition} reflects the hierarchical nesting of dynamical scales on $\mathbf{Q}[N]$:
\begin{itemize}
    \item The outer subflows $\Phi_{h/2}^{\mathrm{drift}}$ advance the free barycentric and composite coordinates by a half-step $h/2$.
    \item The intermediate subflows $\Phi_{h/2}^{\mathrm{super}}$ and $\Phi_{h/2}^{\mathrm{leaf}}$ evolve the tightly bound internal degrees of freedom under their respective regularized local clocks.
    \item The central subflow $\mathcal{K}_h^{\mathrm{apex}}$ delivers the full screened tidal interaction kick at the midpoint of configuration space, where positions are evaluated with second-order midpoint centering.
\end{itemize}

% ------------------------------------------------------------------------------
\subsubsection{Unconditional Extended Symplecticity and Energy Confinement}
\label{sub:strang_symplecticity_proof}

We now demonstrate that $\mathcal{S}_2(h)$ is an exact symplectomorphism on $(\mathcal{P}_{\mathrm{ext}}, \Omega_{\mathrm{ext}})$ and satisfies strict backward error confinement without dissipative energy drift.

\begin{thmbox}[Unconditional Extended Symplecticity of Strang Splitting]{thm:strang_extended_symplecticity}
In the sense of Strang~\cite{strang1968construction} and Marchuk~\cite{marchuk1968applications}, the 7-operator symmetric hierarchical splitting propagator $\mathcal{S}_2(h) : \mathcal{P}_{\mathrm{ext}} \to \mathcal{P}_{\mathrm{ext}}$ satisfies:
\begin{enumerate}
    \item \textbf{Exact Symplecticity}: $\mathcal{S}_2(h)$ is an exact symplectomorphism of $(\mathcal{P}_{\mathrm{ext}}, \Omega_{\mathrm{ext}})$:
    \begin{equation}
        \mathcal{S}_2(h)^* \Omega_{\mathrm{ext}} = \Omega_{\mathrm{ext}}.
        \label{eq:strang_symplectic_pullback}
    \end{equation}
    \item \textbf{Strict Energy Invariance}: The extended energy variable $E$ is an exact discrete first integral:
    \begin{equation}
        E_{k+1} = E_k = E_0 \quad \forall k \in \mathbb{N}.
        \label{eq:strang_energy_conservation}
    \end{equation}
    \item \textbf{Confinement of Physical Energy}: Across any integration interval $[0, T]$, the physical energy defect satisfies:
    \begin{equation}
        \sup_{k \ge 0, \, k h \le T} \big| H(\mathbf{q}_k, \mathbf{p}_k) - E_0 \big| \le C \, h^2,
        \label{eq:strang_energy_defect_bound}
    \end{equation}
    where $C > 0$ is a constant depending only on the universal Lipschitz constant $L_{\mathrm{apex}}$ and the mass spectrum.
\end{enumerate}
\end{thmbox}
\leanrefright{strang_order_two_theorems}{strang_order_two_symplectic_theorem}{01_StrangOrderTwoSymplectic_Theorems.lean}
\pythonrefright{py_integrators}{integrators/{strang, common}.py}{Float64 / Float128}

\noindent The formal certification badge denotes verification of the \emph{Constructive Combinatorial and 1-Jet Algebraic Kernel} (formalized in Lean~4), certifying exact scalar and 3D tensorial symplectic area preservation under Schwarz symmetry, temporal mirror involution, and palindromic micro-reversibility.

\begin{proofbox}[Proof of Theorem~\ref{thm:strang_extended_symplecticity}]
\textbf{Item (1):} Following the classical geometric splitting theory of Hairer, Lubich and Wanner~\cite[Theorem~II.2.3, pp.~23--27]{hairer2006geometric}, each elementary constituent flow of $\mathcal{S}_2(h)$ is an exact Hamiltonian phase flow on $(\mathcal{P}_{\mathrm{ext}}, \Omega_{\mathrm{ext}})$: the kinematic drifts $\Phi_{\delta t}^{\mathrm{drift}}$ are linear shears with symmetric mass tensor $\mathbb{M}^{-1}$; the leaf Keplerian steps $\Phi_{\vec{\Delta \tau}}^{\mathrm{leaf}}$ are exact analytic Sundman--Stumpff flows on disjoint Poisson-commuting pairs, extended across binary collisions $\mathbf{M}_C$ via McGehee's anti-symplectic boundary bounce $\mathcal{B}_C$ applied at entry and exit; the supernode flows are exact symplectomorphisms on compact sub-clusters (Theorem~\ref{thm:atlas_transition_anti_chattering}); and the apex kick $\mathcal{K}_h^{\mathrm{apex}}$ preserves $\Omega_0$ by Schwarz symmetry of the Hessian of $W_{\mathrm{apex}}^{\mathrm{screen}}$. Since time translation trivially preserves $\mathrm{d}E \wedge \mathrm{d}t$, the composite operator $\mathcal{S}_2(h)$ is an exact extended symplectomorphism: $\mathcal{S}_2(h)^* \Omega_{\mathrm{ext}} = \Omega_{\mathrm{ext}}$.

\textbf{Items (2) and (3):} Throughout all seven stages, no operation modifies $E$, ensuring $E_{k+1} = E_k = E_0$ for all $k$. Backward error analysis~\cite[Chapter~IX, pp.~335--370]{hairer2006geometric} guarantees that the numerical trajectory interpolates the exact flow of a modified Hamiltonian $\widetilde{\mathcal{H}}_{\mathrm{ext}}(h) = \mathcal{H}_{\mathrm{ext}} + h^2 \mathcal{H}_3 + \mathcal{O}(h^4)$, which is an exact first integral of the shadow flow. Gevrey-2 stability (following Marco \& Sauzin~\cite{marco2003stability}) confines the physical energy defect $|H(\mathbf{q}_k, \mathbf{p}_k) - E_0| \le C h^2$ over exponentially long times $T \le \exp(c / h^{1/2})$ without secular energy drift.
\end{proofbox}

% (Covered in §6.4 / Appendix A.25)

% ------------------------------------------------------------------------------
\subsubsection{Baker--Campbell--Hausdorff Parity and Lie--Bargmann Conservation}
\label{sub:strang_bch_parity_bargmann}

A defining property of the symmetric Strang splitting is its micro-reversibility, which eliminates all odd orders of the time step in the formal modified Hamiltonian. Furthermore, following Bargmann~\cite{bargmann1954unitary} and the geometric symmetries established in Proposition~\ref{prop:lie_bargmann_preservation} (Section~\ref{sub:lie_bargmann_and_palindromy}) and Proposition~\ref{prop:lie_bargmann_closure} (Section~\ref{sec:background_conormal}), the scheme preserves all ten Galilean invariants.

\begin{propbox}[Palindromic BCH Parity and Lie--Bargmann Conservation]{prop:strang_bch_parity_bargmann}
Following Bargmann~\cite{bargmann1954unitary} for Galilean central extensions, and Yoshida~\cite{yoshida1990construction} and Hairer, Lubich and Wanner~\cite{hairer2006geometric} for time-reversible palindromic compositions, let $\mathcal{R}_{\mathrm{ext}} : (\mathbf{q}, \mathbf{p}, t, E) \mapsto (\mathbf{q}, -\mathbf{p}, -t, E)$ be the unitary time-reversal involution on $\mathcal{P}_{\mathrm{ext}}$, satisfying $\mathcal{R}_{\mathrm{ext}}^2 = \mathrm{id}$ and anti-symplectic reflection $\mathcal{R}_{\mathrm{ext}}^* \Omega_{\mathrm{ext}} = -\Omega_{\mathrm{ext}}$. Then:
\begin{enumerate}
    \item \textbf{Palindromic Micro-Reversibility}: The Strang propagator $\mathcal{S}_2(h)$ satisfies:
    \begin{equation}
        \mathcal{R}_{\mathrm{ext}} \circ \mathcal{S}_2(h) \circ \mathcal{R}_{\mathrm{ext}} = \mathcal{S}_2(-h) = \mathcal{S}_2(h)^{-1}.
        \label{eq:strang_palindromic_identity}
    \end{equation}
    \item \textbf{Even Shadow Hamiltonian}: The formal Lie generator $Z(h) = \log \mathcal{S}_2(h)$ is strictly odd:
    \begin{equation}
        Z(-h) = -Z(h),
        \label{eq:odd_lie_generator}
    \end{equation}
    so that the formal shadow Hamiltonian $\widetilde{H}(h) = h^{-1} Z(h)$ is strictly even in $h$:
    \begin{equation}
        \widetilde{H}(-h) = \widetilde{H}(h),
        \label{eq:shadow_h_even_box}
    \end{equation}
    eliminating all odd dissipation orders ($\mathcal{O}(h), \mathcal{O}(h^3), \dots$).
    \item \textbf{Exact Lie--Bargmann Invariance}: Under each macro-step $\mathcal{S}_2(h)$, the ten generators of the centrally extended Lie--Bargmann algebra $(\mathbf{P}, \mathbf{L}, \mathbf{G}, \mathcal{H}_{\mathrm{ext}})$ satisfy exact discrete conservation:
    \begin{equation}
        \Delta \mathbf{P} = \mathbf{0}, \qquad \Delta \mathbf{L} = \mathbf{0}, \qquad \Delta \mathbf{G} = \mathbf{0},
        \label{eq:discrete_bargmann_invariance}
    \end{equation}
    and the Galilean boost Poisson bracket compensation $\{\mathbf{G}, \mathcal{H}_{\mathrm{ext}}\} = \mathbf{0}$ holds identically.
\end{enumerate}
\end{propbox}
\leanrefright{lie_bargmann_theorems}{lie_bargmann_conservation_theorem}{02_LieBargmannConservation_Theorems.lean}

\noindent The formal certification badge verifies the algebraic closure of the Lie--Bargmann invariants, 3D vector action-reaction, vanishing total apex torque, and exact Galilean boost step-by-step invariance in Lean~4.

\begin{proofbox}[Proof of Proposition~\ref{prop:strang_bch_parity_bargmann}]
\textbf{Items (1) and (2):} Under the time-reversal involution $\mathcal{R}_{\mathrm{ext}}$, momenta and coordinate time reverse sign ($\mathbf{p} \mapsto -\mathbf{p}, t \mapsto -t$). Since every elementary stage is time-reversible ($\mathcal{R}_{\mathrm{ext}} \circ \Phi_{\delta t} \circ \mathcal{R}_{\mathrm{ext}} = \Phi_{-\delta t}$), conjugating the symmetric composition~\eqref{eq:strang_7_stage_definition} inverts the stage sequence and negates each sub-step, establishing $\mathcal{R}_{\mathrm{ext}} \circ \mathcal{S}_2(h) \circ \mathcal{R}_{\mathrm{ext}} = \mathcal{S}_2(-h) = \mathcal{S}_2(h)^{-1}$. By the Baker--Campbell--Hausdorff theorem for symmetric compositions (Hairer, Lubich and Wanner~\cite[Theorem~V.3.1, pp.~143--150]{hairer2006geometric}), the Lie generator $Z(h) = \log \mathcal{S}_2(h)$ is strictly odd ($Z(-h) = -Z(h)$), which forces the formal modified Hamiltonian $\widetilde{H}(h) = h^{-1} Z(h)$ to be strictly even in $h$, eliminating all odd dissipation orders.

\textbf{Item (3):} Following Bargmann~\cite{bargmann1954unitary}, spatial translation invariance of $W_{\mathrm{apex}}^{\mathrm{screen}}(\mathbf{q})$ enforces $\sum_{i=1}^N \nabla_{\mathbf{q}_i} W_{\mathrm{apex}}^{\mathrm{screen}} = \mathbf{0}$, ensuring exact total linear momentum conservation $\Delta \mathbf{P} = \mathbf{0}$. Rotational $\mathrm{SO}(3)$ invariance ensures that the total apex torque vanishes identically ($\boldsymbol{\tau}_{\mathrm{tot}} = \sum_{i=1}^N \mathbf{q}_i \times \nabla_{\mathbf{q}_i} W_{\mathrm{apex}}^{\mathrm{screen}} \equiv \mathbf{0}$), ensuring $\Delta \mathbf{L} = \mathbf{0}$. Under center-of-mass drift, $\mathbf{R}_{\mathrm{cdm}}(t + \delta t) = \mathbf{R}_{\mathrm{cdm}}(t) + \frac{\delta t}{M_{\mathrm{tot}}} \mathbf{P}$, which exactly cancels the coordinate time shift in $\mathbf{G}(t) = M_{\mathrm{tot}} \mathbf{R}_{\mathrm{cdm}}(t) - t \mathbf{P}$, ensuring $\Delta \mathbf{G} = \mathbf{0}$ across every stage.
\end{proofbox}

\paragraph{Transition to Higher-Order Schemes and the Suzuki Barrier (\S\ref{sec:suzuki_barrier_mcgehee}).}
While the second-order Strang scheme $\mathcal{S}_2(h)$ provides unconditional extended symplecticity, exact palindromic parity, and exact discrete Lie--Bargmann invariance, modern high-precision astrophysics (such as long-term solar system ephemerides and dense star cluster simulations) demands integrators of order $4$ or higher. In classical Hamiltonian integration, higher-order symmetric schemes are constructed via composition:
\begin{equation}
    \mathcal{S}_{2k}(h) = \mathcal{S}_{2k-2}(w_1 h) \circ \dots \circ \mathcal{S}_{2k-2}(w_m h).
    \label{eq:composition_teaser}
\end{equation}
However, the celebrated theorem of Suzuki~\cite{suzuki1991general} and Sheng~\cite{sheng1989solving} establishes that any symplectic composition of order $p \ge 3$ must necessarily include at least one negative sub-step $w_i < 0$. In classical coordinates, a negative time step facing a singular Newtonian collision $\rho_C \to 0$ causes catastrophic failure (diverging inward instead of receding). In Section~\ref{sec:suzuki_barrier_mcgehee}, we show how the desingularized Fulton--MacPherson configuration geometry $\mathbf{Q}[N]$ combined with the unitary McGehee collision bounce involution $\mathcal{B}_C$ effortlessly overcomes the Suzuki barrier.