Compensating Probe Misplacements in On-Wafer S-Parameters Measurements

As the maximum frequency of electronics is rising, on-wafer measurements play an important role in modeling of integrated devices. Most of the time, due to the lack of measurement accuracy beyond 110 GHz, such models are usually extracted at frequencies much below their working frequencies and are subsequently extrapolated. The validity of such models is then mostly verified after fabrication of the complete chip, with a simple pass and fail test. This is stating the necessity of enhancing measurement results by any means possible, i.e., to reduce the overall uncertainty in such measurements. It is widely accepted that one of the main sources of uncertainty in such measurements is probe contact repeatability, since it is difficult to reach position accuracy below a few micrometers. We are presenting in this article a method to model the $S$ -parameter variation with probe position on the pads, which can then be used to either estimate contact repeatability uncertainty or further enhance measurement results. The approach is validated based on the measurements performed at 500 GHz.

The performance of vector network analysis was significantly enhanced in recent years, thanks to the introduction of zero bias Schottky diodes in mixers and multiplier chains' part of the frequency extension modules [1]. In parallel, the probing technology also got further enhanced [2], reaching nowadays up to 1.1 THz, owing to the micromachining process.
Since most calibration algorithms suppose that only one or two modes propagate in the transmission lines, the unshielded environment of on-wafer measurements makes such calibration sensitive to many more uncertainty sources. These include additional unexpected propagating modes [3] and evanescent modes around the probe transition [4]. In addition, any change in the probe behavior between calibration standards and device under test (DUT) will also degrade the measurement accuracy. Such changes may be caused by parasitic coupling to neighboring structures [5], or more obviously, by probe contact variation on the pads making the transition, i.e., probe contact repeatability.
A lot of research has been carried out recently to further improve the accuracy of probe positioning [6] and probe planarity [7], permitting a clear enhancement of measurement accuracy at high frequencies [8]. Automated techniques require a well-calibrated x yz stage and may not be applicable to manual probe stations.
In addition, probe contact repeatability is usually estimated by repeated measurements on various calibration standards [9]. To obtain a statistically valid uncertainty evaluation, at least a dozen of touchdowns are usually necessary. This describes the type A approach to evaluate contact repeatability uncertainty [10].
By evaluating the variation in the transition behavior with the probe position, we are aiming to model the systematical part of contact repeatability, composed of the probe mispositioning error. This is not meant to replace the statistical study which is still necessary to understand the random effects, linked to contact resistance, tilt, and probe height variation. As optical microscopes are usually used to visualize the contact, it is possible with simple image processing techniques to obtain the final position on the pads.
In Section II, we first present a method to evaluate the probe transition variation with respect to the probe position on the pads. Then, in Section III, we introduce a calibration algorithm permitting to further enhance the measurement results by compensating the probe position in the measurements. Finally, in Section IV, we present the measurement results when This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ applying this new method, before discussing its validity in Section V.

II. MODEL DESCRIPTION
Based on the multiline thru reflect line (TRL) calibration, S-parameters correspond to the traveling wave effectively propagating in the transmission line [11]. Therefore, the Sparameters are normalized to the line impedance, which is a priori unknown. Several methods exist to determine the transmission line impedance, either based on the propagation constant and the hypothesis that the capacitance remains constant over the whole frequency band [12], [13] or using model-enforced calibration comparison [14]. While the first method may only be valid in specific cases, the second method heavily relies on well-defined error terms. As a consequence, we would benefit from a general method without the necessity to specify any impedance, i.e., by only using S-parameters.
In [15], we were aiming to model such variation based on lumped elements. While such model requires a known line impedance, it is quite limited to the number of parameters that would need to be fit, easily falling into Occam's razor caveat, and still remains an approximation. Another method described in [16] also proposed the use of a simple stub model based on the line propagation constant and impedance.
The following model only presumes that there is one mode propagating and that the change in transition behavior with the probe position is a local phenomenon. While the unimodality hypothesis is similar for most calibration algorithms, the locality of the phenomena should be verified with a sufficient distance to the calibration plane. From the multiline TRL calibration, we first obtain the following relationship between switch terms' corrected measurements T raw , left/right error boxes T A /T B , and calibrated results T cal : We can then introduce the variation in the transition behavior as two cascade matrices localized at the left/right probe pads R a /R b . Using the propagation constant, we can numerically relocate the phenomena by deembedding them from the calibration plane to the pads with L. This permits the reflection components' phases to have a more stable cross-frequency behavior and provides much better agreement with the simulated results. In addition, since no correction was applied yet, T A and T B of (1) are only approximating the actual error boxes T A and T B of the following equation: (2)

