Frequency Noise Characterization of Narrow-Linewidth Semiconductor Lasers: A Bayesian Approach

We describe a Bayesian estimation approach to infer on the frequency noise characteristics of narrow-linewidth semiconductor lasers from delayed self-heterodyne beat note measurements. Our technique is grounded in a statistical model of the measurement process that accounts for both the impact of the interferometer and the detector noise. The approach yields accurate results, even in scenarios where the intrinsic linewidth plateau is obscured by detector noise. The analysis is carried out using a Markov-chain Monte Carlo method in the frequency domain and exploits prior knowledge about the statistical distribution of the data. The method is validated using simulated time series data from a stochastic laser rate equation model incorporating $1/f$-type non-Markovian noise.


I. INTRODUCTION
N ARROW linewidth lasers are key components of many applications in modern science and technology, ranging from precision metrology, e.g., optical atomic clocks [1] and gravitational wave interferometers [2], to coherent optical communication systems [3] and quantum computers based on trapped ions [4] or neutral atoms [5].To perform well in these technologies, the laser requires a high degree of spectral coherence, i.e., low phase or frequency noise corresponding to a narrow spectral linewidth.The laser linewidth quantifies the width of the optical power spectrum, which consists of two parts: The intrinsic (Lorentzian) linewidth, which is dominated by the rate of spontaneous emission into the laser mode, that is typically broadened by an additional 1/f -like technical noise component (also flicker noise) [6], [7], [8].As the width of the optical power spectrum is dominantly determined by frequency noise, the corresponding frequency noise power spectral density (FN-PSD) provides an almost complete characterization of the spectral quality of the laser, where the effects of non-Markovian noise can be well separated from that of white noise.The white noise component manifests itself as a plateau in the high-frequency The authors are with the Weierstrass Institute for Applied Analysis and Stochastics (WIAS), 10117 Berlin, Germany (e-mail: mertenskoetter@ wias-berlin.de;kantner@wias-berlin.de).
Digital Object Identifier 10.1109/JPHOT.2024.3385184part of the FN-PSD, from which the intrinsic laser linewidth [9], [10] (that is of major interest for most of the aforementioned applications) can be deduced.In ultra-narrow linewidth lasers, the determination of the white noise plateau can be challenging, since it often sets in only at very high frequencies and may be obscured by 1/f -type noise (at low frequencies) or by detector noise (at high frequencies).
The standard technique for the experimental measurement of the laser linewidth is the delayed self-heterodyne (DSH) method [11], [12], [13], [14], [15], which involves measuring the beat note signal between the optical field with a delayed and frequency-shifted copy of itself.The method is attractive because it does not require any external frequency standards or active frequency stabilization.A drawback of the technique, however, is the need for advanced data analysis, as both the footprint of the interferometer as well as the detector noise must be removed in order to reconstruct the FN-PSD [16].
In this paper, we present a Bayesian estimation approach to infer on the laser's FN-PSD from time series data.Our method is based on statistical modeling of the DSH measurement process and allows to extract accurate estimates of key parameters such as the intrinsic linewidth, even when the white noise plateau is obscured by detector noise.Other methods do not fully exploit the information inherent in the data, as they often require approximations neglecting 1/f noise, detector noise or fluctuations of the reference frequency [17].Our method is demonstrated for simulated time series data based on a stochastic rate equation model for a single-mode semiconductor laser including non-Markovian noise.The estimation results are critically discussed and assessed by comparison with analytical reference values and input parameters of the numerical simulation.

A. Experimental Setup
The DSH method is a standard technique for experimental characterization of the laser's FN-PSD [13].It involves splitting of the laser beam into two paths using an acousto-optic modulator (AOM), where one of the two beams is also frequencyshifted, see Fig. 1.The frequency-shifted beam is delayed using a long fiber before being finally superimposed with the other beam at a photo detector.The optical field E(t) = Re(E(t)) received by the detector reads E(t) ∝ P (t) e i(ω 0 t+φ(t)) where where is the phase jitter, η is the detector efficiency and ΔΩ is the harmonic difference frequency after beating of the laser and digital downshifting in the spectrum analyzer.The detector noise is assumed to be Gaussian white noise that is not correlated with fluctuations of the optical field.
The measured phase jitter, from which the fluctuations of the instantaneous frequency can be inferred, is obtained from the I and Q data as The detector noise in (2) affects the measured phase jitter and leads to the effective phase measurement noise ξ φ (t) that is approximately white [16] where σ = σ det /(ηP ).This white noise approximation requires that the average optical power P is much larger than the detector noise level σ det .

