% --- Section Content : The Sheng--Suzuki Barrier and McGehee's Involutive Collision Bounce ---
\subsection{The Sheng--Suzuki Barrier and McGehee's Involutive Collision Bounce}
\label{sec:suzuki_barrier_mcgehee}

The second-order palindromic Strang scheme $\mathcal{S}_2(h)$ established in~\S\ref{sec:strang_propagator} provides an unconditional, structure-preserving foundation on the compactified phase space $\mathcal{P}_{\mathrm{ext}} = T^*\mathbf{Q}[N] \times \mathbb{R}_t \times \mathbb{R}_E$. However, in precision celestial mechanics (such as secular planetary ephemerides over multi-millennial horizons or dense globular cluster cores), second-order integrators exhibit an intrinsic limitation: while energy and first integrals remain bounded without secular drift, the accumulated orbital phase error grows as $\Delta \theta(t) \sim \mathcal{O}(h^2 \cdot t)$. Eliminating this secular phase error demands higher-order geometric integrators ($p \ge 4$).

In classical geometric numerical integration, higher-order symmetric integrators are generated by composing symmetric second-order base propagators $\mathcal{S}_2$. Yet, the celebrated Sheng--Suzuki impossibility theorem~\cite{sheng1989solving, suzuki1990fractal, suzuki1991general} establishes that any real-coefficient splitting or composition scheme of order $p \ge 3$ must contain at least one strictly negative time sub-step ($w_i < 0$). In standard gravitational $N$-body codes, backward time steps are fatal: near an unsoftened collision $\rho_C \to 0$, integrating backward accelerates approaching bodies into singular division-by-zero or numerical deadlock. 

In this section, we show that the Fulton--MacPherson blow-up desingularization on $\mathbf{Q}[N]$ effortlessly surmounts the Suzuki barrier: the universal Sundman--Stumpff two-body flow is an exact, real-analytic diffeomorphism on the whole real line $\mathbb{R}$, while McGehee's collision boundary bounce operator is a unitary involution ($\mathcal{B}_C^2 = \operatorname{Id}$) guaranteeing exact, time-reversible collision traversal during negative sub-steps. We formalize this resolution through Yoshida's canonical order-$4$ symmetric composition.

% ------------------------------------------------------------------------------
\subsubsection{The Sheng--Suzuki Obstruction to Positive-Coefficient Splittings}
\label{sub:suzuki_obstruction}

To elevate a symmetric second-order integrator $\mathcal{S}_2(h)$ to order $p \ge 4$, one considers an $s$-stage palindromic composition:
\begin{equation}
    \mathcal{S}_p(h) = \mathcal{S}_2(w_s h) \circ \mathcal{S}_2(w_{s-1} h) \circ \dots \circ \mathcal{S}_2(w_1 h), \qquad \sum_{i=1}^s w_i = 1,
    \label{eq:general_symmetric_composition}
\end{equation}
with symmetric weights $w_{s+1-i} = w_i$. The formal Baker--Campbell--Hausdorff expansion of $\mathcal{S}_p(h)$ expresses the modified Lie generator $Z_p(h) = \log \mathcal{S}_p(h)$ as a series of nested Lie brackets. Achieving order $p \ge 3$ requires eliminating the lower-order commutators. The fundamental algebraic barrier governing this elimination is formulated as follows:

\begin{thmbox}[The Sheng--Suzuki Impossibility Barrier]{thm:sheng_suzuki_barrier}
Following Sheng~\cite{sheng1989solving} and Suzuki~\cite{suzuki1990fractal, suzuki1991general}, let $\mathcal{A}$ and $\mathcal{B}$ be non-commuting operators in a free Lie algebra. Any real-coefficient exponential splitting:
\begin{equation}
    \exp\bigl(h(\mathcal{A} + \mathcal{B})\bigr) = \prod_{i=1}^m \exp(c_i h \mathcal{A}) \exp(d_i h \mathcal{B}) + \mathcal{O}(h^{p+1})
    \label{eq:suzuki_splitting_form}
\end{equation}
or any symmetric composition of second-order schemes:
\begin{equation}
    \mathcal{S}_p(h) = \prod_{i=1}^s \mathcal{S}_2(w_i h) + \mathcal{O}(h^{p+1})
    \label{eq:suzuki_composition_form}
\end{equation}
achieving convergence order $p \ge 3$ must necessarily possess at least one strictly negative stage coefficient:
\begin{equation}
    \exists i \in \{1, \dots, s\} \quad \text{such that} \quad w_i < 0.
    \label{eq:negative_step_necessity}
\end{equation}
\end{thmbox}

\begin{proofbox}[Proof of Theorem~\ref{thm:sheng_suzuki_barrier}]
For any symmetric composition $\mathcal{S}_p(h) = \prod_{i=1}^s \mathcal{S}_2(w_i h)$, eliminating the order-3 Baker--Campbell--Hausdorff error term requires $\sum_{i=1}^s w_i = 1$ and $\sum_{i=1}^s w_i^3 = 0$ (Hairer, Lubich and Wanner~\cite[Section~II.4, p.~46]{hairer2006geometric}). If all $w_i \ge 0$, then $\sum w_i^3 > 0$, an immediate contradiction requiring $\min_i w_i < 0$. The general impossibility for arbitrary exponential splittings was established in full generality by Sheng~\cite{sheng1989solving} and Suzuki~\cite{suzuki1990fractal, suzuki1991general} through free Lie algebra commutator bases.
\end{proofbox}

\begin{rembox}[Scope of the Sheng--Suzuki Barrier and Normal-Cone Desingularization]\label{rem:sheng_suzuki_scope}
It is essential to circumscribe the exact mathematical domain of Theorem~\ref{thm:sheng_suzuki_barrier}: the necessity of negative sub-steps ($w_i < 0$ or $c_i < 0$) is an algebraic consequence of operating in a \emph{free Lie algebra} $\mathcal{L}(\mathcal{A}, \mathcal{B})$. In that unconstrained setting, the two operators $\mathcal{A}$ and $\mathcal{B}$ possess comparable magnitudes ($\varepsilon \sim 1$), and all independent Lie brackets $[\mathcal{A}, [\mathcal{A}, \mathcal{B}]]$, $[\mathcal{B}, [\mathcal{B}, \mathcal{A}]]$, etc., must be eliminated simultaneously without scale differentiation.

In general-purpose celestial mechanics compositions where no asymptotic hierarchy is presupposed ($\varepsilon \sim 1$), such as the 3-stage Yoshida propagator $\mathcal{S}_4$ below and the optimal 7-stage Blanes--Moan scheme $\mathrm{BM}_4$ studied in~\S\ref{sec:bm4_seven_stages}, negative sub-steps remain algebraically unavoidable. It is precisely for these schemes that the Fulton--MacPherson desingularization provides an unconditional guarantee: backward time steps are rendered completely singularity-free through the unitary McGehee collision bounce involution ($\mathcal{B}_C^2 = \operatorname{Id}_{\mathbf{M}_C}$), as established below.
\end{rembox}

In classical celestial mechanics, Theorem~\ref{thm:sheng_suzuki_barrier} represents an insurmountable barrier for standard unsoftened Cartesian codes. In standard Newtonian coordinates, the gravitational force $\mathbf{F}_{ij} = -G m_i m_j (\mathbf{q}_i - \mathbf{q}_j) / \|\mathbf{q}_i - \mathbf{q}_j\|^3$ diverges as $\rho_{ij} \to 0$. A negative sub-step $\delta t = w_0 h < 0$ inverts the velocity vector $\dot{\mathbf{q}} \mapsto -\dot{\mathbf{q}}$ with respect to the forward marching direction, driving approaching bodies backward along their orbit toward the singularity. In adaptive step-size codes, this backward motion produces step-size collapse, numerical stiffness, and loss of symplecticity (Skeel and Gear~\cite{skeel1992does}).

% ------------------------------------------------------------------------------
\subsubsection{Sundman--Stumpff Analytical Diffeomorphism under Retrograde Steps}
\label{sub:sundman_stumpff_diffeomorphism}

The multiscale splitting architecture on $\mathbf{Q}[N]$ resolves the negative time step deadlock by decoupling physical time from the internal integration of close encounters. Under the Fulton--MacPherson splitting on $\mathbf{Q}[N]$, every tight binary pair $C \in \mathcal{V}_{\mathrm{leaf}}(\mathcal{T})$ is integrated using the Sundman conformal clock $\mathrm{d}t_C = \rho_C \, \mathrm{d}\tau_C$. We establish that this transformation remains a global diffeomorphism under backward time steps:

\begin{propbox}[Backward Real-Analytic Invertibility of Universal Kepler Flows]{prop:sundman_stumpff_diffeomorphism}
In the Sundman--Stumpff universal two-body representation~\cite{sundman1912memoire, stumpff1959himmelsmechanik}, the physical time advance $\Delta t(\Delta \tau)$ generated by an isolated leaf pair $C \in \mathcal{V}_{\mathrm{leaf}}(\mathcal{T})$ with reduced mass $\mu_C$ and total mass $M_C$ is:
\begin{equation}
    \Delta t(\Delta \tau) = \rho_0 \Delta \tau \, c_1\bigl(\alpha \Delta \tau^2\bigr) + \sigma_0 \Delta \tau^2 \, c_2\bigl(\alpha \Delta \tau^2\bigr) + G M_C \Delta \tau^3 \, c_3\bigl(\alpha \Delta \tau^2\bigr),
    \label{eq:sundman_stumpff_time_relation}
\end{equation}
where $\rho_0 = \|\mathbf{r}_0\|$, $\sigma_0 = \mathbf{r}_0 \cdot \boldsymbol{\pi}_0 / \mu_C$, $\alpha = 2 G M_C / \rho_0 - \|\boldsymbol{\pi}_0\|^2 / \mu_C^2$, and $c_k(z) \coloneqq \sum_{n=0}^\infty \frac{(-1)^n z^n}{(2n+k)!}$ are the universal Stumpff transcendental series. For any Keplerian orbit with non-zero angular momentum ($L > 0$), the mapping $\Delta \tau \mapsto \Delta t(\Delta \tau)$ is a strictly monotonic, real-analytic ($\mathcal{C}^\omega$) diffeomorphism from $\mathbb{R}$ to $\mathbb{R}$. For exact radial collision orbits ($L = 0$), $\Delta t(\Delta \tau)$ is a strictly monotonic global homeomorphism from $\mathbb{R}$ to $\mathbb{R}$, smooth on $\mathbb{R} \setminus \{0\}$ with a classical $1/3$ Hölder singularity ($\Delta t \sim \Delta \tau^3$) at collision. For every physical time step $\Delta t \in \mathbb{R}$ (positive or negative), there exists a unique Sundman time increment $\Delta \tau \in \mathbb{R}$ satisfying:
\begin{equation}
    \operatorname{sgn}(\Delta \tau) = \operatorname{sgn}(\Delta t),
    \label{eq:sign_preservation_diffeo}
\end{equation}
and the universal state transition matrix preserves exact symplectic unimodularity $f \dot{g} - \dot{f} g = 1$ identically.
\end{propbox}

\begin{proofbox}[Proof of Proposition~\ref{prop:sundman_stumpff_diffeomorphism}]
Differentiating the Sundman relation~\eqref{eq:sundman_stumpff_time_relation} yields $\frac{\mathrm{d}(\Delta t)}{\mathrm{d}(\Delta \tau)} = \rho(\Delta \tau) \ge 0$, which is strictly positive everywhere for $L > 0$ and vanishes only at isolated collision instants for $L = 0$ (Stumpff~\cite{stumpff1959himmelsmechanik}; Hairer, Lubich and Wanner~\cite[p.~25]{hairer2006geometric}). Hence $\Delta t(\Delta \tau)$ is a strictly monotonic global homeomorphism on $\mathbb{R}$, smooth away from collision with Hölder exponent $1/3$ ($\Delta t \sim \Delta \tau^3$). In particular, $\operatorname{sgn}(\Delta \tau) = \operatorname{sgn}(\Delta t)$ holds for any step $\Delta t \in \mathbb{R}$, and the state transition coefficients satisfy the exact symplectic unimodularity identity $f \dot{g} - \dot{f} g = 1$ to floating-point precision floor.
\end{proofbox}

% ------------------------------------------------------------------------------
\subsubsection{Unitary McGehee Bounce Involution under Backward Sub-Steps}
\label{sub:mcgehee_bounce_unitary}

When an unsoftened binary pair experiences a near-head-on encounter or exact radial collision ($e \to 1$, $L \to 0$), the physical separation reaches the collision divisor $\rho_C = 0$. On the compactified Fulton--MacPherson space $\mathbf{Q}[N]$ (Section~\ref{ch:ch04_fulton_macpherson_splitting}), this locus corresponds to the exceptional boundary manifold $\mathbf{M}_C$. In blow-up coordinates $(\rho_C, \mathbf{s}_C, \mathbf{w}_C)$, where $\rho_C = \|\mathbf{r}_C\| \ge 0$, $\mathbf{s}_C = \mathbf{r}_C / \rho_C \in \mathbb{S}^2$ is the collision unit direction, and $\mathbf{w}_C = \rho_C^{1/2} \boldsymbol{\pi}_C / \sqrt{\mu_C} = w_r \mathbf{s}_C + \mathbf{w}_\theta$ partitions the scaled momentum into radial velocity $w_r = \mathbf{w}_C \cdot \mathbf{s}_C$ and tangential momentum $\mathbf{w}_\theta = \mathbf{w}_C - w_r \mathbf{s}_C$, the boundary manifold is $\mathbf{M}_C = \{\rho_C = 0\} \times \mathbb{S}^2 \times \mathbb{S}^2_{\sqrt{2 G M_C}}$, with energy condition $w_r^2 + \|\mathbf{w}_\theta\|^2 = 2 G M_C$ (and $\|\mathbf{w}_\theta\| = 0$ for radial collisions). While the regular continuation of binary collisions by elastic reflection traces back to Sundman~\cite{sundman1912memoire} and Levi-Civita~\cite{levicivita1904, levicivita1920regularisation}, McGehee's blow-up framework~\cite{mcgehee1974triple} formalizes this continuation on the collision manifold. We now establish that McGehee's collision boundary bounce operator regularizes collision crossings symmetrically under negative time sub-steps:

\begin{lemmabox}[Unitary Involutivity and Microscopic Reversibility of Collisions]{thm:mcgehee_bounce_unitary}
Following McGehee~\cite{mcgehee1974triple}, on the compactified collision boundary manifold $\mathbf{M}_C = \{\rho_C = 0\} \times \mathbb{S}^2 \times \mathbb{S}^2_{\sqrt{2 G M_C}}$, the specular radial bounce operator $\mathcal{B}_C : (\rho_C, \mathbf{s}_C, w_r, \mathbf{w}_\theta) \mapsto (\rho_C, \mathbf{s}_C, -w_r, \mathbf{w}_\theta)$ is an exact unitary involution:
\begin{equation}
    \mathcal{B}_C \circ \mathcal{B}_C = \operatorname{Id}_{\mathbf{M}_C}.
    \label{eq:bounce_involution_identity}
\end{equation}
Consequently, crossing $\mathbf{M}_C$ under any sequence of forward and backward sub-steps returns the orbit identically to its pre-collision state on the incoming sheet without phase discontinuity, dissipative loss, or numerical hysteresis.
\end{lemmabox}
\leanrefright{suzuki_barrier_theorems}{yoshida_order_four_suzuki_theorem}{01_YoshidaOrderFourSuzuki_Theorems.lean}
\pythonrefright{py_kernels}{solvers/{kepler_sundman, optimized_kepler_sundman, stumpff}.py, integrators/yoshida.py}{Float64 / Float128}

\noindent The formal certification badge denotes verification of the \emph{Constructive Combinatorial and 1-Jet Algebraic Kernel} in Lean~4, certifying the algebraic necessity of negative sub-steps ($w_0 < 0$), the cubic order-3 Lie cancellation $2w_1^3 + w_0^3 = 0$, and the unitary involutivity of the radial boundary bounce $\mathcal{B}_C^2 = \operatorname{Id}_{\mathbf{M}_C}$.

\begin{proofbox}[Proof of Lemma~\ref{thm:mcgehee_bounce_unitary}]
Following McGehee~\cite{mcgehee1974triple}, the specular radial bounce $\mathcal{B}_C : (\rho_C, \mathbf{s}_C, w_r, \mathbf{w}_\theta) \mapsto (\rho_C, \mathbf{s}_C, -w_r, \mathbf{w}_\theta)$ satisfies $\mathcal{B}_C^2(\rho_C, \mathbf{s}_C, w_r, \mathbf{w}_\theta) = (\rho_C, \mathbf{s}_C, -(-w_r), \mathbf{w}_\theta) = \operatorname{Id}_{\mathbf{M}_C}$. A retrograde sub-step $h_0 < 0$ following a forward crossing $h_1 > 0$ drives the system back along the time-reversed orbit by Proposition~\ref{prop:sundman_stumpff_diffeomorphism}; applying $\mathcal{B}_C$ upon return restores the pre-collision velocity $w_r^{\mathrm{in}}$ without dissipation or numerical hysteresis (see Figure~\ref{fig:ch05_mcgehee_bounce_involution}).
\end{proofbox}

\begin{figure}[htbp]
\centering
\resizebox{0.98\textwidth}{!}{%
\begin{tikzpicture}[>=stealth, font=\small]

    % ==========================================================================
    % PANEL (a): Blown-Up Configuration Space (\rho_C >= 0, \mathbf{s}_C \in S^2)
    % ==========================================================================
    \begin{scope}[shift={(-4.45, 0)}]
        \node[font=\footnotesize\bfseries] at (0.15, 2.80) {(a) Blown-up configuration $(\rho_C, \mathbf{s}_C)$};

        % Faint coordinate axes
        \draw[->, line width=0.45pt, black!35] (-1.95, 0) -- (2.55, 0) node[right=2pt, font=\scriptsize, black!75] {$\rho_C \mathbf{s}_C$};
        \draw[->, line width=0.45pt, black!35] (0, -1.85) -- (0, 2.05);

        % Blown-up collision sphere E_C = {\rho_C = 0} (radius 1.05)
        \fill[black!10] (0, 0) circle (1.05);
        \draw[line width=1.15pt, black] (0, 0) circle (1.05);
        \node[font=\scriptsize] at (-0.28, 0.24) {$E_C$};
        \node[font=\scriptsize] at (-0.28, -0.24) {$\rho_C = 0$};

        % Incoming and outgoing radial trajectories touching the boundary point (1.05, 0)
        \draw[->, line width=1.05pt, black] (2.40, -0.22) -- (1.13, -0.22);
        \node[below=2pt, font=\scriptsize] at (1.78, -0.22) {$w_r < 0$};

        \draw[->, line width=1.05pt, black] (1.13, 0.22) -- (2.40, 0.22);
        \node[above=2pt, font=\scriptsize] at (1.78, 0.22) {$w_r > 0$};

        % Specular boundary bounce point s_C on {\rho_C = 0}
        \filldraw[fill=white, draw=black, line width=0.95pt] (1.05, 0) circle (2.5pt);
        \node[above left=2pt, font=\scriptsize] at (0.98, 0.25) {$\mathbf{s}_C$};

        % Curved specular reflection arrow at the boundary circle
        \draw[<->, line width=0.8pt, black] (0.78, -0.25) arc (-135:135:0.35);
        \node[left=1pt, font=\scriptsize] at (0.62, 0) {$\mathcal{B}_C$};
    \end{scope}

    % Morphism arrow between (a) and (b)
    \draw[->, line width=0.9pt, black] (-0.95, 0) -- (0.15, 0)
        node[midway, above=3pt, font=\scriptsize] {$w_r = \rho_C^{1/2}\pi_r/\sqrt{\mu_C}$};

    % ==========================================================================
    % PANEL (b): McGehee Cotangent Phase Plane (\rho_C >= 0, w_r) & Suzuki Bounce
    % ==========================================================================
    \begin{scope}[shift={(3.55, 0)}]
        \node[font=\footnotesize\bfseries] at (0.45, 2.80) {(b) McGehee phase plane $(\rho_C, w_r)$};

        % Shaded forbidden region \rho_C < 0
        \fill[black!7] (-1.35, -2.05) rectangle (-0.80, 2.05);

        % Collision manifold wall M_C = {\rho_C = 0} at x = -0.80
        \draw[line width=1.2pt, black] (-0.80, -2.05) -- (-0.80, 2.05)
            node[above=2pt, font=\scriptsize] {$\mathbf{M}_C = \{\rho_C = 0\}$};

        % Horizontal axis \rho_C >= 0
        \draw[->, line width=0.6pt, black!65] (-1.35, 0) -- (2.75, 0)
            node[right=2pt, font=\scriptsize, black] {$\rho_C$};
        \node[below left=1pt, font=\scriptsize] at (-0.80, 0) {$0$};

        % Unscaled Cartesian momentum divergence \pi_r ~ -\rho_C^{-1/2} (dotted curve plunging down)
        \draw[dotted, line width=0.8pt, black!60, domain=-0.38:2.15, samples=40]
            plot (\x, {-0.55 / sqrt(\x + 0.65) - 0.55});
        \node[right=2pt, font=\scriptsize, black!70] at (-0.28, -1.92) {$\pi_r \sim -\rho_C^{-1/2} \to -\infty$};

        % Exact McGehee regular energy curve w_r^2 - 2\rho_C E_C = 2 G M_C
        % Lower incoming branch (w_r < 0) approaching (-0.80, -1.35)
        \draw[line width=1.15pt, black, domain=-0.80:2.35, samples=40]
            plot (\x, {-sqrt(1.8225 - 0.34*(\x + 0.80))});
        \draw[->, line width=1.0pt, black] (1.30, {-sqrt(1.8225 - 0.34*2.10)}) -- (0.45, {-sqrt(1.8225 - 0.34*1.25)});

        % Upper outgoing branch (w_r > 0) departing from (-0.80, +1.35)
        \draw[line width=1.15pt, black, domain=-0.80:2.35, samples=40]
            plot (\x, {sqrt(1.8225 - 0.34*(\x + 0.80))});
        \draw[->, line width=1.0pt, black] (0.45, {sqrt(1.8225 - 0.34*1.25)}) -- (1.30, {sqrt(1.8225 - 0.34*2.10)});

        % Retrograde negative sub-step arrows (w_0 h < 0, dashed) reversing the flow
        \draw[->, dashed, line width=0.75pt, black!75] (0.75, {sqrt(1.8225 - 0.34*1.55) - 0.18}) -- (-0.15, {sqrt(1.8225 - 0.34*0.65) - 0.18})
            node[midway, below=1pt, sloped, font=\scriptsize] {$w_0 h < 0$};
        \draw[->, dashed, line width=0.75pt, black!75] (-0.15, {-sqrt(1.8225 - 0.34*0.65) + 0.18}) -- (0.75, {-sqrt(1.8225 - 0.34*1.55) + 0.18})
            node[midway, above=1pt, sloped, font=\scriptsize, fill=white, inner sep=0.8pt] {$w_0 h < 0$};

        % Forward sub-step labels (w_1 h > 0)
        \node[below=2pt, font=\scriptsize] at (1.55, -1.02) {$w_1 h > 0$};
        \node[above=2pt, font=\scriptsize] at (1.55, 1.02) {$w_1 h > 0$};

        % Involutive McGehee bounce arrow along the collision wall \rho_C = 0
        \draw[<->, line width=1.05pt, black] (-0.52, -1.22) -- (-0.52, 1.22)
            node[pos=0.64, right=2pt, font=\scriptsize, fill=white, inner sep=1pt] {$\mathcal{B}_C\;(\mathcal{B}_C^2 = \operatorname{Id})$};

        % Boundary collision points w_r^- = -\sqrt{2GM_C} and w_r^+ = +\sqrt{2GM_C}
        \filldraw[fill=white, draw=black, line width=0.95pt] (-0.80, -1.35) circle (2.5pt)
            node[left=3pt, font=\scriptsize] {$-\sqrt{2GM_C}$};
        \filldraw[fill=white, draw=black, line width=0.95pt] (-0.80, 1.35) circle (2.5pt)
            node[left=3pt, font=\scriptsize] {$+\sqrt{2GM_C}$};
    \end{scope}

\end{tikzpicture}%
}
\caption{Geometric resolution of binary collisions and the Sheng--Suzuki negative-timestep barrier via McGehee's blow-up and unitary bounce involution (Lemma~\ref{thm:mcgehee_bounce_unitary}). \textbf{(a)} In blown-up configuration coordinates $(\rho_C, \mathbf{s}_C)$, the singular collision origin is replaced by the spherical boundary divisor $E_C = \{\rho_C = 0\}$, where the specular bounce $\mathcal{B}_C$ reflects the radial velocity $w_r \mapsto -w_r$ at $\mathbf{s}_C$. \textbf{(b)} In the McGehee phase plane $(\rho_C, w_r)$, whereas the unscaled Cartesian momentum $\pi_r \sim -\rho_C^{-1/2}$ diverges (dotted curve), the rescaled energy trajectory $w_r^2 - 2\rho_C E_C = 2GM_C$ reaches the collision manifold $\mathbf{M}_C = \{\rho_C = 0\}$ at finite endpoints $\mp\sqrt{2GM_C}$. Because $\mathcal{B}_C : -\sqrt{2GM_C} \leftrightarrow +\sqrt{2GM_C}$ satisfies $\mathcal{B}_C^2 = \operatorname{Id}_{\mathbf{M}_C}$, forward sub-steps ($w_1 h > 0$, solid arrows) and retrograde Suzuki/Yoshida sub-steps ($w_0 h < 0$, dashed arrows) traverse $\mathbf{M}_C$ reversibly without hysteresis.}
\label{fig:ch05_mcgehee_bounce_involution}
\end{figure}

% ------------------------------------------------------------------------------
\subsubsection{The 3-Stage Symmetric Yoshida-4 Propagator}
\label{sub:yoshida_propagator}

Armed with the global backward invertibility of Sundman--Stumpff flows (Proposition~\ref{prop:sundman_stumpff_diffeomorphism}) and the unitary collision bounce (Lemma~\ref{thm:mcgehee_bounce_unitary}), we construct the canonical fourth-order symmetric integrator:

\begin{propbox}[Order-4 Convergence and Invariance of the Yoshida Scheme]{prop:yoshida_order_four}
Following Yoshida~\cite{yoshida1990construction}, let $\mathcal{S}_2(h)$ be the second-order palindromic Strang scheme on $\mathcal{P}_{\mathrm{ext}}$ established in Proposition~\ref{prop:strang_bch_parity_bargmann}. The 3-stage symmetric composition:
\begin{equation}
    \mathcal{S}_4(h) = \mathcal{S}_2(w_1 h) \circ \mathcal{S}_2(w_0 h) \circ \mathcal{S}_2(w_1 h)
    \label{eq:yoshida_composition}
\end{equation}
with the unique real algebraic weights ($w_1 \approx 1.3512, w_0 \approx -1.7024$):
\begin{equation}
    w_1 = \frac{1}{2 - 2^{1/3}}, \qquad w_0 = -\frac{2^{1/3}}{2 - 2^{1/3}},
    \label{eq:yoshida_weights}
\end{equation}
satisfies the following structural properties:
\begin{enumerate}
    \item \textbf{Order-4 Convergence:} The second-order truncation error vanishes identically, and the formal shadow Hamiltonian expands exclusively in even powers:
    \begin{equation}
        \widetilde{H}_4(h) = H + h^4 H_5 + \mathcal{O}(h^6).
        \label{eq:yoshida_shadow_hamiltonian}
    \end{equation}
    \item \textbf{Unconditional Extended Symplecticity:} $\mathcal{S}_4(h)^* \Omega_{\mathrm{ext}} = \Omega_{\mathrm{ext}}$ on $\mathcal{P}_{\mathrm{ext}}$.
    \item \textbf{Exact Lie--Bargmann Invariance:} All $10$ generators $(\mathbf{P}, \mathbf{L}, \mathbf{G}, E)$ are conserved exactly across all three sub-steps: $\Delta \mathbf{P} = \mathbf{0}$, $\Delta \mathbf{L} = \mathbf{0}$, $\Delta \mathbf{G} = \mathbf{0}$, $\Delta E = 0$.
\end{enumerate}
\end{propbox}

\begin{proofbox}[Proof of Proposition~\ref{prop:yoshida_order_four}]
\textbf{Item (1):} Applying the Baker--Campbell--Hausdorff formula to the symmetric composition $\mathcal{S}_4(h) = \mathcal{S}_2(w_1 h) \circ \mathcal{S}_2(w_0 h) \circ \mathcal{S}_2(w_1 h)$, the order-3 error term is proportional to $(2w_1^3 + w_0^3) h^3 H_3$, which vanishes identically under the algebraic weights $w_1 = (2 - 2^{1/3})^{-1}$ and $w_0 = -2^{1/3} w_1$ (Yoshida~\cite[eq.~(2.11)]{yoshida1990construction}; Hairer, Lubich and Wanner~\cite[pp.~44--45]{hairer2006geometric}). Palindromic symmetry $\mathcal{S}_4(-h) = \mathcal{S}_4(h)^{-1}$ eliminates all remaining odd orders, leaving leading error $\mathcal{O}(h^4)$.

\textbf{Items (2) and (3):} By Proposition~\ref{prop:strang_bch_parity_bargmann}, each elementary sub-step $\mathcal{S}_2(w_i h)$ is an exact symplectomorphism preserving $\Omega_{\mathrm{ext}}$ and all ten Lie--Bargmann generators $(\mathbf{P}, \mathbf{L}, \mathbf{G}, E)$. Because the composition of symplectomorphisms is a symplectomorphism and symmetry invariants are conserved stage by stage, $\mathcal{S}_4(h)$ preserves $\Omega_{\mathrm{ext}}$ and the symmetry algebra unconditionally.
\end{proofbox}

\paragraph{Transition to Optimal 7-Stage Compositions and Boundary Correctors.}
Although Yoshida's 3-stage composition $\mathcal{S}_4(h)$ achieves the optimal formal order $p = 4$ with a minimal stage count ($s = 3$), it suffers from a significant practical drawback in celestial mechanics: its negative weight $w_0 \approx -1.702$ is large in absolute magnitude, exceeding the total macro-step ($|w_0| > 1$). In hierarchical systems featuring fast Keplerian inner binaries or eccentric orbits, this large backward excursion forces substantial coordinate deflections, increasing the $L_2$ norm of the residual order-5 and order-7 Lie commutators in the shadow Hamiltonian. 

To overcome this limitation without compromising symplecticity or Bargmann invariance, multi-stage symmetric compositions optimize the coefficient norm: following Blanes and Moan~\cite{blanes2002practical}, their optimal 7-stage scheme $\mathrm{BM}_4$ (studied in~\S\ref{sec:bm4_seven_stages}) reduces the negative step to a mild $a_4 \approx -0.4344$ while compressing the commutator error constant by over two orders of magnitude (a factor $> 400$). In Section~\ref{sec:symplectic_correctors}, this engine is further augmented on the regular hierarchical domain $\mathcal{D}_{\mathrm{reg}}$ by the partial boundary corrector $\mathrm{PBM}_4$.