## 3 pLEO Time Scale Development
### The Basic Time Scale Equations  
In a pLEO constellation, where multiple clocks are present, the effects of clock errors or failures can be detected and mitigated by analyzing the behavior of clock phase discrepancies across the different relative phase measurements. Specifically, after completing the clock difference process for all clocks \(j=1,\ldots,N\) (where \(j\
eq i\)), we can investigate the individual clock phase discrepancies by solving a system of \(N\) simultaneous equations for the \(N\) phase discrepancies.  
\[\sum_{i=1}^{N}a_{i}\hat{e}_{i}(t_{k}) =0 \tag{56}\] \[\hat{e}_{i}(t_{k})-\hat{e}_{j}(t_{k}) =k_{ij}^{p}(t_{k+1})[z_{ij}(t_{k+1})-H_{ij}\Phi_{\delta}\hat{x}_{ ij}(t_{k})]\,.\]  
Thus, it leads to  
\[\hat{e}_{j}(t_{k})=\sum_{i=1}^{N}a_{i}k_{ij}^{p}(t_{k+1})[H_{ij}\Phi_{\delta} \hat{x}_{ij}(t_{k})-z_{ij}(t_{k+1})]. \tag{57}\]  
Now, it should be realized that of all the phase shocks estimated as in (57), the individual clock phases are given by  
\[\hat{p}_{j}(t_{k+1}) =\hat{p}_{j}(t_{k})+\delta\hat{f}_{j}(t_{k})+\frac{\delta^{2}}{2} \hat{p}_{j}(t_{k})+\hat{e}_{j}(t_{k}) \tag{58}\] \[=\hat{p}_{j}(t_{k})+\delta\hat{f}_{j}(t_{k})+\frac{\delta^{2}}{2} \hat{p}_{j}(t_{k})+\sum_{i=1}^{N}a_{i}k_{ij}^{p}(t_{k+1})[H_{ij}\Phi_{\delta} \hat{x}_{ij}(t_{k})-z_{ij}(t_{k+1})]. \tag{59}\]  
Interestingly, the above result highlights another compelling aspect of the proposed time scale problem: during each iteration, there is inadequate information to determine the phase states of individual clocks, as the prior frequency state, \(\hat{f}_{j}(t_{k})\), and the frequency drift state, \(\hat{r}_{j}(t_{k})\), are unknown. However, as shown in (55), the difference between frequency shocks is known. Therefore, it is essential to estimate the individual frequency shocks for each clock. A similar requirement applies to the individual frequency-drift random shocks for each clock.  
It is important to note that there are only \(N-1\) simultaneous equations for \(N\) unknowns. This ambiguity can be resolved by applying the same assumption about the frequency and frequency drift shocks that was previously made about the phase shocks, as shown in (56)  
\[\sum_{i=1}^{N}b_{i}(t_{k})\hat{\eta}_{i}(t_{k}) =0 \tag{60}\] \[\sum_{i=1}^{N}c_{i}(t_{k})\hat{a}_{i}(t_{k}) =0\,. \tag{61}\]  
Using the result (55), it is easy to see that both frequency and frequency drift random shocks can be obtained for each clock  
\[\hat{\eta}_{j}(t_{k}) =\sum_{i=1}^{N}b_{i}k_{ij}^{f}(t_{k+1})[H_{ij}\Phi_{\delta}\hat{x }_{ij}(t_{k})-z_{ij}(t_{k+1})] \tag{62}\] \[\hat{a}_{j}(t_{k}) =\sum_{i=1}^{N}c_{i}k_{ij}^{r}(t_{k+1})[H_{ij}\Phi_{\delta}\hat{x }_{ij}(t_{k})-z_{ij}(t_{k+1})]\,. \tag{63}\]  
Consequently, the solutions for frequency and frequency aging can be obtained by inserting the values for the random shocks from Equations (61) and (62) into Equation (46)  
\[\hat{f}_{j}(t_{k+1})=\hat{f}_{j}(t_{k})+\delta\hat{r}_{j}(t_{k})+\sum_{i=1}^{ N}b_{i}k_{ij}^{f}(t_{k+1})[H_{ij}\Phi_{\delta}\hat{x}_{ij}(t_{k})-z_{ij}(t_{k+1})] \tag{64}\]  
and  
\[\hat{r}_{j}(t_{k+1})=\hat{r}_{j}(t_{k})+\sum_{i=1}^{N}c_{i}k_{ij}^{r}(t_{k+1}) [H_{ij}\Phi_{\delta}\hat{x}_{ij}(t_{k})-z_{ij}(t_{k+1})]\,. \tag{65}\]  
The algorithm for the clock ensemble, as described in Equations (59), (64) and (65), begins with initial values for three states: \(p_{j}(t_{0})\), \(f_{j}(t_{0})\), and \(r_{j}(t_{0})\). For each iteration, the algorithm generates the phase, frequency, and frequency drift states by adjusting the weights associated with these parameters. The weights are set in such a way that they satisfy the limit theorem, which states that the sum of the weighted random shocks approaches zero as the number of clocks increases.