Performing noise characterizations of lasers and optical frequency combs on sampled data offers numerous advantages compared to analog measurement techniques. One of the main advantages is that the measurement setup is greatly simplified. Only a balanced detector followed by an analog-to-digital converter is needed, allowing all the complexity to be moved to the digital domain. Secondly, near-optimal phase estimators are efficiently implementable, providing accurate phase noise estimation in the presence of measurement noise. Finally, joint processing of multiple comb lines is feasible, enabling the computation of the phase noise correlation matrix, which includes all information about the phase noise of the optical frequency comb. This tutorial introduces a framework based on digital signal processing for phase noise characterization of lasers and optical frequency combs. The framework is based on the extended Kalman filter (EKF) and automatic differentiation. The EKF is a near-optimal estimator of the optical phase in the presence of measurement noise, making it very suitable for phase noise measurements. Automatic differentiation is key to efficiently optimizing many parameters entering the EKF framework. More specifically, the combination of EKF and automatic differentiation enables the efficient optimization of phase noise measurement for optical frequency combs with arbitrarily complex noise dynamics that may include many free parameters. We show the framework’s efficacy through simulations and experimental data, showcasing its application across various comb types and in dual-comb measurements, highlighting its accuracy and versatility. Finally, we discuss its capability for digital phase noise compensation, which is highly relevant to free-running dual-comb spectroscopy applications.

## I. INTRODUCTION

Optical frequency combs are indispensable tools in many scientific and engineering disciplines.^{1,2} Over the past two decades, the field of optical frequency combs has witnessed significant advancements, leading to the development of various types of frequency combs with distinct characteristics and applications.^{3–6} The development of broad and stable combs has led to their integration in various fields, notably playing vital roles in the establishment of ultra-stable frequency references for metrology applications, high-capacity coherent telecommunication, and high-precision spectroscopy.^{7–10}

While extensive research efforts have been made to realize various types of optical frequency combs, the development of novel noise characterization techniques has not followed a similar trend. Performing noise characterization of optical frequency combs poses a challenge due to their complex phase noise dynamics and the corresponding scaling with comb-line number.^{11–15} Additionally, different comb lines may have different signal-to-noise ratios (SNRs), requiring measurement techniques to estimate phase noise over a wide range of SNRs accurately.^{16}

Understanding and characterizing the noise properties of optical frequency combs is paramount to harnessing their full potential. The magnitude of their phase noise is a crucial property, as its presence degrades the frequency stability of the individual comb lines, resulting in limited system performance.^{1,17}

Various phase noise measurement techniques have been proposed in the literature.^{18–20} As these methods are implemented in the analog domain, they introduce additional complexity to the experimental setup. More importantly, analog methods usually rely on optical bandpass filters to isolate single comb lines from the spectrum. In the aforementioned approaches, the measurement accuracy can be limited by the SNR of single comb lines. This is because the combination of low-power comb lines and measurement noise may lead to low SNR, leading to an inaccurate phase noise estimate. While cross-correlation methods can reduce the impact of measurement noise on the phase noise estimate, they require additional calibration and introduce complexity.^{21} Additionally, optically filtering out single lines can be infeasible for very narrowly spaced combs.

Another drawback of performing phase noise characterization of frequency combs based on line-by-line measurement is the inability to extract the time-domain correlations between neighboring lines. To remedy this, analog measurement methods have been developed to measure multiple comb lines simultaneously.^{22,23} However, these methods again introduce more experimental complexity while still being limited to measuring only a few lines.

The complementary approach of using balanced detection in combination with digital signal processing for phase noise estimation shifts most of the measurement complexity into the digital domain. The availability of ultra-wideband photodetectors and analog-to-digital (ADC) converters at the order of 100 GHz bandwidth allows the detection of many comb lines simultaneously in multi-heterodyne measurements. The sampled data are stored offline, allowing advanced digital signal processing (DSP) for joint phase noise estimation of multiple frequency comb lines to be efficiently implemented.^{7,24}

Applying DSP techniques thereby allows simultaneous phase noise characterization of multiple comb lines while retaining a minimal measurement setup.^{25} The extracted phase noise can then be used to compute the phase noise correlation matrix, which can then be used to identify the underlying phase noise sources and their comb line number dependent scaling.^{11,16}

However, some of the aforementioned DSP-based phase noise characterization methods employ digital bandpass filtering of individual comb lines, thereby inheriting some of the limitations of analog techniques. This implies that the SNR of the individual lines limits the estimation accuracy of each comb line’s phase noise, and narrowly spaced combs still pose an issue.

Bayesian filtering, implemented as an extended Kalman filter (EKF), offers joint digital signal processing of multiple comb lines. It eliminates the shortcomings of digitally filtering individual comb lines while maintaining the reduced experimental complexity. In particular, the EKF can provide near-optimal phase noise estimation in the presence of measurement noise.^{26} Given an analytical description (model) of the detected signal and the associated phase noise sources, the EKF can distinguish between phase and measurement noise. This distinction allows the digital suppression of measurement noise and an improved phase noise estimate over state-of-the-art conventional methods.^{16} The EKF inherently provides confidence intervals (CIs) with the phase noise estimates, allowing the estimate’s uncertainty to be quantified. In a laser phase noise characterization, the EKF was able to predict phase noise up to 7.6 dB below the measurement noise floor with 95% confidence.^{27}

In Ref. 16, a framework based on the EKF is presented to characterize the phase noise of an electro-optic (EO) comb. The EKF operates in the time domain using a joint model of the detected comb lines, so no prior bandpass filtering is required. The estimation accuracy is improved compared to the conventional DSP methods, and the phase noise correlations can be recovered. In this framework, however, the computational complexity scaling of the EKF is $O(N3)$, where *N* is the number of detected comb lines. This scaling is unfavorable and practically disqualifies the method for combs with many lines.

In Ref. 28, the EKF was used to compensate for phase noise in a dual-comb spectroscopy setup. Here, the complexity scaling issue was avoided by estimating the carrier-envelope frequency and the repetition rate phase noise of the comb instead of estimating each line. The complete reconstruction of a noise-free comb spectrum from a free-running comb is demonstrated. However, this method is tailored toward applying noise compensation in dual-comb spectroscopy rather than characterization.

This tutorial introduces an improved EKF-based framework for phase noise characterization of optical frequency combs. The framework’s significant improvement compared to Ref. 28 is the inclusion of automatic differentiation (AD). As the model within the framework may contain many free parameters that need to be determined, AD allows for gradient-based static parameter optimization. This approach to model parameter optimization significantly reduces complexity and enables efficient optimization in high-dimensional search spaces. Moreover, AD allows for the introduction of a latent phase noise space, a compressed representation of comb line phase noise. The structure of the latent space is learned directly from the measurement data. The introduction of an optimizable latent space not only solves the complexity scaling problem but also enables the search for phase noise modes beyond the commonly agreed upon noise sources in combs.^{11}

We apply AD and adaptive gradient-based optimization to learn the structure of the latent space. AD allows the optimization of an analytical model that describes the signal. Critically, AD enables the calculation of all model parameter gradients within any differentiable EKF model expression, eliminating the need for almost all auxiliary parameter estimation techniques. In this joint optimization process, the models of the measured signal and the statistical properties of the noise sources are learned. This versatility makes the method very general and, in principle, applicable to any comb.

In this work, we demonstrate in detail how to use this framework on different types of combs in simulation and experimental data. The remainder of this tutorial is structured as follows: In Sec. II, we introduce the required experimental setup and derive analytical expressions to describe the detected signal. We supply the necessary theoretical background for applying the EKF and optimizing a model. The framework is applied to a synthetic laser beat note measurement to detail every step in the process. Additionally, the state-of-the-art conventional phase estimation method is explained and applied. In Sec. III, the framework is applied to an EO comb on simulated and experimental data. Here, we demonstrate the versatility of the optimization process, the accuracy of the phase noise estimation, and the extraction of phase noise correlations. In Sec. IV, we show how the framework can be used in dual-comb measurements. Using two mode-locked lasers, we present the phase noise estimation process based on an experimental measurement. Furthermore, we discuss how this framework can be applied for digital noise compensation, similar to the method presented in Ref. 28. Finally, we give some concluding remarks and a brief outlook in Sec. V.

## II. METHODOLOGY

This section describes the complete framework in detail. First, the experimental setup common to all multi-heterodyne measurements considered in this work is described. A conventional phase noise measurement method using simple digital phase noise estimators is briefly explained. Finally, the EKF-based method is detailed. We describe how AD is applied in an efficient and versatile optimization framework. At the end of the section, we step through the example of a laser phase noise characterization, as it is conceptually simpler than the frequency comb case while still showcasing the intricacies of the method.

### A. Setup for numerical simulations and experimental measurements

The experimental setup described in this work is shown in Fig. 1(a). The same setup is also used for simulation studies. It is used to generate heterodyne beat note signals for use in frequency comb and laser characterizations. The setup requires minimal analog components, including an optical coupler, a BPD, and an ADC, making it easy to reproduce.

Throughout this article, the same experimental setup is used for all characterization studies, whether in simulation or experiment. All optical components in the setup are assumed to be fiber-bound, but a free-space setup would also work.

