## 1 Introduction  
The development of microsatellites leads to a general trend in Networking with LEO mega-constellations. Many countries and companies have proposed building this type of constellation, such as OneWeb and Starlink in the United States [1]. Of these, the Starlink communications constellation of 12,000 satellites is expected to be constructed around 2025, and accumulated to 42,000 satellites afterward [2]. Compared with the traditional middle and high Earth orbit navigation constellations, the LEO mega-constellations are characterized by low orbit height, dense satellite distribution, large perturbation, and high requirements for station position maintenance. Although the relative phase between satellites will not drift in constellations with the same value of mean elements, due to errors in the measurement data, dynamic model, solution algorithm, actuator, and other factors, the constellation configuration will be out of order in the uncontrolled state, leading to a high risk of collisions within the constellations. To investigate how much the relative phase between satellites in the same orbital plane in the LEO mega-constellations will be drifted by the perturbations, it is critical to first point out the nominal value of the constellation design. From the semi-analytic theory [3], the satellite orbits can be described by their initial mean elements. If the initial mean elements of the constellation are identical, the relative phases within the constellation will not theoretically diverge. So, the constellation nominal values can be considered a set of mean elements. Accordingly, the problem of phase evolution affected by perturbations can be viewed as the initial uncertainties propagation of the constellation: how to establish the mapping relationship of the terminal state to the initial state. The initial state is the deviation of the position velocity of the two satellites from the nominal design value, and the terminal state is the deviation of the relative phase of the two satellites from the design nominal value after a while. As widely recognized, orbit propagation is a fundamental problem in astrodynamics. There are three main approaches for modeling orbit propagation [4]: the numerical approach, analytical approach, and semi-analytical approach. The numerical approach [5; 6] has high precision and is suitable for scenarios with complex mechanical environments, but consumes more computational resources. The analytical approach [7] is developed on [PERSON]'s perturbation solutions [8; 9], using an approximate analytical solution instead of the original ordinary differential equations of motion, preserving the main features of the perturbation while neglecting higher-order small terms. The analytical approach makes the calculation much more efficient while losing some accuracy. The semi-analytical approach [10; 11] mixes the analytical and numerical approaches' advantages, classification of perturbation solutions into long-term, long-period, and short-period terms according to perturbation characteristics. [PERSON] [12] presents the computational efficiency strength of the semi-analytical approach in long-term dynamics. [PERSON] [13] presents a dynamic approach of satellite formation flights based on the THEONA semi-analytical theory with arbitrary values of the eccentricity. For general stochastic dynamical systems, the evolution equation of the probability density over time satisfies the Fokker-Planck equation [14] (FPE), while for orbital dynamical systems with high dimensionality (6 dimensions) and nonlinear perturbations, the FPE solution becomes very tough. Thus, the probability density function (PDF) evolution often requires the assumption of local linearization and only the lower order moments (e.g., mean and covariance matrices) are analyzed. Currently, nonlinear uncertainty propagation approaches include polynomial chaos expansions [15], state transition tensor approaches [16; 17], differential algebra approaches [18; 19], and Gaussian mixture model approaches [20].