% --- Section 3.1 : Canonical Synodic Formulation of the PCR3BP ---
\subsection{Canonical Synodic Formulation of the PCR3BP}
\label{sec:pcr3bp_canonical}



To benchmark our normal-cone-inspired splitting methodology on a concrete celestial mechanics problem before tackling the full $N$-body collision stratification, we examine the planar circular restricted three-body problem (PCR3BP). This classic Hamiltonian system provides a rigorous, low-dimensional testing ground possessing isolated binary singularities, non-trivial Coriolis coupling, and an exact first integral.

\subsubsection{Rotating Synodic Coordinates and Canonical Phase Space}
Let two massive bodies (the primaries) of masses $m_1 = 1-\mu$ and $m_2 = \mu$ (with $0 < \mu \le 1/2$) revolve in circular Keplerian orbits about their common center of mass with normalized mean motion $\Omega = 1$ and separation distance $R = 1$. A third massless body moves in the orbital plane of the primaries. 

In the co-rotating (synodic) reference frame, the primaries are fixed at coordinates $\mathbf{q}_{(1)} = (-\mu, 0)$ and $\mathbf{q}_{(2)} = (1-\mu, 0)$. Denoting the Cartesian coordinates of the third body by $\mathbf{q} = (x, y) \in \R^2$ and its conjugate canonical momenta by $\mathbf{p} = (p_x, p_y) \in \R^2$, the phase space is $T^*(\R^2 \setminus \{\mathbf{q}_{(1)}, \mathbf{q}_{(2)}\})$, equipped with the Darboux symplectic form:
\begin{equation}
    \omega_0 = \mathrm{d}x \wedge \mathrm{d}p_x + \mathrm{d}y \wedge \mathrm{d}p_y.
\end{equation}

\subsubsection{Unsoftened Synodic Hamiltonian and Jacobi Energy Surface}
The unsoftened autonomous Hamiltonian governing the third body in canonical variables is:
\begin{equation}
    H_{\mathrm{syn}}(x, y, p_x, p_y) = \frac{1}{2}(p_x^2 + p_y^2) + (y p_x - x p_y) - \frac{1-\mu}{r_1} - \frac{\mu}{r_2},
\end{equation}
where $r_1 = \sqrt{(x + \mu)^2 + y^2}$ and $r_2 = \sqrt{(x - 1 + \mu)^2 + y^2}$ denote the pairwise Euclidean distances to the primaries. The kinematic coupling term $y p_x - x p_y = -\mathbf{q} \times \mathbf{p}$ accounts for the Coriolis acceleration of the rotating synodic frame. The kinematic velocity vector $\mathbf{v} = (\dot{x}, \dot{y})$ is related to canonical momenta by $\mathbf{v} = \mathbf{p} - \mathbf{\Omega} \times \mathbf{q} = (p_x + y, p_y - x)$, yielding the squared velocity relation $v^2 = \|\dot{\mathbf{q}}\|^2 = (p_x^2 + p_y^2) + 2(y p_x - x p_y) + (x^2 + y^2)$.

\begin{propbox}[Autonomous Conservation of the Jacobi First Integral]{prop:pcr3bp_jacobi}
Following Jacobi~\cite{jacobi1836} and Szebehely~\cite[Chapter~1, \S1.2--\S1.5]{szebehely1967theory}, because the synodic Hamiltonian $H_{\mathrm{syn}}$ is autonomous ($\partial H_{\mathrm{syn}} / \partial t = 0$), the Jacobi energy observable:
\begin{align}
    C_J(x, y, p_x, p_y) &= -2 H_{\mathrm{syn}} \nonumber \\
    &= -(p_x^2 + p_y^2) - 2(y p_x - x p_y) + 2\left(\frac{1-\mu}{r_1} + \frac{\mu}{r_2}\right) \label{eq:jacobi_canonical} \\
    &= 2 U_{\mathrm{eff}}(x, y) - (\dot{x}^2 + \dot{y}^2), \label{eq:jacobi_velocity}
\end{align}
where $U_{\mathrm{eff}}(x, y) = \frac{1}{2}(x^2 + y^2) + \frac{1-\mu}{r_1} + \frac{\mu}{r_2}$ is the augmented effective potential, is an exact first integral along every unsoftened phase space trajectory:
\begin{equation}
    \frac{\mathrm{d}C_J}{\mathrm{d}t} = \{C_J, H_{\mathrm{syn}}\}_0 \equiv 0.
