Characteristics of a real-time ensemble time scale corrected by the KRISS-PSFS using a reduced Kalman filter (Kred) algorithm

An ensemble time, which is indicated by TA (temps atomique, in French), was generated every minute by combining four active hydrogen masers with the reduced Kalman filter (Kred) algorithm. The TA frequency was averaged over 12 h and then the mean frequency was steered to the primary and secondary frequency standards developed by KRISS with a micro-phase stepper (MPS). The steered and physically realized TA frequency is called KRISS ensemble time (KET). We investigated the frequency-state estimate outputted from the Kred algorithm and the control signal applied to the MPS, according to the abnormal temperature change of the room where the hydrogen masers were placed. When the frequency of the TA was corrected by the KRISS-F1, the Allan deviation of KET was 1.72 × 10−15 at the average time τ = 1 d, and the variation of [UTC–KET] was stable within 2 ns for 75 d.


Introduction
Many countries around the world have time standard laboratories that maintain the standard time for that country or region [1]. In the lab lab, commercial atomic clocks (e.g. cesium atomic clocks, hydrogen masers) are usually used to generate the local standard time UTC(lab). The lab lab contributes to the creation of Coordinated Universal Time * Author to whom any correspondence should be addressed.
Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.
(UTC) and International Atomic Time (TAI) by continuously reporting UTC(lab) data to the International Bureau of Weights and Measures (BIPM). Then BIPM announces, once a month through the FTP Server [2], information such as [UTC-UTC(lab)], weights of individual clocks contributing to TAI generation, frequency differences between individual clocks and TAI, frequency drifts of clocks participating in the computation of TAI, and so on.
Some laboratories have developed PSFS such as cesium atomic fountain clocks and/or optical clocks and use them as a reference for generating their local time UTC(lab) or TA(lab) [3][4][5][6][7]. These PSFS are used to measure the difference in frequency with a clock (mainly hydrogen maser) participating in TAI generation. When the measured values are submitted to BIPM, BIPM announces the frequency difference (indicated by d) between the PSFS and TAI through Circular T once a month [8]. This allows lab lab to have a frequency standard that contributes to generating the TAI from the EAL (free atomic time scale).
KRISS measures the frequency difference between the KRISS-Yb1 optical frequency standard and the active hydrogen maser (KRISS code: H28) used as a master clock of UTC(KRIS), and periodically submits it to BIPM. In addition, the time difference data between the other three hydrogen masers (KRISS codes: H26, H25, H24) and the master clock are also reported periodically. Another PSFS developed by KRISS, the cesium atomic fountain KRISS-F1 (hereafter F1), also provides a reference frequency while the KRISS-Yb1 (hereafter Yb1) is not operating.
When generating the ensemble time scale, an algorithm based on AT1 [9,10] or a Kalman filter is often used. The Kalman filter-based algorithm is used to generate Global Positioning System Time [11]. It has also been used to create UTC [12] or UTC(k) [13], which was created using AT1-based algorithms. The Kalman filter-based algorithm does not have a separate process of assigning weights unlike the AT1-based algorithm. However, according to [14], the Kalman filter-based algorithm is functionally identical with the AT1-based algorithm. Some review papers and commentary papers describing the Kalman filter-based algorithm have been published [15][16][17][18].
In this study, using the reduced Kalman filter (Kred) algorithm developed by Greenhall [19][20][21][22], we combined four hydrogen masers to create an ensemble time (TA), as named in [23]. We used one of the hydrogen masers (KRISS code: H26) as a reference clock for generating TA and measured the time differences of each member clock in the ensemble with respect to the reference clock before entering the result into the Kred algorithm as input data. Among the state vectors of the clocks outputted every minute from the Kred algorithm, the frequency offset (marked y (REF−TA)1min ) of the reference clock with respect to the TA was essential in this study. The Kred algorithm is characterized by high short-term frequency stability of TA compared to other Kalman filter-based time scale algorithms [19]. In order to improve the long-term frequency stability of the TA, it was necessary to steer the TA frequency to PSFS. For this, we estimated the frequency difference between TA and Yb1 (or F1) using frequency difference data measured between the reference clock and Yb1 (or F1) operating intermittently. This signal (marked y Steer ) made the TA frequency close to the frequency of Yb1 (or F1). Finally, by controlling the MPS with the dimensionless fre- , a frequency close to the PSFS could be physically generated from the MPS.
The structure of this paper is as follows. In section 2.1, we define the phase offset and frequency offset between two clocks, and introduce the notation used in this paper. In section 2.2, we briefly introduce the Kred algorithm. In section 2.3, we describe how to steer the TA frequency to the primary and secondary frequency standards (PSFS) using an MPS. In section 3, we introduce a system for realizing a realtime TA scale named KET. In section 4.1, we show how to predict the standard frequency at the present moment from PSFS data. In section 4.2, we show and discuss the changes in y (H26−TA)1min and y Steer due to the abnormal temperature change in the room in which hydrogen masers are placed. In section 4.3, we show the time difference of KET with respect to rapid UTC (UTCr) or UTC. In section 4.4, we compare the Allan deviations (ADEVs) of the phase states and the frequency states outputted from the Kred algorithm. Finally, we conclude this paper in section 5.

