Analytical Uncertainty Evaluation of Material Parameter Measurements at THz Frequencies

Material parameter extraction algorithms are studied and simplified both for transmission-reflection and transmission-only methods. The simplified relations, which are closed-form in some cases, are analyzed to establish the uncertainty sensitivity coefficients and therefore, to clarify the main uncertainty contributions and reduce the systematic and random errors. Simple closed-form expressions presented in this paper show the sensitivity of the extracted permittivity to each input parameter such as S21 (phase and amplitude), frequency, and the material thickness. Results are presented for several material slabs for three waveguide frequency ranges 75–110 GHz, 140–220 GHz, and 500–750 GHz using VNA-based free-space technique in the THz domain. Comparison of results (and the associated uncertainties) between different algorithms can help to choose the optimal one suitable for lossy or low-loss materials, and thin or thicker slabs. This can explain why the same set of S-parameters data usually gives different final results (permittivity and permeability) with different algorithms and verify the reliability of the calibration and extraction process.


Introduction
Measurement of complex permittivity and permeability is very important in THz and millimeter-wave technology (substrates and dielectric components), biomedical and space applications, material engineering, and others. There exist different measurement techniques for different frequency ranges and materials which can cover lossy to low-loss slabs with various sizes and thicknesses [1]. Closed and open transmission cells and waveguides [2] or resonant cavities [3] are usually used in the microwave region. Waveguide or cavity methods, however, seem to be less applicable in the mm-Wave/THz domain because of mechanical problems and dimensional restrictions. Free-space non-contact methods are therefore used largely at these frequencies [4,5]. Laser-based optical techniques have been, for a long time, the only solution for THz spectroscopy [6,7], but recent progress in high-frequency semiconductor technology made available electronic-based components in the THz domain. Works [8,9] discuss fundamental problems in reliable data analysis with known measurement uncertainty. A vector network analyzer (VNA) can be used as the core measurement unit with quasi-optical associated devices to set the free-space non-contact measurements up. The setup has to be adequately calibrated to derive S-parameters on the surface of the material-under-test (MUT) [4]. Monte Carlo simulation is also used to evaluate the measurement uncertainty of the VNA [10]. Based on the measured S-parameters and the MUT thickness and positioning, the permittivity and permeability can be extracted.
Material parameter extraction methods are based on reflection [11] and transmission parameters (or only one of them [1,12,13]). The algorithms usually suffer from instability or complicated iterative process and use genetic algorithms [14] or other numerical/optimization techniques [1,12,13,15]. Also some methods require certain a priori knowledge of the material [4,12,16]. It is often observed that with the same set of S-parameters data, the extracted permittivity and permeability are not the same with different algorithms and therefore, the reliability of the full process of the material characterization gets questioned. We believe that the best solution here is the "sensitivity analysis" of the extraction methods. The abovementioned classic algorithms and many others are usually complicated and the sensitivity/error analysis may need a more complicated numerical process. By establishing reliable simplified (mostly closed-form expressions) extraction models, we will be able to show the sensitivity in a comprehensive and clear manner. The "simplification" and "sensitivity analysis" are the key points and main novelties of this work. This paper presents simplified relations for both transmission-reflection (Section 2) and transmission-only (Section 3) methods with analytical closed-form expressions of uncertainty sensitivity coefficients. Measurement results are based on non-contact free-space techniques by using VNA and the associated frequency extensions for the bands 75-110 GHz, 140-220 GHz, and 500-750 GHz. Non-magnetic materials with relative permeability μ r = 1 are assumed in this paper. The quasioptical device here is a commercially available mode-convertor (Swissto12@ MCK [17,18] set of two corrugated horn antennas). The system is supposed to convert waveguide propagation modes to TEM free-space mode on its aperture where the MUT is installed. We have been working on the setup optimization and in-depth study of its electromagnetic characteristics (propagation modes, matching, losses, aperture/gap scattering, etc.) and the system calibration. The study of the setup and its optimization is mainly based on theoretical and experimental works of [19,20]. This paper is focused on extraction methods, only, whereas both the VNA calibration and the associated frequency extender optimization are performed to minimize the uncertainties of the input parameters (S 21 , amplitude, and phase).
The results are presented at higher frequencies and in the mm-Wave/THz domain, but we believe the methods and analysis process presented here can be used for other frequency ranges and therefore helpful for the microwave and THz community, in general.

