Characterizing Accuracy of Normal Forms to Study Trajectories in Cislunar Space

Abstract

Qualitative understanding of cislunar trajectories is increasingly important as lunar missions become more commonplace. Researchers have used the Lie series method to reduce the Circular Restricted Three-Body Problem (CR3BP) to its normal form up to some approximation order and in the vicinity of the five libration points. This approximation allows for analytical propagation in proximity to the libration points by defining action-angle variables. These variables are such that, up to the approximation order, the actions are constant and the angles are linear in time. An objective of this work is to examine how the normal form coordinates characterize trajectories in the vicinity of the libration points and maintain accuracy of propagation. Normal form coordinates qualitatively separate periodic, quasiperiodic, transit, and reflective trajectories. Another objective of this work is to examine the accuracy of the approximate normal form centered at $L_1$ and $L_2$ at various approximation orders, distances, and energy levels. At higher approximation orders, the normal form is able to accurately propagate trajectories in a ball around the libration point of origin. Finally, two example applications of this method are then examined, including maneuver characterization and Halo orbit identification.

Appendix A Preliminary Transformations

1.1 A.1 Polynomial Approximation and Resulting Hamiltonian

New equations of motion for the system defined from (6) may be approximated by taking series expansions of the two gravitational potentials:

$$\begin{aligned} \frac{1}{r_i(x,y,z)} = \frac{1}{\Vert \varvec{\rho } - \textbf{d}_i\Vert } = \sum _{n\ge 0}^{\infty } \frac{\rho ^n}{d_i^{n+1}} P_{n} \left( \frac{\varvec{\rho } \cdot \textbf{d}_i}{\rho d_i} \right) \end{aligned}$$

(A1)

under the assumption $\rho < d_i$, where $\rho ^2 = x^2 + y^2 +z^2$, $\textbf{d}_1 = [\frac{-1\pm \gamma }{\gamma }, 0, 0]^T$ for the Earth’s potential, $\textbf{d}_2 = [\pm 1, 0, 0]^T$ for the Moon’s potential, and $P_{n} \left( \cdot \right)$ are the Legendre polynomials of the first kind. As shown by Richardson [17], these new equations of motion are

$$\begin{aligned} \begin{aligned} \ddot{x} - 2\dot{y} - (1+2c_2)x&= \frac{\partial }{\partial x} \sum _{n\ge 3} c_n \rho ^n P_{n} \left( \frac{x}{\rho } \right) \\ \ddot{y} + 2 \dot{x} + (c_2 -1) y&= \frac{\partial }{\partial y} \sum _{n\ge 3} c_n \rho ^n P_{n} \left( \frac{x}{\rho } \right) \\ \ddot{z} +c_2 z&= \frac{\partial }{\partial z} \sum _{n\ge 3} c_n \rho ^n P_{n} \left( \frac{x}{\rho } \right) , \\ \end{aligned} \end{aligned}$$

(A2)

where $c_n$ captures the coefficients of the expansion

$$\begin{aligned} c_n(\mu ) = \frac{1}{\gamma _j^3}\left( (\pm 1)^n \mu + (-1)^n (1-\mu ) \left( \frac{\gamma _j}{1 \mp \gamma _j}\right) ^{n+1} \right) , \end{aligned}$$

where again, the upper sign is for $L_1$, the lower sign is for $L_2$, $\gamma _j$ is the distance from the libration point of interest to the Moon, and $j=1,2$ indexes the libration point of interest. Transforming to the conjugate momenta–defined in the same way as in the synodic reference frame, but now in the libration-centered coordinates–the new Hamiltonian–first shown in (7)–is restated here:

$$\begin{aligned} \begin{aligned} H(x, y, z, p_x, p_y, p_z)&= \frac{1}{2} (p_x^2 + p_y^2 +p_z^2) + yp_x-xp_y\\&\hspace{3em}-c_2 x^2 + \frac{c_2}{2} (y^2 +z^2) - \sum _{n\ge 3} c_n \rho ^n P_{n} \left( \frac{x}{\rho } \right) . \end{aligned} \end{aligned}$$

(A3)

1.2 A.2 Diagonalization

Letting the quadratic terms be collected in $H_2$, the objective now is to define a symplectic change of variables to eliminate the cross terms from $H_2$ and decouple the linearized system. The z-direction is already a decoupled harmonic oscillator with a frequency of $\nu _0=\sqrt{c_2}$ in the linearized system [37], therefore the planar case is examined. Under the planar assumption, the EOM may be written as

$$\begin{aligned} \dot{\varvec{\xi }} = {{\textbf {J}}} \nabla H_2(\varvec{\xi }) = {{\textbf {M}}}\varvec{\xi } \end{aligned}$$

(A4)

where $\varvec{\xi }=[x,y,p_x,p_y]^T$ and

$$\begin{aligned} {{\textbf {J}}} = \begin{bmatrix} \textbf{0} &{} \textbf{I}_2 \\ -\textbf{I}_2 &{}\textbf{0} \end{bmatrix} \qquad {{\textbf {M}}}=\begin{bmatrix} 0&{}1&{}1&{}0 \\ -1&{}0&{}0&{}1 \\ 2c_2&{}0&{}0&{}1 \\ 0&{}-c_2&{}-1&{}0\end{bmatrix}. \end{aligned}$$