A. Estimation of R Matrix
It was demonstrated in [15] that during measurements, it is possible to track down the variation in the transition's S-parameters with respect to the probe position. This can be achieved by letting one of the probes fixed on any line standard, and moving the second probe at various positions p i . For reference, herein, we chose to move the probe on the right side. By forming (3), the steady term R a cancels. Under the locality and unimodality assumptions, we can therefore conclude that the exact position of probe port 1 has no influence on the extraction method The reference position p 0 can be chosen arbitrarily, where we set R b0 = I . Yet, T A and T B are not exactly known, and as a consequence, we can only approximate them in the first instance. We can relate the estimates and actual error boxes by the following equation where the position p c designates the approximate position of the first calibration, and R bc designates the bias introduced in this experiment when p 0 = p c : While the position p c is only approximate, we can expect it to be close to the probe position on the Thru line of the initial calibration. Hence, to ensure that R bc → I , a good choice of reference position p 0 is the closest one to the probe position on the Thru.
We now need to interpolate R est (x, y) with the probe position on the pads. A first method consists of simply imposing reciprocity condition, fixing the number of complex variables to interpolate at three. We can then fit a polynome to the nth order to each of these variables, i.e., S 11 (x, y), S 12 (x, y), and S 22 (x, y).
While this general approach seems to already give good results, it may be useful to further reduce the number of variables to fit by imposing additional conditions to this matrix, to increase its robustness to uncertainty.

B. Lossless Condition
Connector repeatability error is usually represented as a small lossless symmetric perturbation [17]. Thus, retaking our previous formulation, S 11 = S 22 . While this may be true in connectors, it is possible to find a unitary matrix factorization much closer to the actual phenomena of probe movements along the pads, as illustrated in Fig. 1. Such factorization may be written as in the following equation: With φ(x), θ(x), φ a (x), and φ b (x) being fit parameters according to the position. More specifically, for this matrix, lossless condition additionally enforces φ a = φ b . The simulation results in Fig. 2 imply that with small displacement, the lossless model is very close to the actual phenomena. Though in the case of significant displacement, to further reduce the interpolation error, it is possible to release this last constraint by considering φ a = φ b , introducing our lossy model. Imposing lossless condition also avoids additional variation with the position of this matrix's transmission. As we will see in Section V, we can observe a dependence of the extracted model to the length of the line standard on which it was extracted on. We believe that this is caused by the nonlocal effects, certainly linked to the parasitic probe effects described in [18] and [4], which mainly affects the transmission coefficients.

C. Extension to Lateral Position
While current factorization may accurately model any displacement along the pads, it can also capture the variation in behavior with lateral movements if they remain lossless. However, slotline (SL) mode is expected to propagate when coplanar waveguide (CPW) grounds are not connected by bridges. In these conditions, the CPW mode may start to weakly couple to SL modes, with lateral mispositions.
The transition itself will then not consist of only a two-port matrix, but most likely of a four-port matrix, linking on both sides the modes propagating in the probe tips and the modes propagating in the transmission lines Ideally, we can expect the conditions described by (8) on the different 2 × 2 matrices forming the four-port S-parameters. While such relationships can be verified thanks to 3-D full-wave simulations, transmission between CPW and SL modes where expected only to be found below −35 dB in the worst condition, i.e., most forward position, most misaligned. While this may be the case in simulations, it is possible to find probes that certainly couple to the SL modes directly due to introduced asymmetry at the probe tips. This can be explained by various factors, from the most probable to the least: impurities sticking to the probe, probe aging, or even manufacturing inaccuracies.