Transmission/Reflection Method
Reflection (S 11 , S 22 ) and transmission (S 21 , S 12 ) data can be processed to obtain real and imaginary parts of both permittivity ε = ε 0 ε r and permeability μ = μ 0 μ r of the MUT (the relative permittivity is expressed as ε r = ε r − jε r ). Various iterative and non-iterative algorithms have been presented for this goal [1,12,21]. The classic Nicolson-Ross-Weir (NRW) method [1] deals with Γ and T (T is the transmission coefficient and Γ is the reflection coefficient) to extract μ and ε, simultaneously. The method gives unstable results for low-loss and high-refractive materials when S 11 ∼ 0 and |S 21 | ∼ 1 at specific frequencies where and d is the material slab thickness, c is the free-space speed of light, and ω is the angular frequency. A comprehensive method has been presented [22] to correct the instability. Using the classic (1), (2), and (3) to simplify the final relations relevant to ε and μ as a closed-form formula, we obtain The expansion of Eq. 6 leads to The material extraction method using Eq. 8 is independent of the material thickness (i.e., the material thickness d is not a variable here and the extraction process is not affected by improper input data d) and can give reliable and stable results if S 11 is not very small. For some materials, e.g., body tissue or substrates for wearable antennas, the thickness d is difficult to measure precisely. Yet, the instability problems appear for low-loss materials when S 11 ∼ 0 and |S 21 | ∼ 1 at specific frequencies which is equivalent to the classic NRW method [1,21]. Similar relations have been already used to determine the impedance of transmission lines with dielectric-filled sections [24]. Closed-form Eq. 8 provides a more comprehensive analysis of the uncertainties with a clear contribution of each measured S-parameter. By applying the partial derivatives of Eq. 8, sensitivity coefficients of ε r can be derived as a function of S 11 and S 21 . Simple and closed-form structure of Eq. 8 makes its derivatives closed-form, as well (10) Figure 1 shows the sensitivity coefficients of the permittivity for a 3-mm Plexiglass slab. As shown, for both S 11 and S 21 , there are nearly simultaneous minimum sensitivity points which are in agreement with the results of "instability corrected" curves in [22]. Given nominal |ΔS 11 | and |ΔS 21 | approximately 0.015 (the overall measurement uncertainty of S 11 and S 21 due to repeatability, VNA stability, noise, drift, etc.), the associated uncertainty of |ε r | is estimated around 0.1 (∼ 4%, relative) when ε r = 2.6. Therefore, the simple closed-form (yet unstable) relation (8) can give reliable results with calculable uncertainties around the frequencies 140 GHz, 170 GHz, and 200 GHz (Fig. 1). Figure 2 shows the sensitivity coefficients of the permittivity for a high-refractive low-loss 1.2-mm sample YZrO 2 . As shown, for both the S 11 and S 21 , there are nearly simultaneous minimum sensitivity points. In this case, the |ΔS 11 | and |ΔS 21 | are approximately 0.015, and the associated uncertainty of |ε r | is estimated around 7% when ε r = 33.7. The best uncertainties are observed at the frequency points around 157 GHz, 179 GHz, 201 GHz, and 220 GHz (Fig. 2). The minima and maxima of the uncertainty can be explained by analysis of the denominator of Eqs. 9 and 10. The uncertainty is minimal when the |S 21 | is minimal and the |S 11 | is maximal and vice versa.  (10) Comparing with the results of Fig. 1, we can see the relative uncertainty is higher and this shows the difficulties of extracting parameters for high-refractive thin materials with the transmission-reflection method. To demonstrate these challenges, we have calculated the permittivity (real and imaginary parts) of the same sample of Fig. 2 with two well-known extraction methods: stable non-iterative [21] and NISTiterative [12] (Fig. 3). METAS VNA-Tools II [25] has been used to show the extracted and ε r (bottom), with two well-known classic methods, NIST (black) and stable non-iterative (red). The methods show diverged results and partially "unphysical" for ε r material parameters with classic methods. VNA-Tools II is a metrology software to compute uncertainties of the measured S-parameters with various post-processing features including material parameter extraction with associated uncertainties.
As shown in Fig. 3, the real part of ε r varies between 33.5 and 34 and the results of imaginary part ε r are much more deviated or even meaningless for some frequencies.
In the next section dedicated to the transmission-only method, we show how we can find more reliable extracted parameters by simplifying the algorithm and analyzing the associated uncertainties.
Here at the end of Section 2, we will show an original application of the simple closed-form formula (8) and its partial derivatives (9), (10) [22]. Actually, reliable results of ε r (real part) even for a limited number of points in a given frequency range are of interest. Most of the "normal" homogeneous materials have nearly constant ε r in the mm-Wave/sub-THz domain. Besides, the calculation of ε r and its associated uncertainties is independent of the slab thickness. The thickness of the slab and its positioning can be one of the major uncertainty contributions (cause a systematic error) in reflection-and-transmission methods [26] and transmission-only methods [9], as well.
As mentioned, closed-form Eq. 8 is independent of the slab's thickness and can correct the systematic errors caused by this sensitive parameter. Figure 4 shows the Fig. 4 A 2-mm fused silica (F.Si) slab, ε r , and sensitivity coefficients of ε r (140-220 GHz). Extracted ε r with thicknesses 2.00 mm (top), 2.04 mm (middle), and 2.07 mm (bottom). Non-iterative stable method (black) and Eq. 8 (red) are used, optimal match (ε r = 3.8) with 2.07 mm thickness at the "best" uncertainty points ∼ 167 GHz and ∼ 205 GHz high sensitivity of the extracted ε r (real part) of a 2-mm (nominal thickness) fused silica slab.
The parameter extracted by the non-iterative stable method [21] varies about 5%, which is a considerable systematic error. This is regarding ε r = 4.05 (for nominal 2 mm thickness), to ε r = 3.9 (2.04 mm thickness, measured with high-precision mechanical devices) and finally ε r = 3.8 (for 2.07 mm thickness, which gives reasonable match to Eq. 8 at the best uncertainty points ∼ 167 GHz and ∼ 205 GHz).
It can be shown that the ε r uncertainty minimum points occur when the measured |S 11 | and |S 21 | reach their maximum and minimum (Fabry-Perot multiple reflection inside the slab), consequently. This fact may still make the extraction and uncertainty evaluation process easier.
Measurement of reflection parameters using a VNA is problematic at higher frequencies, because of the lack of suitable short standard and positioning challenges [4] with micrometer precision. The short standard must be perfectly conducting and ideally reflective (i.e., polished), which is not easy at sub-millimetre frequencies. The positioning (straightness and tilt, exact and repeatable position) is crucial especially for the measurement of the S 11 phase. For non-magnetic materials (μ r = 1 − j0), the transmission-only method (the S 21 parameter, only) can be more reliable and easier in terms of setup calibration and parameter extraction. The |ΔS 11 | is usually bigger than |ΔS 21 | and the total uncertainty of the method using S 11 is then bigger (note that the sensitivity coefficients must be multiplied by the |ΔS 11 | or |ΔS 21 | to get the final measurement uncertainty).

