## Abstract

Dual-comb spectroscopy allows for high-resolution spectra to be measured over broad bandwidths, but an essential requirement for coherent integration is the availability of a phase reference. Usually, this means that the combs’ phase and timing errors must be measured and either minimized by stabilization or removed by correction, limiting the technique’s applicability. We demonstrate that it is possible to extract the phase and timing signals of a multiheterodyne spectrum completely computationally, without any extra measurements or optical elements. These techniques are viable even when the relative linewidth exceeds the repetition rate difference and can tremendously simplify any dual-comb system. By reconceptualizing frequency combs in terms of the temporal structure of their phase noise, not their frequency stability, we can greatly expand the scope of multiheterodyne techniques.

- Frequency combs
- spectroscopy
- semiconductors
- computational
- quantum cascade lasers

## INTRODUCTION

Dual-comb spectroscopy is essentially the simplest possible version of multiheterodyne spectroscopy, and thus, they are usually considered synonymous. A frequency comb is a broadband coherent source that requires only two frequencies to fully describe it, the offset and the repetition rate (*1*). In dual-comb spectroscopy, two frequency combs are shined onto a common detector, and the heterodyne beating between different pairs of lines manifests at different radio frequencies (RF) (*2*–*8*). Although the idea is simple, the implementation is complicated by the fact that carrier-phase drift of the combs precludes coherent averaging (*5*, *7*, *9*). When this drift is known, its effect can be removed, but even measuring the absolute frequency of a comb line can be challenging. The most rigorous approach to measuring the carrier-envelope offset (CEO) directly is by means of the *f −* 2*f* technique, which requires that the comb be octave-spanning (*10*). Another approach is to beat the comb with a stable continuous wave (CW) laser (*11*–*13*). As long as the CW laser in question is moderately stable, this will provide a measurement of the comb’s relative CEO. Yet another approach is to use a narrowband optical filter, such as a fiber Bragg grating, to select only a portion of the combs’ optical spectra and to extract out the RF dual-comb beating of different portions (*14*). This approach requires no extra frequency references and allows for the extraction of the phase and timing fluctuations from multiheterodyne spectra, although it does require reasonably stable combs and narrowband optical filters.

Traditionally, optical frequency combs have been based on near-infrared mode-locked lasers, whose output consists of a uniform train of pulses. More recently, there has been a flurry of activity surrounding the development of novel chip-scale comb sources, including combs based on microresonators (*15*, *16*), combs based on interband cascade lasers, and combs based on quantum cascade lasers (QCLs), which can operate in the spectroscopically interesting mid-infrared (*17*, *18*) and terahertz (*19*–*21*) wavelengths. Dual-comb spectroscopy based on these lasers (*22*–*25*) is promising for enabling compact spectroscopic systems. However, performing phase and timing correction at these wavelengths is much more challenging. All of the previously discussed approaches are viable for combs based on near-infrared mode-locked lasers, which have a well-defined time-domain profile and operate in relatively technologically mature wavelengths. By contrast, long wavelengths have their own challenges. For example, every additional reference channel requires an additional optical detector, but the highest-sensitivity high-speed detectors are often cryogenically cooled (for example, mercury cadmium telluride in mid-infrared wavelengths and hot-electron bolometers in terahertz wavelengths). In addition, the lasers themselves are often cryogenically cooled, particularly at terahertz wavelengths, meaning that every extra CW laser and reference channel can greatly increase the size and complexity of the optical system.

To overcome these challenges, we show here that nearly all of the information needed to coherently correct a multiheterodyne spectrum is contained within the RF spectrum itself, excluding information about the average absolute frequency of the two combs. These techniques are broadly applicable, applying to all sorts of combs, and, as we will show, they apply even when the phase and timing errors are extremely large and when little clear time- or frequency-domain structure is present.

## RESULTS

For this demonstration, we consider the case of dual terahertz QCLs biased to a regime of marginal stability (nearing negative differential resistance). The lasers themselves are heterogeneous QCLs (*21*) that lase around 2.8 THz and are dispersion-compensated (*19*), whereas the detectors used are hot-electron bolometers and Schottky diode mixers. This paper is focused on coherent correction [for information on the terahertz spectroscopy performed with similar lasers, refer to the study of Yang *et al*. (*24*)].

