Self-calibration technique for characterization of integrated THz waveguides

Emerging high-frequency accelerator technology in the terahertz regime is promising for the development of compact high-brightness accelerators and high resolution-power beam diagnostics. One resounding challenge when scaling to higher frequencies and to smaller structures is the proportional scaling of tolerances which can hinder the overall performance of the structure. Consequently, characterizing these structures is essential for nominal operation. Here, we present a novel and simple self-calibration technique to characterize the dispersion relation of integrated hollow THz-waveguides. The developed model is verified in simulation by extracting dispersion characteristics of a standard waveguide a priori known by theory. The extracted phase velocity does not deviate from the true value by more than $ 9 \times 10^{-5} ~\%$. In experiments the method demonstrates its ability to measure dispersion characteristics of non-standard waveguides embedded with their couplers with an accuracy below $ \approx 0.5~\% $ and precision of $ \approx 0.05~\% $. Equipped with dielectric lining the metallic waveguides act as slow wave structures, and the dispersion curves are compared without and with dielectric. A phase synchronous mode, suitable for transverse deflection, is found at $ 275~\text{GHz} $.

Alternatively, recent developments in high-power multicycle THz production [13] support the direct use of THz waveforms in non-resonant structures, e.g., a dielectriclined waveguide (DLW). In DLWs, the structure characteristics are entirely determined by the transverse dimensions and material parameters of the structure: the inner radius a, outer radius b, dielectric permittivity of the layer r and conductivity of the metallic shell [14][15][16] -see Fig. 1(a). In a DLW, an efficient interaction between the driving field and charged bunch requires synchronicity between the phase (v ph ) and bunch (v e ) velocities. The phase velocity is determined by the characteristic equation of the DLW and depends on the frequency of the incident wave. In laser-based approaches, the central frequency of the produced waveforms depend on the crystal temperature [17], providing a way to control the phase velocity inside the structure. However, lower temperatures are beneficial to reduce THz absorption within the lithium niobate crystal. Therefore, with limited tunability, it is necessary to properly characterize a structure for efficient interactions.
Several well-established techniques exist for conventional radio frequency (RF) structures, such as the beadpull method [18,19] for measuring field profiles of cavities, * max.kellermeier@desy.de shorting planes for determining dispersion curves [20], and the coaxial-wire method [21,22] or the recently applied Goubau-line method [23] for impedance measurements of vacuum components. The small dimensions of mm-scale DLWs make these techniques challenging. Alternatively, optical methods have been used to measure the dispersion characteristics of a structure in the THz regime [24]. Here a photoconductive antenna was used to sample the electric field after the structure, providing the dispersive characteristics of the structure for various mode excitations. A limitation of this technique resides in the fact that it measures the integrated effect of the structure, i.e., it is not capable of providing local information within the waveguide, and includes coupler effects.
To eliminate the coupler contribution, two waveguides of different length, but equal cross section and same coupler geometry [25,Sec. 7.1.2] can be compared. Considering the phase difference between these two structures instead of the absolute phase shift, the phase velocity is deduced. This approach requires a precise knowledge of the length difference, and high machining precision for identical cross sectional and horn geometry. The waveguides cannot be characterized independently.
In this paper, we introduce an error network model to characterize the phase shift in integrated high frequency structures (∼ 200 GHzto300 GHz). The technique relies on a movable reflecting device within the waveguide, in analogy to a non-contacting variable short. The measured phase shift is used to deduce the phase velocity. The paper is organized as follows: first, the self-calibration method is introduced. The following section demonstrates the applicability in simulations. Then, the method is applied to experimental data. Finally, the results and limits of the error network model are discussed.