Transmission-Only Method
The transmission-only method utilizes the S 21 (or S 12 ) parameter on the material slab to derive T and consequently, to extract the permittivity ε r = ε r − jε r . Here T is a function of the material thickness d and its permittivity and permeability (5). For the materials which are not very lossy (ε r /ε r = tan δ < 0.1), T can be separated directly to phase and magnitude parts [4] This means the phase and magnitude of T give real and imaginary parts of permittivity, consequently ln |T | = − ωd 2c ε r ε r .
Here, we suggest a very easy and original way to evaluate φ (T ) from φ (S 21 ) which are related by the Fabry-Perot relation (2) with the unknown Γ .

Simple Expression for the Phase of T
First, we look in the phase of 1 − Γ 2 in Eq. 2, if the material is not very lossy (ε r /ε r = tan δ < 0.1). In this case, the phase is equal to zero with a very small error. Then, we search the points for which the phase of 1 − T 2 Γ 2 is zero. If φ T 2 = nπ (n ∈ Z), then φ 1 − T 2 Γ 2 = 0 and In order to better understand (14), suppose that the phase φ (T ) is equal to nπ/2, then calculate the phase φ (S 21 ) with the assumption of φ Γ 2 ≈ 0. This means at the points where the measured φ (S 21 ) is equal to nπ/2, the φ (T ) has also the same value and therefore, the complicated equation (2) does not need to be solved and used as the extraction parameter basis. Figure 5 shows the phase of S 21 and T (approximated) for the material from Fig. 2 (1.2-mm YZrO 2 sample) and their equality at "nπ/2" points. In this figure, we have also calculated the real part of permittivity from Eq. 12 at the points where φ (S 21 ) = nπ Here, ε r is intentionally calculated only for "nπ " points rather all "nπ/2" points, as we will show later that at "nπ " points the uncertainty (sensitivity coefficient) is minimum in the transmission-only method. It can be shown that around "nπ " points the |S 21 | reaches its maximum (Fig. 5). Figure 5 shows the reliability of the simple method (14), (15) with a constant ε r = 33.7 comparing with Fig. 3 in which well-known classic methods are used. Once ε r is calculated at "nπ " points with minimum uncertainty, Γ 2 is calculated (by definition) Note that the calculated Γ 2 in this step is a constant for the whole studied frequency range. Meanwhile, from Eq. 2, T will be determined (as a complex number for all the frequency points) with a quadratic relation where the condition |T | < 1 is imposed. The calculated T is the key for the extraction of permittivity from Eqs. 12 and 13.
To show the reliability of the simple method demonstrated in Eqs. 14 and 15, here we present an interesting and challenging case, i.e., results for a high-refractive lowloss 0.6-mm alumina slab in the 500-750-GHz band, for which the S-parameters are not time gated, see Fig. 6.
In general spectroscopy, materials with peaks in the extracted permittivity (spectral fingerprint) show sharp changes both in the amplitude and phase of S 21 and Fig. 7 A 1.2-mm YZrO 2 sample, 140-220 GHz: uncertainty sensitivity coefficient with transmission-only method (bottom) and calculation of ε r with best uncertainties at φ (S 21 ) = nπ = φ (T ) points (top), and the best fit (black solid line) for the entire frequency range 0.13 < ε r < 0.17 Fig. 8 A 3-mm plexiglass slab, 140-220 GHz: uncertainty sensitivity coefficient with transmission-only method (bottom, solid red curve is |S 21 |, dashed blue curve is the sensitivity) and calculation of ε r with best uncertainties at φ (S 21 ) = nπ = φ (T ) points (middle), the best fit (black solid line) for the entire frequency range 0.030 < ε r < 0.034. Comparing with transmission-reflection method (top), where the well-known classic methods give diverged results 0.025 < ε r < 0.039 are usually observed at higher THz frequencies, although some energetic materials exhibit peaks at frequencies below 500 GHz [23]. The material parameter extraction around these points may be directly based on S 21 phase and amplitude and then corrected for the Fabry-Perot effect and multiple reflections.
In the next part of this section, the sensitivity analysis of transmission-only equations will be presented to calculate real and imaginary parts of permittivity with estimations of associated uncertainties.