In the measurement process, the comb or laser under test (CUT/LUT) and the local oscillator (LO) are combined using a balanced fiber coupler. The combined signal is then detected by a balanced photodetector, which suppresses intensity noise and reduces shot noise contributions, resulting in an improved signal-to-noise ratio (SNR) of the signal.^{29}

The detection process produces one or multiple beat notes, depending on the number of optical tones in the light sources and their detuning with respect to each other. This setup is known as (multi-) heterodyne detection.

*ω*

_{ceo}, and an integer multiple of a repetition frequency

*ω*

_{rep}that defines the spacing of the comb lines. The complex electric fields for both sources can be expressed as

*P*

_{l}+ Δ

*P*

_{l}(

*t*) is the power of the

*l*-th line plus its time-dependent fluctuation, and

*ϕ*

_{l}(

*t*) is the time-dependent phase noise of the

*l*-th line. As the infinite sums suggest, an optical frequency comb may have arbitrarily many comb lines. In practice, however, only a subset of these lines is of interest. Often, experimental thresholds like a minimal detectable optical power or the optical frequency of given lines, limit which lines are detectable.

This general description can model any combination of CUT/LUT and LO. Mathematically, a LUT is equivalent to a CUT with only a single line. As this is a very general signal description, it can be unwieldy and impractical. We will, therefore, examine multiple special cases of this expression throughout this article and detail how this expression simplifies.

*t*

_{k}=

*k*Δ

*T*

_{s}indicates time discretization by the ADC with sampling period Δ

*T*

_{s},

*η*is a detection proportionality constant, and

*ξ*(

*t*

_{k}) is detector measurement noise.

*ξ*(

*t*

_{k}) contains contributions originating from thermo-electric noise in the detector, shot noise, and ADC quantization noise.

^{25,29}While bandwidth limitations in the detection system do not allow for an infinite number of detected comb lines, we keep the notation for simplicity.

If a single photodiode were used instead of a BPD, the detected signal would be proportional to |*E*_{a}|^{2}. This would result in constant power offset terms and no amplitude noise suppression. Additionally, single photodiode detection has a 3 dB SNR penalty compared to BPD. While BPD is the superior detection method, single PD detection can be used to the same effect if amplitude noise and SNR are not a concern.

*P*

_{l/m}(

*t*

_{k}) enter Eq. (3) only in square root order. In addition, they are typically significantly less prominent compared to phase noise, so they can often be neglected. Under this assumption, the expression simplifies to

The quantities of interest for the phase noise characterization of the CUT are $\varphi lCUT(tk)$. However, as both $\varphi lCUT(tk)$ and $\varphi mLO(tk)$ are randomly fluctuating phase noise terms, it is not possible to separate their difference $\varphi lmtotal(tk)$. As we cannot rely on additional prior knowledge about their statistics, they become indistinguishable. To still be able to estimate $\varphi lCUT(tk)$, one of two approaches should be chosen. For the first option, a LO needs to be selected whose phase noise is negligible compared to the CUT, in which case $\varphi lmtotal(tk)\u2248\varphi lCUT(tk)$. The other option is to use a reference with close to identical phase noise properties as the CUT. It is still impossible to separate their contributions in the time domain, as the time evolution of their phases is uncorrelated. However, both noise contributions will have the same statistical behavior, which leads to them having approximately identical phase noise power spectral densities (PSDs) $S\varphi CUT,l(f)\u2248S\varphi LO,l(f)$ per line *l*. As the total phase noise $\varphi lmtotal(tk)$ is the difference between the individual noise sources, the total PSD can be expressed as $S\varphi total,ll(f)=S\varphi CUT,l(f)+S\varphi LO,l(f)\u22482S\varphi CUT,l(f)$. Here it is immediately clear that in this specific case, a division by two yields the PSD of the *l*th line of the CUT. Notably, this argument is only effective when investigating the beat notes that originate from comb lines with the same properties, as indicated by the identical line index *l*. This scenario is useful when multiple copies of the CUT are available.

If neither of these options can be chosen due to a lack of suitable references, the measured phase noise will contain an unknown mixture of both sources. While this is not ideal for the characterization of the CUT, extracting the total phase noise per comb line can still be useful in certain applications. In dual-comb spectroscopy, for example, the amount of total phase noise per detected comb line can be a limitation.^{30} Here, the exact composition of $\varphi lmtotal(tk)$ from the CUT and LO is less important in practice.

### B. Conventional phase noise estimation algorithm

In this section, a brief summary of the conventional digital phase noise estimation algorithm is given.^{11,16,26} It is based on the discrete Hilbert transform $H$ to obtain the second quadrature of the detected multi-heterodyne beat, followed by an inverse tangent operation to obtain the phase. A schematic illustration of the algorithm is given in Fig. 1(c).

Starting from the detected multi-heterodyne signal as in Eq. (4), each beat note needs to be isolated. This is accomplished using a bank of digital bandpass filters. Each filter is centered around the beat note frequencies *ω*_{lm}. The bandwidth is determined by the frequency difference between the detected notes. Applying the filter bank produces a time-domain trace *y*_{lm}(*t*_{k}) for each beat note. Each *y*_{lm}(*t*_{k}) now only contains one beat note and represents one of the quadratures of the signal.

*π*/2 phase shift. The phase noise $\varphi lmtotal(tk)$ of each beat note can now be estimated using

*ω*

_{lm}

*t*

_{k}represents a linear detrending operation. The approach is comparatively low in complexity and does not require any prior knowledge about the signal.

Notably, this method does not attempt to remove measurement noise present in the signal. The measurement noise term *ξ*(*t*_{k}) in Eq. (4) is still present in each *y*_{lm}(*t*_{k}) after the bandpass filtering operation. In high SNR regimes where measurement noise is negligible, this method has been shown to be an optimal estimator of $\varphi lmtotal(tk)$.^{26} In moderate to low SNR regimes, however, measurement noise can dominate such that the estimator is no longer optimal. The resulting phase noise estimate then contains measurement noise contributions, as the estimator cannot distinguish between the two. The resulting phase noise estimate is biased and overestimates the phase noise magnitude. In regimes with low SNR and low phase noise, the method can even fail completely to produce a reasonable estimate.

Additionally, the bandpass filtering operation limits the maximum offset frequency at which the phase noise can be estimated. This frequency is half of the bandpass filter’s bandwidth, which is half the frequency distance of a given beat note to its neighboring beat note. This limitation becomes relevant for narrowly spaced combs. In dual-comb spectroscopy, the multi-heterodyne beat notes can have spacings of less than 100 Hz. This method would, therefore, only allow a phase noise characterization of up to 50 Hz.

Throughout this work, this method will be called the conventional method and serve as a benchmark for the EKF framework.

### C. EKF-based phase noise estimation framework

In signals with low SNR or low phase noise beat notes, a method that can reduce the impact of measurement noise and isolate phase noise is necessary. For this purpose, Bayesian filtering in the form of the *extended Kalman filter* (EKF) has been proposed. The EKF has been shown to provide near-optimal phase noise estimates in (multi-) heterodyne measurements in the presence of measurement noise.^{16,26} After formulating a *state-space model* that describes the detected signal in Eq. (4), it can be used to track dynamic quantities inside the model over time that cannot be measured directly. These are the phase noise terms and are often called the “hidden” state of the state space model. Simultaneously, the EKF attempts to isolate the measurement noise present in the signal based on its statistical properties. By separating the measurement noise from the phase noise, an accurate phase noise estimate can be produced even in low SNR regimes.

If the analytic expression underlying the detected signal was linear, the regular Kalman filter can be shown to provide the optimal estimate with regard to the mean squared error.^{31} In our case, however, the expression in Eq. (4) is highly non-linear. This is known as the non-linear filtering problem, which does not have an analytical solution. This implies that no optimal filter exists for our purpose. The EKF provides an approximate solution by locally linearizing the state-space model. While the optimality property is lost, this first-order approximation is typically sufficient to produce accurate estimates. We, therefore, call the estimate “near-optimal.”

The EKF operates directly on the detected signal in the time domain given in Eq. (4). All beat notes in the signal are simultaneously part of the state-space model, which makes this a joint estimation approach. This means that no prior bandpass filtering around the individual combines is required as a joint estimation of the phase noise of each beat note line is performed. Consequently, the method does not exhibit the same limitations induced by bandpass filtering as the analog methods or the conventional method.

#### 1. State space model

*state-transition function*

**f**and a

*measurement function*

**h**, which have the form

**y**

_{k}is the measured signal at time

*t*

_{k}, which is modeled by

**h**(

**x**

_{k},

**). In multi-heterodyne measurements,**

*θ***h**is described by the analytical form in Eq. (4). The arguments of

**h**are the vector

**x**

_{k}, which contains the values of the hidden variables at sample

*k*.

**x**

_{k}corresponds to $\varphi lmtotal(tk)$ in Eq. (4). The second argument

**is a vector of time-independent parameters of the model, which we will refer to as static parameters. Here,**

*θ***contains all other quantities in Eq. (4) that are not phase noise, like the beat note amplitudes and their frequencies. The additional term**

*θ***r**

_{k}describes an additive measurement noise contribution to the measured signal, which corresponds to the detector noise

*ξ*(

*t*

_{k}).

**r**

_{k}is a stochastic variable that is distributed according to a zero-mean Gaussian probability distribution $N$ with a covariance matrix

**R**.

The state-transition function **f**(**x**_{k}, ** θ**) models how the hidden variables evolve over time. In addition to the deterministic function