Figure 1A shows a simplified experimental setup. Comb lasers A and B are shined onto a common detector, and the multiheterodyne beating between them is recorded. The beatnotes associated with each laser—the detected signal that contains the beatings between a comb’s adjacent lines—are shown in Fig. 1B, as is the recorded dual-comb signal. Ideally, because the beatnotes are separated by about Δ*f*_{A} − Δ*f*_{B} = 35 MHz, their mutual beating should lead to an RF comb with a repetition rate of 35 MHz. Unfortunately, because the lasers are operated at a bias in which they are only marginally stable, it is to some extent a stretch to call them combs. The beatnotes associated with these lasers are quite broad, and in addition, laser A has very clear sidebands (spaced by 140 kHz). As a consequence, the multiheterodyne signal obtained from these devices, shown in the time domain in Fig. 1B, is of poor quality and has none of the periodicity that one would expect from dual combs. Consequently, in the frequency domain, shown in Fig. 1C, the multiheterodyne signal is broad and has very little evidence of a downconverted dual-comb structure: Over 100 μs, the duration of the recording time, all features are completely washed out. One can imagine that this is simply an issue of long-term stability and that some structure might be obtained by processing the signal over shorter time intervals; spectra over 1 μs are shown in Fig. 1D. Although the comb structure is now evident on some spectra, it is not the case for all of them. There remain many instances in which phase instabilities completely spoil the spectrum, no matter how short the signal trace is cropped. As a result, no correction procedure that relies on interferogram cropping and alignment will succeed here (*26*): The signal must be corrected within the duration of an interferogram by the instantaneous phase and timing signals.

Extracting the phase and timing errors from the observed multiheterodyne signal is essentially a nonlinear estimation problem. Nevertheless, although we have no a priori knowledge of these errors, we have a model of what the RF comb should look like. Specifically, we expect it to take the form(1)where *y*(*t*) is the measured signal, is the dual-comb amplitude of the *n*th line, and φ_{0} and Δφ are the phase corresponding to the offset and repetition rate signals, and , respectively. In addition, the signal itself is corrupted by additive detector noise, and the parameters are all perturbed by multiplicative amplitude noise and phase noise. Although this estimation problem may seem intractable, it has been known for some time that if the measurement was a linear function of the parameters, it would be exactly solvable by the celebrated Kalman filter (*27*). In the case of a nonlinear measurement, one must linearize, resulting in an inexact solution. Nevertheless, good results can still be obtained. Essentially, we are fitting the measured multiheterodyne signal to the dual-comb model with the regularization constraint that the dual-comb amplitudes vary slowly (see Materials and Methods for more detailed information). With this approach, we can continuously update our estimates of the offset and repetition rates without any form of cropping; this, in principle, makes it very amenable to real-time processing (*11*). Because the data have been recorded in this case, we additionally perform Rauch-Tung-Striebel smoothing (*28*), which uses future knowledge to refine the state estimate and to correct for the group delay introduced by the standard filter. The resulting correction is superior to the correction obtained from the causal Kalman filter.

The physics of the comb enter primarily in the form of the multiplicative noise. Specifically, we assume that the comb complex amplitudes are perturbed only slightly at each time step (giving them a long time constant), whereas the phase and timing errors are perturbed much more (giving them a short time constant). In other words, we assume that the comb’s phase noise covariance is approximately rank −2. The Kalman filter quite naturally provides a way to test the validity of this assumption, because it makes a prediction about what the next measurement will be at every time step. By simply comparing the measured signal to the predicted signal, we can verify the efficacy of the prediction. For all of the data shown in this paper, the prediction residual is less than 8% of the signal power.