B. Periodogram and Power Spectral Density
The power spectral density (PSD) S y,y (ω) of a stationary stochastic process y(t) is given by the Fourier transform of its auto-correlation function C y,y (τ ) = y(t)y(t + τ ) (Wiener-Khinchin theorem) [18].Given only a sample of the trajectory (in discrete time), the PSD is typically estimated from the periodogram [19], which is given by the absolute square of the (discrete) Fourier transform of the time series The PSD then follows as the expectation value of the periodogram (ensemble mean) [19] S y,y (ω) = S y,y (ω) .
In the following, we use the notation where x, the hidden variable, is the instantaneous laser frequency and z, the observed variable, is the noisy measurement of the time-differentiated phase jitter.We work with frequencies rather than phases, as the phase itself is not a stationary process.Our goal is to estimate the FN-PSD of the free-running laser from the noisy DSH-measurements.Note that the FN-PSD S x,x (ω) is not directly observed in the experiment, but must be reconstructed from the measurements z, which can be deduced from the slow beat note of the intensity time series via (4).
The measured signal of the DSH experiment (in terms of frequency fluctuations) is described by where z(t) is the observed time series, h(t) is a convolution kernel given by the transfer function of the interferometer [16] and x(t) is the hidden time series of the instantaneous frequency fluctuations of the laser.Finally, ξ(t) is (colored) additive measurement noise (not correlated with the hidden signal).A sample simulated time series of frequency fluctuations x(t), which exhibits characteristic frequency drifts induced by non-Markovian noise as commonly observed in semiconductor lasers, is shown in Fig. 2.
Our goal is to characterize the statistical properties of the time series x(t).Fourier transform of (9) yields a relation between the PSDs of the observed and the hidden signal where the Fourier transformed transfer function reads Since z(t) follows from ( 4) by differentiation with respect to time, and thus ξ(t) = ξφ (t), the PSD of ξ(t) is related to that of ξ φ (t) by Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.as time-differentiation is equivalent with multiplication with iω in the frequency domain.The model for the FN-PSD can be derived from the stochastic laser rate ( 19) via linearresponse theory (see Appendix B) and includes both 1/f -type noise (described by C and ν) and a white noise plateau, where S ∞ quantifies the intrinsic linewidth.

C. Computation of the Periodograms
Care has to be taken in the numerical computation of the periodograms for two reasons.First, the PSD only exists for stationary processes.The optical phase, however, is not stationary.Even without the colored noise contributions, it performs to leading order a random walk (i.e., a Brownian motion), cf.(19c), so that its variance is increases over time.Computation of periodograms via (12) is thus not possible.While phase-noise PSDs are often reported in literature, it should be noted that these are necessarily given only for a certain measurement time and are not PSDs in the strict sense [20].The PSD of the differentiated process, i.e., the instantaneous frequency, does however exist.The stationarity of the frequency fluctuations is preserved even in the presence of 1/f -type colored noise in the phase (which is inherently non-stationary).Consequently, the FN-PSD must be computed from the time-differentiated time series of the optical phase as opposed to first calculating the phase-noise PSD and then multiplying it by ω 2 .Due to the non-stationarity of the optical phase fluctuations, both approaches are not equivalent, which is a common source of error.
Here, a second issue arises however.When computing the periodogram of differentiated white noise this way, one does not recover the expected ω 2 -dependency over the whole frequency range, but rather obtains a plateau at very low frequencies.As the detector noise is essentially differentiated white noise in the frequency time-series, we hence over-estimate its contribution to S z,z (ω) at low frequencies.This is only an issue very close to the first few roots of the transfer function where these inaccuracies dominate the calculated value of S z,z (ω).Hence we exclude points, where S z,z (ω) < 100 Hz from our analysis, which amounts to about 0.01% of the data.Typically, these kinds of numerical effects in spectral estimation are mitigated via windowing techniques, which, however, do not preserve the amplitude of the periodogram and thereby might shift the intrinsic linewidth plateau.