The characteristic polynomial of ${{\textbf {M}}}$ is $\lambda ^4 +(2-c_2)\lambda ^2 + (1+c_2-2c_2^2)=0$, giving

$$\begin{aligned} \eta _{1,2} = \frac{c_2-2 \mp \sqrt{9c_2^2-8c_2}}{2} \end{aligned}$$

(A5)

where $\eta =\lambda ^2$. Given $\mu \in [0,\frac{1}{2})$ forces $c_2>1$ thus causing $\eta _1 < 0$ and $\eta _2>0$, the collinear libration points are type saddle$\times$center$\times$center when the z-component is also included [37]. The variables $\omega _0$ and $\lambda _0$ are defined such $\omega _0 = \pm \sqrt{-\eta _1}$ and $\lambda _0 = \pm \sqrt{\eta _2}$, where, the correct sign will be discussed later. Let the eigenvectors of the matrix ${{\textbf {M}}}$, which are defined and carried out in reference [21, 37], initially define the transformation. This resulting change of variables may be cast as

$$\begin{aligned} \textbf{C} = \begin{bmatrix} {2\lambda _0} &{} 0 &{} {-2\lambda _0} &{} {2\omega _0} \\ {\lambda _0^2-2c_2-1} &{} {-\omega ^2-2c_2-1} &{} {\lambda _0^2-2c_2-1} &{} 0 \\ {\lambda _0^2+2c_2+1}&{} {-\omega ^2+2c_2+1} &{} {\lambda _0^2+2c_2+1} &{} 0 \\ {\lambda _0^3 +(1-2c_2)\lambda _0} &{} 0 &{} {-\lambda _0^3 - (1-2c_2)\lambda _0} &{} {-\omega _0^3 + (1-2c_2)\omega _0} \end{bmatrix}. \end{aligned}$$

A symplectic transform must satisfy $\textbf{C}^T{{\textbf {J}}}\textbf{C}={{\textbf {J}}}$, but in this case

$$\begin{aligned} \textbf{C}^T{{\textbf {J}}}\textbf{C} = \begin{bmatrix} 0&{}0&{} d_{\lambda _0}&{} 0 \\ 0&{}0&{}0&{} d_{\omega _0} \\ d_{\lambda _0}&{}0 &{}0&{}0 \\ 0&{}d_{\omega _0}&{}0&{}0\end{bmatrix}. \end{aligned}$$

Letting $s_{\lambda _0}=\sqrt{d_{\lambda _0}}$ and $s_{\omega _0}=\sqrt{d_{\omega _0}}$, the first and third columns of $\textbf{C}$ must be scaled by

$$\begin{aligned} s_{\lambda _0}= \sqrt{2\lambda _0 ((4+3c_2)\lambda _0^2 +4 + 5c_2 - 6c_2^2)},\\ \end{aligned}$$

while the second and fourth columns of $\textbf{C}$ must be scaled by

$$\begin{aligned} s_{\omega _0}= \sqrt{\omega _0((4+3c_2)\omega _0^2 - 4 - 5c_2 +6c_2^2)}. \end{aligned}$$

In order for these scaled quantities to be real, then $\lambda _0, \omega _0 >0$ as stated in reference [21, 37].

Given these scaling factors and reintroducing the z-direction under the transformation

$$\begin{aligned} z \rightarrow \frac{1}{\sqrt{\nu _0}} z, \qquad p_z \rightarrow \sqrt{\nu _0} p_z, \end{aligned}$$

(A6)

the symplectic change of variables from the libration-centered coordinate, $\varvec{\xi } = [x, y, z, p_x, p_y, p_z]^T$, to the diagonal coordinates, $\tilde{\varvec{\xi }} = [\tilde{x}, \tilde{y}, \tilde{z}, \tilde{p}_x, \tilde{p}_y, \tilde{p}_z]^T$, is given by

$$\begin{gathered} \mathcal {T}_{diag}:(x, y, z, p_x, p_y, p_z) \longrightarrow (\tilde{x}, \tilde{y}, \tilde{z}, \tilde{p}_x, \tilde{p}_y, \tilde{p}_z),\nonumber \\ \varvec{\xi } = \textbf{C} \tilde{\varvec{\xi }} \end{gathered}$$

(A7)

where the scaled matrix $\textbf{C}$ with the z-transform included is

