# Near-Real-Time Global Thermospheric Density Variations Unveiled by Starlink Ephemeris
## 2 Data and Methods
The precise force model incorporates the following components: 120-order gravity field model from EGM2008 [41], the solid Earth tide model from IERS2010 [42], and the ocean tide model from EOT11a [43]. Figure 2: Semi-major axis (SMA) variations of STARLINK-1009. (**a**) Osculating SMA computed from the reference batch. (**b**) Mean SMA derived from the reference batch and the subsequent batch. (**c**,**d**) The osculating SMA and mean SMA computed from continuous ephemeris formed by removing overlapping arcs from multiple batches, respectively, and color-coded data represent different batches. Between consecutive epochs, the satellite trajectory is approximated as a short arc. To ensure the accuracy of the gravitational work integration, an eight-point piecewise Lagrange interpolation scheme is employed to calculate both the Earth gravitational force \(\overset{\rightarrow}{\delta}_{ECI}\) and the satellite's geocentric position vectors \(\overset{\rightarrow}{r}\). Then, the mechanical energy loss (\(\Delta E\)) can be derived from the Starlink ephemeris data using Runge-Kutta methods. Notably, the mechanical energy loss of Starlink satellites at 500-600 km will be small during brief observational epochs; thus, errors introduced from imperfect force modeling approximations and cumulative numerical integration inaccuracies are magnified [44]. Nevertheless, our analysis focuses primarily on relative variations in mechanical energy loss. Through strict adherence to the consistent mechanical modeling throughout all computations, the method in this paper can mitigate the effect of the systematic errors while maintaining the computational efficiency required for near-real-time processing of massive satellite constellation data. It should be noted that Starlink satellites periodically perform orbital maneuvers using their onboard electric propulsion systems, which significantly increase mechanical energy and thus affect the accuracy of mechanical energy loss calculations. To address this, we implement the aforementioned maneuver detection method to exclude data points associated with maneuvers. After calculating the mechanical energy loss (\(\Delta E\)) caused by atmospheric drag, a 9-point moving window median filter is applied to suppress anomalous data points and enhance the statistical robustness. Then, the mean thermospheric density between two epochs (\(\Delta t\)) can be computed using the following equation [29]:  
\[\rho=\frac{\Delta E}{v\Delta t}/\left(\frac{C_{d}A}{2m}v_{r}^{2}\right) \tag{3}\]  
where \(C_{d}\) is the drag coefficient, \(A\) is the effective cross-sectional area, \(m\) is the satellite mass, and \(v_{r}\) represents velocity of relative motion between satellite and atmosphere. The accurate calculation of \(C_{d}A\) requires detailed knowledge of satellite attitude information and panel design specifications, including surface area and material properties. However, these parameters of Starlink satellites remain unavailable due to SpaceX's non-disclosure. Therefore, we utilize the mechanical energy loss of Starlink satellites as a proxy indicator for thermospheric density variations. In this study, we utilize Swarm-B's thermospheric density data computed by the European Space Agency (ESA), which can be obtained through ESA's data portal ([[https://swarmhandbook.earth.esa.int/catalogue/SW_DNSxPOD_2](https://swarmhandbook.earth.esa.int/catalogue/SW_DNSxPOD_2)]([https://swarmhandbook.earth.esa.int/catalogue/SW_DNSxPOD_2](https://swarmhandbook.earth.esa.int/catalogue/SW_DNSxPOD_2)), accessed on 13 April 2025), to validate the results of the mechanical energy loss derived from Starlink ephemerides. We adopt an inverse approach by converting the Swarm-derived thermospheric density into equivalent mechanical energy loss (\(\Delta E\)) through parameter substitution in Equation (3).