## 1 Introduction  
In the last decade, the number of proposed formation flying missions has increased significantly [1, 2, 3], led by the growing number of low-cost rideshare opportunities to low Earth orbits (LEO) with medium inclination (50-55 degrees) and to sun-synchronous orbits (SSO). Satellite formations operating in these orbits are affected by perturbations associated with the non-uniform Earth gravitational field [4, 5], which depend on the position of the spacecraft and therefore have a different impact on each satellite of the formation [6, 7]. Satellite formations designed using the traditional approach based on the Hill-Clohessy-Wiltshire (HCW) equations show an unstable behavior over time, characterized by marked variations in the amplitude of the in-plane and out-of-plane oscillations, and by the drift of the relative orbit center [8, 9]. These effects have a dramatic impact on the performance of the formation, since the relative distance between the satellites can increase and exceed the operative range or decrease down to potential collision. In such a scenario, periodic reconfiguration maneuvers are required and they can eventually involve expensive coordination activity from the ground segment [10, 11, 12, 13, 14]. An effective alternative is to identify configurations which are stable with respect to the gravitational perturbations. Several research works address the identification of satellite formations that are invariant under gravitational perturbations, most of them considering only the effects associated with the dominant \(\mathbb{J}2\) term. The first work on this topic is by [PERSON] and [PERSON], who derived two first-order conditions to design \(\mathbb{J}2\)-invariant relative orbits, which can be applied on a large number of practical cases [15]. [PERSON] and [PERSON] developed a high-fidelity linearized model which takes into account the time-averaged effect of J2 [16], [PERSON], [PERSON], and [PERSON] used the Hamiltonian formalism and the variation of parameters procedure to determine configurations stable under the effects of J2 [17]. [PERSON], [PERSON], and [PERSON] identified some values of the inclination corresponding to minimum drift orbits for the J2 perturbed case [18]. He, [PERSON], and [PERSON] developed a model to design bounded relative orbits in the full zonal problem, therefore including the \(J_{n}\) terms up to the desired degree, based on the Taylor approximation of Poincare maps and numerical refinement [19]. [PERSON] et al. examined the same problem to determine invariant relative orbits (negligible drift and variation in the amplitude of the oscillations) using double-averaged Delaunay orbital elements [20]. In this work, a model providing a compact representation of invariant relative orbits for formation flying in the full zonal problem is introduced. Using the Hamiltonian formalism, canonical transformations allow absorbing the terms associated with the perturbation [21,22] and providing a new representation of the system dynamics in terms of three mutually orthogonal couples of state variables. It is proved that the relative orbit can be characterized in terms of four parameters, corresponding to the value of four state variables at a given time, with two of them defining the boundedness condition (or equivalently the drift rate) and the other two characterizing (independently) the amplitude of the in-plane and out-of-plane oscillations. The model proposed here can be conveniently used to design bounded satellite formations stable under the full zonal harmonic perturbation. Furthermore, it allows developing a simple guidance strategy to inject a spacecraft approaching a target satellite into a stable bounded relative orbit. The effectiveness of the model and of the guidance strategy proposed are verified by means of numerical integration, using gravity field models including harmonic coefficients with high degree and order and real ephemerides.