1 Introduction↩︎

We study polar orbits, i.e., spatial periodic orbits with large vertical components and close approaches to a primary, in the circular restricted three-body problem (CR3BP). Such orbits are of special interest in space mission design. For instance, the halo orbit family is widely used in missions for its advantageous out-of-plane geometry. Near rectilinear halo orbits (NRHOs) [1] offer repeated close approaches to one celestial body while maintaining continuous visibility of the other, a feature underlying their use in NASA’s planned Gateway mission, [2]. NRHOs have also been proposed as science orbits for plume sampling missions to Saturn’s icy moon Enceladus, [3], [4], which is one of the primary motivations for this study.

We will use a Hamiltonian description of the CR3BP in a uniformly rotating frame, \[H =\frac{1}{2}|\mathbf{p}|^2 +p_1q_2-p_2q_1 -\frac{1-\mu}{|\mathbf{q}-(-\mu,0,0)|} -\frac{\mu}{|\mathbf{q}-(1-\mu,0,0)|}.\] Here, \(\mathbf{q}=(q_1, q_2, q_3)\) and \(\mathbf{p}=(p_1, p_2, p_3)\) denote the position and momentum, and \(\mu\in [0,1]\) is the mass ratio of the primaries. We focus on the case \(\mu\) close to \(0\) and on motion near the light primary at \(\mathbf{q}=(1 - \mu, 0, 0)\). By rescaling the coordinates and taking the limit \(\mu\to 0\), one obtains the Hill’s problem (see Section 2.3 for details), whose Hamiltonian is \[H =\frac{1}{2}|\mathbf{p}|^2 +p_1q_2-p_2q_1 -\frac{1}{|\mathbf{q}|} +\frac{1}{2} |\mathbf{q}|^2 -\frac{3}{2}q_1^2 .\]

From an abstract perspective, the limit cases of families of periodic orbits are of particular importance. For example, many periodic orbits in the CR3BP arise as perturbations of periodic orbits in the limiting case of the Hill’s problem, which Hénon referred to as generating orbits, [5]. In the study of polar orbits, an important example is the vertical collision orbit in Hill’s problem, a trajectory confined to the vertical \(q_3\)-axis that reaches a maximum height before ending in a collision, see Figure 6. After regularization, this orbit becomes periodic and serves as an important generating orbit for nearby polar orbits.¹ Recent work by Aydin and Batkhin, [7], provides a comprehensive study of bifurcating families of multiple covers of the vertical collision orbit within Hill’s problem, along with an overview of earlier results, [8]–[10]. One aim of this work is to reveal how this vertical collision orbit continues into the CR3BP and to organize the bifurcation structure of related families of polar orbits.

A detailed numerical study of near-collision periodic orbits requires a regularization framework compatible with high-precision continuation. The Kustaanheimo–Stiefel (KS) regularization, [11], as described for the numerical setting in [1], has traditionally been used for this purpose. In this study, we adopt an alternative approach based on Moser regularization. This framework offers several advantages: it is applicable in arbitrary dimensions (e.g., planar or spatial CR3BP), preserves periodicity of orbits between the regularized and unregularized systems, and provides a direct correspondence of Floquet multipliers (Lemma 4), simplifying stability analysis. We provide a first, complete implementation of Moser regularization for numerical continuation and stability analysis of periodic orbits. Our implementation treats the dynamics as a Hamiltonian system with constraints, whose framework we review in Appendix 7. In particular, the regularized phase space appears as a \(6\)-dimensional constrained submanifold embedded in an \(8\)-dimensional Euclidean space.

Based on this setup, we investigate the bifurcation structure of polar orbits near the light primary in the Moser regularized CR3BP and in Hill’s problem. We identify the pitchfork, period-doubling, and period-tripling bifurcations of the vertical collision family, the associated branching families, and how these structures persist or change under perturbation from Hill’s problem to the CR3BP. Bifurcation graphs are computed for several representative systems, including the Saturn–Enceladus, Earth–Moon, and Copenhagen problems. A summary of the periodic families and their key properties (including symmetries and Conley–Zehnder indices) is provided in the conclusion (Table 2) for convenient reference. Our results provide a global understanding of key orbit families relevant to mission design, such as halo and butterfly families, while providing new observations on their behavior near collision. In particular, this work complements and extends the foundational analysis of halo and near rectilinear halo orbits, [12]–[15], and further connects these findings to the symplectic geometry perspective through a Hamiltonian formulation and the use of Moser regularization.

Notably, we describe the near-collision bifurcation structure of four families of three-lobed polar orbits, which we refer to as tri-fly orbits due to their shape. In Hill’s problem, these families are associated with a period-tripling bifurcation of the vertical collision orbit, while in the Saturn–Enceladus and Earth–Moon systems they are organized by touch-and-go and pitchfork bifurcations which occur near the triple cover of the \(L_2\) halo orbit. Determining the bifurcation mechanism involved in these orbits for small mass parameters, such as in Saturn–Enceladus, requires both high precision and regularization. Members of these families include low-altitude south pole flyby trajectories in the Enceladus system (see Figure 26).

For detailed analysis, we compute the Conley–Zehnder indices of the periodic orbit families in this study; the Conley–Zehnder index is a symplectic invariant that remains constant along non-degenerate continuous families and whose change implies bifurcations. This index provides both a diagnostic tool for numerical detection of bifurcations and a topological classification of families, and has been recently used for numerical bifurcation studies of CR3BP and related systems, [7], [16], [17]. We give details and computation methods in Appendix 9.

The paper is organized as follows. The main text focuses on the underlying celestial mechanics, while the contents of a purely mathematical nature are collected in the appendix. A brief outline is given below. Section 2 introduces the rescaled Hamiltonian formulation for CR3BP and its limiting relation to Hill’s problem. Section 3 presents the numerical framework for Moser regularization and introduces the vertical collision orbit. Section 4 presents the numerical results. In Appendix 7, we provide details on Hamiltonians with constraints and the restricted phase space. We follow that in Appendix 8 with additional details of differential correction in Moser regularization. Appendix 9 discusses computations of Conley–Zehnder indices. Finally, Appendix 10 provides numerical data for the orbit families discussed in this work.

2 Dynamical Models and Symplectic Scaling Framework↩︎

2.1 The Circular Restricted Three-Body Problem (CR3BP)↩︎

We consider the spatial circular restricted three-body problem (CR3BP), which models the motion of an infinitesimal particle in three-dimensional Euclidean space, under the gravitational influence of two massive primary bodies in circular motion. The mass ratio of the primaries is denoted by the parameter \(\mu \in [0, 1]\).

We choose a rotating coordinate frame fixing the two primaries along the \(x\)-axis, the first primary at \(\mathbf{e} = (-\mu, 0, 0)\) and the second at \(\mathbf{m} = (1 - \mu, 0, 0)\), and the barycenter of the system at the origin. For our choice of rotation direction, the Hamiltonian of the CR3BP is given by \[H_\mu(\mathbf{q}, \mathbf{p}) = \frac{1}{2}|\mathbf{p}|^2 + p_1q_2 - p_2q_1 - \frac{\mu}{|\mathbf{q} - \mathbf{m}|} - \frac{1 - \mu}{|\mathbf{q} - \mathbf{e}|},\] where the position vector is \(\mathbf{q} = (q_1, q_2, q_3) \in \mathbb{R}^3\), the momentum vector is \(\mathbf{p} = (p_1, p_2, p_3) \in \mathbb{R}^3\). The motion of the infinitesimal object is described by the Hamilton equations, given by \[\dot{\mathbf{q}} = \partial_{\mathbf{p}} H, \quad \dot{\mathbf{p}} = - \partial_{\mathbf{q}} H.\] We define the Hamiltonian vector field of \(H_\mu\) as \(X_\mu = (\partial_{\mathbf{p}} H, -\partial_{\mathbf{q}} H)\). Due to energy conservation, each trajectory is confined to a five-dimensional energy surface \(\Sigma_c = H^{-1}(c)\). The constant \(c\) is referred to as the Jacobi energy, and the corresponding Jacobi constant is given by \(\Gamma = -2c\). The system has five equilibrium points, called Lagrange points. The three collinear points \(L_1, L_2, L_3\) lie on the \(q_1\)-axis, while the two equilateral points \(L_4, L_5\) form equilateral triangles with the primaries (see Figure 1).

Our main focus is on the dynamics near the light primary, including the region around the collinear points \(L_1\) and \(L_2\), in the regime where \(\mu\) is small. This setting is relevant for most planet–moon systems, and is of particular interest for the exploration of icy moons such as Europa orbiting Jupiter and Enceladus orbiting Saturn. Table 1 lists the mass ratios \(\mu\) of interest in this study, with values from NASA’s Jet Propulsion Laboratory solar system dynamics website: https://ssd.jpl.nasa.gov.

2.2 Hill’s Problem and Symmetries↩︎

