Propagation of Linear Uncertainties through
Multiline Thru-Reflect-Line Calibration

The eigenvalue problem is formulated by multiplying an $N\times N$ weighting matrix $\boldsymbol{W}$ on the right-hand side of (4a). Then, we multiply the new equation on the left-hand side of (4b). This finally results in

In general, the problem presented in (6) is a similarity equation between the matrices $\boldsymbol{F}$ and $\boldsymbol{H}$, where $\boldsymbol{X}$ is the transformation matrix. However, choosing an appropriate $\boldsymbol{W}$, the matrix $\boldsymbol{H}$ can turn into a diagonal matrix. As a result, the similarity problem in (6) transforms into an eigenvalue problem. From [12], the optimal $\boldsymbol{W}$ that maximizes the eigenvalue separation is found to be:

where $(\cdot)^{H}$ is the Hermitian transpose (conjugate transpose). The vectors $\boldsymbol{y}$ and $\boldsymbol{z}$ are given by

Therefore, under the condition that $\boldsymbol{W}$ was chosen according to (7), the eigenvalue problem is now given by

with

The calibration coefficients can be solved as the eigenvectors associated with the eigenvalues $\pm\lambda$. However, eigenvector solutions are only known up to a complex scalar factor. Therefore, the calibration coefficients obtained from the eigenvectors must be normalized to obtain a unique solution. Additionally, as $\boldsymbol{X}$ is constructed as a Kronecker product, it can be shown that the remaining inner columns of $\boldsymbol{X}$ can be obtained directly from the normalized eigenvectors of $\pm\lambda$ [12]. The denormalization of the eigenvectors is accomplished with the measurements of the reflect and the thru standards, while simultaneously solving for the 7th error term $k$. The equations for the denormalization can be found in [12]. Lastly, the propagation constant can be determined as a final step from the normalized calibration coefficients and the measurements of the line standards. In this step, the exponential terms of the line standards are extracted and combined in the least squares sense. This estimate of the propagation constant can be used to shift the reference plan if desired or used in other measurements such as characteristic impedance estimation [16] or permittivity measurement [17].

The weighting matrix $\boldsymbol{W}$ is an essential part of achieving a reliable calibration, as it optimally weights the measurement pairs according to their importance. Measurement pairs that are near singular are given less weight. In general, we can determine $\boldsymbol{W}$ by knowing the propagation constant and the length of the line standards, as given in (7). In [12], we originally proposed to compute the propagation constant from the measurements through non-linear optimization. However, numerical optimization procedures can pose issues in linear uncertainty propagation, as such methods are solved iteratively, and a proper convergence is required.

In this article, we propose a new technique for deriving $\boldsymbol{W}$ without having to compute the propagation constant. We first define a new equation by multiplying (4a) with the right-hand side of (4b). In this way, we get the following equation:

where

with $\boldsymbol{y}$ and $\boldsymbol{z}$ being the same as defined in (8).

From (11), it is clear that the equation does not depend on error-boxes. Additionally, from (12), we can see that its structure is very similar to the equation for $\boldsymbol{W}^{H}$ given in (7). The only difference between (7) and (12) is that the subtraction changes to a summation sign. This structural similarity implies that both equations share the same vector basis. To derive this vector basis from (11) and construct $\boldsymbol{W}$ accordingly, we use the Eckart-Young-Mirsky theorem [18, 19] by performing a singular value decomposition (SVD) and recovering the first two dominant singular values and their associated singular vectors. This process can be written as

where $s_{i}$, $\boldsymbol{u}_{i}$ and $\boldsymbol{v}_{i}$ are the singular values and associated left and right singular vectors of the measurements in (11). Furthermore, as (11) has a symmetric structure, we can use the Takagi decomposition to split the matrix into a symmetric product as

where $\boldsymbol{U}_{t}$ is a unitary matrix that holds the Takagi singular vectors and $\boldsymbol{S}_{t}$ is a diagonal matrix with positive real-valued singular values (equal to the singular values in SVD). As a result, $\boldsymbol{G}=\boldsymbol{U}_{t}\sqrt{\boldsymbol{S}_{t}}$. Takagi decomposition is analogous to the eigen-decomposition of real-valued symmetric matrices but for complex-valued symmetric matrices (without conjugation). In fact, the Takagi decomposition can be computed through SVD or eigen-decomposition [20, 21]. Therefore, we can simultaneously apply the Eckart-Young-Mirsky theorem and perform Takagi decomposition in a single step, either with SVD or eigen-decomposition.

Lastly, the matrix $\boldsymbol{G}$ is the common vector basis between equations (7) and (12). As a result, $\boldsymbol{W}$ is determined by

