High-order analytical solutions of bounded relative motions for Coulomb formation flying

Close-proximity Coulomb formation flying offers attractive prospects in multiple astronautical missions. A semi-analytical method of constructing the series solution up to an arbitrary order for the relative motion near the equilibrium of Coulomb formation systems is proposed to facilitate the design of Coulomb formations based on a Lindstedt–Poincaré method. The details of the series expansion and coefficient solution for arbitrary m:n\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$m\!:\!n$$\end{document}-period orbits are discussed. To verify the effectiveness of the proposed method, the practical convergence domain of the high-order series solution is computed by comparing it with corresponding numerical solutions that satisfy the specified boundary conditions. Simulation results demonstrate the efficacy of the Lindstedt–Poincaré method in constructing series solutions for the relative motion of Coulomb formation systems.


Introduction
Close-proximity spacecraft formation has huge potential in telescope-occulter pairs, space-borne interferometry, and Earth imaging [1]. The fast throttle ability and avoidance of plume impingement issues of Coulomb control make it get attracted sustained attention in close-proximity formation mission design [2]. By actively charging spacecraft, repulsive and attractive forces can be generated to execute Coulomb control over a close formation. Several space missions, including SCATHA [3], ATS [4], and CLUSTER [5], have preliminary explored the feasibility of active spacecraft potential control [6].
Unlike traditional formation flying, Coulomb force control can simply generate more types of static configurations, in which spacecraft remain stationary with respect to the rotating center-of-mass frame [7]. On the contrary, the dynamical configurations of Coulomb formation are more complicated due to its high nonlinearity and strong-coupling feature. Therefore, early researches focused more on the dynamics and control of static Coulomb formations. Natarajan et al. [8] investigated a particular static formation shape of a twocraft Coulomb formation, where the spacecraft aligned the orbit nadir direction. Then, a charge feedback control law is firstly introduced to stabilize the formation. Hogan et al. [9] discovered families of multiple invariant shape solution for a collinear three-craft Coulomb formation in the absence of relevant gravitational forces. More generally, an analytical method for computing the required charge to maintain static N-craft Coulomb formations is proposed by Berryman and Schaub in Ref. [10]. Inampudi et al. [11] analyzed the orbit and attitude stability of a two-craft virtual Coulomb structure at Earth-Moon libration points. Then, the gravity gradient torques are utilized to stabilize the formation in orbit radial direction. Alikhani et al. [12] studied a static triangular three-craft Coulomb formation and presented a fault-tolerance control law to maintain the desired formation.
Another topic of particular interest is the design, maintenance, and reconfiguration of the dynamic Coulomb formations. Coulomb thrust leads to the high nonlinearity of the formation dynamics, which makes the dynamic formations have rich but very complicated dynamic properties. Jones [13] extended the static Coulomb formation research and derived the periodic relative orbit of two-craft Coulomb formation under the action of open-loop charge for the first time by using the dynamic system theory. In Ref [14], the author further studied three-craft collinear equilibria in the presence of a central body gravity field and discussed the possibility of exploiting invariant manifolds to achieve minimum-fuel reconfigurations. Wang [15] developed circular transfer trajectory and the general patched conic section trajectory programming algorithms for the maneuver mission of two-craft formation. Aslanov [16] studied the attitude dynamics of a defunct satellite with flexible appendages under the action of Coulomb forces and found the chaotic motion in Coulomb formation for the first time. Lin et al. [17] also revealed the intermittent chaos generated from a tangent bifurcation point in the Coulomb formation. Ref. [18] studied two different two-craft Coulomb formations in nearcircular reference orbits. The Udwadia-Kalaba formulation is implemented to achieve and maintain the for-mation along the nadir direction and the formation with a constant separation distance.
The traditional spacecraft formation in a circular or near-circular orbit can be designed by solving the Clohessy-Wiltshire (CW) equation and determining the parameters of its linear solution [19,20]. However, the high nonlinearity feature of the dynamical Coulomb formation introduces a variety of dynamical complexities, making a linear solution unable to describe the relative motion well. Thus, tools for analytical descriptions of the nonlinear relative motion are necessary to study the dynamics of the Coulomb formation. Many approximate analytical methods, such as multiple scales method [21,22], Krylov-Bogoliubov-Mitropolski (KBM) method [23,24], method of variation of parameter [25,26] and Lindstedt-Poincaré technique [27,28], can be used to approximately solve nonlinear ordinary equations. However, the first three methods are mainly used to find approximate periodic and quasiperiodic solutions. It is difficult for them to construct approximate analytical solutions of invariant manifolds. These invariant manifolds, including stable/unstable manifolds and homoclinic/heteroclinic orbits, play a key role in the dynamics and reconfiguration of formation flying. Besides, the method of variation of parameter can only be employed to find approximate solutions of invariant manifolds of equilibrium instead of periodic orbits. The Lindstedt-Poincaré method is first developed to compute periodic orbits of weak nonlinear equations based on the solution of the corresponding linearized equation. It has been successfully extended to studies on stronger nonlinear problems [29,30], including the computations of periodic and quasiperiodic orbits (center manifolds) [31,32], hyperbolic manifolds [33], homoclinic and heteroclinic connection [34], and even their bifurcation predictions [35]. Therefore, the Lindstedt-Poincaré method is employed to describe the various relative motions of Coulomb formation flying for the parameterized formation design (center manifolds), maintenance, and reconfiguration (hyperbolic manifolds).
This paper is first devoted to the arbitrary high-order expansion of relative motion for Coulomb formation flying, which greatly facilitates its formation design, maintenance, and configuration. The series expansion of arbitrary m : n-period orbits for the general dynamics system is first given. Finally, the analytical approximations can be regarded as a good replacement of numerical results with a tolerable error and greatly reduces the heavy computation.
The remainder of this paper is organized as follows. First, the dynamical model of a typical eight-charge Coulomb formation [17,36] is introduced in Sect. 2. Then, Sect. 3 constructs high-order series solutions of Lissajous orbits and m:n-period orbits around equilibrium points in Coulomb formation system by employing the Lindstedt-Poincaré method. In Sect. 4, the constructed high-order analytical solutions are applied to derive planar and vertical Lyapunov orbits, general Lissajous orbits, and Halo orbits as well as general m:n-period orbits. Furthermore, the high-order analytical solutions are compared with the corresponding exact numerical solutions to verify their effectiveness. Finally, a practical convergence domain of the analytical solutions is given for applying the analytical solutions to practical formation mission design. Finally, Sect. 5 makes some concluding remarks.