II. REFLECTION METHOD WITH SELF-CALIBRATION
In RF network theory, scattering parameters characterize RF devices as linear time invariant systems. The S-parameters S ij are complex quantities describing the ratio between outgoing traveling waves at port i when port j is excited [26]. By definition, the traveling waves are complex-valued such that S-parameters also include the phase shift. Here, it should be mentioned that the networks can also be described in terms of so called powerwaves [27], in analogy to electrical circuits. However, the two definitions mainly differ in the renormalization of the S-matrix, while the method described here is independent of the renormalization. In terms of RF network theory, a waveguide is modeled as a transmission line, i.e. a 2-port device which is reciprocal, symmetric and has zero reflection coefficients if it is impedance matched to its connected devices, Here, γ = α + iβ is the propagation constant, including the field attenuation coefficient α and the wavenumber β, and L is the length of the transmission line. The quantity of interest here is the wavenumber as it is related to the phase velocity via v ph = ω /β, with ω being the angular frequency. When short-circuiting a waveguide, the phase φ 11 := arg S 11 in reflection relates to the phase velocity via By convention the phase is given in the interval [−π, π) instead of [0, 2π). While so-called phase-unwrapping, which adds 2π at jumps in φ 11 (f ), allows for absolute phase relations across the frequency band of interest, the lowest frequency point lacks an absolute offset in multiples of 2π. If a certain dispersion shape is already preconditioned, the offset can be determined as a free fit parameter, for instance with the analytical dispersion of a perfect electric conductor waveguide. But without prior knowledge the offset remains unknown. One way to overcome this limitation is measuring an almost identical waveguide with slightly different length L [25]. Here, it is assumed that the waveguides have the same dispersion and that the length difference is on a sub-wavelength scale. Then, the phase φ 11 differs by which allows to deduce the phase velocity. However, manufacturing imperfections could limit the similarities between both structures. In order to characterize a single waveguide independently, it requires the short circuit to be freely movable inside. As this is technically very challenging the short can be replaced by some reflecting object, denoted as an obstacle. This replacement requires that all reflections behind the object are negligible, which limits systematic errors. Fig. 1(b) and Fig. 1(c) demonstrate the measurement principle. The reflecting obstacle causes a total phase shift at the input port of φ = 2l i β + φ offset , where l i is the length of the section from the input to the reflecting device, and φ offset is an arbitrary offset caused by the reflection. When moving the obstacle by a defined step ∆l the measured phase changes by ∆φ = 2∆lβ, providing the unknown β and, thereby, the phase velocity.
However, since the waveguide is embedded with its couplers, it cannot be characterized independently. An error network model is constructed, Fig. 2, which summarizes all elements upstream of the embedded device under test (DUT) as P -network. Then, the measured reflection coefficient S (m) 11 at the test port reads as, S (m) 11 = P 11 + Q 11 where P ij and Q 11 denote the S-parameters of the Pnetwork and the Q-network, respectively. S ij of the DUT are substituted by Eq. (1) in the derivation, but with a shorter effective length l, l < L, of the transmission line due to the obstacle being inserted. The goal of a calibration is to extract the DUT's S-parameters from Eq. (4). In fact, the network model is an extension of the standard three-term error model for the measurement of a 1-port device [28,29], which is calibrated via Open-Short-Match (OSM) standards. The terms P 11 , P 21 P 12 and P 22 correspond to directivity E D , source match E S and reflection tracking E R in the error terms representation, respectively. The DUT from the three-term model is replaced by the reflection coefficient k 2 Q 11 . Due to the additional unknown propagation constant the model is denoted as the four-term error model. However, since the DUT cannot be detached from the couplers, no calibration standards can be applied. Placing the obstacle at four different positions l i with steps smaller than the free space wavelength, is in principle sufficient in order to solve the system of equations for the four unknowns a, b, c and β. Here, the error terms from Eq. (4) are summarized as arbitrary unknowns a, b and c. However, it turns out to be numerically more stable to over-determine the system of equations and, instead, apply a non-linear least squares algorithm to the measured reflection coefficients S 11 (l) depending on the obstacle position. In order to include bounds for the solution of Eq. (5) the trust region reflective algorithm is chosen for the non-linear least-squares problem, as implemented in the scipy package [30,31]. In Eq. (5), the attenuation is constrained to be zero, α = 0. For most DUTs, especially the metallic waveguides, the attenuation is too small to be significant on the length scale of the obstacle position sweep. For the case of lossy waveguides the attenuation is significant, and can be included in the fit routine, as shown in Section IV B.
During a measurement, the Vector Network Analyzer (VNA) sweeps the excitation frequency within the band of interest, between 220 GHz to 330 GHz. This is repeated for each obstacle position such that the final data set includes both dependencies, S 11 (f, l). The least-squares problem is solved for each frequency point f j independently, which results in the phase velocity v ph (f j ) = 2πfj /βj.
In the next section the developed method is validated by simulations, before application to experimental data in the subsequent sections.