III. POSSIBLE RECALIBRATION SCHEME
To find the actual error boxes T A and T B of (2), we can use a nonlinear least-square fitting method following developments made in [19]. After a first multiline TRL calibration, we already obtain good estimates of the propagation constant, left and right error boxes T A and T B , which is useful for the initial conditions of the optimization scheme.

A. Mathematical Approach
The algorithm consists of the minimization of a set of equations constituting the measurement model. Retaking (2) and now that we have an estimate for the displacement function, we are aiming to estimate T A and T B . The overdetermined set of equations are given by the left and right R-matrix the different lines defined by their propagation constant (9), and by the reflects only defined by their symmetrical termination as follows: Finally, we define the cost function to be minimized as the difference in raw measured and raw estimated S-parameters, which are defined by the seven error terms in S a and S b , propagation constant γ , and the symmetrical reflection coefficients of the reflects Lk Such minimization problem is then solved using the Levenberg-Marquardt (LM) algorithm. This is completed with eight cores in 58 s for 200 frequency points. In comparison, more classical minimization techniques approximately perform five times slower.
In Appendix A, we also developed equations that include repeatability correction to the general Thru-circuit-unknown (TCX) algorithm from [20]. The derivation in terms of S-parameters helps greatly to include our new terms to the first order. In comparison, this last algorithm only takes 10 s on a single core to be solved while giving similar results.
Finally, once the new error terms are found, by knowing the position of the probes during the DUT measurement, we can simply apply the inverse of the relationship (2), which leads to the following equation: The diagram of Fig. 3 describes the global behavior of the proposed calibration algorithm, where we represented the different fluxes of information necessary to its completion.

B. Residuals Definition and Uncertainty
To verify the validity of the approach, we first need to introduce few definitions since classical calibration comparison to the multiline TRL will only give a rough idea of the gain in accuracy.
As pointed out in [19], the main advantage of the LM minimization technique is that we can easily obtain a good estimate of the error terms' covariance matrix, based on the calibration residuals r i j k , the sensitivity of the objective function to each parameter J T J with J being the Jacobian matrix at the minimum, and the naive statistical degree of freedom n − m, with n the number of observations, and m the number of parameters. For clarification, each line gives four observations, while each reflect only gives two Using this first definition, we can finally obtain the error terms of the calibration with their uncertainty. With this covariance matrix, we can easily propagate error terms' uncertainty into the calibrated DUT S-parameters using the Monte Carlo simulation. However, to obtain comparable results in terms of propagated uncertainty, we will have to distinguish repeatability uncertainty in case 1, when we do not apply mispositioning compensation, and in case 2 when applying the full compensation algorithm.
In case 1, to include repeatability uncertainty to the uncompensated results, we have to remind that the current R-matrix was extracted with reference to T B . Thus, it is necessary to compensate for the shift of reference plane between T B and T B . We relate the R-matrix of case 1 R b to the R-matrix We can then reuse (4) to obtain the actual repeatability function in case 1 Finally, mispositioning uncertainty is propagated by injecting distributions in the R-matrix and form the following equation: In case 2, we would require a more complete study to model the uncertainty of the displacement function. A possible type A approach to the evaluation of the random part of probe contact repeatability is to perform the complete cycle of "calibrationfunction extraction-compensation" several times. Through the variation in the extracted function and the compensated results, we would consequently obtain an estimate of the uncertainties caused by the function and by the random part of contact repeatability. However, since we did not perform such measurements, we only considered the position indeterminacy linked to the microscope resolution, at around 0.5 μm. Similar to case 1, position indeterminacy is propagated using (16).
Another method, which does not require to define the uncertainty of the extracted function, is to compare the overall residuals of the LM calibration in both the cases. For the TCX approach, we used the same residuals' definition as in (11). The final residuals for both the cases in the LM and TCX calibrations can be found in Fig. 4, where we clearly observe a reduction in overall residuals when applying the correction.

