Frequency comb ptychoscopy

Multiheterodyne techniques using frequency combs—radiation sources whose lines are perfectly evenly-spaced—have revolutionized science. By beating sources with the many lines of a comb, their spectra are recovered. Even so, these approaches are fundamentally limited to probing coherent sources, such as lasers. They are unable to measure most spectra that occur in nature. Here we present frequency comb ptychoscopy, a technique that allows for the spectrum of any complex broadband source to be retrieved using a comb. In this approach, the spectrum is reconstructed by unfolding the simultaneous beating of a source with each comb line. We demonstrate this both theoretically and experimentally, at microwave frequencies. This approach can reconstruct the spectrum of nearly any complex source to high resolution, and the speed, resolution, and generality of this technique will allow chip-scale frequency combs to have an impact in a wide swath of new applications, such as remote sensing and passive spectral imaging.


SUPPLEMENTARY NOTE 1: DETAILED DERIVATION
To show that equation (2) of the main text holds, we will derive the result for the dual comb version first and extend it to the delayed comb version. All electric fields are expressed as a superposition of exponentials as where E ci (t) is the field of the ith comb and E s (t) is the field of the signal to be measured. For convenience the signal is represented as a summation rather than an integral. In a heterodyne measurement, the raw signal that is recorded is where we neglect the intracomb beat terms (which occur only at vanishingly-narrow multiples of the repetition rate) and the intrasignal beat terms (which are assumed to be small). We assume that our detector has a bandwidth larger than f r /2 (half the comb repetition rates) and that the signal is digitized with a sample period ∆t sufficiently small to avoid aliasing. In the spirit of spectrograms, time is divided into two separate time axes-the fast time t i and the slow time T j -such that the total time is given by t = t i + T j . The data is similarly divided into batches that are N t ∆t long (which determines the resolution bandwidth). In other words, T j can be taken as T j = (j − 1)N t ∆t, where j is an integer.

Dual comb version.
First, we calculate the signals in terms of the FFT IF frequencies ω k ≡ 2π Nt∆t k. For the FFT, we use the convention that F[f ](ω k ) ≡ 1 Nt i e −iω k ti f i . For comb 1, we find that the spectrogram F 1 (ω k , T j ) is given by Because our resolution bandwidth is determined by our batch length (RBW= 1 Nt∆t ), in order to proceed we make an approximation in which our signal frequencies are all an integer number of RBWs away from comb 1, i.e. ω (s) Nt∆t l for some integer l. As a result, where in the last line we used the definition of the FFT twice. Due to our finite resolution bandwidth approximation, the spectrogram would appear to be constant in T . However, once this approximation has been made it cannot be modified when computing the same quantity for comb 2. Performing a 1 → 2 substitution, one finds that Thus, while F 1 is stationary in time, F 2 is not, and in fact beats periodically at the dual comb frequencies ∆ n . Provided the data has been recorded long enough to resolve individual beat frequencies (i.e., that the data is at least recorded for 1/∆f r = 1/|f r2 − f r1 |), these beatings can be resolved. For simplicity, we assume that the number of batches N T has been chosen to ensure that N T N t ∆t = 1/∆f r , which will resolve exactly one dual comb beat tooth. Finally, we compute the Fourier transform of the correlation function and its expectation value: Because the number of batches was chosen to be an integer number of 1/∆f r , the summation 1 N T j e i(∆n−∆ l )Tj is a summation over roots of unity and vanishes unless n = l, leaving This result is general for any source. For the vast majority of sources we can make an additional assumption, which is that the long-term correlation between components of frequencies spaced by the repetition rate of the comb vanishes. This assumption is essentially valid for any source that is not a comb of the same spacing as either of the LO combs. With this additional assumption we can eliminate the cross terms, leaving By calculating every C n , the signal's power relative to a comb line can always be extracted.
Delayed comb version. For the delayed comb, the analysis is similar. The analysis for comb 1 is fully identical n ), and the analysis for comb 2 is found by setting E n τ and ∆ n = 0. However, a subtlety arises when delay is a linear function of lab time (i.e., a linear scan is performed rather than a step scan). Because delay changes during the batch, this has the same effect as Doppler shifting the IF frequencies in a manner similar to a nonzero ∆ n . One should therefore proceed as before, but using the explicit mapping τ → 2v c (t i + T j ). This results in By comparing this to the dual comb version of the same result, we find that −ω (c) n 2v c has replaced ∆ n . If we now define τ j ≡ 2v c T j and make the appropriate substitutions into our definition of C n , we find that we must instead calculate which results in Example: single-line reconstruction When only one line is present in the signal, vernier spectroscopy is sufficient for reconstruction. However, it is still instructive to see how ptychoscopy reconstructs such signals. Consider the dual comb case and consider a signal consisting of a single line, E s e iωst , which beats optical power  Figure 1a). Assuming once again that the signal is an integer multiple of the RBW away from the comb 1 line, equation (4) becomes In other words, each spectrogram is a delta function at the beat frequency. Because we have posed the problem discretely, here it is Kronecker. Additionally, they have a mutual phase difference that rotates at the dual-comb difference frequency. Note that for the signal to be real, we must strictly also consider a negative frequency contribution (E * s e −iωst ), which beats at negative frequencies (Supplementary Figure 1b). If we now compute the instantaneous correlation between the two spectograms, we find that As we shifted F 2 towards positive ∆ n , the negative frequency contributions will vanish-the two delta functions at negative frequencies are shifted apart (Supplementary Figure 1c), resolving the double sideband ambiguity. The comb line ambiguity is resolved by Fourier transforming the result in the T direction, obtaining provided that data has been recorded for an integer number of 1/∆f r 's. Thus, we have correctly determined that the signal spectrum is a delta function offset from the pth line by ω s − ω (c1) p (and only the pth line), and we have correctly determined that its magnitude is |E s | 2 .
Note that the prior analysis says nothing about the bandwidth of the detector. For dense combs, multiple pairs of comb lines can beat with the signal. In this case, each pair of combs will actually yield a competing version of the same data, albeit at different intermediate frequencies-a Rashomon effect. In the single line case shown above, the (p + 1)th line would yield a delta function at ω s − ω (c1) p − ω r (a large negative frequency). However, as detector responsivity falls at high IFs, the data corresponding to beating with the nearest comb line will typically be the most reliable.