Figure 2 (A and B) shows the instantaneous repetition rate and offset frequency of the RF comb in Fig. 1. Several features are immediately apparent. The first is that both frequencies suffer a perturbation that recurs every 7 μs, which corresponds to the aforementioned 140-kHz sidebands evident in the beatnote of laser A. In other words, the beatnote undergoes a periodic instability that is imprinted onto the multiheterodyne spectrum. Second, the magnitude of the offset fluctuations (phase error) greatly exceeds the magnitude of the repetition rate fluctuations (timing error). This is not unexpected, because timing fluctuations correspond only to group index, whereas phase fluctuations also depend on phase index (*29*). Note also that the speed of the offset fluctuations—as much as 140 MHz per 220 ns during the perturbations, roughly 0.5Δ*f* per 1/Δ*f*—would be problematic for conventional techniques, because it implies that the short-term linewidth of a comb tooth greatly exceeds the spacing between comb teeth, a situation generally considered incompatible with dual-comb spectroscopy (*5*, *7*, *9*). Stabilization by thermal tuning is out of the question because the perturbation occurs on time scales that are too short, and even techniques that rely on measuring the beating with a CW laser would require an additional fast detector. Figure 2 (C and D) shows the time-domain multiheterodyne signals before and after the phase and timing correction, both during the instability (Fig. 2C) and away from it (Fig. 2D). During the instability, no clear periodicity or structure is obvious in the raw data; away from it, some periodicity is evident. In both cases, the signal predicted by the filter agrees very well with the actual data. As a result, following the phase and timing correction (discussed in Materials and Methods), the periodic comb structure is recovered.

Figure 3 shows the results of the computational correction in the frequency domain. The raw data from before are shown in Fig. 3A and, once again, show no comb structure. The phase-corrected data are shown in Fig. 3B; because phase correction removes the average offset frequency of the signal in addition to its fluctuations, the average offset 〈*f*_{0}〉 has been readded to correspond with the raw data. Phase correction reveals the individual multiheterodyne comb lines, although lines near the center of the comb are better-corrected than lines near the edge because of timing fluctuations that are still present. Finally, Fig. 3C shows the phase- and timing-corrected spectra, with insets showing zoomed views of several lines. All of the lines in the spectrum have been corrected, with full width at half maxima near the uncertainty limit of 10 kHz. Following the correction, some lines that were not apparent in the raw data have appeared out of the noise floor, such as the one shown in the rightmost inset. By filtering the data, one can verify that this is real signal (that is, not a computational artifact), but given detector dynamic range limitations, it may arise from detector nonlinearity rather than from heterodyne beating. Although the laser has a large disparity in mode amplitudes and large phase errors, it is possible to perform spectroscopy with these techniques (*24*).

## DISCUSSION

Computational phase and timing correction has advantages beyond its experimental elegance. As already shown, this approach can deal with extremely large phase-timing fluctuations. Figure S1 shows an even more extreme case in which the laser is biased in an even more unstable regime, causing the comb to chaotically switch between multiple operating conditions. Even here, correction remains possible. As long as the combs are coherent in the weak sense that the lines are evenly-spaced (*30*), with computational correction they become coherent in the strong sense that mutually coherent dual comb spectroscopy (*9*) can be performed. Although we focused here on unstable combs, it is also beneficial for stable combs and even combs operated in pulsed mode (*24*). Because the extended Kalman filter used here makes only weak assumptions about the evolution of the comb parameters (that is, a random walk with a Gaussian prior), almost no assumptions are needed to be made about the stationarity of the multiheterodyne signal. Figure S1 is an example of a dual-comb signal whose signal is nonstationary in three separate ways: the aforementioned instability at microsecond time scales, the chaotic switching between regimes at millisecond time scales, and a periodic fluctuation due to cryocooler fluctuations at 2 kHz.