Phase offset and frequency offset between two clocks
The basic quantity needed to create an ensemble time is the time difference between two clocks. As detailed in [15], the time difference between two clocks is the difference in the clock readings of the two clocks and at the same time the phase difference of the two clocks. The phase offset x(k) with the dimension of time at epoch k represents the offset of a real clock (Clk) from the ideal one. However, in actual measurement, the phase offset of Clk means the time difference with respect to a reference clock REF as expressed in equation (1) x If x(k) is measured as a discrete time series of time interval τ , the normalized frequencyȳ k averaged over time τ can be expressed as equation (2) where the subscripts k and k−1 of x indicate the consecutive epochs whose time difference equals τ . By substituting the time difference between Clk and REF in equation (1) into equation (2), it becomes equation (3) If the nominal frequency of REF and Clk is the same as f REF , and the average of the difference frequency between REF and Clk over τ is ∆f, the normalized frequency can be written as equation (4) Therefore, while equation (3) expresses the normalized frequency as the time difference between two clocks, equation (4) expresses it as the frequency difference.
In this paper, the following notation is established to clarify the expression of equations: • In y representing the normalized frequency, the bar representing the average is omitted and the order of calculating the time difference between two clocks is indicated in the subscript of y as in equation (5) • Normalized frequency expressed as the frequency difference is also expressed in the subscript of y as in equation (6) y If the above notation is used, the frequency offsets between individual clocks (e.g. Clk, REF, TA) measured at the same epoch can be added or subtracted using subscripts as in equation (7) y Here, since all k indicating the measurement epoch is the same, it is omitted in equation (7).

Kred algorithm
The Kalman filter algorithm firstly defines a clock state vector, and then describes the change in the clock state over time by repeating the following three steps: (1) prediction of the state, (2) measurement of the time difference between clocks, and (3) update of the state estimate. In this cyclic iteration process, whenever the time difference measurement is updated, the estimated value of the state is determined from the predicted state.
However, as this process is repeated, the estimation error increases indefinitely. This comes from the fact that the number of measurements is less than the number of clocks. As the estimation error increases, the uncertainty of the estimation value increases, and there is also the problem that it takes up a lot of computational capacity in the computer. To solve this problem, K. Brown used a method that periodically reduces the size of the error covariance matrix without affecting the state vector [11].
Meanwhile, C. Greenhall discovered a simple way to solve this problem in 2003 and named it the Kred (reduced Kalman filter) scale [19]. This method reduces the magnitude of the error covariance by giving zeros to all matrix elements corresponding to the phase components (represented by x) in the error covariance matrix, which is called 'covariance x-reduction' [22]. The covariance x-reduction is performed after the state estimates are updated in the iterations of the algorithm. He mathematically proved that this does not affect the future estimates of frequency and frequency drift.
In this paper, we used a three-state model representing the state of a clock as phase (x), frequency (y), and frequency drift (d) [24][25][26]. If we denote the clock's state vector as S, it can be written as S = [x 1e , y 1e , d 1e ,..., x 4e , y 4e , d 4e ]. Since there are 3 states per clock, a total of 12 states are represented. Here, the number in the subscript is the number representing each clock and e means the ensemble time (TA). Therefore, y 1e means the frequency state of clock number 1 (i.e. the reference clock, H26) relative to the ensemble, which is expressed in our notation as y (H26−TA)1min . The subscript 1 min indicates that the state of the clock is estimated every minute. This signal is used to control the MPS after sign reversal.
In the Kred algorithm, it is necessary to determine the qfactors related to the noises of the clocks in the ensemble. We obtained them in the following way. First, the short-term frequency stabilities of individual clocks were obtained using the three-cornered hat method [27]. The mid-term stabilities were obtained by measuring the clocks' frequencies with respect to H28 which has better stability than other clocks and the longterm stabilities of clocks longer than τ ⩾ 1 d were measured with respect to PSFS. Combining these three frequency stability curves for each clock yields an ADEV graph of the clock over the entire averaging time (τ = 60 s to 1 × 10 6 s). Then, using the ADEV equation, which depends on the q-factors [24,25], we found the values of q-factors, that were the most similar to the ADEV graphs obtained previously.