Formulation of the problem
This section proposed a typical dynamical model of Coulomb formation, and then, the problem description based on the proposed model is addressed. Considering an eight-charge Coulomb formation, the chief spacecraft has an independent Coulomb propulsion system, where eight controllable charged spheres are evenly distributed on eight vertices of a cube around the chief as shown in Fig. 1. The mass of the chief is far larger than that of the deputy, where the deputy is modeled as an ideal charged sphere. The chief can manipulate the relative motion of the deputy by controlling the eight charged spheres.
The dynamics of Coulomb formation is generally derived based on CW equations. The relative motion of the deputy with respect to the chief spacecraft is described in the Hill frame. As shown in Fig. 1, the origin of the Hill frame is at the centroid of the chief. The x-axis is directed from the chief radially outward, the z-axis is normal to the orbit plane and positive in the direction of the angular momentum vector, and the yaxis completes a right-hand triad. r = [x y z] T denotes the position of deputy; R i = [X i Y i Z i ] T denotes the constant position of the ith charged sphere of the chief (i = 1, 2 …, M). M = 8 is the number of controllable charged spheres. Then, the relative motion of the deputy with respect to the chief can be formulated as where the right-hand side of (1) represents the sum of electrostatic forces between two spacecraft, r i = [x i y i z i ] T = r − R i is the position to deputy from the ith charged sphere, ω is the constant orbital rate of the chief, m is the deputy mass, k c = 8.99 × 10 9 Nm 2 /C 2 is the Coulomb's constant, λ d is the Debye length over which the charge carriers are electrically screened, q is the deputy charge, and q i is the ith active charge of the chief. It is denoted as where τ is timelike nondimensional variable and X represents the state variables of (1). It is noticed that Then, (1) can be normalized as where C is called the Jacobi integral and thus (4) represents a conservative system. Defining momenta as p x = x − y, p y = y + x, p z = z , we can formulate the dynamical equations in a Hamiltonian form as It is well known that approximate analytical solutions are essential for the configuration design of formation flying. For traditional formation with chemical propulsion, the linear solution of the CW equation describes the relative motion between spacecraft and thus is the cornerstone of dynamics and control of space formation flying. Equation (4) shows that the linear solution is not accurate enough to describe the relative motion of Coulomb formation due to its high nonlinearity and inter-satellite coupling feature. Considering the complicated phase space of Coulomb formation flying, the objective of this paper is to construct several approximate analytical solutions describing its various relative motions. More precisely, a dynamical system of formation flying with (4) is given. Eight identical charged spheres are located on eight vertices of a cube with a side length of L as shown in Fig. 1. According to [13], the detailed parameters of this formation are defined in Table 1. Then, as a nonlinear extension of the linear solution in the CW equation, a general approximate analytical solution of (4) up to arbitrary order should be found in the following form [19] x where x i jkm is determined coefficient of the solution projection in the x-axis o f Hill frame, ρ xy and ρ z are amplitudes in x y-plane and z-direction, respectively, θ 1 and θ 2 are the corresponding phases. At the same time, the accuracy of these solutions should be high enough to replace numerical integration solutions and satisfy the designated requirement of formation design and reconfiguration.