III. SIMULATION
In order to construct a network reasonably described by the error model from Section II a coupler-to-coupler transition was simulated with the time domain solver of CST Studio Suite 2019 [32]. The overall transition consists of a WR3.4 waveguide port attached to a pyramidal out-coupling horn, a free space section of 10 mm, another horn coupler of the same shape for in-coupling The rectangular waveguide has a standard cross section of 0.864 mm by 0.432 mm while the horn has an aperture cross section of 11.14 mm by 8.95 mm after a transition of 35 mm length. The dimensions are based on the optimal horn condition [33]. Here, for simplicity the out-and in-couplers were considered to be identical and the free space section is kept short for a smaller computational domain. While these two conditions are not met in the experiments in Section IV, the error network is capable of modeling both cases.
At the ports the horn-to-horn transition was excited in the frequency range of interest from 220 GHz to 330 GHz in order to compute the scattering parameters. These were further processed with the RF network modeling package scikit-rf [34][35][36]. The horn-to-horn transition is treated as the P -network from Fig. 2, while the Qnetwork is constructed from a randomly generated Sparameter, which is constant with respect to frequency. Both networks are cascaded with a transmission line of adjustable length l i in between. The transmission line is modeled with the geometrical parameters of the WR3.4 waveguide. Fig. 4 shows an example of the computed reflection coefficient in terms of real and imaginary parts depending on the obstacle position for a fixed frequency. While the periodicity of 1/2β is obvious, one can also observe that S 11 (l) does not follow a simple sine shape. The term a from Eq. (5) leads to an offset, while b and c distort the symmetry around a peak. Especially the real part shows a flatter slope on the left side while the slope to the right is steeper.
The goal of fitting the four-term error model is to extract the phase velocity and verify it with the expected one. In addition, the model outputs the error terms a, b, c from Eq. (5). To find a suitable starting point, β 0 , for the parameter β in the least squares algorithm, a fast Fourier transform is applied to S 11 (l) and the point of maximum amplitude in the reciprocal space is chosen. For the frequency point in Fig. 4 the fit results in a wavenumber β = 5.209 mm −1 which gives v ph = 1.220 c. In Fig. 5 the phase velocity is plotted over the whole frequency range considered, together with its relative deviation from the The extracted values do not differ from the input reference values by more than 9 × 10 −5 %. Based on the scattering parameters of P and Q used in constructing the network, the error terms can also be compared to the prediction, where a = P 11 is expected as well as b = Q 11 P 21 P 12 and c = P 22 Q 11 . The fit parameter a is the most precise one of the three as it describes the offset of S 11 (l). Across the whole frequency range it deviates by 0.04 % on average. At frequencies where P 11 ≈ 0 the relative deviation increases significantly to about 3.4 %. The other two fit parameters b and c have a similar effect on the shape of S 11 (l) such that the fit routine can determine them less precisely. As before, their deviations from the expected value is largest when they approach zero. On average, they deviate by ∼ 7 %. Since the fit parameter of interest is β, deviations up to 10 % in the error terms are acceptable. With the preceding simulation, the four term error model has demonstrated its applicability in measuring the phase velocity of an integrated waveguide. Here, the WR3.4 waveguide supports only the fundamental mode which is why the horns are placed in close proximity. For unknown devices the phase front error may couple power to higher order modes if present. In the experiment this is mitigated by placing the horn's further apart from each other.