SUPPLEMENTARY NOTE 2: SENSITIVITY
For astronomical applications, it is important to analyze the sensitivity of this approach. We do so for the case of additive white Gaussian noise (AWGN) and compare with the sensitivity of a tunable local oscillator. We assume that each detector measurement S i (t) is perturbed by a white noise source with variance σ 2 . This noise source is also taken to be uncorrelated between each detector. By propagating this noise through the reconstruction, one can show that the variance of the extracted signal in the dual comb version is given by As an example, in Supplementary Figure 2 we plot this expression for the data plotted in Figure 3. Because the summation over m is over both positive and negative frequencies, this means that the noise at a single frequency comes from the double-sided power spectral density of all signal components that share an IF with the frequency under consideration. The noise is frequency-dependent and resembles the sum of the two raw power spectral densities. This is practically speaking a dynamic range limitation: if one attempts to measure a weak signal that shares an IF with a much stronger signal, image noise generated by the stronger signal will swamp the weaker one. Whether this is tolerable depends on the source and on the application. For spectra consisting of relatively narrow lines it can be avoided by a proper choice of f r , for example.
The corresponding expression for using a single LO with electric field E 0 e iω0t + c.c. to measure double sideband intensity is Both expressions have a constant term independent of the signal-usually neglected since it is fourth-order in the noise-and both have a term that decreases linearly with the total measurement time. To compare these two expressions, we note that for a uniform comb the measurement time of the tunable LO must be N c times smaller to account for the fact that the comb measurement is multiplexed. In addition, the power per comb tooth for the comb version should be N c times smaller than the tunable LO's power if the mixer is to be optimally pumped. When there are M overlapping lines in the IF, the signal-to-noise ratio (SNR) for both the comb and the tunable LO are respectively given by In the best case scenario, where no lines overlap, the SNRs are identical. In the worst case scenario (a broadband light source), M = N c and the SNR is a factor of √ N c times worse. In addition to additive noise, there is multiplicative noise that arises from the fact that even incoherent signals can have transient frequency domain correlations.
Even when E s (ω)E * s (ω + nω r ) = 0, it will not be the case that |E s (ω)E * s (ω + nω r )| 2 = 0. This manifests as noise, and it effectively limits the dynamic range of the measurement. While the exact form depends on the details of the source-white frequency noise differs from flicker noise, for example-the variance is typically proportional to Supplementary Equation (20) for sufficiently broadband spectra.

SUPPLEMENTARY NOTE 3: PHASE CORRECTION
In the derivation of this approach we assumed that the comb lines being used to probe the signal were free of phase noise. In fact, this assumption can be relaxed if an additional dual comb spectroscopy measurement is performed to measure their mutual phase fluctuations and compensate for them. Suppose that the dual comb beating of pair n has been digitized and is given by n E (c1) * n e i∆nt , and that its magnitude has been divided out to construct p n (t) = e i(φ2n−φ1n) e i∆nt . The correlation function is then given by Since the total time is given by t = t i + T j , the explicit dependence on ∆ n can be removed by substituting in p n e i(φ1n−φ2n) and noting that the summation over i is merely an FFT over the t i axis: This expression is convenient since it eliminates any explicit references to ∆ n by premultiplying the signal before the spectrogram calculation. Not only does it remove phase noise, but it also makes calculation of the power spectrum more convenient, as it removes the global phase of the dual comb lines. By defining the alternative correlation functionĈ n (ω) aŝ we find that This result guarantees that phase noise does not contribute to the reconstruction, essentially by eliminating the phase entirely. Supplementary Figure 3 illustrates this for the data in Figure 3 for free-running combs. The combs have random walk phase noise producing offset fluctuations with a 1 MHz FWHM and repetition rate fluctuations with a 1 kHz FWHM, similar to what is found in free-running QCLs. As a result, the reconstruction produced by the previous approach fails. However, using the alternative correlation defined in (26), the correct result is recovered.