Lindstedt-Poincaré method
In this section, the high-order expansions for a particular eight-charge Coulomb formation are investigated by means of the Lindstedt-Poincaré method. Based on the given parameters, equilibrium points of (4) can be derived by The equilibrium points on the X -axis and Y -axis can be obtained, as shown in Table 2 [17].
Without loss of generality, the Taylor expansion of the Hamiltonian function in (6) around the chosen equilibrium point is given as: Lindstedt-Poincaré method is an iterative calculation technology. The main principle is to observe the fact that the nonlinear terms of the system affect the Table 2 The equilibrium points on the X -axis and Y -axis and the type of their linearized vector field Equilibrium points x y z Linear behavior basic frequency of the corresponding linear system, and then the relationship between the frequency and amplitude of the system is given. Starting from the linear basic solution of the system, it continuously adjusts the relationship between the frequency and amplitude by iterating the known low-order solution step by step to obtain a higher-order series solution.
In order to perform the computations using the Lindstedt-Poincaré method, the motion equations of Coulomb formation are expanded in power series as where c x,i jk , c y,i jk , and c z,i jk are the coefficients of (i, j, k)-th item of the power series of −φ x , −φ y and −φ z , respectively.

Lissajous orbits
In this subsection, the analytical solution of Lissajous orbits is derived. The first step is to find the linear solution. If the higher-order terms on the right side in (10) are ignored, the linearized equation of the system is The part of solution of center manifolds for (11) can be obtained as where 2ω 0 , α and β are the amplitude in XYplane and Z axis, respectively. φ 1 and φ 2 are the corresponding phases, respectively. When the nonlinear terms are considered, the analytical solution around the equilibrium can be expanded in the power series about α and β.
where θ 1 = ωt + φ 1 and θ 2 = νt + φ 2 . Due to the nonlinear terms, the frequency of the system is no longer constant, but should also be expanded in the power series about α and β, i.e., In (13) and (14), i and j ∈ N, k and m ∈ Z. Due to the symmetry of the system, k and m have the same parity as i and j, respectively. Besides, due to the symmetry of sine and cosine functions, it can be assumed that k ≥ 0, and m ≥ 0 when k = 0. The series of ω and v only include those even items.
Our goal is to compute the coefficients x i jkm , y i jkm , (13) and (14) up to finite order N .
The Lindstedt-Poincaré method calculates these coefficients following an iterative scheme from the linear term. Compared to the solution of the linear part (12), we can know x 1010 = 1, y 1010 = κ, z 0101 = 1, ω 00 = ω 0 , v 00 = v 0 . By substituting this linear solution into (10), the coefficients of the second-order solution can be derived. Similarly, when the coefficients up to order N − 1 are obtained, i.e., x(t), y(t), and z(t) are determined up to order N − 1, ω and v are determined up to order N − 2. Substituting them into the right side of (10), we can obtain three power series up to N order, denoted by p, q, and r . Here, what we are interested in is those N -order terms. Without losing generality, the N -order terms of p, q, and r are denoted by p i jkm , q i jkm , and r i jkm (i + j = N ), respectively, and the N -order terms of ω and v are denoted by ω i j and v i j . Next, we analyze the composition of N -order terms on the left side of (10). According to (13), the derivatives of variable x can be given as Similarly, we can obtain the derivatives of y and z. Let f g denote the first-order derivate term, where f represents the power series of the frequencies (ω i j and v i j ) and g represents the coordinate variables (x, y, and z). i f and i g denote their corresponding order. Then, the f g satisfying i f + i g = N is the N -order terms we need. When i f is 0 or N − 1, the corresponding i g is N or 1 and f g is an unknown term that needs to be solved. When i f = 1, 2, …, N − 2, f g is a known term, which needs to be moved to the right-hand side of (10). Table 3 summarizes the unknown and known terms of the first derivatives, where δ i j denotes Kronecker function.
For the second derivatives of x, it can be similarly summarized as shown in Table 4. Then, we move all the known terms to the right-hand side of (10) to add them to p i jkm , q i jkm , and r i jkm , which are re-denoted byp i jkm ,q i jkm , andr i jkm . Besides, we should further discuss the calculation of unknown term f i−1 j and f i j−1 . In fact, they are composed of the unknown term 2ω 00 ω i−1 j (2v 00 v i−1 j ) and the remaining known terms, i.e., f i−1 j = 2ω 00 ω i−1 j + . For the second derivatives of y and z, similar results can be obtained. In summary, the linear equation of N -order unknown terms is yielded: Table 3 Derivatives of x and y with respect to time Table 4 Second derivatives of x, y, and z with respect to time where = kω 00 + mv 00 and When  (16) is zero, and thus z i jkm is set zero.
Remark 1 Due to the effect of small divisors, the series solution (13) by the Lindstedt-Poincaré method will diverge when the frequencies of XY -plane and Z -direction approach or are at resonance. When the frequencies are at m:n resonance, m:n-period orbits bifurcate from Lyapunov periodic orbits. In this case, the series solution (13) cannot accurately describe the motion and another applicable series solution need to be provided.