Frequency correction using PSFS
An ensemble time (TA) created by combining multiple clocks with the Kred algorithm has higher frequency stability, particularly in short-term stability, than the individual clocks. If this TA is realized as a physical signal using an MPS, a hardware clock with high stability is secured although the noise of the MPS is added. However, the frequency may be different from that of PSFS. This frequency difference eventually creates a time deviation from UTC or UTCr.
Thus, the TA frequency was corrected using the optical clock Yb1 [28] or the cesium fountain F1 [29]. Then we compensated for the frequency difference between these PSFS and TAI, because we wanted to generate a TA frequency close to the TAI frequency. Finally, we produced KET by realizing the corrected TA frequency using an MPS. As a result, the longterm frequency stability of KET followed PSFS and TAI. This resulted in an ensemble clock close to UTC or UTCr. For this purpose, we followed the following process.
At KRISS, the average frequency of H28 with respect to PSFS is provided almost daily. However, in this study, we used the H26 as a reference clock when generating TA, so we needed to know the average frequency of H26 with respect to PSFS, which we obtained by equation (8) where y (H26−H28) is the average frequency of H26 relative to H28. We need the average frequency y (H26−PSFS) every 12 h. However, y (H28−PSFS) is measured on a daily basis, so the measurement data for the past 90 d (or 60 d, 30 d) of y (H28−PSFS) are linearly fitted to predict the optimal value every 12 h. Whenever the new measurement is updated, the process of linear fitting is repeated with the data for 90 d (or 60 d, 30 d) excluding the oldest data.
To obtain the frequency of PSFS with respect to TAI as in equation (10), we need to know y (H26−TAI) and y (H26−PSFS) . For We define y Steer as the dimensionless frequency to control the MPS to match the TA frequency to the TAI via the PSFS, and obtain y Steer by equation (11) using the results calculated in equations (9) and (10). Note that in the real experiments, only y (TA−PSFS) was initially used as the y Steer , and then y (PSFS−TAI) was included in the y Steer to see the effect of y (PSFS−TAI) Through this process, we can obtain a frequency close to TAI from the MPS output terminal. Figure 2 gives a schematic diagram of a system for the realization of a real-time KET whose frequency is corrected by PSFS using an MPS. Signals with all nominal frequencies of 5 MHz from four hydrogen masers enter the multi-channel measurement system (MMS, Timing Solutions Corp.). The MMS measures the time differences between the clocks using the dual-mixer time difference method. The hydrogen masers are numbered for simplicity, H26 is the clock number 1, and it is used as the reference clock in the ensemble. The measured time differences (z i1 , i = 2, 3, 4) of other clocks with respect to clock 1 are outputted from the MMS and enter the Kred algorithm. Among the state vector estimates outputted from the algorithm, the frequency estimate of H26 with respect to TA (i.e.,y (H26−TA)1min ) is used for controlling the MPS. An auxiliary output generator (AOG, model 110) manufactured by Symmetricom was used as the MPS. LABVIEW was used for communication with MMS and MPS, and MATLAB was used for the calculation and data analysis. The two programming languages were linked via data files.

Realization of KET
The systematic uncertainty of Yb1 was 1.7 × 10 −17 , and the ADEV was 3.2 × 10 −15 / √ τ , where τ is in seconds [28]. The frequency stability of F1 was reported in [29]. The frequency difference between PSFS and the reference clock was provided as a daily data file.
The KET signal was inputted to the MMS and the characteristics such as frequency stability were analyzed by measuring the time difference with other clocks. Furthermore, the long-term stability and accuracy of KET were analyzed using [UTCr-UTC(KRIS)] and/or [UTC-UTC(KRIS)] data released by BIPM.

Prediction of the standard frequency
We predicted the current standard frequency by linear fitting of the previously measured frequency data of H26 with respect to PSFS. We calculated three predicted values using measurement data for the past 90 d, 60 d, and 30 d. We also determined the current standard frequency as the average of these three predicted values because the slopes of the three fitting lines are different, and the 30 day fitting line represents the recent trend, while the 90 day fitting line represents the long-term trend. Figure 3 shows the error of the predicted value for the current moment (in this case, on modified Julian date (MJD) 59 944) using y (H26−F1) data of the past 90 d, 60 d, and 30 d.   stable among the four predicted values, so we chose the mean value as the current value.

