% ==============================================================================
% chapters/ch06_algorithmic_architecture/sec04_symplectic_subcycle/content.tex
% Section 6.4: The Symplectic Sub-Cycle: Kepler--Sundman Analytic Flow and Screened Apex Kick (Phase 3)
% ==============================================================================

\subsection{The Symplectic Sub-Cycle: Kepler--Sundman Flow and Screened Apex Kick}
\label{sec:symplectic_subcycle}

Each substep $h_k = \Delta t_{\mathrm{macro}} / 2^k$ established by the dyadic scheduler is propagated via a symmetric palindromic composition of Blanes--Moan operators ($\mathrm{BM}_6^{(11)}$ or $\mathrm{BM}_4^{(7)}$). This section formalizes Phase~3: the exact algebraic execution of the kinematic drift, the universal Kepler--Sundman analytic flow on terminal leaves, the synchronous integration of supernodes, and the screened apex $\mathcal{C}^{1,1}$ kick with exact analytical gradients and IEEE~754 hardware safeguards.

% ------------------------------------------------------------------------------
\subsubsection{The Universal Kepler--Sundman Analytic Flow (Phase 3.1)}
\label{sub:kepler_sundman_analytic_flow}

In the Fulton--MacPherson configuration space $\mathbf{Q}[N]$, every terminal leaf node $v \in \operatorname{Leaves}(T)$ corresponds to an isolated binary pair $(i, j) \in \operatorname{LeafPairs}(T)$ whose relative dynamics is governed by the uncoupled two-body Keplerian Hamiltonian:
\begin{equation}
    H_{\mathrm{leaf}}(\mathbf{r}_v, \mathbf{w}_v) = \frac{\|\mathbf{w}_v\|^2}{2\mu_v} - \frac{G m_i m_j}{\|\mathbf{r}_v\|},
    \label{eq:leaf_kepler_hamiltonian}
\end{equation}
where $\mu_v = m_i m_j / (m_i + m_j)$ is the reduced mass, $\mathbf{r}_v = \mathbf{q}_i - \mathbf{q}_j$ is the relative separation vector, and $\mathbf{w}_v = \mu_v \mathbf{v}_{\mathrm{rel}} = \mu_v (\dot{\mathbf{q}}_i - \dot{\mathbf{q}}_j)$ is the conjugate Jacobi momentum.

Rather than subjecting this singular potential to numerical operator splitting (which incurs an unavoidable collision catastrophe in $r_v^{-(2p+1)}$ as diagnosed in Section~\ref{ch:ch05_higher_order_symplectic}), the flow $\exp(h X_{H_{\mathrm{leaf}}})$ is integrated analytically via the universal Sundman--Stumpff diffeomorphism established in Proposition~\ref{prop:sundman_stumpff_diffeomorphism} \cite{sundman1912memoire, stumpff1959himmelsmechanik, stumpff1962himmelsmechanik}. Under the fictitious Sundman time transformation:
\begin{equation}
    \mathrm{d}t = \|\mathbf{r}_v(\tau)\| \, \mathrm{d}\tau = \rho(\tau) \, \mathrm{d}\tau,
    \label{eq:sundman_time_diff}
\end{equation}
the physical trajectory is governed by the universal Sundman time equation:
\begin{equation}
    F(\tau) \coloneqq \rho_0 \tau c_1(z) + \sigma_0 \tau^2 c_2(z) + G M_v \tau^3 c_3(z) - \Delta t = 0,
    \label{eq:sundman_time_equation}
\end{equation}
where $M_v \coloneqq m_i + m_j$, $\rho_0 \coloneqq \|\mathbf{r}_v(0)\|$, $\mathbf{v}_0 \coloneqq \mathbf{w}_v(0)/\mu_v$, $\sigma_0 \coloneqq \mathbf{r}_v(0) \cdot \mathbf{v}_0$, and the energy parameter $\alpha$ and dimensionless Stumpff argument $z$ are:
\begin{equation}
    \alpha \coloneqq \frac{2 G M_v}{\rho_0} - \|\mathbf{v}_0\|^2 = - \frac{2 E_{2b}}{\mu_v}, \qquad z \coloneqq \alpha \tau^2.
    \label{eq:stumpff_z_def}
\end{equation}
The transcendental functions $c_k(z)$ are the universal Stumpff functions \cite{stumpff1959himmelsmechanik, stumpff1962himmelsmechanik} defined by:
\begin{equation}
    c_k(z) \coloneqq \sum_{m=0}^\infty \frac{(-1)^m z^m}{(2m + k)!}, \qquad k \ge 0,
    \label{eq:stumpff_series_def}
\end{equation}
which unify elliptic ($\alpha > 0, z > 0$), parabolic ($\alpha = 0, z = 0$), and hyperbolic ($\alpha < 0, z < 0$) orbits into a single real-analytic manifold without branching singularities.