m:n-periodic orbits
It is known that the system frequencies vary with the amplitudes α and β, and thus, the frequencies will be resonant when ω/v is m/n (m, n ∈ Z + ) at a particular point. In this case, the relative trajectories of deputy spacecraft will become m:n-period orbits. In this subsection, we discuss the analytic solution of m:n-period orbits in the case that the frequencies ω and v are resonant. Since m:n-period orbits appear as a result of the effect of nonlinear frequencies, they cannot be obtained in the linear equation (11). To apply Lindstedt-Poincaré method to the high-order solution of m:n-period orbits, the form of (10) should be adjusted by adding a linear correction term Δ about z. Then, the equation becomes Our goal is to find one-dimensional torus (m:nperiod orbits) while the condition Δ = 0 is satisfied. In the iterative process up to finite order using the Lindstedt-Poincaré method, the correction term coefficient Δ is also expanded into a power series about α and β, i.e., Δ = i, j≥0 d i j α i β j , where its coefficients are also obtained by iterative computation. Analogously to the procedure to find the analytical solution of Lissajous orbits, the iteration starts from the linear solution of the system and the linear solution is where ω 0 : v 0 = m : n, φ is an arbitrary phase, d 00 = d − n 2 m 2 ω 0 2 . It is known that the amplitudes α and β are independent for Lissajous orbits. However, the m : n-period orbits, as one-dimensional torus, are determined by only α or β. α and β are coupled and their relationship is constrained by the condition Δ = 0, i.e., α = α(β). According to Δ = 0, we have d 00 = 0 for linear solution (20). However, this does not meet the necessary condition for one-dimensional torus, which means that the m : n-period orbit has no linear solution. Further considering the nonlinear term in (18) and the symmetry of sine and cosine functions, the analytical solution around the equilibrium can be expanded in the power series as Then, the Norder analytical solution can be obtained by calculating the coefficients x i jk , y i jk , z i jk , Ω i j , and d i j up to N order. Similarly, k has the same parity as (i + j) due to the symmetry of the system. The series of Ω and Δ only include those even items. For the coefficients of linear part in (21), compared to (20), we can know x 101 = 1, y 101 = κ, z 011 = 1, Ω 00 = ω 0 /m, d 00 = d − n 2 m 2 ω 0 2 . By substituting this linear solution into (18), the coefficients of the second-order solution can be solved using the Lindstedt-Poincaré method. Similarly, when the coefficients up to order N − 1 are obtained, i.e., x(t), y(t), and z(t) are determined up to order N − 1, ω and Δ are determined up to order N − 2. Substituting them into the right side of (18), we can obtain four power series up to N order, denoted by p, q, r , and Δz. Similarly, we are interested in the N -order terms of p, q, r , and Δz, which are denoted by p i jkm , q i jkm , r i jkm , and d i j (i + j = N ), respectively, and the N -order terms of Ω denoted by Ω i j .
To compute the N -order power series of m:n-period orbits, we should also analyze the composition of the N -order term on the left side of (18). Different from Lissajous orbits, the N -order unknown term of Δz appears in m : n-period orbits. The unknown and known terms of the first derivatives for m : n-period orbits are summarized in Table 5. The known terms are added to the corresponding series p i jk , q i jk , and r i jk and re-denoted byp i jk ,q i jk , andr i jk . The second derivatives can also be similarly summarized in Table 6, where f represents the series corresponding to Ω 2 and ∂ 2 v/∂θ 2 represents the second derivatives of x, y, and z with respect to θ . Table 5 Derivatives of x and y with respect to time for m : n-period orbits Table 6 Second derivatives of x, y, and z with respect to time for m : n-period orbits