III. PARAMETER INFERENCE
In order to estimate the parameters characterizing the FN-PSD of the laser, we perform a Bayesian regression on the PSD S z,z (ω) of the detected signal using the transfer function (11) and the model PSDs ( 12) and ( 13).In the regression procedure, we regard the periodogram S z,z (ω) (available at a set of discrete frequencies) as a random variable, whose mean is given by (10), i.e, S z,z (ω) is an unbiased estimator of S z,z (ω).In order to conduct a maximum likelihood estimation using frequency-domain data, prior knowledge about the expected statistical distribution of the periodogram S z,z (ω) is required.
The frequency fluctuations x(t) and the measurement noise ξ(t) are assumed to be normally distributed, which is in excellent agreement with experimental data.The detected time series z(t) observed at discrete instances of time is thus a multi-variate Gaussian characterized by its covariance matrix.By random variable transformation, the periodogram of the measured time series is found to be exponentially distributed where the probability distribution function (PDF) is characterized by a parameter λ that is itself a function of the frequency ω and the unknown parameters θ = (C, ν, S ∞ , σ) T .We identify the parameter function λ(ω, θ) with the inverse mean of the periodogram data given by ( 10) such that which has a highly nonlinear frequency dependency.Fig. 3 shows a double-logarithmic plot of this PSD as well as the PSD (13) along with the respective periodograms of sampled data.The likelihood of observing a certain realization of the periodogram S z,z (ω k ) (at a discrete set of frequencies ω k ) given a set of parameters θ is then given by the likelihood function where the PDF p is given by (15).Note that here the underlying joint probability distribution factorizes completely, as the periodogram data S z,z (ω) at different frequencies are statistically independent.Bayesian regression on the parameters θ is performed by maximizing the likelihood function (17) using a Markov chain Monte Carlo (MCMC) approach [21].We employ the Metropolis-Hastings algorithm to sample a Markov chain of parameter Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.sets θ (j−1) → θ (j) → . . .that is distributed according to L(θ).The maximum of this distribution is then located at the parameter set θ * which is most likely to underlie the observed data.
We would like to stress that in the present scenario the regression on spectral data (i.e., on the periodogram S z,z (ω)) as described above has several advantages compared to the direct estimation from time-domain data using a Kalman filter [22], [23].The reason for this is that the construction of the likelihood function is either conceptually challenging due to the long interferometer delay in the equation for the observation process (9) or would require a high-dimensional Markovian embedding.Furthermore, the system exhibits long-range time-correlations due to the presence of 1/f -type noise, which would require very long time series to be considered in the regression and additional overhead in state augmentation.These issues render the estimation procedure in the time domain computationally expensive, but are easily overcome by transforming the problem to the frequency domain.I.

A. Markov Chain Monte Carlo Method
We demonstrate the estimation method described above for simulated data using the stochastic laser rate equations (19) given in Appendix A. The DSH measurement is simulated as described in Ref. [16] with simulation code that is available online [24].
Fig. 4 shows histograms of the sampled Markov chains θ (j) , which are estimators of the marginal distributions of the joint PDF of the parameters that is proportional to the likelihood function L(θ).A notable feature of the MCMC method is that it provides not only an estimate of the most probable set of parameters but rather their entire distribution, allowing to asses the uncertainty of the estimates, see Table I.We observe that the analytic reference values are in good agreement with the results of the estimation procedure.The analytic reference values for C and S ∞ are derived from linear response theory, see Appendix B, and the reference value for the effective noise amplitude σ is Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
obtained from an approximation of the actual detector noise, see [16].It should be noted that the present method allows for a direct observation of the detector noise over a wide frequency range (at the roots of the transfer function) without the need for an additional dedicated noise measurement.The reference value for ν is exact, as it is set as an input parameter.In addition to the histogram, the plot also indicates the position of the maximum likelihood estimator (MLE), which is the parameter set θ * that maximizes the likelihood function (17).
The MLE values partially differ significantly from the means of the marginal distributions calculated via MCMC.The reason for this is that the different parameters are strongly correlated, such that the joint distribution of the parameters does not factorize into the marginal distributions shown in Fig. 4. The mutual correlation of all parameter estimates is quantified by the matrix of Pearson correlation coefficients which is obtained as We observe that the estimate of the intrinsic linewidth is negatively correlated with the detector noise estimate σ, which is expected from (10).
The MCMC routine was run with normally distributed proposal functions for the parameters ν, S ∞ and σ and with a log-normal proposal distribution for C. The latter automatically enforces C > 0 and provides an efficient sampling as it matches the shape of the target distribution.With the variances of the proposal distribution tuned appropriately, we achieved an acceptance rate of the proposed parameter samples of around 45%.