IV. EXPERIMENT
The setup consists of a Rohde & Schwartz ZVA67 VNA which is connected to a ZC330 millimeter-wave converter spanning the frequency range between 220 GHz and 330 GHz. The converter provides a rectangular waveguide output port of type WR3.4 to which a standard gain horn antenna is attached (Anteral, SGH-26). In order to compensate internal systematic errors the output port was calibrated via a standard OSM calibration. Errors up to the waveguide flange, such as directivity errors from the cables attached, are not covered by the above introduced error model. So calibrating up to the flange moves the reference plane for the S-parameter measurement to the connector between waveguide and horn antenna. The monolithic DUT was mounted on an optical post in line with the beam axis of the out-coupling horn. A rail allows different positioning along the axis which helps with in-coupling. The free-space distance of ∼ 5 cm was chosen as compromise between lower phase error at further distances and higher incident power at closer distances.
The obstacle was mounted on a motorized stage on top of the rail, downstream of the structure. The choice for obstacles is restricted by being non-contacting movable inside small waveguide apertures in the order of ≈ 1 mm diameter. Thin stiff devices such as a blunt needle of a dosing syringe, and blue steel wires were applied. The gauge 22 needle has an outer diameter of 0.72 mm and nominal inner diameter of 0.51 mm. The set of blue steel wires covers a range from OD = 0.3 mm to OD = 4 mm while only wires with OD = 0.50 mm, 0.81 mm, 1.21 mm and 1.40 mm have been used in the present study so far. All measurements have been conducted at ambient air at 21°C and relative humidity of 50 %. Dry air contributes to the refractivity with N dry = (n − 1) × 10 6 ≈ 270 [37, Sec. 1.3.5.1], while the frequency independent part of water vapor contributes N H2O ≈ 70 [38]. Although water vapor provides 49 absorption lines in the frequency range of interest [39,40], even the highest one at 325 GHz contributes only by about N 325 GHz ≈ 0.3, based on the phase shift lineshapes [41,42].