In addition, computational correction offers very good performance even in the case of low-signal reference measurements. For example, fig. S2 shows a dual-channel correction based on the multiheterodyne signal from one of the channels, in this case, a Schottky mixer with a raw signal-to-noise ratio (SNR) of less than 25 dB. Even when the mixer’s noise is artificially boosted by 10 dB and practically no signal remains, computational correction remains informative on both channels. One may wonder whether computational correction offers any multiplex advantage analogous to the one present in conventional Fourier transform spectroscopy (*31*, *32*). For additive white noise, there is no advantage: The SNR of the extracted correction signal is identical between a reference comb and a reference CW laser of the same power. Even so, computational correction confers an advantage in the presence of excess phase noise, that is, phase noise that is not offset- or timing-related. Because it uses all of the lines of the comb, it is less susceptible to the excess noise of a particular line. We have observed that it is sometimes the case that even when most of the multiheterodyne spectrum is well-corrected, the weaker comb teeth can have sidebands (see fig. S3). These sidebands result from multiple lines sharing a mode, similar to what has been observed in microresonator combs (*33*), and manifest in QCL combs as a lack of coherence (*30*). If one were to use one of these lines to perform the correction instead of a clean comb line, the correction would be poor, using all of the comb’s lines averages out these effects.

Some caveats apply. Because the combs are not fully referenced, there remains an ambiguity in their average absolute frequency (*9*). That said, this is often unimportant for spectroscopy, because the free-running linewidth of most combs is narrower than typical spectroscopic features, and the absolute frequency can be localized using an etalon or an interferometer if need be. Another consideration is that because the filter needs to track the complex amplitudes of each line, this scheme is easier to apply to combs where the number of lines is not too large. This is easy to accomplish for chip-scale combs because the number of lines is typically in the hundreds but is computationally much more difficult for large mode-locked lasers [where the number of lines approaches 100,000 (*34*)]. One way this issue can be remedied is by intentionally tracking just a few of the multiheterodyne lines. For example, fig. S4 shows a comparison between a multiheterodyne signal corrected by a Kalman filter that incorrectly assumes that only two lines are present in the signal and one that correctly assumes that no more than 30 lines are present (*N* = 2 and *N* = 30, respectively). Although the corrected amplitudes are similar, because the 30-line model can use information from all of the signal’s lines, the resulting correction is more accurate. (An extra heuristic is needed to correct the two-line model; see “Removal of the coherent artifact” in Supplementary Materials and Methods for more information.)

In addition, because the linearized Kalman filter is suboptimal, the correction is not perfect and always leaves a noise pedestal. For the data shown in Fig. 3C, the pedestal is approximately 40 dB down from the peak of each line. Although we subsequently implemented nonlinear unscented Kalman filters (*35*), we found that the improvement was only marginal and that the result diverged more frequently. The overall robustness of the linearized filter in this case can be interpreted using an analogy to phase-locked loops. It has long been known that when a linearized Kalman filter is used to track a single oscillator, it is equivalent to an ideal phase-locked loop, one whose gain perfectly reflects the amount of noise in the system (*36*). Because a comb is simply a system of oscillators whose phase noise is strongly correlated, a similar interpretation can apply here. The linearized filter essentially constructs a correlated bank of ideal loops and adjusts the phase of each loop while ignoring signals resulting from other lines. The chief deficiency of the linearized filter is the possibility of mislocking, which is described in “Nonconvexity of the error function” in Supplementary Materials and Methods. Although nonconvexity is corrected heuristically here, it could be resolved without heuristics by using a more general tracking algorithm, such as a particle filter (*37*).

In conclusion, we have shown that coherent multiheterodyne spectroscopy can effectively be performed using a computationally enhanced method. Nearly all of the relevant phase and timing information traditionally acquired by separate reference channels is buried within the dual-comb spectrum and can be used to self-correct or to cross-correct a spectroscopic signal channel. These approaches are particularly useful for dual-comb systems based on semiconductor lasers or on microresonators, on account of their large free spectral range. Although we have used these techniques to correct comb spectra, they are very general and can apply to any coherent spectrum where a robust model of the phase-locking exists. The light source does not need to be comb-like, merely deterministic.

## MATERIALS AND METHODS

### Phase and timing correction

