Linear Propagation of Uncertainty in Probe Position Compensated Multiline-TRL

1 Linear Propagation of Uncertainty in Probe Position Compensated Multiline-TRL Robin Schmidt dt dtdt dt , Marco Garelli, Member, IEEE, Andrea Ferrero, Fellow, IEEE, Michael Dieudonn´ e, Dominique Schreurs, Fellow, IEEE Abstract—The quantification of uncertainty sources in onwafer S-parameters measurements is crucial to clarify the quality of the extracted models used in Process Design Kits. In fact, a clear uncertainty budget expands our understanding of where first actions can be taken to improve such measurements. In this paper, we propose a simplified approach to uncertainty quantification without the requirement of any waveguide or coaxial connection, e.g., without removing the probes. Then, we develop the uncertainty propagation algorithm applicable to these uncertainty sources. Index Terms—Broadband communications, coplanar waveguides, instrumentation and measurement techniques, on-wafer calibration, planar transmission lines, I. INTRODUCTION ON-WAFER S-parameters measurements allow accurate characterization of devices and eliminate the inaccurate and often cumbersome process of de-embedding test fixtures. However, some issues may arise when using test probes, mainly due to contact repeatability, instrumentation drift, and signal leakage into parasitic modes [1]. Many of these undesirable effects can be minimized as demonstrated in [2], [3], [4], improving the accuracy of these measurements. A good example can be found in [5], where Arz was able to obtain traceable on-wafer measurements with a comprehensive uncertainty budget. More particularly, the uncertainties of the test fixture were evaluated with Euramet’s guidelines [6], which requires probe disconnection and calibration at the connector. In a similar context, we were able in [7] to significantly improve measurements suffering from probe misplacements on contact pads by quantifying the modification of error boxes with the probe position. However, the uncertainty quantification performed in [7] was limited to calibration residuals using expressions from [8], noise, and probe position uncertainty. In the current paper, our aim is to complete the uncertainty study from [7] and include test fixture drift and substrate proManuscript received March 29, 2025; revised May 23, 2025; accepted June 8, 2025. The work was in part supported by the 23IND10 OnMicro project. The project (23IND10 OnMicro) has received funding from the European Partnership on Metrology, co-financed from the European Union’s Horizon Europe Research and Innovation Programme and by the Participating States. (Corresponding authors: Robin Schmidt) Robin Schmidt, Andrea Ferrero, Marco Garelli, and Michael Dieudonn´ e are with Keysight Technologies Inc. (e-mail: [email protected]) Robin Schmidt and Dominique Schreurs are with KU Leuven. Color versions of one or more of the figures in this paper are available online at http://ieeexplore.ieee.org. Digital Object Identifier: 10.1109/TMTT.2025.3580938 cessing variations, before propagating them in our versions of the multiline-TRL algorithm derived from [9], [10] and in the extracted displacement functions of the probes. Additionally, we propose an approach to noise quantification less sensitive to drift and a simplified method for drift quantification, both possible when probes are still connected to the Vector Network Analyzer. Receiver linearity errors were however not considered in the frame of this research. In the first part of this paper, we clarify how noise, drift, and wafer inhomogeneity are estimated. The second part consists of Jacobian-based first-order propagation of uncertainty in both the calibration algorithm and the displacement function used for probe misplacement compensation. Finally, we compare uncertainties propagated in the different implementations of the multiline-TRL, including the most recent versions from [11], [12]. II. UNCERTAINTY QUANTIFICATION A. Uncertainty Due to Noise S-parameter measurements are based on the ratio of complex numbers coherently measured by I/Q receiver chains that receive a signal from directional couplers. These receiver chains are subject to both thermal noise and oscillator phase noise, which are responsible for trace noise and noise floor in S-parameters. The modeling of these first uncertainty sources can be done based on [13], where multiplicative and additive random variables are introduced in each receiver of the VNA, following the model (1), with ˆ xithe actual value of the reading, δnHi and δnLi the complex random variables for highand low-level noise. xi=ˆ xi(1 + δnHi) + δnLi (1) In [13], a simple experiment is proposed to estimate the contribution of highand low-level noise to the measurements, based on several receiver measurements made at two different signal levels. Following this method, for a system of at least two ports, no additional load connection is required, as a low signal level in the receivers can be achieved when the opposite port is active. This allows us to write the expression (2), where the measurements of a single standard with a high reflection coefficient and a transmission coefficient close to the noise floor are sufficient to recover the variance of the different receiver noises. In on-wafer measurements, it is in general possible to reach these requirements by lowering the chuck and separating open probes by several mm. © 2025. Personal use of this material is permitted. However, permission to reprint/republish this material for advertising or promotional purposes or for creating new collective works for resale or redistribution to servers or lists, or to use any copyrighted component of this work in other works must be obtained from the IEEE. Link to publisher version with DOI: 10.1109/TMTT.2025.3580938 2 xi,i xi,j=xi,i 1 xi,j 1∗δnHxi δnLxi with j6=i(2) Since these measurements may take some time, the receiver measurements xi,i of (2) can be subject to instrumentation drift, as it may be the case for older versions of frequency extenders, or for extremely low high-level noise (e.g., like for the new PNA DDS sources [14]). Therefore, the total distribution of xi,i will also be affected by drift, effectively increasing the variance of the distribution. In such case, local differences prove to be much more robust against drift when quantifying noise. This introduces equations (3) that use the distribution of local differences instead. diff(xi,i) diff(xi,j )=xi,i 1 xi,j 1∗√2δnHxi √2δnLxi with j6=i(3) In Fig. 1, the validity of the new equations was tested by extracting noise levels with repeated receiver measurements with a short connected to both ports, while the RF power level was slowly drifting. We clearly observe that drifts had much less impact on this new extraction methodology as the acquisition time increases and the high-level noise term reduces with the IF bandwidth. Fig. 1: Port 1 Receiver noise levels extracted for various IF Bandwidth (dotted: from [14], continuous: This Paper) Noise models found in the literature usually decompose each S-parameter with a multiplicative part and an additive part [6], [15]. Using the formulation (1), we can derive, to the first order, the expression (4) where we easily identify the composition of the multiplicative and additive parts in any receiver ratio, forming the raw S-parameters or the switch terms [16]. xi yj =ˆ xi ˆ yj (1+δnHxi −δnHyj −1 |ˆ yj|δnLyj)+ 1 |ˆ yj|δnLxi (4) In (4), we understand that the eventual correlations of the high-level noise in the receivers must also be taken into account if one wants to retrieve the actual noise levels in the ratio measurements. In fact, the phase noise of the VNA oscillators is mainly responsible for the high-level noise, and it is consequently clear that the high-level noise will be correlated. Meanwhile, we assume that the noise floor essentially originates in thermal noise and should therefore be uncorrelated. To estimate such correlations, we can simply study the product (5). xiy∗ j−xiy∗ j≃ˆ xiˆ yj ∗δnHxiδn∗ Hyi(5) Then, using local differences in (6), it is also possible to observe cross-correlations, while minimizing drift impact. diff(xi)diff(y∗ j)≃2∗ˆ xiˆ yj ∗δnHxiδn∗ Hyi(6) After the extraction of noise levels and correlations based on repeated measurements of a reflect standard, we could only reliably determine the correlations of high-level noises on the same port, between the reflect and the reference receivers. Fig. 2: Standard Deviation of noise propagated to S-parameters (30dB attenuator) In Fig. 2, we compare the actual standard deviation of a 30dB attenuator’s S-parameters due to noise with the predicted standard deviation using our model. Noise error was propagated to the first order using (4) with and without correlations, respectively labeled ”Jacobian” and ”Jacobian no corr”), which proved similar to the results obtained using METAS UncLib [17]. By taking into account the cross-correlations between the receiver chains, the model predicts well the expected noise variances for reflection coefficients. However, since we neglected correlations in transmission parameters, the impact of noise remains overestimated. B. Uncertainty Due to Drift Drift errors in on-wafer measurements originate from various sources, from the temperature-dependent behavior of Vector Network Analyzers’ receivers and frequency extenders’ mixers to stability of the probe station. Time-dependent drift depends on many factors, such as room-temperature stability, the power at which multiplier chains are driven, and the airflow around the probe station. As a consequence, an additional experiment is required to understand the magnitude of such deviation under conditions similar to those of the measurements. At lower frequencies, corresponding to the use of coaxial connectors, it is possible to monitor drifts thanks to an electronic calibration (e-cal) by performing a repeatable calibration regularly over a period of several hours. By this means, it is possible to obtain the variance of each of the error terms, eventually including their covariance by defining a multivariate normal Wiener random process [18]. However, as is the case in on-wafer measurements, when no e-cals are available and making a repeatable connection can be an issue, measurement drift can only be assessed with the regular measurement of a unique standard. The description of a possible method can be found in [6], where after an initial calibration, the measurement of a simple flush-through standard (thru) is used. The consequence of drift can be then expressed as in [15] with 3 a perturbation matrix positioned at the calibration plane, using uncorrelated multiplicative terms in transmission and additive terms in reflection. On-wafer, a thru can only be realized when probes are placed on the wafer, which means that these measurements will also be subject to vibration of the probe station, therefore progressively affecting the initial good contact. Thus, we propose the use of the more reliable probe open for drift quantification, after an initial probe tip calibration. As we only access the drift of the reflection coefficient at each port, it is not possible to understand exactly which error term is responsible for the drift. Therefore, we chose to use conservative estimates for the drift of the instrument as described for port pby equations (7), (8) and (9), similar to those for a flush-through in [6]. However, we should be aware that such a reflection-based estimate will not be sensitive to LO phase drift present in the transmission coefficient phase, which can be tracked based on the measurement of any reciprocal structure. ∆pp =max(|ℜ(Spp(tj)−Spp(ti))|,|ℑ(−)|)(7) ∆app =arg(Spp (tj)/Spp (ti))(8) ∆mpp =|Spp (tj)/Spp (ti)|−1(9) Based on these values, it is possible to define several Wiener random processes Wbased on a time-dependent normal distribution Nas described by Eq. (10). A classical approach to the estimation of αuses local differences applied to the process studied. W(α, ∆t) = N(0,Σ = α∆t)(10) Yet, as seen in the previous part, the reflection coefficient is subject to trace noise and so will the local differences. To reduce the influence of high-level noise, we first convolve the delta functions of size Nwith a Keizer window of size M, obtaining f ∆k. Finally, a better estimate for αis found in Eq. (11), which uses observed differences between any time tiand tj. For clarity, the formulation describes the probability that the delta function is under a certain time-dependent value described by the Wiener process, independently of the starting time. ˆαk=1 ΣiΣjwij N−M X i=0 N−M−i X j=M wij f ∆k(ti, tj)·f ∆k(ti, tj) ti+j−tj (11) wij =1 N−M−j(12) The drift error is then expressed at the calibration plane using the perturbation matrix (13), where δW represents independent random variables obtained from the Wiener process of (13) with αextracted on each of the delta functions defined in (7), (8) and (9). Dp(t) = δWpp1(1 + δWmpp1)ejδWapp1 (1 + δWmpp2)ejδWapp2δWpp2(13) To verify that the different methods are indeed equivalent, we performed a drift study using an e-cal, where we apply the different methodologies mentioned above: using a pair of shorts, a pair of opens, a thru line and finally, by directly tracking the error terms through recalibration. A circuit board was also introduced at one of the ports to affect the error box and be able to observe a more significant drift. The results are plotted in Fig. 3 for various verification structures. Under these conditions, we clearly observe an overlap of the curves for different extraction methodologies, which gives us some confidence in their equivalence. Finally, our estimate ˆαgives satisfactory results, enclosing the observed drift error within the 95% uncertainty bounds. Fig. 3: Estimated drift 95% bounds (continuous) against actual error (dashed) on e-cal verification standards C. Wafer Processing Inaccuracy The manufacturing of calibration structures and devices onwafer is also subject to non-idealities, responsible for part of the measured S-parameters’ uncertainty. To understand their impact on calibration, we performed a dimensional characterization with a confocal microscope on some of the calibration standards. The first results are already described in Table I of our original paper [7]. On the basis of these measurements, a covariance matrix representing the dimensional variation across the entire wafer was determined, which, with the CPW model of [19], could be propagated to the propagation constant γand line impedance Z0of the transmission line. These were then used in order to understand the uncertainty that affects the transmission lines of the calibration kit and the access lines that connect probes to the different standards and DUTs. Finally, for the sake of simplicity, we did not consider the dimensional variation of the pads and tapers. Fig. 4: Measured propagation constant of ”L1” and ”L3” samples, 95% error bars indicate the portion of error originating from geometry variation As shown in Fig. 4, a significant spread can be observed in the electrical properties of the transmission line between the different samples, where the propagation loss for ”L1” sample 4 Fig. 5: Measured impedance of ”L1” and (current) ”L3” samples, 95% error bars indicate the geometry variation on model fitted with method of [24] is more than twice the loss of the currently studied ”L3” sample, for the exact same dimensions. In fact, the .39µm gold deposited on the InP wafer for L1sample was measured with an average conductivity of 11.1 MS/m, while L3sample was measured at 21.0MS/m. This indicates that the leak in the reactor observed during fabrication had a significant impact on the properties of conductors. A more detailed description of the fabrication procedure can be found in [20]. In Fig. 5, the line impedance extraction using [21] and [22] did not match the method of [23]. This can be explained by the presence of a continuous gradient of conductivity within the metal strips, which, similar to surface roughness [24], affects local current densities and abnormally increases conductor inductance outside of what is expected from a uniform material. Proper modeling of the effect however requires information on the conductivity profile [25]. We therefore chose to model the uncertainty as a variation around the value extracted using [21]. III. UNCERTAINTY PROPAGATION For first-order propagation of uncertainty, we computed a Jacobian matrix around the solution for each step of the multiline-TRL algorithm, where uncertainty sources are considered as perturbations either located at the calibration plane or at the measurement plane. At first, we derived the equations for the propagation constant from the NIST’s implementation of [9] and [26]. Then, the error terms are determined based on the general algorithm of [10]. Finally, the displacement function that we introduced in [7] is included in the complete calibration scheme. A summary of the propagation of uncertainty is presented in Fig. 6, where uncertainty sources and measurements are presented at the input of each of the calibration steps. A. Propagation Constant To compute the uncertainty of the propagation constant, we first introduce the different parameters that affect the accuracy of each line’s measurement. As in [13], we apply noise directly at the receiver level in Eq. (14) with the matrix δBmand δAm the band a-wave matrix perturbated by noise, e.g. δBm= Bm+Nb,δAm=Am+Nawith Naand Nbthe matrices containing noise errors applied on each receiver as introduced in Eq. (1). Fig. 6: Procedure for uncertainty propagation The repeatability of probe contacts is introduced at the calibration plane level in Eq. (15) with the cascade matrices R1and R2described in [7] and s2tthe function that converts Sto T-parameters. The influence of wafer processing on the propagation constant will be considered later when computing error terms. Mi=s2t(Sraw) = s2t(δBmδA−1 m)(14) Ni=R1(x1i)·e−γli0 0eγli·R2(x2i)(15) In [26], a best common line is found and, forming a pair with each of the other (N−1) lines, we obtain (N−1) eigenvalues from which an optimal weighted average is found. In the original algorithm, the optimal weighting was found on the basis of uncertainty expressed by a perturbation matrix introduced at the calibration plane. It was already demonstrated there that at the first order, the determined propagation constant is only sensitive to changes in the transmission coefficient. In our conditions, since repeatability is also expressed at the calibration plane, only the transmission of our displacement matrix is of importance. Mk=Mij =MiM−1 j(16) Nk=Nij =NiN−1 j(17) We therefore define the problem of N−1eigenvalues based on the measurements of each pair of lines Miand Mj, which contain noise error, and the modeled lines Ni and Nj, which contain repeatability error. As measured Mij and modeled Nij expressed in (16) and (17) both share the same eigenvalues, we can express (18), where λMdesignate the eigenvalue after the correct root choice. Λ(Mk) = Λ(Nk) = λMk =1 2(λ1+1 λ2 )(18) Given our notation of (19), we reformulated the firstorder expressions obtained by Marks in [26] by including the random part of the transmission coefficient of the displacement matrix located at port p of each line i,δSrpi 21, giving (20). gNk = (lj−li)γk=ln(λNk)∼ˆ gNk +δgNk + ∆gNk (19) 5 δgNk =δSr1j 21 +δSr2j 21 −δSr1i 21 −δSr2i 21 (20) This formulation explicitly gives the uncertainty formulated at the calibration plane. We can therefore write the Jacobian matrix of (21) linking gNk to all real and imaginary terms of the left and right repeatability matrix. We did not include the reciprocity directly here for convenience, as they are taken care of by the correlations of the repeatability terms. JN(32×2) =1 20−I2−I20... I2I20T(21) δgNk =JN·δxLNi δxLNj (22) Also, knowing the position of the probes for each touchdown, it is possible to correct each of the gNk by subtracting the bias introduced in Eq. (23), where ∆Srpi 21 designates the systematic part of the transmission coefficient of the displacement matrix. ∆gNk = ∆Sr1j 21 + ∆Sr2j 21 −∆Sr1i 21 −∆Sr2i 21 (23) For the second step, it is now necessary to account for uncertainty expressed at the measurement plane, e.g., on the eigenvalues of Mij. For this purpose, we numerically compute the Jacobian matrix linking gMk to the perturbations accounted for in (14), giving the Jacobian matrix JMk used in (24), where δxLMi consists of a vector of 16 terms representing the impact of noise on the raw S-parameters of the line Li. δgMk =JMk ·δxLMi δxLMj (24) Combining the different equations results in Eq. (25). gk≡gMk +δgMk +δgNk −∆gNk (25) As there is no correlation between the variances of the measured and modeled propagation constants, as suggested by Hatab in [12], we can directly sum these two contributions to the total uncertainty in (26), with (27) and (28). Vgk =VMgk +VNgk (26) VMgk =JMk VMLi0 0VMLjJT Mk (27) VNgk =JNk VNLi0 0VNLjJT Nk (28) with VMLi=δxLMi δxT LMi =Vn16×16 (29) VNLi=δxLNi δxT LNi =Vr18×80 0Vr28×8(30) The propagation constant is then estimated on the basis of the weighted average of the different results gk. In (31), we reiterate the clear formulation from [9], where the optimal solution is calculated based on the vector of line length differences L and the covariance matrix V is obtained based on uncorrelated noise with similar variance introduced on the real and imaginary part of the transmission coefficient of the lines. However, this first formulation is not compatible with the covariance matrix that we computed so far, as we are manipulating a covariance matrix of complex values. γ=LTV−1G LTV−1L=C1×N−1·G(31) A convenient way to represent a complex number z, as introduced in Appendix C of [13], is the 2x2 matrix (32). z=ℜ(z)−ℑ(z) ℑ(z)ℜ(z)(32) By substituting this matrix for each term of the Land G vectors, introducing the L′and G′matrices, it is therefore possible to recompute with (33) the best linear unbiased estimator, though this time including the complete covariance matrix VG. ˆγ=C′ opt ·G′= (L′TV−1 GL′)−1·L′V−1 GG′(33) The covariance matrix of the propagation constant gamma can then simply be obtained with (34), by substituting the representation of (32) for each element of the original C, obtaining C’, or by using the updated version C′ opt. VγM=C′·  Vg1... Vg1/gN−1 ... ... ... VgN−1/g1 ... VgN−1  ·C′T(34) Within this matrix, since a common line is used to determine the covariance matrix VGof the vector G′, it is necessary to account for the correlation between the different terms, introducing the matrices Vgk/gm . We can again split the calculation since the modeled and measured contributions are not correlated, giving (35). Vgk/gm =JMk ·VMkm ·JT Mn +JNk ·VNkm ·JT Nn (35) With VMkm and VNkm computed as in (36) by placing the parameters of the common line Lisystematically first. VMkm =VMLi 0 0 0(36) B. Error Coefficients Now that uncertainty is expressed for the propagation constant, we apply a similar methodology within the context of [10], though this time using the complete covariance of the lines described in (37), with the matrix Vtthe covariance of the transmission line parameters obtained in section II-C. VLi =    Vn16×16 0 0 0 0Vr18×80 0 0 0 Vr28×80 0 0 0 Vt4×4    (37) Uncertainty contributions are introduced on both sides of the calibration equations (38) coming from [10], [7] where E contains the first six error terms, and k the symmetry factor determined during the next algorithm step of section III-C. Repeatability and topology are applied to the modeled response of the standards contained in the matrix N, while 6 noise terms are applied at the measurement plane, thus on terms appearing in N,Gand H. NE =G+kH (38) Before including the displacement terms, the ideal response of each line (39) is subject to the propagation constant error, composed of the sum of a measured part δγM, obtained in section III-A, and the process-related variation δγTof section II-C. Finally, we also introduce δΓLn = (δZLn − ZL)/(δZLn +ZL)that describes the line mismatch caused by changes in the cross section of the line, and fRthe function cascading the displacement matrices to S-parameters equivalent to Eq. (15). SLn =fR δΓLn e−(γ+δγM+δγT)ln e−(γ+δγM+δγT)lnδΓLn  (39) The covariance of the propagation constant Vγis then included in the complete covariance matrix (40) next to those of each line standard VLi . Vγ+L=    VγVL1/γ . VLN/γ Vγ/L1 VL1.0 . . . . Vγ/LN 0. VLN    (40) Since each line measurement is reused in this section, there is an existing correlation Vγ/Ln between the mean propagation constant and the modeled line transmission parameters. These can be calculated using (41), where we introduce the correlation between the extracted gkand the parameters of each line in (42). We then determined the expression at the measurement plane XMand at the calibration plane in XNin (43). VLn/γ2×52 =VT γ/Ln =δγδxT Ln =C′  δg1δxT Ln ... δgN−1δxT Ln  (41) δgkδxT Li =XM2×16 XN2×16 02×4(42) XM2×32 =δgMkδxT LMn =JMk "δxLMiδxT LMn δxLMj δxT LMn#(43) Finally, we calculate the Jacobian matrix JElinking the parameters considered with the vectors Xand Y, solutions of the least squares problem of (38) as described in [10] and [7]. We therefore obtain the covariance matrix of the terms of the Xand Yvectors in (44). VE24×24 =JEVγ+LJT E(44) C. Symmetry Factor k To finally obtain the covariance of all the error terms, we again separate the terms acting at the measurement plane and at the calibration plane. In addition to the ideal reflect response, we consider the mismatch and change of propagation constant of the transmission line sections connecting the reflect standards to pads of length lofs, using the local process covariance described in section II-C. As the exact reflect behavior ΓRis approximately known after calibration, we introduce the ratio (45) that defines the variation of the reflection coefficient with respect to changes in this access line. For each port, δΓ1,δΓ2represents the line offset mismatch and δγ1,δγ2the change in propagation constant. ∆Γ=Γb Γa =δΓ2+ΓRe−δγ2lofs (1 + δΓ2e−δγ2lofs +o(δΓ2 2)) δΓ1+ΓRe−δγ1lofs (1 + δΓ1e−δγ1lofs +o(δΓ2 1)) (45) The ideal response of the reflects is then modified in (46) by considering probe contact repeatability and asymmetry, which also modifies the matrices ∆G,∆Hand N′in (22), (23) and (24) of [7]. C11 =La11 +L2 a12Γa[1 + La22Γa+o(LkiiLkjj)] C22 =Lb11 +L2 b12Γa·∆Γ[1 + Lb22Γa·∆Γ+o(LkiiLkjj)] (46) Noise is added to the scheme similarly to the lines in the previous section. We therefore obtain the covariance matrix for each of the reflect standards VRi56×56 , as an additional processing accuracy covariance matrix of shape P4×4was included to consider the change in propagation characteristics of the left and right offsets. Finally, a Jacobian matrix is again computed tracing back the sensitivity of the final error terms to the previously computed X and Y vectors and to the uncertainty sources considered in the reflect standards. With (47) the complete covariance containing noise, drift, process variation, and repeatability for each of the reflects, we finally obtain the covariance of the error terms with (48). VE+R=    VE0.0 0VR1.0 . . . . 0 0 . VRM    (47) VK14×14 =JKVE+RJT K(48) D. Displacement Function In [7], we demonstrated the extraction of a displacement function linking the variation of S-parameters of the transitions to the position of the probes. This displacement function is also subject to the uncertainty sources previously discussed, essentially reducing possible result enhancement. We now propagate uncertainty in the parameters of this function. For a port 2 probe, these parameters are fitted on the basis of the matrix (49), with T′ raw designating the switch terms corrected measurements of any 2-port transmitting device when the probe is placed on its pads at a position pi. R2(xi) = L−1T2T′−1 raw0T′ rawiT−1 2L(49) Because the extraction method is insensitive to the exact S-parameters of the on-wafer device, and the pad geometry variation was not taken into account, process variation was not considered. Then, we chose the microscope resolution of σx=.5µm as the indeterminacy of the position of the probe. Vf=  Vp0.0 . . . 0. Vpn (50) 7 The complete covariance of the input parameters can therefore be written as in (50), where Vpi designate the covariance of the measurement taken for each position of the probe (51). Vpi =Vn16×16 0 0σ2 x(51) These input parameters propagate into the fitting parameters βof function using (52). Vβ=JβVfJT β(52) Finally, we obtain the covariance matrix Vrof the displacement function in (53) considering its parameters βand their covariance Vβ, the measured position xof the probe and its indeterminacy σx. Vr8×8=f(β, x, Vβ, σx)(53) E. Correction The last step to obtain the uncertainty of the DUT Sparameters is the propagation during correction of systematic errors, which can be equivalently described in the classical 8-terms error model [16], or more easily in [13]. Using the latter, we obtain Eq. (54) that accounts for noise, drift, and displacement function. Finally, δK,δM,δL and δH describe the error terms with covariance, fRand fDare functions cascading repeatability and drift matrices at the calibration plane. fD(fR(SDUT )) = (δKδBm−δMδAm)(δLδBm−δHδAm)−1 (54) Finally, we propagate these uncertainties on the basis of a 8×62 Jacobian matrix and the complete 62 ×62 covariance matrix of the uncertainties considered. IV. RESULTS AND DISCUSSION The different algorithms were implemented in Python with the scikit-rf library [27]. At first, we validated the approach by comparing the results with other available methods [12]. Then we applied the algorithm to the measurements from [7], where a more complete uncertainty evaluation was performed. Finally, since we computed the Jacobians during uncertainty propagation, complete uncertainty budgets are quickly obtained. A. Comparison Between Different Implementations To verify the validity of our approach, we compare the results obtained with our approach with a second implementation that uses automatic differentiation with METAS unclib [15], and with the version introduced in [12]. In Fig. 7 we plot the difference in the propagation constant between the common line pairing method from the MultiCal implementation of [26] and the method introduced in [11]. Additionally, the standard deviations obtained after first-order propagation of the noise error are sensitively the same for all three versions. This shows that in terms of first-order moment and hence for a small error, there is no difference between the two approaches. Fig. 7: Standard Deviation of Propagation Constant due to Noise and Difference of Estimates Obtained with TUG and NIST Algorithms The more interesting part concerns the computation of error terms using the general formulation of Silvonen [10] (TCX “transmitting circuit - any circuit - unknown circuit”) rather that using eigenvector-based methods, as in [26] and [12]. In [7] we already observed that the results obtained with the TCX algorithm were closer to the results obtained from the optimal mTRL of [8], where they already demonstrated a reduced uncertainty compared to the MultiCal implementation. Fig. 8: Calibrated 288um Line behavior (continuous) and standard deviation due to noise only (dotted) In Fig. 8, we compare the results obtained with the algorithm of [12] with our version of the multiline-TRL algorithm. Again, our results were also tested against automatic differentiation. The results confirm a slight reduction in standard deviation of the reflection coefficients when using the TCX approach instead of eigenvector-based methods of [28] and [7]. As expected, there is a more significant reduction in uncertainty of the transmission coefficient, since the Thru line was not considered ideal in our approach. This is different for other versions of the multiline-TRL, as they are all enforcing the calibration plane to be placed exactly in the middle of the Thru. It has the advantage of clearly specifying the calibration plane in reference to the probe position but at the cost of transmission measurement accuracy. This is especially true in on-wafer measurements, where parasitic modes participate in the degradation of the transmission of short lines [29]. B. Displacement Function Next, we analyze the accuracy of the extracted displacement function. As described in section III-D, we only considered position indeterminacy and noise influence on the parameters of the interpolation function. In the plot Fig. 9, the total uncertainty bounds of the displacement function’s S-parameters for a probe displacement of 8µm were obtained for each of the following scenarii: when error terms are obtained by (i.) applying the nonlinear optimization method (omtrl) or (ii.) the TCX approach, and when interpolation is based on (a.) the polynomial method or (b.) the lossless method. 8 Fig. 9: Displacement Function S-parameters (x= 8µm) with 95% error bars plotted for ii.a. and ii.b. The results first show that there is no significant difference between the results obtained with nonlinear optimization (omtrl) and the general algorithm (tcx) both in value and in 1st-order moments. A notable difference between the two interpolation methods is that a bias appears for the reflection coefficient angle at low frequency. This is related to the wellknown change in behavior of coplanar waveguide lines at low frequencies when the skin depth approaches the conductor thickness. Although, since this phenomenon occurs at relatively low frequencies, the corresponding magnitude of the reflection coefficient introduced by a displaced probe is quite low, and thus has a limited impact on compensated results. A more important difference between the two interpolation approaches is the magnitude of the transmission coefficient, which is clearly more uncertain for the polynomial approach, as the reflection coefficient in the lossless approach is mainly dependent on reflection coefficients, i.e. on 1−|s11|2. Fig. 10: Uncertainty Budget on function (extracted from TCX) (continuous: polynomial method, dashed: lossless method) Finally, in the budget of uncertainty contributions in both polynomial and lossless functions Fig. 10, we separated the impact on the function parameters of the probe position variance caused by the image processing, labeled ”position”, from its direct impact during compensation, labeled ”indeterminacy”. Concerning the polynomial approach, we can observe that the variation seen on the transmission magnitude at high frequency is indeed not coming from noise, but is more likely due to non-local effects that were not taken into account. Additionally, the use of simplified models only marginally improves the transmission phase because the position indeterminacy remains the same, even though the function parameters improve. C. Uncertainty of the Propagation Constant We now compute the propagation constant and its uncertainty bounds estimated on the measurements from [7], this time with a different approach that is clarified by Eq. (19). The results shown in Fig. 11 first show a reduction in the effective permittivity interband discontinuities after compensation. Fig. 11: Propagation constant with 95% bounds Although the eigenvalue and optimal approaches are giving similar results, we observe more significant ripples in the optimal approach. This signifies that the optimal approach is more prone to uncertainty than the eigenvalue one. In addition, in the optimal version, residual-based uncertainty quantification using (55) seems to underestimate the uncertainty of the effective permittivity, showing the limits of such a general formulation. cov(β)∼P|rij k|2 n−m∗(JTJ)−1(55) Finally, since we did not account for probe-to-probe crosstalk, parasitic probe effects, and the random part of probe contact repeatability, the obtained budget is still incomplete. This is quite visible in the real part of the propagation constant, where the eigenvalue approach underestimates uncertainty. This is particularly visible in the WR3.4 band, where one of the probes was damaged and was responsible for significant crosstalk and parasitic mode injection. Consequently, only the general residual-based estimation was able to capture the effect. D. Uncertainty of Calibrated Measurements The test of the proposed propagation algorithm was verified with numerous verification standards already described in [7]. To remain concise, we focus on one of the offset load standards, and compare uncertainty propagated through the TCX and the optimal multiline-TRL (omtrl) approaches. The calibrated measurements with and without probe position compensation are shown in Fig. 12, where the 95% confidence intervals include noise, probe placement, and geometry. Drift was only taken into account in the 0-220GHz frequency band since the experiments were not yet properly setup at the time of the measurements. We realize that compared to the values obtained in [7] for the optimal calibration, the standard deviation slightly 9 Fig. 12: Calibrated offset load with 95% bounds (top to bottom plots: mag uncompensated; mag compensated; angle uncompensated; angle compensated) increased, explained by the displacement function, which now accounts for the uncertainty in its parameters and includes additional sources other than the indeterminacy of the probe position. There is a relative agreement between the two approaches apart from the following issues: the underestimation of uncertainty in transmission magnitude for similar reasons than for the propagation constant’s real part; and a more significant variance at low frequency for the TCX approach, since we accounted for process variation. Finally, this is a passive standard, and hence the difference between the forward and reverse transmission coefficient’s angles is attributable to a relative change of the LO signal phase between port 1 and port 2. This could have been caused by the movements of the cables connecting the frequency extenders or drifts that were not taken into account in the WR3.4 and WR2.2 bands. Fig. 13: Uncertainty budget of S-parameters’ magnitude (top: uncompensated, bottom: compensated) In Fig. 13, we finally retrieve the uncertainty budget for the offset load for transmission and reflection magnitudes, both for the optimal and TCX versions. Clearly, the agreement between the two methods in the uncompensated case comes from the displacement function, labeled ”rfunc”. After compensation, most of the error in reflection at lower frequency is attributable to process variation, labeled ”topo”. In transmission, the optimal algorithm residuals, labeled ”residuals”, indicate that other uncertainty sources are present and not accounted for, like probe crosstalk for example. Also, it is noise and drift that are responsible for ripples seen at the edge of each waveguide band, when the VNA port power is reduced and error-terms deteriorate. Finally, we understand that each methodology has its strengths, but a complete uncertainty study is only possible if all parasitic modes are avoided, or modeled and quantified. In that context, a partially damaged probe or bad cable management can lead to an unpredicted measurement error. This makes quantification of uncertainty in on-wafer measurements a particularly difficult task, especially on a dense wafer, often encountered in the industry. Fig. 14: Extracted capacitance of an open standard with 95% bounds (green: port 2 uncompensated, blue: port 2 compensated) Meanwhile, the methodology we propose permits a more accurate characterization of small devices. In fact, as demonstrated by Fig. 14, there is a clear reduction in the spread of extracted lumped elements across frequencies. V. CONCLUSION In on-wafer S-parameter measurements, comprehensive uncertainty budgets are of importance to help understand the main sources of uncertainties. In the case where calculated bounds are obviously not representative of the error observed on verification structures, it can help identify that the source of error is not originated in test fixtures or probe contact repeatability. Thus, it is most likely an issue linked to the design of the calibration set, a too dense probe environment, or that probes are injecting unwanted propagation modes. Such a comprehensive study is necessary in the realization of primary calibration, and thus in on-wafer measurements, it is absolutely critical for obtaining reliable data. It has been seen more than once that a portion of yield loss is caused by unreliable measurement data. For that purpose, we enhanced and verified a few methodologies found in the literature to permit the in situ quantification of test fixture noise and drift, e.g., when on-wafer probes are still mounted on the frequency extenders, simplifying the extraction of such quantities. We then proposed an approach to the quantification of uncertainty in displacement functions used in the compensated algorithms of [7]. Although it is possible to understand the effect of probe misplacements