A. Results for metallic waveguides
Three circular metallic waveguides have been tested, and compared to predictions based on the known analytical dispersion relation [43], where f c is the cutoff frequency, a the radius of the cross section, and u 11 the first root of the derivative of the first order Bessel function, J 1 (u 11 ) = 0, due to the fundamental mode being the T E 11 mode. Although the measurements are not conducted in vacuum, it is assumed for simplicity that the refractive index n = √ r = 1. The error introduced by the assumption is in the order of 0.03 %, which is below the current limitations in experimental uncertainties. Microscope images of the waveguide apertures were taken and analyzed by a circle detection, using the scikit-image package [44]. The algorithm performs a Hough transform [45] of the detected edges in the image and returns a score for each pair of pixel and radius. The pair with the highest score is identified as the center and the radius of the aperture. The peak width compared to the background of the score is considered as the uncertainty. Fig. 6 shows the three structures, the first one being a copper block with the waveguide being drilled, the second one was also drilled but in brass, and the third one is again a copper structure which was wire electron-discharge machined (EDM). Due to machining limitations, the structure was manufactured in two half shells which were precisely assembled afterwards with alignment pins. This structure is referred to as split block waveguide. The two different machining techniques have significantly different final features. A ridge formation is observable from the bored structures while the EDM provides a much more explicit edge.
The scans of the obstacle position covered a range of 10.00 mm inside the hollow channel, with minimum distance to the free space transition of 5 mm. This may be important for the assumptions on the error model, as will be explained in the discussion section. Measured dispersion relations v ph (f ) for the two block waveguides are shown in Figs. 7 and 8. For both cases, applying a least-squares routine with Eq. (7) and the radius as free parameter led to reasonable fits and effective radii of 498 µm and 495 µm, respectively. The dispersion unravels the actual geometry inside, in contrast to the radii of (511 ± 10) µm and (505 ± 7) µm measured at the apertures. While the measurements do not match within the uncertainty range around the prediction, the almost identical shape indicates a systematic error leading to the offset between prediction and measurement. Due to the drilling artifacts on the channel opening observed on the microscope images, the systematic error is attributed to the measurement of the radii. It is reasonable to assume that the channel widens close to the opening as the drill head was less constraint in space. Both structures were measured with the blunt needle as obstacle as well as an 0.80 mm blue steel wire. Deviations in v ph (f ) from different obstacles were at maximum 0.3 %.
In contrast, the measured dispersion v ph of the split block waveguide, Fig. 9, matches very well the prediction based on (657 ± 7) µm for frequencies below 280 GHz. The fit of the measurement with the analytical model, Eq. (7), results in a radius of a = 654 µm which is well within the uncertainty range of the radius measured by microscopy. The standard deviation of the fit error is below 0.1 µm which shows how well the measurement and model match. However, above 280 GHz the measured phase velocity clearly deviates from the fit, as well as from the prediction. The reason for that became clear when looking at the scan data for a single frequency point, for instance at 283 GHz. Although not shown here explicitly, S 11 (l) clearly follows a beat pattern, indicating a second mode being present. When applying a discrete Fourier transform, two peaks of almost identical amplitude were identified in the reciprocal space. This explains why the model of Eq. (5) fails, as it assumes a single mode propagation in the waveguide. The discrete Fourier transform is insufficient in resolving a smooth change in β(f ) since the reciprocal space is sampled in discrete points of 1/(N ∆l), where N refers to the number of obstacle positions and ∆l is the step size. Zero-padding shifts the spectral sampling points, but does not improve the resolution and is still limited by the truncation artifacts. Filter diagonalization method (FDM) [46,47] is not limited to discrete steps and therefore used to study the multi-mode behavior of the waveguide, using the harminv [48] implementation. For interested readers, section 2 of Ref. [49] provides a more rigorous comparison of the discrete Fourier transform and the FDM. As the FDM is modeling the input signal as a sum of sinusoids, it is not able to properly account for the errors introduced in Sec. II. But as seen before in Fig. 4, the 1/(2β) periodicity in S 11 (l) is at least similar to a sinusoid if the error terms b and c are small. Fig. 10 compares the modes found by FDM with analytical expected dispersion relations of the next higher order modes. First, the FDM finds several modes which are not following any dispersion line and seem to be spurious. This may be due to the sinusoids not properly representing the error model. Second, the majority of modes follow dispersion lines with respect to frequency. As expected, the fundamental T E 11 mode dominates and is present over the entire frequency range. Above 280 GHz, a second mode is present which ranges from v ph ≈ 2.2 c to ≈ 1.4 c. Even a third dispersion line is observed above 295 GHz. Finally, the second dispersion line does not match with an analytical mode. In contrast to the first and the third line, which coincide very well with T E 11 and T E 01 modes, the discrepancy between the second line and the expected T E 21 mode has not yet been understood.   Fig. 6(b). Prediction based on Fig. 6(e).