B. Comparison With Welch's Method
The MCMC estimation results are benchmarked against Welch's method [25], which is widely used in PSD estimation.Here we focus on estimation of the intrinsic linewidth parameter S ∞ , which must be extracted from the flat section of the FN-PSD S x,x (ω).This approach thus necessitates, first, a clearly visible white noise plateau, that is not obscured by 1/f -noise or detector noise and, second, a reconstruction of the FN-PSD from the available S z,z (ω) periodogram.Although it is possible to compute an artifact-free reconstruction of S x,x (ω) (i.e., the periodogram of the time series of the hidden signal) from S z,z (ω) using parametric Wiener filters [16], much simpler methods are often used in practice.Most commonly, one disregards the detector noise in (10) such that a simple approximation is obtained by merely inverting the transfer function: It is well known that this approach features a series of reconstruction artifacts (at the roots of the transfer function) and systematically overestimates the signal if the detector noise is not negligible [16].The key step in Welch's method is partitioning of the time series z(t) into multiple segments, for each of which the periodogram is then computed individually and finally averaged.This results in an estimate of the FN-PSD S x,x (ω) with reduced variance, from which eventually the linewidth can be estimated.Due to the reconstruction artifacts, however, only small frequency domains inbetween the roots of the transfer function can be used for estimation.
For the simulated data set considered here, we obtained S Welch ∞ = (565 ± 6) Hz using the approach described above.The analysis was limited to data at frequencies far from the roots of the transfer function characterized by |H(2πf )| 2 ≥ 3.9 within the frequency window f ∈ [50, 80] MHz.While the value captures the correct order of magnitude of the reference value, it is less accurate than the MCMC result.In fact, the value is systematically too large (about twice the reference value) because the data were not corrected for detector noise that is not negligible in this realization, cf.Fig. 3.We see a significantly better agreement between the two methods for vanishing detector noise.

V. CONCLUSION
Bayesian inference methods allow for an accurate extraction of the parameters characterizing the FN-PSD (including the intrinsic laser linewidth) along with its uncertainties from DSH measurement data.Based on statistical modeling of the underlying measurement process, the method enables a reliable characterization of the spectral coherence properties of the laser even when the intrinsic linewidth plateau is obscured by detector noise.By incorporating explicit assumptions about the statistical distribution of the measured data, the method extracts the information encoded in the available data very efficiently, so that no averaging over many realizations or long times is required.The method has been demonstrated for simulated data, where the estimation results were found to be in excellent agreement with analytical reference values.

APPENDIX A STOCHASTIC LASER RATE EQUATIONS
The evolution of the photon number P , the optical phase φ and the carrier number N in the active region of a generic singlemode semiconductor laser is described by a set of Langevin equations modeling the fluctuation dynamics: Here, γ is the optical loss rate, P th is the thermal photon number (Bose-Einstein factor), Γ is the optical confinement factor, v g is the group velocity, Ω 0 is a detuning from the nominal CW frequency, α H is the linewidth enhancement factor (Henry factor), η inj is the injection efficiency, I is the pump current and q is the elementary charge.The net-gain is modeled as where g 0 is the gain coefficient, ε is the inverse saturation photon number (modeling gain compression) and N tr is the carrier number at transparency.Following [10], the rate of spontaneous emission into the lasing mode is modeled as which does not require any additional parameters and shows the correct asymptotics at low and high carrier numbers.The rate of stimulated absorption entering the noise amplitudes follows as g abs (P, N ) = g sp (P, N ) − g(P, N ).Finally, non-radiative recombination and spontaneous emission into waste modes is described by where V is the active region volume.We refer to Table II for a list of parameter values used in the simulation.The model (19) include several zero-mean Gaussian white noise F i (t), i ∈ {P, N, φ}, and Gaussian colored noise F i (t) sources as driving sources of the fluctuation dynamics.At a stationary point (P , N ), the white noise is characterized by the frequency-domain correlation functions The diffusion matrix has been obtained from the Itô-type stochastic differential equation model presented in [16] using algebraic conditions for the noise-free steady state and the high-power limit P 1 P th .The colored noise sources F i (t) phenomenologically describe the technical flicker noise and obey the frequency-domain correlation functions where the matrix elements have a 1/ω ν -type frequency dependence The colored noise is numerically generated using the Davies-Harte method [26].White and colored noise sources are assumed to be uncorrelated.