\end{equation}
\end{propbox}
\leanrefright{pcr3bp_hamiltonian}{PhaseState.jacobi_conservation_identity, PhaseState.coriolis_step_symplectic}{01_PCR3BPHamiltonian_Theorems.lean}

\begin{proofbox}[Proof of Proposition~\ref{prop:pcr3bp_jacobi}]
Hamilton's equations for $H_{\mathrm{syn}}$ with respect to $\omega_0 = \mathrm{d}x \wedge \mathrm{d}p_x + \mathrm{d}y \wedge \mathrm{d}p_y$ read $\dot{x} = \partial_{p_x} H_{\mathrm{syn}} = p_x + y$, $\dot{y} = \partial_{p_y} H_{\mathrm{syn}} = p_y - x$, $\dot{p}_x = -\partial_x H_{\mathrm{syn}} = p_y - \partial_x V_{\mathrm{grav}}$, and $\dot{p}_y = -\partial_y H_{\mathrm{syn}} = -p_x - \partial_y V_{\mathrm{grav}}$, where $V_{\mathrm{grav}}(x,y) = -\frac{1-\mu}{r_1} - \frac{\mu}{r_2}$. Differentiating $C_J = -2 H_{\mathrm{syn}} = 2(x p_y - y p_x) - (p_x^2 + p_y^2) - 2 V_{\mathrm{grav}}(x,y)$ along any solution curve yields:
\[
    \frac{\mathrm{d}C_J}{\mathrm{d}t}
    = 2(\dot{x} p_y - \dot{y} p_x + x \dot{p}_y - y \dot{p}_x) - 2(p_x \dot{p}_x + p_y \dot{p}_y) - 2(\partial_x V_{\mathrm{grav}}\,\dot{x} + \partial_y V_{\mathrm{grav}}\,\dot{y}) = -2\{H_{\mathrm{syn}}, H_{\mathrm{syn}}\}_0 \equiv 0.
\]
Substituting $p_x = \dot{x} - y$ and $p_y = \dot{y} + x$ into~\eqref{eq:jacobi_canonical} expands $-(p_x^2+p_y^2) - 2(y p_x - x p_y) = (x^2+y^2) - (\dot{x}^2+\dot{y}^2)$, recovering~\eqref{eq:jacobi_velocity}.
\end{proofbox}

The energy hypersurfaces $C_J = \mathrm{const}$ determine Hill's forbidden regions whose topological bifurcations at the Lagrange equilibrium points $L_1, \dots, L_5$ govern capture, resonance transitions, and escape~\cite{szebehely1967theory}. In classical celestial mechanics, Szebehely~\cite{szebehely1967theory} adds the gauge constant $\frac{1}{2}\mu(1-\mu)$ to the effective potential $\Omega(x,y) = U_{\mathrm{eff}}(x,y) + \frac{1}{2}\mu(1-\mu)$, which doubles to $+\mu(1-\mu)$ in $C_{\mathrm{Sz}} = 2\Omega(x,y) - (\dot{x}^2 + \dot{y}^2) = C_J + \mu(1-\mu)$, leaving the physical trajectories identical.

% (removed forced pagebreak)
To synthesize this synodic Hamiltonian architecture into an operational reference, Figure~\ref{fig:ch03_pcr3bp_synodic_polarization} illustrates the co-rotating synodic base geometry alongside the canonical cotangent polarization, and we summarize the resulting structure in the dedicated specification below:

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

