## 2 Hamiltonian Formalism for Formation Flying under Zonal Harmonic Perturbation  
The relative motion between two spacecraft \(SC_{1}\) and \(SC_{2}\) can be conveniently described using an orthogonal reference system \(\mathcal{F}_{r}=\left[ \begin{array}{cc}\hat{x}&\hat{y}&\hat{z}\end{array}  \right]^{T}\) centered in \(SC_{1}\), named the chief satellite, with axes \(\hat{x}\) pointing from the center of the Earth to \(SC_{1}\), \(\hat{z}\) parallel to the angular momentum of the chief orbit, and \(\hat{y}\) completing the orthogonal frame. The relative position and velocity of \(SC_{2}\), named the deputy satellite, with respect to \(SC_{1}\) are given by \(\rho=\left[ \begin{array}{cc}x&y&z\end{array} \right]^{T}\) and \(\dot{\rho}=\frac{d\rho}{dt}=\left[ \begin{array}{cc}v_{x}&v_{y}&v_{z}\end{array}  \right]^{T}\). The dynamics of the two satellites evolve under the effects of the Earth gravitational potential, which can be modeled using spherical harmonics [4]:  
\[V_{i}=-\frac{\mu}{r_{i}}\left[1-\sum_{n=2}^{\infty}\sum_{m=0}^{n}\left(\frac{ R_{E}}{r_{i}}\right)^{n}P_{nm}(\cos\phi)(C_{nm}\cos m\lambda+S_{nm}\sin m \lambda)\right] \tag{1}\]  
where \(\mu=3.986\times 10^{5}\) km\({}^{3}/s^{2}\) and \(R_{E}=6378\) km are, respectively, the Earth gravitational parameter and equatorial radius, \(r_{i}\) is the orbit radius of the \(i\)-th spacecraft, \(P_{nm}\) are the Legendre associated functions of degree \(n\) and order \(m\), \(C_{nm}\) and \(S_{nm}\) are the spherical harmonics coefficients, \(\phi\) is the geocentric co-latitude, and \(\lambda\) is the longitude. To develop the model proposed in Section 3, Equation (1) is simplified considering only the zonal harmonic terms (\(m=0\)). \[V_{i}=-\frac{\mu}{r_{i}}\Bigg{[}1-\sum_{n=2}^{\infty}J_{n}\bigg{(}\frac{R_{E}}{r_ {i}}\bigg{)}^{n}P_{n}(\cos\phi)\Bigg{]} \tag{2}\]  
It shall be noted that for satellite formations operating in Earth orbits, this approximation is acceptable, based on the fact that the secular perturbation associated with the zonal harmonic terms is larger by more than one order of magnitude than the periodic perturbations associated with the tesseral and sectorial harmonic ones (\(m>0\)) of the same degree [30]. Nevertheless, this is not a general case and, therefore, the model developed under this approximation is verified by means of numerical integration, discussed in Section 4, using orbital propagators which include all the harmonic coefficients. The relative dynamics are expressed here using the Hamiltonian formalism in terms of the dimensionless state variables \(\mathbf{q}\) and \(\mathbf{p}\):  
\[ \begin{cases}\dot{\mathbf{q}}=\frac{\partial H}{\partial\mathbf{p}}\\ \dot{\mathbf{p}}=-\frac{\partial H}{\partial\mathbf{q}}\end{cases}  \tag{3}\]  
For the sake of compactness, the derivation of the model is presented in Appendix A and the expressions of the conjugate momenta and of the Hamiltonian function \(H\) are reported below:  
\[ \begin{cases}p_{x}=v_{x}-y+v_{x}^{(1)}\\ p_{y}=v_{y}+x+1+v_{y}^{(1)}\\ p_{z}=v_{z}+v_{z}^{(1)}\end{cases}  \tag{4}\]  
\[H=H^{(0)}(\mathbf{q},\mathbf{p})+H^{(1)}(\mathbf{q},\mathbf{p},J_{n})+H^{(2)} \Big{(}\mathbf{q},\mathbf{p},J_{n}^{2}\Big{)} \tag{5}\]  
where  
\[H^{(0)}=\frac{1}{2}(p_{x}+y)^{2}+\frac{1}{2}(p_{y}-1-x)^{2}+\frac{p_{z}^{2}}{ 2}-\frac{3}{2}x^{2}+\frac{z^{2}}{2} \tag{6}\]  
\[H^{(1)}=-v_{x}^{(1)}p_{x}-v_{y}^{(1)}p_{y}-v_{z}^{(1)}p_{z}+K_{10}x+K_{20}y+K_ {30}z \tag{7}\]  
\[H^{(2)}=K_{11}x^{2}+K_{22}y^{2}+K_{33}z^{2} \tag{8}\]  
The reader will notice that the expression of \(H^{(0)}\) is equal to that of the Hamiltonian function for the unperturbed case, while \(H^{(1)}\) and \(H^{(2)}\) depend on, respectively, the parameters \(J_{n}\) and \(J_{n}^{2}\), collected in the coefficients \(v_{i}^{(1)}\) and \(K_{ij}\). As detailed in Appendix A.2, the parameter \(J_{n}\) is small (i.e., \(J_{n}<<1\)).