% ============================================================================== % chapters/ch01_introduction/sec01_background/content.tex % Chapter 1: Introduction % ============================================================================== \section{Introduction} \label{sec:background_conormal} \label{ch:ch01_introduction} \subsection{The Gravitational $N$-Body Problem: Contemporary Relevance and Open Questions} The Newtonian $N$-body problem models the mutual gravitational interaction of $N \ge 3$ point masses $m_1, \dots, m_N > 0$ in three-dimensional Euclidean space, governed on the cotangent bundle $T^*(\R^{3N} \setminus \Delta)$ by the classical Hamiltonian: \begin{equation} H(\mathbf{q}, \mathbf{p}) = \sum_{i=1}^N \frac{\|\mathbf{p}_i\|^2}{2 m_i} - \sum_{1 \le i < j \le N} \frac{G m_i m_j}{\|\mathbf{q}_i - \mathbf{q}_j\|}, \label{eq:newton_hamiltonian} \end{equation} where $\mathbf{q}_i \in \R^3$ denote Cartesian positions, $\mathbf{p}_i = m_i \dot{\mathbf{q}}_i \in \R^3$ are conjugate momenta, and $\Delta = \bigcup_{1 \le i < j \le N} \{\mathbf{q} \in \R^{3N} \mid \mathbf{q}_i = \mathbf{q}_j\}$ represents the singular collision divisor. In contemporary astrophysics, the demand for accurate, long-term multi-body integration of this Hamiltonian has never been more acute. Investigations into the secular stability and chaotic diffusion of the Solar System over gigayear timescales, the orbital migration and resonant capture within tightly packed exoplanetary systems, the core-collapse dynamics of globular star clusters, and the relativistic coalescence of intermediate-mass black holes in galactic nuclei all rely upon numerical integrations spanning millions to billions of dynamical times. Within these dense multi-body environments, close encounters play a decisive dynamical role. In continuous phase space, seminal theorems of Saari \cite{saari1971improbability, saari1973improbability} establish that initial conditions leading to true collisions in the Newtonian $N$-body problem form a set of Lebesgue measure zero in the full phase space $\R^{6N}$ (and in the center-of-mass reduced phase space $\R^{6N-6}$). Among collisional configurations, simultaneous multiple collisions ($k \ge 3$) form a lower-dimensional stratum of measure zero, while non-collision singularities are physically pathological configurations conjectured to have measure zero (Saari \cite{saari1977singularities}; see also Xia \cite{xia1992existence}). Consequently, in continuous time, the dynamics is non-singular almost everywhere with respect to the physical Liouville measure. However, in discrete numerical integration with finite time step $h > 0$, generic trajectories (which constitute open sets of strictly positive Lebesgue measure) inevitably penetrate compact quasi-collision spheres of radius $\delta(h) \ll 1$, triggering catastrophic step-size collapse or artificial numerical ejections in conventional codes. While exact multiple collisions remain absent almost surely, resolving binary encounters and multiscale cluster hierarchies symplectically without potential softening ($\varepsilon = 0$) desingularizes these generic positive-measure quasi-collisional encounters and ensures robust, unsoftened propagation across all physical regimes (see Figure~\ref{fig:ch01_dyadic_trajectories_showcase}). \begin{figure}[htbp] \centering \resizebox{\textwidth}{!}{% \begin{tikzpicture}[>=stealth, font=\small] % ========================================================================= % TOP DYADIC STEP-LEVEL BEAD SCALE BAR (k = 0, 1, 2, 3) % ========================================================================= \begin{scope}[shift={(0, 5.45)}] \draw[line width=0.4pt, black!35] (-6.50, -0.30) rectangle (6.50, 0.30); % k = 0 (Macro-step h: White Bead) \filldraw[fill=white, draw=black, line width=0.55pt] (-5.35, 0.0) circle (1.75pt); \node[right=3pt, font=\scriptsize] at (-5.30, 0.0) {$k=0\; (\Delta t = h)$}; % k = 1 (Sub-step h/2: Light Gray Bead) \filldraw[fill=black!25, draw=black, line width=0.5pt] (-2.05, 0.0) circle (1.60pt); \node[right=3pt, font=\scriptsize] at (-2.00, 0.0) {$k=1\; (h/2)$}; % k = 2 (Sub-step h/4: Medium-Dark Gray Bead) \filldraw[fill=black!55, draw=black, line width=0.45pt] (1.05, 0.0) circle (1.45pt); \node[right=3pt, font=\scriptsize] at (1.10, 0.0) {$k=2\; (h/4)$}; % k = 3 (Sub-step h/8: Near-Black Bead) \filldraw[fill=black!85, draw=black, line width=0.4pt] (4.05, 0.0) circle (1.35pt); \node[right=3pt, font=\scriptsize] at (4.10, 0.0) {$k \ge 3\; (h/2^k)$}; \end{scope} % ========================================================================= % PANEL (a): CLOSED ECCENTRIC KEPLERIAN BINARY (e = 0.75) % ========================================================================= \begin{scope}[shift={(-3.55, 2.85)}] \node[font=\small\bfseries] at (0, 1.96) {(a) $\mathbf{Q}[2]\; (e = 0.75)$}; \draw[line width=0.4pt, black!35] (-2.95, -1.58) rectangle (2.95, 1.55); % Axes \draw[->, line width=0.3pt, black!30] (-2.55, 0) -- (2.55, 0) node[right=1pt, font=\scriptsize, black!65] {$q_x$}; \draw[->, line width=0.3pt, black!30] (0, -1.35) -- (0, 1.24) node[above=0.5pt, font=\scriptsize, black!65] {$q_y$}; % Concentric dyadic threshold rings around focus (1.170, 0.0) \fill[black!7] (1.170, 0.0) circle (0.62); \draw[line width=0.4pt, dashed, black!40] (1.170, 0.0) circle (0.62); \draw[line width=0.35pt, dotted, black!35] (1.170, 0.0) circle (1.28); \node[font=\tiny, black!65] at (1.82, -0.74) {$k=2$}; \node[font=\tiny, black!55] at (0.35, -1.32) {$k=1$}; \node[font=\tiny, black!55] at (-1.75, 1.18) {$k=0$}; % Exact Keplerian Ellipse Curve (e = 0.75) \draw[line width=0.72pt, black] plot[smooth cycle] coordinates { (-1.560,-0.000) (-1.552,-0.103) (-1.529,-0.205) (-1.490,-0.306) (-1.435,-0.405) (-1.363,-0.502) (-1.274,-0.595) (-1.168,-0.684) (-1.043,-0.767) (-0.898,-0.844) (-0.731,-0.911) (-0.541,-0.968) (-0.325,-1.009) (-0.079,-1.031) (0.201,-1.023) (0.520,-0.973) (0.887,-0.849) (1.087,-0.740) (1.292,-0.578) (1.475,-0.335) (1.560,0.000) (1.475,0.335) (1.292,0.578) (1.087,0.740) (0.887,0.849) (0.520,0.973) (0.201,1.023) (-0.079,1.031) (-0.325,1.009) (-0.541,0.968) (-0.731,0.911) (-0.898,0.844) (-1.043,0.767) (-1.168,0.684) (-1.274,0.595) (-1.363,0.502) (-1.435,0.405) (-1.490,0.306) (-1.529,0.205) (-1.552,0.103) }; % Direction arrows \draw[->, >=stealth, line width=0.65pt, black] (-0.325,-1.009) -- (-0.18,-1.026); \draw[->, >=stealth, line width=0.65pt, black] (-0.18,1.026) -- (-0.325,1.009); % Dyadic beads: k = 2 (black!55) \foreach \pt in {(0.395,-0.998), (1.241,-0.625), (1.241,0.625), (0.395,0.998)} { \filldraw[fill=black!55, draw=black, line width=0.45pt] \pt circle (1.45pt); } % Dyadic beads: k = 1 (black!25) \foreach \pt in {(-0.490,-0.980), (0.791,-0.889), (0.791,0.889), (-0.490,0.980)} { \filldraw[fill=black!25, draw=black, line width=0.5pt] \pt circle (1.60pt); } % Macro-step boundary beads: k = 0 (white) \foreach \pt in {(-1.560,0.000), (-1.401,-0.454), (-0.898,-0.844), (0.056,-1.031), (1.560,0.000), (0.056,1.031), (-0.898,0.844), (-1.401,0.454)} { \filldraw[fill=white, draw=black, line width=0.55pt] \pt circle (1.75pt); } % Gravitational Focus \fill[black] (1.170, 0.0) circle (2.2pt); \node[left=2pt, font=\scriptsize\bfseries] at (1.15, 0.16) {$\mathbf{q}_{\mathrm{attr}}$}; \node[above right=1pt and 2pt, font=\scriptsize] at (1.58, 0.04) {$\mathbf{q}_{\mathrm{peri}}$}; \node[left=2pt, font=\scriptsize] at (-1.62, 0.22) {$\mathbf{q}_{\mathrm{apo}}$}; \end{scope} % ========================================================================= % PANEL (b): HYPERBOLIC CLOSE-ENCOUNTER SLINGSHOT (e = 1.25) % ========================================================================= \begin{scope}[shift={(3.55, 2.85)}] \node[font=\small\bfseries] at (0, 1.96) {(b) $\mathbf{Q}[2]\; (e = 1.25,\; \Delta\theta = 106^\circ)$}; \draw[line width=0.4pt, black!35] (-2.95, -1.58) rectangle (2.95, 1.55); % Axes \draw[->, line width=0.3pt, black!30] (-2.55, -0.35) -- (2.55, -0.35) node[right=1pt, font=\scriptsize, black!65] {$q_x$}; \draw[->, line width=0.3pt, black!30] (0, -1.35) -- (0, 1.24) node[above=0.5pt, font=\scriptsize, black!65] {$q_y$}; % Asymptotes \draw[line width=0.35pt, dashed, black!35] (-2.05, 1.30) -- (0.0, -0.945) -- (2.05, 1.30); % Concentric dyadic threshold rings around focus (0.0, -0.35) \fill[black!7] (0.0, -0.35) circle (0.58); \draw[line width=0.4pt, dashed, black!40] (0.0, -0.35) circle (0.58); \draw[line width=0.35pt, dotted, black!35] (0.0, -0.35) circle (1.15); \node[font=\tiny, black!65] at (0.0, -1.14) {$k_{\mathrm{sym}}=3\; (h/8)$}; \node[font=\tiny, black!55] at (1.42, -0.82) {$k=1$}; \node[font=\tiny, black!55] at (-1.78, 0.32) {$k=0$}; % Exact Hyperbolic Slingshot Curve \draw[line width=0.72pt, black] plot[smooth] coordinates { (-1.928,1.128) (-1.828,1.004) (-1.725,0.878) (-1.619,0.751) (-1.511,0.621) (-1.399,0.489) (-1.282,0.354) (-1.159,0.217) (-1.029,0.075) (-0.888,-0.071) (-0.731,-0.222) (-0.549,-0.376) (-0.441,-0.453) (-0.317,-0.526) (-0.170,-0.584) (0.000,-0.609) (0.170,-0.584) (0.317,-0.526) (0.441,-0.453) (0.549,-0.376) (0.731,-0.222) (0.888,-0.071) (1.029,0.075) (1.159,0.217) (1.282,0.354) (1.399,0.489) (1.511,0.621) (1.619,0.751) (1.725,0.878) (1.828,1.004) (1.928,1.128) }; % Direction arrows \draw[->, >=stealth, line width=0.65pt, black] (-1.619,0.751) -- (-1.511,0.621); \draw[->, >=stealth, line width=0.65pt, black] (1.511,0.621) -- (1.619,0.751); % Dyadic beads: k = 3 (black!85) \foreach \pt in {(-0.506,-0.408), (-0.200,-0.575), (0.200,-0.575), (0.506,-0.408)} { \filldraw[fill=black!85, draw=black, line width=0.4pt] \pt circle (1.35pt); } % Dyadic beads: k = 2 (black!55) \foreach \pt in {(-0.368,-0.498), (0.368,-0.498)} { \filldraw[fill=black!55, draw=black, line width=0.45pt] \pt circle (1.45pt); } % Dyadic beads: k = 1 (black!25) \foreach \pt in {(-1.000,0.044), (0.000,-0.609), (1.000,0.044)} { \filldraw[fill=black!25, draw=black, line width=0.5pt] \pt circle (1.60pt); } % Macro-step boundary beads: k = 0 (white) \foreach \pt in {(-1.824,0.999), (-1.303,0.378), (-0.624,-0.315), (0.624,-0.315), (1.303,0.378), (1.824,0.999)} { \filldraw[fill=white, draw=black, line width=0.55pt] \pt circle (1.75pt); } % Gravitational Focus \fill[black] (0.0, -0.35) circle (2.2pt); \node[above=2pt, font=\scriptsize\bfseries] at (0.0, -0.32) {$\mathbf{q}_{\mathrm{attr}}$}; \node[left=3pt, font=\scriptsize] at (-1.88, 1.02) {$\mathbf{v}_\infty^-$}; \node[right=3pt, font=\scriptsize] at (1.88, 1.02) {$\mathbf{v}_\infty^+$}; \end{scope} % ========================================================================= % PANEL (c): ARENSTORF SYNODIC THREE-BODY MULTI-LOOP ORBIT (PCR3BP) % ========================================================================= \begin{scope}[shift={(-3.55, -1.15)}] \node[font=\small\bfseries] at (0, 1.96) {(c) $\mathrm{PCR3BP}\; (\mu = 0.01228,\; \mathbf{Q}[3])$}; \draw[line width=0.4pt, black!35] (-2.95, -1.58) rectangle (2.95, 1.55); % Axes \draw[->, line width=0.3pt, black!30] (-2.55, 0) -- (2.45, 0) node[right=1pt, font=\scriptsize, black!65] {$x_{\mathrm{syn}}$}; \draw[->, line width=0.3pt, black!30] (0, -1.35) -- (0, 1.24) node[above=0.5pt, font=\scriptsize, black!65] {$y_{\mathrm{syn}}$}; % Lunar Hill Sphere & Terrestrial Encounter Zone \fill[black!7] (0.106, 0.0) circle (0.36); \draw[line width=0.4pt, dashed, black!40] (0.106, 0.0) circle (0.36); \fill[black!7] (1.226, 0.0) circle (0.28); \draw[line width=0.4pt, dashed, black!40] (1.226, 0.0) circle (0.28); \node[font=\tiny, black!65] at (1.226, -0.44) {$R_{\mathrm{Hill}}$}; % Exact Arenstorf Orbit Curve \draw[line width=0.72pt, black] plot[smooth cycle] coordinates { (1.233,0.000) (1.088,-0.045) (0.970,-0.007) (0.850,0.077) (0.710,0.201) (0.530,0.350) (0.299,0.482) (0.049,0.561) (-0.171,0.596) (-0.345,0.621) (-0.475,0.656) (-0.566,0.707) (-0.624,0.774) (-0.650,0.851) (-0.648,0.935) (-0.618,1.019) (-0.566,1.098) (-0.494,1.167) (-0.408,1.222) (-0.311,1.259) (-0.210,1.278) (-0.111,1.275) (-0.020,1.253) (0.057,1.213) (0.115,1.155) (0.148,1.085) (0.152,1.003) (0.123,0.912) (0.056,0.812) (-0.047,0.697) (-0.182,0.557) (-0.327,0.381) (-0.450,0.176) (-0.538,-0.029) (-0.602,-0.203) (-0.661,-0.338) (-0.726,-0.433) (-0.800,-0.491) (-0.881,-0.515) (-0.964,-0.507) (-1.046,-0.471) (-1.121,-0.410) (-1.185,-0.327) (-1.233,-0.228) (-1.264,-0.117) (-1.274,0.000) (-1.264,0.117) (-1.233,0.228) (-1.185,0.327) (-1.121,0.410) (-1.046,0.471) (-0.964,0.507) (-0.881,0.515) (-0.800,0.491) (-0.726,0.433) (-0.661,0.338) (-0.602,0.203) (-0.538,0.029) (-0.450,-0.176) (-0.327,-0.381) (-0.182,-0.557) (-0.047,-0.697) (0.056,-0.812) (0.123,-0.912) (0.152,-1.003) (0.148,-1.085) (0.115,-1.155) (0.057,-1.213) (-0.020,-1.253) (-0.111,-1.275) (-0.210,-1.278) (-0.311,-1.259) (-0.408,-1.222) (-0.494,-1.167) (-0.566,-1.098) (-0.618,-1.019) (-0.648,-0.935) (-0.650,-0.851) (-0.624,-0.774) (-0.566,-0.707) (-0.475,-0.656) (-0.345,-0.621) (-0.171,-0.596) (0.049,-0.561) (0.299,-0.482) (0.530,-0.350) (0.710,-0.201) (0.850,-0.077) (0.970,0.007) (1.088,0.045) }; % Dyadic beads: k = 2 (black!55) \foreach \pt in {(0.985,-0.014), (0.390,0.438), (0.390,-0.438), (0.985,0.014)} { \filldraw[fill=black!55, draw=black, line width=0.45pt] \pt circle (1.45pt); } % Dyadic beads: k = 1 (black!25) \foreach \pt in {(0.748,0.167), (-0.587,-0.163), (-0.587,0.163), (0.748,-0.167)} { \filldraw[fill=black!25, draw=black, line width=0.5pt] \pt circle (1.60pt); } % Macro-step boundary beads: k = 0 (white) \foreach \pt in {(1.233,0.000), (-0.067,0.582), (-0.618,1.019), (0.021,1.235), (-0.182,0.557), (-0.840,-0.507), (-1.274,0.000), (-0.840,0.507), (-0.182,-0.557), (0.021,-1.235), (-0.618,-1.019), (-0.067,-0.582)} { \filldraw[fill=white, draw=black, line width=0.55pt] \pt circle (1.70pt); } % Earth and Moon Primaries \fill[black] (0.106, 0.0) circle (2.2pt); \node[below=2pt, font=\scriptsize\bfseries] at (0.106, -0.02) {$m_1$}; \fill[black] (1.226, 0.0) circle (1.6pt); \node[above=2pt, font=\scriptsize\bfseries] at (1.226, 0.04) {$m_2$}; \node[font=\tiny, black!60] at (-1.85, 1.22) {$\mathcal{T}_{13} \leftrightarrow \mathcal{T}_{23}$}; \end{scope} % ========================================================================= % PANEL (d): THREE-BODY FIGURE-8 CHOREOGRAPHY (CHENCINER-MONTGOMERY) % ========================================================================= \begin{scope}[shift={(3.55, -1.15)}] \node[font=\small\bfseries] at (0, 1.96) {(d) $\mathbf{Q}[3]\; (m_1 = m_2 = m_3 = 1)$}; \draw[line width=0.4pt, black!35] (-2.95, -1.58) rectangle (2.95, 1.55); % Axes \draw[->, line width=0.3pt, black!30] (-2.55, 0) -- (2.55, 0) node[right=1pt, font=\scriptsize, black!65] {$q_x$}; \draw[->, line width=0.3pt, black!30] (0, -1.35) -- (0, 1.24) node[above=0.5pt, font=\scriptsize, black!65] {$q_y$}; % Central mutual-encounter Voronoi switching zone \fill[black!7] (0.0, 0.0) circle (0.68); \draw[line width=0.4pt, dashed, black!40] (0.0, 0.0) circle (0.68); \node[font=\tiny, black!65] at (0.0, -0.92) {$k_{\mathrm{sym}} \ge 1\; (\mathcal{T}_{ij})$}; % Exact Figure-8 Lemniscate Curve \draw[line width=0.72pt, black] plot[smooth cycle] coordinates { (0.000,0.000) (-0.173,-0.158) (-0.352,-0.308) (-0.542,-0.438) (-0.742,-0.540) (-0.946,-0.604) (-1.143,-0.627) (-1.321,-0.611) (-1.474,-0.566) (-1.599,-0.501) (-1.698,-0.425) (-1.772,-0.343) (-1.827,-0.258) (-1.864,-0.172) (-1.885,-0.086) (-1.892,-0.000) (-1.885,0.086) (-1.864,0.172) (-1.827,0.258) (-1.772,0.343) (-1.698,0.425) (-1.599,0.501) (-1.474,0.566) (-1.321,0.611) (-1.143,0.627) (-0.946,0.604) (-0.742,0.540) (-0.542,0.438) (-0.352,0.308) (-0.173,0.158) (-0.000,0.000) (0.173,-0.158) (0.352,-0.308) (0.542,-0.438) (0.742,-0.540) (0.946,-0.604) (1.143,-0.627) (1.321,-0.611) (1.474,-0.566) (1.599,-0.501) (1.698,-0.425) (1.772,-0.343) (1.827,-0.258) (1.864,-0.172) (1.885,-0.086) (1.892,-0.000) (1.885,0.086) (1.864,0.172) (1.827,0.258) (1.772,0.343) (1.698,0.425) (1.599,0.501) (1.474,0.566) (1.321,0.611) (1.143,0.627) (0.946,0.604) (0.742,0.540) (0.542,0.438) (0.352,0.308) (0.173,0.158) }; % Direction arrows on lobes \draw[->, >=stealth, line width=0.65pt, black] (-1.321,0.611) -- (-1.143,0.627); \draw[->, >=stealth, line width=0.65pt, black] (1.143,-0.627) -- (1.321,-0.611); % Dyadic beads: k = 2 (black!55) \foreach \pt in {(-0.330,-0.290), (-1.071,-0.623), (-1.071,0.623), (-0.330,0.290), (0.330,-0.290), (1.071,-0.623), (1.071,0.623), (0.330,0.290)} { \filldraw[fill=black!55, draw=black, line width=0.45pt] \pt circle (1.45pt); } % Dyadic beads: k = 1 (black!25) \foreach \pt in {(-0.691,-0.518), (-0.691,0.518), (0.691,-0.518), (0.691,0.518)} { \filldraw[fill=black!25, draw=black, line width=0.5pt] \pt circle (1.60pt); } % Macro-step boundary beads: k = 0 (white) \foreach \pt in {(0.000,0.000), (-1.401,-0.591), (-1.892,0.000), (-1.401,0.591), (1.401,-0.591), (1.892,0.000), (1.401,0.591)} { \filldraw[fill=white, draw=black, line width=0.55pt] \pt circle (1.75pt); } % Three Choreography Masses at t = 0 \fill[black] (0.000, 0.000) circle (2.2pt); \node[above=3pt, font=\scriptsize\bfseries] at (0.0, 0.02) {$m_3$}; \fill[black] (-1.698, 0.425) circle (2.2pt); \node[above left=0.5pt and 2pt, font=\scriptsize\bfseries] at (-1.71, 0.43) {$m_1$}; \fill[black] (1.698, -0.425) circle (2.2pt); \node[below right=0.5pt and 2pt, font=\scriptsize\bfseries] at (1.71, -0.43) {$m_2$}; \end{scope} % ========================================================================= % PANEL (e): GRAZING ENCOUNTER WITH BINARY PAIR EXPULSION (Q[3]) % ========================================================================= \begin{scope}[shift={(-3.55, -5.15)}] \node[font=\small\bfseries] at (0, 1.96) {(e) $\mathbf{Q}[3]\; (\mathcal{T}_{13} \to \mathcal{T}_{12}\text{ Binary Expulsion})$}; \draw[line width=0.4pt, black!35] (-2.95, -1.58) rectangle (2.95, 1.55); % Axes \draw[->, line width=0.3pt, black!30] (-2.75, 0) -- (2.55, 0) node[right=1pt, font=\scriptsize, black!65] {$q_x$}; \draw[->, line width=0.3pt, black!30] (0, -1.35) -- (0, 1.24) node[above=0.5pt, font=\scriptsize, black!65] {$q_y$}; % Grazing resonant exchange / capture sphere \fill[black!7] (0.28, -0.05) circle (0.38); \draw[line width=0.4pt, dashed, black!40] (0.28, -0.05) circle (0.38); \node[font=\tiny, black!65] at (0.68, -0.58) {$k_{\mathrm{sym}}=3$}; % Curve 1 (Body m_1: enters from bottom, captured into expelled binary spiraling right) \draw[line width=0.72pt, black] plot[smooth] coordinates { (0.054,-1.139) (0.043,-1.059) (0.032,-0.949) (0.027,-0.799) (0.041,-0.595) (0.114,-0.316) (0.317,0.065) (0.365,0.195) (0.405,0.106) (0.561,-0.042) (0.821,-0.100) (1.071,-0.019) (1.223,0.146) (1.251,0.312) (1.184,0.393) (1.107,0.338) (1.104,0.265) (1.217,0.105) (1.446,0.043) (1.670,0.120) (1.798,0.282) (1.803,0.440) (1.718,0.512) (1.667,0.497) }; % Curve 2 (Body m_2: grazes from top-right, captured into expelled binary with m_1) \draw[line width=0.72pt, black] plot[smooth] coordinates { (0.960,1.053) (0.836,0.964) (0.702,0.861) (0.555,0.734) (0.395,0.568) (0.233,0.338) (0.158,0.181) (0.141,-0.053) (0.336,-0.217) (0.519,-0.180) (0.663,-0.107) (0.771,-0.009) (0.835,0.100) (0.854,0.200) (0.838,0.272) (0.805,0.308) (0.771,0.306) (0.749,0.273) (0.750,0.218) (0.779,0.150) (0.837,0.082) (0.923,0.024) (1.031,-0.013) (1.154,-0.020) (1.279,0.008) (1.389,0.070) (1.469,0.161) (1.509,0.262) (1.507,0.354) (1.473,0.419) (1.424,0.448) (1.376,0.441) (1.343,0.402) (1.334,0.342) (1.354,0.271) (1.405,0.201) (1.485,0.142) (1.587,0.106) (1.704,0.099) (1.820,0.128) }; % Curve 3 (Body m_3: grazed and violently recoils to bottom-left) \draw[line width=0.72pt, black] plot[smooth] coordinates { (-0.174,-0.094) (-0.039,-0.085) (0.107,-0.092) (0.259,-0.115) (0.403,-0.153) (0.493,-0.202) (0.478,-0.221) (0.382,-0.192) (0.150,-0.136) (-0.199,-0.240) (-0.575,-0.337) (-0.980,-0.434) (-1.417,-0.534) (-1.884,-0.639) (-2.312,-0.733) (-2.648,-0.806) }; % Outgoing direction arrows (recoil vs expelled binary) \draw[->, >=stealth, line width=0.65pt, black] (-2.243,-0.718) -- (-2.550,-0.785); \draw[->, >=stealth, line width=0.55pt, dashed, black!60] (0.45, 0.02) -- (2.22, 0.36); \node[font=\tiny, black!70] at (1.98, 0.72) {$(m_1, m_2)$}; \node[font=\tiny, black!70] at (-1.95, -0.36) {$m_3$ recoil}; % Dyadic beads: k = 3 (black!85) \foreach \pt in {(0.335,0.119), (0.182,-0.145), (0.324,-0.154)} { \filldraw[fill=black!85, draw=black, line width=0.4pt] \pt circle (1.35pt); } % Dyadic beads: k = 2 (black!55) \foreach \pt in {(0.046,-0.564), (0.375,0.543), (0.478,-0.193), (0.838,0.272), (1.184,0.393), (1.424,0.448), (1.798,0.282)} { \filldraw[fill=black!55, draw=black, line width=0.45pt] \pt circle (1.45pt); } % Dyadic beads: k = 1 (black!25) \foreach \pt in {(0.027,-0.841), (0.593,0.768), (-0.400,-0.294), (0.821,-0.100), (1.154,-0.020), (1.446,0.043)} { \filldraw[fill=black!25, draw=black, line width=0.5pt] \pt circle (1.60pt); } % Macro-step boundary beads: k = 0 (white) \foreach \pt in {(0.054,-1.139), (0.960,1.053), (-0.174,-0.094), (-0.980,-0.434), (-1.884,-0.639), (-2.648,-0.806), (1.667,0.497), (1.820,0.128)} { \filldraw[fill=white, draw=black, line width=0.55pt] \pt circle (1.75pt); } \node[right=2pt, font=\scriptsize] at (0.08, -1.16) {$m_1$}; \node[above right=1pt, font=\scriptsize] at (0.96, 1.04) {$m_2$}; \node[above left=1pt, font=\scriptsize] at (-0.18, -0.06) {$m_3$}; \end{scope} % ========================================================================= % PANEL (f): SIMULTANEOUS TRIPLE COLLISION & MCGEHEE SUPERNODE BOUNCE (k=3) % ========================================================================= \begin{scope}[shift={(3.55, -5.15)}] \node[font=\small\bfseries] at (0, 1.96) {(f) $\mathbf{Q}[3]\; (\rho_S \to 0,\; \mathcal{B}_S^2 = \operatorname{Id})$}; \draw[line width=0.4pt, black!35] (-2.95, -1.58) rectangle (2.95, 1.55); % Axes \draw[->, line width=0.3pt, black!30] (-2.55, 0) -- (2.55, 0) node[right=1pt, font=\scriptsize, black!65] {$q_x$}; \draw[->, line width=0.3pt, black!30] (0, -1.35) -- (0, 1.24) node[above=0.5pt, font=\scriptsize, black!65] {$q_y$}; % Outer dyadic ring and central McGehee Supernode Blow-Up Sphere \rho_S = 0 \draw[line width=0.35pt, dotted, black!35] (0,0) circle (0.72); \fill[black!10] (0,0) circle (0.26); \draw[line width=0.55pt, dashed, black!65] (0,0) circle (0.26); \draw[<-, >=stealth, line width=0.35pt, black!55] (0.23, -0.13) -- (0.68, -0.48); \node[right=1pt, font=\tiny, black!75] at (0.68, -0.50) {$\partial\mathbf{Q}[3]\; (\rho_S \to 0)$}; \node[font=\tiny, black!55] at (-1.22, 0.68) {$k_{\mathrm{sym}}=3$}; % Body 1 Trajectory (Symmetric C_3 Petal 1) \draw[line width=0.72pt, black] plot[smooth] coordinates { (-0.940,-1.300) (-0.843,-1.252) (-0.741,-1.195) (-0.633,-1.130) (-0.519,-1.054) (-0.398,-0.965) (-0.270,-0.857) (-0.134,-0.725) (0.008,-0.554) (0.079,-0.444) (0.145,-0.310) (0.188,-0.134) (0.145,0.084) (-0.022,0.230) (-0.196,0.280) (-0.346,0.290) (-0.476,0.283) (-0.591,0.267) (-0.790,0.221) (-0.959,0.167) (-1.106,0.108) (-1.236,0.048) (-1.352,-0.013) (-1.457,-0.074) (-1.596,-0.164) }; % Body 2 Trajectory (Symmetric C_3 Petal 2) \draw[line width=0.72pt, black] plot[smooth] coordinates { (1.596,-0.164) (1.457,-0.074) (1.352,-0.013) (1.236,0.048) (1.106,0.108) (0.959,0.167) (0.790,0.221) (0.591,0.267) (0.476,0.283) (0.346,0.290) (0.196,0.280) (0.022,0.230) (-0.145,0.084) (-0.188,-0.134) (-0.145,-0.310) (-0.079,-0.444) (-0.008,-0.554) (0.134,-0.725) (0.270,-0.857) (0.398,-0.965) (0.519,-1.054) (0.633,-1.130) (0.741,-1.195) (0.843,-1.252) (0.940,-1.300) }; % Body 3 Trajectory (Symmetric C_3 Petal 3) \draw[line width=0.72pt, black] plot[smooth] coordinates { (-0.656,1.420) (-0.662,1.356) (-0.665,1.239) (-0.663,1.113) (-0.654,0.976) (-0.637,0.827) (-0.607,0.663) (-0.561,0.479) (-0.527,0.378) (-0.483,0.270) (-0.424,0.154) (-0.341,0.030) (-0.210,-0.096) (0.000,-0.168) (0.210,-0.096) (0.341,0.030) (0.424,0.154) (0.483,0.270) (0.527,0.378) (0.561,0.479) (0.607,0.663) (0.637,0.827) (0.654,0.976) (0.663,1.113) (0.665,1.239) (0.662,1.356) (0.656,1.420) }; % Dyadic beads: k = 3 (black!85) \foreach \pt in {(0.189,-0.121), (0.145,0.084), (-0.010,0.225), (0.010,0.225), (-0.145,0.084), (-0.189,-0.121), (-0.200,-0.103), (0.000,-0.168), (0.200,-0.103)} { \filldraw[fill=black!85, draw=black, line width=0.4pt] \pt circle (1.35pt); } % Dyadic beads: k = 2 (black!55) \foreach \pt in {(-0.108,-0.696), (0.152,-0.290), (-0.175,0.277), (-0.657,0.254), (0.657,0.254), (0.175,0.277), (-0.152,-0.290), (0.108,-0.696), (-0.549,0.442), (-0.327,0.014), (0.327,0.014), (0.549,0.442)} { \filldraw[fill=black!55, draw=black, line width=0.45pt] \pt circle (1.45pt); } % Dyadic beads: k = 1 (black!25) \foreach \pt in {(-0.237,-0.827), (-0.835,0.208), (0.835,0.208), (0.237,-0.827), (-0.598,0.619), (0.598,0.619)} { \filldraw[fill=black!25, draw=black, line width=0.5pt] \pt circle (1.60pt); } % Macro-step boundary beads: k = 0 (white) \foreach \pt in {(-0.880,-1.271), (-0.474,-1.023), (0.026,-0.528), (-0.445,0.286), (-1.123,0.101), (-1.541,-0.127), (1.541,-0.127), (1.123,0.101), (0.445,0.286), (-0.026,-0.528), (0.474,-1.023), (0.880,-1.271), (-0.660,1.398), (-0.648,0.922), (-0.470,0.242), (0.470,0.242), (0.648,0.922), (0.660,1.398)} { \filldraw[fill=white, draw=black, line width=0.55pt] \pt circle (1.75pt); } \node[left=2pt, font=\scriptsize] at (-0.92, -1.28) {$m_1$}; \node[right=2pt, font=\scriptsize] at (1.56, -0.18) {$m_2$}; \node[left=2pt, font=\scriptsize] at (-0.68, 1.38) {$m_3$}; \end{scope} \end{tikzpicture}% } \caption{\textbf{Multi-regime trajectory showcase of unsoftened normal-cone symplectic integration with palindromic dyadic substepping.} Along every trajectory, white beads ($\circ$) mark uniform macro-step boundaries ($t_n = n h, k=0$), while progressively darker gray beads denote internal dyadic sub-steps ($\Delta t_n = h / 2^{k_{\mathrm{sym}, n}}$ for $k = 1, 2, \ge 3$) triggered automatically inside the shaded close-encounter spheres without potential softening ($\varepsilon = 0$). \textbf{(a)}~Closed eccentric Keplerian binary ($e = 0.75$) showing palindromic refinement $k_{\mathrm{sym}} : 0 \to 1 \to 2 \to 1 \to 0$ around pericenter $\mathbf{q}_{\mathrm{peri}}$. \textbf{(b)}~Unbound hyperbolic slingshot ($e = 1.25$, deflection $\Delta\theta = 106^\circ$) reaching $k_{\mathrm{sym}} = 3$ ($h/8$) at closest approach. \textbf{(c)}~Synodic multi-loop Arenstorf orbit in the Earth--Moon PCR3BP ($\mu = 0.012277471$), combining local Levi-Civita / Sundman leaf regularization with Hill-sphere chart switching ($\mathcal{T}_{13} \leftrightarrow \mathcal{T}_{23}$). \textbf{(d)}~Equal-mass three-body Figure-8 choreography (Chenciner--Montgomery~\cite{chenciner2000remarkable}), where democratic Voronoi chart transitions ($\mathcal{T}_{12} \to \mathcal{T}_{23} \to \mathcal{T}_{13}$) and symmetric dyadic sub-steps resolve the central crossings. \textbf{(e)}~Grazing three-body resonant encounter ($\mathcal{T}_{13} \to \mathcal{T}_{12}$) triggering the high-velocity recoil of $m_3$ (coasting at $k=0$) and the expulsion of a tightly bound binary pair $(m_1, m_2)$ tracked at $k \in \{1, 2\}$. \textbf{(f)}~Simultaneous three-body Lagrange collapse onto the blown-up McGehee supernode manifold ($\rho_S \to 0$), where the synchronous cluster hierarchy and unitary bounce involution ($\mathcal{B}_S^2 = \operatorname{Id}$) resolve the triple quasi-collision at $k_{\mathrm{sym}} = 3$.} \label{fig:ch01_dyadic_trajectories_showcase} \end{figure} \subsection{Survey of Contemporary Numerical Approaches and the Pathology of Invariant Drifting} To address this challenge, computational celestial mechanics has historically explored three distinct avenues, each of which entails substantive mathematical trade-offs: \begin{enumerate}[label=(\alph*), noitemsep, topsep=2pt] \item \textbf{High-Order Adaptive Non-Symplectic Integrators:} Modern astrophysical packages frequently deploy high-order non-symplectic schemes with adaptive step-size control, such as 15th-order Gauss--Radau quadrature with compensated summation~\cite{rein2015ias15, pham2024new}, embedded Runge--Kutta pairs, or Bulirsch--Stoer rational extrapolation techniques. In scenarios featuring chaotic scattering or highly eccentric orbits, adaptive step-size controllers automatically contract the integration time step to resolve steep acceleration gradients, achieving remarkable local accuracy. However, because these integrators do not preserve the canonical symplectic 2-form $\omega_0 = \sum_{i=1}^N \mathrm{d}\mathbf{q}_i \wedge \mathrm{d}\mathbf{p}_i$ nor the phase space volume, they do not possess an exact modified shadow Hamiltonian. Consequently, while standard non-symplectic integrators suffer from secular linear energy drift ($\mathcal{O}(t)$), even compensated high-order quadrature schemes that suppress linear energy drift down to Brouwer's random-walk roundoff limit ($\propto t^{1/2}$ in energy~\cite{brouwer1937measures}) inevitably suffer from a secular orbital phase error growing as $\mathcal{O}(t^{3/2})$ and lack an exact shadow Hamiltonian $\widetilde{H}$, eroding the long-term qualitative fidelity of secular resonances and semi-major axes. \item \textbf{Classical Analytical Regularization Techniques:} In the early twentieth century, Sundman \cite{sundman1912memoire} achieved the analytical regularization of isolated binary collisions by introducing a local conformal time reparameterization $\mathrm{d}t = r_{ij}\,\mathrm{d}u$ (Section~V, Eq.~71), and extended it globally to the 3-body problem via the product clock $\mathrm{d}t = \Gamma\,\mathrm{d}\omega$ with $\Gamma = \prod_{k=0}^2 (1 - e^{-r_k/l})$ (Section~XI, Eq.~179--180). This was extended to three dimensions via the spinor transformation of Kustaanheimo and Stiefel (KS) \cite{kustaanheimo1965perturbation} and subsequently generalized into chain regularization algorithms by Mikkola and Aarseth \cite{mikkola1993implementation, aarseth2003gravitational}. While exceptionally powerful for isolated binary pairs, classical regularizations are burdened by a fundamental drawback: the \emph{global time-freezing dilemma}. Because the regularizing clock $\mathrm{d}t = g(\mathbf{q})\,\mathrm{d}\tau$ is defined as a global scalar field over configuration space (such as $\Gamma \to 0$ in Sundman's global clock or $g = 1/(T+U) \to 0$ in chain methods), whenever any two bodies approach a collision ($r_{ij} \to 0$), the physical time increment $\mathrm{d}t$ collapses toward zero for the entire system. Consequently, the motion of distant, non-colliding bodies is artificially frozen, rendering global Sundman regularization computationally intractable for general $N \ge 3$ systems. \item \textbf{Phenomenological Potential Softening:} To avert singular gradient explosions in large-scale $N$-body simulations, numerical practitioners frequently replace the physical Newtonian potential with a softened kernel, such as the Plummer potential $V_\varepsilon(\mathbf{q}) = -G m_i m_j / \sqrt{\|\mathbf{q}_i - \mathbf{q}_j\|^2 + \varepsilon^2}$ with $\varepsilon > 0$, or impose an artificial force cutoff distance $d_0 > 0$. Softening is physically well-justified and indispensable in collisionless galactic dynamics and cosmological simulations ($N \sim 10^6 - 10^{11}$ particles), where it suppresses artificial two-body relaxation and converges toward the smooth mean-field Vlasov--Poisson kinetic limit. However, in collisional, point-mass systems such as compact planetary architectures, dense globular cluster cores, or multiple stellar systems, introducing an ad hoc length scale $\varepsilon > 0$ fundamentally alters the physical point-mass Hamiltonian: it suppresses large-angle gravitational deflections, introduces unphysical apsidal precession, shifts resonant commensurabilities, and precludes the faithful study of close gravitational encounters. \end{enumerate} Importantly, while smooth, stationary, radially symmetric softening potentials (such as $V_\varepsilon$) preserve the kinematic Galilean invariants ($\mathbf{P}, \mathbf{G}, \mathbf{L}$) via standard spatial Noether symmetries, they corrupt the physical energy spectrum and orbital geometry. The \emph{Lie--Bargmann invariant algebra} is then directly violated whenever practitioners introduce discontinuous force cutoffs ($d_0 > 0$), dynamical time-dependent softening scales $\varepsilon(t)$, or non-symplectic dissipations, permanently eroding the 10 universal first integrals of an isolated Newtonian system (total energy $H$, linear momentum $\mathbf{P}$, center-of-mass boost $\mathbf{G}$, and angular momentum $\mathbf{L}$) and degrading the long-term qualitative credibility of the simulated phase space trajectories. \subsection{The Geometric Promise and Singular Dilemmas of Symplectic Integration} Symplectic integrators \cite{hairer2006geometric} offer a rigorous remedy to secular dissipation. By construction, a numerical map $\Phi_h : (\mathbf{q}_n, \mathbf{p}_n) \mapsto (\mathbf{q}_{n+1}, \mathbf{p}_{n+1})$ is symplectic if its Jacobian matrix $J$ satisfies identically: \begin{equation} J^T \mathbb{J}_{6N} J = \mathbb{J}_{6N}, \qquad \mathbb{J}_{6N} = \begin{pmatrix} \mathbb{O}_{3N} & \mathbb{I}_{3N} \\ -\mathbb{I}_{3N} & \mathbb{O}_{3N} \end{pmatrix}. \end{equation} The decisive theoretical advantage of symplectic schemes is elucidated by \emph{Backward Error Analysis} \cite{hairer2006geometric, marsden2001discrete, benettin1994hamiltonian}: the discrete sequence generated by a symplectic integrator does not approximate the original continuous system with secularly drifting energy; rather, it tracks \emph{exactly} the trajectory of a nearby formal Hamiltonian, known as the \emph{shadow} or \emph{modified Hamiltonian}: \begin{equation} \widetilde{H}(\mathbf{q}, \mathbf{p}) = H(\mathbf{q}, \mathbf{p}) + h^2 H_2(\mathbf{q}, \mathbf{p}) + h^4 H_4(\mathbf{q}, \mathbf{p}) + \dots \end{equation} Because $\widetilde{H}$ is an exact invariant of the discrete flow, the energy error $H(\mathbf{q}_n, \mathbf{p}_n) - H(\mathbf{q}_0, \mathbf{p}_0)$ does not grow linearly with time, but remains confined within an adiabatic oscillatory envelope of amplitude $\mathcal{O}(h^p)$ over astronomical durations. Furthermore, whenever the underlying splitting respects the spatial symmetries of the potential, discrete Noether theorems guarantee that linear momentum $\mathbf{P}$ and angular momentum $\mathbf{L}$ are exact first integrals ($\Delta_{\mathrm{trunc}} \equiv \mathbf{0}$) without secular discretization drift. However, applying symplectic integrators to the unsoftened gravitational $N$-body problem encounters three formidable theoretical bottlenecks: \begin{enumerate}[label=(\roman*), noitemsep, topsep=2pt] \item \textbf{The Constant Time-Step Constraint:} As demonstrated by Skeel and Gear \cite{skeel1992does} and Hairer et al. \cite{hairer2006geometric}, standard adaptive step-size mechanisms based on local error estimators destroy the symplectic 2-form. Introducing a variable step $h_n = h(\mathbf{q}_n, \mathbf{p}_n)$ converts the discrete map into a non-Hamiltonian system, reinstating secular energy drift comparable to that of non-symplectic integrators. \item \textbf{Singularity of the Newtonian Potential:} Under a fixed time step $h$, when two bodies experience a close encounter ($r_{ij} \to 0$), the physical acceleration $\|\nabla V\| \sim r_{ij}^{-2}$ and its Hessian $\|\operatorname{Hess}(V)\| \sim r_{ij}^{-3}$ diverge rapidly. The modified Hamiltonian series $\widetilde{H}$ fails to converge, leading to numerical divergence or unphysical ejections. \item \textbf{The Sheng--Suzuki Barrier for High-Order Compositions:} To achieve high precision, one constructs higher-order integrators via symmetric composition of lower-order stages. However, a rigorous non-existence theorem proved independently by Sheng \cite{sheng1989solving} and Suzuki \cite{suzuki1991general} establishes that \emph{any symplectic partitioned Runge--Kutta or splitting method of order $p \ge 3$ with real coefficients must necessarily include at least one negative intermediate time step} ($w_k < 0$). In the presence of singular gravitational forces, propagating backward in time near a collision boundary triggers catastrophic numerical instabilities unless the singularity has been regularized in an analytically reversible manner. \end{enumerate} In the regular planetary regime where orbital eccentricities remain small and a single central mass dominates ($H = A + \varepsilon B$ with $\varepsilon \ll 1$), the celebratory perturbative splitting school pioneered by Wisdom and Holman~\cite{wisdom1991symplectic}, equipped with symplectic correctors~\cite{wisdom1996correctors, blanes2000processing} and generalized into high-order symmetric families by Laskar and Robutel~\cite{laskar2001large}, Blanes et al.~\cite{blanes2002practical, blanes2013optimized}, and Farrés, Laskar et al.~\cite{farres2013new_families}, achieves unmatched efficiency at fixed macro-steps. Yet these perturbative splittings rely fundamentally on the uniform analyticity of the interaction potential $B$ away from collision diagonals: whenever a subsystem enters an extreme eccentricity excursion ($e \to 1$), undergoes a grazing planetary flyby, or forms a democratic multi-body encounter ($r_{ij} \to 0$), the Lie commutator $\operatorname{ad}_A^n(B)$ explodes as $\mathcal{O}(r_{ij}^{-(n+1)})$ and the fixed-chart Kepler--drift hierarchy breaks down. Rather than competing with fixed-step planetary integrators on quasi-circular regular orbits, the theory developed in this paper constructs their \emph{universal, unsoftened singular extension} on the Fulton--MacPherson compactification $\mathbf{Q}[N]$: it seamlessly reduces to an exact hierarchy of commuting Keplerian leaves and cotangent-corrected $\mathrm{PBM}_4$ propagators in regular hierarchical regimes while traversing arbitrary binary collisions and dense multi-body encounters in a finite number of physical steps. \subsection{The Deformation to the Normal Cone and Canonical Hamiltonian Splittings} A promising mathematical bridge across this impasse has emerged from modern symplectic geometry. In a groundbreaking memoir, Damien Calaque and Pavel Safronov \cite{calaque2024deformation} formalized the \emph{deformation to the normal cone} within the framework of shifted cotangent structures and symplectic groupoids on derived stacks. While their foundational work operates in the abstract setting of derived algebraic geometry, the present treatise adopts this canonical geometric degeneration as a conceptual guiding paradigm, translating its core insight into a classical Hamiltonian operator splitting on the ordinary cotangent bundle $T^*M$. Let $M = \R^{3N}$ denote the ambient configuration space of the $N$ bodies, $X = M \setminus \Delta$ the regular non-collision locus (defined as the set complement of the collision divisor $\Delta$), and let $Z \subset \Delta \subset M$ be an isolated collision submanifold (such as a pairwise collision diagonal $\Delta_{ij} = \{\mathbf{q} \in M \mid \mathbf{q}_i = \mathbf{q}_j\}$). The deformation space $\mathcal{D}_Z M \to \mathbb{A}^1$ over the real affine line $\mathbb{A}^1 \cong \R$ (parameterized by a continuous deformation parameter $t \in [0, 1]$) interpolates geometrically between two distinct symplectic regimes: \begin{itemize}[noitemsep, topsep=2pt] \item \textbf{Generic Fiber ($t = 1$):} The fiber $(\mathcal{D}_Z M)_{t=1}$ is canonically isomorphic to the ambient configuration space $M$, and on the regular locus $X$ its cotangent bundle carries the standard physical Darboux form $\omega_0$. \item \textbf{Degenerated Normal Fiber ($t = 0$):} As $t \to 0$, the space degenerates onto the normal bundle $N_Z M$, and the cotangent space lifts to $T^*(N_Z M)$, which is endowed with a canonical, linear, flat symplectic structure. \end{itemize} Transposed to celestial mechanics, this normal cone perspective suggests a principled, non-heuristic method for Hamiltonian operator splitting. Rather than modifying the potential with ad hoc cutoffs, one decomposes the dynamics into: \begin{enumerate}[noitemsep, topsep=2pt] \item An exact groupoid exponential along the base manifold $X$, governed by kinetic inertia and generating unconstrained Galilean drift; \item A family of transverse impulses acting in the flat normal cotangent fibers $T^*(N_Z M)$ at $t = 0$. \end{enumerate} Crucially, in the flat normal cotangent fiber, the pairwise gravitational force between particles $i$ and $j$ acts strictly along the central direction $(\mathbf{q}_i - \mathbf{q}_j)$. By virtue of Newton's third law and the symmetry of the central interaction, the transverse impulse field is \emph{strictly orthogonal to the level sets of the 9 kinematic first integrals}: \begin{equation} \label{eq:1.4_orthogonality} \sum_{i=1}^N \nabla_{\mathbf{q}_i} V \equiv \mathbf{0}, \qquad \sum_{i=1}^N \mathbf{q}_i \times \nabla_{\mathbf{q}_i} V \equiv \mathbf{0}. \end{equation} Consequently, this normal cone splitting preserves total linear momentum $\mathbf{P}$, center of mass $\mathbf{G}$, and total angular momentum $\mathbf{L}$ unconditionally and identically at each stage of the integration, while energy $H$ is preserved through the shadow Hamiltonian $\widetilde{H}$ without secular dissipation. \subsection{Fulton--MacPherson Compactification: Resolving the Sundman Time-Freezing Dilemma} While binary collisions $Z = \Delta_{ij}$ form smooth submanifolds of codimension 3, general $N$-body systems ($N \ge 3$) introduce multiple collision strata $\Delta_I = \bigcap_{i, j \in I} \{\mathbf{q}_i = \mathbf{q}_j\}$ with non-transverse intersections. As Siegel \cite{siegel1941der} proved the non-analytic regularizability of triple collisions and McGehee \cite{mcgehee1974triple} uncovered chaotic dynamics on the blown-up collision manifold, the intrinsic normal cone along multiple collisions in $\R^{3N}$ is not a regular bundle, but a singular, reducible cone with intersecting branches, causing naive Cartesian splittings to generate unphysical "ghost poles" at cluster barycenters. To resolve this obstruction and decouple trivial center-of-mass motion, one passes to relative configurations via the algebraic quotient $M_{\mathrm{rel}} = M / \R^3 \cong (\R^3)^N / \R^3$. We then invoke the real Fulton--MacPherson compactification $\mathbf{Q}[N] = \operatorname{Bl}_{\Pi(N)}\left( M / \R^3 \right)$ \cite{fulton1994compactification, axelrod1994chern} (see Axelrod \& Singer~\cite[Section~5]{axelrod1994chern}), constructed by performing canonical spherical blow-ups along collision diagonals ordered by the partition lattice $\Pi(N)$ of $\{1, \dots, N\}$. On the compactified manifold with corners $\mathbf{Q}[N]$, the singular collision arrangement $\Delta$ becomes a \emph{simple normal crossing (SNC) divisor with corners}, whose boundary strata are smooth hypersurfaces with regular normal line bundles, and whose maximal corners of codimension $N-1$ are bijectively dual to binary Jacobi trees $\mathcal{T} \in \mathrm{Trees}(N)$. This compactification provides the definitive geometric resolution to the Sundman time-freezing dilemma and multi-body singular barriers: \begin{enumerate}[label=(\arabic*), noitemsep, topsep=1pt] \item \textbf{Decentralized Cluster Clocks:} Each hierarchical node $v$ of the Jacobi tree $\mathcal{T}$ carries its own \emph{local cluster clock} $\mathrm{d}t_v = \rho_v\,\mathrm{d}\tau_v$, where $\rho_v$ is the internal hyper-radius of cluster $C_v$. \item \textbf{Asynchronous Multiscale Evolution:} Fast internal orbital or collisional motions within tight sub-clusters evolve in their own local time, crossing pairwise collisions analytically via universal Stumpff functions $c_k(z)$ \cite{stumpff1959himmelsmechanik} and regularized elastic bounces $\mathcal{B}_C : v_r \mapsto -v_r$ in McGehee coordinates, while distant bodies drift symplectically without artificial deceleration. \item \textbf{Resolution of the Suzuki Barrier:} Because McGehee's regularized bounce is a unitary involution ($\mathcal{B}_C^2 = \mathrm{Id}$) combined with odd Stumpff transformations, the flow is strictly time-reversible across collisions, enabling the safe deployment of high-order symplectic compositions ($\mathcal{S}_4, \mathrm{BM}_4^{(7)}, \mathrm{BM}_6^{(11)}$) with negative intermediate time steps. \item \textbf{Elimination of the Multipolar Ghost Pole and Universal $\mathcal{C}^{1,1}$ Regularity:} Naive subtraction of a barycentric monopole creates an unphysical ``ghost pole'' $+G M_L M_R / \|\mathbf{r}_v\|$ whenever two cluster barycenters coincide ($\|\mathbf{r}_v\| \to 0$) without any physical collision. Introducing the scale-adapted screening function $\Theta(u_v)$ ($u_v = \|\mathbf{r}_v\|/(\bar{\rho}_L + \bar{\rho}_R)$) cancels both the monopole and dipole identically, leaving a pure quadrupolar tidal interaction governed by a symmetric trace-free tensor $\mathbf{T}_{\mathrm{apex}} \in \operatorname{Sym}_0(3)$ that is uniformly $\mathcal{C}^{1,1}$ Lipschitz across all of $\mathbf{Q}[N]$ without any softening cutoff ($d_0 = 0$, Theorem~\ref{thm:universal_c11_apex}). \item \textbf{Synchronous Anti-Zeno Super-Nodes for Multi-Body Encounters:} Whenever a dense sub-cluster of cardinality $k \ge 3$ crosses the scale-invariant Hill/Marchal--Sundman tidal threshold ($\eta_S \le \eta_{\mathrm{crit}} = 3^{-1/3}$, Proposition~\ref{prop:intrinsic_compactness_threshold}, stabilized by a spectral Schmitt trigger $[2^{-1/3}, 2^{-1/6}]$, Corollary~\ref{cor:hill_hysteresis_anticipation}), asynchronous pairwise splitting would trigger a catastrophic $\mathcal{O}(\delta_{\mathrm{cluster}}^{-5/2})$ Lie commutator divergence (Proposition~\ref{prop:excluded_middle_dilemma}). The sub-cluster dynamically fuses into a \emph{synchronous super-node} integrated under its local unfolded Sundman clock $\mathrm{d}t_S = \chi_S^2\,\mathrm{d}\tau_S$, guaranteeing traversal of arbitrary close multi-body encounters in a uniformly bounded, finite number of sub-steps (universal Anti-Zeno property, Theorem~\ref{thm:anti_zeno_regularization}). \end{enumerate} \subsection{Core Claims and Practical Numerical Schemes} To bridge the persistent chasm between algebraic geometry and computational astrophysics, this paper demonstrates that deep geometric structures directly govern high-performance symplectic integrators. Our overarching claim is anchored in a decisive measure-theoretic foundation: because collisions form a set of Lebesgue measure zero \cite{saari1971improbability, saari1973improbability} with multiple collisions occupying a lower-dimensional sub-stratum of measure zero, resolving pure binary close encounters via normal cone desingularization and decoupling cluster hierarchies via Fulton--MacPherson compactification yields an exact, unsoftened symplectic integration valid on \emph{every event set of non-zero measure} in phase space. We formulate our contributions along two complementary axes: foundational geometric structures and operational numerical schemes. \paragraph{1. Foundational Geometric and Symplectic Framework.} \enlargethispage{3\baselineskip} From an analytical and geometric perspective, this paper establishes four fundamental structural pillars: \begin{enumerate}[label=(\alph*), noitemsep, topsep=1pt] \item \textbf{Normal Cone Hamiltonian Splitting on $T^*M$:} Formalizing the graded Rees deformation to the normal cone $\pi : \mathcal{D}_Z M = \operatorname{Spec}_M \mathcal{R}(\mathcal{I}_Z) \to \mathbb{A}^1$ (Theorem~\ref{thm:calaque_safronov_degen}) decomposes dynamics into unconstrained groupoid/coadjoint base drifts $\Phi_h^{\mathcal{G}} = \exp_{\mathcal{G}}(h\mathbf{v}(\mathbf{p}_N))$ on the generic fiber $t = 1$ ($X = M \setminus \Delta$) and flat abelian cotangent fiber translations $\Phi_s^{\mathrm{trans}}(\mathbf{p}_N) = \mathbf{p}_N - s\nabla_n V(n)$ on the special fiber $t = 0$ ($T^*(N_{Z/M})$, endowed with the additively split symplectic $2$-form $\omega_{\mathrm{normal}} = \mathrm{d}x\wedge\mathrm{d}p + \mathrm{d}n\wedge\mathrm{d}w$, Proposition~\ref{prop:conormal_fiber_linearity}). This unifies Lie--Poisson Casimir-preserving splittings on $\mathfrak{se}(3)^*$ and $(\mathbb{S}^2)^K$ (Theorem~\ref{thm:casimir_preservation}, Propositions~\ref{prop:strang_time_reversibility} and~\ref{prop:heavy_top_casimirs}) with planar synodic and conformal regularizations (Proposition~\ref{prop:pcr3bp_jacobi}, Theorem~\ref{thm:conormal_bridge}), while guaranteeing pure Newtonian motion without softening ($\varepsilon = 0, d_0 = 0$) and unconditional exact conservation ($\Delta_{\mathrm{trunc}} \equiv \mathbf{0}$) of the 9 kinematic Lie--Bargmann first integrals ($\mathbf{P}, \mathbf{G}, \mathbf{L}$) without external projection. \item \textbf{Fulton--MacPherson Multiscale Desingularization on $T^*\mathbf{Q}[N]$:} Blowing up $M / \R^3$ by iterated real oriented spherical blow-ups along the irreducible building set $\mathcal{G}_{\mathrm{irr}}$ of the Whitney partition lattice $\Pi(N)$ transforms multiple collision strata into a simple normal crossing divisor with corners covered by $(2N-3)!!$ binary Jacobi trees, where disjoint sub-clusters commute identically on the standard cotangent bundle $T^*\mathbf{Q}[N]$ without the logarithmic boundary poles of Melrose's $b$-cotangent bundle ${}^b T^*\mathbf{Q}[N]$ (Theorem~\ref{thm:fulton_macpherson_snc}, Proposition~\ref{prop:disjoint_cotangent_commutation}). Associating decentralized clocks $\mathrm{d}t_v = \rho_v\,\mathrm{d}\tau_v$, $\mathcal{C}^{1,1}$ screened apex potentials (with dual $\mathcal{C}^\infty$ Gevrey-2 Cauchy mollifier and order-matched $\mathcal{C}^{2k}$ Hermite spline realizations, Theorem~\ref{thm:universal_c11_apex} and Propositions~\ref{prop:lie_bargmann_preservation}--\ref{prop:extended_palindromic_symplecticity}), and recursive synchronous Anti-Zeno super-nodes triggered at the parameter-free Hill/Marchal--Sundman threshold $\eta_{\mathrm{crit}} = 3^{-1/3}$ (Proposition~\ref{prop:excluded_middle_dilemma}, Theorem~\ref{thm:anti_zeno_regularization}, Proposition~\ref{prop:intrinsic_compactness_threshold}, Corollary~\ref{cor:hill_hysteresis_anticipation}) resolves both the Sundman time-freezing dilemma and the $\mathcal{O}(\delta_{\mathrm{cluster}}^{-5/2})$ commutator barrier. \item \textbf{Analytic Continuation and Time-Reversible Crossing of Binary Collisions:} Coupling McGehee's radial elastic bounce $\mathcal{B}_C : v_r \mapsto -v_r$ (Lemma~\ref{thm:mcgehee_bounce_unitary}) with the real-analytic Sundman--Stumpff diffeomorphism $\Delta\tau \mapsto \Delta t(\Delta\tau)$ (Proposition~\ref{prop:sundman_stumpff_diffeomorphism}) and Weierstrass pericenter conditioning (Proposition~\ref{prop:battin_unimodular_projection}) integrates pairwise collisions analytically across the singular divisor without cutoff. Because $\mathcal{B}_C$ is an exact unitary involution ($\mathcal{B}_C^2 = \mathrm{Id}$), the regularized trans-collisional flow is strictly time-reversible ($t \mapsto -t$), providing the indispensable mathematical stability required for negative intermediate time steps ($\beta\tau < 0$) in high-order symplectic compositions (Theorem~\ref{thm:sheng_suzuki_barrier}). \item \textbf{Universal Double Obstruction Theorem and Incompressible Quadratic Tidal Floor:} Within the bi-graded Poisson--Lie algebra $\mathcal{L}_{n,m}$ on $T^*\mathbf{Q}[N]$---where exact Jacobi--Kepler leaf integration collapses the unperturbed column $\{A_{\mathrm{JK}}, A_{\mathrm{JK}}\} \equiv 0 \implies \mathcal{L}_{n,0} \equiv \{0\}$ ($n \ge 2$) and spatial commutativity $\{W_{\mathrm{apex}}^{\mathrm{screen}}, V_{\mathrm{leaf}}\}_{\Omega_0} \equiv 0$ annihilates $\operatorname{ad}_{W_{\mathrm{apex}}^{\mathrm{screen}}}^3(A_{\mathrm{JK}}) \equiv 0$ (Proposition~\ref{prop:bigraded_bch_conormal}, Theorem~\ref{thm:bm6_order_conditions})---we prove that (i) momentum-parity alternation ($\operatorname{Im}(D_A|_{\mathcal{P}_0}) \subseteq \mathcal{P}_1$ versus $H_5^{\mathrm{linear}} \in \mathcal{P}_2 \oplus \mathcal{P}_0$) strictly forbids any pure configuration-space corrector $K(\mathbf{r}) \in \mathcal{P}_0$ from elevating $\mathrm{BM}_4^{(7)}$ to order 6, and (ii) a non-vanishing Hall cohomology class $[D_B D_A^3 B] \in H^1(D_A, \mathcal{L}_{5,2})$ ($\dim\mathcal{L}_{5,2} = 2 > \dim\operatorname{Im}(D_A|_{\mathcal{L}_{4,2}}) = 1$) with strictly positive Poincaré orbital average induces an incompressible quadratic tidal error floor $\mathcal{O}(\varepsilon^2 h^4)$ (Theorem~\ref{thm:universal_order6_obstruction}, Theorem~\ref{thm:c4_corrector_solvability}). \end{enumerate} \paragraph{2. Operational Numerical Schemes.} To serve practical computational needs in celestial mechanics and stellar dynamics, this theoretical architecture is realized as a concrete hierarchy of geometric integrators: \begin{enumerate}[label=(\alph*), noitemsep, topsep=1pt] \item \textbf{Strang Normal Cone Propagator $\mathcal{S}_2$:} A symmetric 7-operator second-order splitting on $(\mathcal{P}_{\mathrm{ext}}, \Omega_{\mathrm{ext}})$ interleaving barycentric drift, exact Kepler--Sundman leaf solvers, synchronous super-nodes, and $\mathcal{C}^{1,1}$ screened apex kicks, bounding energy errors within an $\mathcal{O}(h^2)$ envelope without secular dissipation while conserving all 9 kinematic Lie--Bargmann invariants exactly (Definition~\ref{def:strang_7_stage_propagator}, Theorem~\ref{thm:strang_extended_symplecticity}, Proposition~\ref{prop:strang_bch_parity_bargmann}). \item \textbf{Symmetric 4th-Order Yoshida Integrator $\mathcal{S}_4$:} A 3-stage symmetric composition $\mathcal{S}_4(h) = \mathcal{S}_2(\alpha h)\circ \mathcal{S}_2(\beta h)\circ \mathcal{S}_2(\alpha h)$ ($\alpha = 1/(2 - 2^{1/3}), \beta = -2^{1/3}/(2 - 2^{1/3})$)~\cite{yoshida1990construction}. Benefiting from the involutive time-reversibility of $\mathcal{B}_C$, it overcomes the Sheng--Suzuki barrier \cite{sheng1989solving, suzuki1991general}, executing trans-collisional steps with negative intermediate durations $\beta h < 0$ without instability (Proposition~\ref{prop:yoshida_order_four}). \item \textbf{Optimized 4th-Order 7-Stage Scheme $\mathrm{BM}_4^{(7)}$:} Higher-order symmetric composition of Blanes \& Moan~\cite{blanes2002practical} ($s=6$ non-trivial FSAL evaluations in direct ABA/RKN form, or $7$ FSAL evaluations in symmetric $\mathcal{S}_2^*$ composition, Table~\ref{tab:comparative_symplectic_synthesis}) with moderate retrograde sub-step ($a_4 \approx -0.4344$), compressing the principal BCH order-5 commutator error norm by a factor exceeding $400$ relative to Yoshida ($|C_5| \le 0.01236$ vs $5.3719$) and ensuring exponential Benettin--Giorgilli shadow confinement (Definition~\ref{def:bm4_seven_stages}, Theorem~\ref{thm:bm4_error_minimization}, Theorem~\ref{thm:unsoftened_shadow_confinement}, Proposition~\ref{prop:secular_phase_kpis}). \item \textbf{Optimized 6th-Order 11-Stage Scheme $\mathrm{BM}_6^{(11)}$:} An 11-stage symmetric palindromic composition ($s=10$ FSAL in direct RKN form $\mathrm{SRKN}_{11}^b$, or $11$ FSAL in $\mathcal{S}_2^*$ composition) solving the 7 independent order-6 RKN algebraic conditions on $\mathcal{L}_1, \mathcal{L}_3, \mathcal{L}_5$, bounding negative sub-steps ($|a_i| \le 0.36$), and variationally minimizing the residual order-7 commutator norm $\|H_7\|$ (Definition~\ref{def:bm6_eleven_stages}, Theorem~\ref{thm:bm6_order_conditions}, Theorem~\ref{thm:bm6_shadow_confinement}, Proposition~\ref{prop:bm6_secular_phase_kpis}). \item \textbf{Cotangent-Lifted Corrector $K_4 \in \mathcal{P}_1$ and Processed Scheme $\mathrm{PBM}_4$ (Pre-Asymptotic Order 6):} Bypassing the first parity obstruction via the cotangent lift $K_4(\mathbf{r}, \mathbf{w}) = c_4^{\mathrm{BM4}} h^4 D_A^3 B = c_4^{\mathrm{BM4}} h^4 \sum_\lambda \operatorname{Tr}(\mathbf{T}_{\mathrm{apex},\lambda}\,\mathbf{M}_{\mathrm{sym}}(\mathbf{r}_\lambda, \mathbf{w}_\lambda)) \in \mathcal{L}_{4,1} \cap \mathcal{P}_1$ ($c_4^{\mathrm{BM4}} \approx +1.7600 \times 10^{-4}$) realized by the 2-stage palindromic Lie conjugator $\mathcal{Z}_4^{\mathrm{BM4}}(h)$ (with antipodal half-shell factorization, Theorem~\ref{thm:c4_corrector_solvability}, Theorem~\ref{thm:shadow_alignment}, and Proposition~\ref{prop:corrector_lie_bargmann}), $\mathrm{PBM}_4$ annihilates the linear tidal error $\mathcal{O}(\varepsilon h^4)$ with zero internal torque ($\boldsymbol{\tau}_{\mathrm{int}} \equiv \mathbf{0}$) and $\mathcal{O}(r^{7/2})$ free-fall extinction across binary collisions (Theorem~\ref{thm:asymptotic_conormal_scaling}), achieving pre-asymptotic effective order 6 ($p_{\mathrm{eff}} = 6.25$, up to $2{,}537\times$ energy accuracy gain) down to the $\mathcal{O}(\varepsilon^2 h^4)$ floor at zero additional interior cost ($s=6$, Theorem~\ref{thm:c4_corrector_solvability}). \item \textbf{Covariant Symplectic Atlas Gauge Transition Operator $T_{(A,h_A)\to(B,h_B)}^{(\mathrm{PBM4})}$ and Canonical Voronoi Admissibility:} Dynamic restructuring across the $(2N-3)!!$ Jacobi tree charts via unimodular $\mathrm{SL}(2,\R) \hookrightarrow \operatorname{Sp}(6N-6,\R)$ tree rotations under the three-regime Hill-weighted Delaunay--Jacobi admissibility criterion at canonical Voronoi threshold $\gamma = 1.0$ and time-reversal-covariant adjoint hysteresis $\mathcal{A}_- = \mathcal{I}_T \circ \mathcal{A}_+^{-1} \circ \mathcal{I}_T$ (Theorem~\ref{thm:atlas_transition_anti_chattering}, Proposition~\ref{prop:delaunay_jacobi_admissibility}, Theorem~\ref{thm:anti_chattering_stability}, Remark~\ref{rem:adjoint_time_reversal_hysteresis}), glued in shadow space across both chart flips and dyadic step-size junctions by the non-commutative cotangent-lifted gauge operator $T_{(A,h_A)\to(B,h_B)}^{(\mathrm{PBM4})} = \mathcal{Z}_4^{\mathrm{BM4},(B)}(-h_B) \circ \mathbf{M}_{\mathrm{coord}} \circ \mathcal{Z}_4^{\mathrm{BM4},(A)}(+h_A)$ which closes the \v{C}ech 1-cocycle ($5.30 \times 10^{-36}$ algebraically, $4.08 \times 10^{-35}$ dynamically) and suppresses interfacial shadow energy jumps to machine precision ($2.36 \times 10^{-35}$ in \texttt{Float128} and $0.00$ to $4.19 \times 10^{-16}$ in \texttt{Float64}, Theorem~\ref{thm:pbm4_gauge_operator_exact} and Theorem~\ref{thm:pbm4_gauge_operator_sec6}). \item \textbf{Skeel--Gear Synchronized Two-Endpoint Dyadic Multirate Scheduler with Leaf/Apex Decoupling:} A time-reversible dyadic step hierarchy $h_k = \Delta t_{\mathrm{macro}}/2^{k_{\mathrm{sym},n}}$ governed by the two-endpoint symmetric selector $k_{\mathrm{sym},n} = \max(k_{\mathrm{dyn}}(\mathbf{z}_n), k_{\mathrm{dyn}}(\Phi_{\Delta t_n}^{(0)}(\mathbf{z}_n)), k_{\mathrm{sched}})$ with bit-exact scaled-integer timestamping (Remark~\ref{rem:ieee754_bit_exact_timestamps}), enforcing intra-block chart freezing (Proposition~\ref{thm:skeel_gear_boundary_sync}, Proposition~\ref{prop:causal_macro_boundary_ordering}), bilateral pericenter anticipation, and the Leaf/Apex Decoupling Theorem with two-regime quadrupolar tidal modulation ($\tau_{\mathrm{pair}}(C) = \min(\tau_{\mathrm{dyn}}(C), \max(\tau_{\mathrm{reg},1}, \tau_{\mathrm{reg},2}))$, Theorem~\ref{thm:leaf_apex_decoupling}, Proposition~\ref{prop:dyadic_scale_invariance}, Proposition~\ref{prop:dyadic_shadow_continuity}, and Theorem~\ref{thm:operational_leaf_apex_decoupling}), preventing dyadic step-size collapse ($k \le 2$, $k_{\mathrm{sym}} \equiv 0$ on isolated binaries) during extreme pericenter passages ($e \to 1$). \end{enumerate} \newpage \subsection{Outline of the Paper} The remainder of this treatise is organized into eight principal sections, progressing systematically from differential-geometric foundations and canonical Hamiltonian splittings to high-order processed symplectic schemes, asynchronous multirate algorithms, dual-precision empirical metrology, and overarching perspectives, complemented by a unified formal and computational appendix: \begingroup \renewcommand{\descriptionlabel}[1]{\hspace\labelsep\upshape\bfseries\rlap{#1:}}% \begin{description}[leftmargin=1.8em, itemsep=4pt, topsep=2pt] \item[\textbf{Section~\ref{ch:ch01_introduction} (Introduction and Common Mathematical Foundations)}] \leavevmode\par\nopagebreak Establishes the astrophysical motivation, the measure-theoretic foundation grounded in Saari's theorems guaranteeing unsoftened integration on all positive-measure phase-space events, the core theoretical and numerical claims, and the common mathematical primitives ($\S\ref{sec:common_foundations}$: relative configuration manifold $M_{\mathrm{rel}}$, singular collision arrangement $\Delta$, canonical Darboux symplectic structure $\omega_0$, and the 11-dimensional centrally extended Lie--Bargmann Poisson algebra $\mathfrak{bargmann}(3)$, Proposition~\ref{prop:lie_bargmann_closure}). \item[\textbf{Section~\ref{ch:ch02_conormal_foundations} (Calaque--Safronov Normal Cone Methodology and Algorithmic Recipes)}] \leavevmode\par\nopagebreak Develops the geometric theory of the deformation to the normal cone $\pi : \mathcal{D}_Z M = \operatorname{Spec}_M \mathcal{R}(\mathcal{I}_Z) \to \mathbb{A}^1$ via the graded Rees algebra and establishes the additive splitting of the degenerated symplectic $2$-form $\omega_{\mathrm{normal}}$ on $T^*(N_{Z/M})$ into base coordinates and flat abelian cotangent fiber translations ($\S\ref{sec:rees_deformation}$, Theorem~\ref{thm:calaque_safronov_degen}, Proposition~\ref{prop:conormal_fiber_linearity}). Formulates the canonical four-step normal-cone Strang splitting $\mathcal{S}_2(h) = \Phi_{h/2}^{NL} \circ \Phi_h^{\mathcal{G}} \circ \Phi_{h/2}^{NL}$ coupling $t=0$ transverse impulses with the $t=1$ symplectic groupoid exponential map $\Phi_h^{\mathcal{G}}$, proving unconditional Casimir conservation and exact time-symmetry with odd Baker--Campbell--Hausdorff generator parity ($\S\ref{sec:algorithmic_recipe}$, Theorem~\ref{thm:casimir_preservation}, Proposition~\ref{prop:strang_time_reversibility}). Demonstrates how this framework retrospectively unifies and explains the structural fidelity of explicit Lie--Poisson integrators for the $\mathrm{SE}(3)$ heavy rigid top on $\mathfrak{se}(3)^*$ and Bogomolov--Kirchhoff point-vortex dynamics on $(\mathbb{S}^2)^K$ ($\S\ref{sec:retrospective_cases}$, Proposition~\ref{prop:heavy_top_casimirs}). \item[\textbf{Section~\ref{ch:ch03_three_body_paradigm} (The Planar Three-Body Prototype and Classical Regularizations)}] \leavevmode\par\nopagebreak Investigates the planar circular restricted three-body problem (PCR3BP) in synodic coordinates as a foundational low-dimensional archetype, establishing the exact invariance of the Jacobi constant $C_J = -2 H_{\mathrm{syn}}$ (alongside Szebehely's gauge shift $C_{\mathrm{Sz}} = C_J + \mu(1-\mu)$) and the canonical decoupling between the exact linear Coriolis-kinematic flow $\Psi_h^{\mathrm{kin}} \in \operatorname{Sp}(4,\R)$ and transverse gravitational kicks ($\S\ref{sec:pcr3bp_canonical}$, Proposition~\ref{prop:pcr3bp_jacobi}). Proves that classical analytical regularizations---Sundman's conformal clock $\mathrm{d}t = r\,\mathrm{d}\tau$, Levi-Civita's branched $2:1$ spinor covering $z = w^2$ on $\mathbb{P}(N_Z M) \cong \mathbb{S}^1$, Birkhoff's global primary regularization, Kustaanheimo--Stiefel's 4D spinor lift, and McGehee's real oriented spherical blow-up $\operatorname{Bl}_Z M \cong \mathbb{S}^1 \times \R_{\ge 0}$---arise as local coordinate slices of the special fiber $t = 0$ of $\mathcal{D}_Z M$ while exposing their global time-freezing, spinor-branching, and non-transversality obstructions for $N \ge 3$ ($\S\ref{sec:classical_regularizations}$, Theorem~\ref{thm:conormal_bridge}, Remark~\ref{rem:sundman_conormal_bridge}). Validates the dyadic normal-cone solver governed by the harmonic Sundman scale $\chi(q) = r_1 r_2 / (r_1 + r_2)$ on multi-loop Arenstorf periodic orbits ($|\Delta C_J| < 10^{-14}$, $\S\ref{sec:arenstorf_validation}$, Remark~\ref{rem:conormal_subsumption}, Appendix~\ref{app:arenstorf_solver_sec}). \item[\textbf{Section~\ref{ch:ch04_fulton_macpherson_splitting} (Fulton--MacPherson Compactification and Multiscale Splitting)}] \leavevmode\par\nopagebreak Constructs the real oriented Fulton--MacPherson compactification $\mathbf{Q}[N] = \operatorname{Bl}_{\mathcal{G}_{\mathrm{irr}}}(M_{\mathrm{rel}})$ along the irreducible building set of the Whitney partition lattice $\Pi(N)$, turning $\Delta$ into a simple normal crossing divisor with corners covered by $(2N-3)!!$ binary Jacobi tree charts, contrasting the smooth standard cotangent bundle $T^*\mathbf{Q}[N]$ against Melrose's logarithmic $b$-cotangent bundle ${}^b T^*\mathbf{Q}[N]$, and proving exact cotangent commutation between disjoint sub-clusters ($\S\ref{sec:whitney_lattice}$, Theorem~\ref{thm:fulton_macpherson_snc}, Proposition~\ref{prop:disjoint_cotangent_commutation}). Neutralizes multipolar ``ghost poles'' via the scale-invariant screening function $\Theta(u_v)$ (in both $\mathcal{C}^\infty$ Gevrey-2 Cauchy mollifier and order-matched $\mathcal{C}^{2k}$ Hermite spline realizations) to establish universal $\mathcal{C}^{1,1}$ Lipschitz regularity, trace-free quadrupolar tidal structure $\mathbf{T}_{\mathrm{apex}} \in \operatorname{Sym}_0(3)$, exact preservation of the 10-generator Lie--Bargmann algebra on $(\mathcal{P}_{\mathrm{ext}}, \Omega_{\mathrm{ext}})$, and palindromic BCH odd-parity symplecticity ($\S\ref{sec:screened_apex}$, Theorem~\ref{thm:universal_c11_apex}, Propositions~\ref{prop:lie_bargmann_preservation} and~\ref{prop:extended_palindromic_symplecticity}). Identifies the $\mathcal{O}(\delta_{\mathrm{cluster}}^{-5/2})$ asynchronous commutator barrier and resolves it via synchronous super-nodes with universal Anti-Zeno regularization, parameter-free Hill/Marchal--Sundman compactness thresholds ($\eta_{\mathrm{crit}} = 3^{-1/3}$), spectral Schmitt triggers ($[2^{-1/3}, 2^{-1/6}]$), and unimodular $\mathrm{SL}(2,\R) \hookrightarrow \operatorname{Sp}(6N-6,\R)$ atlas transitions with positive dwell time under canonical Voronoi partitioning ($\gamma = 1.0$, with optional $\gamma = 0.6$ anti-chattering hysteresis; $\S\ref{sec:synchronous_supernodes}$, Proposition~\ref{prop:excluded_middle_dilemma}, Theorem~\ref{thm:anti_zeno_regularization}, Theorem~\ref{thm:atlas_transition_anti_chattering}, Proposition~\ref{prop:intrinsic_compactness_threshold}, Corollary~\ref{cor:hill_hysteresis_anticipation}). \item[\textbf{Section~\ref{ch:ch05_higher_order_symplectic} (Symplectic Discretization and High-Order Schemes)}] \leavevmode\par\nopagebreak Formulates the 7-operator palindromic Strang propagator $\mathcal{S}_2(h)$ on $(\mathcal{P}_{\mathrm{ext}}, \Omega_{\mathrm{ext}})$ with exact Noether conservation of the 9 kinematic Lie--Bargmann invariants ($\S\ref{sec:strang_propagator}$, Definition~\ref{def:strang_7_stage_propagator}, Theorem~\ref{thm:strang_extended_symplecticity}, Proposition~\ref{prop:strang_bch_parity_bargmann}), resolves the Sheng--Suzuki negative-timestep barrier via the real-analytic Sundman--Stumpff diffeomorphism and McGehee's unitary elastic bounce involution $\mathcal{B}_C^2 = \mathrm{Id}$ to certify the 4th-order Yoshida composition $\mathcal{S}_4(h)$ across exact binary collisions ($\S\ref{sec:suzuki_barrier_mcgehee}$, Theorem~\ref{thm:sheng_suzuki_barrier}, Proposition~\ref{prop:sundman_stumpff_diffeomorphism}, Lemma~\ref{thm:mcgehee_bounce_unitary}, Proposition~\ref{prop:yoshida_order_four}), and deploys the optimal 7-stage 4th-order scheme $\mathrm{BM}_4^{(7)}$ ($>400\times$ BCH error compression and exponential Benettin--Giorgilli shadow confinement at $\varepsilon = 0$, $\S\ref{sec:bm4_seven_stages}$, Definition~\ref{def:bm4_seven_stages}, Theorem~\ref{thm:bm4_error_minimization}, Theorem~\ref{thm:unsoftened_shadow_confinement}, Proposition~\ref{prop:secular_phase_kpis}) alongside the 11-stage 6th-order scheme $\mathrm{BM}_6^{(11)}$ and the bi-graded Lie algebra $\mathcal{L}_{n,m}$ governed by $\{A_{\mathrm{JK}}, A_{\mathrm{JK}}\}\equiv 0$ and $\operatorname{ad}_{W_{\mathrm{apex}}^{\mathrm{screen}}}^3(A_{\mathrm{JK}})\equiv 0$ ($\S\ref{sec:bm6_eleven_stages}$, Definition~\ref{def:bm6_eleven_stages}, Theorem~\ref{thm:bm6_order_conditions}, Theorem~\ref{thm:bm6_shadow_confinement}, Propositions~\ref{prop:bm6_secular_phase_kpis} and~\ref{prop:bigraded_bch_conormal}). Section~\ref{sec:symplectic_correctors} proves the \emph{Universal Double Obstruction Theorem} at order 6 (Theorem~\ref{thm:universal_order6_obstruction}), constructs the cotangent-lifted partial corrector $K_4 \in \mathcal{L}_{4,1} \cap \mathcal{P}_1$ and 2-stage palindromic conjugator $\mathcal{Z}_4^{\mathrm{BM4}}(h)$ characterized by the incompressible quadratic floor $\mathcal{O}(\varepsilon^2 h^4)$ and zero internal torque (Theorem~\ref{thm:c4_corrector_solvability}), establishes the 5-strata regularity partition of $T^*\mathbf{Q}[N]$ with $\mathcal{O}(r^{7/2})$ free-fall extinction of $K_4$ at binary collisions (Theorem~\ref{thm:asymptotic_conormal_scaling}), and derives the covariant cotangent-lifted gauge transition operator $T_{(A,h_A)\to(B,h_B)}^{(\mathrm{PBM4})}$ closing the \v{C}ech 1-cocycle (Theorem~\ref{thm:pbm4_gauge_operator_exact}). Finally, $\S\ref{sec:dyadic_multirate_scheduling}$ establishes the Skeel--Gear synchronized dyadic multirate lattice (Proposition~\ref{thm:skeel_gear_boundary_sync}), the Fundamental Leaf/Apex Decoupling Theorem with two-regime quadrupolar tidal modulation ($\tau_{\mathrm{pair}}(C) = \min(\tau_{\mathrm{dyn}}(C), \max(\tau_{\mathrm{reg},1}, \tau_{\mathrm{reg},2}))$, Theorem~\ref{thm:leaf_apex_decoupling}, Proposition~\ref{prop:dyadic_scale_invariance}), and dyadic shadow Hamiltonian continuity (Proposition~\ref{prop:dyadic_shadow_continuity}). \item[\textbf{Section~\ref{ch:ch06_algorithmic_architecture} (Computational Architecture and Algorithmic Realization)}] \leavevmode\par\nopagebreak Translates the mathematical theory into a deterministic, zero-allocation computational architecture: details the six-phase operational pipeline and 16-component state vector ($\S\ref{sec:global_pipeline}$), $\mathcal{O}(N^2\log N)$ agglomerative Jacobi tree construction under the unilateral Delaunay--Hill intruder guard ($\gamma_{\mathrm{build}} = 1.0$) and antipodal half-shell shadow pre-processing $\mathcal{Z}_4^{\mathrm{BM4}}(h)$ ($\S\ref{sec:initialization_shadow}$, Proposition~\ref{prop:tree_builder_algorithm}, Theorem~\ref{thm:shadow_alignment}, Proposition~\ref{prop:corrector_lie_bargmann}), the time-reversible two-endpoint predictive dyadic scheduler $k_{\mathrm{sym},n} = \max(k_{\mathrm{dyn}}(\mathbf{z}_n), k_{\mathrm{dyn}}(\Phi_{\Delta t_n}^{(0)}(\mathbf{z}_n)), k_{\mathrm{sched}})$ with bit-exact scaled-integer timestamping ($\S\ref{sec:dyadic_scheduler}$, Proposition~\ref{prop:macroblock_palindromic_symplecticity}, Theorem~\ref{thm:operational_leaf_apex_decoupling}, Remark~\ref{rem:ieee754_bit_exact_timestamps}), the analytical Kepler--Sundman leaf solver with Weierstrass pericenter conditioning and Battin unimodular projection alongside recursive normal-cone super-nodes and $\mathcal{C}^{2p}$ Hermite-screened apex kicks with exact monopole cancellation $V_{\mathrm{kick}} \equiv V_{\mathrm{inter}}$ ($\S\ref{sec:symplectic_subcycle}$, Propositions~\ref{prop:halley_bisection_fallback}, \ref{prop:battin_unimodular_projection}, and~\ref{prop:apex_galilean_invariance}), the strict four-step causal macro-boundary sequence, three-regime Hill-weighted Delaunay--Jacobi chart admissibility, time-reversal-covariant adjoint hysteresis $\mathcal{A}_-$, \v{C}ech 1-cocycle closure under $T_{(A,h_A)\to(B,h_B)}^{(\mathrm{PBM4})}$, and Battin Wronskian / Lie--Bargmann metrology ($\S\ref{sec:atlas_transitions_diagnostics}$, Propositions~\ref{prop:delaunay_jacobi_admissibility} and~\ref{prop:causal_macro_boundary_ordering}, Theorem~\ref{thm:anti_chattering_stability}, Remark~\ref{rem:adjoint_time_reversal_hysteresis}, Theorem~\ref{thm:pbm4_gauge_operator_sec6}, Proposition~\ref{prop:lie_bargmann_conservation_spectrum}), and high-performance vectorized and LLVM JIT-compiled micro-kernels ($0.0\text{ ULP}$ bit-exact equivalence and up to $62.54\times$ apex and $59.83\times$ super-node speedups across $N \in [3, 64]$) alongside 35-digit \texttt{Float128} reference solvers ($\S\ref{sec:computational_optimizations}$, Propositions~\ref{prop:apex_acceleration_invariance} and~\ref{prop:supernode_jit_performance}). \item[\textbf{Section~\ref{ch:ch07_benchmarks_and_lean4} (Numerical Benchmarks)}] \leavevmode\par\nopagebreak Conducts a comprehensive empirical validation following the spectral metrology of the Laskar school across dual-precision (53-bit \texttt{Float64} and 116-bit \texttt{Float128}) implementations. Section~\ref{sec:unsoftened_benchmarks} evaluates Laskar's inverse-cost metric $\tau/s$ on the Chenciner--Montgomery Figure-8 choreography~\cite{chenciner2000remarkable} and 3D eccentric hierarchical triples, isolates Brouwer's $\propto N^{1/2}$ roundoff regime from true asymptotic slopes ($p = 2.000, 3.999, 4.000, 4.000, 5.999$ down to $3.74 \times 10^{-31}$ in energy and $1.89 \times 10^{-31}$ in position, $\S\ref{sub:laskar_dual_convergence}$), experimentally confirms the two-channel $p_{\mathrm{eff}} = 6.25 \to 4.00$ crossover ($2{,}537\times$ energy and $1{,}335\times$ position gain) and $100\times$/decade $\mathcal{O}(\varepsilon^2 h^4)$ quadratic floor of Theorem~\ref{thm:c4_corrector_solvability} ($\S\ref{sub:pbm4_quadratic_floor_proof}$), and benchmarks four unsoftened singular limit cases ($\S\ref{sub:unsoftened_limit_cases}$: L1 analytic pericenter conditioning down to $1-e = 10^{-12}$ and $e = 1.0$ with $k_{\mathrm{sym}} \equiv 0$, L2 $\mathcal{O}(r^{7/2})$ asymptotic free-fall extinction of the $\mathrm{PBM}_4$ cotangent corrector across $r \in [10^{-12}, 10^0]\text{ AU}$, L3 exact McGehee central-collision traversals for $k \in \{2, 3, 4\}$, and L4 chaotic 3-4-5 Burrau scattering and $i_0 = 84.0^\circ$ 3D Kozai--Lidov secular resonance). Section~\ref{sec:long_term_stability} validates historical astronomical encounters ($\S\ref{sub:historical_encounters_validation}$) and long-term secular stability ($\S\ref{sub:secular_stability_laplace}$): Comet Lexell's 1770/1779 Jovian encounters ($3.031\,R_J$, demonstrating interfacial shadow-jump suppression down to $2.36 \times 10^{-35}$ in \texttt{Float128} and $0.00$ in \texttt{Float64} under $T_{(A,h_A)\to(B,h_B)}^{(\mathrm{PBM4})}$), Shoemaker--Levy~9's 1992 sub-Roche tidal disruption ($1.336\,R_J < 2.45\,R_J$, $\Lambda_S^{\max} = 16{,}690$), Asteroid Apophis's 2029 Aten $\to$ Apollo deflection and full-flyby palindromic time-reversibility, the 4:2:1 Jovian Laplace resonance under active dyadic substepping ($k_{\mathrm{sym}} \equiv 3$, $17.6\times$ max / $32.6\times$ RMS shadow energy reduction in \texttt{Float128} [$5.8\times$ / $7.5\times$ in \texttt{Float64}] and $179.85^\circ \pm 0.80^\circ$ resonant locking), and Inner Solar System secular integrations across $264$ dyadic junctions ($22.4\times$ \texttt{Float64} / $27.3\times$ \texttt{Float128} RMS shadow energy reduction, unbiased $\propto t^{1/2}$ Brouwer Noether diffusion, and $2{,}120\times$ Laskar NAFF frequency convergence down to the $3.70 \times 10^{-13}$ spectral floor with KAM symplectic conjugacy invariance~\cite{laskar1993frequency}). Finally, Section~\ref{sub:rebound_head_to_head} ($\S 7.3.1$--$\S 7.3.4$, including the cross-benchmark synthesis in $\S\ref{sub:rebound_cross_benchmark_synthesis}$) conducts three head-to-head comparative benchmark suites against the compiled C99 integrators of \texttt{REBOUND~5.2.2} (\texttt{IAS15}~\cite{rein2012rebound,rein2015ias15,pham2024new}, \texttt{TRACE}~\cite{lu2024trace}, \texttt{MERCURIUS}~\cite{rein2019mercurius}, and \texttt{WHFast}~\cite{rein2015whfast,rein2019highorder}) across star-grazing pericenters ($1-e_1 \in [10^{-9}, 10^{-1}]$ with $29{,}423\times$ tighter round-trip reversibility than \texttt{IAS15}) and exact collision $e_1 = 1.0$ ($9.13 \times 10^{11}\times$ more accurate than \texttt{TRACE}) with exact $0.00$ index-permutation invariance (\texttt{B1}), deep Jovicentric Hill-sphere crossings ($3.98 \times 10^{-4}\,R_{\mathrm{Hill},J}$, $17{,}082\times$ interfacial jump suppression via $T_{A\to B}^{(\mathrm{PBM4})}$, overcoming the Wisdom~\cite{wisdom2017symplectic} second-order hybrid switching barrier with $3.69 \times 10^7\times$ lower energy error and $61{,}734\times$ tighter reversibility than \texttt{TRACE}) and repeated giant-planet Hill crossings (\texttt{B2}), and democratic Pythagorean $3\text{--}4\text{--}5$ chaotic scattering across $27$ Voronoi chart transitions ($67.1\times$ dyadic and $1{,}834\times$ fixed-step tighter pre-horizon reversibility in \texttt{Float64}) and across the \texttt{Float64} Lyapunov horizon via 35-digit \texttt{Float128} ($1.80 \times 10^{20}\times$ tighter post-horizon reversibility at $T = 100$, \texttt{B4}). \item[\textbf{Section~\ref{ch:ch08_conclusion_outlook} (Conclusion and Outlook)}] \leavevmode\par\nopagebreak Synthesizes the theoretical, algebraic, and algorithmic achievements of the monograph alongside its demonstrated numerical capabilities and intrinsic mathematical boundaries ($\S\ref{sec:conclusion_synthesis}$), and formulates four open research questions regarding ultra-high-order processed cotangent schemes beyond $\mathrm{PBM}_4$, anisotropic multi-body collision regularizations ($k \ge 3$), $\mathcal{O}(N\log N)$ fast-multipole coupling for dense star clusters, and momentum-dependent post-Newtonian or stratified potentials ($\S\ref{sec:conclusion_perspectives}$). \item[\textbf{Appendix~\ref{sec:master_appendix} (Certified Lean~4 Formal Artifacts and Numerical Schemes)}] \leavevmode\par\nopagebreak Throughout Sections~\ref{ch:ch01_introduction} to~\ref{ch:ch07_benchmarks_and_lean4}, definitions, propositions, theorems, algorithms, and empirical benchmarks are anchored by verification badges cross-referencing Appendix~\ref{sec:master_appendix} ($\S\text{A.1}$--$\S\ref{app:github_benchmark_suite}$), which co-locates the complete, zero-\texttt{sorry} mechanized Lean~4 formal proofs (35 modules comprising 125 formal declarations covering the Constructive 1-Jet Algebraic Kernel, graded Rees degeneration, planar synodic and Levi-Civita regularizations, Whitney partition lattice, $\mathcal{C}^{1,1}$ apex bounds, Anti-Zeno super-nodes, Universal Double Obstruction Theorem, cotangent corrector $K_4 \in \mathcal{P}_1$, and \v{C}ech 1-cocycle gauge transitions) and the complete directory tree and module-by-module locator of the dual-precision companion GitHub repository (\texttt{normal\_cone\_nbody\_float64/} and \texttt{normal\_cone\_nbody\_float128/}, $\S\ref{app:github_benchmark_suite}$). \end{description} \endgroup \newpage \subsection{Common Mathematical Foundations and Notational Primitives} \label{sec:common_foundations} To bridge continuous differential geometry with mechanized verification, our accompanying Lean~4 development~\cite{moura2015lean} formalizes the \emph{Constructive 1-Jet Algebraic Kernel}: by operating on algebraic 1-jets of observables over an arbitrary commutative ring rather than infinite-dimensional Fréchet manifolds, it certifies the finite-dimensional tangent-space symplectic structure, the constructive duality of collision loci, and the closed Lie algebra of observables that underpin our numerical integrators. \subsubsection{Ambient and Relative Configuration Manifolds} Let $N \ge 3$ be the number of gravitationally interacting point masses $m_1, \dots, m_N > 0$. The \emph{ambient configuration manifold} models the unconstrained spatial positions of all $N$ bodies in physical three-dimensional space: \begin{equation} M = (\R^3)^N = \bigl\{ \mathbf{q} = (\mathbf{q}_1, \dots, \mathbf{q}_N) \bigm| \mathbf{q}_i \in \R^3 \bigr\} \cong \R^{3N}. \end{equation} The manifold $M$ is naturally equipped with the mass-weighted Riemannian kinetic metric: \begin{equation} \langle \mathbf{v}, \mathbf{w} \rangle_{\mathbb{M}} = \sum_{i=1}^N m_i \mathbf{v}_i \cdot \mathbf{w}_i = \mathbf{v}^T \mathbb{M} \mathbf{w}, \end{equation} where $\mathbb{M} = \operatorname{diag}(m_1 \mathbb{I}_3, \dots, m_N \mathbb{I}_3)$ denotes the positive-definite mass tensor, and the total mass is $\Mtot = \sum_{i=1}^N m_i > 0$. \leanrefright{ambient_config}{AmbientConfiguration}{01_ConfigurationSpace_Defs.lean} The diagonal linear subspace of uniform spatial translations is defined by: \begin{equation} \Delta_{\mathrm{trans}} = \bigl\{ \mathbf{q} \in M \bigm| \mathbf{q}_1 = \dots = \mathbf{q}_N \bigr\} \subset M \cong \R^3. \end{equation} Under an arbitrary spatial translation $\mathbf{q}_i \mapsto \mathbf{q}_i + \mathbf{a}$ with $\mathbf{a} \in \R^3$, the center of mass vector $\mathbf{R}_{\mathrm{cdm}}(\mathbf{q}) = \frac{1}{\Mtot} \sum_{i=1}^N m_i \mathbf{q}_i$ undergoes a pure rigid translation $\mathbf{R}_{\mathrm{cdm}} \mapsto \mathbf{R}_{\mathrm{cdm}} + \mathbf{a}$. To decouple this trivial Galilean drift, we identify the \emph{relative configuration manifold} with the orthogonal barycentric slice $M_0 = \Delta_{\mathrm{trans}}^\perp = \{ \mathbf{q} \in M \mid \sum_{i=1}^N m_i \mathbf{q}_i = \mathbf{0} \}$: \begin{equation} M_{\mathrm{rel}} = M / \Delta_{\mathrm{trans}} \cong M_0 \cong (\R^3)^{N-1}. \end{equation} \leanrefright{relative_config}{RelativeConfiguration}{01_ConfigurationSpace_Defs.lean} \subsubsection{The Singular Collision Arrangement and Locus} For each distinct pair of particles $1 \le i < j \le N$, the pairwise binary collision locus is the linear submanifold: \begin{equation} \Delta_{ij} = \bigl\{ \mathbf{q} \in M \bigm| \mathbf{q}_i = \mathbf{q}_j \bigr\} \subset M, \end{equation} which is of codimension~3 in $M$. It is algebraically cut out by the pairwise vanishing ideal: \begin{equation} \mathcal{I}(\Delta_{ij}) = \langle q_{i,a} - q_{j,a} \mid a = 1, 2, 3 \rangle \subset \mathcal{O}_M. \end{equation} \leanrefright{pairwise_collision}{BinaryCollision}{02_CollisionDivisor_Defs.lean} The global singular collision arrangement $\Delta$ is the union of all pairwise collision diagonals across all possible particle pairs: \begin{equation} \Delta = \bigcup_{1 \le i < j \le N} \Delta_{ij} \subset M, \qquad \mathcal{I}_\Delta = \bigcap_{1 \le i < j \le N} \mathcal{I}(\Delta_{ij}) \subset \mathcal{O}_M. \end{equation} \leanrefright{collision_divisor}{InCollisionLocus}{02_CollisionDivisor_Defs.lean} The physical domain of collisionless Newtonian motion is the \emph{regular configuration locus}, defined by excising the collision arrangement from the configuration manifold: \begin{equation} X = M \setminus \Delta, \qquad X_{\mathrm{rel}} = M_0 \setminus (M_0 \cap \Delta). \end{equation} \leanrefright{regular_locus}{NonCoincident}{02_CollisionDivisor_Defs.lean} The regular locus $X$ represents configurations where all particles remain pairwise distinct. The certified formal duality establishes that this condition is constructively and logically equivalent to the negation of membership in the collision locus: \leanrefright{regular_duality}{regular_iff_not_collision}{03_CollisionDivisor_Theorems.lean} \subsubsection{Symplectic Phase Space and Lie--Bargmann Algebra} The physical phase space of the system is the cotangent bundle $T^*X \subset T^*M \cong \R^{3N} \times \R^{3N}$, parameterized by the canonical position and momentum coordinates $(\mathbf{q}, \mathbf{p}) = (\mathbf{q}_1, \dots, \mathbf{q}_N, \mathbf{p}_1, \dots, \mathbf{p}_N)$ with $\mathbf{p}_i = m_i \dot{\mathbf{q}}_i \in \R^3$. The cotangent bundle carries the canonical Liouville 1-form $\theta_0 = \sum_{i=1}^N \mathbf{p}_i \cdot \mathrm{d}\mathbf{q}_i$ and the exact Darboux symplectic 2-form: \begin{equation} \omega_0 = -\mathrm{d}\theta_0 = \sum_{i=1}^N \sum_{a=1}^3 \mathrm{d}\mathbf{q}_{i,a} \wedge \mathrm{d}\mathbf{p}_{i,a}. \end{equation} \leanrefright{symplectic_form}{omega0}{04_SymplecticForm_Defs.lean} The symplectic form $\omega_0$ is closed ($\mathrm{d}\omega_0 = 0$), skew-symmetric, and alternating on all tangent variations: \leanrefright{symplectic_skew}{omega0_skew}{05_SymplecticForm_Theorems.lean} It satisfies constructive non-degeneracy on the configuration tangent space: \begin{equation} \forall v \in T(T^*X), \quad \bigl( \forall w \in T(T^*X), \, \omega_0(v, w) = 0 \bigr) \implies v = 0. \end{equation} \leanrefright{symplectic_nondeg}{omega0_tangent_nondeg}{05_SymplecticForm_Theorems.lean} On smooth observables $f, g \in \mathcal{C}^\infty(T^*X)$, the canonical Poisson bracket induced by $\omega_0$ is given by: \begin{equation} \{f, g\}_0 = \sum_{i=1}^N \left( \frac{\partial f}{\partial \mathbf{q}_i} \cdot \frac{\partial g}{\partial \mathbf{p}_i} - \frac{\partial f}{\partial \mathbf{p}_i} \cdot \frac{\partial g}{\partial \mathbf{q}_i} \right). \end{equation} \leanrefright{bargmann_obs}{poissonBracketVal}{06_BargmannAlgebra_Defs.lean} The unsoftened Newtonian dynamics is governed by the kinetic energy $T(\mathbf{p}) = \frac{1}{2}\mathbf{p}^T \mathbb{M}^{-1}\mathbf{p}$ and the singular pairwise gravitational potential $V(\mathbf{q}) = -\sum_{1 \le i < j \le N} \frac{G m_i m_j}{\|\mathbf{q}_i - \mathbf{q}_j\|}$, yielding the fundamental autonomous first integrals: \begin{equation} \mathbf{P}_{\mathrm{tot}} = \sum_{i=1}^N \mathbf{p}_i, \quad \mathbf{L}_{\mathrm{tot}} = \sum_{i=1}^N \mathbf{q}_i \times \mathbf{p}_i, \quad \mathbf{K} = \sum_{i=1}^N m_i \mathbf{q}_i = \Mtot \mathbf{R}_{\mathrm{cdm}}, \quad H = T(\mathbf{p}) + V(\mathbf{q}). \end{equation} \begin{propbox}[Central Poisson Closure of the Lie--Bargmann Algebra in $T^*X$]{prop:lie_bargmann_closure} Following Bargmann~\cite{bargmann1954unitary}, the autonomous observables $(\mathbf{P}_{\mathrm{tot}}, \mathbf{L}_{\mathrm{tot}}, \mathbf{K}, H, \mathbf{1})$ generate an 11-dimensional closed Poisson Lie algebra on $T^*X$ isomorphic to the centrally extended Galilean algebra $\mathfrak{bargmann}(3) \cong \mathfrak{gal}(3) \oplus_\omega \R \mathbf{1}$, satisfying the canonical brackets: \begin{equation} \begin{aligned} \{K_a, P_b\}_0 &= \Mtot \delta_{ab} \mathbf{1}, & \qquad \{L_a, K_b\}_0 &= \varepsilon_{abc} K_c, \\ \{L_a, L_b\}_0 &= \varepsilon_{abc} L_c, & \qquad \{L_a, P_b\}_0 &= \varepsilon_{abc} P_c, \\ \{H, \mathbf{P}_{\mathrm{tot}}\}_0 &\equiv \mathbf{0}, & \qquad \{H, \mathbf{L}_{\mathrm{tot}}\}_0 &\equiv \mathbf{0}, & \qquad \{K_a, H\}_0 &= P_a. \end{aligned} \end{equation} \end{propbox} \leanrefright{bargmann_so3}{BargmannAlgebra.poisson_L_X_L_Y}{07_BargmannAlgebra_Theorems.lean} \leanrefright{bargmann_charge}{BargmannAlgebra.poisson_P_G_same}{07_BargmannAlgebra_Theorems.lean} \pythonrefright{py_core}{core/invariants.py, diagnostics/metrics.py}{Float64 / Float128} \begin{proofbox}[Proof of Proposition~\ref{prop:lie_bargmann_closure}] By the canonical Darboux brackets $\{q_{i,a}, p_{j,b}\}_0 = \delta_{ij}\delta_{ab}$ and $\{q_{i,a}, q_{j,b}\}_0 = \{p_{i,a}, p_{j,b}\}_0 = 0$, the kinematic brackets follow from direct componentwise contraction of the Levi-Civita tensor $\varepsilon_{abc}$: $\{K_a, P_b\}_0 = \sum_{i,j} m_i \{q_{i,a}, p_{j,b}\}_0 = \bigl(\sum_i m_i\bigr)\delta_{ab} = \Mtot \delta_{ab}\mathbf{1}$, $\{L_a, K_b\}_0 = \sum_{i,c,d} \varepsilon_{acd} m_i q_{i,c} \delta_{db} = \varepsilon_{abc} K_c$, $\{L_a, P_b\}_0 = \varepsilon_{abc} P_c$, and $\{L_a, L_b\}_0 = \varepsilon_{abc} L_c$. For the Hamiltonian $H = T(\mathbf{p}) + V(\mathbf{q})$, since $T(\mathbf{p}) = \sum_{i=1}^N \frac{\|\mathbf{p}_i\|^2}{2m_i}$ depends solely on momenta and $V(\mathbf{q})$ depends solely on pairwise Euclidean distances $r_{ij} = \|\mathbf{q}_i - \mathbf{q}_j\|$, skew-symmetry $\nabla_{\mathbf{q}_i} r_{ij} = -\nabla_{\mathbf{q}_j} r_{ij}$ and collinearity $(\mathbf{q}_i - \mathbf{q}_j) \times \nabla_{\mathbf{q}_i} r_{ij} = \mathbf{0}$ yield $\sum_{i=1}^N \nabla_{\mathbf{q}_i} V = \mathbf{0}$ and $\sum_{i=1}^N \mathbf{q}_i \times \nabla_{\mathbf{q}_i} V = \mathbf{0}$ as in~\eqref{eq:1.4_orthogonality}, while $\sum_{i=1}^N \mathbf{p}_i \times \nabla_{\mathbf{p}_i} T = \sum_{i=1}^N m_i^{-1}(\mathbf{p}_i \times \mathbf{p}_i) = \mathbf{0}$; hence $\{H, \mathbf{P}_{\mathrm{tot}}\}_0 = \sum_i \nabla_{\mathbf{q}_i} V \equiv \mathbf{0}$ and $\{H, \mathbf{L}_{\mathrm{tot}}\}_0 \equiv \mathbf{0}$. Finally, since $\{K_a, V(\mathbf{q})\}_0 = 0$, one obtains $\{K_a, H\}_0 = \{K_a, T(\mathbf{p})\}_0 = \sum_{i=1}^N m_i \frac{\partial T}{\partial p_{i,a}} = \sum_{i=1}^N p_{i,a} = P_a$. \end{proofbox} In extended phase space-time, the non-autonomous secular Galilean boost $\mathbf{G}(t) = \mathbf{K} - t \mathbf{P}_{\mathrm{tot}}$ completes the ten classical Galilean first integrals, satisfying $\{P_a, G_b\}_0 = -\Mtot \delta_{ab} \mathbf{1}$ and vanishing total convective derivative along the Hamiltonian flow: \begin{equation} \frac{\mathrm{d}\mathbf{G}}{\mathrm{d}t} = \frac{\partial \mathbf{G}}{\partial t} + \{\mathbf{G}, H\}_0 = -\mathbf{P}_{\mathrm{tot}} + \mathbf{P}_{\mathrm{tot}} \equiv \mathbf{0}. \end{equation}