where the sign ambiguity results from not knowing the order of the sum from the matrix decomposition. This sign ambiguity can be resolved by choosing the answer which has the smallest Euclidean distance to a known estimate, for example, an estimate extracted from material properties. The complex number $j$ in (15) is to compensate for the square-root matrix of the permutation matrix in (12), which is implicit in the matrix $\boldsymbol{G}$.

Measurement uncertainties originate from the instrument itself. Such uncertainties could come in the form of additive or multiplicative noise. The general approach to model measurement uncertainties of S-parameters is through sample statistics. The covariance matrix of measured S-parameters can be determined from a sample covariance. If we assume a multivariate Gaussian distribution, then the unbiased sample covariance is given by

where $n$ is the total number of samples, $\boldsymbol{h}_{i}$ is the vectorized real-valued representation of the measurements of the $i$th S-parameters, and $\boldsymbol{\mu_{h}}$ is the corresponding mean value. $\boldsymbol{h}_{i}$ is defined similarly to (16),

In many modern VNAs, the wave parameters of a two-port device can be directly measured for a given port source. Using the wave parameters, we can compute the S-parameters by combining the measurements of both source directions [25], as stated below,

where $a_{ij}^{\prime}$ and $b_{ij}^{\prime}$ represent the measured incident and reflected waves at the receiver port $i$ when sourced by port $j$, respectively. The formulation in equation (21) accounts for switch terms [14], eliminating the possibility of correlation between measurements through switch term correction equations.

In the multiline TRL calibration, as we deal with T-parameters, the sample covariance can be propagated after converting from S- to T-parameters. Knowing the covariance of the measurements allows for the direct propagation of uncertainties at every step of the calibration through the chain rule of Jacobian matrices, ensuring accurate and consistent uncertainty calculations throughout the calibration process.

Forward model uncertainties are parameters that are considered during the calculation of calibration coefficients. An example of such a parameter is the length of line standards used to determine the propagation constant. Even though the propagation constant itself is not directly used in the calibration process, it can affect the reference plane by specifying its location. Additionally, when the calibration coefficients are denormalized using reflection measurements, it is possible to include the uncertainty of the unknown symmetric reflect standard by treating it as an independent variable in the denormalization equation. This can be done by substituting its estimated value into the corresponding Jacobian matrix.

The flow diagram of the forward uncertainty propagation in multiline TRL calibration is shown in Fig. 2. The Jacobian matrix is updated at each step, and the covariance of the calibration coefficients is determined using (18b). If the reference plane is shifted, the Jacobian matrix of the calibration coefficients shall be updated to account for the dependency on the shifted reference plane.

Inverse model uncertainty quantification is the process of identifying the uncertainty contributions arising from the model’s limitations used to estimate calibration coefficients. Various sources of uncertainty cannot be forward propagated in the multiline TRL algorithm, such as line mismatch, connecting/probing repeatability, and cable movement. For instance, the eigenvalue problem in multiline TRL assumes that all line standards have the same cross-section (i.e., the same characteristic impedance and propagation constant) and that error-boxes remain unchanged for each measurement. In reality, these assumptions are not always accurate due to randomness. However, on average, they are considered valid if the underlying phenomena are unbiased (i.e., their effects average to zero).

Various techniques can be used to address inverse model uncertainty quantification, such as numerical methods based on Bayesian modeling [26, 27]. In this work, a more straightforward approach is taken by reformulating the inverse problem as a forward uncertainty propagation by updating the covariance matrix of the measurements. It is important to note that the measurement noise and the model uncertainties are two independent statistical processes. Measurement noise originates from the VNA, while model uncertainties originate from the calibration standards. By assuming that each process follows a zero-mean multivariate Gaussian distribution, the overall distribution of the measurement remains multivariate Gaussian with the same mean value but with a total covariance comprising the sum of the measurement noise and model uncertainties as represented in the equation:

where $\boldsymbol{\Sigma}_{\mathrm{N}}$ is the covariance due to noise, $\boldsymbol{\Sigma}_{\mathrm{F}}$ represents the model uncertainties that can be directly accounted for in the forward uncertainty propagation, and $\boldsymbol{\Sigma}_{\mathrm{I}}$ represents every source of uncertainty that is addressed with inverse uncertainty propagation. In this work, we focus on the line mismatch, as it can be easily modeled in relation to the measurements. Still, the presented technique can also be applied to other uncertainty types, such as cable movement, connection repeatability, or probing uncertainties, if they can be adequately modeled in relation to the measurements.

To quantify the line mismatch uncertainty, we can use the generalized mismatched line model in the error-box model given by

