Ultra-precise determination of thicknesses and refractive indices of optically thick dispersive materials by dual-comb spectroscopy

Precise measurements of the geometrical thickness of a sample and its refractive index are important for materials science, engineering, and medical diagnosis. Among the possible non-contact evaluation methods, optical interferometric techniques possess the potential of providing superior resolution. However, in the optical frequency region, the ambiguity in the absolute phase-shift makes it difficult to measure these parameters of optically thick dispersive materials with sufficient resolution. Here, we demonstrate that dual frequency-comb spectroscopy can be used to precisely determine the absolute sample-induced phase-shift by analyzing the data smoothness. This method enables simultaneous determination of the geometrical thickness and the refractive index of a planar sample with a precision of five and a half digits and an ultra-wide dynamic range. The thickness and the refractive index at 193.414 THz of a silicon wafer determined by this method are 0.52047(3) mm and 3.4756(3), respectively, without any prior knowledge of the refractive index.

The geometrical thicknesses d and the complex refractive indices of layers (n(ω)+ ik(ω), where ω is the optical angular frequency) are fundamental parameters for materials science, and thus ultraprecise determination of these parameters in a wide range of materials is important. Due to advances in optical technologies, several techniques can be used to implement precision metrology 1,2,3,4,5,6 . For example, laser interferometry has been used for the dimensional metrology of optically thin samples (whose thickness is comparable with the wavelength of light, λ) 7 . In this method, d and n(ω) of a planar sample can be determined by analyzing the absolute sample-induced phase change ΔΦabs(ω).
However, in conventional interferometry of optically thick samples (i.e. d>λ), we need to consider the ambiguity in the value of ΔΦabs(ω) 8 , because in this case, ΔΦabs(ω) is usually not identical to the measured phase change Δϕexp(ω) but can differ by multiples of 2π: where M(ω) is an (unknown) integer. Because M(ω) can be used to evaluate d and n(ω) of a planar sample, various interferometric methods have been proposed to precisely determine M(ω) 8 , for instance interferometric methods that employ rotating 9,10,11 or moving 12,13 the sample. However, since the mechanical system that is required for such sample translations, causes additional measurement uncertainties, an ultra-precise determination of d and n(ω) by these methods is still difficult. Moreover, optical ellipsometry is usually only applicable in the case of optically thin samples 14,15,16,17 , and low-coherence interferometry can only be used to measure the group refractive index (ng(ω)) of the sample 18,19,20,21 . Thus, a more sophisticated optical technique is required.
The development of the optical frequency comb (OFC) technology has revolutionized the field of precision metrology. Dual-comb spectroscopy (DCS), which uses two OFCs 22,23,24,25,26,27 , enables us to unambiguously determine ΔΦabs(ω) in open-air ranging applications 28,29,30,31 by combining timeof-flight and interferometric measurements, resulting in an ultra-wide dynamic range of the determination of the absolute distance. This concept has been applied to materials science in order to determine d and ng(ω) of samples like wafers 19,20,32,33 . However, because the distortion of the optical wave-packet obstructs the determination of ΔΦabs(ω) due to the dispersion of n(ω), the same precision as in open-air ranging has not been achieved.
Here, we demonstrate the ultra-precise determination of d, n(ω) of an optically thick dispersive material with an ultra-wide dynamic range. We determine ΔΦabs by analyzing the smoothness of the n(ω) and k(ω) spectra derived from the experimental data (i.e. the spectra should exhibit no artificially induced oscillating behavior with respect to frequency). We show that d and n(ω) of a silicon wafer can be simultaneously determined with a precision of five and a half digits: d=520473.7±1.5 nm and n=3.47563±0.00001 at 193.414 THz with a relative uncertainty of ~3×10 -6 . This precision of n(ω) is the best value that has so far been reported for a planar sample, and is comparable with the precision of the interferometric measurement of prism-shaped samples 12 .
Moreover, the value of n(ω) is consistent with the literature value within the previously reported uncertainty.