B. Characterization of unknown waveguides
While the previous section compared measured dispersions of known DUTs to the analytical predictions, the method was also applied to devices with unknown dispersion. Regarding phase synchronous mode propagation with particles, purely metallic waveguides are unsuitable as v ph > c. By introducing a dielectric lining the modes mainly become hybrid, but with reduced phase velocity. With a given set of the dielectric permittivity ε r , the inner hollow core radius a, and outer radius b of the metallic wall, the dispersion can also be computed analytically [15]. However, this requires a precise knowledge of the three parameters which is often not feasible, and a direct measurement of the dispersion is necessary. As a first unknown device, the copper block waveguide from Fig. 7 was oxidized on the inside by browning (Ballistol Nerofor). The thickness of the oxide layer is unknown as well as its permittivity in the frequency range of interest. The measurement was conducted with the same settings as before oxidization. Fig. 11 shows the dispersion in comparison to the one of the purely metallic waveguide. As expected, the phase velocity has decreased. The drop in v ph is about 0.04c at 220 GHz and decreases towards higher frequencies, i.e. about 0.01c at 290 GHz. Here, the measurement above 290 GHz was clipped. Significant loss above this threshold led to a poor fit of the model in Eq. (4). This loss was already present at low frequencies, while being small enough to allow for a good fit. In Fig. 12  where the exponential decay is observed. Note that, here, the position l obs refers to the distance the obstacle has been moved in from the backside, which is why the attenuation is represented as an exponential growth. The fit of the four term error model includes explicitly a non-zero attenuation coefficient, β = (2.8280 ± 0.0009) mm −1 , α = (0.0250 ± 0.0009) mm −1 . The envelope does not properly cover the maxima of the fit since setting β = 0 does not account for the increase and shift of the amplitude due to the error terms b and c. But the envelope follows very well the trend of the measured data. In addition to attenuation, at high frequencies the error model also fails due to higher order modes. While the oxide layer decreased the phase velocity, it is not sufficient to reach the phase synchronous case of v ph ≤ c. A dielectric loading was created in the split-block waveguide by inserting a 3D-printed polymer capillary. The relative permittivity of the resin used for printing (Moiin Tech clear) was measured in advance [50], being ε r ≈ 2.78 − i0.07 in the frequency range of interest. The wall thickness was estimated as 150 µm. Both values are too inprecise in order to predict the dispersion relation, and the printing process was not fully optimized such that support residuals on the capillary distort the ideal shape, see the inlet of Fig. 13.
The wire with a diameter of 0.5 mm was used as an obstacle and scanned over a 10 mm range inside the waveguide. The resulting phase velocity dispersion is shown in Fig. 13. v ph ranges between 1.10 c at 220 GHz and 0.98 c at 300 GHz and crosses v ph = c at 275 GHz.

C. Uncertainty analysis
Although not visible on the chosen axis scale, dispersion lines presented in the previous section include an uncertainty band around v ph (f ) in the plots. This demonstrates The fit model is extended with a non-vanishing attenuation coefficient α. The envelope is based on β = 0. Note that the generator is to the right side of the plot which is why the attenuation appears as an exponential growth. already the level of accuracy of the measurements.
To assess the uncertainties, each frequency sweep was repeated 10 times and the sample mean S 11 was gathered together with a 99 % confidence interval. Due to the low sampling number, the sample standard error was weighted with the proper Student's t of 3.25. The absolute error in |S 11 | did not exceed 1 × 10 −3 across the entire frequency band of interest, while the phase error was about 0.03°.
To assess the uncertainty for β and v ph in an obstacle scan, the dataset S 11 (l) with its uncertainty was Monte Carlo sampled, assuming normal distribution. Each randomized copy of the dataset was applied to the fit routine of Eq. (5) to extract β and its standard deviation. The sampling and fitting was repeated for each frequency point of S 11 (l, f ). For all measurements presented, relative uncertainties in the phase velocity are below 0.05 %. However, they are not shown in the plots. The standard deviation in the parameter estimate of the fit was about an order of magnitude larger. This error is not a figure of merit for the measurement uncertainty, but a measure of how well the model describes the dataset. Systematic errors due to the limits of the four-term error model dominate over statistical uncertainties.
In cases in which the mismatch of the model and the data was more severe, the fit error was substantially larger. In Fig. 9, above 280 GHz the error band enveloping the dispersion line is clearly visible and covers a range of ∆v ph ≈ 0.02c. The following section discusses the limits of the experimental setup and the model.

V. DISCUSSION
Here, the limitations of the underlying model are discussed first. Potential improvements in the experimental design and the DUTs are considered in the subsequent section.