To perform the phase and timing correction, a general scheme that could take the large fluctuations present into account was needed. Phase correction was fairly simple, because phase fluctuations were common to all comb lines and manifested as a pure multiplication. Timing fluctuations presented a more substantial challenge, because the effect of timing fluctuations was to nonlinearly stretch the time axis. Defining an effective time , the multiheterodyne signal can be written as . The phase-corrected signal was found by calculating , and the phase-timing–corrected signal was found by using a nonuniform fast Fourier transform (*38*) to interpolate *y*_{0}(*t*) onto a linear grid, effectively calculating *y*_{0Δ}(*t*) ≡ *y*_{0}(τ^{−1}(*t*)). No other correction was performed.

### Optimization problem formulation

As mentioned in the text, the extended Kalman filter could be viewed as fitting the measured data to the model with a regularization constraint (*39*). This was an optimization problem that sought to minimize the error(2)where *x*_{k} is the state of our system at time *k*, *y*_{k} is the measurement at time *k*, *h*(*x*_{k}) is the measurement function of the state, *f*(*x*_{k}) represents how the state evolves, *R* is the measurement noise covariance, *Q* is the process noise covariance, and . The first term represents how closely the predicted measurement matches the observed measurement, and the second is a regularization term that takes the role of a time constant, in this case controlling how much multiplicative noise is present in the system. Equation 2 is extremely general, applying to any nonlinear system modeled by a Kalman filter, and all of the physics of the dual-comb signal are present in the form of the model *h*(*x*_{k}) and in the process noise covariance *Q*.

### Model of the dual-comb spectrum

As discussed in the text, the measurement model is a function of the complex amplitudes, phase signal, and timing signal(3)where we have separated the complex amplitudes into their amplitudes and phases. Alternative equivalent formulations existed (for example, considering the quadratures of the complex amplitudes), and they differed somewhat, on account of the nonlinearity of the filter. This particular formulation allowed us to place amplitude and phase noise into the system very naturally. If *N* is the number of lines under consideration, the state of the system is described by a vector length of 2*N* + 4, which contains the offset and repetition rate (*f*_{0} and Δ*f*), the corresponding phases (φ_{0} and Δφ), the mode amplitudes ({*r*_{n}}), and the mode phases ({φ_{n}}). At every time step, the two frequencies, modal amplitudes, and phases were assumed to be left unchanged—perturbed only by Brownian noise—whereas the offset and timing phases were updated by the frequencies(4)

Similarly, the process noise covariance, *Q*, was constructed in a way that was consistent with the spirit of a comb. It contained relatively large amounts of phase and timing noise, relatively small amounts of multiplicative amplitude noise, and relatively small amounts of additional phase noise. Note that the additional phase noise was constructed in such a way that it did not contribute any extra phase/timing error, that is, it has rank *N* – 2. The rest of the processing involved a fairly standard Kalman filter, with some modifications made to account for the particulars of dual-comb spectroscopy. (Both extended and unscented filters were used, with very little performance difference in this case.) These are discussed in Supplementary Materials and Methods.

## SUPPLEMENTARY MATERIALS

Supplementary material for this article is available at http://advances.sciencemag.org/cgi/content/full/2/11/e1601227/DC1

Supplementary Materials and Methods

fig. S1. Dual-comb spectroscopy in a chaotic regime.

fig. S2. Cross-correction of multiple spectra.

fig. S3. Incoherent sidebands in a weak dual-comb tooth.

fig. S4. Multiheterodyne spectrum corrected by a two-line model.

fig. S5. Correction of artificial dual-comb data.

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is **not** for commercial advantage and provided the original work is properly cited.

## REFERENCES AND NOTES

**Funding:**The work was supported by the Defense Advanced Research Projects Agency Spectral Combs from UV to THz program through grant number W31P4Q-16-1-0001 from the Aviation and Missile Research, Development, and Engineering Center and NSF.

**Author contributions:**D.B. developed the formalism and performed the analysis. Y.Y. fabricated the devices and performed the measurements. Q.H. oversaw the project. All authors wrote the manuscript.

**Competing interests:**D.B., Y.Y., and Q.H. has a U.S. patent, 15/259,687, “Computationally assisted multi-heterodyne spectroscopy,” filed 8 September 2016. The remaining authors declare that they have no competing interests.

**Data and materials availability:**All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from the authors.

- Copyright © 2016, The Authors