Variations of y MPS and y Steer
The frequency y MPS that controls the MPS is obtained as −y (H26−TA)1min + y Steer as shown in figure 2. The steering frequency y Steer adjusts TA frequency to the reference frequency (PSFS or TAI) and is obtained by equation (11). In this study, only −y (TA−PSFS) was first applied to y Steer , and then −y (PSFS−TAI) was added to investigate the effect of the latter. Figure 4 shows the abnormal temperature variation of the clock room where the hydrogen masers are placed. Normally, it maintains a stable state at (21.8 ± 0.2) • C. However, between MJD 59 710 and 59 730, the temperature rose by up to 1 • C due to an abnormality in the air conditioning equipment. We investigated changes in the output values of the Kred algorithm and the y Steer according to the temperature variation between MJD 59 702 and MJD 59 731 because the phenomena that occurred during that period helped us understand how the Kred algorithm works.   +9.0 × 10 −15 . Since a new value is added every 12 h (i.e. one value is held for less than 12 h), there are flat parts in the graph as shown in figure 6. To ensure the frequency stability while the value of y Steer changes, we made the frequency slowly change by increasing or decreasing at the rate of 5 × 10 −18 /min. The overall shape of y Steer is similar to that of y (H26−TA)1min in figure 5. This means that the y Steer compensates for the frequency change due to −y (H26−TA)1min applied to the MPS in response to the temperature change while simultaneously adjusting TA to PSFS.   To determine why the amount of change of UTC (KRIS) and KET was different according to the temperature change, we compared the frequency differences of H28 and H26 with respect to Yb1 for about 60 d including the period in which the abnormal temperature change occurred. Figure 8 shows only the residuals after removing the linearly fitted values from the frequency changes of the two clocks. The linearly fitted slopes (i.e. frequency drifts) of H28 and H26 were −8.6 × 10 −21 s −1 and −1.7 × 10 −22 s −1 , respectively. In the case of H28, it can be seen that the frequency decreased from around MJD 59 720 when the temperature rise started. On the other hand, in the case of H26, a frequency step slightly larger than that of the other days occurred only on MJD 59 730. The standard deviations of the residuals of H28 and H26 for about 60 d were similar, 1.5 × 10 −15 and 1.6 × 10 −15 , respectively. However, for 60 d before MJD 59 720, the standard deviation of the residual of H28 was significantly reduced to less than onethird as 4.3 × 10 −16 , but that of H26 was almost unchanged. In the end, we inferred that the main reason for the disparity in the time difference of the two time scales (UTC(KRIS) and KET) against UTCr was that the responses of the clocks to the temperature change were different. That is, H28 was more sensitive than H26 in frequency sensitivity to the temperature change.

Time differences with respect to UTCr or UTC
BIPM announces the frequency difference (indicated by d) and its uncertainty (indicated by u) between TAI and PSFS reported by laboratories around the world through Circular T every month. Measurement data of Yb1 have also been submitted to BIPM since January 2020. From then until recently, the d values of Yb1 were between −6.5 × 10 −16 ∼ 8.7 × 10 −16 , and the u values were between 2.5 × 10 −16 ∼ 1.04 × 10 −15 . Sometimes PSFS fail to submit data to BIPM for a variety of reasons, resulting in unknown d-value for that PSFS. Therefore, we calculated y (PSFS−TAI) values corresponding to the d values by the method described in equation (10), and used them in this study.
In order to calculate y (Yb1−TAI) , we need to know the value of y (H26−TAI) as in equation (10). We predicted the current value of y (H26−TAI) by linear fitting of the [TAI-H26] data published by BIPM. The effect of y (Yb1−TAI) appears in the longterm in terms of frequency stability and in the time difference graph with respect to UTCr, so it was necessary to observe at least several weeks. Figure 9 shows the changes in [UTCr-KET] as results of the following two tests. UTC (KRIS) is also drawn for comparison. The first result is the graph from MJD 59 765 to MJD 59 800 in the figure. During this period the frequency drift estimate (d 1e ) of H26 outputted from the Kred algorithm was replaced by the frequency drift of H26 against Yb1 (i.e. d (H26−Yb1) ) in the algorithm. The original purpose of doing this was to improve the long-term stability of KET, but the actual result was that the time difference graph continued to increase smoothly, so from MJD 59 801, we did not use d (H26−Yb1) , but d 1e which was produced by the Kred algorithm.
The second test included the y (Yb1−TAI) obtained by equation (10) in y Steer as equation (11). Once or twice a week, we calculated the y (Yb1−TAI) and inputted the value into the program through the front panel of the LABVIEW program in order to add to the MPS control signal. The MJDs marked with triangles in the lower part of figure 9 represent the days when this value was applied. Contrary to expectations, the variation of the time difference did not improve significantly. We presume this was because of the following reasons. As explained in figure 8, in the normal condition, the standard deviation of the residual of H26 with respect to Yb1 was about 1.6 × 10 −15 however, y (Yb1−TAI) was as small as 10 −16 level, so its effect could not be seen in the [UTCr-KET] graph. In addition, the frequency difference between H28 and Yb1 was directly measured, but the frequency difference between H26 and Yb1 was calculated using equation (8). As a result, KET made based on H26 could not have a smaller deviation than UTC(KRIS) made based on H28. Since UTC(KRIS) was also steered to Yb1, it can be seen that the changing shape of the two graphs was quite similar after MJD 59 820. The variation of [UTCr-KET] was observed less than 5 ns for 142 d when only y (TA−Yb1) was used as y Steer .
An experiment was performed using F1 as the PSFS under the same condition as above, i.e. not including y (F1−TAI) in y Steer . Figure

Frequency stabilities of state estimates
We calculated the overlapping ADEV of the phase-state estimates (Sx) and the frequency-state estimates (Sy) which are outputs of the Kred algorithm. Among the four hydrogen masers in the clock ensemble, the ADEVs of H26 and H28 are shown in figure 11. One of the characteristics of the Kred algorithm is that white frequency noise is greatly reduced by the covariance x-reduction [19,21]. Thanks to this, the short-term ADEVs of Sy (H26−TA) and Sy (H28−TA) were very small in the 10 −18 area for average time (τ ) less than 10 3 s. In contrast, the ADEVs of Sx (H26−TA) and Sx (H28−TA) show values that are normally expected of hydrogen masers. The ADEV of H28 was superior to that of H26 in the short term. However, around τ = 1 d (86 400 s), H28 showed drift, and at the same time, the ADEV of Sx (H28−TA) became equal to that of Sy (H28−TA) for τ > 1 d. This means that the white frequency noise of H28 disappeared around τ = 1 d and the random walk frequency noise became  and Sy (H26−TA) were equal at τ = 1 × 10 6 s.
The frequency of the KET was measured against F1, and the ADEV is also shown in figure 11, where 1.72 × 10 −15 at τ = 1 d and gradually decreased as τ increased.

Conclusion
Four hydrogen masers were combined to generate an ensemble time (TA) every minute with the Kred algorithm, and the KRISS ensemble time KET was made by steering the TA frequency to the KRISS-PSFS, i.e. the KRISS-Yb1 optical frequency and/or the KRISS-F1 cesium fountain.
Comparing the changes of [UTCr-KET] and [UTCr-UTC(KRIS)] according to the abnormal temperature rise of the clock room, we found that the former was more stable. This was caused by the difference in the sensitivity of clock frequency to temperature change between H26, the reference clock of KET, and H28, the master clock of UTC(KRIS).
To generate a KET closer to UTC or UTCr, we included y (PSFS−TAI) in addition to y (TA−PSFS) in the steering signal y Steer that controls the MPS as in equation (11), but it did not result in an improvement because the standard deviation of the frequency fluctuation of H26 against Yb1 was about 1.6 × 10 −15 for 60 d, whereas the frequency difference between Yb1 and TAI was small at the level of 10 −16 , so the effect was not shown. In other words, as long as H26 is used as a reference clock, there is no need to compensate for the frequency difference between Yb1 and TAI.
In the experiment using F1 as PSFS, [UTCr-KET] was stable within 5 ns, and [UTC-KET] was stable within 2 ns for 75 d when only y (TA−F1) was used for y Steer . In comparison, [UTC-UTC(KRIS)] showed a change of about 7 ns during the same period.
The Allan deviation was calculated for the phase-state estimate and frequency-state estimate outputted from the Kred algorithm. The ADEV for the phase state was worse for H26 than H28 in the short term, but better for H26 for τ longer than 1 d. In the ADEV for the frequency state, the short-term stability was in the 10 −18 region as the white frequency noise disappeared by the covariance x-reduction, which is a feature of the Kred algorithm.
UTC(lab) must be generated without interruption under any circumstances to supply accurate time as required by NMIs as well as navigation systems such as a Global Navigation Satellite System or Regional Navigation Satellite System. Therefore, it is necessary to prepare and have various options for emergency situations [3,31]. The method proposed in this paper can be used as an option to generate UTC(lab) for such a case.