Finally, the linear equation of N -order unknown term is yielded
Known term 1,2,..., When k = 1, the coefficients x i jk , y i jk , and z i jk can be solved immediately. When k = 1, x i jk , y i jk are coupled, and the coefficient of z i jk is zero. Thus, x i jk and z i jk can be set zero, and then y i jk and Ω i−1 j are solved by the first two equations of (22). d i j−1 are solved by the last equation of (22).
gives the approximate analytical condition of a bifurcation of m : n-period orbits from planar or vertical Lyapunov periodic orbit. The higher the order of series, the higher the accuracy of the bifurcation point. Based on this analytical condition, the bifurcation point and bifurcation direction can be quickly calculated in the continuation of periodic orbit families.

Remark 3
The proposed algorithm for approximate analytical computation of Lissajous and Halo orbits is mainly implemented using MATLAB code. To improve the computation efficiency, the part of the expansion of the motion equation (4) is implemented using differential algebraic techniques with C++14 in Visio studio 2017 [37].

Simulation results
This section implements the numerical simulations to verify the accuracy of the proposed method, and to determine the magnitude of errors when describing typical orbits for Coulomb formation flight. For this reason, the numerical simulation of Lissajous orbits is first carried out. Then, a numerical simulation of m:n-period orbits is provided for a supplement of the part of relative motion that the series solution of Lissajous orbits cannot describe.