APPENDIX B LINEAR RESPONSE THEORY
Linearization of the system (19) at the noise-free steady state X = (P , N, 0) yields a linear relation for the response of the system δX = (δP, δN, δφ) to the noise, which is expressed in the frequency domain as −iωδ X (ω) = Jδ X (ω) + F (ω) + F (ω) .
Here, J = J(X) is the Jacobian of the deterministic part of the evolution equations (19) at the steady state.The matrix of noise power spectral densities is obtained via where K(ω) = (−J − iωI) −1 is the transfer matrix.We are interested in the FN-PSD S φ, φ(ω), which follows from the phase noise PSD as S φ, φ(ω) = ω 2 S φ,φ (ω).Note that the structure of the model PSD (13) can already be recognized from (20), which is a sum of a white noise and a colored noise part.Employing the usual approximations for semiconductor lasers [27], [28], the phenomenological model of the FN-PSD (13) can be derived from the present model in the limit of low frequencies far below the relaxation oscillation frequency.The analytical results for the parameters are S ∞ ≈ Γv g g sp P , N where (21a) is dominated by the colored noise in the optical field and (21b) is essentially the (modified) Schawlow-Townes formula for the intrinsic laser linewidth.

Manuscript received 14
March 2024; accepted 29 March 2024.Date of publication 4 April 2024; date of current version 18 April 2024.This work was supported by the German Research Foundation (DFG) through Germany's Excellence Strategy -The Berlin Mathematics Research Center MATH+ (EXC-2046/1) under Project 390685689.(Corresponding author: Lutz Mertenskötter.)

Fig. 1 .
Fig. 1.Experimental setup of the DSH method, where an AOM separates the incoming laser beam into two beams.One part of the signal is delayed by a long fiber and frequency shifted.The two beams are superimposed at a photo detector, which captures only the slow beat note signal.

Fig. 2 .
Fig.2.Simulated time series of frequency fluctuations ω(t) = φ(t) around the nominal continuous wave frequency ω 0 .The instantaneous laser frequency exhibits characteristic drifts, which are commonly observed in semiconductor lasers.The time series has been simulated using the Langevin model(19), where frequency drifts are induced by the 1/f -type colored noise contributions.The trajectory in the plot shows a moving average (over 50 ns) to improve the visibility of the effect.

Fig. 3 .
Fig. 3. (a) Periodogram S z,z (ω) of the observed time series z(t) along with the PSD S z,z (ω) estimated using the MCMC method.The signal features sharp dropouts at the roots of the transfer function (11) at frequencies f n = n/τ d , n ∈ Z.The detector noise PSD S ξ,ξ (ω) has a quadratic frequency dependency, see (12).(b) Estimated FN-PSD S x,x (ω) and periodogram S x,x (ω) of the hidden time series x(t).Both time series stem from the same numerical simulation of the DSH experiment.

Fig. 4 .
Fig. 4. Histograms of estimated parameters characterizing the FN-PSD obtained using the MCMC method.The sampled distributions are shown along with normally (Gaussian) and log-normally distributed probability density distributions.The blue line indicates the position of the maximum likelihood estimate (MLE) and the red line shows the analytical reference value.The estimation results are summarized in TableI.

TABLE II PARAMETER
VALUES USED IN STOCHASTIC TIME SERIES SIMULATION