Results
The criterion for the correct value of M. Figure 1 describes the underlying problem in the DCS measurement; the phase ambiguity of the experimental Δϕexp(ω) data (Fig. 1a) and the resulting possible candidates for ΔΦabs(ω) (Fig. 1b).
Here, we explain how to solve this problem of ambiguity in ΔΦabs(ω) obtained by DCS. As a general approach, we first construct a continuous Δϕ(ω) curve ( Here, the subscript i identifies the discrete frequency values of the FT data in the range from the lowest to the highest frequency in the observable spectral region. The correct M and d values are those that minimize En(M,d)+Ek(M,d) (Fig. 2d).
Note that a similar analytical procedure has been proposed to determine d in the terahertz frequency range 34,35 , and this idea has been applied to the characterization of various materials 36,37,38,39 .
However, in the terahertz frequency range, the value of M is uniquely and easily determined by extrapolating Δϕ(ω) to ω=0, where ΔΦabs should be always zero. This distinct difference between the determination of M and d at terahertz and at optical frequencies makes the determination of N(ω) in the optical frequency range very difficult.
Measurement system. We mounted the sample on a mirror mount and fixed this holder on a motorized translation stage.
By moving the translation stage, we can easily obtain the interferograms with and without the sample.
The direction of the optical beam was carefully aligned to achieve normal incidence on the sample surface within an uncertainty of ≈5 mrad. All measurements were performed in ambient atmosphere at room temperature (23±2 °C).  spectra. Then, we averaged T(ω) and Δϕ(ω) over each subset for a selected Nave to determine the mean and the standard deviation of M and d (for example, Nave=1375 allows us to determine the standard deviation using four data points for M). Figure 6a and 6b show the Nave dependences of M and d, respectively. As shown in Fig. 6a, the averaged M values are always closest to 831, independent of Nave. In addition, the uncertainty in M is smaller than ±1 for Nave=1100 and the uncertainty in M vanishes for Nave=1375 because four determined M have the same value of 831.

Analysis of the interferogram data.
This means that we can determine M uniquely when Nave ≥ 1375, and thus we can precisely determine ΔΦabs(ω). In Fig. 6b, the mean value of d is 520.473-520.474 µm and its error bar is less than 5 nm for all considered Nave values. The error bar gradually becomes smaller for larger Nave values, and it becomes less than 1.5 nm for Nave≥1100. Our obtained value is consistent with the nominal value of the sample (525±25 µm), supporting the validity of our method.

The complex refractive index of silicon
Since d in our experiment was determined precisely when Nave=1833, we consider the n(ω) and k(ω) calculated from the three pairs of T(ω) and Δϕ(ω) obtained for Nave=1833. The averaged n(ω) and k(ω) are plotted in Fig. 6c and 6d, respectively. It is found that both n(ω) and k(ω) are rather smooth despite the oscillating behavior of the underlying data shown in Fig. 4c and 4d. This smoothness indicates that we correctly accounted for the effect of multiple internal reflections. In addition, n(ω) increases with frequency (normal dispersion), while k(ω) is almost zero. Because these behaviors are consistent with the literature 40 , it is proven that our evaluation method is valid.
As shown in the inset of Fig. 6c and in Fig. 6d, even though the amplitudes of the oscillating component in the n(ω) and k(ω) spectra are minimized, a small residual oscillating component with an amplitude on the order of ~10 -5 is still observed. The origin of this small component is considered to be the uncertainty of T(ω). In contrast to k(ω), the uncertainty of T(ω) has only a weak influence on the uncertainty of n(ω), since n(ω) does not directly depend on T(ω) as shown in equation (10).
The red dots in Fig. 6e show the values of the standard deviation of n(ω) at 193.414 THz, σn, obtained for different Nave values using the corresponding standard deviation of M, σM. σn gradually decreases in the range of Nave=125-1100 and it dramatically decreases between Nave=1100 and 1375 (more than 100 times smaller) because the value of M is uniquely determined for Nave≥1375. The blue dots in Fig. 6e correspond to the value of σn calculated using M=831 and σM=0. In this case, σn gradually decreases in the entire range of Nave and seems to approach a constant value for Nave≥1375. precisely.

Discussion
Firstly, we discuss the impact of the unique determination of M on the precision of n(ω). Because the leading term of n in equation (10) is 2πM, it is considered that σn is mainly determined by σM as long as M has not been determined uniquely; σn~(2πσMc)/dω. For the presently used experimental conditions, σn is evaluated to be 3×10 -3 σM, which governs the degree of uncertainty of the red dots in Fig. 6e for Nave≤1100. Once M has been uniquely determined, σn dramatically decreases by more than two orders of magnitude. This drastic decrease in σn is the main impact of the unique determination of M.
Secondly, we discuss the Nave dependence of σn using the uniquely determined M (Fig. 6e; blue dots). When σM=0, σn is proportional to the uncertainty of ∆ϕexp(ω), σ∆ϕ(ω) (See Supplementary   Information). In our experiment, σ∆ϕ(ω) is proportional to ave −0.5 for Nave≤1100 and tends to be constant for Nave>1100 (σ∆ϕ(ω)~10 mrad, data not shown). This lower limit of σ∆ϕ(ω) determines the lower limit of σn, which is 1.2×10 -5 . Furthermore, the standard deviation of d, σd, does also not decrease further when Nave>1100 (see Fig. 6b). This indicates that the lower limits of σn and σd in our experiment are not determined by the stability of our DCS system, but other experimental uncertainties due to changes in environmental conditions (such as temperature, pressure, and humidity), which change nair and n(ω). For instance, a temperature variation of 0.3 K causes a change in the refractive index of silicon on the order of 5.6×10 -5 , which may limit σn. A discussion of the contributions to the experimental uncertainty is provided in the Supplementary Information.
Thirdly, we discuss why we were able to precisely determine d even at Nave≤1100 where σM ≠ 0. As mentioned above, the local minima of E(M,d) appear at an interval of λ0/2nair when plotted as a function of d. As the slope of E(M,d) is quite steep around each local minimum (see Fig. 5b), σd becomes much smaller than λ0/2nair (~0.7 µm) even when we choose one of these local minima instead of the global minimum. In addition, because λ0/2nair is larger than the uncertainty of d for Nave=125 in our experiment, we can choose the correct minimum, i.e., the global minimum, for all values of Nave as shown in Fig. 6b.
Finally, we discuss the experimental requirements for a unique determination of M by our method. This means that the period of the artificial oscillating component in the derived N(ω) spectrum is determined by the optical-path difference between two adjacent optical pulses where one of these pulses is the result of multiple internal reflection. This optical-path difference is equal to 2n(ω)d.
Because σd is small even at smaller Nave values as mentioned above, the uncertainty of 2n(ω)d is governed by σn when σM ≠ 0. When the value of M changes by +1, the value of n(ω) in our experiment changes on the order of 0.003, which causes a change in the period of the artificial oscillating component in n(ω) by 0.09%. As a result, the period of the oscillating behavior of Δϕ(ω) with respect to frequency that is inversely estimated from the calculated n(ω), is not consistent with the measured oscillation period of Δϕ(ω). Because of this difference, the difference between the phases of the inversely calculated Δϕ(ω) and the measured Δϕ(ω) increases by ~0.6 mrad per oscillation. Because the number of oscillations in the measured Δϕ(ω) spectrum is ~60 within the measured frequency region (see Fig. 4d), the maximum phase difference becomes ~36 mrad. Thus, when σΔϕ(ω) is smaller than ~36 mrad, it is possible to uniquely determine M. By accumulating the data at each frequency, we verified that σΔϕ(ω)≈10 mrad. This suggests that an ultra-precise determination of n(ω) requires a precise determination of Δϕexp(ω), where σΔϕ(ω) is small enough to distinguish the small phase change associated with a change in M by ±1.
In summary, we have demonstrated a method that is able to simultaneously determine d and n(ω) of an undoped silicon wafer with σd=1.5 nm and σn=1.2×10 -5 , respectively. We determined the absolute phase spectrum by DCS and a simple criterion with respect to data smoothness. DCS can be used to measure both amplitude and phase with high frequency resolution and accuracy. The dramatic improvement in the precision of the absolute sample-induced phase change ΔΦabs(ω) that is achieved by our method, allows us to measure d and n(ω) simultaneously with degrees of precision that were previously unachievable.

Dual comb spectroscopy measurement
In this work we use DCS as described in ref 41 . We used two OFCs based on erbium(Er)-dopedfiber-based mode-locked lasers with slightly different repetition rates (fr and fr-Δfr). The two OFCs are referred to as the signal (S) comb and the local (L) comb. In our experiment, fr was ~48 MHz and Δfr was set to about 69 Hz. From each laser, one of the two output beams was amplified by an Erdoped fiber amplifier (EDFA) and its spectral bandwidth was broadened by a highly nonlinear fiber (HNLF). The obtained two output beams with the broadened spectrum were used to measure the interferograms. The other two output beams were used to detect the carrier-envelope offset frequencies of the S-and the L-combs, fS,ceo and fL,ceo, by f-2f interferometers. We also detected the beat notes between a 1.54-μm continuous-wave (CW) laser and one of the comb modes for both the S-comb and the L-comb (fS,beat and fL,beat). fS,ceo and fL,ceo were phase-locked to the radio frequency (RF) reference signals generated by function generators. Each carrier-envelope offset frequency was phase-locked via a feedback to the current of the corresponding pump LD. fS,beat and fL,beat were also phase-locked to the RF reference signals via feedback loops (to the electro-optic phase modulator, to the piezo-electronic transducer, and to the Peltier element in each oscillator) with different time constants. We controlled the values of fS,ceo-fL,ceo and fS,beat-fL,beat to obtain multiples of Δfr in order to satisfy the condition for coherent averaging 42 . As the CW laser is stabilized to an ultra-stable cavity, the relative line width between the S-comb and the L-comb is below 1 Hz, which is the resolution limit of the spectrum analyzer.
For the data acquisition process we first moved the sample to the measurement position. Then, we measured about 60 interferograms (recording time: 60/Δfr ≈ 0.87 s). The averaged interferogram is referred to as U1,w/(t) and contains the sample-related information of the first dataset. As 0.87 s is shorter than the coherence time of our DCS system, it is possible to average these 60 interferograms in the time domain. Then we removed the sample from the measurement position by moving the translation stage and again measured about 60 interferograms, and the averaged interferogram is referred to as U1,w/o(t). We repeated the whole procedure for about 10 hours to obtain 5500 pairs of

Uw/(t) and Uw/o(t).
To identify each single dataset, we use the notation Ui,w/(t) and Ui,w/o(t) with i = 1 to 5500.

Calculation of the transmittance and phase difference spectra
To use the sample-path and ref.
-path information (e.g. the peaks at 11.2 ns and 10.3 ns in Fig. 4a) separately, we extracted 7500 data points around the maximum intensity of each interference-signal peak and replaced the other data points with zeros. When the sample is located at the measurement Here, F means Fourier transform. As shown in equations (4) and (5), we use the ref.
-path information to cancel fluctuations in our optical setup, and we also use it as the reference phase.
tas(tsa) is the amplitude transmittance coefficient for transmittance from air to the sample (from the sample to air). rsa is the amplitude reflectance coefficient at the sample/air interface for light inside the sample. According to the Fresnel's law, these coefficients can be written as: as ( ) = 2 air ( )+ air (7) sa ( ) = 2 ( ) ( ) + air (8) sa ( ) = ( ) − air ( ) + air (9) Using equation (6) Here, we added the term for multiple internal reflections, 1 − sa 2 exp � 2 �, to the expressions of n(ω) and k(ω) derived by Tripathi et al. 44 . In our calculations, we simultaneously solve equations (10) and (11) by using the fsolve function of the software MATLAB.

Characterization of the levels of smoothness of n(ω) and k(ω)
For a correct characterization of the smoothness, it is important that the evaluation function only considers the artificially induced oscillating component in n(ω) and k(ω), that is, the oscillating behavior as a function of frequency that is caused by an incorrect analysis of the effect of multiple internal reflections in the T(ω) and Δϕ(ω) spectra. In undoped silicon, n(ω) linearly increases with ω in the measure frequency range. If we use incorrect values for M and d, n(ω) increases with ω but also contains a superimposed oscillating component, i.e. n(ω) can be expressed as n(ω)=α+βω+f(ω), where α and β are constants and f(ω) is the oscillating component. Therefore, by using the secondorder derivative, we can evaluate f(ω). Here, we consider that the second-order derivative is sufficient to evaluate the levels of smoothness of n(ω) and k(ω), because the second-order derivatives of n(ω) and k(ω) oscillate around zero. This means that all components except the oscillating component have been successfully removed from the data.