Lissajous orbits
The series expansion of Lissajous orbits around L 4 in the Coulomb formation system is computed in this subsection. Table 7 gives their coefficients up to order 3. When α = 0 and β = 0, the series solution (13) corresponds to a vertical Lyapunov orbit near the equilibrium. Figure 2 shows a family of vertical Lyapunov . Taking the series solution as the initial value, the exact numerical solution with the corresponding amplitude can be obtained by using the differential correction [38]. Figure 3 shows the comparison of the time-varying positions of the linear analytical solution, the 10-order analytical solution, and the accurate numerical solution in an orbit period. Here, the numerical solution of a periodic orbit is defined as follows. If the numerical solution satisfies the boundary condition where x is the numerical solution, t 0 is the initial time, and T is the orbital period. Considering the floatingpoint precision of the numerical computation, the value of ε is set to 1 × 10 −12 .
It can be observed that the linear solution ignores the motion in the X and Y directions, and thus, it can only approximate the vertical Lyapunov when the amplitude is very small. However, the 10-order analytical solution almost coincides with the exact numerical solution. Finally, Fig. 4 shows the error distribution of the analytical solution of the vertical Lyapunov orbits with the order from 1 to 20 and of the amplitude from [0.001, 0.7]. The higher-order analytical solution can also be well approximated to the vertical Lyapunov with large amplitude.
When α = 0 and β = 0, the series solution (13) corresponds to a planar Lyapunov orbit near the equilibrium. Figure 5 shows a family of planar Lyapunov orbits with the order N = 30 and the amplitude α ∈  Fig. 6. It is shown that the error of the linear solution is very large, while the 10-order analytical solution almost coincides with the accurate numerical result, and its error magnitude is 10 −3 m. Although the analytical solution of planar orbits has high accuracy, it is nearly three orders of magnitude lower than the analytical solution of vertical Lyapunov orbit (10 −6 m). Figure 7 shows the error distribution of the analytical solution of the planer Lyapunov orbits  Fig. 4 with Fig. 7, we can conclude that for the Coulomb formation system, the accuracy of the vertical Lyapunov orbit is much higher than that of the planar Lyapunov orbit in the case of the same order of the analytical solution and the same amplitude. This is because, for an unstable periodic orbit, the deviation between the series solution and the numerical solution will enlarge rapidly with time. On the contrary, the vertical Lyapunov orbit is equivalent to a simple harmonic motion in the Z direction, and its stability is high, so its analytical solution has the highest accuracy.  When α = 0 and β = 0, (13) corresponds to a general Lissajous orbit near the equilibrium. Figure 8 shows a Lissajous orbit near L 4 with the order N = 20, and the amplitudes α = 0.08, β = 0.1. The comparison with the corresponding exact numerical solution and its position error in a time interval of dimensionless length π are shown as in Fig. 9. It should be noted that the Lissajous orbit of the exact numerical solution is very difficult to calculate due to the hyperbolic characteristics of L 4 . The numerical solution here is obtained by taking the initial value provided by the analytical solution with N = 10 as the initial condition, and then numerically integrating the Coulomb formation system equation. Figure 9b indicates that even if the initial error of the analytical solution is very small, it will increase exponentially with the evolution of the orbit, and the error expansion speed in the X and Y directions is much greater than that in the Z direction.
Finally, by comparing the higher-order analytical solution with the exact numerical solution, the convergence region of the higher-order analytical solution can be given. Specifically, for an analytical solution with the given order, its initial condition is obtained after the amplitudes α and β are given, and then the dynamic equation is numerically integrated with the normalized integration time T = π (about one orbital period). By comparing the position of analytical solution and exact numerical solution at the final time T f = π , the ana- lytical solution up to order N is regarded as practically convergent if the Euclidean norm of position error is less than 10 −3 m. It should be noted that due to the small divisor problem, the analytical solution obtained by the Lindstedt-Poincaré method is not convergent in fact. "Practical convergence" refers to that the error of analytical solution during one orbital period is less than a certain threshold, which is used to describe the degree of the high-order analytical solution approaching the exact numerical solution. The reason why the integration time of one orbital period is selected here is that the orbit of the equilibrium is unstable, and even the accurate numerical solution diverges rapidly after several periods. Therefore, comparing the error during one period is enough to illustrate whether the analytical solution converges practically. For the selection of the threshold value, other values can also obtain similar results. Figure 10 shows the practical convergence domain of the analytical solution of Lissajous orbits up to different orders concerning the amplitudes α and β. The results indicate that the convergence region of the series solution obtained by the Lindstedt-Poincaré method enlarges with the increasing order. On the one hand, the convergence domain shows the validity of the construction of the high-order analytical solution. On the other hand, it also provides a reference benchmark for when to use the analytical solutions. If the required accuracy is less than or equal to the specified accuracy of the convergence domain, the orbit can be approximated by the analytical solution in the convergence domain. For the case outside the convergence domain, the analytical solution cannot be directly applied, but at this time, it can still provide a good initial value to obtain the exact numerical solution by employing the differential correction.