IV. MEASUREMENTS A. Setup
The measurement setup consisted of a CM300x probe station (FormFactor) and two RPP504 x yz stages for the probe positioners. To reduce the drift impact on the measurements, very short 45 • sections of waveguide were used to connect the probes to the Virginia Diodes' frequency extenders, themselves placed at 45 • on Formfactor's mini extender arms. The network analyzer and the probe station were controlled via python scripts running on a laptop. For measurements in the WR3.4 and WR2.2 bands, 50-μm pitch Dominion MicroProbe It is important to specify that the setup was regularly mounted and dismounted, and the positioner axes were compensated using our own algorithms while the probe position was considered as steady. Unfortunately, due to limited time, we could not fully identify why the probe position was slightly drifting, even after long warm-up time. Several external sources of drifts can be the cause of such variation, including small variation in room temperature, where the probe position can have a dependence of several μm/ • C. This justifies Formfactor's automatic positioning algorithms (Contact Intelligence) that we, however, did not use during these measurements. As a result, we could only reach position accuracy of ±3 μm, which could have been reduced under ±1 μm with Contact Intelligence.

B. Position Detection
During each measurement, a snapshot from both the landed probes was taken. Based on these pictures, it was possible to use simple image processing techniques to find the relative positions of the probes with respect to the signal launchers. We first performed a bilinear interpolation to numerically augment the number of pixels, thus enhancing the accuracy of the detection below microscope pixel level. Then, to reduce the influence of the signal launcher structure on the detected probe position, we used different masks by filtering the image based on color ranges.
Finding the optimal position consisted of minimizing a cost function composed of the differences between the positioned feature and the corresponding mask. Since we did not require too much speed in postprocessing, we simply went through all the possible positions within a rectangle approximately enclosing the final position. An example of the output can be found in Fig. 6.

C. Calibration Substrate
We designed an extended multiline TRL set that was implemented with evaporated gold on a 400-μm-thick indium phosphide (InP) substrate. The CPW transmission line has dimensions of 6-μm signal linewidth, 5-μm gap, and 90-μm total width. The total width was chosen to avoid the additional radiation effects at high frequencies [22]. The signal and gap sizes were chosen based on the manufacturing constraints, while respecting the conclusions made in [4], i.e., W gnd ≥ 2W gap + W sig . Finally, the pads and tapers were designed to permit the use of various 50-μm pitch probes on such narrow coplanar structures. The size of the signal pad was set to 18 per 30 μm permitting a significant probe skating in case of tip planarity issues.
The calibration set consisted of 12 lines going from a very short Thru of only 180-μm length to a long line of 2484 μm, as well as offset reflects of various lengths. Next to this, we also implemented several verification structures mainly composed of mismatched lines of various lengths. Following the guidelines of [5], a significant spacing was used to avoid coupling to neighboring structures, thus using 750 μm of lateral spacing and 900 μm of longitudinal spacing. Table I describes the dimensions of some of the structures, which includes uncertainty caused by substrate's inhomogeneity, where W gap designates the gap width, W tot is the total width, W sig is the signal conductor width, T cond is the thickness of the conductors, L tot is the total length from edge to edge, and L eff is the effective length from the initial calibration plane. These values were determined by characterizing the topology of several structures across the wafer with a confocal microscope. It is still unclear at this point whether the obtained uncertainty values are inherent to the measurement methodology or to actual wafer inhomogeneity.

D. Crosstalk Effects and Anomalies
The use of shorter lines serves the verification of nonlocal effects when extracting the repeatability model and permits to visualize anomalies in short CPW transmission lines as discussed in [18]. Thus, after several trials, the shortest line effectively used as a Thru was the one of 258 μm, i.e., "Line108" by following our notation. We plotted the calibrated line transmission in Fig. 8, where we can clearly observe the aforementioned anomalies. Interestingly, unlike the Multical implementation of the mTRL algorithm, both the LM and TCX approaches do not assume the chosen Thru as being perfect.
Comparatively, we also measured similar structures with air bridges to understand whether such effects could be associated with SL mode propagation, but such phenomenon was sensibly similar. On the contrary, on the lower portion of the WR3.4 waveguide band, we were able to identify a significant injection of SL modes when measuring open standards with the probe on port 1, which mainly had an influence on the computation of the symmetry factor k (described in Appendix).
Astonishingly, even the 756-μm-long line is experiencing a similar behavior. The fact that such phenomenon is present at k × f c frequencies with k an integer may suggest that it is coming from some resonant mode. In fact, the probe tips of ≈400 μm are reaching one wavelength at frequencies around 250 GHz, thus giving maximum radiation efficiency. In Fig. 8, we were effectively able to measure a significant transmission of more than −30 dB in various reflects' standards. Isolation of such structures being normally much lower, their transmission is quite representative of the magnitude of probe crosstalk effects. We can clearly observe the correlation between crosstalk magnitude on the reflects' transmission and frequencies at which anomalies occur. This strongly suggests that probe crosstalk is at least partially responsible for such anomalies.