The Hill’s problem, [18], is a limiting case of the CR3BP where the first primary is assumed to be much heavier than the second, and the infinitesimal object to be very close to the light primary. The Hamiltonian for the Hill’s problem is given as follows. \[\label{eqn:hill-hamiltonian} \hat{H}_0(\mathbf{q}, \mathbf{p}) = \frac{1}{2}|\mathbf{p}|^2 - \frac{1}{|\mathbf{q}|} + p_1 q_2 - p_2 q_1 + \frac{1}{2}|\mathbf{q}|^2 - \frac{3}{2}q_1^2.\tag{1}\] We denote the Hamiltonian by \(\hat{H}_0\) because in the next section we will introduce a rescaled Hamiltonian \(\hat{H}_\mu\) for the CR3BP, which converges to \(\hat{H}_0\) as \(\mu\to 0\). Among the five equilibrium (Lagrange) points in the CR3BP, only the two collinear points \(L_1\) and \(L_2\) remain in the Hill’s problem (see Figure 1). Like the CR3BP, the Hill’s problem is non-integrable, [19], so numerical methods are essential for studying its detailed dynamics, [20].

The Hill’s problem comes with additional symmetries compared to the CR3BP. In the CR3BP, the Hamiltonian is invariant under the following reversing symmetries, which reverse the direction of the Hamiltonian flow: \[\begin{align} r_y(q_1, q_2, q_3, p_1, p_2, p_3) &= (q_1, -q_2, q_3, -p_1, p_2, -p_3), \tag{2} \\ r_{yz}(q_1, q_2, q_3, p_1, p_2, p_3) &= (q_1, -q_2, -q_3, -p_1, p_2, p_3). \tag{3} \end{align}\] In the Hill’s problem, two additional reversing symmetries appear: \[\begin{align} r_x(q_1, q_2, q_3, p_1, p_2, p_3) &= (-q_1, q_2, q_3, p_1, -p_2, -p_3), \tag{4} \\ r_{xz}(q_1, q_2, q_3, p_1, p_2, p_3) &= (-q_1, q_2, -q_3, p_1, -p_2, p_3). \tag{5} \end{align}\] These symmetries generate a symmetry group of order \(8\), [21], including the pure symmetry \(r_z = r_y \circ r_{yz} = r_x \circ r_{xz}\) which preserves the direction of the flow: \[r_z(q_1, q_2, q_3, p_1, p_2, p_3) = (q_1, q_2, -q_3, p_1, p_2, -p_3).\] Note that the additional symmetries 4 and 5 also appear in the CR3BP in the special case \(\mu=1/2\).

2.3 Rescaled CR3BP and the Hill’s Limit↩︎

We now introduce a rescaled CR3BP Hamiltonian \(\hat{H}_\mu\) which converges to Hill’s problem as \(\mu \to 0\). The derivation follows the descriptions in [22], [23] and proceeds through two coordinate transformations. First, we translate the origin to the position of the light primary: \[T(\mathbf{q}, \mathbf{p}) = (\mathbf{q}, \mathbf{p}) + \mathbf{m}_\mu, \quad \mathbf{m}_\mu = (1 - \mu, 0, 0, 0, 1 - \mu, 0).\] Next, we apply the symplectic scaling from [24] \[\varphi(\mathbf{q}, \mathbf{p}) = (\mu^{1/3} \mathbf{q}, \mu^{1/3} \mathbf{p}),\] which rescales the Hamilton equations by a constant factor; the rescaled CR3BP Hamiltonian is then \[\hat{H}_\mu(\mathbf{q}, \mathbf{p}) = \mu^{-2/3}\left( (H_\mu \circ T)(\varphi(\mathbf{q}, \mathbf{p})) + (1 - \mu) + \frac{(1 - \mu)^2}{2} \right),\] which expands to \[\label{eqn:rescaled-CR3BP-expanded} \hat{H}_\mu(\mathbf{q}, \mathbf{p}) = \frac{1}{2}|\mathbf{p}|^2 - \frac{1}{|\mathbf{q}|} + p_1 q_2 - p_2 q_1 - \frac{1 - \mu}{\mu^{2/3}} \left( \frac{1}{\sqrt{(\mu^{1/3} q_1 + 1)^2 + \mu^{2/3}(q_2^2 + q_3^2)}} + \mu^{1/3} q_1 - 1 \right).\tag{6}\]

We use the hat (\(\,\hat{}\,\)) notation to indicate the Hamiltonian in scaled coordinates, and \(h\) to denote the corresponding energy levels. The energies of the two systems are related by \[c = \mu^{2/3} h - (1 - \mu) - \frac{(1 - \mu)^2}{2},\] so that the energy surface \(\hat{H}_\mu = h\) corresponds to \(H_\mu = c\) under the transformation \(\varphi\circ T\) (see Figure 2). As \(\mu \to 0\), we have \(c \to -3/2\), and the Jacobi constants converge to \(\Gamma = 3\).

In the rescaled Hamiltonian 6 , the term in parentheses is \(O(\mu^{2/3})\), and \(\hat{H}_\mu\) is analytic in \(\nu=\mu^{1/3}\) at \(\mu = 0\). Thus, the limit \(\mu \to 0\) recovers the Hamiltonian for the Hill’s problem 1 .

This rescaling is useful because it provides a perturbative framework where solutions of the Hill’s problem continue into the CR3BP. While the expression 6 is not numerically practical due to the singular denominator \(\mu^{2/3}\), we provide in Appendix 6 an equivalent formulation suitable for numerical computations.

2.4 Non-degeneracy of Periodic and Symmetric Orbits↩︎

In the most general setting, we consider a \(T\)-periodic solution \(\gamma(t)\) of an autonomous Hamiltonian \(H\), where \(\gamma\) lies on the energy level \(H = h\). Let \(\Sigma=H^{-1}(h)\) denote the energy hypersurface containing \(\gamma\). Let \(\psi:S\to S\) be the Poincaré return map on a surface of section \(S\subset \Sigma\), with \(\psi(\gamma(0)) = \gamma(T)\).

We define the appropriate non-degeneracy condition for such orbits as follows.

This condition can be checked numerically using the eigenvalues of the monodromy matrix \(\Phi = d\varphi^T(\gamma(0))\), where \(\varphi^t\) denotes the time \(t\)-flow of the Hamiltonian vector field \(X_H\). These eigenvalues are referred to as the Floquet multipliers.

We now introduce a suitable non-degeneracy condition for symmetric periodic orbits. Let \(r\) be a reversing symmetry of the Hamiltonian, satisfying \[H \circ r = H, \qquad r \circ \varphi^t = \varphi^{-t} \circ r,\] where \(\varphi^t\) denotes the Hamiltonian flow. We consider the case where \(r\) is an involution, i.e.,\(r^2 = \mathrm{id}\). The fixed-point set \[\mathrm{Fix}(r) = \{ x \in \mathbb{R}^6 \mid r(x) = x \}\] has half the dimension of the phase space. An \(r\)-symmetric periodic orbit of period \(T\) satisfies \[r(\gamma(t)) = \gamma(-t), \qquad \text{so that} \quad \gamma(0),\, \gamma(T/2) \in \mathrm{Fix}(r).\]

Let \(S\subset \Sigma\) be a surface of section such that \(r(S)=S\) and \(\mathrm{Fix}(r)\cap \Sigma \subset S\), and assume \[\psi(\gamma(0)) = \gamma(T/2).\] At each \(x \in \mathrm{Fix}(r)\cap \Sigma\), the tangent space decomposes into the eigenspaces of \(dr\): \[T_x S = E_1 \oplus E_{-1}, \qquad \left(dr|_S\right)|_{E_{\pm 1}} = \pm \mathrm{id},\] where \(E_1\) corresponds to directions tangent to \(\mathrm{Fix}(r)\cap \Sigma\).

2.5 Perturbation of Periodic Orbits from Hill’s Problem to CR3BP↩︎

Applying Poincaré’s perturbation theory of periodic orbits, we describe how non-degenerate periodic orbits of the Hill’s problem persist in the CR3BP for small mass ratios. A closely related approach is the method of generating solutions studied in the works of Perko, Hénon, Bruno, and Varin, [5], [26], [27], which continues planar symmetric orbits of Hill’s problem to the CR3BP within the perturbative framework. Here, we state a general perturbative result (Theorem 2) on fixed energy surfaces, with leading-order approximations which will be used for numerical continuation of these solutions.

2.6 Numerical Continuation with respect to Mass Parameter↩︎

The leading-order approximations of Theorem 2 provide initial guesses for numerically continuing periodic orbits when the mass ratio is small. This continuation can be carried out from Hill’s problem to the CR3BP, or vice versa, or between CR3BP systems with different small mass ratios. We make the latter precise in the following corollary.

This scaling provides a practical method for predicting periodic orbits at a nearby mass ratio \(\mu = \mu_{\text{prev}} + \Delta\mu\) from a known orbit at \(\mu_{\text{prev}}\). The Jacobi energy \(c\) is updated via the intermediary Hill energy \(h\), computed by \[h = \mu_{\text{prev}}^{-\frac{2}{3}} \left( c_{\text{prev}} + (1 - \mu_{\text{prev}}) + \frac{(1 - \mu_{\text{prev}})^2}{2} \right), \quad c = \mu^{\frac{2}{3}} h - (1 - \mu) - \frac{(1 - \mu)^2}{2}.\] The initial point is rescaled as \[X - m_\mu = (X_{\text{prev}} - m_{\mu_{\text{prev}}}) \cdot \left( \frac{\mu}{\mu_{\text{prev}}} \right)^{1/3},\] and the period remains unchanged to leading order: \(T = T_{\text{prev}}\). The resulting prediction is then refined using a shooting method with Newton-type correction that ensures periodicity at the updated parameters. This process is summarized in Algorithm 3. We use this continuation procedure at several stages of our computations to transfer known periodic orbits from the CR3BP or Hill’s problem to different mass ratios, providing initial data for orbit continuation. It is also used in the computation of bifurcation surfaces, in particular the bifurcation graphs of halo and vertical collision orbits over the range \(\mu \in [0,0.5]\) in Section 4.4.