m : n-periodic orbits
In this subsection, the series expansion of m:n-period around L 4 in the Coulomb formation system is studied. First, when (m, n) = (1, 1), the analytical solution (21) represents a typical Halo orbit. Figure 11 shows the family of Halo orbits with N = 20 and α ∈ [α min , 0.3], where α min = 0.137099841192699. α min can be obtained by solving the equation Δ =  Fig. 12 when N = 20. The convergence of the series solution of Halo orbits is also investigated. Fixed amplitude α = 0.2, when the order N is from 3 to 20, the corresponding position error of the series solution of the Halo orbit is shown in Fig. 13. It is seen that when the order is small, the position error is large and reaches the magnitude of 1 m, while the error is less than 0.1 mm when the order is N = 20, which can fully meet the requirements for the practical formation mission. To further evaluate the accuracy of the analytical solution, the initial residual acceleration (IRA) of Halo orbits is computed as follows. First, the spacecraft moves along the orbit determined by the series solution in a given time interval, and the corresponding analytical acceleration can be obtained by second-order derivation of the series solution. Then, the numerical acceleration is obtained by substituting the state of the spacecraft into (4). Finally, the average residual acceleration is the average value of the deviation between the numerical acceleration and the analytical acceleration of the spacecraft during this time. Figure 13 shows that the higher the order of the series solution, the smaller the average residual acceleration, which further demonstrates the validity of the high-order analytical solution. Finally, the case of general m : n-period orbits is investigated. Figure 14 shows the relationship between the amplitudes α and β of m : n-period orbits around L 4 in Coulomb formation system under different resonance ratios. It can be observed that whether the resonance ratio is greater than or less than 1, there exists a periodic orbit with the minimum amplitudes. The red curves are those cases that the resonance ratio m : n > 1 : 1. The minimum β indicates that m : nperiodic orbits bifurcate from a vertical Lyapunov orbit. At this time, the amplitude α can be arbitrarily small. On the contrary, the blue curves are those cases that the resonance ratio m : n < 1 : 1, where the m : nperiodic orbits bifurcate from the planar Lyapunov orbit. Besides, we can see that when the basic frequency ratio is close to 1:1, only a very large α can produce m : n-periodic orbits, and the greater the resonance ratio m : n, the smaller the minimum amplitude α min where the periodic orbit bifurcates. Figure 15 shows the periodic orbits with resonance ratios of m : n = 4 : 3, 5: 4, 4: 5, and 5: 6, respectively. It shows that the projection of m : n-periodic orbits in XY plane is actually a planar Lyapunov periodic orbit, and has the properties of Lissajous orbit in the direction of Z axis. The resonance ratios in Fig. 15a and b are greater than 1:1. In these cases, the amplitudes of the orbits in XY plane are very small, and the amplitudes in Z axis are up to 50 m. When the resonance ratios in Fig. 15c and d are less than 1:1, the projection in XY plane is a typical large-amplitude planar Lyapunov orbit, and the projection in the X Z plane is no longer symmetrical.

Conclusion
This paper presents a semi-analytical approach based on the Lindstedt-Poincaré method, which is capable of constructing series solutions up to arbitrary finite order for the relative motion of Coulomb formation systems near equilibrium. By applying this method, we obtain detailed analytical expressions for a range of bounded orbits, including Lyapunov periodic orbit, Lissajous orbit, halo orbit, and arbitrary m:n-periodic orbit near the equilibrium L 4 . The series solutions also provide an approximate analytical condition for the bifurcation of m:n-period orbits from planar or vertical Lyapunov periodic orbits. The effectiveness of our proposed method is verified by comparing the series solutions with corresponding numerical solutions, which demonstrates that the proposed method yields highly accurate results. In particular, the practical convergence domain of the high-order series solution of Lissajous orbits is obtained, which is crucial for the design and analysis of close-proximity Coulomb formation systems in various astronautical missions.
The proposed analytical series solutions provide a useful tool for the design and analysis of Coulomb formation systems, and they offer numerous advantages over traditional numerical approaches. With the series solutions, the formation configuration for actual missions with specific accuracy requirements can be designed quickly and efficiently. The results of this study are expected to have practical implications for the development of future space formation missions.
In the further work, the proposed high-order analytical method can be extended to describe the more complicated relative motions of Coulomb formation systems, such as formation reconfiguration or timedependent cases with elliptical reference orbits.