% =========================================================================
% LEFT PANEL: Rotating Synodic Configuration Plane R^2 \ {q_(1), q_(2)}
% =========================================================================
\begin{scope}[shift={(0,0)}]
    % Background Frame
    \draw[draw=BorderGray, fill=BadgeGray, rounded corners=3pt, line width=0.7pt] (-3.5,-2.3) rectangle (3.7,2.6);
    \node[font=\footnotesize\bfseries, text=NavyBlue] at (0.1, 2.25) {Synodic Configuration Base $\R^2 \setminus \{\mathbf{q}_{(1)}, \mathbf{q}_{(2)}\}$};

    % Axes
    \draw[->, draw=DarkSlate!70, thin] (-3.0,0) -- (3.2,0) node[right, font=\scriptsize] {$x$};
    \draw[->, draw=DarkSlate!70, thin] (-0.5,-1.5) -- (-0.5,1.65) node[above, font=\scriptsize] {$y$};

    % Primaries q_(1) and q_(2)
    \fill[AcademicPurple] (-1.9,0) circle (3.2pt);
    \node[below=5pt, font=\scriptsize, text=AcademicPurple] at (-1.9,0) {$\mathbf{q}_{(1)} = (-\mu,0)$};
    \node[above left=2pt and -2pt, font=\scriptsize\bfseries, text=AcademicPurple] at (-1.9,0) {$m_1 = 1-\mu$};

    \fill[AcademicPurple] (1.6,0) circle (2.2pt);
    \node[below=5pt, font=\scriptsize, text=AcademicPurple] at (1.6,0) {$\mathbf{q}_{(2)} = (1-\mu,0)$};
    \node[above right=2pt and -2pt, font=\scriptsize\bfseries, text=AcademicPurple] at (1.6,0) {$m_2 = \mu$};

    % Test Particle q = (x,y)
    \coordinate (Q) at (0.35, 1.15);
    \fill[NavyBlue] (Q) circle (2.2pt);
    \node[above=3pt, font=\scriptsize\bfseries, text=NavyBlue] at (Q) {$\mathbf{q} = (x,y)$};

    % Distances r_1 and r_2
    \draw[dashed, thick, draw=LeanBlue] (-1.9,0) -- (Q) node[pos=0.52, above left=0pt, font=\scriptsize, text=LeanBlue] {$r_1$};
    \draw[dashed, thick, draw=LeanBlue] (1.6,0) -- (Q) node[pos=0.48, above right=0pt, font=\scriptsize, text=LeanBlue] {$r_2$};

    % Synodic Rotation Indicator Omega = 1
    \draw[->, thick, draw=AcademicPurple] (2.6,1.0) arc (15:140:0.42);
    \node[font=\scriptsize\bfseries, text=AcademicPurple] at (2.35, 1.62) {$\Omega = 1$};

    % Bottom Formula Note
    \node[font=\scriptsize, text=DarkSlate] at (0.1, -1.9) {$\Delta = \{\mathbf{q}_{(1)}, \mathbf{q}_{(2)}\}, \quad \mathbf{v} = \mathbf{p} - \boldsymbol{\Omega} \times \mathbf{q}$};
\end{scope}

