## 4 Numerical Simulations and Results  
We will now demonstrate the method for a specific problem. Because the Starlink constellation is a typical LEO mega-constellation with great research significance, we choose the orbit and mechanical environmental parameters of Starlink for simulation. Take into consideration the satellite with mean Keplerian elements as shown in Table 1 and mechanical environmental parameters as shown in Table 2 to conduct Monte Carlo simulations. Orbital and structural parameters of Starlink satellites can be found in SpaceX's FCC documents (such as 36 FCC Rcd 7995 (11)(2021)). In our simulations, the semi-analytical approach shows a strong advantage in computing efficiency. On a PC equipped with i7-1065G7 CPU 1.30 GHz, the numerical and the semi-analytical approach respectively take 4.258 s and 0.103 s under \(J_{2}\) and atmospheric drag for the 7-day orbit prediction. Without loss of generality, the initial position and velocity uncertainty of the satellite is assumed to be Gaussian distributed, and are defined in the geocentric inertial coordinate system. The mean is zero and the covariance matrix is  
\[ \begin{cases}\Sigma_{r}=\text{diag}(100,100,100)(\text{m})^{2},\\ \Sigma_{\
u}=\text{diag}(0.01,0.01,0.01)(\text{m}/\text{s})^{2}.\end{cases}  \tag{25}\]  
Four thousand points are distributed to be consistent with the initial uncertainty, generated and propagated using the averaged \(J_{2}\) and atmosphere drag dynamics, referring to the second equation of Equation (3). The secular terms of \(J_{2}\) and atmosphere drag are the same as Equations (4), (5), and (8). The conversion between osculating elements and mean elements abides by the first equation of Equation (3), and short-period terms and long-period terms are described in Equations (6) and (7). Figure 2 shows initial points distribution before and after conversion from osculating elements to mean elements. The mean track propagated in averaged \(J_{2}\) and atmosphere drag dynamics is marked as a green line, while the osculating track propagated in full dynamics--which includes short-period and long-period terms besides secular terms--is marked as an orange line. To better display the size of the uncertainty ellipsoid, the absolute position of the points is subtracted from the nominal position, and denoted by \(\delta\). It can be seen that it is the conversion between osculating elements and mean elements that makes the distribution of sample points show strong nonlinearity at initial moments. Figure 3 shows points distribution after different revolutions up to 100 orbits, which is about 6.2 days. The samples are gradually correlated with the mean trajectory as the propagation time progresses. After 1 orbit period, the shape of points takes on a shuttle shape. Hence, it is reasonable to assume that the position nonlinearity becomes weaker after several revolutions. The argument of latitude of a satellite is Gaussian distribution after going through the same revolutions. Table 3 shows the first two order moments of the relative phase of two co-plane satellites in the constellation, \(z\) value of hypothesis test Equation (20) and whether \(H_{0}\) should be accepted are also presented in this table. There is only a \(15^{\circ}\) difference in relative position between the two satellites, and other initial mean elements of them are the same. When the significance level of hypothesis test \(\alpha\) is given as 0.05, the \(z\) value of the two-sided test is 1.96. Thus, the hypothesis \(H_{0}\) is rejected in the beginning and is accepted after 1 orbit, because the standard deviation is growing faster than the mean over time. As a consequence, the phase of two co-plane satellites can also be considered as a Gaussian distribution with zero mean as propagation time goes by.