## 2 Methods
### Numerical Simulation Method  
To validate our approach and identify the optimal modulation parameters (\(\alpha\) and \(\beta\)) for minimizing the estimated uncertainty, we employ numerical simulations. The flow chart of this simulation is shown in Figure 1, which involves three main steps. In step 1, we set the values of \(T_{\text{P}}\), \(T_{0}\), \(N_{0}\), \(\gamma_{0}\), \(\delta_{\text{N}}\), \(\alpha\), and \(\beta\), and calculate \(S_{0}\), \(\sigma_{\text{a}}\), \(\sigma_{\text{b}}\), and \(\sigma_{\text{c}}\). The noise-induced frequency fluctuation is represented by the discrete normal random numbers with the standard deviation of \(\sigma_{\text{a}}\) for case 1, \(\sigma_{\text{b}}\) for case 2, and \(\sigma_{\text{c}}\) for case 3. In step 2, the clock comparison process between two clocks is executed with a simulated time of approximately 2.2 h for each case. The total simulated time is obtained by multiplying the total clock cycle number by \(3T_{0}\), where the reason of the use of \(3T_{0}\) is that three measurements are requiredin real-world scenarios to determine the values of \(\sigma_{a}\), \(\sigma_{b}\), and \(\sigma_{c}\). In this work, the Dick effect is cancelled by setting the clock laser noise to be zero [27]. The cancellation of the Dick effect can be realized in experiment by synchronous frequency comparison between two clocks [23; 28], or by using the in situ imaging technique to compare two regions of cold ensembles in a clock [18; 19; 22; 25]. Regarding the simulation aspect of the in-loop clock operation, the excitation fractions at half-height points are determined by adjusting the center frequency of the Rabi spectrum by \(\pm 0.4/T_{\text{p}}\). Subsequently, the corresponding generated discrete normal random numbers from step 1 are added to the determined excitation fraction, and the frequency corrections of each clock are calculated. These corrections are then utilized to update the center frequency of the Rabi spectrum. In step 3, upon completion of all clock cycles, the comparison stability is derived by computing the Allan variance of the frequency difference between the two clocks. It is worth noting that the Allan variance of a single clock should be divided by 2 [28]. Subsequently, the parameters of \(\sigma_{\text{QPN}}\), \(\sigma_{\text{Shot}}\), and \(\sigma_{\text{Det}}\) can be determined and the value of \(N_{0}\) is determined by \(T_{0}/4S_{0}^{2}\sigma_{\text{QPN}}^{2}\).  
In this work, we maintain \(T_{0}\) at 1 s and utilize Rabi detection with a fixed interrogation time of \(T_{\text{p}}=0.1\) s, which can be easily implemented in the experiment. It is important to note that the specific choice of \(T_{\text{p}}\) has minimal impact on the numerical results obtained. The collision between atoms influences the atomic-density-dependent \(S_{0}\)[13; 14]. However, in our simulation, where \(T_{\text{p}}\) equals 0.1 s and the maximum atom number is lower than 6000, the system operates in the weak interaction region. Consequently, the variation in \(S_{0}\) can be reasonably disregarded as the atom number changes [13; 14]. Nevertheless, if the collision effects become prominent in a regime characterized by high density and strong interactions, the value of \(S_{0}\) will be reduced due to density broadening. This reduction in \(S_{0}\) indicates a higher level of instability for the clock.  
Figure 1: The flow chart of the numerical simulation.