## Appendix A Derivation of the Hamiltonian Formalism
### Lagrangian for the Relative Motion under the Zonal Harmonic Perturbation  
The reader shall recall from Section 2 that the relative coordinates of the deputy satellite with respect to the chief are expressed in the reference frame \(\mathcal{F}_{r}\), using the symbol \(\boldsymbol{\rho}\) to indicate the relative position. \(\mathcal{F}_{r}\) rotates with respect to an Earth-centered inertial (ECI) frame with an angular velocity \(\boldsymbol{\omega}\) depending on both the orbit angular speed \(n\) of the chief satellite and that of its orbit plane caused by the gravitational perturbations. The expression of \(\mathcal{L}_{r}\) is here developed considering only the secular effect of \(J_{2}\) on the rotation of the orbit plane, which causes the time drift in the RAAN and in the argument of latitude expressed below:  
\[\dot{\Omega}=-\frac{3}{2}nJ_{2}\bigg{(}\frac{R_{E}}{\bar{r}}\bigg{)}^{2}\cos i \tag{1}\]  
\[\dot{u}=n+\delta n \tag{2}\]  
with  
\[n=\sqrt{\frac{\mu}{\bar{a}^{3}}} \tag{3}\]  
\[\delta n=\frac{3}{4}nJ_{2}\bigg{(}\frac{R_{E}}{\bar{r}}\bigg{)}^{2}\bigg{(}3 -\frac{7}{2}\sin^{2}i\bigg{)} \tag{4}\]  
\[\bar{r}=\bar{a}+\frac{3}{4}J_{2}\bigg{(}\frac{R_{E}}{\bar{a}}\bigg{)}^{2} \bigg{(}3\sin^{2}i-2\bigg{)} \tag{5}\]  
where \(\bar{r}\) is the mean radius of the \(J_{2}\)-perturbed orbit, \(\bar{a}\) is the semi-major axis, and \(i\) is the inclination of the orbital plane. Based on Equations (1)-(5), the expression of \(\boldsymbol{\omega}\) is given by  
\[\boldsymbol{\omega}=\left[ \begin{array}{c}\dot{\Omega}\sin i\sin u\\ \dot{\Omega}\sin i\cos u\\ \dot{\Omega}\cos i+\dot{u}\end{array} \right] \tag{6}\]  
Equation (6) is now used to derive the expression for the relative velocity in the inertial frame, reported below, and that of the relative kinetic energy \(K_{r}\). \[\mathbf{v}=\boldsymbol{\omega}\times\bar{\mathbf{r}}+\frac{d\boldsymbol{\rho }}{dt}+\boldsymbol{\omega}\times\boldsymbol{\rho} \tag{7}\]For the sake of compactness, it is worth isolating from Equation (A7) the terms \(\mathbf{v}^{(1)}\), used in Equation (4), depending on the perturbation (i.e., on \(J_{2}\)). \[\mathbf{v}=\mathbf{v}^{(0)}+\mathbf{v}^{(1)}\] (A8)  
where  
\[\mathbf{v}^{(0)}=\left[ \begin{array}{c}\dot{x}-ny\\ \dot{y}+nx+n\bar{r}\\ \dot{z}\end{array} \right]\] (A9)  
and  
\[\mathbf{v}^{(1)}=\left[ \begin{array}{c}\dot{\Omega}(z\sin i\cos u-y\cos i)- y\delta n\\ \left(\dot{\Omega}\cos i+\delta n\right)(x+\bar{r})-z\Omega\sin i\sin u\\ \dot{\Omega}\sin i(y\sin u-(x+\bar{r})\cos u)\end{array} \right]\] (A10)  
Finally, the relative kinetic energy can be expressed in compact form based on Equations (A8)-(A10):  
\[K_{r}=\frac{1}{2}|\mathbf{v}^{(0)}|^{2}+\mathbf{v}^{(0)^{T}}\mathbf{v}^{(1)}+ \frac{1}{2}|\mathbf{v}^{(1)}|^{2}\] (A11)  
Once the expression of \(K_{r}\) is derived, an equivalent representation of the relative potential \(V_{r}\) is necessary to obtain the Lagrangian. As suggested by [PERSON] et al. [17], the expression of \(V_{r}\), including the zonal harmonic coefficient, can be obtained after expressing Equation (2) for \(SC_{2}\) and \(SC_{1}\) in \(\mathcal{F}_{r}\) and then subtracting the two, resulting in the following form:  
\[V_{r}=-\frac{\mu}{\bar{r}+\rho}\left[1-\sum_{n=2}^{\infty}J_{n}\left(\frac{R_{ E}}{\bar{r}+\rho}\right)^{n}P_{n}(\cos\phi)\right]\] (A12)  
The geocentric co-latitude \(\phi\) can be obtained from the projection of \(\mathbf{r}=\bar{r}_{1}+\rho\) along the \(\xi_{3}\) axis of the ECI frame, expressed through the 3-1-3 rotation \(\mathcal{R}_{ri}=\mathcal{R}_{3}^{T}(\Omega)\mathcal{R}_{1}^{T}(i)\mathcal{R }_{3}^{T}(u)\):  
\[\cos\phi=\frac{(\bar{\mathbf{r}}+\rho)\mathcal{R}_{ri}^{T}\xi_{3}}{\bar{r}+ \rho}=\left[ \begin{array}{cc}\sin i\sin u&\sin i\cos u&\cos i\end{array}