## 2 Materials and Methods
### Orbital Dynamics Model
#### 2.1.2 Reduced Dynamics Model for Climbing Starlink Satellites  
To investigate the orbital characteristics of Starlink satellites during the orbit-raising phase, a representative set of 22 satellites from the same launch on 24 January 2024 was selected from the TLE catalog. The variations in their orbit altitudes and inclinations with time during the orbit-raising process are shown in Figure 1. Among the selected satellites, 18 experienced two distinct orbital parking stages--approximately at 350 km and 430 km--before reaching their final operational altitude near 480 km. The duration of each parking phase varied among the satellites. In contrast, four satellites--identified by NORAD IDs 58837, 58844, 58845, and 58847--followed a different orbit-raising pattern. These satellites exhibited parking phases at approximately 430 km and 450 km before transitioning to their final orbital altitudes near 480 km and 550 km, respectively. Further analysis of historical TLE data from multiple Starlink batches indicates that, during continuous orbit-raising phases, satellite altitudes tend to increase approximately linearly over time, although at varying rates of ascent. This trend is consistent with the patterns observed in Figure 1a. Given the unknown timing of thruster activation and deactivation, along with the quasi-linear altitude increase observed over multi-day periods, the propulsion system is modeled as delivering a time-averaged, constant thrust acceleration. As illustrated in Figure 1b, the inclination variation among all 22 satellites remains within 0.01\({}^{\circ}\), suggesting that out-of-plane thrust components are negligible during the orbit-raising phase. Moreover, theoretical thrust optimization shows that the efficiency of altitude adjustment is maximized by the along-track thrust. Based on these considerations, the adopted dynamics model assumes a constant equivalent thrust acceleration in the along-track direction.  
Figure 1: _Cont._  
For LEO constellations like Starlink, the primary perturbative forces arise from the Earth's non-spherical gravitational field and atmospheric drag. Among these, atmospheric drag is the dominant factor causing semi-major axis decay, offsetting the orbit-raising effect of the along-track thrust. Considering this dynamic interplay, this study incorporates both atmospheric drag and along-track thrust into a unified estimation framework. Under this formulation, the dynamics model is simplified to three primary components: central gravitational acceleration, perturbations due to the Earth's non-spherical gravitational field, and a constant equivalent low-thrust acceleration, as defined in Equation (4). The perturbation due to the Earth's non-sphericity \(\mathbf{a}_{\text{NS}}\) is calculated using the EGM96S geopotential model, truncated at degree and order 30. The combined effect of thrust and atmospheric drag is modeled as a net equivalent acceleration in the along-track direction. Most of the numerical simulations in this study (Sections 3.2.1-3.2.3) are conducted using the simplified model, Equation (4). Only in Section 3.2.4 did we use Equation (21), which extends Equation (4) by incorporating atmospheric drag perturbations.  
\[\left\{\begin{array}{l}\frac{d\mathbf{r}}{dt}=\mathbf{v}\\ \\ \frac{d\mathbf{v}}{dt}=-\frac{\pi}{r^{3}}\mathbf{r}+\mathbf{a}_{\text{NS}}+M_{RTN\to ECI} \left[\begin{array}{c}0\\ a_{T}\\ 0\end{array}\right]\\ \end{array}\right. \tag{4}\]