## 3 Results and Discussion  
To verify the accuracy of our numerical calculation code, we compare the stabilities obtained from numerical simulations at three sets of parameters to the theoretical results shown in Figure 2. The parameters used for this comparison are \(T_{0}=1\) s, \(N_{0}=500\), \(\gamma_{0}=1\), \(\delta_{\rm N}=3\)[23], \(\alpha=0.2\), \(\beta=0.1\). The excellent agreement observed in Figure 2 between the numerical and theoretical results confirms the correctness of our code and the validity of our method. We investigate the influence of \(\alpha\) and \(\beta\) on the estimated uncertainty of \(N_{0}\) by studying the standard deviation of 50 independent simulations at different combinations of \(\alpha\) and \(\beta\), as shown in Figure 3a for \(N_{0}=500\) and Figure 3b for \(N_{0}=2000\). The white regions in both figures indicate errors larger than 500 for Figure 3a and 2000 for Figure 3b, which occurs when \(\alpha\) or \(\beta\) approaches 1. In such cases, the modulation amplitudes of the parameters become close to zero, indicating that the different noises cannot be distinguished. Similar uncertainty distributions are observed for both cases. The smallest uncertainty is achieved at \(\alpha=5.71\) and \(\beta=0.1\) for \(N_{0}=500\), and \(\alpha=3.84\) and \(\beta=0.27\) for \(N_{0}=2000\). We also found that the maximum difference in estimated uncertainty within the ranges of \(\alpha=0.1\)-\(0.44\) and \(\beta=1.97\)-\(6.73\) is below 4% at the current total simulated time of 6.7 h.  
Figure 3: Numerical results of estimated uncertainties at different modulation parameters of \(\alpha\) and \(\beta\). (**a**) The standard deviation of 50 independent simulations at \(N_{0}=500\). (**b**) The case of \(N_{0}=2000\). Figure 2: Comparisons of numerical and theoretical results of stabilities. The Allan variances of frequency fluctuations at three-group parameters. The points represent the numerical results and dashed lines indicate the corresponding theoretical results. All error bars represent the \(1\sigma\) standard error. Therefore, we choose the parameters of \(\alpha=3.84\) and \(\beta=0.27\) to numerically simulate the measurements of \(N_{0}\) using our approach. Figure 4 demonstrates good agreement between numerical and theoretical results, as \(N_{0}\) exceeds approximately 200. However, it should be noted that larger uncertainties are observed when \(N_{0}\) is lower than 200, due to larger noise and reduced stability. The larger noise can correspond to stronger excitation fraction fluctuation, which may cause deviations in the half-height points of the Rabi spectrum more frequently. Although the frequency-sensitivity slope \(S_{0}\) is not constant for the Rabi spectrum, we use a constant \(S_{0}\) to infer \(N_{0}\), leading to deviations for small \(N_{0}\). We verify our inference by conducting a comparison of \(N_{0}\) simulated results using the Rabi and triangular spectra (shown in Figure 5), respectively. The values of \(S_{0}\) remain unaffected by the frequency detuning on both the left and right sides in the case of a triangular spectrum. As a result, the stability calculation utilizing Equation (1) will not be influenced by fluctuations in the excitation fraction. Figure 4 demonstrates good agreement between theory and numerical result as the triangular spectrum is used, which is the powerful evidence of our hypothesis. Nevertheless, to achieve the triangular spectrum is challenging, this phenomenon suggests that our method can effectively work when \(N_{0}\) exceeds about 200. Exploring the relationship between total time consumption and estimated uncertainty using our approach presents an intriguing avenue for further study, given their crucial role in experiments. Figure 6 shows the relationship between relative uncertainty and total simulated time. It is evident that as the time increases, the uncertainties decrease following a slope of \(-0.5\).