Characterizing Periodic Variations of Atomic Frequency Standards via Their Frequency Stability Estimates

The onboard atomic frequency standard (AFS) is a crucial element of Global Navigation Satellite System (GNSS) satellites. However, it is widely accepted that periodic variations can influence the onboard AFS. The presence of non-stationary random processes in AFS signals can lead to inaccurate separation of the periodic and stochastic components of satellite AFS clock data when using least squares and Fourier transform methods. In this paper, we characterize the periodic variations of AFS using Allan and Hadamard variances and demonstrate that the Allan and Hadamard variances of the periodics are independent of the variances of the stochastic component. The proposed model is tested against simulated and real clock data, revealing that our approach provides more precise characterization of periodic variations compared to the least squares method. Additionally, we observe that overfitting periodic variations can improve the precision of GPS clock bias prediction, as indicated by a comparison of fitting and prediction errors of satellite clock bias.


Introduction
The onboard atomic frequency standard (AFS), commonly referred to as the satellite clock, is a crucial component of the Global Navigation Satellite System (GNSS) satellites. This is due to the close relationship between their characteristics and the system's geodetic performance [1][2][3]. While ground-based AFS processing methodologies such as the power-law noise model [3], the offset and drift model for atomic clocks [4], Allan variance (AVAR) [5], and Hadamard variance (HVAR) [6] are useful in dealing with onboard clock data, they may not be sufficient. The behavior of satellite AFS differs from that of ground AFS due to an integrated non-dispersive effect affecting the broadcast timing signals observed from the ground [2,5,[7][8][9][10]. For instance, sub-daily periodic components can cause a clock bias up to 7.80 ± 1.61 ns [5] (about 2.33 ± 0.48 m equivalent light travel distance). Therefore, a more comprehensive modeling of periodic variations in GNSS satellite clocks is necessary to optimize the current system and plan future improvements.
GNSS satellite clocks are aligned with terrestrial time (TT) [11][12][13][14], which is annually realized by the Bureau International des Poids et Mesures (BIPM) as "TT BIPM * * " (" * * " indicates the year of calculation) [15]. Although TT is currently the most accurate and stable timescale, the post-processed timescale is not suitable for real-time applications. The International Atomic Time (TAI), established by the Consultative Committee for Time and Frequency (CCTF) and maintained by BIPM after the adoption of the atomic definition of second in the 13th General Conference on Weights and Measures (CGPM), is probably the most reliable "real-time" timescale in application. It is generated by more than 400 atomic clocks spread globally as an EAL (échelle atomique libre, or free atomic timescale), and the frequency of the EAL is steered by measurements of 12 primary or secondary frequency standards [15]. While GNSS timescales (e.g., GPS [11], GLONASS [12], Beidou [13], Galileo [14], etc.) are aligned with TAI in the long-term frequency, each satellite system has its own timescale. Therefore, the International GNSS Service (IGS) developed its own timescale, the International GNSS Service timescale (IGST), for GNSS research and application. IGST is also aligned with TAI in the long-term frequency [16]. Additionally, both the timescale algorithm of merely onboard AFS [4] and the IGST generated from both onboard and ground AFS [16] use the same time offset-frequency offset-drift (or aging) clock model as TAI [15,17]. Although time and frequency metrology techniques are essential to GNSS applications, it should be mentioned that the GNSS technique plays a significant role in resolving demands for universal timescales following scientific and industrial progress [18].
The behavior of space AFS is different from that of ground AFS [3]; however, electrical, gravitational, atmospheric, and thermal variations can cause magnetic field [19,20], microwave energy [21], light [22], thermally induced frequency shifts [23] of the clock. For instance, oblateness of geopotential causes periodic variations whose fundamental frequency is 2 cpd (cycles per day) in GPS satellite AFS [7]. Therefore, a satellite clock should be considered as the integration of AFS and "non-dispersive effects of all satellite components that affect the broadcast timing signals as observed from the ground" [5].
In fact, Senior et al. reported the detection of periodic signals at n × [2.0029 ± 0.0005] cpd for n = 1, 2, 3, and 4 [5] in all GPS Block II and Block IIA cesium and rubidium and Block IIR and IIR-M rubidium clocks. The variations were also reported in Block IIF rubidium clock and Block IIF cesium clock by Montenbruck et al. [9] and Fan et al. [2], respectively. While Senior et al. point out that 12-hour variation has a satellite-type-dependent relation with the sun-spacecraft-earth angle [5] and Montenbruck et al. suggest that solar illumination is the root cause of 12-hour period Block IIF Rb variations, they can not explain the "robust" [5] difference (58 ± 22 s) between the fundamental frequency (a weighted mean) of sub-daily GPS periodics (n × [2.0029 ± 0.0005] cpd) and the GPS mean orbit period (2.0057 cpd) [24]. However, the reason for this discrepancy remains unexplained. In addition, these studies are based on Fourier transform and least squares methods: • Satellite clock bias includes both deterministic and random signals: while the spectra of deterministic signals are estimated from its direct Fourier transform, the spectra of stochastic processes are calculated from the Karhunen-Loeve transform of its covariance [25].  [27].
Instead of the frequency-domain Fourier transform and least squares methods, we study the characterization of GPS onboard AFS periodic variations using frequency stability in this paper. In Section 2, we demonstrate that if the AFS signal comprises sinusoidal and random signals, the AVAR and HVAR of the sinusoidal signal are independent of those of the random signals. Additionally, AVAR and HVAR of GPS AFS periodic variations are formulated, and the method of fitting periodic variations from AVAR and HVAR estimates are presented. In Section 3, the effectiveness of the proposed model is tested against both simulated and real clock data.