The correction at fixed Jacobi energy is implemented using a single-shooting scheme on a chosen Poincaré section. For clarity, we describe the case of the section \(q_2=0\), but an analogous scheme on the section \(p_3=0\) is used in Section 4.5.

Let \(q_1,q_3,p_1,p_3\) be the coordinates of a predicted initial condition and \(T\) be a predicted period. The missing momentum \(p_2\) is obtained by solving the quadratic equation \(H(\mathbf{q},\mathbf{p})=c\) for \(p_2\). We also compute its derivatives \(\frac{\partial p_2}{\partial q_1},\cdots,\frac{\partial p_2}{\partial p_3}\). The periodicity condition is expressed as finding zeros of the map \[F(q_1, q_3, p_1, p_3, T) = \pi_{(q_1, q_2, q_3, p_1, p_3)} \circ \varphi^T(q_1, 0, q_3, p_1, p_2, p_3) - (q_1, 0, q_3, p_1, p_3)\] where \(\varphi^t\) denotes the Hamiltonian flow and \(\pi\) is the projection onto the listed coordinates². Starting from an initial guess \(X_0=(q_1, q_3, p_1, p_3, T)\), we apply Newton iterations \[X_{i+1} = X_i - [DF(X_i)]^{-1} F(X_i)\] until convergence, where a vector \(X_i\) is considered to be a periodic orbit if \(|F(X_i)|\) is sufficiently small. The computation of the Jacobian \(DF(X_i)\) involves the linearized flow \(d\varphi^T\), which is obtained by integrating the first-order variational equations along the trajectory.

Note that for symmetric (and doubly symmetric) periodic orbits, we can exploit the standard simplification of integrating only half (or a quarter) of the period. The reduction for symmetric orbits is used throughout our computations. For brevity, we present the details of this simplification only in the Moser regularization setting in Section 3.4.

3 Moser Regularization↩︎

In this section, we develop a numerical framework for detection and continuation of periodic orbits in the Moser regularized CR3BP. These methods will be essential for the numerical study of families of collision and near-collision spatial orbits and their bifurcation graphs.

3.1 Moser’s Transformation↩︎

We describe Moser’s transformation of coordinates \((\mathbf{q}, \mathbf{p})\in \mathbb{R}^3\times \mathbb{R}^3\) to the cotangent bundle \(T^*S^3\) of the unit \(3\)-sphere: \[T^*S^3 = \{(\boldsymbol{\xi}, \boldsymbol{\eta}) \in \mathbb{R}^4 \times \mathbb{R}^4 \, | \, |\boldsymbol{\xi}|^2 = 1,\;\boldsymbol{\xi}\cdot \boldsymbol{\eta} = 0 \}.\] To define the transformation, we first switch the roles of position and momentum coordinates through the switch map \((\mathbf{x}, \mathbf{y}) = (-\mathbf{p}, \mathbf{q})\). Moser’s transformation is then given by: \[\label{eqn:moser95transf} \begin{align} &\xi_0 = \frac{|\mathbf{x}|^2 - 1}{|\mathbf{x}|^2 + 1}, \qquad \xi_k = \frac{2x_k}{|\mathbf{x}|^2 + 1}, \qquad k=1,2,3. \\ &\eta_0 = \mathbf{x} \cdot \mathbf{y}, \qquad \eta_k = \frac{|\mathbf{x}|^2+1}{2} y_k -(\mathbf{x} \cdot \mathbf{y}) x_k, \qquad k=1,2,3. \end{align}\tag{7}\] Note that here the momentum \(\mathbf{x}\in \mathbb{R}^3\) is mapped to the unit \(3\)-sphere \(S^3\subset\mathbb{R}^4\) using the inverse stereographic projection, where the singularity \(p=\infty\) corresponds to \(\xi_0 = 1\), i.e., the north pole \(N\) of the sphere. Figure 4 shows a schematic picture of the transformation. The inverse of this transformation is given as follows: \[q_k = y_k = \eta_0 \xi_k + (1 - \xi_0)\eta_k, \qquad p_k = -x_k = -\frac{\xi_k}{1 - \xi_0}, \qquad k=1,2,3.\]

3.2 Moser Regularization of CR3BP and Hill’s Problem↩︎