SUPPLEMENTARY NOTE 4: DETECTOR NONLINEARITY
Because we are using many frequencies to reconstruct our signal, it is possible that detector/mixer nonlinearity could negatively impact the results of the measurement, for example by distorting the signal. For example, this can be a challenge in dual comb spectroscopy when highly nonlinear mixers such as hot electron bolometers are used to measure the spectrum [1], as nonlinearity will cause different comb lines to mix. However, in this case we do not expect such effects to be significant. When operating in the heterodyne limit, the measured signal is given by and the |E ci | 2 term is much larger than the heterodyne term. Even so, it is easy to ignore because it only beats at multiples of the repetition rate (i.e., is periodic). Any nonlinearity will act upon it only to produce another signal that beats at the repetition rate, and we can therefore expect to continue to be able to ignore it. We therefore expect that the lone effect of nonlinearity on the heterodyne measurement is to reduce its sensitivity. For example, if a standard two-level saturation model is used to model the detector response, then the output signal will be related to the input power by where P sat is the saturation power. The heterodyne responsivity is essentially the differential response of this model, which is squared since both detectors suffer this decrease: Once again, we have verified this numerically. Supplementary Figure 4 shows the result of passing the narrowband data from before through the two-level saturation model, choosing the saturation power to coincide with the average power P sat = P (t) . As expected, there is essentially no effect on the reconstruction, except that all measured reconstructed power spectral densities have been reduced by a factor of 2 4 = 16. As with conventional heterodyne measurements, there will always be an optimal pumping power where the mixer is not yet saturated, but where the conversion gain is maximal. Supplementary Figure 5. Microwave experimental schematic. Microwave schematic outlining the experiment used to produce Figure 6. Unique components are labeled.
A schematic of the apparatus used to demonstrate broadband ptychoscopy is shown in Supplementary Figure 5. Two microwave frequency combs with bandwidths exceeding 4 GHz (High-Sierra Microwave) are high-pass filtered at 700 MHz to remove any low-frequency comb lines from the generated spectrum. Comb 1 has a tunable repetition rate that is set according to the desired difference frequency, while comb 2 has a fixed frequency of 5 or 25 MHz (depending on the measurement). The signal is then directed into two branches: the dual comb reference branch and the detection branch. The two branches utilize individual amplifiers since amplifiers are high directionality devices and reduce the amount of cross-talk between the detectors. In the reference branch, the two combs are summed together before being detected by a square-law detector (LTC5535ES6). This provides the reference spectrum that is used to correct the phase and amplitudes of the individual comb teeth. In the signal branch, each comb is amplified and summed with the signal to be measured and detected by a square-law detector. For simulating broadband continuous spectra, a frequency-modulated (FM) voltage-controlled oscillator (VCO) (CRYSTEK Microwave CRBV55BE) was used. To generate the FM voltage waveform a digital to analog converter (DAC) was used (Arduino Due). The DAC output was randomly generated from a designed distribution in order to ensure incoherence. Additionally, the tunable comb was set to 24.95 MHz and the 25 MHz fixed comb was used. By having a higher difference frequency a smoother spectrum is obtained. For our discrete spectrum the combs were chosen to be 4.997 MHz and 5 MHz, respectively. This was done in order to be able to achieve over 200 comb lines within a frequency octave.
Each microwave spectra displayed in the main text contains an "actual" and a "reconstructed" curve. The actual curve is generated by measuring the spectra via an RF spectrum analyzer (Agilent E4440A) at a sweep time of 1s and 100 averages. This allows us to treat this spectrum as the ground truth. The reconstructed curve was generated by measuring the voltage signal from the square-law detectors via a high-speed oscilloscope (Picoscope 6404D) and processing them according to the method derived in Supplementary Note 1. If our reconstruction method and the spectrum analyzer are set to have to have the same resolution, number of aver-ages, and measurement time, ptychoscopy outperforms the spectrum analyzer in terms of SNR, as the nonsimultaneous nature of a spectrum analyzer creates measurement artifacts that require many averages to eliminate.