Methods
To characterize the periodic variations of GNSS satellite AFS using its frequency stability estimates, we derive the AVAR and HVAR of periodic variations in this section. During the process, we demonstrate that the frequency stability of periodic variations is independent of the frequency stability of the stochastic component of AFS signals.
Finally, we present the numerical method utilized in this paper for characterizing periodic variations of AFS through frequency stability estimates.

Allan and Hadamard Variances of Periodic Variations
Suppose the satellite AFS clock bias, x(t), is composed of stochastic and periodic signals, i.e., then the theoretical values of AVAR and HVAR of the frequency standard are: and respectively, wherex(t) is the random part of x(t), a superposition of power-law noise (PLN) processes. Next,x(t) denotes the periodic variations of the clock: For the sake of brevity, we define two auxiliary functions, ∆ 2 H(x, τ) and ∆ 3 H(x, τ): respectively. The averaging interval τ = mτ 0 , m should be a positive integer, and τ 0 is the sampling interval of the data. By the theory of statistics, for a fixed averaging interval τ, Given a PLN process, f α , where f denotes the Fourier frequency, f α shapes the power spectral distribution (PSD) of the PLN process. When α > −3, E[ ∆ 2 H(x(t), τ)|τ] = 0; that is: For signals consistin of PLN processes, f α , α > −3, and sinusoidal variations, the AVAR of the periodic variations is independent of the AVAR of the random signals.
By applying Theorem 1 to Equation (1), we can separate the AVAR of stochastic fluctuations from the one of sinusoidal variations: Since any continuous periodic signals in 2 -space can be decomposed as the summation of a series of sinusoidal signals, AVAR of periodic variations of the satellite AFS is independent of those of the stochastic vibrations of the clock. There are plenty of documents on the characterization of AFS clock random vibrations; we will focus on characterizing periodic variations of AFS clock using AVAR.
After expanding E σ 2 AVAR (x, τ) with the definitions ofx(t) and ∆ 2 H(x(t), τ) and Equation (1), we have: Since t goes from −∞ to +∞, the value of the initial phase, ϕ 0 , makes no difference in the computation of E σ 2 AVAR (x, τ) . By replacing ωt + ϕ 0 with ωt and taking the following relations into account, we formulate the AVAR of periodic variations as following: Similarly, the mathematical expectation of HVAR for the periodic variations can be formulated as following: Numerical computation of Equations (4) and (5) may result in negative values due to truncation errors. To prevent numerical evaluations of Equation (4) from being negative values, we recast Equation (4) as a summation of squares: Similarly, the mathematical expectation of HVAR for periodic variations, Equation (5), can be reformulated as following:

Characterizing Periodic Variations Using Frequency Stability Estimates
To evaluate periodic variations from frequency stability estimates of GNSS satellite AFS, we modify the stochastic model proposed in [28]: where "s.t." is the abbreviation of "subject to"; σ is a column vector of frequency stability estimates (when the word "vector" is used, we mean column vector by default): Subscript k is used as a generic form that indicates different kinds of frequency stability: The variable h is a vector of noise intensity coefficients, a is a vector of amplitudes of sinusoidal variations, and Φ b is a matrix comprised of Φ k (α i , τ): The independent variables of Equation (6) are h and a. To solve the optimization problem, we set ε = 0.0025 and Q = [diag(σ)] −2 .
Since ε is nonzero, the feasibility of Equation (6) depends on the values of frequency stability estimates. When Equation (6) is infeasible, which means the optimization problem has no solution, we use the following alternative model: The independent variables of Equation (9) are µ, ν, h, and a; µ and ν are auxiliary variables defined such that nonzero components of µ and ν indicate violations of the inequalities respectively. Details about the optimization models, Equations (6) and (9), and their numerical solutions can be found in references [28,29].

Results
In this section, we test the proposed model, Equations (4) and (5), with both simulated and GPS clock data using the method described by Equations (6) and (9).

Simulated Data
To verify Equations (4) and (5), we generate sinusoidal signals with frequencies, i × 2.0029, i = 1, 2, 3 cpa, and list their amplitudes in the second column of Table 1.
(The frequencies and amplitudes are set according to reference [5]). As demonstrated in Figure 1a-c, the AVAR computed from Equation (4) (represented by red dots in Figure 1) aligns with the one estimated from the simulated sinusoidal signals (represented by black solid lines in Figure 1). The same holds true for Hadamard variances (HVAR), though these results are not presented here in the interest of brevity. In Figure 1d, AVARs estimated from simulated clock data and their "theoretical" values are represented as a black solid line and red dots, respectively. The "theoretical" AVAR is computed as the sum of Equations (7) and (4), and the simulated satellite clock biases are comprised of white PM, flicker PM, white WHFM, FLFM, RWFM, and sinusoidal variations of 12-, 6-, and 4-hour cycles (whose AVARs are shown in Figure 1a-c, respectively). As shown in Figure 1d, the two AVARs are close to each other, except for the values around averaging times 10 4 s. It can be observed from Figure 1a-c that the AVAR of 12-, 6-, and 4-hour sinusoidal variation reaches its maximum near averaging times 5 × 10 3 , 10 4 , and 2 × 10 4 s, respectively. Discrepancies between the two AVARs can be explained by the interactions between stochastic and periodic variations due to the finite data set. Furthermore, to compare the effectiveness of detecting periodic variations from satellite clock biases directly using the least squares method and from Equations (4) and (5) using Equations (6) and (9), the periodic variations in the simulated data are fitted and removed using both methods. The amplitudes given by the least squares method and Equations (6) and (9) are listed in the third and fourth columns of Table 1, and AVARs of the resulting periodic variation-free clock biases are plotted as a red dash-dot line and blue dots, respectively, in Figure 2. In the least squares estimation process of this article, the satellite clock biases are are firstly detrended by fitting and removing a second-order polynomial. Then, a 12-hour sinusoidal signal was fitted and removed, followed by a 6-hour sinusoidal signal, a 4-hour frequency signal, and, finally, a 3-hour sinusoidal signal. While the least squares method gives greater periodic variation amplitudes than the set values of the simulated variation amplitudes and the ones estimated by Equations (6) and (9), AVARs estimated from simulated satellite clock bias with periodics fitted and removed by the least squares method does not seem to be periodic-variation-free. It can be observed in Figure 2 that in comparison to the AVARs estimated from the simulated data (plotted as a black solid line in Figure 2), the "lump" near averaging intervals 5 × 10 3 suggests an extra 4-hour variation. In other words, the least squares method overfits the periodic variations. On the other hand, AVARs computed from simulated data with periodic variations fitted and removed using Equations (6) and (9) are the most periodic-variation-free AVARs among the three. To have a closer look at the overfitting problem, the root mean square (RMS) of the detrended simulated clock bias, detrended simulated data with periodic variations fitted and removed using the least squares method, and detrended simulated data with periodic variations fitted and removed using Equations (6) and (9) are computed: the result are 8.02 × 10 −9 , 7.80 × 10 −9 , and 7.95 × 10 −9 ns, respectively. While the least squares method provides the lowest fitting errors, it seems to overfit flicker noise with a 4-hour periodic variation.