where $\boldsymbol{L}_{i}^{\prime}$ is the $i$th mismatched line model given by [28],

with $\Gamma_{i}=\delta_{Z_{i}}/(2\mu_{Z}+\delta_{Z_{i}})$ and $\gamma_{i}=\mu_{\gamma}+\delta_{\gamma_{i}}$, where $\delta_{Z_{i}}$ and $\delta_{\gamma_{i}}$ are random processes that describe uncertainties in the impedance and propagation constant of the $i$th line standard. Under the assumption that $\delta_{\gamma_{i}}$ and $\delta_{Z_{i}}$ are jointly distributed after zero-mean Gaussian distribution, then the joint distribution of $\Gamma_{i},\gamma_{i}$ is also Gaussian, $\Gamma_{i},\gamma_{i}\sim\mathcal{N}\left([\mu_{\Gamma},\mu_{\gamma}],\boldsymbol{\Sigma}_{\Gamma_{i},\gamma_{i}}\right)$. Following that $\delta_{Z_{i}}$ is assumed to be zero-mean Gaussian, then it follows that $\mu_{\Gamma}=0$. Additionally, for the value of $\mu_{\gamma}$ we use the estimated value computed by the multiline TRL calibration. The joint covariance $\boldsymbol{\Sigma}_{\Gamma_{i},\gamma_{i}}$ is determined through propagating uncertainties in the cross-section of the transmission line (e.g., martial properties and dimensional quantities) by using an appropriate model (e.g., analytical or electromagnetic simulation).

Lastly, to determine $\boldsymbol{\Sigma}_{\mathrm{I}}$, we only need to propagate $\boldsymbol{\Sigma}_{\Gamma_{i},\gamma_{i}}$ through the error-box model in (23), using the estimated values of the calibration coefficients and the propagation constant from the multiline TRL as mean values. Thus, $\boldsymbol{\Sigma}_{\mathrm{I}}$ for the $i$th line standard is given by

With knowledge of $\boldsymbol{\Sigma}_{\mathrm{I}}$, we update the covariance matrix of the measurements and propagate it in the forward uncertainty propagation, as discussed in the previous subsection. The flow diagram in Fig. 3 summarizes the steps for inverse uncertainty propagation of the line mismatch.

It is important to note that the definition of uncertainties in the transmission line model can be rephrased in alternative terms. For instance, instead of using the reflection coefficient $\Gamma$, the analysis could be based on the characteristic impedance $Z$, and instead of using the propagation constant $\gamma$, the relative effective permittivity $\epsilon_{\mathrm{r,eff}}$ could be used. It is also important to remember that the Jacobian matrix needs to be adjusted accordingly.

The measurement setup includes a VectorStar VNA from Anritsu, with mm-wave heads up to 150 GHz, used together with a SUMMIT200 probe station from FormFactor. The CPW calibration standards were implemented on a commercial impedance standard substrate (ISS) from FormFactor. The measurements were conducted using $150\,\mu\mathrm{m}$-pitch ground-signal-ground (GSG) ACP probes with a 0.8 mm connector interface from FormFactor. Fig. 5 depicts the measurement setup. The calibration standards comprise six CPW lines with edge-to-edge length $\{200,450,900,1800,3500,5250\}\,\mu\mathrm{m}$ and an open standard, which is measured by setting the probes to float.

The measurements were conducted in the frequency range 1 to 150 GHz with a step of 1 GHz, resulting in 150 points. A low number of frequency points was necessary to reduce the measurement time, as 500 frequency sweeps of the wave parameters were collected for each standard at an IF-bandwidth of 1 kHz and power level of -10 dBm. It took less than 10 minutes to collect all sweeps for each standard, resulting in a measurement time of under 70 minutes for all standards. The S-parameters were calculated and the sample covariance matrix for each standard was estimated. Fig. 6 shows scatter plots and the corresponding marginal histograms of the uncalibrated measurements of $S_{11}$ and $S_{21}$ from the $1800\,\mu\mathrm{m}$ line standard at 120 GHz.

To model the CPW structures used in the measurements, we used the dimensional and material parameters along with their uncertainties. Since publicly available data on these parameters are limited, we conducted optical measurements for dimensional parameters, and obtained material properties from publicly available sources. Gold was used for the conductor layer, and lossless Alumina with a relative permittivity of 9.9 was used for the substrate. To demonstrate the propagation of uncertainty in multiline TRL calibration, we assumed $\pm 10\%$ standard uncertainty for the smallest planar dimensional feature, $\pm 10\%$ for both the thickness and the conductivity of the Gold layer, and $\pm 0.2$ for the relative permittivity of the Alumina substrate. These uncertainty assumptions are for demonstration purposes only and are not intended to be considered accurate. Table I summarizes the parameters and their uncertainties.

