## 3 Characterization of Invariant Relative Orbits
### Canonical Transformations Absorbing the Zonal Harmonic Perturbation  
The canonical transformations \(\mathcal{T}_{i}:(\mathbf{q},\mathbf{p})\rightarrow(\mathbf{Q},\mathbf{P})\) are developed to rearrange the Hamiltonian function from the form given by Equation (5) to the following one:  
\[\widehat{H}=H^{(0)}(\mathbf{Q},\mathbf{P})+o(J_{n}^{k}) \tag{9}\]  
where \(o(J_{n}^{2})\) represents the negligible higher-order terms associated with the zonal harmonic perturbation.  
As shown by Equations (6)-(8), the dynamics onto the \([\hat{x},\hat{y}]\) plane are not coupled to the one along the \(\hat{z}\) direction, so the canonical transformations for the in-plane \((x,y,p_{x},p_{y})\) and the out-of-plane \((z,p_{z})\) state variables can be developed separately. The perturbation terms in the out-of-plane dynamics can be absorbed by introducing a type-3 generating function \(S_{3}(z,P_{z})\) such that  
\[\begin{cases}p_{z}=\frac{\partial S_{3}}{\partial z}\\ Q_{z}=\frac{\partial S_{3}}{\partial P_{z}}\end{cases} \tag{10}\]  
Setting \(S_{z}=\frac{\partial S_{3}}{\partial z}\), introducing Equation (11) into Equations (6)-(8), and isolating the out-of-plane terms results in  
\[\bar{H}_{z}=\frac{S_{z}^{2}}{2}+\frac{1+K_{33}^{2}}{2}z^{2}-v_{z}^{(1)}S_{z}+ K_{30}z \tag{11}\]  
The desired transformation \(S_{3}\) can be found after determining the expression of \(S_{z}\) such that  
\[\bar{H}_{z}=\frac{P_{z}^{2}}{2}+\frac{1+K_{33}^{2}}{2}z^{2} \tag{12}\]  
Equations (11)-(12) can be solved by quadrature  
\[S_{z}^{2}-2v_{z}^{(1)}S_{z}+2K_{30}z-P_{z}^{2}=0\rightarrow\]  
\[\to S_{z}=v_{z}^{(1)}\pm\sqrt{{v_{z}^{(1)}}^{2}-(2K_{30}z-P_{z}^{2})} \tag{13}\]  
and Equation (13) can be integrated to produce  
\[S_{3}(z)=v_{z}^{(1)}z\mp\frac{\left[{v_{z}^{(1)}}^{2}-2K_{30}z+P_{z}^{2} \right]^{\frac{3}{2}}}{3K_{30}}+C_{3} \tag{14}\]  
where \(C_{3}\) is the (inessential) constant of integration. Introducing Equation (13) into Equation (11) results in the following canonical changes in coordinates:  
\[\mathcal{T}_{3}:z=\frac{\left[{p_{z}^{2}+v_{z}^{(1)}}^{2}-K_{30}^{2}\frac{Q_ {z}^{2}}{P_{z}^{2}}\right]}{2K_{30}} \tag{15}\]  
\[\mathcal{T}_{6}:p_{z}=v_{z}^{(1)}\pm K_{30}\frac{Q_{z}}{P_{z}} \tag{16}\]  
The in-plane terms are too complex to be processed and the use of normal forms is preferred. As proved in Appendix B, using the type-3 generating function \(S(x,y,P_{x},P_{y})\) results in the following canonical change in coordinates:  
\[\mathcal{T}_{1}:x=Q_{x}-g_{1} \tag{17}\]  
\[\mathcal{T}_{2}:y=Q_{y}-g_{2} \tag{18}\]\[\mathcal{T}_{3}:p_{x}=P_{x}+f_{1} \tag{19}\]  
\[\mathcal{T}_{4}:p_{y}=P_{y}+f_{2} \tag{20}\]  
where  
\[\begin{cases}f_{1}=\frac{v_{x}^{(1)}(1+2K_{22})+K_{20}}{2K_{22}}\\ f_{2}=\frac{K_{10}}{2(1-K_{11})}+\frac{v_{y}^{(1)}(3-2k_{11})}{2(1-K_{11})}\\ g_{1}=v_{y}^{(1)}-f_{2}\\ g_{2}=f_{1}-v_{x}^{(1)}\end{cases} \tag{21}\]  
The Hamiltonian function obtained after applying the \(\mathcal{T}_{i}\) transformations is reported below  
\[\tilde{H}=\frac{1}{2}\big{(}P_{x}+Q_{y}\big{)}^{2}+\frac{1}{2}\big{(}P_{y}-1- Q_{x}\big{)}^{2}+\bigg{(}K_{11}-\frac{3}{2}\bigg{)}Q_{x}^{2}+K_{22}Q_{y}^{2}+ \frac{P_{z}^{2}}{2}+\frac{1+K_{33}^{2}}{2}Q_{z}^{2} \tag{22}\]  
The reader can verify that the Hamilton Equation (3) obtained using \(\tilde{H}\) correspond to the following set of second-order ordinary differential equations:  
\[\begin{cases}\dot{Q}_{x}=(2K_{11}-3)Q_{x}+2\dot{Q}_{y}\\ \dot{Q}_{y}=-2\dot{Q}_{x}-2K_{22}Q_{y}\\ \dot{Q}_{z}=-\big{(}1+K_{33}^{2}\big{)}Q_{z}\end{cases} \tag{23}\]  
As declared at the beginning of this section, Equation (23) is equivalent to the HCW equations except for the negligible terms \(K_{ii}\).