% =========================================================================
% RIGHT PANEL: Cotangent Bundle Polarization T*(R^2 \ Delta)
% =========================================================================
\begin{scope}[shift={(7.6,0)}]
    % Background Frame
    \draw[draw=BorderGray, fill=BadgeGray, rounded corners=3pt, line width=0.7pt] (-3.5,-2.3) rectangle (3.7,2.6);
    \node[font=\footnotesize\bfseries, text=NavyBlue] at (0.1, 2.25) {Normal-Cone Splitting on $T^*(\R^2 \setminus \Delta)$};

    % Vertical Cotangent Fibers T*_q R^2
    \draw[draw=BorderGray, thick] (-2.1,-1.25) -- (-2.1,1.45);
    \draw[draw=BorderGray, thick] (0.0,-1.25) -- (0.0,1.45);
    \draw[draw=BorderGray, thick] (2.1,-1.25) -- (2.1,1.45);
    \node[above, font=\scriptsize, text=DarkSlate] at (-2.1,1.45) {$T^*_{\mathbf{q}_0}\R^2$};
    \node[above, font=\scriptsize, text=DarkSlate] at (2.1,1.45) {$T^*_{\mathbf{q}_1}\R^2$};

    % Base Manifold Curve
    \draw[thick, draw=NavyBlue] (-2.8,-1.25) -- (2.9,-1.25);
    \node[below=2pt, font=\scriptsize, text=NavyBlue] at (-2.1,-1.25) {$\mathbf{q}_0$};
    \node[below=2pt, font=\scriptsize, text=NavyBlue] at (1.55,-1.25) {$\mathbf{q}_1 = \mathbf{R}_{-h}(\mathbf{q}_0 + h\mathbf{p}_{1/2})$};

    % Step 1: Half Vertical Kick Phi_{h/2}^normal
    \coordinate (P0) at (-2.1, -0.4);
    \coordinate (P05) at (-2.1, 0.85);
    \fill[DarkSlate] (P0) circle (2pt) node[left=2pt, font=\scriptsize] {$(\mathbf{q}_0, \mathbf{p}_0)$};
    \draw[->, very thick, draw=PythonGreen] (P0) -- (P05)
        node[midway, left=2pt, font=\scriptsize\bfseries, text=PythonGreen] {$\Phi_{h/2}^{\mathrm{normal}}$};
    \fill[PythonGreen] (P05) circle (2pt);

    % Step 2: Exact Linear Coriolis-Kinematic Flow Psi_h^base in Sp(4, R)
    \coordinate (P15) at (2.1, 0.85);
    \draw[->, very thick, draw=AcademicPurple] (P05) to[out=20, in=160]
        node[midway, above=2pt, font=\scriptsize\bfseries, text=AcademicPurple] {$\Psi_h^{\mathrm{base}} \in \operatorname{Sp}(4,\R)$}
        node[midway, below=2pt, font=\scriptsize, text=DarkSlate] {Exact Coriolis Flow} (P15);
    \fill[AcademicPurple] (P15) circle (2pt);

    % Step 3: Half Vertical Kick Phi_{h/2}^normal
    \coordinate (P1) at (2.1, -0.35);
    \draw[->, very thick, draw=PythonGreen] (P15) -- (P1)
        node[midway, right=2pt, font=\scriptsize\bfseries, text=PythonGreen] {$\Phi_{h/2}^{\mathrm{normal}}$};
    \fill[NavyBlue] (P1) circle (2pt) node[right=2pt, font=\scriptsize] {$(\mathbf{q}_1, \mathbf{p}_1)$};

    % Bottom Conservation Invariant
    \node[font=\scriptsize, text=DarkSlate] at (0.1, -1.9) {$\mathcal{S}_2(h)^*\omega_0 = \omega_0, \quad \widetilde{C}_{J,h} = C_J + \mathcal{O}(h^2)$};
\end{scope}

\end{tikzpicture}%
}
\caption{Left: Co-rotating synodic configuration space $\R^2 \setminus \{\mathbf{q}_{(1)}, \mathbf{q}_{(2)}\}$ of the planar circular restricted three-body problem (PCR3BP) with primaries $m_1 = 1-\mu$ and $m_2 = \mu$. Right: Canonical normal-cone Strang splitting $\mathcal{S}_2(h) = \Phi_{h/2}^{\mathrm{normal}} \circ \Psi_h^{\mathrm{base}} \circ \Phi_{h/2}^{\mathrm{normal}}$ on $T^*(\R^2 \setminus \Delta)$, alternating vertical cotangent fiber momentum translations with the exact linear Coriolis-kinematic groupoid flow $\Psi_h^{\mathrm{base}} \in \operatorname{Sp}(4,\R)$.}
\label{fig:ch03_pcr3bp_synodic_polarization}
\end{figure}