E. Results
Initially, after performing a multiline TRL calibration, we directly extracted the displacement model on several transmission lines using the method described in Section II. A first verification of the model validity was then performed by compensating measurements of various reflects' standards when the probe was voluntarily displaced on the pads (not shown here). While this seemed to give satisfying results on each individual standard by significantly reducing the S-parameters' variation, we could not see any enhancement in the behavior of extracted lumped element models of similar standards with various offset lengths. Such unsatisfactory results are explained by the undefined calibrated position p c in (4), i.e., by the approximate error terms of the initial calibration T A and T B . This yielded our new approach to the calibration, permitting to properly define error terms and reference positions.
To demonstrate the output of different calibration algorithms, we plotted in Fig. 9 the calibrated results of a mismatched line in which reflection coefficients are quite sensitive to probe displacements. In addition, in Fig. 10 we plotted the calibrated results of an offset serial load after moving the calibration plane to its center. These results clearly demonstrate an enhancement in terms of continuity between waveguide bands and reduce the discrepancies caused by inaccurate probe positioning. Finally, we can observe that the TCX algorithm results are very close to the LM approach in both the cases.
The error bars were computed based on the Monte Carlo simulations by only considering high-and low-level noises, calibration residuals, mispositioning in case 1, and position indeterminacy in case 2. Noises were considered uncorrelated from the other uncertainty sources. From (13), we obtained estimates of the error term covariance matrix based on the residuals, which could be directly used in a multivariate normal distribution. Finally, the random position obtained from the distributions of Fig. 7 in case 1, or from position indeterminacy in case 2, can be fed to the extracted displacement function, which consequently keeps the correlations between its output S-parameters. Table II summarizes the contribution of each considered sources to the computed error bars. In case 2, we computed the values for position indeterminacy of ±0.5 and ±1 μm. Yet, some important contributions are still missing, such as drifts, the random part of probe contact repeatability, and wafer inhomogeneity.
Finally, in Fig. 11, we are comparing the displacement function obtained from the measurements and from 3-D fullwave simulations (Keysight EMpro). The simulation results were obtained with an approximate probe EM model obtained from [16] that we updated based on microscope snapshots of the probes. All the probe materials were taken as lossless to reduce the differences in behavior between simulations. The results are astonishingly close to the measurements, proving that such set of simulations could reproduce the phenomena occurring in the measurements. This gives us an additional tool for the design of signal launchers less sensitive to probe displacements.

V. LIMITATIONS
When extracting the R-matrix on lines of different lengths using the method presented in Section II, we were able to identify a strong dependency of its terms with the line length, as observed in Fig. 12. In fact, on short lines, the locality hypothesis on which we base ourselves will be invalidated by the presence of evanescent modes and probe crosstalk, as the transitions are closer to each other. The fact that our displacement function is able to capture these effects is suggesting that the anomalies described in [18] are also probeposition-dependent. On excessively short lines, e.g., Thru and Line72 in our case, even the angles of the different terms are affected. The former suggests that the method presented here requires sufficiently long lines to respect the locality hypothesis.
Yet, even under these conditions, concerning Line216, the variation in transmission magnitude with the probe position shown in Fig. 13 is still difficult to capture accurately. In fact, it is first subjected to the anomalies mainly affecting the  transmission of lines, and to the random contact resistance. Thus, it is preferable to let the lossless model transmission only depend on the R-matrix reflection terms, i.e., θ(x) is estimated based on the reflection terms.
The current extraction method also relies on one measurement at the reference position, e.g., on T raw0 of (5). As a consequence, any error introduced in this unique measurement will then be propagated in the model itself. After a clarification of the uncertainty sources affecting this measurement, we would then benefit from a probabilistic approach to the estimation of such function. Main random effects that are expected to have a significant impact on the extraction of the function are probe height, nonlocal effects, and test fixture drifts.
Finally, as it was the case in the lower portion of the WR3.4 waveguide band, SL mode propagation also impacted this method, thus reducing the enhancement of error terms' definition.