$$\begin{aligned} \textbf{C} = \begin{bmatrix} \frac{2\lambda _0}{s_{\lambda _0}} &{} 0 &{} 0 &{} \frac{-2\lambda _0}{s_{\lambda _0}} &{} \frac{2\omega _0}{s_{\omega _0}} &{} 0 \\ \frac{\lambda _0^2-2c_2-1}{s_{\lambda _0}} &{} \frac{-\omega _0^2-2c_2-1}{s_{\omega _0}} &{} 0 &{} \frac{\lambda _0^2-2c_2-1}{s_{\lambda _0}} &{} 0 &{} 0\\ 0 &{} 0 &{} \frac{1}{\sqrt{\nu _0}} &{}0 &{}0 &{}0 \\ \frac{\lambda _0^2+2c_2+1}{s_{\lambda _0}}&{} \frac{-\omega _0^2+2c_2+1}{s_{\omega _0}} &{} 0 &{} \frac{\lambda _0^2+2c_2+1}{s_{\lambda _0}} &{} 0 &{} 0\\ \frac{\lambda _0^3 +(1-2c_2)\lambda _0}{s_{\lambda _0}} &{} 0 &{} 0 &{} \frac{-\lambda _0^3 - (1-2c_2)\lambda _0}{s_{\lambda _0}} &{} \frac{-\omega _0^3 + (1-2c_2)\omega _0}{s_{\omega _0}} &{} 0 \\ 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} \sqrt{\nu _0} \end{bmatrix}. \end{aligned}$$

Applying this transform to the Hamiltonian in (7)

$$\begin{aligned} H(\varvec{\xi }) \circ \mathcal {T}_{diag} = \lambda _0 \tilde{x}\tilde{p}_x+ \frac{\omega _0}{2}(\tilde{y}^2 + \tilde{p}_y^2) + \frac{\nu _0}{2}(\tilde{z}^2 + \tilde{p}_z^2) + \sum _{n\ge 3} H_n(\tilde{x}, \tilde{y}, \tilde{z}, \tilde{p}_x, \tilde{p}_y, \tilde{p}_z), \end{aligned}$$

(A8)

where $H_n(\tilde{x}, \tilde{y}, \tilde{z}, \tilde{p}_x, \tilde{p}_y, \tilde{p}_z)$ collects higher order terms.

1.3 A.3 Complexification

Additionally, using complex variables, $\textbf{q}$ and ${{\textbf {p}}}$, when defining the Lie series method’s generating functions is advantageous. Therefore, the complexifying transformation is introduced such that

$$\begin{gathered} \mathcal {T}_{comp}: (\tilde{x}, \tilde{y}, \tilde{z}, \tilde{p}_x, \tilde{p}_y, \tilde{p}_z) \longrightarrow ({\textbf{q}}, {{{\textbf {p}}}})\\ \begin{bmatrix} \tilde{x}\\ \tilde{p}_x\end{bmatrix} = \mathbf {I_{2\times 2}} \begin{bmatrix} q_1\\ {p}_1\end{bmatrix}, \qquad \begin{bmatrix} \tilde{y}\\ \tilde{p}_y\end{bmatrix} = \textbf{Q} \begin{bmatrix} q_2\\ {p}_2\end{bmatrix}, \qquad \begin{bmatrix} \tilde{z}\\ \tilde{p}_z\end{bmatrix} = \textbf{Q} \begin{bmatrix} q_3\\ {p}_3\end{bmatrix}, \end{gathered}$$

(A9)

where $\mathbf {I_{2\times 2}}$ is the identity maxtrix and

$$\begin{aligned} \textbf{Q} = \frac{1}{\sqrt{2}} \begin{bmatrix} 1 &{} \sqrt{-1} \\ \sqrt{-1} &{} 1 \end{bmatrix}. \end{aligned}$$

1.4 A.4 Action-Angle Transformation

Finally, the action-angle transformtion is defined:

$$\begin{gathered}\mathcal {T}_{AA} : (\bar{\textbf{q}}, \bar{{{\textbf {p}}}}) \longrightarrow (I_1, I_2, I_3, \phi _1, \phi _2, \phi _3)\nonumber \\ \begin{aligned} &\quad \bar{q}_1= \sqrt{I_1} \exp (\phi _1) \qquad \qquad&\bar{p}_1&= \sqrt{I_1} \exp ({-\phi _1}) \nonumber \\&\quad \bar{q}_2= {\sqrt{I_2}} \exp ({\sqrt{-1}\phi _2})&\bar{p}_2&= -\sqrt{-1}{\sqrt{I_2}} \exp ({-\sqrt{-1}\phi _2}) \nonumber \\&\quad \bar{q}_3= {\sqrt{I_3}} \exp ({\sqrt{-1}\phi _3})&\bar{p}_3&= -\sqrt{-1}{\sqrt{I_3}} \exp ({-\sqrt{-1}\phi _3}). \end{aligned} \end{gathered}$$

(A10)

The Hamiltonian of (7) expressed in action-angle coordinates is given by

$$\begin{aligned}{} & {} H(\varvec{\phi }, \textbf{I}) = \left( \left( H(x, y, z, p_x, p_y, p_z) \circ \mathcal {T}_{diag} \right) \circ \mathcal {T}_{comp} \right) \circ \mathcal {T}_{AA} \nonumber \\{} & {} \quad H(\varvec{\phi }, \textbf{I}) = \lambda _0 I_1+ \omega _0 I_2+ \nu _0 I_3+ \sum _{n\ge 3} H_n(\varvec{\phi }, \textbf{I}). \end{aligned}$$

(A11)