**f**, the state space model includes an additive stochastic noise term

**q**

_{k}. This so-called process noise is a random variable distributed as a zero-mean Gaussian with a covariance matrix

**Q**. This allows the modeling of hidden variables based on their statistical properties. In the multi-heterodyne case, Eq. (6) models the phase noise evolution.

^{26}Intuitively, phase noise accumulates random fluctuations over time. Modeling these additive fluctuations to be normally distributed is equivalent to modeling the phase noise evolution as a random walk. A beat note with phase noise that behaves as a random walk has a Lorentzian lineshape. This is the same lineshape profile a quantum fluctuation limited laser cavity produces, which further motivates this approximation.

^{32}It yields a very simple state-transition function

#### 2. Filtering equations

Once a state space model of the form given in Eqs. (6) and (7) that describes the detected signal has been defined, the EKF can be applied. Applying the EKF consists of iteratively applying the so-called filtering equations. Given an initial value for the phase noise **x**_{0}, applying the EKF to the first signal sample *y*_{1} produces an estimate for the phase noise **x**_{1} at the first sample. Recursively stepping through the equations produces a vector of phase noise estimates **x**_{1:T} for each time step, where *T* is the number of samples in the signal.

Furthermore, the EKF not only produces the point estimates **x**_{1:T}. The EKF estimate for each sample is a Gaussian distribution; therefore, a vector of covariances **P**_{1:T} is produced alongside the point estimates. For each time step *t*_{k}, the EKF produces a Gaussian probability distribution $N(xk,Pk)$, which describes the likelihood for a certain value of the estimate. This means that the phase noise prediction inherently includes a quantification of the uncertainty of the estimate. At each *t*_{k}, the standard deviation of the Gaussian estimate can be interpreted as a confidence interval.

**f**with respect to

**x**, evaluated at (

**x**

_{k},

**).**

*θ*Iterating through these equations and saving the estimates **x**_{k} and their covariances **P**_{k} then produces a series of phase noise estimates, including uncertainties that are based on the presence of measurement noise. The equations are given here for completeness; for an in-depth introduction and derivations, see Ref. 31.

#### 3. Static parameter estimation using automatic differentiation enhanced adaptive optimization

Thus far, we have assumed that the values of the static parameters ** θ** are known. This is generally not the case, and their values need to be determined. This step is crucial, as accurate values for

**are a strict prerequisite for accurate phase noise estimation.**

*θ*In the EKF framework, this is most commonly performed by maximizing the posterior distribution *p*(** θ**|

**y**

_{1:T}) of

**given the observations**

*θ***y**

_{1:T}. It represents the probability of a specific value of

**being responsible for producing the observations**

*θ***y**

_{k}. If there is no prior information about the optimal value of

**, the prior probability distribution**

*θ**p*(

**) is constant. According to Bayes’s rule, maximizing the posterior distribution then becomes equivalent to maximizing the likelihood**

*θ**p*(

**y**

_{1:T}|

**).**

*θ**a posteriori*(MAP) parameter $\theta \u0302MAP$. Working in logarithmic space increases numerical stability as values of

*p*(

**y**

_{1:T}|

**) can become very large.**

*θ**LL*

^{−}can be approximated with the expression

^{31}

**y**

_{1:T}.

A common technique to solve the optimization problem in Eq. (11) is the expectation–maximization algorithm. It is an iterative optimization scheme in which convergence to a (local) optimum is theoretically guaranteed. For certain state-space models, including the one used in Ref. 16, Eq. (11) can even be solved using a closed-form analytic expression.^{33}

For state-space models with many static parameters ** θ**, however, an analytic form to solve the optimization problem is not available. For high-dimensional

**, an efficient approach to solve Eq. (11) requires the explicit computation of the gradients**

*θ**∂LL*

^{−}/

*∂*

**. While it is possible to calculate the gradients manually, it requires the symbolic propagation of the gradients through the filtering equations (9). This quickly becomes infeasible for state-space models with many parameters**

*θ***.**

*θ*^{16}To alleviate this complication, the use of automatic differentiation (AD) to calculate

*∂LL*

^{−}/

*∂*

**has been proposed.**

*θ*^{34}Using AD allows the accurate and efficient evaluation of partial derivatives of any differentiable function specified in a computer program.

^{35,36}It keeps track of all elementary operations inside the function to build a graph of the function composition. By defining the differentiation rules of these elementary operations, the gradient of a function can be calculated by propagating through the composition graph and recursively applying the chain rule of differentiation. The key advantage of AD is its ability to differentiate arbitrarily complex functions. Instead of the error-prone derivation of symbolic derivatives, AD handles the computation automatically.