VI. CONCLUSION
In this article, we successfully developed a new approach to evaluate probe mispositioning effects in on-wafer S-parameters' measurements. After studying the properties of the extracted matrix with accurate 3-D full-wave simulations, we were able to develop a model that does not require the definition of impedance or propagation constant at pad level. Thanks to the generality of the approach, with only small modifications, we expect such model to be applicable to a quite extensive variety of signal launchers. Furthermore, we demonstrated a possible application of this model by compensating probe misplacement in measurements performed up to 500 GHz. While measurements yielded very promising results by compensating measurement errors introduced by quite significant positioning error of ±3 μm, the method still requires a deeper uncertainty study to understand how far measurements could be enhanced.

APPENDIX ALTERNATIVE CALIBRATION ALGORITHM
Based on the calibration algorithm described in [20], it is possible to introduce the terms on the left and right displacement matrices R a and R b . To simplify calculus, we used the following notation: By following the signal flow charts when adding the displacement matrix on each side of an ideal line, it is first possible to obtain the following system, where A corresponds to the S-parameters of a displacement corrected ideal line:

A. Computation of Line Transmission
The following equation is equivalent to the eigenvalue problem from mTRL calibration, thus falling in the same singularity when the angle of the transmission between the Thru and the Line standards are within ±kπ: where M A and M B correspond to the raw measured S-parameters of lines of length l a and l b after switch terms' correction, respectively, and M a and M b their determinant. Manipulations similar than in the mTRL algorithm of [23] can be realized for choosing the correct root and obtaining the weighted least-square solution.

B. Error Terms Computation
To compute the first error terms, we deviate from the multiline TRL algorithm by correcting each of the lines. In fact, the eigenvalue problem on which this algorithm is based is this time ill-conditioned due to the mispositioning of the probes. Recognizing that the ideal lines are corrected by the repeatability function, cascade matrices cannot be considered as diagonal anymore.
Thanks to the expressions clearly derived in [20], a much simpler formulation for our lines can be found based on the S-parameters. Equations [20, eqs. (1)- (4)] can be reimplemented, this time using the corrected expressions (17) for each of the ideal lines. The propagation constant having been already calculated in the previous part, we can consider each line to be fully known. We can then simply use S 12 = e −γ l k .
Hence, this forms an overdetermined set of equations, and a least-square solution is found describing error boxes' ports 1 and 2 apart from the well-known symmetry factor. This set of equations benefits from being much less sensitive to errors than the expressions used in the original mTRL algorithm.

C. Symmetry Factor
The next step consists of solving the last symmetry factor, and this is where expressions had to be the most readapted. By considering that errors in reflection are sufficiently low, we can safely use the first-order expressions of (19), where C designates the S-parameters of a displacement corrected ideal symmetrical reflect To solve the system (20)-(26), an additional step compared with the original article has to be followed. First, we can solve the linear part of the equation using least-square method for both k = 0 and 1, giving, respectively, U and W vectors as solutions for E(k = 0) and E(k = 1). Introducing V = U − W , we can obtain a first estimate of using (27). Root choice is made based on the previously declared offset lengths and approximate values for the different reflect standards, e.g., min(|1 − roots(P( ))/ est e −2γ l ofs |) Once a valid estimate is found for , we can modify the second part of (25) where we use this last estimate to correct G and H using (28). Since both G, H 1, the final solution converges in very few loops Provided that a more close solution is found to estimate T A and T B , it is possible to further refine the estimate of the repeatability function, thus ensuring that R bc → I in (5). Again here, due to the same reason, the final solution converges in few loops.

ACKNOWLEDGMENT
The authors especially wish to thank Gavin Fisher (Formfactor) for his help in setting up the probe station and Virginia Diodes Inc. (VDI) for the frequency extender heads.