\begin{conormalbox}{Planar Circular Restricted 3-Body Problem (PCR3BP)}
    \textbf{Geometric Base \& Normal Fiber:} The physical system describes the coplanar motion of an infinitesimal third body under the Newtonian gravitational attraction of two massive primaries revolving in circular orbits about their common barycenter. The state is parameterized by the synodic position vector $\mathbf{q} = (x, y) \in \R^2 \setminus \{\mathbf{q}_{(1)}, \mathbf{q}_{(2)}\}$, the conjugate canonical momentum vector $\mathbf{p} = (p_x, p_y) \in \R^2$, and the normalized mass parameter $\mu = m_2/(m_1 + m_2) \in (0, 1/2]$.
    
    The geometric phase space is the canonical cotangent bundle of the punctured plane equipped with the standard Darboux 2-form:
    \begin{equation}
        \mathcal{M} = T^*\bigl(\R^2 \setminus \{\mathbf{q}_{(1)}, \mathbf{q}_{(2)}\}\bigr), \qquad \omega_0 = \mathrm{d}\mathbf{q} \wedge \mathrm{d}\mathbf{p} = \mathrm{d}x \wedge \mathrm{d}p_x + \mathrm{d}y \wedge \mathrm{d}p_y.
    \end{equation}
    The collision locus is the discrete codimension-2 subvariety formed by the two primary masses:
    \begin{equation}
        \Delta = \{\mathbf{q}_{(1)}, \mathbf{q}_{(2)}\} = \{(- \mu, 0), (1 - \mu, 0)\}, \qquad \mathcal{I}_\Delta = \langle y, \, (x + \mu)(x - 1 + \mu) \rangle \subset \R[x, y].
    \end{equation}
    \vspace{1pt}

    \textbf{Hamiltonian Polarization:} The synodic kinetic-potential model originating in the synodic integral of \textbf{Jacobi (1836)}~\cite{jacobi1836}, cast in canonical Hamiltonian form by \textbf{Poincaré (1892)}~\cite{poincare1892}, and systematized by \textbf{Szebehely (1967)}~\cite{szebehely1967theory} identifies with the canonical degree-$0$ cotangent lift of the deformation to the normal cone, polarizing the autonomous synodic energy into:
    \begin{equation}
        H_{\mathrm{syn}}(\mathbf{q}, \mathbf{p}) = \underbrace{\frac{1}{2}\|\mathbf{p}\|^2 + (y p_x - x p_y)}_{H_{\mathrm{base}}} + \underbrace{\left(-\frac{1-\mu}{r_1} - \frac{\mu}{r_2}\right)}_{H_{\mathrm{normal}}}.
    \end{equation}
    \begin{itemize}[leftmargin=1.5em, itemsep=1.5pt, topsep=1.5pt]
        \item \emph{Base Kinematic Hamiltonian $H_{\mathrm{base}}$}: Quadratic kinetic-Coriolis generator governing uniform rotating frame kinematics and free inertial drift in $\operatorname{Sp}(4, \R)$, integrated exactly without truncation error via the linear groupoid action $\Psi_h^{\mathrm{base}}(\mathbf{q}, \mathbf{p}) = (\mathbf{R}_{-h}(\mathbf{q} + h\mathbf{p}), \mathbf{R}_{-h}\mathbf{p})$.
        \item \emph{Normal Potential Hamiltonian $H_{\mathrm{normal}}$}: Singular Newtonian gravitational potential depending strictly on base configuration $\mathbf{q} \in \R^2 \setminus \Delta$, integrated as an exact abelian transverse kick $\mathbf{p} \mapsto \mathbf{p} - h \nabla V_{\mathrm{grav}}(\mathbf{q})$ along the vertical cotangent fibers $T^*_{\mathbf{q}} \R^2$, corresponding to the cotangent lift of the normal cone coordinates.
    \end{itemize}
    \vspace{1pt}

    \textbf{Preserved Invariants:} Under the symmetric Strang splitting $\mathcal{S}_2(h) = \Phi_{h/2}^{\mathrm{normal}} \circ \Psi_h^{\mathrm{base}} \circ \Phi_{h/2}^{\mathrm{normal}}$, the symplectic 2-form is preserved identically at any finite step $h$, while the continuous flow preserves $C_J = -2H_{\mathrm{syn}}$ exactly and the discrete propagator preserves the modified shadow Jacobi integral $\widetilde{C}_{J,h} = -2\widetilde{H}_{\mathrm{syn},h} = C_J + \mathcal{O}(h^2)$ without secular drift:
    \begin{equation}
        \mathcal{S}_2(h)^* \omega_0 = \omega_0, \qquad \{C_J, H_{\mathrm{syn}}\}_0 \equiv 0 \implies \frac{\mathrm{d}C_J}{\mathrm{d}t} \equiv 0,
    \end{equation}
    guaranteeing exact conservation of the phase-space Liouville volume $\mathrm{d}\mathbf{q} \wedge \mathrm{d}\mathbf{p}$ and bounded Jacobi energy oscillations with zero secular drift.\par\vspace{3pt}
    \leanrefright{pcr3bp_hamiltonian}{PhaseState.jacobi_conservation_identity, PhaseState.coriolis_step_symplectic}{01_PCR3BPHamiltonian_Theorems.lean}
    \pythonrefright{py_interfaces}{interfaces/pcr3bp_canonical_solver.py}{Planar PCR3BP Float64}
\end{conormalbox}