Uncertainty Analysis of the Transmission-Only Method
Partial derivative process on Eq. 2 gives sensitivity coefficients of T as a function of At all the points and especially at "nπ " points, ∂Γ 2 is much smaller than ∂S 21 (and negligible); thus, Eq. 18 can be shortened as (20) Now, by using closed-form equations for real and imaginary parts of permittivity (12), (13) and performing partial derivatives, the sensitivity coefficients are determined ωd 2c ε r ∂ε r = ∂|T | |T | ≤ |∂ S 21 |  Figure 7 shows ∂T /∂S 21 for the YZrO 2 sample with thickness 1.2 mm with the minimum sensitivities at the points where |S 21 | is maximum (i.e., "nπ " points of φ (S 21 ), the phase of S 21 ). These are the points where the imaginary part of ε r has the minimum uncertainty and the real part of ε r is directly calculated from φ (S 21 ).
By comparing Figs. 5 and 7 and Fig. 3, the advantages of the simplified transmission-only equations, presented in this paper, and comprehensive uncertainty analysis can be demonstrated. Another set of results is presented in Fig. 8 for a thicker 3-mm plexiglass slab to compare simplified transmission-only method with well-known classic algorithms (NIST-iterative and stable non-iterative) based on transmission-reflection data. The S-parameters data in the 140-220 GHz range are Fig. 11 Float-glass 1.5 mm, 500-750 GHz: ε r and ε r time gated (smoothening) and have been provided from Swissto12@ [26]. They have been processed in this paper to compare different extraction algorithms and evaluate/analyze the uncertainties.
This method needs to observe a couple of minima and maxima on the |S 21 | curves ("nπ/2" points of φ (S 21 )). This occurs if the material is not very thin (electrical length in the THz domain) and very low-refractive. Examples here can be an alumina slab (ε r ∼ 9.5) with 1 mm thickness or more, and plexiglass (ε r ∼ 2.7) thicker than 2 mm for the 75-110-GHz band. This can change to thin slabs, i.e., alumina with the thickness 0.2 mm and plexiglass with the thickness 0.4 mm (relevant to the electric-length) for the 500-750-GHz band.
Note that the S-parameters and their measurement uncertainty depend on the MUT conditions (positioning, surface quality), environmental conditions (VNA noise, drift, stability), VNA calibration and measurement repeatability, etc. Most of these influences can be demonstrated by |ΔS 11 | and |ΔS 21 |. In this paper, Fig. 12 Manufactured lossy material, 4 mm sample, 75-110 GHz: ε r and ε r we focused mostly on the material parameter extraction and the sensitivity coefficients.