In Moser’s work, [6], the transformation 7 was introduced to regularize the singularity of the Kepler problem. Here, we consider the regularization of Hamiltonian systems of the form \[H(\mathbf{q},\mathbf{p}) = \frac{1}{2}|\mathbf{p}|^2 + p_1q_2 - p_2q_1 - \frac{g}{|\mathbf{q}|} + \tilde{V}(\mathbf{q}),\] at the singularity at \(\mathbf{q}=0\), where \(g\neq 0\) is a constant and \(\tilde{V}\) is analytic in a neighborhood of \(\mathbf{q}=0\). For the CR3BP and Hill’s problem, the corresponding values of \(g\) and \(\tilde{V}\) are given by \[\begin{align} \text{CR3BP:}\quad g &= \mu, \quad \tilde{V}(\mathbf{q}) = -\frac{1-\mu}{\sqrt{(q_1+1)^2 + q_2^2 + q_3^2}} - (1-\mu) q_1 - \frac{(1-\mu)^2}{2}, \notag \\ \text{Hill's problem:}\quad g &= 1, \quad \tilde{V}(\mathbf{q}) = \frac{1}{2}|\mathbf{q}|^2 - \frac{3}{2} q_1^2 . \end{align}\] In the CR3BP case, the Hamiltonian has been translated by \[q_1 \mapsfrom q_1 - (1-\mu), \qquad p_2 \mapsfrom p_2 - (1-\mu),\] so that the collision singularity is located at the origin. For a general description of the Moser regularization in Stark-Zeeman systems, see [28].

To define the regularized Hamiltonian, we first set \[K(\mathbf{q},\mathbf{p}) = (H - c)|\mathbf{q}| = \left( \frac{1}{2}(|\mathbf{p}|^2 + 1) - (c + 1/2) + p_1q_2 - p_2q_1 + \tilde{V}(\mathbf{q}) \right)|\mathbf{q}| - g.\] The energy level set \(H=c\) now corresponds to \(K=0\), and the dynamics on this energy surface is reparametrized by \(X_K = |\mathbf{q}|X_H\). In the regularized coordinates \((\boldsymbol{\xi}, \boldsymbol{\eta})\), we have \[K(\boldsymbol{\xi}, \boldsymbol{\eta}) = \left( 1 - (1 - \xi_0)(c + 1/2) + (1 - \xi_0)(\xi_2\eta_1 - \xi_1\eta_2) + (1 - \xi_0)\tilde{V}(\boldsymbol{\xi}, \boldsymbol{\eta}) \right)|\boldsymbol{\eta}| - g.\] Finally, we smooth this to define the Moser regularized Hamiltonian \[Q:T^*S^3\to \mathbb{R},\qquad Q(\boldsymbol{\xi}, \boldsymbol{\eta}) = \frac{1}{2g}(K(\boldsymbol{\xi}, \boldsymbol{\eta}) + g)^2.\] The energy level set \(K=0\) corresponds to \(Q=\frac{g}{2}\), where the dynamics is preserved via \(X_Q = X_K\). The following summarizes the result.

The Moser regularization gives dynamics on the restricted phase space \(T^*S^3\subset\mathbb{R}^8\) defined by the two constraints: \[f_1 = |\boldsymbol{\xi}|^2 - 1 = 0, \qquad f_2 = \boldsymbol{\xi}\cdot \boldsymbol{\eta} = 0.\] The corresponding Hamiltonian vector field on \(T^*S^3\) is obtained as follows.

In Appendix 7, we give a general formulation of restricted phase spaces with \(2k\) constraints, and the proof of the formula for \(X_Q\) in this generalized setting.

3.3 Floquet Multipliers and Bifurcations↩︎

A key advantage of the Moser regularization is that the Floquet multipliers of periodic orbits can be computed in a simple way. Let \(\gamma\) be a periodic orbit of period \(\tau\) in the unregularized system, and let \(\tilde{\gamma}\) be the corresponding periodic orbit in the regularized system on \(T^*S^3\subset \mathbb{R}^8\).

Consider the monodromy matrix \(\tilde{\Phi} = d\tilde{\varphi}^\tau(\tilde{\gamma}(0))\) of the flow of \(X_Q\) in \(\mathbb{R}^8\). Because this flow has two additional integrals \(f_1\) and \(f_2\), we can take a basis adapted to these two integrals and to the tangent space of \(T^*S^3\), for which the monodromy matrix takes the block form \[\left( \begin{array}{c|c} \text{id}_2 & 0 \\\hline * & S \end{array} \right),\] where \(S\) represents the derivative of the flow restricted to \(T^*S^3\) (see Lemma 2.3 of [25]). The matrix \(S\) is symplectic and is conjugate to the monodromy matrix of \(\gamma\) in the original coordinates via the derivative of Moser’s transformation 7 . Thus, the nontrivial (i.e., not equal to \(1\)) Floquet multipliers of \(\gamma\) and \(\tilde{\gamma}\) are the same.

If more than four eigenvalues are equal to \(1\), the orbit is degenerate and bifurcations may occur. Let \(X_s\) be a smooth family of vector fields, and let \(s \mapsto \gamma_s\) be a smooth one-parameter family of periodic orbits of \(X_s\). The linearized return map \(R_s\) then depends smoothly on \(s\), and its eigenvalues can be chosen to vary continuously with \(s\): the latter claim follows from Theorem 5.1 and 5.2 in [29]. In Hamiltonian systems, these eigenvalues follow the symplectic eigenvalue theorem: if \(\lambda\) is an eigenvalue, then so are \(\bar \lambda\), \(1/\lambda\) and \(1/\bar \lambda\). They either lie on the real axis \(\mathbb{R}\) or on the unit circle \(S^1 \subset \mathbb{C}\), or appear as complex quadruples. Hence, changes in stability may occur when a pair of eigenvalues merges at \(1\) (tangent bifurcation), at \(-1\) (period-doubling bifurcation), or when two pairs collide on \(S^1\) and leave the unit circle as a complex quadruple (Krein collision), [30]. Below we describe the local bifurcations relevant to our study, following the terminology of [31].

See Figure 5 for a visual description. Note that the modified secondary Hopf bifurcation does not involve a change in stability.

None

3.4 Differential Correction and Continuation↩︎

We present a differential correction scheme formulated directly in the Moser regularized coordinates 7 , which allows periodic orbits to be continued smoothly through collision configurations. The Poincaré section is chosen as \(\xi_3 = 0\), corresponding to \(p_3 = 0\) in the unregularized coordinates. This section naturally captures turning points in the \(z\)-direction, such as the periapsis and apoapsis of halo orbits, and contains the fixed-point locus of the symmetries. All numerical integrations are performed using the high-order Taylor integrator implemented in the heyoka library, [32].

We choose the vector of free variables \[X = (\xi_1, \xi_2, \eta_1, \eta_2, \tau),\] where \(\tau\) is the (predicted) period of the orbit. Given \(X\), we consecutively determine the \(\xi_0\) and \(\eta_0\) variables based on the constraints \(|\boldsymbol{\xi}|^2 = 1\) and \(\boldsymbol{\xi} \cdot \boldsymbol{\eta} = 0\), respectively. Then, \(\eta_3\) is determined by solving the energy constraint \(Q = g/2\), for which we use numerical root finding based on Brent’s method, [33]. Define the constraint function: \[F(X) = p \circ \varphi^\tau(\xi_0, \xi_1, \xi_2, 0, \eta_0, \eta_1, \eta_2, \eta_3) - (\xi_1, \xi_2, 0, \eta_1, \eta_2),\] where \(\varphi^\tau\) denotes integration of the flow of \(X_Q\) for time \(\tau\), and \(p\) is the projection to \((\xi_1, \xi_2, \xi_3, \eta_1, \eta_2)\) coordinates. The differential correction scheme is given by Newton’s method as follows: \[X_{i + 1} = X_i - DF(X_i)^{-1} F(X_i), \qquad i=0,1,2,\cdots.\] We iterate until convergence, where a vector \(X_i\) is considered to be a periodic orbit if \(|F(X_i)|\) is sufficiently small: we chose the threshold to be less than \(10^{-10}\) for most computations, and switched to \(10^{-30}\) for high-precision computations, which were needed for the investigation of touch-and-go bifurcations in Saturn–Enceladus. To compute the Jacobian matrix \(DF\), we evaluate the linearized flow, as well as symbolically compute the derivatives of the dependent variables \(\xi_0, \eta_0\) and \(\eta_3\) with respect to the free variables. See Appendix 8 for the details of this derivation.

We describe the correction scheme for symmetric periodic orbits with respect to the reflection \(r_y\) across the plane \(y=0\); the formulation for the reflection \(r_x\) across \(x=0\) is analogous. For this symmetry, the fixed point locus is given by \(q_2 = p_1 = p_3 = 0\), which in Moser coordinates corresponds to \(\xi_1 = \xi_3 = \eta_0 = \eta_2 = 0\). We can therefore reduce the free variable vector to \[X_{\textrm{sym}} = (\xi_2, \eta_1, \tau).\] The remaining variables \(\xi_0\) and \(\eta_3\) are determined as before. The constraint function is then defined as \[F_{\textrm{sym}}(X_{\textrm{sym}}) = p_{\textrm{sym}} \circ \varphi^{\tau/2}(\xi_0, 0, \xi_2, 0, 0, \eta_1, 0, \eta_3),\] where \(p_{\textrm{sym}}\) projects to \((\xi_1, \xi_3, \eta_0)\) coordinates. This gives a reduced \(3\)-dimensional correction scheme. The halved integration time also improves numerical robustness. For doubly symmetric orbits, the integration interval can be further reduced to a quarter period; however, this additional reduction was not required for the computations in this study.

For the numerical continuation of periodic orbits, we apply pseudo-arclength continuation in the unregularized system, and switch to natural parameter continuation with respect to the Hamiltonian energy in the Moser-regularized system as trajectories approach the collision singularity. Initial periodic orbits are obtained either from the JPL Three-Body Periodic Orbit Catalog (https://ssd.jpl.nasa.gov/tools/periodic_orbits.html) or from internal data computed using the cell-mapping method, [34], and are transferred between systems and mass ratios using Algorithm 3. For several families studied in Sections 4.5 and 4.7, branch switching at local bifurcation points is applied. We refer to [35] for the formulation of the continuation schemes and the branch-switching method.

3.5 Vertical Collision Orbits↩︎

Vertical collision orbits in the Hill’s problem arise from dynamics confined to the \(z\)-axis, i.e., the submanifold \(q_1 = q_2 = p_1 = p_2 = 0\), which is invariant under the Hamiltonian flow. Restricting to the \(z\)-axis, we get a reduced \(2\)-dimensional Hamiltonian system on the \((q_3, p_3)\)-plane: \[\dot{q_3} = \frac{\partial \hat{H}_0}{\partial p_3} = p_3, \qquad \dot{p_3} = -\frac{\partial \hat{H}_0}{\partial q_3} = - \frac{q_3}{|q_3|^3} - q_3.\] This system is singular at \(q_3=0\), corresponding to collision with the primary at the origin. For initial conditions with \(q_3\neq 0\), the trajectory reaches a collision after a finite time, which is computed in [10].

To continue these trajectories through collision, we apply Moser regularization. In this framework, the dynamics restricted to the vertical subspace corresponds to the invariant set \[\begin{align} & \xi_1 = \xi_2 = \eta_1 = \eta_2 = 0, \\ & \xi_0^2 + \xi_3 ^2 = 1, \\ & \xi_0\eta_0 + \xi_3\eta_3 = 0, \end{align}\] giving a one degree-of-freedom Hamiltonian system. The regularized Hamiltonian on this subset takes the form \[Q = \frac{1}{2}\left( 1 - (1 - \xi_0)(c + 1/2) + (1 - \xi_0)\cdot \frac{1}{2}(\eta_0\xi_3 + (1 - \xi_0)\eta_3)^2 \right)^2(\eta_0^2 + \eta_3^2).\] This function is smooth across the collision locus \(p_3=\pm \infty\), which corresponds to \(\xi_0 = 1\). Thus, the collision trajectories extend to smooth periodic orbits in the regularized phase space, as shown in Figure 6. We call these orbits vertical collision orbits.

For each energy level \(\hat{H}_0 = h\), there exists a unique pair of vertical collision orbits. These can be characterized explicitly using the Poincaré section \(p_3 = 0\), at which the corresponding position \(q_3\) is given by \[\label{eqn:collsion95h95to95q3} \hat{H}_0 = - \frac{1}{|q_3|} + \frac{1}{2}q_3^2 = h.\tag{8}\] Solving for \(q_3\), we get two solutions for each \(h\) with opposite signs, corresponding to two orbits reflected about the origin. We refer to the orbit with \(q_3 > 0\) as the northern vertical collision orbit, and the other with \(q_3 < 0\) southern. The positive solution \(q_{3,\textrm{max}}\) represents the maximal height of these orbits, and increases monotonically with energy \(h\), starting from zero and increasing to infinity.

4 Numerical Results↩︎

4.1 Bifurcations of Vertical Collision Orbits and Halo Orbits in Hill’s Problem↩︎

The bifurcation diagram of the northern vertical collision family is shown in Figure 7. In all bifurcation diagrams presented in this work, including Figure 7 and those in subsequent sections, the labeled integers along each branch denote the associated Conley–Zehnder indices of the periodic orbits. To highlight regions of linear stability, we additionally shade the corresponding energy intervals for which all Floquet multipliers lie on the unit circle with a green background.

As a one-parameter family parameterized by the Hamiltonian energy \(\hat{H}_0=h\), the family begins as a linearly stable periodic orbit possessing two pairs of elliptic Floquet multipliers. As \(h\) increases, five distinct bifurcation events are detected:

(1) First period-doubling bifurcation (\(h \approx -1.02\)) A period-doubling bifurcation occurs as an elliptic pair of Floquet multipliers transitions to a hyperbolic pair, at which point a doubly symmetric butterfly family of periodic orbits emerges.

(2) First pitchfork bifurcation (\(h \approx -0.85\)) The remaining elliptic pair transitions to a positive hyperbolic pair at a point of a pitchfork bifurcation. This bifurcation gives rise to two branches of \(r_y\)-symmetric orbits, which appear in symmetry with respect to the reflection \(r_x\). These correspond to the northern \(L_1\) and \(L_2\) halo families.

(3) Second pitchfork bifurcation (\(h \approx 0.04\)) A second pitchfork bifurcation occurs, producing two branches of \(r_x\)-symmetric orbits that appear in symmetry with respect to the reflection \(r_y\). After perturbation to the CR3BP, the \(r_x\)-symmetry of the orbits is broken, giving families of non-symmetric orbits which connect to vertical families of orbits emanating from the \(L_4\) and \(L_5\) Lagrange points. These correspond to the W4/W5 families described in [36].

(4) Second period-doubling bifurcation (\(h \approx 0.09\)) A second period-doubling bifurcation occurs, giving rise to the “moth” family of doubly symmetric orbits. This family also appears in the bifurcation graph presented in [7], where its continuation is shown to connect to the planar direct family of periodic orbits (family \(g\)).

(5) Secondary Hopf bifurcation (\(h \approx 0.11\)) Finally, a secondary Hopf bifurcation is observed, in which two pairs of Floquet multipliers collide on the unit circle and move off into the complex plane. This transition marks the onset of strong instability, with the modulus of one pair increasing far beyond unity while the other decreases toward zero.

The branches described in (1)-(4) are illustrated in Figure 8. Table 3 lists numerical data for selected orbits on these branches. We note that an anonymous reviewer pointed out that some of the results in this section were also described by Lidov and Liakhova in [8], however, we were unable to access this source.

The bifurcation graphs of the two halo branches are plotted together in Figure 7. The halo families originate from the planar Lyapunov family at \(h \approx -2.0\). As the energy increases, each halo branch undergoes two successive period-doubling bifurcations followed by a fold bifurcation. In the literature, this fold bifurcation is often associated with the onset of the near rectilinear halo orbit (NRHO) regime, [37]. Beyond the fold, the family undergoes a narrow interval of linear stability that terminates at another period-doubling bifurcation. Continuing further in \(h\), the family connects with the vertical collision orbit at the pitchfork bifurcation near \(h \approx -0.85\). In Figure 9 we visualize the northern \(L_1\) halo family as it bifurcates from the planar Lyapunov orbit and terminates in the vertical collision orbit.

4.2 Perturbation of Vertical Collision Orbits to CR3BP↩︎

We study the continuation of the vertical collision orbits of Hill’s problem under perturbation to CR3BP, with respect to the parameter \(\nu=\mu^{1/3}\). The discussion here extends the similar description given for the KS regularization in [15].

Fix a Hamiltonian energy \(\hat{H}_\nu = h\) and let \(X_\nu\) denote the Hamiltonian vector field on the Moser regularized energy surface. Denote by \(\varphi_\nu(x, t)\) the time-\(t\) flow of \(X_\nu\). Let \(x_0\) be a point of the vertical collision orbit of Hill’s problem on the Poincaré section \(\xi_3=0\), with period \(t_0\). We seek perturbed solutions \((x_\nu, t_\nu)\) satisfying \[\varphi_\nu(x_\nu, t_\nu) - x_\nu = 0.\] The constraints on \(x_\nu\) include the energy surface \(\hat{H}_\nu = h\), the section \(\xi_3=0\), and the constraints \(f_1 = |\boldsymbol{\xi}|^2 = 1\) and \(f_2 = \boldsymbol{\xi}\cdot\boldsymbol{\eta} = 0\) from the Moser regularization.

Linearizing with respect to \(\nu\) at \(\nu=0\) gives \[\label{eqn:perturbed-linearized} (d\varphi_0 - \text{id}) \Delta x + X_0(x_0) \Delta t + \Delta \varphi= 0,\tag{9}\] where \(\Delta x = \frac{d x_\nu}{d\nu}\big|_{\nu=0}\), \(\Delta t = \frac{d t\nu}{d\nu}\big|_{\nu=0}\), and \[\Delta \varphi= \int_0^{t_0} \frac{\partial X_0}{\partial \nu}(\varphi_0(x_0, s)) ds.\]

Since the vertical collision orbit lies in the invariant “vertical” subspace \(\xi_1=\xi_2=\eta_1=\eta_2=0\), we project 9 onto the “horizontal” coordinates \((\xi_1, \xi_2, \eta_1, \eta_2)\). In this projection, the \(\Delta t\) term vanishes, yielding \[A \begin{bmatrix} \Delta \xi_1 \\ \Delta \xi_2 \\ \Delta \eta_1 \\ \Delta \eta_2 \end{bmatrix} + b = 0,\] where \[A = \frac{\partial (\pi_{(\xi_1, \xi_2, \eta_1, \eta_2)} \circ \varphi_0)}{\partial (\xi_1, \xi_2, \eta_1, \eta_2)} - \text{id}, \quad b = \pi_{(\xi_1, \xi_2, \eta_1, \eta_2)} (\Delta \varphi).\] The matrix \(A\) is invertible except at two energy levels \(h_1\approx -0.85\) and \(h_2\approx 0.04\) corresponding to degeneracies. Numerical evaluation shows that the solutions \(\Delta x = -A^{-1}b\) diverges near \(h=h_1\), but is continuous near \(h=h_2\), with \(\Delta \xi_1 = \Delta \eta_2 = 0\), implying perturbation to symmetric orbits.

Next, we consider the \(r_y\)-symmetric continuation of the vertical collision orbits. We look for perturbed \(r_y\)-symmetric solutions \((x_\nu, t_\nu)\) satisfying \[\pi_{(\xi_1, \xi_3, \eta_0, \eta_2)} \circ \varphi_\nu (x_\nu, t_\nu) = 0\] while restricting to the fixed point locus \[x_\nu \in \text{Fix}(r_y) = \{(\xi,\eta) \,|\, \xi_1 = \xi_3 = \eta_0 = \eta_2 = 0\}.\] The linearized perturbation projected onto \((\xi_2, \eta_1)\) coordinates give \[\label{eqn:perturbed-symmetric-linearized-projected} A_{\text{sym}} \cdot \begin{bmatrix} \Delta \xi_2 \\ \Delta \eta_1 \\ \end{bmatrix} + b_{\text{sym}} = 0,\tag{10}\] where \[A_{\text{sym}} = \frac{\partial (\pi_{(\xi_1, \eta_2)} \circ \varphi_0)}{\partial (\xi_2, \eta_1)}, \quad b_{\text{sym}} = \pi_{(\xi_1, \eta_2)} (\Delta \varphi).\] Numerical evaluation shows that \(A_{\text{sym}}\) is singular only at \(h=h_1\), whereas the orbit at \(h=h_2\) remains non-degenerate as a symmetric orbit.

The solution to the system 10 at the two intersections of the vertical collision orbit with \(\xi_3=0\) – the apoapsis (\(q_3\) maximal) and periapsis (\(q_3=0\)) – is shown in Figure 10. Since \(\Delta \xi_2\neq 0\), we have \(|\xi_0| < 1\) after perturbation, indicating that vertical collision orbits are continued to non-collision orbits. Using the relation \(\Delta q_1 = (1 - \xi_0) \Delta \eta_1\), the displacement of the unregularized coordinate \(q_1\) has the same sign as \(\Delta \eta_1\). We summarize the observations as follows:

1. The vertical collision orbit is degenerate at two energy levels, \(h=h_1\) and \(h=h_2\). As an \(r_y\)-symmetric orbit, the orbit at \(h=h_2\) is non-degenerate. Hence, the vertical collision orbit continues for \(h\neq h_1\) to an \(r_y\)-symmetric orbit after perturbation to the CR3BP by Theorem 2.

2. For \(h < h_1\), the vertical collision orbit continues under perturbation to a symmetric, non-collision \(L_2\) halo orbit, with apoapsis displaced outward (\(q_1 > 1-\mu\)) and periapsis inward (\(q_1 < 1-\mu\)) for small \(\mu\).

3. For \(h > h_1\), the vertical collision orbit continues under perturbation to a symmetric, non-collision \(L_1\) halo orbit, with apoapsis displaced inward (\(q_1 < 1-\mu\)) and periapsis outward (\(q_1 > 1-\mu\)) for small \(\mu\).

4.3 Near-Collision Behavior of Halo Orbits in the Saturn–Enceladus System↩︎

As an explicit case of the perturbation from Hill’s problem to the CR3BP, we examine the near-collision behavior of halo orbits in the Saturn–Enceladus system (\(\mu = 1.901109735892602 \times 10^{-7}\)), computed via continuation in the Moser-regularized system.

The bifurcation diagrams of the northern \(L_1\) and \(L_2\) halo families are shown in Figure 11, with numerical data for selected orbits given in Table 4. Compared with Hill’s problem, the breaking of the \(r_x\)-symmetry in the CR3BP removes the pitchfork bifurcation connecting the vertical collision orbit and the halo families. Consequently, the low-energy branch of the vertical collision orbit connects to the \(L_2\) halo family, while the high-energy branch connects to the \(L_1\) halo family. Note that this behavior is similar to the bifurcation structure of the planar direct periodic orbit families in the CR3BP and Hill’s problem (families \(g\) and \(g'\)), as presented in [16] and earlier in [38]. We refer to the near-vertical continuations of the halo branches as collisional halo orbits.

Figure 12 illustrates the transitions from planar to collisional of the \(L_1\) and \(L_2\) halo families, respectively. Prior to the collisional branch, the orbits appear in a near symmetric configuration with respect to the plane \(x = 1 - \mu\), due to the approximate symmetry \(r_x\) continued from the Hill’s problem. Beyond the NRHO regime, the orbits become increasingly confined in the \(yz\)-plane and acquire a near vertical geometry. In this regime, the \(z\)-coordinate of the apoapsis grows rapidly for the \(L_1\) family and decreases for the \(L_2\) family, while both retain their vertical structure.

4.4 Bifurcation Surface of Halo Families from \(\mu=0\) to \(\mu=0.5\)↩︎

To understand the global evolution of the halo families and their relation to the vertical collision orbits, we compute the bifurcation diagrams of the \(L_1\) and \(L_2\) halo orbits for mass ratios \(\mu \in [0, 0.5]\) of the rescaled CR3BP Hamiltonian \(\hat{H}_\mu\). The Jacobi energy is varied over the range \(h \in [-1, 0.4]\) to construct a two-parameter bifurcation surface that continuously connects Hill’s problem to the equal-mass case.

The resulting bifurcation surface is shown in Figure 13. For small mass ratios, the \(L_1\) family shows two distinct fold bifurcations, bounding an interval which includes the linearly stable regime corresponding to the near-rectilinear halo orbits (NRHOs). A qualitative transition occurs at a critical mass ratio \(\mu^*\) between \(\mu=0.05\) and \(\mu=0.06\), where the two folds disappear and the Jacobi energy becomes a monotone function along the family. Thus, for most planet–moon systems, the bifurcation structure of the \(L_1\) family includes the same set of degeneracies, with exceptions such as the Pluto–Charon system (\(\mu \approx 0.12\)). The critical value \(\mu^*\) corresponds to the upper limit \(\mu^* \approx 0.0573\) identified by Howell, [13], for the existence of a linearly stable region of the \(L_1\) halo family. Meanwhile, the third degeneracy at higher energy, corresponding to the pitchfork bifurcation of the vertical collision family described in Section 4.1, persists for \(\mu\in[0,0.5]\).

4.4.1 Earth–Moon Case↩︎

The bifurcation diagram of halo orbits in the Earth–Moon system is shown in Figure 14, with numerical data for selected orbits given in Table 4. We observe the same topological structure as in the Saturn–Enceladus system, with the graphs displaying the same sequence of bifurcations but at different parameter scales. The perturbation from Hill’s problem is more pronounced, with less symmetry between \(L_1\) and \(L_2\) halo families and more distorted bifurcation diagrams.

4.4.2 Copenhagen Case \((\mu=0.5)\)↩︎

We examine the equal-mass case \(\mu=0.5\), known as the Copenhagen problem, named after the series of studies carried out at the Copenhagen Observatory under the direction of Strömgren, [39]. This case possesses additional symmetries: reflection across both the plane \(x=0\) 4 and the \(y\)-axis 5 leaves the equations of motion invariant.

The bifurcation diagram for the Copenhagen problem is shown in Figure 15, with numerical data for selected orbits given in Table 4. In contrast to the small-\(\mu\) regime, the \(L_1\) family no longer goes through fold bifurcations, and its Jacobi energy increases monotonically. The \(L_2\) family retains a fold bifurcation, but the order and position of the period-doubling bifurcations are modified.

Figure 16 shows the \(L_1\) and \(L_2\) halo families of the Copenhagen problem. In this case, the \(L_1\) halo orbits are doubly symmetric with respect to both reflections \(r_y\) across the \(xz\)-plane and \(r_{xz}\) across the \(y\)-axis. The behavior near collision is visibly distinct from the small mass ratio case: the \(z\)-coordinate of the apoapsis of the \(L_1\) family increases consistently from the start of the continuation, while for the \(L_2\) family, the orbits contract smoothly toward collision with the light primary.

4.5 W4/W5 families in Hill’s Problem and Earth–Moon System↩︎

We study the pair of \(r_x\)-symmetric periodic orbits that bifurcate from the second pitchfork bifurcation of the vertical collision orbit in Hill’s problem at energy \(h \approx 0.04\), as described in Section 4.1. These orbits correspond after continuation to the Earth–Moon system to the W4/W5 families studied in [36]. For consistency, we refer to these orbits in the Hill’s problem also as the W4/W5 families.

Figure 17 shows the bifurcation graph in Hill’s problem, with numerical data for selected orbits of the W5 branch given in Table 3. From the pitchfork bifurcation, the two branches emerge in symmetry in the direction of increasing energy. Along each branch, we observe three period-doubling bifurcations: one pair of Floquet multipliers alternates between negative hyperbolic, elliptic, negative hyperbolic again, and finally elliptic, while the other pair remains positive hyperbolic. The continuation of the W4 family as it bifurcates from the vertical collision orbit is shown in Figure 18.

As explained in Section 2.5, this pitchfork bifurcation persists under perturbation to the CR3BP because the \(r_y\)-symmetry of the pitchfork bifurcation persists. The resulting branches are no longer symmetric orbits in the CR3BP, but instead form two distinct non-symmetric families that connect to the vertical periodic orbits emanating from the equilateral Lagrange points \(L_4\) and \(L_5\).

Figure 19 shows the corresponding bifurcation graph in the Earth–Moon system, which is continued using the Poincaré section \(p_3 = 0\). Table 5 lists numerical data for selected orbits of the W4 branch. As in Hill’s problem, the two branches bifurcate from the \(L_1\) halo family in the direction of increasing energy and pass through three period-doubling bifurcations. As the energy increases further, the apoapsis \(y\)-coordinate starts to decrease until the family connects at energy \(c \approx -0.97\) to the vertical periodic orbit associated with \(L_4\). Figure 20 shows the continuation of the W4 family from the \(L_1\) halo orbit to its termination in the \(L_4\) vertical orbit.

4.6 Period-Doubling Bifurcation and Butterfly Orbit Family↩︎

The butterfly orbit family is a well-known class of periodic orbits cataloged in the JPL Three-Body Periodic Orbit Catalog. These orbits share several characteristics with halo orbits, including a vertical geometry, close proximity to the light primary, and extended coverage of the polar region, [40]. Recently, butterfly orbits in the Saturn–Enceladus system have been proposed as candidate science orbits for a mission to Enceladus due to their coverage of the south polar plume region, which is of high scientific interest because of its active geysers, [41]. In the literature, butterfly orbits have been primarily studied in the Earth–Moon system, where they arise through a period-doubling bifurcation of the \(L_2\) halo family, [14].

By continuation in the Moser regularized rescaled Hamiltonian \(\hat{H}_\mu\), we obtain the following characterization of the butterfly orbits:

1. In the Hill’s problem, the butterfly family arises as a doubly symmetric orbit from a period-doubling bifurcation of the vertical collision orbit at energy level \(h\approx-1.02\).

2. Under perturbation to the CR3BP, the butterfly family emerges as a symmetric orbit from the period-doubling bifurcation of the collisional \(L_2\) halo orbit.

Figure 21 shows the bifurcation diagrams of the butterfly family in Hill’s problem and in the Saturn–Enceladus CR3BP, with numerical data for selected orbits given in Table 3 and Table 6, respectively. In both systems, the butterfly family emerges from a period-doubling bifurcation of either the collisional \(L_2\) halo or vertical collision family, where the Floquet multipliers change from elliptic-elliptic (stable) to elliptic-negative hyperbolic (unstable). The butterfly branch itself begins with elliptic-positive hyperbolic type and evolves in the direction of decreasing Hamiltonian energy. The family then undergoes a fold bifurcation, after which the energy increases. In both systems, the Floquet multipliers transition to elliptic-elliptic type right after the fold, then immediately transition through a secondary Hopf bifurcation to four complex conjugate multipliers.

Figure 22 shows the emergence of the butterfly family in the Saturn–Enceladus system from a period-doubling bifurcation of the collisional \(L_2\) halo orbit. At the bifurcation point, the orbit is nearly vertical, with a diameter in the \(xy\)-plane of order \(10^{-5}\) length units versus a vertical height of 0.00437 LU (1042 km). The perilune radius is approximately 0.05 km, well below the surface of Enceladus (radius 252.1 km). This is significantly closer to collision than in the Earth–Moon system, where the analogous bifurcation occurs at a perilune radius of about 1830 km, [42]. The emerging butterfly orbit preserves the vertical height of the parent halo orbit while developing two symmetric “wings” extending in the \(x\)-direction at apoapsis. After a fold bifurcation, the wings tilt inward and intersect in the \(xz\)-projection.

4.7 Tri-Fly Orbits from Period-Tripling Bifurcations↩︎

We examine the period-tripling bifurcation of the vertical collision and collisional \(L_2\) halo families, which occurs at an energy level lower than the period-doubling bifurcation from which the butterfly family emerges. This period-tripling bifurcation gives rise to a set of four periodic orbit families with symmetric configurations. Examples of each of the four orbits are shown in Figure 23. These orbits exhibit a three-lobed structure in the vertical direction and form a triangular pattern in the \(xy\)-projection. We refer to these as tri-fly orbits to reflect this characteristic geometry. The appearance of these families from the third cover of the vertical collision orbits also appears in [7].

The four tri-fly orbits are characterized by the orientation of their triangular structure in the \(xy\)-projection: “right”, “left”, “up” and “down”. The “right” and “left” orbits (plotted in green) are symmetric with respect to the reflection \(r_y\), while the “up” and “down” orbits (plotted in purple) are \(r_x\)-symmetric in the Hill’s problem but become non-symmetric after perturbation to CR3BP. The “up” and “down” families are exact mirror images of one another with respect to the plane \(y=0\).

Figure 24 shows the bifurcation structure of the tri-fly families in Hill’s problem, the Saturn–Enceladus system, and the Earth–Moon system. In Hill’s problem, all four tri-fly families emerge simultaneously from the triple cover of the vertical collision orbit, with the four branches exhibiting a \(\mathbb{Z}_2 \times \mathbb{Z}_2\) symmetry, [7]. After perturbation to the CR3BP, this symmetry is broken. In the Saturn–Enceladus system, the “up” and “down” families bifurcate via a pitchfork bifurcation from the “left” tri-fly branch, whereas in the Earth–Moon system they bifurcate from the “right” branch. As a result, the four tri-fly families emerge sequentially rather than simultaneously. Table 7 lists numerical data for selected orbits from the branches in the Saturn–Enceladus and Earth–Moon systems.

Figure 25 shows a zoomed view near the period-tripling bifurcation of the \(L_2\) halo family for the Saturn–Enceladus and Earth–Moon systems. In both cases, two branches intersect transversely in the direction of increasing Jacobi energy, with both branches remaining on both sides of the energy at the bifurcation point. This is characteristic of a touch-and-go bifurcation, [43], also referred to as a phantom kiss, [44], in Meyer’s classification of generic bifurcations, [45]. Immediately after the intersection, the “left” tri-fly branch undergoes a fold bifurcation, after which the Jacobi energy decreases. These observations are consistent with the computer-assisted proofs of touch-and-go bifurcations by Walawska and Wilczak, [46], for \(L_2\) halo orbits in the Earth–Moon and Sun-Jupiter systems. Such bifurcations occur on very small energy scales (on the order of \(10^{-12}\) in the Saturn–Enceladus system). To reliably identify the bifurcation type at this scale, we employed for the Saturn–Enceladus system arbitrary-precision arithmetic from MPFR, [47], together with a high-order Taylor integrator implemented in the CAPD library, [48]. The result confirms the presence of the same touch-and-go bifurcation in the Saturn–Enceladus system.

Finally, we note that the “right” tri-fly orbit corresponds to the three-petal tulip-shaped orbit studied in [49] for the Earth–Moon system. The remaining tri-fly families exhibit similar three-lobed spatial structures, differing primarily by orientation, and the corresponding orbits in Saturn–Enceladus involve close flybys of the south pole region of Enceladus (see Figure 26), making them potentially relevant for future mission design. From a dynamical viewpoint, the bifurcation mechanisms identified here provide insight into the organization of these orbit families.

In the Earth–Moon system, we additionally observe a short-lived, non-symmetric “spin” tri-fly branch, shown in brown in Figure 24, which connects the “left” and “right” branches over a narrow parameter range. The brown circle in the bifurcation diagram represents the evolution of the initial condition located at the top of one of the three petals as this point makes a full turn around the orbit. Because the orbit has three petals, this curve traverses the same family evolution three times. Representative orbits from the “left”, “down”, and “spin” branches in the Earth–Moon system are shown in Figure 27.

5 Conclusion↩︎

In this work, we investigated the bifurcation structure of polar orbits near the light primary in the CR3BP using a Hamiltonian formulation together with Moser regularization, using the vertical collision orbit in Hill’s problem as a central organizing structure for the emergence and persistence of orbit families under perturbation. We identified and organized the pitchfork, period-doubling, and period-tripling bifurcations of the vertical collision family, along with the associated halo, butterfly, and tri-fly families, and constructed bifurcation graphs for representative systems including Saturn–Enceladus, Earth–Moon, and the Copenhagen problem. The use of Conley–Zehnder indices provided a consistent topological classification of these families and a diagnostic tool for detecting bifurcations. For convenience of reference, we summarize the principal orbit families, their symmetries, and indices in Table 2.

We conclude by noting several limitations and directions for future work. The present study focuses on small mass ratios and near-collision regimes, and is based primarily on high-precision numerical continuation, while a rigorous analytical treatment of the observed bifurcation structures remains to be developed. Extensions to more general settings, such as elliptic or multi-body problems, are also left for future investigation. Nevertheless, the framework and orbit families identified here provide a foundation for systematic study of near-collision dynamics and for future work integrating symplectic methods, numerical continuation, and mission design, particularly for low-altitude trajectories near small bodies.

6 Numerical Implementation of Rescaled CR3BP↩︎

The rescaled CR3BP Hamiltonian \(\hat{H}_\mu\) in 6 contains a singular denominator \(\mu^{2/3}\). For numerical integration at small \(\mu\), we let \(\nu = \mu^{1/3}\) and rewrite the Hamiltonian in a convenient form using the identity: \[\frac{1}{\sqrt{1 + x}} - 1 + \frac{x}{2} = x^2 \left( \frac{2 + \sqrt{1 + x}}{2\sqrt{1 + x} \, (1 + \sqrt{1 + x})^2} \right), \quad x > -1.\] For \(x = 2q_1 \nu + |\mathbf{q}|^2 \nu^2 > -1\), i.e., away from the first primary located at \(\mathbf{q} = (-\mu^{-\frac{1}{3}}, 0, 0)\), we obtain the numerically convenient form: \[\hat{H}_\mu(\mathbf{q},\mathbf{p}) = \frac{1}{2}|\mathbf{p}|^2 - \frac{1}{|\mathbf{q}|} + p_1q_2 - p_2q_1 + \frac{1-\mu}{2}|\mathbf{q}|^2 - \frac{(1-\mu) (2q_1 + |\mathbf{q}|^2\nu)^2 (2 + \sqrt{1 + 2q_1\nu + |\mathbf{q}|^2\nu^2})}{2\sqrt{1 + 2q_1\nu + |\mathbf{q}|^2\nu^2} \, (1 + \sqrt{1 + 2q_1\nu + |\mathbf{q}|^2\nu^2})^2} .\] This expression is analytic in \(\nu = \mu^{1/3}\), and allows accurate numerical evaluation for small \(\mu\).

7 Restricted phase space↩︎

The Moser regularized Hamiltonian described in Section 3.2 fits into the general setting of Hamiltonian systems with constraints, [50], [51]. Here, we give a brief outline of the method together with a motivating example. We consider a Hamiltonian \(\tilde{H}\) on a phase space \({\mathbb{R}}^{2n}\) with coordinates \((\mathbf{q},\mathbf{p})\), which we want to restrict to constraints \(f_1(\mathbf{q},\mathbf{p})=\ldots=f_{2k}(\mathbf{q},\mathbf{p})=0\). Let \(\omega=\sum_{i=1}^n dp_i\wedge dq_i\) be the standard symplectic form and consider the values \[\label{eqn:poisson-bracket-matrix} A_{ij} = \omega(X_{f_i}, X_{f_j}) = \omega(\nabla f_i, \nabla f_j) = \{f_i, f_j\}, \quad i,j=1,\cdots,2k,\tag{11}\] where \(\{,\}\) denotes the Poisson bracket.

Our goal is to study the dynamics on \(R\) given by the restricted Hamiltonian \(H := \tilde{H}|_{R}\).

8 Details of Differential Correction in Moser Regularization↩︎

We derive the Jacobian matrix used in the differential correction scheme formulated in Moser regularized coordinates (Section 3.4). We consider a Moser regularized Hamiltonian \(Q\) corresponding to the original Hamiltonian \(H\), and apply corrections for a fixed energy level \(H = c\). The procedure described here is valid for general (non-symmetric) periodic orbits; for symmetric cases, the reduced formulation in Section 3.4 can be applied analogously.

We work with the Poincaré section \(\xi_3 = 0\) and define the vector of free variables as \[X = (\xi_1, \xi_2, \eta_1, \eta_2, \tau),\] where \(\tau\) denotes the (predicted) period of the orbit. The remaining variables \(\xi_0\), \(\eta_0\), and \(\eta_3\) are determined from the three constraint equations: \[f_1 = |\boldsymbol{\xi}|^2 - 1 = 0, \qquad f_2 = \boldsymbol{\xi} \cdot \boldsymbol{\eta} = 0, \qquad Q(\boldsymbol{\xi}, \boldsymbol{\eta}) = \frac{g}{2}.\] From the first two constraints we obtain \[\xi_0 = \sqrt{1 - \xi_1^2 - \xi_2^2}, \qquad \eta_0 = \frac{1}{\xi_0}(\xi_1 \eta_1 + \xi_2 \eta_2),\] and \(\eta_3\) is determined by numerically solving the energy constraint \(Q(\xi, \eta) = g/2\) near a predicted value, for instance using Brent’s method, [33].

Periodic orbits satisfy the constraint function \[F(X) = p \circ \varphi^{\tau}(\xi_0, \xi_1, \xi_2, 0, \eta_0, \eta_1, \eta_2, \eta_3) - (\xi_1, \xi_2, 0, \eta_1, \eta_2),\] where \(\varphi^{\tau}\) denotes the flow of the Hamiltonian vector field \(X_Q\) for time \(\tau\), and \(p\) is the projection onto \((\xi_1, \xi_2, \xi_3, \eta_1, \eta_2)\). The correction proceeds via Newton iterations: \[X_{i+1} = X_i - DF(X_i)^{-1} F(X_i), \qquad i = 0, 1, 2, \dots,\] until convergence, i.e., when \(|F(X_i)| < 10^{-10}\). For some high-precision computations in Saturn–Enceladus system in Section 4.7, we used the MPFR library and employed a similar strategy with an order 50 Taylor integrator from the CAPD library with \(|F(X_i)| < 10^{-30}\) as stop criterion.

To compute the Jacobian \(DF\), we first evaluate the derivatives of the dependent variables with respect to the free variables. The expressions are: \[\begin{align} & \frac{\partial \xi_0}{\partial \xi_i} = - \frac{\xi_i}{\sqrt{1 - \xi_1^2 - \xi_2^2}}, \qquad i = 1, 2, \\ & \frac{\partial \eta_0}{\partial \xi_0} = - \frac{\xi_1 \eta_1 + \xi_2 \eta_2}{\xi_0^2}, \qquad \frac{\partial \eta_0}{\partial \xi_i} = -\frac{\eta_i}{\xi_0}, \qquad \frac{\partial \eta_0}{\partial \eta_i} = -\frac{\xi_i}{\xi_0}, \qquad i = 1, 2, \\ & \frac{\partial \eta_3}{\partial \xi_i} = \left( \frac{\partial Q}{\partial \eta_3} \right)^{-1} \frac{\partial Q}{\partial \xi_i}, \qquad \frac{\partial \eta_3}{\partial \eta_i} = \left( \frac{\partial Q}{\partial \eta_3} \right)^{-1} \frac{\partial Q}{\partial \eta_i}, \qquad i = 0, 1, 2. \end{align}\]

The columns of the Jacobian matrix \(DF\) can be computed from these derivatives together with the linearized flow \(d\varphi^\tau\), which we obtain by integrating the first-order variational equations along the trajectory. For example, the derivative with respect to \(\xi_1\) is given by \[\frac{\partial F}{\partial \xi_1} = dp \circ d\varphi^{\tau} \circ \begin{bmatrix} \frac{\partial \xi_0}{\partial \xi_1} & 1 & 0 & 0 & 0 & \frac{\partial \eta_0}{\partial \xi_1} & 0 & \frac{\partial \eta_3}{\partial \xi_1} \end{bmatrix}^T - \begin{bmatrix} 1 & 0 & 0 & 0 & 0 \end{bmatrix}^T,\] and similarly for the variables \(\xi_2\), \(\eta_1\), \(\eta_2\). The derivative with respect to the time variable \(\tau\) is \[\frac{\partial F}{\partial \tau} = dp \circ X_Q(\varphi^{\tau}(\xi_0, \xi_1, \xi_2, 0, \eta_0, \eta_1, \eta_2, \eta_3)).\]

9 Bifurcation Analysis via Conley–Zehnder Indices↩︎

To a non-degenerate periodic Hamiltonian orbit \(\gamma\) one can associate an integer valued invariant \(\mu_{CZ}(\gamma)\) called the (transverse) Conley–Zehnder index, [52]–[54]. Intuitively, this invariant measures the winding number of the linearized flow along \(\gamma\). We recall some important properties most relevant for our purposes:

1. (invariance along non-degenerate families) Let \(\gamma_s\) be a continuous family of non-degenerate periodic orbits of Hamiltonians \(H_s\). Then the Conley–Zehnder index is constant along the family: \[\mu_{CZ}(\gamma_s) = \mu_{CZ}(\gamma_{s_0}) \quad \text{for all } s.\]

2. (index jump implies bifurcation) Let \(\gamma_s\) be a continuous family of isolated periodic orbits of Hamiltonians \(H_s\). If \(\gamma_{s_0 \pm \varepsilon}\) are non-degenerate and satisfy \[\mu_{CZ}(\gamma_{s_0 - \varepsilon}) \neq \mu_{CZ}(\gamma_{s_0 + \varepsilon}),\] then there is a bifurcation occurring in the interval \((s_0-\varepsilon,s_0+\varepsilon)\) generating at least one additional orbit. This follows from the behavior of local Floer homology; see Section 3.2 of [55].

3. (parity) If \(\gamma\) is a non-degenerate periodic Hamiltonian orbit in a \(2n+2\)-dimensional phase space with reduced monodromy \(M_\gamma\), i.e., the differential of the Poincaré map restricted to the energy surface, then \[\text{sign}(\det(\mathrm{id}-M_\gamma))=(-1)^{\mu_{CZ}(\gamma)-n}.\]

A direct corollary of these properties is that if the type of precisely one pair of eigenvalues changes from positive hyperbolic to elliptic (or vice versa), then another periodic orbit family must be involved. These properties make the Conley–Zehnder index an effective tool for classification of periodic orbit families and bifurcation analysis. Recent studies have worked to develop practical and rigorous methods for computing the Conley–Zehnder index, [16], [56], as well as applying it to the analysis of bifurcation networks, [7], [17].

Here we recall a precise definition of the Conley–Zehnder index following the paper of Hofer, Wysocki, Zehnder, [57], see also [53] for an older reference. A detailed exposition of index theory can be found in the book of Long, [58]. First, we consider a path \(\psi: [0,1] \to Sp(2n)\) of symplectic matrices beginning at the identity \(\psi(0) = \mathrm{id}\), which is non-degenerate, meaning \(\psi(1)\) has no eigenvalue equal to \(1\). Extend \(\psi\) to a path \(\tilde{\psi}: [0,2]\to Sp(2n)\) ending at either \(W_+=-\text{id}\) or \(W_- = \text{diag}(2, 1/2, -1, \cdots, -1)\) in such a way that the extension avoids intersecting with the Maslov cycle \[V = \{\psi\in Sp(2n) \,|\, \det (\psi - \text{id}) = 0\}.\]

In other words, the extended path \(\tilde{\psi}\) retracts to a loop of unitary matrices based at the identity after taking the square, and the Conley–Zehnder index measures the winding number of the loop in the circle obtained by taking the complex determinant.

For a periodic Hamiltonian orbit \(\gamma\), the Conley–Zehnder index \(\mu_{CZ}(\gamma)\) is defined as the index of the symplectic path representing the transverse linearized flow along \(\gamma\). This requires choosing a suitable symplectic frame along the orbit. In this work, we follow the computation methods, including the choice of frames and extension procedures, described in [16]. For vertical collision orbits, we take as transverse frame \(\{\partial_{\xi_1},\partial_{\xi_2},\partial_{\eta_1},\partial_{\eta_2}\}\), following [59]. By Property [itm:cz-invariance], the index is constant along any non-degenerate family of periodic orbits, so it suffices to compute the index at a single orbit in the family. As an additional consistency check, we track index changes along branches via Floquet multipliers corresponding to eigenvalues of the first kind (see [53]), or equivalently, via Krein signatures of elliptic multipliers as suggested in Section 1.2.5 of [60] and in Section 3 of [7]. All bifurcation graphs are consistent with local Floer homology invariance.

The resulting Conley–Zehnder indices are shown as labels in the bifurcation graphs in Figures 7, 11, 13, 14, 15, 17, 19, 21, 24, 25 and summarized in Table 2. In exactly two instances, namely the indices 10 and 11 along the “left” tri-fly branch of the Saturn–Enceladus system in Figure 25, where the orbits lie very close to collision, we inferred the values of these indices from the bifurcation graph and the invariance of local Floer homology.

10 Tables of Numerical Data↩︎

This section lists representative numerical initial conditions and parameters for the periodic orbit families discussed throughout the paper. Tables 3, 4, 5, 6, 7 correspond, respectively, to families in Hill’s problem, \(L_1\) and \(L_2\) halo families, the W4 family in Earth–Moon system, the butterfly family in Saturn–Enceladus system, and tri-fly families in the Saturn–Enceladus and Earth–Moon systems. The velocity coordinates can be computed from the momentum coordinates using \[\dot{x} = p_x + y, \qquad \dot{y} = p_y - x, \qquad \dot{z} = p_z .\]

Acknowledgements↩︎

Chankyu Joung and Otto van Koert were supported by the National Research Foundation of Korea (NRF), grant number (MSIT) (RS-2023-NR076656), funded by the Korean government. Chankyu Joung was partially supported by BK21 SNU Mathematical Sciences Division. The authors thank Cengiz Aydin for many helpful discussions during the course of this project. Chankyu Joung and Dayung Koh are grateful to Damon Landau at the Jet Propulsion Laboratory for his kind support and for valuable conversations related to Enceladus orbits. Chankyu Joung thanks researchers at the Jet Propulsion Laboratory for valuable interactions during his visit, as well as Dongho Lee and Beomjun Sohn for their assistance. Finally, we are grateful to the referees for their helpful suggestions.

References↩︎

---

1. Another limit to consider occurs as the Hamiltonian energy decreases to \(-\infty\). After an appropriate rescaling, the dynamics of the energy surface converges to the Kepler problem. In this setting, Moser, [6], showed using the averaging method that periodic families bifurcate from a set of four “basic” Keplerian orbits: the retrograde and direct planar circular orbits, and a pair of vertical collision orbits.↩︎

2. The period \(T\) is included in the map \(F\), because the \(k\)-th arrival time at the Poincaré section may be discontinuous in the perturbation parameter.↩︎