GPS SVN63 Clock Data
To verify the effectiveness of Equations (4) and (5), we computed the AVAR of the GPS PRN01 (SVN63) onboard rubidium AFS, which is shown as a black solid line in Figures 3 and 4. The IGS final combined precise clock and orbit data from MJD 56739.0 to MJD 56745.9965 [30,31] are used in the computation. In this paper, the preprocessing of all IGS satellite clock data includes: • Removing the J2 relativistic effect using Kouba's method [7] and cubic spline interpolation of IGS final combined orbit data; • Removing the day boundary discontinuities using Yao et al.'s methods [32].
For comparison, AVARs of 12-, 6-, and 4-hour periodic variations are plotted as a purple dash-dot line, blue dots, and a red dash line, respectively, in Figure 3. The values of these AVARs are normalized to be tangent to the AVAR of SVN63 satellite AFS at the maxima of their AVARs.  It can be observed from Figure 3 that the local maxima and minima of AVARs estimated from SVN63 clock bias and the 12-hour sinusoidal coincide for averaging time τ ≥ 10 4 s, suggesting the influence of 12-and 6-hour periodic variations. This consistency reinforces the conclusion that GPS onboard atomic frequency standards are influenced by 12-hour periodic variations. To quantify the influence of periodic variations, the variations are fitted and removed using both the least squares method and Equations (6) and (9); AVARs computed from the resulting clock biases are represented as a red dash-dot line and blue dots, respectively, in Figure 4. It seems that the least squares method may overfit the 12-, 6-, 4-, and 3-hour periodic variations of the GPS Block IIR rubidium AFS with frequency noise: gaps between the black solid line (AVAR estimated from IGS final combined SVN63 clock data) and the red dash-dot line (AVAR resulting from the least squares processing) keep increasing after averaging time 2 × 10 4 . On the other hand, after processing with Equations (6) and (9), the sigma-tau plot of the resulting AVAR estimates is the closest to the classical sigma-tau plot of the AVAR calculated from the summation of white PM, flicker PM, white WHFM, FLFM, and RWFM noise processes.   (6) and (9), respectively.

Other GPS Clock Data
To verify the general applicability of Equations (4) and (5), we compute the AVARs of GPS PRN02∼32 onboard AFS using IGS final combined precise clock data, clock biases with periodic variations removed by the least squares method, and Equation (9); the resulting AVARs are drawn as a black solid line, a red dash-dot line, and blue dots, respectively, in Figures 5 and 6. During the time-span of IGS clock data used in this study (from 23 March 2014 to 27 December 2020), several satellites, such as PRN04, PRN14, PRN18, PRN23, and PRN32, were replaced. In addition, some PRN satellites, such as PRN04, have been substituted several times. For brevity, we only chose the oldest and most recent SVN satellites of these satellites from the dataset.  MJD.57451~57452 clock bias periodic variations removed using least square periodic variations removed using Eqn (9) Figure 6. AVARs of GPS PRN18∼32 onboard clock bias 12-, 6-, 4-, and 3-hour sinusoidals fitted and removed using least squares method and Equation (9), respectively.
By comparison with Figure 1, we can make the following observations from Figures 5 and 6: • For more than half of the satellites, the sigma-tau plot of the AVAR computed from clock bias with periodic variations removed using Equation (9) is closer to the standard sigma-tau plot of AVARs than the least squares method. In other words, the optimization model Equation (9) is "conservative": it tends to underestimate periodic variations of GPS satellite AFS. In effect, a large portion of the numerical solutions of 4-and 3-hour variations returned by solving Equation (9) are zero. According to the theory of statistics, the least squares method makes huge assumptions about data structure and produces stable results. On the other hand, Equation (9) only requires the least violation of the 95% confidence intervals by the AVAR estimates; it does not require optimal estimation of the amplitudes for periodic variations. By comparison with AVARs computed using the least squares method, as shown in Figures 5 and 6, Equation It can also be observed from Figures 4-6 that the AVARs computed from clock biases with periodic variations fitted and removed using Equation (9) are always smaller than the AVARs computed from IGS final combined clock biases, and the two coincide at averaging intervals greater than 4 × 10 4 s. For most cases, this also holds for AVARs estimated from clock bias with periodic variations fitted and removed using the least squares method. There are only two exceptions: • The gap between the AVAR estimated from PRN01 clock bias with periodic variations fitted and removed using the least squares method and the AVAR computed from PRN01 clock bias enlarges with increasing averaging interval. Since the tail of AVAR estimated from PRN01 clock bias with periodic variations removed using the least squares method has a similar shape to the AVARs of 12-hour sinusoidal variations, and the AVAR computed from PRN01 clock bias with periodic variations fitted and removed using Equation (9) suggests strong frequency noise, the discrepancies between the AVAR estimated from PRN01 clock bias with periodic variations fitted and removed using the least squares method and the AVAR computed from PRN01 clock bias is caused by overfitting the periodic variations by taking a portion of frequency noises as 12-hour variation. • The AVAR estimated from PRN22 clock bias with periodic variations fitted and removed using least squares method is greater than the AVAR computed from IGS final combined PRN22 clock bias around averaging time 2 × 10 4 s. Since AVARs computed from the three PRN22 clock biases increase with the averaging interval for τ ≥ 2 × 10 4 s, PRN22 AFS is influenced by strong FM noise processes. It seems that the least squares method overfits the periodic variations of PRN22 by taking a portion of frequency noises as 12-hour variation.
Therefore, Theorems 1 and 2 hold for GPS onboard AFS.
Overfitting of periodic variations can sometimes help improve GNSS onboard clock prediction, as shown in Table 2. By taking frequency noises as sinusoidals, the RMS of one-day GPS onboard atomic clock prediction under two-day observation is decreased. However, removing periodic variations based on Equations (6) and (9) can increase the prediction RMS. This could be due to several reasons:

•
Overfitting of periodic variations can reduce the clock residuals caused by powerlaw noise processes. When periodic variations are removed, the interaction between random clock behaviors and periodic variations is suppressed, leading to an increase in clock residuals and prediction RMS. • Only high-variability estimates (HVAR) are used in solving Equations (6) and (9), which may not capture all the periodic variations present in the data. The maximum averaging time of HVARs estimated from two-day GPS clock bias is 1.44 × 10 4 s, while the first local minimum of Equation (5) appears at averaging interval τ = 2 × 10 4 s. This means that some periodic variations may not be captured by HVAR estimates and could contribute to an increase in the RMS prediction when removed. Table 2. One-day GPS clock bias prediction RMS (root mean square) of fitting and removing first-or second-order polynomial 12-, 6-, 4-, and 3-hour sinusoidals against two-day observations by using least squares and solving Equations (6) and (9), respectively. (6) and (9)

Conclusions and Discussion
The detection of periodic variations in GNSS satellite AFS presents a dilemma: while the spectra of deterministic sinusoidals can be estimated directly from the Fourier transform, those of power-law noise processes, which are random, must be calculated from the Karhunen-Loeve transforms of their covariance. Currently, time domain methods are generally used in the evaluation of AFS and timescales. In this paper, we characterize the periodic variations of GNSS satellite AFS clock bias using AVAR and HVAR. In addition, we demonstrate that the AVAR and HVAR of sinusoidal signals present in the clock bias are independent of the AVAR and HVAR of the random vibrations of the satellite AFS. The method to detect periodic variations using these characterizations is given, and it is test against both simulated and real data.
The proposed method characterizes the periodic variations of satellite AFS more accurately than the least squares method, as demonstrated in both simulated and real clock data tests. Therefore, the method developed in this article can serve as a criterion for detecting periodic variations. The least squares method tends to overfit the 12-, 6-, 4-, and 3-hour frequency variations, resulting in sinusoidals that may have no physical meaning. However, the overfitting demonstrated by the least squares method suggests that by treating some of the frequency noise as periodic variations, we can improve GPS prediction precision in practice.
Author Contributions: Conceptualization, methodology, validation, and writing-original draft preparation, W.C.; writing-review and editing, J.Z.; data collection, writing-review and editing, and supervision, G.N. All authors have read and agreed to the published version of the manuscript.

Conflicts of Interest:
The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript; or in the decision to publish the results.

Abbreviations
The following abbreviations are used in this manuscript: