June 24, 2026
Rapid and robust laser-frequency auto-locking is essential for the field deployment of quantum communications, quantum computing, and precision-measurement technologies; however, achieving this remains a considerable challenge. Here, we propose and demonstrate an auto-locking scheme employing Bayesian optimization and discrete biorthogonal wavelet transformation. First, the reference is rapidly sought by making intelligent use of historical observations, eliminating the inherent blindness of the traditional parameter-scanning method. Second, the frequency reference is robustly identified by pinpointing transition signals with the discrete biorthogonal wavelet transformation and analyzing their immutable frequency differences and relative magnitudes, which are determined by the inherent atomic structure and remain resistant to environmental disturbances. This proposed approach achieves a fivefold acceleration in reference searching compared to conventional scanning methods in the case where the laser frequency drifts far away from the reference. Crucially, it achieves an identification accuracy of more than 99.5%, even under severe 50% laser-intensity fluctuations, \(9.95^\circ\) photodiode misalignment, and \(18^\circ\)C Rb cell temperature elevation. Finally, locking the laser frequency to the identified reference with a lead zirconate titanate-current double-servo loop narrows the linewidth to 20 kHz. We believe that this rapid, robust, and high-performance auto-locking technique will be pivotal towards the deployment of the next generation of practical quantum technologies in demanding field environments.
Laser-frequency locking plays a critical role in a variety of quantum technologies, including quantum computing [1], [2], quantum communication [3], [4], and quantum precision-measurement systems [5]–[8]. As these technologies progressively transition from the laboratory to field applications [9]–[11]—for example through deployment in space stations [12], [13] and sounding rockets [14]—laser-frequency auto-locking will become increasingly important. Although various auto-locking methods have been developed [14]–[24], their practical use is constrained by slow reference search speeds and limited robustness in reference identification. To enable broader field deployment, there is an urgent need for laser-frequency auto-locking methods that can quickly locate and reliably lock on to a reference, thereby ensuring the immediate readiness of quantum systems and improving their adaptability to changing environmental conditions.
First, rapidly searching for the reference is an essential capability for field-deployable quantum technologies, particularly in applications such as quantum inertial navigation, where continuous sensor operation is preferable. Compared with conventional gradual-scanning methods that use a fixed scanning interval and fixed step size, the mean-shift algorithm, which can dynamically calculate the initial scanning value [25], and the coarse-to-fine search strategy, which applies a large step size followed by a refined small step size [21], can improve search efficiency. Despite these advances, such approaches have not fundamentally changed the gradual-scanning nature of the process, and this leads to inherent limitations. In particular, when the search range is large, the search efficiency is still low.
Second, in addition to searching of the reference, accurately and robustly identifying that reference constitutes another indispensable requirement for quantum systems employed in field applications. Threshold comparison methods [15]–[22] and cross-correlation calculation algorithms [14], [23], [24] can be used to accurately identify the reference within specific contexts, but they may misidentify the reference when the laser intensity is affected by environmental variations such as temperature fluctuations [26]. Gaussian continuous-wavelet-transform-based noise reduction [27] can improve robustness, but it does so at the expense of altering the signal. In a similar trade-off, a deep-neural-network identifier [28] has been shown to provide robust reference identification, but it requires 80 s of processing time. To our best knowledge, no existing method can both rapidly search for and robustly identify the reference signal simultaneously.
In this paper, we present a rapid and highly robust laser-frequency auto-locking scheme using Bayesian-optimization and discrete-wavelet-transformation algorithms. First, a Gaussian-process-based Bayesian-optimization method is employed to speed up the reference search. By purposefully selecting the next test value based on historical observations of the laser-frequency control parameters, this approach avoids the blindness inherent in conventional gradual-parameter-scanning methods. A discrete biorthogonal wavelet transformation is then used to accurately locate atomic transition signals within the collected spectra. By assessing the frequency differences and relative magnitudes of these atomic transition signals, the robustness of reference identification is strengthened. Our experimental results show that the Gaussian-process-based Bayesian-optimization algorithm achieves a fivefold increase in reference searching speed compared to the conventional gradual-scanning method, especially when the laser frequency drifts far away from the reference, and the discrete biorthogonal wavelet-transformation method achieves a 99.5% reference-identification rate, even under 50% laser-intensity variation, \(9.95^\circ\) photodiode misalignment, and \(18^\circ\)C Rb cell temperature elevation. Finally, the laser frequency is auto-locked by rapidly searching for and accurately identifying the atomic spectrum, reducing the laser linewidth to 20 kHz using a combined current and lead-zirconate-titanate (PZT) double-servo feedback control loop. This rapid and robust laser-frequency auto-locking method with narrow-linewidth performance is thus beneficial for advancing quantum technologies towards field applications.
Figure 1 shows the framework of an advanced scheme for rapidly searching for and accurately identifying reference signals. The reference signal is rapidly sought using Gaussian-process-based Bayesian optimization. As illustrated in Fig. 1 (a), the observed laser-frequency control parameters are used to calculate the covariance with unobserved parameters. The predicted signal values and uncertainties of the unobserved parameters are provided by the multivariate Gaussian posterior. With this posterior, an expected-improvement (EI) sampling function is used to estimate the next parameter at which the maximum or minimum signal value, corresponding to the atomic transition, may occur. Let \(x\) denote the laser-frequency control parameter and its set be \(\boldsymbol{\chi}\). The functional relationship between the reference signal and the parameter \(x\) is \(g(x)\), which is unknown because the laser frequency may drift due to variations in environmental conditions. Gaussian processes offer a good balance between modeling accuracy and computational complexity; thus, the Gaussian process \(f(x)\) is used as an alternative model for \(g(x)\). Let \(n\) observed parameter values \(x_{i}\in \boldsymbol{\chi}\) form a vector \(\boldsymbol{X} = \{x_i\}_{i = 1}^n\). The corresponding observed signal values \({f_i}{\buildrel \Delta \over =}f(x_i)\) form the vector \(\boldsymbol{f} = \{f_i\}_{i = 1}^n\). The covariance matrix between the two parameter vectors \(\boldsymbol{X}^1\) and \(\boldsymbol{X}^2\) is defined as \[\label{eq:1} \boldsymbol{K}\left( {\boldsymbol{X^1},\boldsymbol{X^2}} \right) = \left[ {\begin{array}{*{20}{c}} {k\left( {x_1^1,x_1^2} \right)}&{k\left( {x_1^1,x_2^2} \right)}& \cdots &{k\left( {x_1^1,x_n^2} \right)}\\ \vdots & \vdots & \ddots & \vdots \\ {k\left( {x_m^1,x_1^2} \right)}&{k\left( {x_m^1,x_2^2} \right)}& \ldots &{k\left( {x_m^1,x_n^2} \right)} \end{array}} \right],\tag{1}\] where \(m\) and \(n\) are the numbers of elements in the vectors \(\boldsymbol{X^1}\) and \(\boldsymbol{X^2}\), and \(k\left(x_p,x_q\right)=\exp\left(-\left|x_p-x_q\right|^2/2\right)\). Let \(\boldsymbol{X_*}\) denote a vector of unobserved parameters and \(\boldsymbol{f_*}\) the corresponding unobserved function values. Since \(g(x)\) is modeled by the Gaussian process \(f(x)\), the random variables \(\boldsymbol{f}\) and \(\boldsymbol{f_*}\) follow a multivariate Gaussian distribution [29], [30]: \[\label{eq:2} \left[ {\begin{array}{*{20}{c}} \boldsymbol{f}\\ \boldsymbol{f_*} \end{array}} \right] \; = {\cal N}\left( {\begin{array}{*{20}{c}} \boldsymbol{0}, & {\left[ {\begin{array}{*{20}{c}} \boldsymbol{K}(\boldsymbol{X},\boldsymbol{X}) & \boldsymbol{K}(\boldsymbol{X},\boldsymbol{X_*})\\ \boldsymbol{K}(\boldsymbol{X_*},\boldsymbol{X}) & \boldsymbol{K}(\boldsymbol{X_*},\boldsymbol{X_*}) \end{array}} \right]} \end{array}} \right).\tag{2}\] According to the conditional probability formula of multivariate Gaussian functions, the posterior probability density function of the Gaussian process is given by \[\label{eq:3} \begin{align} \boldsymbol{f_*} &\mid \boldsymbol{X_*}, \boldsymbol{X}, \boldsymbol{f} \sim \mathcal{N}\Big( \boldsymbol{K}(\boldsymbol{X_*},\boldsymbol{X}) \, \boldsymbol{K}(\boldsymbol{X},\boldsymbol{X})^{-1} \, \boldsymbol{f}, \\ &\boldsymbol{K}(\boldsymbol{X_*},\boldsymbol{X_*}) - \boldsymbol{K}(\boldsymbol{X_*},\boldsymbol{X}) \, \boldsymbol{K}(\boldsymbol{X},\boldsymbol{X})^{-1} \, \boldsymbol{K}(\boldsymbol{X},\boldsymbol{X_*}) \Big) \end{align}\tag{3}\]
To find the maximum or minimum signal value, the next test value \(x_{n+1}\) is purposefully selected by the EI sampling function based on the mean and variance of \(\boldsymbol{f_*}|\boldsymbol{X_*},\boldsymbol{X},\boldsymbol{f}\) and the currently observed maximum or minimum value of \(f(x)\). The following formula shows the maximum signal value search; the minimization case is analogous: \[\label{eq:4} x_{n + 1} = \mathop{\text{arg}\max}\limits_{x \in \chi} E{\left[ {\left( f(x) - \mathop{\max}\limits_{i = 1,\ldots,n} f(x_i) \right) \left| \boldsymbol{f} \right.} \right]^+ },\tag{4}\] where \([\cdot]^+ = \max(\cdot,0)\). Based on Eqs. (1 )\(\--\)(4 ), the next laser-frequency control parameter can be purposefully selected using observed parameters and their corresponding signal values, rather than through traditional step-by-step scanning. This significantly improves the efficiency of the search, as will be experimentally validated in Section 3.
Accurately identifying a desired reference signal, such as an atomic transition frequency or the characteristic frequency of an optical cavity, is indispensable in laser-frequency auto-locking. As shown in Fig. 1 (b), in our approach, the reference signal is accurately identified using discrete wavelet transformations. Because the frequency differences and relative amplitudes of the atomic transition signals within the spectrum are determined by the atomic energy-level structure and transition probabilities, the reference signal can be robustly identified, even in the presence of environmental variations. When the laser frequency coincides with an atomic transition frequency, the collected signal will change abruptly, indicating a specific atomic transition. The discrete biorthogonal wavelet transformation can accurately locate these abrupt changes without artificial shifting, owing to the perfect symmetry of its filter coefficients. Mathematically, let \(\widetilde{\phi}_{j,k}(t) = 2^{j/2}\widetilde{\phi}(2^j t - k)\) denote the biorthogonal scaling function at scale \(j\) and position \(k\). Here, a B-spline function of degree \(n\), denoted as \(B_n(t)\), is chosen as the fundamental scaling function \(\widetilde{\phi}(t)\) [31]. The real-time collected signal \(f(t)\) can then be expressed as an expansion in this scaling space: \[f(t) = \sum_{k \in \mathbb{Z}} \widetilde{a}_{j,k} \widetilde{\phi}_{j,k}(t),\] where \(\widetilde{a}_{j,k}\) are the scaling coefficients, and \(\mathbb{Z}\) is the set of integers. For a chosen degree \(n\), the discrete spline scaling filter \(\widetilde{h}\) can be derived analytically using the following formula: \[\widetilde{h}_k = \frac{\sqrt{2}}{2^{n+1}} \begin{cases} \binom{n+1}{\frac{n+1}{2} - k}, & k = -\frac{n+1}{2}, \dots, \frac{n+1}{2} \quad \text{for odd } n \\[2ex] \binom{n+1}{\frac{n}{2} + k}, & k = -\frac{n}{2}, \dots, \frac{n}{2} + 1 \quad \text{for even } n \end{cases}\] with all other coefficients \(\widetilde{h}_k = 0\), where \(\binom{n}{k}\) denotes the standard binomial coefficient. The symmetric dual scaling filter \(h\), which pairs with the spline filter \(\widetilde{h}\), can be subsequently constructed to satisfy the biorthogonality symbol condition: \[\widetilde{H}(\omega)\overline{H(\omega)} + \widetilde{H}(\omega + \pi)\overline{H(\omega + \pi)} = 1\] for all \(\omega \in \mathbb{R}\), where \(H(\omega) = \frac{1}{\sqrt{2}}\sum_{k \in \mathbb{Z}} h_k e^{-ik\omega}\), \(\widetilde{H}(\omega) = \frac{1}{\sqrt{2}}\sum_{k \in \mathbb{Z}} \widetilde{h}_k e^{-ik\omega}\), and \(\overline{H\left(\cdot\right)}\) denotes the complex conjugate of \({H\left(\cdot\right)}\). To accurately locate atomic transition signals within the collected spectrum, a discrete wavelet decomposition is performed on \(f(t)\) using the derived dual filters \(h\) and \(g\): \[\label{waveletDecomposition} \begin{align} f(t) =& \sum_{k \in \mathbb{Z}} \left( \sum_l h_{l - 2k} \widetilde{a}_{j,l} \right) \widetilde{\phi}_{j-1,k}(t) \\ &+ \sum_{k \in \mathbb{Z}} \left( \sum_p g_{p - 2k} \widetilde{a}_{j,p} \right) \widetilde{\psi}_{j-1,k}(t). \end{align}\tag{5}\] Here, \(\widetilde{\psi}_{j-1,k}(t) = 2^{(j-1)/2}\widetilde{\psi}(2^{j-1}t-k)\) represents the wavelet function, defined by \(\widetilde{\psi}(t) = \sqrt{2}\sum_{k \in \mathbb{Z}} \widetilde{g}_k \widetilde{\phi}(2t-k)\), and \(g_k = (-1)^k\widetilde{h}_{1-k}\) serves as the dual wavelet filter. Iterative wavelet decomposition can be applied to the first term, i.e., the low-frequency approximation, on the right-hand side of Eq. 5 , resulting in a multiscale representation: \[f(t) = \sum_{m \in \mathbb{Z}} b_m \widetilde{\phi}_{r,m}(t) + \sum_{n = r}^{j - 1} \sum_{i \in \mathbb{Z}} c_{n,i} \widetilde{\psi}_{n,i}(t),\] where \(r < j\) is the final coarse resolution level, and \(b_m\) and \(c_{n,i}\) represent the resulting scale and detail coefficients, respectively. When an atomic transition occurs, the values in the detail function \(\sum\nolimits_{i \in \mathbb{Z}} {{c_{n,i}}{\widetilde{\psi} _{n,i}}(t)}\) increase. The \(N\) largest values in the detail function are selected, and the density-based spatial clustering of applications with noise (DBSCAN) method [32], which does not require predefined numbers of clusters, clusters the \(x\) coordinates of these points. For each cluster, zero-crossing or local-maximum points of the signal are sought. The frequency differences of the atomic transition signals are calculated based on the \(x\) coordinates of these points. By also comparing the maximum and minimum values within the signal range corresponding to the clusters in the detail function, the desired spectrum and reference are robustly identified, as will be experimentally validated in the next section.
Algorithm 1 Rapid and robust laser-frequency auto-locking
\(S(v)\), search space \(\mathcal{V}\), threshold \(\epsilon\), max iters \(n\). Laser locked to the target transition frequency.
Initialize observation set \(\mathcal{O} \gets \emptyset\); Randomly sample initial parameter \(v_{next} \in \mathcal{V}\);
Set parameter to \(v_{next}\) and acquire signal \(S(v_{next})\);
Wavelet Decomposition: Decompose \(S(v_{next})\) using dual filters \((h,g)\); Compute detail functions at target scales;
Feature Extraction: Extract local maxima from absolute detail functions; Cluster maxima coordinates using DBSCAN; Locate candidate transition frequencies per cluster; Extract signal extrema (maxima/minima) per cluster;
Reference Identification: Compute transition frequency differences \(\Delta f\); Compute relative amplitudes \(\Delta A\) of the extrema; Verify \(\Delta f\) and \(\Delta A\) against atomic parameters;
Engage PID lock at target transition frequency; break;
\(\mathcal{O} \gets \mathcal{O} \cup \{ (v_{next}, \max|S(v_{next})|) \}\); Update Gaussian Process (GP) model using \(\mathcal{O}\); Select next \(v_{next} \in \mathcal{V}\) by maximizing the EI function;
To demonstrate the effectiveness of the proposed scheme, a laser-frequency auto-locking experimental apparatus based on the modulation transfer spectrum (MTS) was set up as shown in Fig. 2. In this apparatus, the output of the laser (an external cavity diode laser) is split into three beams after passing through two polarizing beam splitters (PBSs). One beam is directed into a wavelength meter to judge whether the laser frequency is locked to the desired frequency reference. Another beam serves as the input to the subsequent optical path, and the third beam, which is reflected by a mirror, is again split into the probe and pump beams after passing through a half-wave plate and a PBS. The polarization of the pump beam is controlled by two half-wave plates located next to the electro-optic modulator (EOM). After phase modulation by the EOM, the pump beam propagates in a direction opposite to and overlapping with the probe beam in a rubidium cell. The beat signal of the probe beam is detected by a fast photodiode. After amplification, this signal is mixed with the modulation signal of the pump beam. The output of the mixer is filtered by a low-pass filter, and the MTS error signal is finally obtained, as shown in Fig. 3 (b), with the corresponding saturated absorption spectrum (SAS) shown in Fig. 3 (a). Figures 3 (c)–3 (e) show enlarged sections of the transition signals marked in Fig. 3 (b). The MTS spectrum is composed of \(F=2 \to F^\prime=3\), \(F=2 \to F^\prime=\text{CO} 2,3\), and \(F=2 \to F^\prime=\text{CO} 1,3\) transition signals of \(^{87}\)Rb atoms in Fig. 3 (c); \(F=3 \to F^\prime=4\), \(F=3 \to F^\prime=\text{CO} 3,4\), and \(F=3 \to F^\prime=\text{CO} 2,4\) transition signals of \(^{85}\)Rb atoms in Fig. 3 (d); and \(F=2 \to F^\prime=3\), \(F=2 \to F^\prime=\text{CO} 2,3\), and \(F=2 \to F^\prime=\text{CO} 1,3\) transition signals of \(^{85}\)Rb atoms in Fig. 3 (e). Here, \(\text{CO}\) indicates crossover resonance. Note that the spectrum of transitions from the ground state of \(^{87}\)Rb \(F=1\) to the excited states is not shown due to its small signal amplitude in the current experimental apparatus. Figure 3 (c) is taken as the desired spectrum, and the transition frequency of \(F=2 \to F^\prime=3\) is used as the frequency reference in the subsequent experiments. The reference-identification algorithms should be able to distinguish the desired spectrum from the others, especially when their shapes are very similar.
The appropriateness of the algorithm parameter settings determines the performance of the proposed reference-identification method, and the implementation details of the algorithms are presented in Algorithm 1. To select the best parameters, the incident laser power of the MTS apparatus was first adjusted to 17.6 mW to ensure that the MTS error signal exhibited a good signal-to-noise ratio (SNR). A total of 1000 MTS error signals, as shown in Fig. 3 (c), were collected and then used to optimize the algorithm parameters to achieve the highest desired spectrum-identification rate. Specifically, a frequency window of 310 MHz was first set to encompass the three atomic transition signals of interest. Additionally, to facilitate a fast search for the desired spectrum, spectrum identification was performed only when the maximum absolute value of the signal within the frequency window exceeded 0.2 V. Selecting the degree of B-spline function in the biorthogonal wavelet transformation [31] involves a critical mathematical trade-off: higher-degree splines provide more vanishing moments to suppress noise, but their wider compact support smears the narrow transition signal.
| Filter pair | S1a | S2b | S3c |
|---|---|---|---|
| Linear (5,3) | 100% | 98.7% | 50.1% |
| Quadratic (8,4) | 99.9% | 99.7% | 53.6% |
| Cubic (7,5) | 86.5% | 84.0% | 38.5% |
a Normal alignment condition.
b PD misaligned by \(9.95^\circ\).
c PD misaligned by \(11.54^\circ\)
.
Table 1 demonstrates this trade-off using photodiode misaglignments of \(0^\circ\) (S1), \(9.95^\circ\) (S2), and \(11.54^\circ\) (S3) as low-SNR tests (which will be further detailed later in this section). The shorter (5,3) filter lacks sufficient vanishing moments to suppress severe detector noise, causing its accuracy to degrade in the S2 and S3 conditions. Conversely, the wider (7,5) filter over-smooths the signal, drastically reducing accuracy even under the clean S1 condition. Because the quadratic (8,4) filter perfectly balances noise suppression and sharp feature localization to yield the highest overall accuracy, the (8,4) biorthogonal filter pair was selected to gradually decompose the collected signals six times. As shown in Fig. 4 (b) by black diamonds, 15 maximum-value points were selected from the detail function, which was obtained by wavelet decomposition and is indicated by the red line. The DBSCAN clustering algorithm was used to cluster the \(x\) coordinates of these points. The clustering parameters, including the epsilon neighborhood and minimum number of neighbors required for a core point, were set to 19.4 MHz and 2, respectively, allowing the selected maximum-value points corresponding to the three atomic transition signals to be clustered into three categories. The average of each cluster was then calculated. Using this average as the center, the zero-crossing coordinates marked by the asterisks in Fig. 4 (a), as well as the maxima and minima of the MTS error signals within \(\pm 19.4\) MHz, were recorded, and the frequency differences between each pair of zero-crossing points were calculated. To identify the desired spectrum, the ranges of frequency differences between the \(F=2 \to F^\prime=3\) and \(F=2 \to F^\prime=\text{CO} 2,3\), \(F=2 \to F^\prime=\text{CO} 2,3\) and \(F=2 \to F^\prime=\text{CO} 1,3\), and \(F=2 \to F^\prime=3\) and \(F=2 \to F^\prime=\text{CO} 1,3\) transition signals were set to 113–152, 66–105, and 206–245 MHz, respectively. The calculated frequency differences between the zero-crossing points of the MTS error signal were then checked against these specified ranges. Furthermore, whether the recorded maxima corresponding to the \(F=2 \to F^\prime=3\), \(F=2 \to F^\prime=\text{CO} 2,3\), and \(F=2 \to F^\prime=\text{CO} 1,3\) transition signals were sequentially increasing, and whether the minima were sequentially decreasing, were verified to identify the desired spectrum.
Rapid laser-frequency locking was demonstrated using the Gaussian-process-based Bayesian-optimization and discrete-wavelet-transformation algorithms. The laser-frequency range was scanned by the PZT from 384.2252 to 384.2285 THz, ensuring that no mode hopping occurred within this range. This 3.3 GHz range was discretized into 65 samples with a PZT step of 0.1 V, simulating the challenging field scenarios, such as severe environmental disturbances or long-term operation, which require a wide-range search to recover the reference signal when the laser frequency drifts far away from the reference point. Two spectrum-searching methods were tested 100 times, with the frequency reference set to 384.227995 THz, corresponding to the \(F = 2 \to F' = 3\) transition frequency. The desired spectrum was identified using the wavelet-transformation-based method. Figure 5 shows the difference between the laser frequency read from the wavelength meter and the reference frequency set according to the experimental requirement after the laser frequency was locked. Both methods successfully found the desired spectrum within a 3.3-GHz frequency range in all 100 tests. The red squares show the frequency differences obtained using the Gaussian-process-based Bayesian-optimization method, and the black dots show the frequency differences obtained by the gradual-scanning method. There were seven instances in which the locked laser frequency deviated significantly from the reference. These discrepancies can be attributed to the communication rate of the DigiLock 110 because it took 100 ms to enable the proportional-integral-derivative (PID) controller to lock the laser frequency after defining the PID setpoint. During this interval, the MTS signal may have changed, preventing the PID controller from immediately locking on to the zero-crossing point of the MTS signal. It can be seen from Fig. 5 that the Bayesian-optimization method has superior performance when compared to the gradual-parameter-scanning approach, achieving better results in 99 out of 100 experimental trials. This improved performance results from our method’s ability to make full use of historical observation data within the Gaussian-process-based Bayesian-optimization framework, enabling intelligent prediction of the most promising parameters for reference identification. In contrast, the gradual-scanning method relies on brute-force parameter exploration from initial to target values. The experimental results indicate that the average number of PZT parameter-searching iterations was 5.52 for Bayesian optimization, compared to 44.76 for gradual scanning, quantitatively validating the much higher efficiency of our proposed approach. The average search times for the two methods were 7.08 and 36.15 s, respectively, indicating a fivefold improvement in search efficiency compared to the gradual-scanning method. To further compare the searching efficiency of our method with that of the gradual-scanning method using the optimized PZT step, a 2.6 GHz mode-hop-free range was discretized into 11 samples with a PZT step of 0.5 V for accelerating the searching speed. This step can also safely capture the three atomic transitions without skipping the desired spectrum. The gradual-scanning method under the optimized condition required an average of only 9.12 iterations (8.7 s). Remarkably, the Bayesian optimization method still demonstrated superior efficiency, requiring an average of only 5.02 iterations (3.8 s). This confirms that the Bayesian-optimization method maintains an algorithmic advantage in search efficiency under the conditions of both a fine grid and a coarse grid. While generally being rapid, the Bayesian-optimization algorithm required an extended search time in 6 out of the 100 trials in Fig. 5. These specific instances arise when the algorithm commands a large, instantaneous voltage step, inducing transient mechanical ringing in the PZT actuator. This ringing distorts the collected MTS signal and leads to a false-negative identification. Because the Bayesian-optimization acquisition function is tuned to favor exploitation, the algorithm temporarily searches other suboptimal regions before eventually returning to explore the true reference, thereby increasing the total search time for these specific trials.
It is also worth noting that in Fig. 5, the search time for the gradual-scanning method exhibits an overall decreasing trend across the 100 experimental trials, because the voltage coordinates corresponding to the atomic transitions drift closely to the start of the search range. The proposed auto-locking scheme consistently and accurately identified the desired spectrum. This result demonstrates advantage of the discrete-wavelet identifier, i.e., by evaluating the immutable relative frequency differences and spectral morphology of the transitions rather than relying on absolute voltage lock-points. The laser frequency was locked after the reference was autonomously sought and identified. As the black data points and black fit line in Fig. 6 show, the full width at half height (FWHH) of the beatnote signal is 840 kHz with only the PZT feedback loop, measured by the heterodyne beatnote method using two identical, independently locked diode lasers. Because this 840 kHz linewidth is largely dominated by the free-running laser frequency noise, it is difficult to meet the experimental demand of tens of kilohertz in the sensitivity-limit measurement [33], [34]. Therefore, the servo loop bandwidth must be expanded into the megahertz regime. To further reduce the linewidth, the bandwidth of the feedback loop was expanded by adding a laser-current feedback loop with a bandwidth of 20 MHz. With this double feedback loop, the high-frequency noise was suppressed, and the FWHH was significantly narrowed from 840 to 28 kHz, as the red data points and red fit line in Fig. 6 show, implying that the single laser linewidth is below 20 kHz.
The robustness of the auto-locking scheme was tested by changing the laser power. It was first increased from 17.6 to 19.4 mW, before being decreased to 8.8 mW. Correspondingly, the maximum voltage of the desired spectrum first increased by 1.14% and then decreased by 25.79%. Under each laser-power condition, 1000 MTS error signals were collected for testing. During the testing process, if the desired spectrum was not identified at the preset scale factor, an attempt was made to identify the reference signal at the preceding or subsequent scale factors in the wavelet decomposition. If the desired spectrum still could not be successfully identified, the zero-crossing points of the spectrum obtained at the current and previous scale factors were combined to determine whether the current signal corresponded to the desired spectrum. The test results indicate that the identification rate of the desired spectrum varies with laser power. Nonetheless, the proposed method achieved a minimum accuracy of 99.8%, with the highest identification rate reaching 99.9%, validating the robustness of our approach. There are three possible reasons for false-negative spectrum identification. First, the selected 15 maximum-value points of the detail function either did not contain—or only contained one or two—maximum-value points corresponding to the \(F=2 \to F^\prime=\text{CO} 1,3\) transition signal. Even when two maximum points were present, the frequency difference between them exceeded the preset frequency range and they were thus not clustered together. Second, the maximum-value points corresponding to the \(F=2 \to F^\prime=\text{CO} 1,3\) and \(F=2 \to F^\prime=\text{CO} 2,3\) transition signals were accidentally clustered together due to spectrum noise. Third, the clustering of maximum values of the detail function corresponding to the \(F=2 \to F^\prime=\text{CO} 2,3\) transition signal is also affected by spectrum noise, shifting the search range for zero-crossing points of the transition signal. To further verify the robustness of the proposed method, the algorithm is evaluated by adjusting the PD alignment and the Rb cell temperature, as shown in Fig. 7. First, the PD was deliberately misaligned by angles of \(2.83^\circ\), \(5.67^\circ\), \(9.95^\circ\), and \(11.54^\circ\). Under each condition, 1,000 target MTS signals were recorded and evaluated, as shown in Fig. 8. Although the shape and amplitude of the spectrum were changed, the high identification accuracies of \(99.7\%,\;99.5\%,\;99.7\%\) for the first three angles were maintained, demonstrating high robustness up to a deviation of \(10^\circ\). At the misalignment of \(11.54^\circ\), the accuracy dropped to \(53.6\%\), because severe clipping of the probe beam on the PD active area limits the light exposure, causing the weakest spectral feature, i.e., \(F=2 \to F^\prime=\text{CO} 1,3\) transition, to be buried in the detector noise. Similarly, the robustness can also be verified by changing the spectral shape and amplitude through adjusting the misalignment between the pump and probe beams. Second, the Rb vapor cell was heated from \(22^{\circ}\)C to \(40^{\circ}\)C. As shown in Fig. 8, for the 1000 signals collected under the heated Rb cell condition, the identification rate reached \(100\%\). This excellent performance is attributed to the higher atomic vapor density caused by the increased temperature, which naturally enhances the SNR of the MTS error signal. To evaluate the false-positive rate of the proposed method, 1000 frames of data for each non-desired spectrum were collected and analyzed, as illustrated in Figs. 3 (d) and 3 (e). These results show that there are no false-positive cases. Overall, these test results demonstrate that our method exhibits a highly successful identification rate with the strong robustness against variations of laser power, Rb cell temperature, and PD alignment.
In this paper, we have proposed and demonstrated a rapid and robust laser-frequency auto-locking method with a narrowed linewidth. A Gaussian-process–Bayesian-optimization algorithm was developed to rapidly search for the reference signal by making use of the observed data to purposefully select the next parameter for finding the desired reference, rather than blindly scanning the parameters step by step. A discrete-wavelet-transformation algorithm was used to locate the transition signals in the collected signals in real time. By analyzing the frequency differences and relative magnitudes of the transition signals, which are determined by the atomic energy-level structure and the relative atomic transition probabilities, the reference signal was robustly identified. After the reference signal was rapidly located and robustly identified, the laser frequency was auto-locked using the PZT-current dual-loop servo. The experimental results showed that our method achieves a fivefold improvement in search efficiency compared to the gradual-scanning method when the laser frequency drifts far away from the reference, along with 99.5% identification accuracy, even under 50% laser power fluctuations, \(9.95^\circ\) photodiode misalignment, and \(18^\circ\)C Rb cell temperature elevation, validating its high robustness. Using the heterodyne beatnote method, the laser linewidth was found to be 20 kHz. Although the approach was demonstrated here for MTS lines, its universality enables straightforward extension to other spectral references (e.g., SAS lines). In the future, the system will be tested through vibration, shock, and variations in pressure, humidity, and temperature for field applications. Furthermore, we will explore integrating machine-learning-based noise profiling to autonomously diagnose mode hops.
We acknowledge the financial support provided by the National Key Research and Development Program of China under Grant No. 2025YFF0515200; the Quantum Science and Technology-National Science and Technology Major Project of China under Grant No. 2021ZD0300604; the National Natural Science Foundation of China under Grants No. U25D9005, No. 12104466, and No. 12504571; the Postdoctoral Innovation Research Posts in Hubei Province of China under Grant No. R20R0004; the Natural Science Foundation of Wuhan under Grant No. 2025040601030117; the Shenzhen Fundamental Research (Key Program) under Grant No. JCYJ20241202124931042; the Key Science and Technology Project of the Shenzhen Science and Technology Innovation Commission (SSTIC) under Grant No. KJZD20230923115505011; the cluster special project of Shenzhen Institutes of Advanced Technology, Chinese Academy of Sciences under Grant JQ0101-2025/02.
The authors declare no conflicts of interest.
All relevant data are available from the authors upon request.