To evaluate *∂LL*^{−}/*∂*** θ**, the AD is propagated through the recursive EKF equations (9) that are required to evaluate

*LL*

^{−}. For our implementation of the EKF framework, we used the JAX library, which provides a fast and easy-to-use implementation of AD.

^{37}

With a way to calculate *∂LL*^{−}/*∂*** θ**, we are able to use any gradient-based optimization method. A robust and efficient optimization algorithm is required that performs well in high-dimensional search spaces since

**can potentially have many free parameters. Therefore, we propose the use of modern adaptive gradient-based optimizers, which are commonly used in the optimization of large neural networks. These algorithms are designed for fast convergence in high-dimensional optimization problems. In particular, we will use the AdaBelief optimizer throughout this article, as it converged very quickly and accurately in testing.**

*θ*^{37,38}

**. The new optimization problem then has the form**

*θ*The new parameters $\theta \u0303$ that are exposed to the optimizer are unconstrained real numbers. Choosing appropriate transformation functions Λ then allows the actual model parameters ** θ** to be scaled and constrained.

To illustrate this concept, we explain a possible choice of optimization transform function for an exemplary parameter. The variance of the measurement noise **R** in the measurement Eq. (7) of the state-space model needs to be strictly non-negative. In addition, its value can become so small that numerical accuracy may become an issue. Both of these issues are avoided if Λ is chosen to be the exponential function, which maps the real numbers to non-negative real numbers. On top of that, small values close to zero are mapped to large negative values. This effectively causes the optimization to be performed in logarithmic space. Λ can be a different transformation for every parameter in ** θ**, as long as it is a differentiable operation.

To reduce the computational load during optimization, a truncated version $LLtrain\u2212=LL\u2212(y1:Ntrain,\theta )$ of the full negative log-likelihood can be used, which only considers *N*_{train} < *T* samples of the measurement vector **y**_{1:T}. Under the assumption that the measurement and phase noise are approximately stationary, this truncation is a good approximation of the full negative log-likelihood. A balance between computational load and accuracy has to be found. Considering too few samples will result in a biased estimate of the model parameters, while using a large number of samples will result in unnecessary slow convergence. However, the long-term behavior of the system cannot be captured by the likelihood if too few samples are considered.

*LL*^{−} is generally a non-convex function, which means that the optimization problem potentially has multiple local minima. The optimization algorithm can, therefore, get stuck at a local minimum and cannot guarantee a global optimum. To mitigate this problem, we typically run the optimization multiple times with different initial guesses for the model parameters ** θ** and choose the model parameters that yield the highest likelihood.

To test whether a good optimum was found, the parameters $\theta \u0302MAP$ can be tested on a section of the measured signal **y** that was not used to calculate $LLtrain\u2212$ in the optimization phase. To test a set of parameters after convergence in the optimization phase, we simply evaluate $LLtest\u2212=LL\u2212(yNtrain:Ntrain+Ntest,\theta )$. *N*_{test} is the number of signal samples used and is usually chosen as *N*_{test} = *N*_{train}.

If $LLtest\u2212$ evaluates to a similar value as $LLtrain\u2212$, this is a good indication that a set of parameters was found that is general over the whole signal and a good optimum was reached. A similar value is also a good indication that the choice of *N*_{train} is reasonable. It means that the optimized parameters ** θ** are generalized and represent the stationary equilibrium of the signal.

### D. Application to laser phase noise characterization

To finalize the methodology section, we will step through a sample case of laser phase noise characterization in detail. The setup in Fig. 1 will be simulated for the case where the LUT and LO both are single-frequency lasers. Subsequently, we will apply the conventional and EKF-based phase noise estimation methods and compare their accuracy.

#### 1. Heterodyne laser beat note signal simulation

*ϕ*

^{LO}(

*t*

_{k}) to be negligible, such that

*ϕ*

^{total}(

*t*

_{k}) ≈

*ϕ*

^{LUT}(

*t*

_{k}). As the simulation parameters, we choose the beat note amplitude

*A*= 10 mV, the beat frequency

*ω*

_{beat}= 220 MHz, and a sampling frequency of

*f*

_{sampling}= 1 GHz. Both the measurement noise

*ξ*(

*t*

_{k}) and the phase noise

*ϕ*

^{LUT}(

*t*

_{k}) are random variables. The measurement noise samples

*ξ*

_{k}have zero auto-correlation in time and are independently drawn from a zero mean normal distribution with variance $\sigma R2=5\u22c510\u22126$. This corresponds to a measurement noise PSD level of −140 dBV

^{2}/Hz, which is a typical experimental value for state-of-the-art measurement equipment.

It should be noted that in the heterodyne beat note case, where the signal as expressed in (16) is just a single sinusoidal, it is also possible to demodulate the signal experimentally. This would result in a measured signal *y*(*t*_{k}), which is linear in *ϕ*^{LUT}(*t*_{k}) and would, therefore, not require the full EKF formalism. Instead, the (linear) Kalman filter could be used, simplifying the DSP. However, in the general multi-heterodyne case, this is no longer possible, which is why we introduce the non-linear EKF formalism here.

^{39}For simplicity, however, we consider a phenomenological model for the frequency noise PSD of the phase noise of the laser. Using

_{flicker}

*f*

^{−1}and a Lorentzian contribution Ξ

_{Lorentz}.

^{40}The frequency-independent term Ξ

_{Lorentz}models the theoretical phase noise contribution of an ideal quantum noise limited laser cavity.

^{32}This term causes the line shape to have a Lorentzian profile. In addition, it can be used to define the intrinsic linewidth of a laser as Δ

*ν*=

*π*Ξ

_{Lorentz}. In real-world systems, low-frequency perturbations typically cause additional laser phase noise, which is reflected in the flicker noise term.

To generate a time series realization of *ϕ*^{LUT}(*t*_{k}) that follows the model in Eq. (17), we use the method given in Ref. 41. In this study, we choose Ξ_{flicker} = 10^{6} Hz^{2} and Ξ_{Lorentz} = 100/*π* Hz to generate *ϕ*^{LUT}(*t*_{k}). Substituting into Eq. (16) produces the simulated signal *y*(*t*_{k}), which is displayed in Fig. 2, together with its PSD.

#### 2. EKF-based phase noise estimation

*y*(

*t*

_{k}). The first step is the definition of a state-space model. We first define a measurement function

**h**that describes

*y*(

*t*

_{k}). As the analytical expression in Eq. (16) was used to simulate

*y*(

*t*

_{k}), we can use this expression as the measurement function. Reformatting Eq. (16) to comply with the notation used in Sec. II C yields

Here, *x*_{k} is the hidden state of the EKF and corresponds to the scalar-valued phase noise *ϕ*^{LUT}(*t*_{k}). *r*_{k} represents the scalar-valued measurement noise *ξ*(*t*_{k}) with variance *σ*_{R}.

**f**. It models the time evolution of the phase noise and should ideally be as close as possible to the true time evolution. In this simulation, the true phase noise model is given in Eq. (17). We will, however, assume that the model is unknown to us and use the random walk approximation introduced in Sec. II C 1. Therefore, the state-transition function is given by

*x*

_{0}and its variance $\sigma P0$ as well as the noise variances of the state space model were included. Instead of optimizing the beat note frequency

*ω*

_{beat}directly, an initial estimate $\omega \u0303$ is obtained by running a peak-finding routine on the signal PSD displayed in Fig. 2. We then define $\omega beat=\omega \u0303+\delta \omega $, where $\omega \u0303$ is fixed and

*δω*allows for fine-tuning of the frequency. To ensure good convergence behavior in the optimization and to constrain certain parameters, we define the parameter-wise optimization transformation function

*δω*is constrained to values between −

*δω*

^{max}and

*δω*

^{max}using the hyperbolic tangents function. We choose

*δω*

^{max}= 5 kHz as the initial estimate $\omega \u0303$ is quite accurate.

With the state model and the parameter transformation function defined, the optimization loop illustrated in Fig. 1 is executed using the AdaBelief optimizer. To minimize the chance of sub-optimal convergence by running into a local optimum, 100 sets of randomly drawn initial parameters *θ*_{0} are chosen. Each parameter set is independently optimized until the change in parameters per iteration is negligible. The convergence behavior for a selected subset of parameters is illustrated in Fig. 3. The initial values are distributed over a wide range of multiple decades to demonstrate the optimization routine’s robustness. Some of the optimization runs for which the random initialization of *θ*_{0} is particularly poor do not converge. Most runs, however, converge to one of two local optima, which can be seen in the two distinct final values for the variances in Figs. 3(b) and 3(d). Ultimately, only the final parameters ** θ** of the best-performing run are used.

Using the optimal parameters $\theta \u0302MAP$, we can apply the EKF to estimate the time evolution of the phase noise as the final step. We receive a time series of estimates *x*_{k} with a corresponding time series of variances $\sigma P,k2$. Together, they form a time series of Gaussian distributions $N(xk,\sigma P,k2)$ of near-optimal estimations of the phase noise time evolution. To put the estimation into perspective, we apply the conventional phase estimation algorithm described in Sec. II A to the same simulated signal. Both methods are compared to the true phase noise in Fig. 4. The time domain traces shown in (a) indicate that the long-term behavior of the phase is evaluated accurately by both methods. Zooming in, however, reveals that the EKF estimate is more accurate than the conventional method. The frequency noise PSD picture in Fig. 4(b) shows that the low-frequency dynamics up to 300 kHz are accurately captured by both methods. For higher offsets, the conventional method is limited by the measurement noise floor *ξ*(*t*_{k}), which causes the estimate to have a ∝*f*^{2} frequency dependency. The EKF estimate follows the ground truth below the noise floor, demonstrating that the measurement noise is effectively filtered out. The limits of the estimation accuracy are quantified by the variances $\sigma P,k2$ of the estimates. Using second-order error propagation, the time domain uncertainty can be propagated to the PSD estimate, allowing us to draw a CI around the PSD. A derivation is presented in the Appendix.

In summary, the presented method using the EKF produces a more accurate estimate than the conventional method. The introduction of adaptive optimization algorithms and the use of AD make the method very versatile, as state space models of arbitrary form with arbitrarily many parameters can be efficiently optimized. In Sec. III, we will demonstrate how this framework can be used to estimate phase noise in frequency combs.

## III. DOWN-CONVERTED COMBS

The beat note generation process described in Eq. (21) is illustrated in Fig. 5. The resulting signal is once again a comb with the same spacing as the original line spacing. To simultaneously detect as many lines as possible, we want to maximize the number of down-converted lines within the detection bandwidth. This is achieved by selecting a LO with a suitable frequency *ω*^{LO} that is close to the frequency of the highest power comb lines of the CUT. All lines for which |*ω*_{l}| < 2*πf*_{BW} can then be detected. The maximum number of detectable lines within the detection bandwidth *f*_{BW} is *N*_{detected} ≤ 2⌊2*πf*_{BW}/*ω*_{rep}⌋ + 1. While we can always detect at least one comb line, this relation clearly illustrates that this setup is only suitable for combs that have a repetition frequency significantly below the available bandwidth. Only then can we detect multiple lines at the same time and benefit from the reduced experimental complexity of the setup to infer information about the comb phase noise from only a single measurement.

This condition disqualifies this setup for types of combs that have a repetition frequency *ω*_{rep} in the tens of GHz and above, as even state-of-the-art photodiodes are limited to around 100 GHz of bandwidth.^{42} Even with these extreme bandwidths, only a few comb lines would be detectable. Therefore, combs with large spacing can effectively only be characterized in dual-comb setups, which are discussed in Sec. IV.

Furthermore, this setup only has limited usability for combs with many lines on the order of thousands, which includes femtosecond laser-based combs or generally octave-spanning combs. Their repetition rate is typically on the order of GHz or less, and beating with a single laser will produce many detectable beat notes. However, a single measurement will only ever cover a limited range of lines in the comb under test. For a full characterization of all lines, multiple multi-heterodyne measurements with a wavelength-tuneable LO would be required to access the phase noise of all lines. While not infeasible, a dual-comb setup with a suitable reference comb is preferable in this scenario.

For combs with *ω*_{rep} up to 10 GHz and few comb lines, however, utilizing this single-frequency laser setup is more advantageous compared to employing a dual-comb configuration. With a single-frequency laser acting as the local oscillator (LO), it is still possible to detect numerous lines while the experimental setup is simplified significantly. Furthermore, assessing the contribution of the LO toward the total phase noise per beat note $\varphi ltotal(tk)$ is more straightforward. As derived in Eq. (21), the LO phase noise *ϕ*^{LO}(*t*_{k}) is an additive contribution independent of the comb line index *l*. It, therefore, only adds a common phase noise offset that does not add any index-dependent dynamics. This distinction is generally not as easy in a dual-comb setup, as the LO itself is a comb and its phase noise is, therefore, dependent on the line index.

To show how the EKF-based method can be used to extract the phase noise term $\varphi ltotal(tk)$ in this scenario, we will first demonstrate its application on a simulated signal. An electro-optic (EO) comb is chosen as the CUT, as it is well-studied in the literature and has a known phase noise model. After showcasing the method on simulated data, the method’s application is also demonstrated on experimental data.

### A. Simulating an EO comb

^{43}

*V*(

*t*) is the sinusoidal voltage that is generated by the RF source and applied to the EOM, which has a peak voltage

*V*

_{0}and a frequency

*ω*

^{RF}. The EOM acts as a pure phase modulator.

*K*is a material-dependent proportionality constant of the EOM that determines the phase change per applied voltage. The LO is simply a single-frequency laser, and its electric field

*E*

_{LO}(

*t*) is given by Eq. (15).

*E*

_{CUT}(

*t*) and

*E*

_{LO}(

*t*) are propagated through the setup using Eqs. (2) and (3). After balanced coupling, detection in the BDP, and digitization by the ADC, the final signal then has the form

*ϕ*

^{LO}(

*t*

_{k}) of the LO is insignificant compared to the phase noise of the seed laser

*ϕ*

^{CW}(

*t*

_{k}) for simplicity such that

*ϕ*

^{CW}(

*t*

_{k}) −

*ϕ*

^{LO}(

*t*

_{k}) ≈

*ϕ*

^{CW}(

*t*

_{k}). For the general case where the LO noise cannot be neglected, we will keep

*ϕ*

^{LO}(

*t*

_{k}) as part of the model throughout this section. For

*ϕ*

^{CW}(

*t*

_{k}), we use the same phase noise model as in the laser beat note study in Sec. II D. We use the model described by Eq. (17), in which the phase noise is composed of a random walk with added low-frequency flicker noise. In this study, we choose Ξ

_{flicker}= 10

^{6}Hz

^{2}and Ξ

_{Lorentz}= 10

^{3}/

*π*Hz as the parameters for the model.

The third phase noise term that appears in Eq. (24) is hidden in *V*(*t*_{k}), the output of the RF source. As Eq. (23) states, *ϕ*^{RF}(*t*_{k}) is the phase noise in the RF signal. RF source phase noise generally has different and more complicated characteristics compared to laser phase noise. Therefore, we will use a different phenomenological model to generate *ϕ*^{RF}(*t*_{k}).

^{44}A physically accurate simulation would heavily depend on implementation details and is out of the scope of this work. To still roughly emulate this behavior and to intentionally deviate from the laser phase noise model,

*ϕ*

^{RF}(

*t*

_{k}) is simulated using

_{Lorentz,RF}and Ξ

_{Plateau,RF}. The former models an underlying random walk of the phase noise. The latter models the scale of a sinusoidal deviation that is added to the random walk. The factor

*α*= 1 Hz

^{−1}only ensures that the arguments to the logarithm and sine are unitless. For this numerical study, we choose the free parameters as Ξ

_{Lorentz,RF}= 100 Hz and Ξ

_{Plateau,RF}= 10

^{4}Hz with the corresponding unit, which creates a phase noise PSD model that has multiple plateaus.

Note that this model and its parameters are phenomenologically chosen to emulate the phase noise of a specific device, which we have used as Ref. 45. The secondary purpose of this model is its explicit deviation from a simple phase noise random walk. As discussed in Sec. II D, the EKF state-transition function we use to approximate the phase noise behavior does not consider any deviations from a random walk. An intentional mismatch between the true and approximated model can, therefore, give an indication of whether the random walk state-space model can still capture complicated phase noise dynamics.

To simulate time series realizations *ϕ*^{CW}(*t*_{k}) and *ϕ*^{RF}(*t*_{k}) based on the PSD models in Eqs. (17) and (25), we again use the algorithm given in Ref. 41. Equation (24) can now be used to directly simulate the time domain signal of a down-converted EO comb. The parameters chosen for this numerical study are $2\eta PCWPLO$ = 100 mV, *KV*_{0} = 1.6, *ω*^{CW} = 2*π* × 1.1 GHz, and *ω*^{RF} = 2*π* × 500 MHz. Here, we have combined some individual parameters into effective combined quantities, as only total values, as opposed to individual contributions, matter. The spectrum of the resulting signal *y*(*t*_{k}) is shown in Fig. 6(b). The resulting comb is centered around the line with index *l*_{center} = 2 and features *N*_{alias} = 3 and *N*_{reg.} = 8 aliased and non-aliased comb lines, respectively.

### B. Latent phase noise state-space model

To use the EKF-based phase noise estimation method, a state-space model needs to be defined. The measurement equation of the model, which describes the detected signal, is based on the general analytical expression in Eq. (24).

*ϕ*

_{l}(

*t*

_{k}) of the individual beat notes. They are approximated by a multi-dimensional random walk. The state-space model takes the form

*N*=

*N*

_{alias}+

*N*

_{reg.}dimensional vector of hidden variables

**x**

_{k}as the beat note phase noise terms

**(**

*ϕ**t*

_{k}).

*N*

_{alias}and

*N*

_{reg.}are the number of aliased and non-aliased comb lines with significant amplitude in the recorded spectrum, respectively.

*ω*

_{0}is chosen as the lowest of the beat note frequencies that is not aliased (as indicated in Fig. 5). The static parameters are

*θ*_{indiv.}= {

**A**,

*ω*

_{0},

*ω*

_{rep},

**x**

_{0},

**P**

_{0},

**Q**

_{N},

*σ*

_{R}}. This model has $|\theta individual|=12N2+7N+3$ free parameters that need to be optimized.

Using this state-space model for EKF-based phase noise tracking has been shown to produce accurate estimations that outperform the conventional method in an EO comb.^{16} The approach does, however, have multiple drawbacks that have to be addressed.

In this individual line tracking approach, the correlations between the phase noise terms $\varphi ltotal(tk)$ are found through the optimization of **Q**_{N}. This matrix represents the differential phase noise covariance, from which the differential correlations can be calculated. No prior information about the phase noise in the system is used in this model.

The state vector **x**_{k} has the same dimension as the number of detected comb lines *N*. This causes the associated covariance matrices **P**_{k} as well as the process noise covariance matrix **Q**_{N} of the state vector to be a matrix of size $N\xd7N$. The EKF filtering equations (9) include matrix operations that have an undesirable computational complexity scaling of $O(N3)$. Applying the EKF to a large number of comb lines *N* can make this approach computationally infeasible. Additionally, the covariance matrix **Q**_{N} needs to be optimized, which requires the optimization of $O(N2)$ free parameters.

To solve the complexity scaling problem, we can exploit the fact that the individual comb line phase noise terms $\varphi ltotal$ is a combination of a few underlying phase noise sources. The expression for the signal in Eq. (24) shows that only three phase noise sources contribute. The *N* individual beat note phase noise terms $\varphi ltotal(tk)$ are a combination of the three sources, which means that they are correlated.

^{19,46}The phase noise of each beat note is

**C**

_{true}is zero at the index of the central comb line

*l*

_{center}. We have introduced a mixing matrix

**C**which describes the linear transformation that connects the true phase noise sources with the phase noise per beat note. Independent of the number of beat notes

*N*, each $\varphi ltotal(tk)$ is described by only two independent contributions: A common mode noise term

*ϕ*

^{comm}(

*t*

_{k}) =

*ϕ*

^{CW}(

*t*

_{k}) −

*ϕ*

^{LO}(

*t*

_{k}) that is added to all beat notes independent of

*l*, and a differential mode noise term

*ϕ*

^{RF}(

*t*

_{k}) which has a linear dependence on

*l*.

**ϕ**

^{total}(

*t*

_{k}) can be constructed if

**C**

_{true}is known. In more formal terms, we introduce a linear latent space into the state-space model, which is a lower-dimensional representation of the comb line phase noises

**ϕ**

^{total}(

*t*

_{k}). In general, this can be expressed as

**ψ**(

*t*

_{k}) is a

*L*<

*N*dimensional vector of latent phase noise sources.

**x**

_{k}is the

*L*dimensional vector of latent phase noise sources

**ψ**(

*t*

_{k}). $\varphi l0$ is the global initial phase of the

*l*th line. While not strictly required, as the initial values

**x**

_{0}can also contribute to an initial phase, the inclusion of $\varphi l0$ improves optimization convergence.

*θ*_{lat.}= {

**A**,

**ϕ**

^{0},

*ω*

_{0},

*ω*

_{rep},

**x**

_{0},

**P**

_{0},

**C**,

**Q**

_{L},

*σ*

_{R}}. The dimension

*L*of the latent space is a hyperparameter of the model. This means that it cannot be optimized directly but instead needs to be chosen based on prior information. In the EO comb case, Eq. (29) shows that a two dimensional latent space with

*L*= 2 can explain the phase noise of all comb lines. Furthermore, the form of the mixing matrix

**C**is given by Eq. (30). The only free parameter in this expression is the index of the central comb line, which shifts the zero-crossing point of the second column of the matrix. As

**C**is part of

*θ*_{lat.}, the index offset can be optimized alongside the other parameters. To perform the optimization, we define the optimization transformation function Λ as introduced in Sec. II C 3,

**P**

_{0}and

**Q**

_{L}are matrices, as the hidden state dimension is larger than one.

*δ*

_{ij}is the Kronecker delta, which forces both

**P**

_{0}and

**Q**

_{L}to be diagonal matrices. This explicitly disallows correlations between the latent noise sources

**ψ**(

*t*

_{k}). While it is possible to allow non-zero off-diagonal elements, this complicates the optimization process as

**Q**

_{L}has to be constrained to the space of positive semi-definite matrices as it is a covariance matrix. Only allowing diagonal matrices avoids this problem and additionally reduces the number of parameters to optimize. This is equivalent to constraining the latent phase noise sources to be uncorrelated. For the EO comb, this is a reasonable assumption as the true phase noise sources are also uncorrelated. To optimize the comb frequencies, we again choose

*δω*

^{max}= 5 kHz.

The function Λ_{C} that transforms the mixing matrix is based on Eq. (30). Only a single free parameter, *l*_{center}, is required to define the matrix. The latent phase noise space is constrained to only have a common mode and a differential mode phase noise contribution. We will call this the constrained latent space model, indicated by the subscript **C**_{c.}.

The constrained model has |*θ*_{c.}| = 2*N* + 3*L* + 4 free parameters and *L* tracked variables. For large *N*, the linear scaling of the number of parameters causes it to have significantly fewer free parameters compared to the individual line approach. Notably, the number of parameters no longer scales quadratically with the number of comb lines. Furthermore, the EKF computational complexity reduces to $O(L3)$, which is no longer dominated by scaling with the number of comb lines. This means that it becomes feasible to use this model to apply the EKF-based phase noise estimation method to combs with a large number of comb lines.

Unconstrained latent spaces While the latent phase noise space of an EO comb can be explicitly described by Eq. (30), a similar expression might not be available for all types of combs. In particular, evidence for latent phase noise sources beyond the common mode and differential mode has been reported for mode-locked lasers.^{11,47}

_{C}as defined in Eq. (34), it can be replaced by the identity function

**C**becomes a free parameter of the state-space model. The latent space decomposition becomes fully data-driven as the optimizer learns the entries of

**C**. No assumptions about the comb line phase noise are made, except for the assumption that only a few latent sources contribute to each comb line. We will call this variant the unconstrained model, indicated by the subscript

**C**

_{uc.}.

In this form, the model has |*θ*_{uc.}| = *NL* + 2*N* + 3*L* + 3 parameters, which still has a favorable scaling compared to the individual line tracking approach as *N* becomes large. In addition, the complexity of applying the EKF is still $O(L3)$.

### C. Application on simulated EO combs

To evaluate the latent space phase noise model, we apply it to the simulated EO comb discussed in Sec. III A. We will compare two variants of the latent space model against ground truth phase noise.

The constrained latent state-space model assumes that the comb line phase noise is explained by a common and a differential phase noise mode. This model is derived from the known behavior of EO comb phase noise as expressed in Eq. (28). The state-space model is defined in (32) and (33). Together with the optimization transition function, which is defined in Eq. (34), the parameters *θ*_{c} are optimized, as detailed in Sec. II C 3, and used to estimate the latent phase noise sources $\psi \u0302c(tk)$.

For comparison, we also test the unconstrained state-space model variant on the same simulated signal. This variant uses the same state-space model. The difference lies in the optimization transformation function, where we use the replacement detailed in Eq. (35). While the latent space is still kept two-dimensional, it is not constrained to common and differential noise modes. Using the optimization framework, we find the parameters *θ*_{uc.} and use them to estimate the latent phase noise $\psi \u0302muc.(tk)$.

Figure 7 displays the results of the EKF-based estimation. In Fig. 7(a), the columns of the optimized mixing matrix **C**_{c.} are compared to the true matrix. The single free parameter **C**_{c.} in Eq. (34) that identifies the index of the central comb lines was optimized. Its final value, $lcenterMAP\u22482$, is very close to the true value. Therefore, the columns of **C**_{c.} align perfectly with the true values. In (b), the frequency noise PSDs of the latent phase noise sources are shown. As **C**_{c.} matches **C**_{true}, the latent noise sources map directly to the physical phase noises $\varphi truecomm(tk)=\varphi trueCW(tk)\u2212\varphi trueLO(tk)$ and $\varphi trueRF(tk)$. At high offset frequencies *f* > 100 MHz, the EKF estimates deviate from the true PSDs. The CI around the PSD, however, captures this estimation inaccuracy.

In Fig. 7(c), the optimized mixing matrix **C**_{uc.} of the unconstrained model is shown. The column vectors do not display a well-defined separation into common and differential noise modes. The latent phase noise sources will, therefore, not correspond to physical noise sources. This can be seen in (d), where their PSDs are shown. The shape of the PSDs is similar to the true ones, but they deviate considerably. In addition, both latent phase noise PSDs display strong spikes at high frequencies, which are not present in true phase noise sources. The origin of these high-frequency oscillations is unclear, as they are not part of the state-space model. It is likely that due to the sub-optimal set of parameters *θ*_{uc.}, the EKF produces erratic phase estimates to compensate for the accordingly inaccurate state-space model. This showcases the importance of finding good model parameters since phase noise estimates cannot be trusted otherwise.

While it can be useful to isolate the physical noise sources in a comb, this is only possible if *a prior* model of the comb phase noise exists, as in Eq. (28). In all other cases, especially in the unconstrained mixing matrix case, the optimizer will find any **C** that can approximately decompose the comb line phase noise into latent noise sources.

Here, the latent space is not interpreted as a physical quantity. However, the observable of interest is the phase noise per comb line $\varphi ltotal(tk)$. The latent space simply serves as a compressed representation. Using Eq. (31), $\varphi ltotal(tk)$ can easily be computed from the latent phase noise.

#### 1. Estimation accuracy

In Fig. 8, the estimation accuracy of both the constrained and unconstrained EKF-based approach and the conventional method is compared against the ground truth. The phase noise per comb line can also be estimated using the conventional method discussed in Sec. II A. The accuracy is estimated based on the estimated variance of the phase noise.

In Fig. 8(a), the integrated phase noise PSDs of each individual comb line for the different methods are shown. This measure is similar to the empirical phase noise time domain variance as it has the same unit, except that it considers a wide part of the spectrum instead. The phase noise PSD of each comb line $S\varphi total,l(f)$ is integrated from *f*_{min} = 100 kHz to *f*_{max} = 100 MHz. The value of *f*_{max} is chosen since 100 MHz is the maximum offset frequency we can evaluate using the conventional method in this case.

*l*as

*A*

_{l}is its amplitude, $\sigma R2$ is the variance of the white measurement noise, Δ

*T*

_{s}is the sampling period, and

*f*

_{BW}is the detection bandwidth. While this seems counterintuitive, this impressive measurement noise rejection can be explained by the restrictive state-space model. Since

**C**

_{c.}is determined by only one parameter, there are not enough degrees of freedom in the model to allow for phase noise overestimation. Intuitively, the single free parameter is determined during the optimization based on the high SNR lines. The phase noise of the noisy low SNR lines is then extrapolated based on the constrained phase noise model, such that the low SNR of the lines does not impact the estimation process.

The unconstrained variant with **C**_{uc.} produces accurate estimates for the central comb lines, which have a high SNR. For the outer lines with SNR < 15 dB, the estimated variance is larger than the true value, which indicates that the phase noise magnitude is overestimated. The conventional method estimates are similar, where the high SNR comb line phase noise is estimated accurately. For the outer lines, however, the method completely fails. The very low SNR indicates that the measurement noise dominates over the beat notes, a regime in which the conventional method cannot produce a good estimate.

In Figs. 8(b) and 8(c), two exemplary phase noise PSDs of the comb lines with indices *l* = −1 and *l* = 3 are shown. The former beat note has a low SNR, which is, however, still enough such that the conventional method does not fail completely. Beyond *f* = 2 MHz, its PSD estimate is dominated by measurement noise. Both EKF variants follow the true PSD, which is below the noise floor. The beat note with index *l* = 3 is one of the highest SNR notes. The conventional estimate is considerably closer to the true PSD, as the noise floor is lower than that of line *l* = −1.

#### 2. Correlation matrices

**ϕ**(

*t*

_{k}) is the N-dimensional vector of differential comb line phase noise at time

*t*

_{k}.

**R**are then calculated component-wise from the covariance matrices using

In Fig. 9, the correlation matrix estimates for the different phase estimation methods are shown. Both correlation matrices from the EKF-based variants are similar to the true matrix. In the unconstrained case, however, the correlations for the line with index *l* = 7 deviate from the ground truth. This line has a very low SNR, which causes problems in the optimization process. The optimized parameters that correspond to this line are harder to find, which causes the inaccuracy.

The conventional phase estimations’ correlations match the true correlations for the central lines. As discussed before, the method fails for low SNR lines, resulting in completely uncorrelated estimates for the outer comb lines. In summary, the accuracy of the correlation matrix estimated seems to follow the phase noise variance estimate that is shown in Fig. 8. This is not surprising, as inaccurate phase noise estimates will lead to inaccurate estimates of their correlations.

### D. Application to an experimental EO comb

In this section, we demonstrate the robustness of the method and its applicability to real-world experimental data. The EKF-based method will be applied to experimental traces of a down-converted EO comb with similar characteristics to the simulated comb discussed in Sec. III A.

The experimental setup for the measurement considered in this section is identical to the setup displayed in Fig. 6(a). A *NKT Photonics Koheras BASIK E15* fiber laser is mixed with an EO comb. The EO comb is generated using a separate *NKT Photonics Koheras BASIK E15* laser and an EO modulator, driven by an *Agilent E8247C* signal generator. The mixed signal is detected by a balanced photodetector and digitized using a sampling oscilloscope with *f*_{BW} = 13 GHz and *f*_{sampling} = 40 GHz. After digitization, the recorded signal is downsampled digitally by a factor 3, as no beat notes appear in the high-frequency part of the signal spectrum. This step reduces the bandwidth of the measurement, which increases the SNR. The spectrum of the resulting measured signal is displayed in Fig. 10(a). We identify 13 non-aliased down-converted comb lines and 13 aliased lines, which are centered around the line with index *l*_{center} = 6. Notably, all lines except lines four to eight have low powers and, therefore, low SNRs. Additionally, each aliased line happens to be only about 36 MHz away from a neighboring non-aliased line. Since most of the regular lines have a higher SNR, they dominate over the aliased lines, which makes them hard to characterize. This is a problem when using the conventional method since it relies on bandpass filtering around each line. Because of the narrow spacing, the maximum bandwidth of each filter is limited to 36 MHz. This limits the maximum measurable phase noise PSD offset frequency to half the bandwidth. For this measurement, this means that the conventional method can be used to characterize each comb line only up to 18 MHz. While in practice, higher bandwidths could be used, the resulting phase noise estimate would be biased since it would include an aliased line.

This issue could be mitigated in this particular experimental setup by adjusting the central frequency of the seed laser. Nevertheless, in frequency combs with low repetition rates, the narrow spacing between lines remains an inherent challenge. Consequently, this measurement highlights the benefit of employing the EKF, as it eliminates the need for bandpass filtering and, therefore, circumventing the problem.

As long as the state space model correctly describes the detected signal, it does not matter how close the lines are to each other in frequency. We can, therefore, separate the contribution of each line and, in principle, isolate each line from the signal. However, here we are interested in their phase noise properties.

To apply the EKF to the detected signal, we use the same steps as in Sec. III A. The same state-space model as defined in Eqs. (32) and (33) is used, which describe the latent space phase noise model. As shown in Sec. III C, the constrained model variant produces a more accurate phase noise estimate, especially for very low SNR comb lines. In this section, we will, therefore, only use this variant. The corresponding optimization transformation function Λ is defined in Eq. (34). After optimizing the model parameters ** θ**, the EKF is applied to estimate the phase. As we used the constrained model for the mixing matrix

**C**

_{c.}, the latent phase noise sources can be identified as the common and differential mode noise sources. Here, the common mode corresponds to the difference between the phase noise of the comb seed laser and the LO

*ϕ*

^{CW}(

*t*

_{k}) −

*ϕ*

^{LO}(

*t*

_{k}). Both contribute equally to all down-converted comb lines, as in Eq. (28). The differential mode is the phase noise originating from the RF source

*ϕ*

^{RF}(

*t*

_{k}).

Both the common and differential mode contributions are separately accessible in the experiment, so we can record independent phase noise reference measurements. The common mode phase noise is estimated by creating a single beat note signal between the seed laser and the LO, which has the same total phase noise as the common mode of the EO comb. For the RF source reference phase noise measurement, the source is connected directly to the ADC. The phase noise of both reference measurements is estimated using the conventional method.

The frequency noise PSDs of the estimated latent phase noise sources are shown in Fig. 10(b), alongside the reference measurements. The PSDs of the EKF estimates are accurate over a wide frequency range, as they match the reference measurements very well. The common mode latent phase noise *ψ*_{0} matches the seed laser phase noise over the full frequency range. Up to 1 MHz offset frequency, the differential mode phase noise *ψ*_{1} matches the RF source phase noise. Beyond 1 MHz, the reference measurement is limited by measurement noise, which can be seen by the ∝*f*^{2} dependency of the frequency noise PSD. At the same frequency, the PSD of *ψ*_{1} flattens out and its CI grows larger. Still, the interval predicts the PSD to be significantly below the noise floor of the reference measurement. A flat frequency noise PSD corresponds to the random walk behavior that is part of the EKF state space model in Eq. (32). It is, therefore, not surprising that the EKF PSD estimate takes on this behavior: Once the phase noise is limited by measurement noise, extrapolating from the state-space model is the only remaining information the EKF can use. The CI does not reflect this model assumption. Its calculation, as presented in the Appendix, does not include uncertainty in the parameters of the model itself. It represents uncertainty originating only from the measurement noise while assuming complete certainty about the state-space model. While a random walk is a valid assumption for laser phase noise, this is not necessarily the case for an RF source, as discussed in Sec. III A. The latent phase noise PSDs below the measurement noise floor should only be trusted if the phase noise behavior is well-modeled in the state-space model. Since we cannot assume that the random walk is a good model for the RF source phase noise, we cannot trust the PSD of *ψ*_{1} beyond 1 MHz.

Next, we investigate the estimation accuracy with respect to the phase noise of the individual comb lines. Just as for the simulated EO comb, the phase noise per comb line can be calculated using Eq. (31). They are compared to estimations using the conventional method.

Figure 11(a) shows the integrated phase noise variance per comb line for the EKF and conventional method. The phase noise variance per comb line is calculated by integrating the phase noise PSDs of each comb line between *f*_{min} = 4 kHz and *f*_{max} = 15 MHz. For the central comb lines with SNR > 15 dB, both methods agree on the phase noise variance. As we saw for the simulated EO comb in Fig. 8, both methods produce accurate estimates in high SNR regimes. Therefore, we can assume that both methods produce an accurate estimate for lines four to eight. Beyond these lines, the conventional method fails. Since we have used the constrained EKF model, its estimate follows the theoretical phase noise model of the EO comb in Eq. (28). Just as in the simulation, the variance follows the expected parabolic form. Since this is an experimental measurement, we do not have access to ground truth values. The simulation results from Sec. III C, however, indicate an accurate prediction.

In Figs. 11(b) and 11(c), the estimated correlation matrices of both methods are shown, which are calculated according to Eqs. (37) and (40). As expected, the EKF predicts strongly correlated comb line phase noise. The conventional method can only recover the correlations for the few central lines that exhibit enough SNR.

To summarize this section, we have shown how to use the EKF to characterize the phase noise of down-converted frequency combs using the example of an EO comb. Both in simulation and on experimental data, we demonstrate that the EKF outperforms the conventional DSP method. Especially for very low SNR comb lines, the EKF can still recover the phase noise while the conventional method fails. While the EO comb was used as a benchmark system, the latent space EKF method is much more general. Two variants of the latent space were introduced, which, in principle, allow for its application on any comb.

## IV. DUAL-COMB SETUPS

In this section, we investigate the dual-comb setup. We continue to utilize the same experimental setup as displayed in Fig. 1. Both the CUT and the LO are chosen to be frequency combs. In this case, the detected signal has the most general form as derived in Eq. (4). Each optical comb line of the CUT can, in principle, create a beat note with any of the optical comb lines of the RS. The number of beat notes in the detected signal depends on the repetition frequencies $\omega repCUT$ and $\omega repLO$ of the comb, as well as the available detection bandwidth.

The beat notes in dual-comb setups can be categorized into different orders. One line of each optical comb is arbitrarily chosen to have comb line index zero. The electrical beat note created by their mixing is a zeroth order beat note and is assigned a line index from each generating line. Its frequency will, therefore, be denoted as $\omega 0,0=\omega 0LO\u2212\omega 0CUT$. All beat notes generated from the neighboring comb lines of each comb with matching index are also considered zeroth order. The index difference of the generating lines determines the order number of a beat note. A beat note with frequency $\omega 0,1=\omega 0LO\u2212\omega 1CUT$ would, therefore, be considered a beat note of order one.

In Fig. 12, the beat note generation process is illustrated. The spectrum of the zeroth-order beat notes is again a comb in the electrical domain with repetition frequency $\Delta \omega rep=\omega repCUT\u2212\omega repLO$.

*ω*

_{rep}to be very small compared to the optical comb repetition frequencies [Eq. (42)], the detection bandwidth can be reduced via low-pass filtering such that the first set of inequations (41) are fulfilled.

*f*

_{BW}and will not be part of the detected signal. The analytical form of the detected signal in Eq. (4) is reduced to

*l*as presented in Fig. 12, where the comb line corresponding to

*l*= 0 is arbitrarily chosen for each comb. The beat note frequencies are

### A. Experimental dual mode-locked laser

We will demonstrate the application of the EKF in dual-comb setups using an experimental measurement with two frequency-modulated quantum dot mode-locked lasers (MLL). A description of the fabrication and a general characterization of the MLLs used in the measurement are given in Ref. 48. Both MLLs were fabricated in the same way and should have similar noise properties. The distinction between CUT and LO is, therefore, not useful. We will, however, keep the naming convention to maintain the notation. Their repetition frequencies are $\omega repCUT\u2248$ 59.4 GHz and $\omega repLO\u2248$ 58.4 GHz. Both combs are centered around *λ* = 1298 nm with an average comb line power of −20 dBm. The beat notes have a spacing of Δ*ω*_{rep} ≈ 1 GHz. A single photodetector with a 15 GHz bandwidth was used to detect the signal. The resulting spectrum contains 16 regular and 12 aliased comb lines. With these parameters, the inequations (41) and (42) are fulfilled. The signal spectrum, therefore, only contains zeroth-order beat notes.

While we can use the same state-space model as for the EO comb, we have to decide whether the mixing matrix that defines the latent phase noise space should be constrained. In Sec. III, we used an optimization transformation function Λ [as defined in Eq. (34)], which constrains the latent phase noise sources. In particular, we chose to constrain **C**_{c.} to two latent noise sources: a common mode and a differential mode source. This model was chosen because, in the EO comb, a comb line index-dependent model of the phase noise is known and defined in Eq. (28).

For the MLLs used in this experiment, we do not have an explicit phase noise model available. However, we can use the so-called elastic tape model for phase noise to our advantage. It assumes that the physical noise sources that contribute to the comb line phase noise have a linear dependence on the comb line number.^{47} Each physical noise source may, however, differ in its fixed point, which is the comb line index at which its contribution is minimal. This model has been shown to describe the comb line phase noise dependency in MLLs as well as microring-resonator combs.^{49,50}

*φ*

_{n}(

*t*) is a physical noise source centered at the comb line with index

*l*

_{fix,n}. The total phase noise per comb line is then the sum of the independent physical noise source contributions. As the contribution of each source is linear with respect to the comb line index

*l*, the total phase noise dependence on line number must also be linear. If there are at least two physical noise sources with different fixed points, the total phase noise may also have a line number independent common mode contribution. This means that the total phase noise in each comb can be described by the decomposition of one common mode and one differential mode latent phase noise source. By the same logic, the phase noise in the electrical beat notes can also be described by only two noise modes.

This is the same decomposition we used in the EO comb case. We could, therefore, use the same two-dimensional latent space constraints as before. To, however, showcase the versatility of the method and the validity of the assumptions made earlier, we propose a hybrid latent space model. We choose the latent space dimension *L* = 3, where the first two noise sources are constrained to common and differential modes. Analogous to the EO comb, we constrain the first column of **C** to be constant with line index and the second to be linear, with an optimizable offset *l*_{center}. These two latent dimensions represent the common mode and differential mode noise sources. The third column of **C** is left unconstrained. If any additional noise sources are present that do not follow the elastic tape model, this additional latent space component should capture their behavior.

**C**drastically improves the robustness and stability of the optimization. As we are not interested in interpreting the latent noise sources in this case, changing the scale of

**C**is not a problem. The normalization is applied after the optimization transformation function as

**C**, where ‖·‖

_{2}is the 2-norm. Apart from this change, the same optimization transform function Λ as in Eq. (34) is used.

While we cannot interpret the latent space noise sources as physical quantities, we can still gain insights by investigating their properties. In Figs. 13(a) and 13(b), the column vectors of the optimized **C** matrix are shown, alongside the frequency noise PSDs of the latent sources. The two constrained vectors **C**^{(0)} and **C**^{(1)} look as expected, where the optimized offset of the zero-crossing of the differential mode vector occurs at $lcenterMAP=3$. The third free column of **C** exhibits a mostly smooth shape with smaller values at the central lines and larger values toward the outer lines. The frequency noise PSDs of the latent sources reveal the true scale of each noise mode. The common mode noise is the largest latent source and makes the dominant contribution to the phase noise per comb line. The differential mode is significantly smaller and has a mostly flat frequency noise PSD profile.

The free noise mode *ψ*_{3} has an approximate ∝*f*^{2} dependency on the PSD offset frequency. This makes the third latent noise source a white noise contribution to the comb line phase noise. It is very likely that some of the measurement noise present in the signal is interpreted as phase noise by the EKF. While this is generally undesirable, it also shows that during the optimization process, the EKF was not able to find phase noise contributions that violate the elastic tape model. This is not proof of the non-existence of higher-order noise terms, but it does show that the comb line phase noise can be accurately estimated using only the common and differential noise mode models.

Finally, Fig. 13(c) displays the intergrated phase noise PSD per comb line. The EKF estimate is compared to the conventional method. As we have seen before, both methods agree on the high SNR of the comb lines. For the outer low SNR lines, the conventional method returns a significantly larger variance estimate compared to the EKF.

### B. Digital noise compensation

We have shown how the EKF can be used to characterize the total phase noise of frequency combs. However, the EKF can also be used to digitally compensate for the phase and measurement noise of a multi-heterodyne signal.

As the distinction between the actual signal and perturbing measurement noise is inherent to the method, filtering out the measurement noise is straightforward. After optimizing the parameters ** θ** of the state-space model and applying the EKF, we are left with a time series of the latent phase noise estimates

*ψ*

_{m}(

*t*

_{k}). If we iterate over the samples

*k*and use

*ψ*

_{m}(

*t*

_{k}) as an input to the state-space model, we can, therefore, reconstruct the estimate of the EKF for the signal $y\u0303EKF(tk)$ without the measurement noise

*ξ*(

*t*

_{k}), such that $y(tk)=y\u0303EKF(tk)\u2212\xi (tk)$.

In principle, the phase noise of the detected comb lines can be compensated for in the same way. Using the latent space Eq. (31), we receive a phase noise time series *ϕ*_{l}(*t*_{k}) for each comb line. Applying the inverted phase noise −*ϕ*_{l}(*t*_{k}) as a phase shift to each line in $y\u0303EKF(tk)$ should result in a phase and measurement noise-free estimate of the original noisy signal *y*(*t*_{k}).

In a similar approach, it has been shown that the EKF can be used to recover comb lines in significantly more noisy conditions.^{28,51} In dual-comb spectroscopy experiments, it is common to produce very narrowly spaced multi-heterodyne signals to achieve good frequency resolution. If the combs have significant phase noise, the fluctuations of the comb line frequencies can cause overlap in the spectrum. Even in situations where individual lines were not discernible in the spectrum, using the EKF to compensate for the phase noise was possible.

Digital noise compensation methods can relax the requirements of frequency combs for metrology applications.^{52} Instead of relying on complex experimental setups to lock combs, the phase noise can be removed retroactively. While the presented method involves complex DSP and requires significant computational processing power, further improvements to its efficiency could enable the use of free-running combs in dual-comb spectroscopy applications.

## V. CONCLUSION

In this tutorial, we have shown how to characterize laser and frequency comb phase noise using multi-heterodyne measurements and the extended Kalman filter. By applying automatic differentiation and adaptive optimization methods, a state space model of the detected multi-heterodyne signal with an arbitrary shape can be learned from a measured signal. On simulated and experimental datasets in different experimental setups, we have demonstrated the method’s versatility and its advantages over conventional DSP approaches. Furthermore, we have discussed how the framework can be used for offline digital phase and measurement noise compensation. Wideband noise compensation in free-running dual-comb spectroscopy could reduce the hardware requirements on the combs and eventually lead to more efficient operation.

## ACKNOWLEDGMENTS

This work has been funded by the SPOC Center (Grant No. DNRF 123) and Villum Fonden (Grant No. VI-POPCOM 54486). The authors acknowledge Mario Dummont, Osama Terra, Bozhang Dong, and John E. Bowers from the University of California, Santa Barbara, USA, for providing the dual-comb MLL data.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**Jasper Riebesehl**: Conceptualization (equal); Methodology (equal); Software (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal). **Holger R. Heebøll**: Methodology (equal); Writing – review & editing (equal). **Aleksandr Razumov**: Writing – review & editing (equal). **Michael Galili**: Supervision (lead); Writing – review & editing (equal). **Darko Zibar**: Conceptualization (equal); Funding acquisition (equal); Methodology (equal); Supervision (equal); Writing – review & editing (equal).

## DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding author upon reasonable request.

### APPENDIX: SECOND-ORDER UNCERTAINTY PROPAGATION TO CALCULATE CONFIDENCE INTERVALS OF POWER SPECTRAL DENSITIES

In this section, we derive an expression for the confidence interval of a power spectral density given a sampled time series with time-dependent uncertainties. In our case, the expression is used to propagate the EKF uncertainty estimates of the phase noise time series to the phase noise PSD. This allows the quantification of the EKF’s uncertainty in the PSD picture. The approach is similar to a derivation presented in Ref. 53. Our derivation is, however, generalized for time-dependent uncertainties and exact due to higher-order error propagation.

#### 1. PSD estimation

*x*(

*t*

_{k}), Welch’s method is used.

^{54}

*t*

_{k}=

*k*Δ

*T*is the discretized time variable with sampling period Δ

*T*.

*x*(

*t*

_{k}) is divided into

*N*

_{seg.}potentially overlapping segments of length

*K*, and each segment is multiplied by a window function

*w*(

*k*) of the same length. Using

*l*corresponds to the Fourier frequency,

*f*

_{l}=

*l*/(

*K*Δ

*T*).