The CPW model we are basing our analysis based on is the analytical model presented in [29, 30, 31]. This model includes corrections for radiation loss from references [29, 30]. A comparison between the CPW model based on the data in Table I and the extracted relative effective permittivity and loss per unit length from the multiline TRL measurements is shown in Fig. 7. The measurements and the model are in reasonable agreement. Therefore, the MC analysis given in the following subsection should represent a realistic scenario of the CPW on-wafer measurements. It is worth noting that the offset seen in the relative effective permittivity can be attributed to the fact that the calibration substrate has a thickness of $254\,\mu\mathrm{m}$ placed on an absorbent chuck, whereas the CPW model assumes an infinite extent of the substrate.

The uncertainties considered in our analysis are measurement noise, line mismatch, length offset, and reflect asymmetry. Measurement noise was modeled based on the sample covariance from the measurements. Line mismatch uncertainties were established by perturbing the cross-sectional parameters based on the values in Table I. The length offset was assumed to have $\pm 40\,\mu\mathrm{m}$ standard uncertainty to account for the contact overtravel of the ACP probes. The reflect asymmetry was introduced as an offset of $\pm 40\,\mu\mathrm{m}$ standard uncertainty to emulate the effects of changes in the open capacitance.

For the MC analysis, we generated S-parameter data of the calibration standards using the CPW model and embedded them into the error-boxes obtained from the multiline TRL measurements. We perturbed the relevant parameters in each realization and collected the corresponding output. The uncertainty of the output parameters was estimated as sample covariance. We conducted the MC analysis using 5000 samples to ensure reliable numerical convergence. While it was possible to use more samples, we noted that using as few as 100 samples yielded results comparable to those achieved with 5000 samples. It is worth noting that the MC analysis with 5000 samples took more than an hour and a half to complete.

For linear uncertainty propagation, we estimated the necessary covariance matrices. The covariance for line mismatch was estimated using numerical calculation of the Jacobian matrix of the CPW model, as the analytical model proposed by [29, 30, 31] contains the complete elliptic integral of the first, second, and third kind, which hinders the ability to compute analytical derivatives. The covariance of the asymmetry of the open standard was computed analytically by its Jacobian matrix. We implemented the multiline TRL algorithm in Python and used the package scikit-rf [32] for processing S-parameter data and the package Metas.UncLib [23] to automatically compute Jacobian matrices and propagate uncertainties linearly.

To analyze the impact of uncertainty in calibration, we considered a symmetric, lossless DUT with equal transmission and reflection to be calibrated. In Fig. 8, we compare MC analysis and linear uncertainty propagation for relative effective permittivity and loss per unit length, as well as $S_{11}$ and $S_{21}$ of calibrated DUT. From the plots we can see an excellent agreement between the two methods, validating our approach. Since we are running a finite number of MC trials and the linear propagation remains an approximation, there is an error tolerance between the derived uncertainties from the two methods. On average, over all frequency points, the relative error of the estimated uncertainties from the linear propagation method compared to the MC method is about $0.6\%$ for the relative effective permittivity, $5.33\%$ for the loss per unit length, $4.61\%$ for the magnitude of $S_{11}$, and $4.99\%$ for the magnitude of $S_{21}$. It is worth noting that for the MC method, we perturbed each parameter directly in the simulation, whereas in the linear uncertainty propagation, we only provided the covariance matrices. This means that if the covariance matrices are available (even if it is a rough estimate), the uncertainty evaluation can be performed quickly, simultaneously with the calibration.

Linear uncertainty propagation allows for a quick evaluation of the uncertainty budget caused by different sources of uncertainty. As shown in Fig. 9, the impact of each type of uncertainty can be analyzed. For example, length uncertainty and line mismatch have significant uncertainties in the extracted relative effective permittivity. In addition, noise in measurements affects all quantities significantly, except for the relative effective permittivity. As observed in the uncertainty budget, the noise contribution from the VNA is higher than expected at higher frequencies, most significantly after 100 GHz. After contacting the VNA manufacturer, we discovered that the issue originated from additional statistical error introduced by the VNA firmware. The manufacturer was able to fix this issue with an updated firmware version. This experience emphasizes the usefulness of analyzing the uncertainty budget, as it allows for proper debugging.

Furthermore, by breaking down the uncertainty contribution, we can see how each standard contributes to the overall uncertainty, as shown in Fig. 10, which can be used to optimize calibration standards and identify measurement bottlenecks.