More Results with the Simplified Transmission-Only Equations
We present more results for two other waveguide frequency ranges 75-110 GHz and 500-750 GHz with different types of material. Swissto12@ MCK (Fig. 9) is used and calibrated for transmission data S 21 . For these measurements, it is chosen not to use time-gating option of the VNA in order to have a more clear perspective of the overall uncertainty budget.
For the results here (Figs. 10, 11, 12, and 13), the estimated uncertainty for ε r is 3-4% and 10-15% for ε r . It depends on the material thickness d (as one of the most important contribution to the uncertainties ∂d/d), |T |, and other parameters in Eqs. 21 and 22.  Figure 13 shows an interesting behavior of (manufactured) high-loss materials in THz domain for which the real part of permittivity decreases with frequency.

Conclusion
Transmission-reflection and transmission-only methods have been studied and reviewed. Simplified closed-form equations are presented to correct the systematic errors and analyze the uncertainty sensitivity coefficients in material parameter extraction. Various materials (low-loss, high-refractive, lossy, thin, and thick) are measured in THz domain for three waveguide frequency ranges: 75-110 GHz, 140-220 GHz, and 500-750 GHz. The presented method helps to understand the uncertainty (and reduce it) and the reliability of the relevant algorithms.
It is interesting to mention that the "sensitivity" curves of transmission-reflection method and transmission-only method give minimum uncertainty at different points (often opposite, regarding the maximum and minimum points of |S 21 |). This may help to consider combined methods to reduce the uncertainties at wider frequency bands.