*X*(

*l*)

_{n}is the discrete Fourier transform of the n-th windowed segment, defined as

*S*

_{n}(

*l*), these periodograms are averaged to obtain the final PSD estimate using

#### 2. Uncertainty propagation

We assume we are given a time series *u*_{x}(*t*_{k}) of the uncertainty of *x*(*t*_{k}). In our case, *u*_{x}(*t*_{k}) is the standard deviation of the Gaussian distribution the EKF predicts.

*S*(

*l*),

*x*(

*t*

_{k}) first enters the calculation of the Fourier transform in Eq. (A2). To propagate the uncertainty, the real and imaginary parts of

*X*(

*l*) are treated separately.

^{55}As the Fourier transform is linear in

*x*(

*t*

_{k}), we can use linear error propagation,

*S*

_{n}(

*l*) exactly since the Eq. (A1a) is not linear. The expression is quadratic, which means that second-order error propagation is sufficient as all higher-order sensitivity coefficients vanish. Using the analytical form given in Ref. 56, we arrive at the expression

*κ*= 3, the kurtosis of a Gaussian distribution.

*S*(

*l*) can be calculated according to linear uncertainty propagation,

## REFERENCES

^{−17}fractional frequency instability through a 2220 km optical fibre network

*CLEO 2024*

*,*

*CLEO 2024*

*,*

*CLEO 2024, Technical Digest Series*

*Fiber-Optic Communication Systems*

*Bayesian Filtering and Smoothing*

*2020 IEEE 30th International Workshop on Machine Learning for Signal Processing (MLSP)*

*Advances in Nonlinear Dynamics and Control of Mechanical and Physical Systems*

*Springer Proceedings in Physics*