\medskip
\noindent\textbf{1. Cubic Halley Root Inversion.}
Equation~\eqref{eq:sundman_time_equation} is inverted numerically for $\tau$ via the third-order Halley iteration. Differentiating $F(\tau)$ with respect to $\tau$ using the Stumpff differential identities $c_k'(z) = \frac{1}{2z}(c_{k-1}(z) - k c_k(z))$ yields the exact analytical derivatives:
\begin{align}
    F'(\tau) &= \rho(\tau) = \rho_0 c_0(z) + \sigma_0 \tau c_1(z) + G M_v \tau^2 c_2(z), \label{eq:halley_fprime} \\
    F''(\tau) &= \sigma_0 c_0(z) + (G M_v - \alpha \rho_0) \tau c_1(z), \label{eq:halley_fsecond}
\end{align}
where $G M_v - \alpha \rho_0 = \rho_0 \|\mathbf{v}_0\|^2 - G M_v = \frac{2\rho_0 E_{2b}}{\mu_v} + G M_v \equiv F'''(0) = \rho''(0)$ is the initial Levi-Civita--Sundman curvature invariant from the regularized harmonic oscillator $\rho''(\tau) + \alpha \rho(\tau) = G M_v$. Because $F'(\tau) = \rho(\tau) > 0$ strictly for all non-collisional configurations, $F(\tau)$ is strictly monotonic, ensuring the existence of a unique real root $\tau^*$. The Halley recurrence relation:
\begin{equation}
    \tau_{k+1} = \tau_k - \frac{2 F(\tau_k) F'(\tau_k)}{2 [F'(\tau_k)]^2 - F(\tau_k) F''(\tau_k)},
    \label{eq:halley_recurrence}
\end{equation}
converges cubically, typically attaining relative tolerance $|\Delta t - F(\tau)| / |\Delta t| \le 10^{-15}$ within $2$ to $4$ iterations.

\begin{propbox}[Deterministic Monotone Bisection Fallback]{prop:halley_bisection_fallback}
In the singular vicinity of an exact parabolic grazing encounter ($z \approx 0$) or when denominator degeneration occurs ($2 (F')^2 - F F'' \le 0.1 (F')^2$), the cubic iteration~\eqref{eq:halley_recurrence} is protected by a trust-region clamp $|\delta \tau| \le 0.5 \tau_k$. If convergence is not achieved within $\ell_{\max} = 25$ iterations, the solver falls back to a deterministic monotone bisection search over the guaranteed bracket $[a, b] = [0, \max(\tau_k, |\Delta t|/\rho_0, 1.0)]$.

Because $F'(\tau) > 0$ unconditionally, 60 bisection steps reduce the uncertainty bracket to:
\begin{equation}
    \frac{b - a}{2^{60}} \approx 8.7 \times 10^{-19} (b - a),
    \label{eq:bisection_60_steps}
\end{equation}
guaranteeing convergence to the double-precision machine epsilon without infinite looping or heuristic relaxation.
\end{propbox}

\medskip
\noindent\textbf{2. Universal Lagrange Multipliers and Battin Projection.}
Once $\tau^*$ is determined, the advanced relative position and velocity vectors are evaluated via the universal Lagrange invariants:
\begin{equation}
    \mathbf{r}_v(\Delta t) = f \mathbf{r}_v(0) + g \mathbf{v}_v(0), \qquad
    \mathbf{v}_v(\Delta t) = \dot{f} \mathbf{r}_v(0) + \dot{g} \mathbf{v}_v(0),
    \label{eq:lagrange_fg_propagation}
\end{equation}
where the universal scalar coefficients are:
\begin{align}
    f &= 1 - \frac{G M_v}{\rho_0} (\tau^*)^2 c_2(z^*), &
    g &= \rho_0 \tau^* c_1(z^*) + \sigma_0 (\tau^*)^2 c_2(z^*), \label{eq:lagrange_f_g} \\
    \dot{f} &= - \frac{G M_v}{\rho_0 \rho_{\mathrm{final}}} \tau^* c_1(z^*), &
    \dot{g} &= 1 - \frac{G M_v}{\rho_{\mathrm{final}}} (\tau^*)^2 c_2(z^*), \label{eq:lagrange_fdot_gdot}
\end{align}
with $\rho_{\mathrm{raw}} \coloneqq F'(\tau^*)$ evaluated from~\eqref{eq:halley_fprime}. Near extreme pericenter passages ($\rho_{\mathrm{raw}} < \frac{1}{4}\rho_0$, $1 - e \to 0$), naive summation of~\eqref{eq:halley_fprime} suffers subtractive cancellation of order $(1-e)^{-1}$ in radius and $(1-e)^{-2}$ in total energy; we eliminate this cancellation analytically by evaluating the conditioned terminal radius from the exact angular-momentum/energy quadratic identity (with one fixed-point refinement step substituting $\rho_{\mathrm{final}}$ back into the denominator):
\begin{equation}
    \begin{aligned}
        \rho_{\mathrm{final}} &= \frac{\sigma_{\mathrm{end}}^2 + \|\mathbf{L}_v\|^2 / \mu_v^2}{2 G M_v - \alpha \, \rho_{\mathrm{raw}}}, \\
        \sigma_{\mathrm{end}} &\coloneqq F''(\tau^*) = \sigma_0 c_0(z^*) + (G M_v - \alpha \rho_0) \tau^* c_1(z^*),
    \end{aligned}
    \label{eq:pericenter_quadratic_conditioning}
\end{equation}
coupled with incremental Wronskian renormalization $\delta\dot{g} = (\dot{f}g - \delta f)/f$ ($\delta f = f - 1$, $\dot{g} = 1 + \delta\dot{g}$) whenever $|f| \ge |\dot{g}|$.

\begin{propbox}[Battin Unimodular Projection and Symplectic Preservation~\cite{battin1999introduction}]{prop:battin_unimodular_projection}
The transition matrix $\mathbf{\Phi}_{\mathrm{Kepler}} \in \mathrm{Sp}(6, \mathbb{R})$ satisfies the continuous Wronskian identity:
\begin{equation}
    W(f, g) \coloneqq f \dot{g} - \dot{f} g \equiv 1.
    \label{eq:wronskian_unimodularity}
\end{equation}
To eliminate floating-point subtractive cancellation drift in $W(f, g)$ during extended integrations, the coefficient with the largest numerical sensitivity is re-projected via Battin's algebraic constraint:
\begin{equation}
    \begin{cases}
        \text{if } \max(|f|, |\dot{g}|) \ge \max(|g|, |\dot{f}|): &
        \begin{cases}
            \dot{g}_{\mathrm{proj}} \gets \dfrac{1 + \dot{f} g}{f}, & \text{if } |f| \ge |\dot{g}|, \\
            f_{\mathrm{proj}} \gets \dfrac{1 + \dot{f} g}{\dot{g}}, & \text{if } |\dot{g}| > |f|,
        \end{cases} \\[1.5ex]
        \text{if } \max(|g|, |\dot{f}|) > \max(|f|, |\dot{g}|): &
        \begin{cases}
            \dot{f}_{\mathrm{proj}} \gets \dfrac{f \dot{g} - 1}{g}, & \text{if } |g| \ge |\dot{f}|, \\
            g_{\mathrm{proj}} \gets \dfrac{f \dot{g} - 1}{\dot{f}}, & \text{if } |\dot{f}| > |g|,
        \end{cases}
    \end{cases}
    \label{eq:battin_projection}
\end{equation}
guaranteeing $|f \dot{g} - \dot{f} g - 1| \le 10^{-16}$ at every substep. The updated momentum is $\mathbf{w}_v(\Delta t) = \mu_v \mathbf{v}_v(\Delta t)$, preserving the exact Darboux symplectic 2-form $\mathrm{d}\mathbf{r}_v \wedge \mathrm{d}\mathbf{w}_v$.
\end{propbox}

\medskip
\noindent\textbf{3. Radial Head-On Collisions and McGehee Unitary Bounce.}
For degenerate rectilinear trajectories with vanishing specific angular momentum:
\begin{equation}
    \|\mathbf{L}_v\|^2 = \|\mathbf{r}_v \times \mathbf{v}_v\|^2 < 10^{-24} \max\left( G M_v \rho_0 + \|\mathbf{v}_0\|^2 \rho_0^2, \; 10^{-30} \right),
    \label{eq:radial_collision_threshold}
\end{equation}
the orbital plane collapses onto a 1D ray. In this regime, the system undergoes an exact collision-ejection bounce modeled by the McGehee unitary reflection operator~\cite{mcgehee1974triple}:
\begin{equation}
    \mathcal{B}_C : (r_{\mathrm{scalar}}, w_r) \mapsto (-r_{\mathrm{scalar}}, -w_r), \qquad \mathcal{B}_C^2 = \mathrm{Id},
    \label{eq:mcgehee_bounce_unitary}
\end{equation}
acting as a specular reflection of the radial momentum along the invariant directional ray $\hat{\mathbf{s}}_0 = \mathbf{r}_0 / \rho_0$, as established in Section~\ref{ch:ch05_higher_order_symplectic} (Lemma~\ref{thm:mcgehee_bounce_unitary}) and certified in Lean~4.

Because the Sundman radius $\rho(\tau) = F'(\tau) \ge 0$ attains a quadratic minimum $\rho(\tau_{\mathrm{coll}}) = 0$ at collision and remains non-negative for $\tau > \tau_{\mathrm{coll}}$ while the radial virial $\sigma_{\mathrm{end}}(\tau) = F''(\tau) = \rho'(\tau)$ smoothly reverses sign from $< 0$ (infall) to $> 0$ (rebound) across $\tau_{\mathrm{coll}}$, using the conditioned terminal radius $\rho_{\mathrm{final}} \ge 0$ from~\eqref{eq:pericenter_quadratic_conditioning} directly reconstitutes the post-collision state symplectically along the incident ray $\hat{\mathbf{s}}_0$ via:
\[
    \mathbf{r}_{\mathrm{final}} = \rho_{\mathrm{final}} \, \hat{\mathbf{s}}_0, \qquad
    \mathbf{v}_{\mathrm{final}} = \left(\frac{\sigma_{\mathrm{end}}}{\rho_{\mathrm{final}}}\right) \hat{\mathbf{s}}_0, \qquad
    \mathbf{w}_{\mathrm{final}} = \mu_v \mathbf{v}_{\mathrm{final}}
\]
(with asymptotic limit $v_{\mathrm{radial}} = c_1(z) / (c_2(z)\tau)$ as $\rho_{\mathrm{final}} \to 0$), without any sheet-sign branching or non-symplectic velocity-norm projection.

\begin{rembox}[IEEE 754 Safeguards: $T$-Involution, Hyperbolic Guard, and Horner Expansion]\label{rem:ieee754_kepler_safeguards}
Three critical machine-precision edge cases in universal Kepler--Sundman propagation are resolved as follows:
\begin{enumerate}
    \item \textbf{Microscopic Time-Reversibility under Retrograde Steps} ($\Delta t < 0$): High-order symplectic integrators ($\mathrm{BM}_4^{(7)}, \mathrm{BM}_6^{(11)}$) require negative intermediate stage weights ($a_i < 0, b_i < 0$) to bypass the Suzuki order barrier. Directly supplying $\Delta t < 0$ to root-finding routines can bias the convergence basin. Exact reversibility $\Psi_{-\Delta t} = \Psi_{\Delta t}^{-1}$ is enforced via the kinematic time-inversion involution:
    \begin{equation}
        T : (\mathbf{r}, \mathbf{w}) \mapsto (\mathbf{r}, -\mathbf{w}), \qquad
        \Phi_{\Delta t}(\mathbf{r}, \mathbf{w}) = T \circ \Phi_{|\Delta t|}(\mathbf{r}, - \mathbf{w}) \circ T,
        \label{eq:retrograde_involution}
    \end{equation}
    which maps backward motion onto forward propagation with inverted initial momentum $\sigma_0^{\mathrm{eff}} = \operatorname{sgn}(\Delta t) \sigma_0$.
    \item \textbf{Hardware Hyperbolic Guard Ceiling} ($Z_{\mathrm{MIN}} = -500\,000.0$): In deep hyperbolic plunging encounters ($E_{2b} > 0, \alpha < 0$), the transcendental argument $z = \alpha \tau^2 \ll 0$ requires evaluating $\sinh(\sqrt{-z}) / \sqrt{-z}$. In IEEE~754 double precision, $\sinh(x)$ overflows for $x > 709.7827$, corresponding to $-z \ge 503\,791$. Setting:
    \begin{equation}
        Z_{\mathrm{MIN}} \coloneqq -500\,000.0 \implies \cosh(\sqrt{500\,000}) \approx 1.13 \times 10^{307} < 1.797 \times 10^{308},
        \label{eq:z_min_hyperbolic}
    \end{equation}
    provides a hardware safety ceiling preventing floating-point overflow during initial iteration overshoot, while preserving strict zero softening ($\varepsilon = 0, d_0 = 0$).
    \item \textbf{Horner Polynomial Expansion for Small Arguments} ($|z| < 0.1$): For near-parabolic regimes, direct evaluation of trigonometric formulas suffers catastrophic subtractive cancellation. The Stumpff functions are evaluated via a 6th-order Horner expansion:
    \begin{equation}
        \begin{aligned}
            c_k(z) &= \frac{1}{k!} - \frac{z}{(k+2)!} + \frac{z^2}{(k+4)!} - \frac{z^3}{(k+6)!} \\
            &\quad + \frac{z^4}{(k+8)!} - \frac{z^5}{(k+10)!} + \frac{z^6}{(k+12)!} + \mathcal{O}(z^7),
        \end{aligned}
        \label{eq:horner_stumpff}
    \end{equation}
    yielding truncation errors strictly bounded by $10^{-18}$, well beneath double-precision floor.
\end{enumerate}
\end{rembox}

% ------------------------------------------------------------------------------
\subsubsection{Synchronous Supernodes and Conformal Sundman Time (Phase 3.2)}
\label{sub:synchronous_supernodes_flow}

When a compact multi-body cluster ($k \ge 3$) enters the canonical Schmitt trigger hysteresis band:
\begin{equation}
    \eta_S \coloneqq \frac{\bar{\rho}_S}{r_P} \left( \frac{M_{\mathrm{ext}}}{M_S} \right)^{1/3} \le 2^{-1/3},
    \label{eq:supernode_activation_condition}
\end{equation}
the hierarchical pairwise tree structure experiences the \emph{Excluded Middle Dilemma} (Section~\ref{ch:ch04_fulton_macpherson_splitting}, Section~\ref{sec:synchronous_supernodes}): asynchronous dyadic sub-stepping within the cluster induces an uncompensable Poisson bracket divergence in $\mathcal{O}(\bar{\rho}_S^{-5/2})$.

To suppress this divergence, the active cluster $S \in \mathcal{S}_{\mathrm{sync}}$ is treated as an indivisible multi-body supernode, whose internal degrees of freedom are integrated synchronously under the unfolded conformal Sundman clock:
\begin{equation}
    \mathrm{d}t_S = \chi_S^2 \, \mathrm{d}\tau_S = \bar{\rho}_S \, \mathrm{d}\tau_S, \qquad \chi_S \coloneqq \sqrt{\bar{\rho}_S} > 0.
    \label{eq:unfolded_supernode_clock}
\end{equation}
Under this conformal time rescaling, binary collisions within $S$ correspond to finite traversal times in $\tau_S$. By Saari's Singularity Theorem (Theorem~\ref{thm:anti_zeno_regularization}), the minimum hyper-radius $\bar{\rho}_{\min} > 0$ is strictly positive almost everywhere on $T^* \mathbf{Q}[N]$, bounding the number of required internal substeps:
\begin{equation}
    N_{\mathrm{sub}}(S) \le \left\lceil \frac{\Delta t_n}{\bar{\rho}_{\min} \, \Delta \tau_S} \right\rceil < \infty,
    \label{eq:anti_zeno_bound}
\end{equation}
where the mass-weighted Jacobi hyper-radius $\rho_S$ of Section~\ref{sec:synchronous_supernodes} (Eq.~\eqref{eq:unfolded_sundman_clock}), the mass-normalized hyper-radius $\bar{\rho}_S = \rho_S / \sqrt{M_S}$, and the unweighted barycentric Euclidean norm $\|\mathbf{x}_S\|_2 = \bigl(\sum_{i \in S} \|\mathbf{q}_i - \mathbf{R}_S\|^2\bigr)^{1/2}$ satisfy the uniform norm equivalence:
\[
    \sqrt{\frac{m_{\min}(S)}{M_S}}\,\|\mathbf{x}_S\|_2 \;\le\; \bar{\rho}_S \;\le\; \sqrt{\frac{m_{\max}(S)}{M_S}}\,\|\mathbf{x}_S\|_2.
\]
Operationalizing~\eqref{eq:anti_zeno_bound}, \texttt{\_compute\_local\_substep\_bound} in \texttt{solvers/supernode.py} evaluates the time-reversal-even internal step bound:
\[
    h_{\mathrm{bound}}(\mathbf{x}_S, \mathbf{v}_S) = 0.05 \min\!\left(2\pi\sqrt{\frac{r_{\mathrm{second}}^3}{G M_S}},\; \frac{r_{\mathrm{second}}}{v_{\max}},\; \frac{\|\mathbf{x}_S\|_2}{v_{\max}}\right),
\]
where $r_{\mathrm{second}}$ denotes the second-smallest pairwise separation outside the analytically desingularized closest pair $(i^*, j^*)$ and $v_{\max} = \max_{i \in S} \|\mathbf{v}_i - \mathbf{V}_S\|$.

Throughout the synchronous supernode sub-cycle, the internal Hamiltonian $H_S^{\mathrm{int}} = T_S + V_S$ is advanced in conformal Sundman time $\tau_S$ through an \emph{adaptive symplectic inheritance}: rather than prescribing a static scheme (such as the 4th-order Yoshida composition $\mathrm{Y}_4$, whose leading error constant is over $400\times$ larger than that of Blanes--Moan schemes), the internal solver dynamically inherits the palindromic stage coefficients $(a_s, b_s)$ of the active parent integrator ($\mathrm{BM}_4^{(7)}$ with 7 drifts and 8 kicks, $\mathrm{BM}_6^{(11)}$ with 11 drifts and 12 kicks, or Strang $\mathrm{S}_2$). This eliminates order mismatch across the hierarchical interface, ensuring that the local shadow Hamiltonian preserves the exact formal order ($2p$) and Lie bracket error structure of the global integrator.

Crucially, the normal cone splitting is applied \emph{recursively} inside the supernode itself. In conformal Sundman time $\tau_S$, the internal degrees of freedom are partitioned into:
\begin{equation}
    H_S^{\mathrm{int}} = H_{\mathrm{leaf}}^{\mathrm{int}} + H_{\mathrm{base}}^{\mathrm{int}} + V_{\mathrm{apex}}^{\mathrm{int}},
    \label{eq:supernode_conormal_split}
\end{equation}
where:
\begin{enumerate}
    \item $H_{\mathrm{leaf}}^{\mathrm{int}}$ desingularizes the tightest internal binary pair $(i^*, j^*) \coloneqq \operatorname{argmin}_{i < j \in S} \|\mathbf{q}_i - \mathbf{q}_j\|$ via the exact analytical Kepler--Sundman flow (governed by Proposition~\ref{prop:battin_unimodular_projection}, with Battin unimodular projection~\eqref{eq:battin_projection} and McGehee radial bounce~\eqref{eq:mcgehee_bounce_unitary} upon near-rectilinear collapse);
    \item $H_{\mathrm{base}}^{\mathrm{int}}$ advances the remaining spectator bodies $k \in S \setminus \{i^*, j^*\}$ through uncoupled linear kinetic drift;
    \item $V_{\mathrm{apex}}^{\mathrm{int}} \coloneqq V_S - V_{\mathrm{Kepler}}(i^*, j^*)$ represents the non-singular internal cross-interaction potential, applied as momentum kicks along the internal cotangent fibers.
\end{enumerate}
This internal splitting is composed through the inherited stages $(a_s, b_s)$, guaranteeing that binary close encounters occurring within a chaotic multi-body encounter are analytically resolved without step collapse or energy divergence. This architectural property embodies the exact fractal self-similarity of the normal cone desingularization across all hierarchical scales of $\mathbf{Q}[N]$.

Throughout this evolution, the internal center of mass remains decoupled from internal gravitational stresses:
\begin{equation}
    \sum_{i \in I_S} \mathbf{F}_i^{\mathrm{int}} \equiv \mathbf{0}, \qquad
    \sum_{i \in I_S} \mathbf{q}_i \times \mathbf{F}_i^{\mathrm{int}} \equiv \mathbf{0},
    \label{eq:supernode_conservation}
\end{equation}
ensuring exact discrete preservation of internal linear and angular momentum, with floating-point roundoff bounded by IEEE~754 precision limits ($< 10^{-15}$ and $< 10^{-14}$ respectively).

% ------------------------------------------------------------------------------
\subsubsection{Screened Apex Kick and Exact Analytical Gradient (Phase 3.3)}
\label{sub:screened_apex_kick_flow}

The non-integrable cross-cluster interactions in the Fulton--MacPherson splitting are governed by the hierarchical screened apex potential $W_{\mathrm{apex}}^{\mathrm{screen}} \in \mathcal{C}^{1,1}(\mathbf{Q}[N])$, whose universal Lipschitz regularity is demonstrated in Theorem~\ref{thm:universal_c11_apex} and whose Galilean invariance is certified by Proposition~\ref{prop:lie_bargmann_preservation}. In this subsection, we formalize the exact analytical evaluation of the symplectic kick operator, which acts as an exact translation along the vertical cotangent fibers $T^*_{\mathbf{q}}\mathbf{Q}[N]$ leaving base configuration coordinates on $\mathbf{Q}[N]$ unperturbed:
\begin{equation}
    \exp\left( b_s h_k X_{W_{\mathrm{apex}}^{\mathrm{screen}}} \right) : \mathbf{p}_i \mapsto \mathbf{p}_i + b_s h_k \mathbf{F}_i^{\mathrm{apex}}, \qquad 1 \le i \le N.
    \label{eq:apex_kick_operator}
\end{equation}

\medskip
\noindent\textbf{1. Scale-Adapted Screening Cutoff and Closed Cauchy Mollifier.}
For every internal composite node $v \in \operatorname{Internal}(T) \setminus \operatorname{Leaves}(T)$, the dimensionless scale ratio is:
\begin{equation}
    u_v \coloneqq \frac{\|\mathbf{r}_v\|}{\bar{\rho}_{L(v)} + \bar{\rho}_{R(v)}},
    \label{eq:dimensionless_u_v}
\end{equation}
where $\mathbf{r}_v \coloneqq \mathbf{X}_{L(v)} - \mathbf{X}_{R(v)}$ is the vector separating the barycenters of the left and right child clusters, and $\bar{\rho}_C \coloneqq \rho_C / \sqrt{M_C}$ are their respective mass-normalized hyper-radii.

The scale-adapted screening function $\Theta(u_v)$ smooths the transition between the compact near-field core and the asymptotic multipolar regime:
\begin{equation}
    \Theta(u) \coloneqq \begin{cases}
        0, & \text{if } u \le 1, \\
        \sigma(u - 1), & \text{if } 1 < u < 2, \\
        1, & \text{if } u \ge 2,
    \end{cases}
    \label{eq:screening_theta_regimes}
\end{equation}
where $\sigma(s)$ is a smooth monotonic transition on $s \in [0, 1]$. To avoid numerical quadratures of transcendental functions while guaranteeing sufficient smoothness for high-order splitting, the symplectic integrator supports either the smooth partition-of-unity Cauchy mollifier:
\begin{equation}
    \sigma_{\mathrm{Cauchy}}(s) \coloneqq \frac{g_1(s)}{g_1(s) + g_2(s)}, \quad
    g_1(s) \coloneqq \exp\left(-\frac{1}{s}\right), \quad
    g_2(s) \coloneqq \exp\left(-\frac{1}{1-s}\right),
    \label{eq:closed_cauchy_mollifier}
\end{equation}
or the exact order-matched Hermite polynomial spline family $\sigma_{C^{2p}}(s)$ ($\sigma_{C^2}(s) = 10s^3 - 15s^4 + 6s^5 \in \mathcal{C}^{2,1}$ for Strang order $2$, $\sigma_{C^4}(s)$ of degree $9$ for $\mathrm{Y}_4/\mathrm{BM}_4^{(7)}$, $\sigma_{C^6}(s)$ of degree $13$ for $\mathrm{BM}_6^{(11)}/\mathrm{PBM}_4$, and centered $\sigma_{C^8}(s)$ of degree $17$) whose derivatives up to order $2p$ vanish identically at $s = 0$ and $s = 1$.

\medskip
\noindent\textbf{2. Analytical Monopole Cancellation.}
In the mathematical definition of the blown-up Hamiltonian on $T^*\mathbf{Q}[N]$ (Section~\ref{ch:ch04_fulton_macpherson_splitting}, Eq.~\eqref{eq:total_screened_apex_def}), the regularized apex potential subtracts the screened monopole:
\begin{equation}
    W_{\mathrm{apex}}^{\mathrm{screen}} = V_{\mathrm{Newton}} - \sum_{(i, j) \in \mathcal{L}_{\mathrm{pairs}}} V_{\mathrm{Kepler}}(i, j) - \sum_{v \in \mathrm{Comp}} \Theta(u_v) \frac{G M_L M_R}{\|\mathbf{r}_v\|} - \sum_{S \in \mathcal{S}_{\mathrm{sync}}} V_{\mathrm{int}}(S).
    \label{eq:v_apex_hier_definition}
\end{equation}
Because extended composite clusters lack closed-form Keplerian solutions, their screened interactions cannot be evolved analytically in $H_0$; they are integrated numerically within the kick stage alongside $W_{\mathrm{apex}}^{\mathrm{screen}}$. In the operator-splitting kick step, recombining the regularizing monopole $+W_v$ with $W_{\mathrm{apex}}^{\mathrm{screen}}$ produces an exact analytical cancellation ($+W_v - W_v \equiv 0$), so that the effective kick potential reduces identically to the exact inter-cluster Newtonian potential $V_{\mathrm{inter}}$:
\begin{equation}
    V_{\mathrm{kick}} \coloneqq W_{\mathrm{apex}}^{\mathrm{screen}} + V_{\mathrm{screened}}^{\mathrm{monopole}} \equiv V_{\mathrm{inter}},
    \label{eq:monopole_cancellation_identity}
\end{equation}
where $V_{\mathrm{inter}}$ is the exact inter-cluster Newtonian potential:
\begin{equation}
    V_{\mathrm{inter}} = - \sum_{\substack{1 \le i < j \le N \\ (i, j) \notin \mathcal{L}_{\mathrm{pairs}} \\ (i, j) \notin \mathcal{S}_{\mathrm{sync}}}} \frac{G m_i m_j}{\|\mathbf{q}_i - \mathbf{q}_j\|}.
    \label{eq:v_inter_exact}
\end{equation}
The terms $+ \Theta(u_v) \frac{G M_L M_R}{\|\mathbf{r}_v\|}$ and $- \Theta(u_v) \frac{G M_L M_R}{\|\mathbf{r}_v\|}$ cancel analytically to machine zero ($+W_v - W_v \equiv 0$). Evaluating $V_{\mathrm{inter}}$ directly eliminates virtual fictitious monopole kicks and uncompensated reaction forces.

\medskip
\noindent\textbf{3. Coupled Hyper-Radius Deformation Gradient.}
When the screened composite potential $W_v = \Theta(u_v) \frac{G M_L M_R}{\|\mathbf{r}_v\|}$ is retained in theoretical diagnostics, its exact Cartesian gradient must account not only for the displacement of cluster barycenters $\mathbf{X}_L, \mathbf{X}_R$, but also for internal hyper-radius deformation:
\begin{equation}
    \nabla_{\mathbf{q}_i} W_v = \frac{\partial W_v}{\partial \mathbf{r}_v} \cdot \frac{\partial \mathbf{r}_v}{\partial \mathbf{q}_i} + \frac{\partial W_v}{\partial R_{\mathrm{scale}}} \cdot \frac{\partial R_{\mathrm{scale}}}{\partial \bar{\rho}_L} \cdot \nabla_{\mathbf{q}_i} \bar{\rho}_L, \qquad i \in \operatorname{bodies}(L(v)),
    \label{eq:coupled_hyperradius_gradient}
\end{equation}
where $R_{\mathrm{scale}} = \bar{\rho}_L + \bar{\rho}_R$ and:
\begin{equation}
    \nabla_{\mathbf{q}_i} \bar{\rho}_L = \frac{m_i (\mathbf{q}_i - \mathbf{X}_L)}{M_L \bar{\rho}_L}.
    \label{eq:nabla_rho_bar_formula}
\end{equation}

\begin{propbox}[Exact Galilean Invariance of the Apex Kick]{prop:apex_galilean_invariance}
Let $\mathbf{F}_i^{\mathrm{apex}} \coloneqq - \nabla_{\mathbf{q}_i} V_{\mathrm{kick}}$ be the net Cartesian force acting on particle $i$. Because $V_{\mathrm{kick}}$ is an invariant scalar function under the Euclidean group $\mathrm{SE}(3)$ (translations and rotations of $\mathbb{R}^3$), the resulting forces satisfy the exact discrete balance identities:
\begin{equation}
    \sum_{i=1}^N \mathbf{F}_i^{\mathrm{apex}} \equiv \mathbf{0}, \qquad
    \sum_{i=1}^N \mathbf{q}_i \times \mathbf{F}_i^{\mathrm{apex}} \equiv \mathbf{0},
    \label{eq:apex_force_conservation}
\end{equation}
with numerical evaluation under IEEE~754 floating-point arithmetic preserving total linear and angular momentum to machine precision floor ($< 10^{-15}$ and $< 10^{-14}$ respectively).
\end{propbox}

% ------------------------------------------------------------------------------
\subsubsection{Palindromic Multi-Stage Splitting Architectures ($\mathrm{BM}_4^{(7)}$ and $\mathrm{BM}_6^{(11)}$)}
\label{sub:blanes_moan_splitting_architectures}

The elementary operators $\Phi_h^{\mathrm{drift}}$, $\Phi_h^{\mathrm{leaf}}$, $\Phi_h^{\mathrm{super}}$, and $\Phi_h^{\mathrm{apex}}$ are assembled into a single, high-order, symmetric substep propagator via the processed Blanes--Moan splitting families~\cite{blanes2002practical}:
\begin{equation}
    \Psi_h = \exp\left( b_s h X_{W_{\mathrm{apex}}^{\mathrm{screen}}} \right) \exp\left( a_{s-1} h X_{H_0} \right) \cdots \exp\left( a_1 h X_{H_0} \right) \exp\left( b_1 h X_{W_{\mathrm{apex}}^{\mathrm{screen}}} \right),
    \label{eq:blanes_moan_general_composition}
\end{equation}
where $H_0 \coloneqq H_{\mathrm{base}} + H_{\mathrm{leaf}} + H_{\mathrm{super}}^{\mathrm{sync}}$.
\begin{enumerate}
    \item \textbf{Scheme $\mathrm{BM}_4^{(7)}$} (7 drift stages, 8 kick stages, pre-asymptotic effective order $6$ on $\mathcal{D}_{\mathrm{reg}}$ for $\varepsilon \ll h^2$ under boundary conjugator $\mathcal{Z}_4^{\mathrm{BM4}}$ generated by $K_4 \in \mathcal{P}_1$): Proven in Theorem~\ref{thm:bm4_error_minimization} and Theorem~\ref{thm:c4_corrector_solvability} to minimize the leading Lie error commutator norm on $\mathbf{Q}[N]$ while bounding secular shadow energy dissipation. Requires 8 force evaluations per standalone substep under generic 7-stage Strang composition $\mathcal{S}_2^*(a_j h)$. By virtue of palindromic symmetry ($a_i = a_{s-i}, b_i = b_{s+1-i}$ with $b_s = b_1$), the first and last kick stages between consecutive substeps coalesce via the First-Same-As-Last (FSAL) property:
    \begin{equation}
        \exp(b_s h X_V) \circ \exp(b_1 h X_V) = \exp((b_s + b_1) h X_V),
        \label{eq:fsal_coalescence}
    \end{equation}
    reducing the amortized evaluation cost to $m = 7$ force evaluations per chained $\mathcal{S}_2^*(a_j h)$ step (and $s = 6$ effective non-trivial evaluations under direct ABA/RKN stage merging~\cite[Tables~2--3]{blanes2002practical}, Definition~\ref{def:bm4_seven_stages}).
    \item \textbf{Scheme $\mathrm{BM}_6^{(11)}$} (11 drift stages, 12 kick stages, unconditional direct order $6$ without boundary corrector): Proven in Theorem~\ref{thm:bm6_order_conditions} and Proposition~\ref{prop:bm6_secular_phase_kpis} to solve the complete 6th-order RKN Lie algebraic system on $\mathbf{Q}[N]$, bypassing the Universal Obstruction Theorem~\ref{thm:universal_order6_obstruction} at the cost of 12 kick stages per standalone substep (reduced to $m = 11$ evaluations under continuous FSAL chaining, or $s = 10$ effective evaluations under direct ABA/RKN stage merging). Provides an unconditional 6th-order integrator across all dynamical regimes without the $\mathcal{O}(\varepsilon^2 h^4)$ quadratic floor.
\end{enumerate}

\medskip
\noindent\textbf{Vectorized and Zero-Allocation CPU Sub-Cycle Execution.}
In contrast to GPU SIMD architectures (which suffer acute warp divergence when different leaf pairs require unequal Halley iteration counts $\ell_v \in \{2, 3, 4\}$ or branch into the McGehee bounce), the execution of Phase~3 on CPUs exploits the geometric decoupling of the phase space:
\begin{itemize}
    \item All terminal leaf pairs $(i, j) \in \operatorname{LeafPairs}(T)$ are disjoint in particle indices; their Kepler--Sundman flows commute exactly and operate on independent coordinate slices.
    \item The inter-cluster apex kick matrix $\mathbf{F}_{ij}$ is accumulated via vectorized NumPy reductions over pre-allocated buffers.
    \item The global phase state update operates with an $\mathcal{O}(N)$ memory footprint ($< 250$~kB for $N \le 1000$, residing entirely in L2/L3 CPU cache).
\end{itemize}

\pythonrefright{py_kernels}{solvers/{kepler_sundman, supernode}.py, screened_apex.py, common.py}{F64 / F128}