A. Limits of the error model
The most limiting assumption in the network model of Fig. 2 is the obstacle being a single port device. Modes propagating further downstream of the obstacle are supposed to be absorbed or radiated away at the end of the transmission line, such that no reflections behind the obstacle interfer with the reflections from the Q-network. An ideal obstacle for which the assumption is always true, would be a short circuit. No power would propagate behind the obstacle. Such devices exist in phase shifters for standard waveguides, e.g. WR-3.4. But for non-standard arbitrary shape waveguides this is technically challenging as the short device has to be non-contacting and movable.
Considering propagation behind the obstacle, the wire and the wall form another transmission line of different wave impedance and propagation constant. The transition at the waveguide exit causes another reflection and is considered as a one-port device. In the network model, Q extends to a two-port device, including the impedance mismatch between the two lines. In total, this makes 5 additional unknown parameters which would have to be included in the measured S of Eq. (4). The extended nine-term error model is a better representation of the experimental setup, but infeasible to solve. The least squares algorithm would have to optimize 16 real parameters, as the error terms are complex numbers.
To suppress mode propagation downstream, the diameter of the obstacle should be as large as possible without scratching the wall. In Fig. 13 the dispersion line is not as smooth as in the other measurement, which is attributed to interference with downstream reflections. Here, the wire is only 0.5 mm thick which is only half the size of the hollow core. In the other measurements the remaining air gap between obstacle and walls was half the size, about ≈ 100 µm.
Another limiting aspect has been shown in Figs. 9 and 10, the excitation of higher order modes. The beat pattern was observed in several obstacle scans, but often small enough to be negligible, e.g. as in Fig. 12. The fourterm error model is designed for single mode operation. Small higher order mode contributions in the reflection coefficient lead to systematic errors in the fit. An extension requires a summation over all modes excited by the current frequency, including error terms a i , b i , c i and wavenumber β i where i = 0, 1, 2, . . . refers to the mode order. Since all the terms are unknown the sum cannot be fitted anymore in general. For the special case of only two modes, one could approach the problem by fitting the fundamental mode first, subtracting the fit from the data, and apply the model again to the remaining signal for the second mode. Subtracting the fit result from the original data and fitting the remaining set again refines the parameter estimates for the fundamental mode. By iterating between the two modes, the error terms could be determined more accurately.
However, extending this scheme for more than two modes is not straight-forward. Assuming a sinusoidal shape of S 11 (l), as it was done with applying the FDM, is useful for finding a first estimate. It is sufficient to compare the dispersions of different modes. But as it does not properly include the asymmetric distortion due to the terms b and c, it may fail in precision. Comparing the analysis presented in Fig. 9 and Fig. 10 up to a frequency of 250 GHz, they deviate by only 0.08 % on average. But at few selected points the deviation is on the order of ∼ 1 %, showing the discrepancy in accuracy between the error model analysis and the FDM analysis.
Transitions also introduce systematical errors which is why the obstacle position is supposed to be scanned at least three wavelengths away from the coupler and the waveguide exit. As a test, the capillary loaded waveguide was studied with a scan covering partially the metalliconly section, the transition, and the dielectric loaded section. The phase velocity was determined with different truncations of the scan S 11 (l) relative to transition, one directly starting from the transition, one being 2 mm away, and one being 3 mm away from the metallic section. v ph differed about 0.01c between each case. With further increase in distance no significant change was observed. The found transition region is only valid for the certain case of step-transition and cannot be generalized to arbitrary transitions. But it gives already a zeroth order estimate for the impact on the measurement.

B. Experimental and device limits
Besides the methodological limitations, the setup and the structures under test are not ideal. The setup is not optimized for maximum coupling. Usually, horncouplers are often designed for ideal plane wave coupling, i.e. maximum aperture efficiency as considered for the optimum horn condition [51,52]. Further improvement in coupling is achieved by extending it to Gaussian beam optics [53]. However, there is no universal condition. The ideal horn design depends on the incident Gaussian beam, or the quasi-optical layout has to match the given horn design.
Due to the lack of Gaussian optics in the current setup the first test devices, Fig. 6(a) and Fig. 6(b), were based on the optimal horn condition. The third one is manufactured with a different horn geometry to match the optical layout of another setup. In both cases, the horn not only captures the incident power, but also causes a phase shift. The effect was simulated with CST Microwave Studio for an exemplary case. A waveguide was excited by an incident Gaussian beam, once with and once without horn antenna. The phase at the output port was compared for both cases, showing a difference of ≈ 110°, almost constant over the excited frequency band from 265 GHz to 275 GHz. For a 20 mm long waveguide this leads to an error of ≈ 1.5 % in the phase velocity. The current scanning reflection method overcomes this systematic error by attributing the horn-related phase shift to the error terms.
With the aim for phase synchronous mode propagation and efficient THz-particle interaction in an accelerator based application, the dielectric lining requires significant manufacturing improvements. On the one hand, mechanical residues from 3D-printing have to be reduced. The circular shape of the hollow core is significantly distorted along the longitudinal axis, causing scattering of the propagating wave. On the other hand, the polymer is lossy and requires further treatment to reduce absorption. Different materials such as fused silica, are more desirable, also because of their higher rf-breakdown limits. But they also increase the machining complexity for an integrated coupling device.

VI. OUTLOOK
Due to coupling to the fundamental mode, the current setup characterizes modes suitable for transverse deflection of bunches [8]. For characterizing the accelerating T M 01 mode the setup will be equipped with a segmented half-wave plate to match the incident field pattern to the waveguide pattern. PTFE lenses are added for focusing the beam into the coupler. However, both optical elements reduce the spanned frequency bandwidth due to being rather monochromatic.
Future structures to be studied include another polymer based dielectric lined waveguide, but instead of being inserted into a metallic structure, the combined polymer horn capillary will be coated with a metallic layer. In addition, fused silica capillaries will be coated and inserted in metallic horn structures.
Potentially more complex devices, such as a tapered dielectric waveguide for synchronous acceleration of low-energy particles [1], or a sectioned waveguide for longitudinal phase space synthesis [54], could be explored with the help of windowing the obstacle scan, resulting in a spectrogram representing the change of β as it varies with position.

VII. CONCLUSION
A new error model is developed with the purpose of selfcalibrating reflection measurements from monolithic integrated hollow waveguides via free-space coupling. Reflections from unmatched couplers, from surrounding mounting elements as well as free-space losses and phase jumps at the reflecting movable obstacle are included in the error-terms. The model is applied to a simulation to investigate the deviations from the input error-terms and phase velocity.
Experiments with metallic waveguides of known cross section are in good agreement with analytical predictions. While the measurements provide a precision on the level of ≈ 0.05 %, systematic errors not covered by the error model currently limit the accuracy to ≈ 0.5 %. One of the systematic errors is introduced by reflections downstream of the obstacle. Comparison measurements between a cannula and a blue steel wire as obstacles deviate up to < 0.3 %. A dedicated choke-type piston as movable short may reduce the downstream reflections in future. In addition, the error model breaks down in case of higher order mode excitation.
The method was also applied to dielectric-lined metallic waveguides with both a copper oxide and a polymer lining, yielding reduced phase velocities. In the case of the thicker 3d-printed polymer, we measured a phase velocity equal to the speed of light at a frequency of ∼ 275 GHz, compatible for efficient beam acceleration and manipulation with relativistic beams.

SUPPLEMENTAL MATERIAL
For reproducibility a python-based implementation of the 4-term error model used for the analysis is attached to the paper, as well as the raw data sets of the presented experimental results.
Appendix A: Derivation of the error model The waves between the networks in the flow graph Fig. 2 are denoted as (a 1 , b 1 ) before the P -network, (a 2 , b 2 ) between P -and S-network, and (a 3 , b 3 ) between S-and Qnetwork. For each interconnection the reflection coefficient from the right-sided device is denoted as ρ i = bi /ai. When going from right to left in the network, each subnetwork represents a Möbius transform of the reflection coefficient. Applying the transform from interconnection 3 to 2, The transform from interconnection 2 to 1 results in b 1 a 2 = P 11 P 12 P 21 P 22 a 1 b 2 ⇒ ρ 1 = P 11 a 1 + P 12 b 2 a 1 (A3) ρ 1 = P 11 + ρ 2 P 21 P 12 (1 − P 22 ρ 2 ) = P 11 + Q 11 P 21 P 12 k −2 − P 22 Q 11 (A4)